Shrinkage Strategies for Structure Selection and Identification of Piecewise Affine Models

Many real-world systems and phenomena are characterized by different operating conditions, among which they can commute smoothly and/or abruptly. In such scenarios, linear models are usually insufficient to accurately describe the behavior of the underlying system and one must resort to more complex model structures to attain satisfactory approximations. However, many learning frameworks lead to models that are rather hard to analyze and use, especially for control purposes.

A good trade-off between complexity, approximation power and intelligibility can be attained by considering piecewise affine (PWA) maps. Indeed, it is well known that PWA functions have universal approximation properties [13, 7], while featuring a rather simple structure, since they are characterized by a finite collection of affine sub-models defined over a polyhedral partition of the regressor space. Thus, PWA models can be analyzed and exploited for control design with state-of-the-art techniques (see e.g., [5]).

Although appealing because of the aforementioned properties, PWA models are challenging to learn from data. Indeed, to identify the local models and estimate the associated regions of the regressor space one should know the operating condition of the real systems at each time step, which is usually not available a priori. The learning problem thus entails not only the identification of sub-models parameters and the computation of a polyhedral partition of the regressor space, but it also requires one to solve an unsupervised clustering problem, since each sample has to be assigned to one of the local models.

Several heuristics have been proposed over the years to deal with this demanding learning task, among which we mention the bounded-error method proposed in [4], the clustering-based procedures in [2, 8, 9, 15], the mixed-integer programming based approaches in [14, 17], the Bayesian method in [11] and the sparse optimization techniques proposed in [1, 16]. However, most of these methods rely on the prior knowledge on the number of operating conditions of the true system and of its order, which might not be available in practice. In [4], the number of sub-models is chosen by solving a minimum number of feasible subsystems (MIN PFS) problem, which is, in turn, formulated by relying on a bound for the estimation error. This parameter is unlikely to be available in practice, thus becoming a tuning knob for the algorithm. A refinement strategy is also proposed therein, where the number of sub-models is further reduced by merging local models that share the same parameters and discarding the ones that are associated to “few” data-points. By exploiting sparse optimization, the methods presented in [1, 16] both allows for the selection of the number of modes. Specifically, the approach proposed in [16] considers an over-parameterized model and, then, exploits $L_{1}$-regularization to enforce the local parameter estimates to be similar. Instead, in [1], $L_{1}$-regularization is used to iteratively search for the sparsest vector of fitting errors to cluster the data-point, estimate local models and decide on their number. To the best of our knowledge, only in [1], it is provided a possible solution for the selection of the sub-models structure, which relies on the introduction of an additional $L_{1}$-regularization term on the parameters.

As an alternative, in this paper we present two regularization-based strategies that allow the partial selection of the PWA model structure along with parameter estimation. These techniques exploit the flexibility of the general purpose approach for jump model fitting presented in [3], which can easily incorporate regularization-based shrinking strategies. In particular, we propose a heuristic to estimate the number of local models by using $L_{\infty}$ regularization on the sub-models parameters, while accounting for the cardinality of the associated clusters to refine the optimization-driven choice. In the second approach, the local model structure is estimated via the elastic net. We remark that the proposed strategy is inspired by the one briefly discussed in [1], but it inherits the benefits of elastic net over simple lasso regularization. The overall identification procedure with partial structure selection requires the alternation between a learning stage and a step in which the model is refined. Once the model structure is determined, the final model is re-estimated by neglecting the shrinkage regularization terms. Although the proposed strategies are quite general and can be straightforwardly extended to other classes of PWA models, in this work we focus on piecewise affine autoregressive models with exogenous inputs (PWARX).

The paper is organized as follows. The problem of learning a PWARX model with partially unknown structure is introduced in Section II and the core algorithm used to address this learning task is summarized in Section III. The proposed regularization-based strategies are then discussed in Section IV, and their effectiveness is shown in Section V on a numerical example and an experimental case study. Conclusions and directions for future work are finally discussed in Section VI.

Consider the following dynamical data-generating system, whose instantaneous response $y_{t}\in\mathcal{Y}\subseteq\mathbb{R}$ to an input $u_{t}\in\mathcal{U}\subseteq\mathbb{R}$ is given by

with $n_{a},n_{b}\in\mathbb{R}^{+}$ denoting the possibly unknown order of the system. The unknown function $f_{\mathrm{o}}:\mathcal{X}\rightarrow\mathbb{R}$ characterizing the underlying regressor/output relationship is assumed to be nonlinear and possibly discontinuous, but not necessarily piecewise affine.

We aim at finding a PWA approximation $f$ of the unknown function $f_{\mathrm{o}}$, based on a set of inputs $U=\{u_{t}\}_{t=1}^{T}$ and output samples $Y=\{y_{t}\}_{t=1}^{T}$, generated by (1a)-(1b). We remark that, due to the universal approximation properties of PWA maps [13, 7], PWA function $f$ can approximate $f_{\mathrm{o}}$ with an arbitrary accuracy.

The approximating PWA function $f:\mathcal{X}\rightarrow\mathbb{R}$ is defined as

with $\tilde{x}_{t}=\left[\begin{smallmatrix}x_{t}^{\prime}&1\end{smallmatrix}\right]^{\prime}$. This piecewise affine autoregressive model with exogenous inputs (PWARX) is a collection of $K\in\mathbb{N}$ local affine models, characterized by the parameters $\Theta_{y}=(\theta_{y,1},\ldots,\theta_{y,K})$ and defined over a complete polyhedral partition¹¹1$\{\mathcal{X}_{k}\}_{k=1}^{K}$ is a complete polyhedral partition of $\mathcal{X}$ if $\bigcup_{k=1}^{K}\mathcal{X}_{k}=\mathcal{X}$ and $\overset{\circ}{\mathcal{X}_{i}}\cap\overset{\circ}{\mathcal{X}_{j}}=\emptyset$, for $i\neq j$, $i,j=1,\ldots,K$, with $\overset{\circ}{\mathcal{X}_{i}}$ denoting the interior of the $i$-th polyhedron $\mathcal{X}_{i}$. $\{\mathcal{X}_{k}\}_{k=1}^{K}$ of the regressor space $\mathcal{X}$. We stress that all the local models share the same structure by definition.

We characterize the partition through a piecewise affine separator function $\phi:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}$, defined as

where $\theta_{x,k}^{\prime}\in\mathbb{R}^{n_{x}+1}$. Thus, a polyhedron $\mathcal{X}_{i}$ can be described by the following linear inequalities

Therefore, learning a PWARX model from data amounts at: $(i)$ choosing the number of affine sub-models $K$; $(ii)$ estimating the parameters $\Theta_{y}$ of the local models, $(iii)$ finding the parameters $\Theta_{x}=(\theta_{x,1},\ldots,\theta_{x,K})$ characterizing the PWA separator in (3), and, eventually, $(iv)$ selecting the order of the sub-models, namely $n_{a}$ and $n_{b}$, when it is unknown.

In this work, we split model structure selection into two separate tasks by assuming that either the number of local models $K$ or the sub-model order is fixed a priori. In the first task, we estimate $(\Theta_{y},\Theta_{x},K)$ for fixed $n_{a},n_{b}$, while in the second one we fix $K$ and retrieve $(\Theta_{y},\Theta_{x},n_{a},n_{b})$ from data.

It is well known that learning PWA approximating map is quite challenging, since this problem is NP-hard [12] even for an a-priori fixed model structure. Indeed, it requires the concurrent solution of local linear regressions, a classification and a clustering problem. To cope with this, we recall the modular optimization problem presented in [3, Section 2.2.1.] for PWA regression and we summarize the steps of the coordinate-descent approach proposed therein to solve the problem at hand.

Let $X=\{x_{t}\}_{t=1}^{T}$ be the sequence of regressors constructed from the available input/output dataset $\{U,Y\}$. Let $s_{t}\in\{1,\ldots,K\}$ be the latent variable, that indicates which local model is associated with the $t$-th regressor/output pair $\{\tilde{x}_{t},y_{t}\}$, such that the sequence $\mathcal{S}=\{s_{t}\}_{t=1}^{T}$ comprises the label associated to each data point.

Within the optimization-based framework of [3], learning a generic jump model entails the solution of the following optimization problem:

From (5a), it is clear that the discrete sequence $\mathcal{S}$ couples the regularized linear regressions and the multi-category discrimination problems, that have to be solved to find $\Theta_{y}$ and $\Theta_{x}$, respectively. Nonetheless, once $\mathcal{S}$ is fixed, the parameters $\Theta_{y}$ and $\Theta_{x}$ can be found independently from one another. On the other hand, for a given $\Theta_{y}$ and $\Theta_{x}$, the discrete sequence $\mathcal{S}$ can be computed using dynamic programming (DP) [6] technique in a computationally efficient manner.

Based on the above considerations, we exploit the coordinate-descent approach proposed in [3] by alternatively solving problem (5a) with respect to $\Theta_{y}$ and $\Theta_{x}$, for a fixed sequence $\mathcal{S}$, and optimizing with respect to $\mathcal{S}$ for fixed $\Theta_{y}$ and $\Theta_{x}$. The three step procedure tailored for the solution of problem (5a) is summarized in Algorithm 1 where, with a slight abuse of notation, we denote

Given the structure of the local models, the problem ought to be solved at step 0.0..0..0. can be split into $K$ distinct regularized linear regression problems, one for each local models. We remark that each sub-model is updated by using only the data points that are associated to it, according to the estimated mode sequence $\mathcal{S}^{i-1}$. The estimation of the PWA separator at step 0.0..0..0. entails the solution of a multi-category discrimination problem, which can be efficiently computed with the Newton-like method proposed in [8]. The mode sequence is finally updated at step 0.0..0..0. via DP [6], as detailed in [3], by neglecting the regularization terms as they are independent of $\mathcal{S}$. Note that, the tunable parameter $\rho\in\mathbb{R}^{+}$ is not considered when optimizing with respect to $\Theta_{y}$ and $\Theta_{x}$ (steps 0.0..0..0.-0.0..0..0.), but it is exploited when updating the mode sequence to trade-off between the fitting error and the cost of choosing a certain polyhedral region.

Suppose that either the number of sub-models characterizing the PWA map in (2) or their local structure is unknown and assume that bounds on the values of these unknowns are given by the maximum tolerable complexity of the estimated model.

To determine the unknown PWARX model structure (i.e., either the number of sub-models $K$ or their model order $n_{a}$, $n_{b}$), the regularization terms in the objective function of (5a) can be shaped so to use well-known regularization-based shrinkage techniques to choose the structure of the estimated PWARX model in (2), directly from data. This allows to avoid grid searches, which usually require an exhaustive exploration of the unknown parameter space.

In this work, starting from an over-parameterized model, we aim at devising iterative structure selection strategies that rely on runs of Algorithm 1 with tailored regularization terms. In particular, we focus on the selection of the group regularization $r_{y}(\Theta_{y})$, which acts on the parameters of the local models.

Two different choices of $r_{y}(\Theta_{y})$ are made to tackle the problems of estimating the number of sub-models $K$ and the model order $n_{a},n_{b}$ respectively. For both problems, we fix the regularizer $r_{x}(\Theta_{x})$ acting on the parameters $\Theta_{x}$ of the PWA seperator as follows

with $\lambda\in\mathbb{R}^{+}$ being a hyper-parameter to be tuned. This choice makes the multi-category discrimination optimization at step 0.0..0..0. of Algorithm 1 strictly convex, thus reducing the complexity of the resulting model. In the following, we introduce two different strategies based on the choice of $r_{y}(\Theta_{y})$, respectively tailored to handle the selection of $K$, for fixed local model structure, and the choice of $n_{a}$ and $n_{b}$, for a given number of local models. In the latter scenario, we assume that $n_{a},n_{b}\geq 1$, namely that the response of the system at a given instant always depends at least on the previous input/output pair. In both cases, we assume that the dataset is relatively balanced, namely we expect that the local models are associated to a reasonable number of samples for the overall model to be identifiable with sufficient accuracy.

In this section, the performance of the proposed model selection strategies is assessed on two benchmark examples, taken from [4]. The first academic example involves the identification of the PWARX system considered in [4], while the second case study involves modeling the dynamics of a placement process in a pick-and-place machine using experimental data [11]. In both case studies, the model structure selection and identification procedures are performed starting from over-parametrized models with $K_{max}=10$ or $n_{a,max}=n_{b,max}=10$, respectively. The tuning parameters $\mu$ and $\nu$ in the regularization terms (9) and (13) are selected so that $\mu\in(0,1)$ and

respectively. In this way, we reduce the number of tunable parameters, by considering a convex combination of the $L_{2}$ and lasso regularizations.

All computations are carried out an a $i7$ $2.8$ GHz Intel core processor running MATLAB 2019a. The local regression problem in Algorithm 1 are solved by using the CVX package [10], while the PWA separator is estimated via the computationally efficient batch approach proposed in [8], by exploiting the auxiliary parameters used in the example reported therein. To improve the quality of the solution and reduce its dependence on the initial conditions, Algorithm 1 is always executed for $N=20$ different initial mode sequences $\mathcal{S}^{0}$ and the PWARX model is selected as the one leading to the minimum cost.

The paper has presented a preliminary study on the application of regularization-based strategies to select the model structure of a PWARX model. The proposed approaches combine the selective power of existing shrinkage techniques with heuristics, thus providing a tool for data-driven model structure selection. The results obtained in the considered examples highlight the potential of these strategy.

Future research will be devoted to find less conservative criteria for the selection of the local model structure and to address a complete structure selection problem. The sensitivity of the presented algorithms w.r.t. a set of hyper-parameters will be analyzed and auto-tuning procedures for this parameters will be investigated. Thanks to the generality of the coordinate-descent approach on which our strategies rely, future work will also include their extension to other classes of switching models.