# Exploiting Model Sparsity in Adaptive MPC: A Compressed Sensing Viewpoint

This paper proposes an Adaptive Stochastic Model Predictive Control (MPC) strategy for stable linear time-invariant systems in the presence of bounded disturbances. We consider multi-input, multi-output systems that can be expressed by a Finite Impulse Response (FIR) model. The parameters of the FIR model corresponding to each output are unknown but assumed sparse. We estimate these parameters using the Recursive Least Squares algorithm. The estimates are then improved using set-based bounds obtained by solving the Basis Pursuit Denoising [1] problem. Our approach is able to handle hard input constraints and probabilistic output constraints. Using tools from distributionally robust optimization, we reformulate the probabilistic output constraints as tractable convex second-order cone constraints, which enables us to pose our MPC design task as a convex optimization problem. The efficacy of the developed algorithm is highlighted with a thorough numerical example, where we demonstrate performance gain over the counterpart algorithm of [2], which does not utilize the sparsity information of the system impulse response parameters during control design.

## Authors

• 10 publications
• 4 publications
• ### Learning to Satisfy Unknown Constraints in Iterative MPC

We propose a control design method for linear time-invariant systems tha...
06/09/2020 ∙ by Monimoy Bujarbaruah, et al. ∙ 0

• ### RLO-MPC: Robust Learning-Based Output Feedback MPC for Improving the Performance of Uncertain Systems in Iterative Tasks

In this work we address the problem of performing a repetitive task when...
10/01/2021 ∙ by Lukas Brunke, et al. ∙ 0

• ### Input Convex Neural Networks for Building MPC

Model Predictive Control in buildings can significantly reduce their ene...
11/26/2020 ∙ by Felix Bünning, et al. ∙ 0

• ### Learning Robustness with Bounded Failure: An Iterative MPC Approach

We propose an approach to design a Model Predictive Controller (MPC) for...
11/22/2019 ∙ by Monimoy Bujarbaruah, et al. ∙ 0

• ### Online Learning Based Risk-Averse Stochastic MPC of Constrained Linear Uncertain Systems

This paper investigates the problem of designing data-driven stochastic ...
11/20/2020 ∙ by Chao Ning, et al. ∙ 0

• ### Design of multiplicative watermarking against covert attacks

This paper addresses the design of an active cyberattack detection archi...
10/01/2021 ∙ by Alexander J. Gallo, et al. ∙ 0

• ### Reduction of the Number of Variables in Parametric Constrained Least-Squares Problems

For linearly constrained least-squares problems that depend on a vector ...
12/18/2020 ∙ by Alberto Bemporad, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

If the true dynamics of a system are uncertain, Adaptive Control strategies have been applied for meeting control objectives and ensuring system stability [krstic1995nonlinear, ioannou1996robust, sastry2011adaptive], typically under no output and input constraints. The uncertainty in these systems can be primarily attributed to two factors: model uncertainty (e.g. modeling mismatch and inaccuracies), and exogenous disturbances (e.g. sensor noise). More recently, utilizing tools from classical Adaptive Control, Adaptive Model Predictive Control (MPC) [tanaskovic2014adaptive, lorenzenAutomaticaAMPC, bujarbaruahAdapFIR, kohlerNonlAMPC] has established itself as a promising approach for control of uncertain systems subject to input and output constraints. For linear systems specifically, the literature on Adaptive MPC has extensively focused on either robust or probabilistic satisfaction of such imposed constraints on the system, using either state-space or input-output modeling.

In [soloperto2018learning, bujarArxivAdap] additive model uncertainties are considered with known system matrices, and imposed state and input constraints are robustly satisfied for all such realizable uncertainties. In [lorenzenAutomaticaAMPC, hernandez2018stabilizing, kohler2019linear] robust state-input constraint satisfaction is extended to include both unknown system dynamics matrices and additive disturbances. The approach introduced in [bujarArxivAdap] is also suited for satisfaction of probabilistic chance constraints on system states. Furthermore, the work in [hewing2017cautious, koller2018learning, ostafew2014learning] uses Gaussian Process (GP) regression for real-time learning of an uncertain model and satisfies probabilistic state constraints with a traditional stochastic MPC [farina2016stochastic] controller. Although such state-space modeling based Adaptive MPC controllers have proven to be effective, they involve construction of positive invariant sets [kolmanovsky1998theory, blanchini1999set], which can become computationally cumbersome. As a consequence, input-output modeling of systems has been opted in literature for proposing computationally efficient Adaptive MPC algorithms for certain applications (e.g. for stable, slow systems).

measured in terms of closed-loop cost, owing to the allowance of output constraint violations with a certain (low) probability.

In this paper, we build on the work of [tanaskovic2014adaptive, bujarbaruahAdapFIR], and propose an Adaptive Stochastic MPC algorithm that considers probabilistic output constraints and hard input constraints for a Multi Input Multi Output (MIMO) system. Similar to [bujarbaruahAdapFIR]

, we consider an uncertain FIR model of the system that is subject to bounded disturbances with known mean and variance. The support for the set of all possible models, which we call the Feasible Parameter Set (FPS), is adapted at each timestep using a set membership based approach. In contrast to previous work

[bujarbaruahAdapFIR], we additionally consider that the impulse response parameters corresponding to each output are sparse. Such sparse impulse response modeling can be motivated by [benesty2006adaptive, etter1985identification] for MIMO systems. Our goal is to utilize this additional sparsity information to demonstrate performance improvement over the algorithm in [bujarbaruahAdapFIR]. Our contributions can be summarized as follows:

• Offline before the control process, we compute a set containing all possible values of the unknown sparse FIR vectors corresponding to each output, with a very

high probability. This set, which we call the Feasible Sparse Parameter Set (FSPS) is computed using the Basis Pursuit Denoising [chen2001atomic] problem.

• Online during the control process, we obtain a point estimate of the unknown system inside the intersection of the FPS and the FSPS, using a Recursive Least Squares (RLS) estimator. Using this estimated system, we propagate our nominal predicted outputs used in the MPC controller objective function to improve performance. Simultaneously, we ensure satisfaction of the output chance constraints for the unknown true system.

• Through numerical simulations, we demonstrate that our algorithm exhibits better performance than the algorithm presented in [bujarbaruahAdapFIR].

## 2 Problem Description

### 2.1 System Modeling and Control Objective

We consider stable linear time-invariant systems described by a Finite Impulse Response (FIR) model of the form

 y(t)=HaΦ(t)+w(t), (1)

where the number of inputs and outputs considered is and , respectively. The FIR regressor vector of length is denoted by , where denotes the input at time . The matrix is a matrix comprising of the impulse response coefficients that relate inputs to the outputs of the system. The disturbance vector

is assumed to be a zero-mean random variable with a known variance, component-wise bounded as

 |wj(t)|≤¯wj,∀j=1,2,…ny, (2)

where the are assumed known. Finally, is the measured output of the system.

Our goal is to control the output while satisfying input and output constraints of the form

 Cu(t) ≤g,t=0,1,…, (3a) P{Ey(t)≤p} ≥1−ϵ,t=0,1,…, (3b)

where is the maximum allowed probability of output constraint violation. Following [bujarbaruahAdapFIR], we consider a single linear output chance constraint, meaning is a row vector and . Notice that joint (linear) chance constraints (3b) can be reformulated into a set of individual (linear) chance constraints using Bonferroni’s inequality, and can therefore also be addressed by our proposed framework.

###### Assumption 1.

We assume that each row of impulse response matrix is sparse. Without loss of generality, we further assume that the sparsity index for each row, i.e. the number of nonzero entries, is at most .

### 2.2 Method Outline

We assume in this paper that the system matrix in (1) is unknown. This paper proposes a method for identifying this unknown system matrix , and using the estimate in a robust control formulation to safely regulate the constrained uncertain system. Our proposed method uses the following steps:

1. Offline before the control process begins: Use number of collected input sequence regressors to compute a set , called the Feasible Sparse Parameter Set (FSPS). We compute this set via the Basis Pursuit Denoising problem [chen2001atomic]. The FSPS contains the true unknown model , with a high probability.

2. Online during the control process: At each timestep ,

1. Obtain the current output measurement and, using the known disturbance bounds (2), update the time-varying set , which we call the Feasible Parameter Set (FPS). The FPS is a set guaranteed to contain the true model .

2. Use the previous applied control inputs and measured outputs to construct an estimate of , lying in the intersection of and the offline-computed FSPS. The estimate is constructed using the Recursive Least Squares method.

3. Compute the input sequence that minimizes the objective function obtained with while satisfying the input and output constraints (3) for all models in the FPS . Apply the first computed control input and continue to step 2a.

In the following Section 3, we discuss steps 1-2(b). Step 2(c) is detailed in Section 4.

## 3 Model Estimation and Adaptation

We approximate system (1) with the form

 y(t)=H(t)Φ(t)+w(t), (4)

where our model is a random variable whose support we estimate online during the control process from the output measurements. The support for the set of all possible models consistent with the recorded system data, which we call the Feasible Parameter Set (FPS), is guaranteed to contain the true model . Based on the knowledge that system (1) has sparse impulse response properties, we also construct a Feasible Sparse Parameter Set (FSPS) offline by solving the Basis Pursuit Denoising Problem. During control run-time, a point estimate of is then computed to lie in the intersection of the offline-computed FSPS and the online-updated FPS. This estimate is then used in the control design. We decouple the offline and online phases of this design process and delineate the steps in detail in the following sections.

### 3.1 Offline: Construct the Feasible Sparse Parameter Set

The Feasible Sparse Parameter Set (FSPS), denoted by , is a function of the (unknown) true system response , and is synthesized offline utilizing the sparsity aspect of system responses from Assumption 1. Proposition 1 clarifies how this set is synthesized. We first introduce the following definition.

###### Definition 1 (Restricted Isometry Property (RIP) [candesTao]).

A matrix satisfies the Restricted Isometry Property (RIP) of order , with constant , if

 (1−δ)∥x∥22≤∥Ax∥22≤(1+δ)∥x∥22, ∀x ¯k sparse.

The order- restricted isometry constant is the smallest number such that the above inequality holds.

###### Proposition 1.

Suppose we collect output measurements offline. Suppose for , where each and , with , and denotes the row of . If , then any solution to the Basis Pursuit Denoising optimization problem

 min∥x∥1s.t.∥Ax−yi∥2≤√q¯wi, (denoted as B(Hai)) (5)

satisfies for some constant , and for all .

###### Proof.

See [YiMaBook, Theorem 3.5.1]. The proof follows from the known result of [candes2006stable, Theorem 1.1], with relaxation of to . ∎

The set is obtained as . Note that this set is synthesized offline, as the regressor vectors

are required to come from a Gaussian distribution in order to ensure the RIP property of each matrix

for . Such Gaussian inputs are not always allowable during control with system constraints (3). Details on our choice of these offline regressors to ensure with high probability are described in the Appendix.

###### Remark 1.

As pointed out in the proof of [YiMaBook, Theorem 3.5.1], a possible choice of the numerical constant in Proposition 1 is , with . However, since is not exactly known, computing and hence accurately is not possible. We see that as . Therefore offline regressor vectors should be chosen to ensure , with . Under such choice of the offline regressors as shown in the Appendix, we pick the constant .

### 3.2 Online: Update the Feasible Parameter Set

Following [tanaskovic2014adaptive, bujarbaruahAdapFIR], a set-membership identification method is used for updating the time-varying FPS, denoted by . The initialization of is done considering the fact that the true system (1) is stable. We update the FPS as given by

 F(t)=F(t−1)∩{H(t):H(t)Φ(t)≤y(t)+¯w}∩{H(t):−H(t)Φ(t)≤−y(t)+¯w}, (6)

where is the bound of the additive disturbance given by (2). In order to bound the computational complexity of (6) over time, an alternative algorithm to compute (6) is presented in [tanaskovic2014adaptive].

### 3.3 Online: Obtain Point Estimate μa(t)

We rewrite (4) as

 y(t) =Φ(t)H(t)+w(t),

where and are shown in the Appendix. Furthermore, let be the variance of the disturbance . Let the initial prior mean and variance estimates for true system be and respectively. Now, the conditional mean and variance estimates, given measurements up to , can be obtained using the Recursive Least Squares method [anderson1979, Sec. (3.1)]. These estimates may differ from the true conditional distribution parameters, as is not necessarily assumed to be Gaussian.

We ensure that the mean point estimate at any time instant is chosen as a point contained in a set , that is, , with

 Fp(t)=F(t)∩B(Ha). (7)

One way of obtaining such a constrained estimate of in is to project the mean. As shown in [bujarbaruahAdapFIR], this is achieved by solving

 μa(t)=argminX∈Fp(t)(X−μa(t))⊤M(X−μa(t)), (8)

where . The mean in matrix form, that is, is obtained by reorganizing into columns. This provides the minimum mean squared error estimate (with any linear estimator) of the true system . Note that (8) is a convex optimization problem, as the set in (7) is an intersection of two convex sets. We can further obtain a polytopic by outer approximation of with the infinity norm.

###### Remark 2.

In [bujarbaruahAdapFIR] the set in (7) is set as the FPS for all timesteps . Thus [bujarbaruahAdapFIR] does not utilize the sparsity information of rows of to construct the nominal point model estimate . The construction of in (7) with our proposed approach utilizes such sparsity information from Assumption 1. However, this comes with additional offline computation of , as shown in Proposition 1.

## 4 Control Synthesis

### 4.1 Prediction Model

Let be the prediction horizon for a predictive controller for system (1). We denote the predicted system outputs at time by , for some . Similarly, denotes the predicted regressor vector, for , and is computed as

 Φ(k|t)=WΦ(k−1|t)+Zu(k−1|t), (9)

where, the matrices and are as reported in the Appendix. With these predicted regressor vectors, the estimated system obtained in Section 3.3 is used to propagate the nominal predicted outputs as , for all . This is shown in the optimization problem presented in Section 4.3.

### 4.2 Reformulation of Chance Constraints

Within each prediction horizon we enforce , where is a function of some . Therefore, to ensure satisfaction of (3b) despite uncertainty in the true system, we must satisfy the constraint for all . Using the theory of distributionally robust optimization [calafiore2006distributionally, zymler2013distributionally], we can conservatively approximate the output chance constraints (3b) as

 κϵ√¯Φ⊤(k|t)Γ¯Φ(k|t)+Φ⊤(k|t)¯EH(t)−p≤0,  ∀H(t)∈F(t), (10)

where we have , and . Here, is an appended covariance matrix shown in the Appendix. As is a polytope, (10) can be succinctly written as

 κϵ√¯Φ⊤(k|t)Γ¯Φ(k|t)+Φ⊤(k|t)¯Efj(t)−p≤0, (11)

where denote the vertices of the polytope .

###### Remark 3.

The reformulated chance constraints (11) are robustly satisfied for all , and only the choice of the point model estimate depends on the set , as shown in (7). Hence, satisfaction of (3b) is guaranteed by (11), despite .

### 4.3 MPC Problem

We solve the following optimization problem for given weight matrices :

 minU(t)t+N−1∑k=t[^y⊤(k|t)Q^y(k|t)+u⊤(k|t)Su(k|t)]+^y⊤(t+N|t)Q^y(t+N|t)s.t.^y(k+1|t)=μa(t)Φ(k+1|t),^y(t|t)=y(t),Cu(k|t)≤g,Φ(t+N|t)=WΦ(t+N|t)+Zu(t+N−1|t),κϵ√¯Φ⊤(k+1|t)Γ¯Φ(k+1|t)+Φ⊤(k+1|t)¯Efj(t)≤p,∀k=t,…,t+N−1,∀fj(t)∈vertex(F(t)), μa(t)∈F(t)∩B(Ha), (12)

where , and the regressor is as in (9). We have included the terminal constraint on the regressor vector as given in [tanaskovic2014adaptive]

 Φ(t+N|t)=WΦ(t+N|t)+Zu(t+N−1|t). (13)

This means the terminal regressor corresponds to a steady state (i.e. the last control inputs in a horizon are kept constant). Problem (12) is a convex optimization problem and can be solved with existing solvers. After solving (12), we apply the first input

 u(t)=u⋆(t|t) (14)

to system (1) in closed-loop. We then re-solve (12) at timestep with new estimated data and . This yields a receding-horizon control scheme. The resulting algorithm is summarized in Algorithm 1.

###### Theorem 1.

Consider Algorithm 1 and the receding horizon closed-loop control law (14) applied to system (1) after solving optimization problem (12). If the optimization problem (12) is feasible at timestep , then it is feasible at all subsequent timesteps .

###### Proof.

Let (12) be feasible at timestep and let the corresponding optimal input sequence be . After applying control (14) in closed-loop to (1), consider a candidate open loop control sequence at the next timestep as

 U(t+1)=[u⋆(t+1|t)⊤,…,u⋆(t+N−1|t)⊤,u⋆(t+N−1|t)⊤]⊤. (15)

This candidate sequence (15) satisfies the input constraints (3a) and the terminal constraint (13) at . So, using candidate sequence (15) and condition (13), we obtain

 Φ(k|t+1)=Φ⋆(k|t),∀k∈{t+2,…,t+N}, (16a) Φ(t+N+1|t+1)=Φ⋆(t+N|t). (16b)

To show feasibility of (12) at , we finally must ensure that (11) satisfied with (15) at . That is, we require for all

 κϵ√¯Φ⊤(k|t+1)Γ¯Φ(k|t+1)+Φ⊤(k|t+1)¯Efj(t+1)−p≤0, (17)

for all vertices. Now, for the chosen input sequence (15), by using (13) and (16), condition (17) can be expressed as

 κϵ√¯Φ⋆⊤(k|t)Γ¯Φ⋆(k|t)+Φ⋆⊤(k|t)¯Efj(t+1)−p≤0, (18)

for all . We can now guarantee that (18) will be satisfied at timestep with (15). This is due to the observation that feasible parameter set follows as new cuts (6) are introduced at each timestep. So vertices at timestep are convex combinations of the ones at . Thus, control sequence (15) is feasible for (12) at timestep . This completes the proof. ∎

## 5 Numerical Simulations

We present simulation results for a simple single-input, single-output system111We choose this model for the sake of simulation simplicity. Sparsity in FIR models is well motivated typically for MIMO systems.. We compare the performance of our Algorithm 1 with that from the adaptive stochastic MPC presented in [bujarbaruahAdapFIR]. This performance is measured in terms of the expected closed-loop cost , where the closed-loop cost of any trajectory which is a function of the disturbance realization , is given by

 V(¯w,Φ(0),Fp(0),μa(0),σa(0))=tend∑t=0y⊤(t)Qy(t)+(u⋆(t))⊤Su⋆(t),

where for Algorithm 1 is obtained as in (7), and for the algorithm in [bujarbaruahAdapFIR], , as Remark 2 points out. For simulating both the algorithms, we use the parameters given in Table 1, with the true system response given as , which is sparse.

We run Monte Carlo simulations with both algorithms for randomly chosen disturbance sequences . We approximate the average closed-loop cost with the empirical average

 ^V(Φ(0),Fp(0),μa(0),σa(0))=1100100∑~m=1V((¯w)⋆~m,Φ(0),Fp(0),μa(0),σa(0)), (19)

where represents the Monte Carlo sample.

### 5.1 Cost Comparison

Fig. 1 shows the comparison of closed-loop cost expressed as for different Monte Carlo draws of trajectories, obtained with Algorithm 1 and the algorithm in [bujarbaruahAdapFIR]. We see that the empirical average closed-loop cost obtained as (19) for Algorithm 1 is around lower than the corresponding value obtained with [bujarbaruahAdapFIR].

This demonstrates performance gain by Algorithm 1 as a consequence of leveraging sparsity information of via the FSPS set . Furthermore, the closed-loop costs obtained for out of the trajectories were lower with Algorithm 1 compared to the corresponding cost obtained with [bujarbaruahAdapFIR]. This improvement is cost is explained in the next section via a closer look at the corresponding trajectories of the system.

### 5.2 Closed-loop Trajectory Comparison

In order to further analyze the performance gain seen above, we consider the difference in absolute value of the outputs along each one of 100 trajectories (Fig. 2). Specifically, we analyze

 Δ|y(t)|=|y1(t)|−|y2(t)|,

for , where and denote the outputs obtained using Algorithm 1 and the algorithm in [bujarbaruahAdapFIR] respectively. We see from Fig. 2 that out of all timesteps simulated, only at about of them is . This implies that at about of all timesteps simulated, Algorithm 1 resulted in outputs closer to the origin, as desired.

Fig. 3 further plots the difference in absolute value of corresponding closed-loop inputs, i.e.,

 Δ|u⋆(t)|=|u⋆,1(t)|−|u⋆,2(t)|,

for along all the 100 trajectories, where and denote the optimal closed-loop inputs obtained using Algorithm 1 and the algorithm in [bujarbaruahAdapFIR] respectively. At only of the 2000 timesteps , which indicates that the improved output regulation in Fig. 2 does not come at the cost of higher magnitudes of inputs. Thus, Fig. 2 and Fig. 3 combined, provides a comprehensive explanation of cost improvement in Fig. 1.

## 6 Conclusions

We proposed an Adaptive Stochastic Model Predictive Control (MPC) algorithm for stable linear time-invariant MIMO systems expressed with an FIR model. The parameters of the FIR model corresponding to each output are sparse, but unknown. We solve a Basis Pursuit Denoising problem offline before the control process, which utilizes the sparsity information of the system model to give set based bounds containing the true system parameters. Online during control process, we estimate these parameters with the Recursive Least Squares algorithm and the estimates are refined using the set based bounds obtained offline. With Set Membership Methods, our MPC controller safely ensures satisfaction of all imposed constraints by the unknown true system. With a thorough numerical example we demonstrated performance gain over the algorithm of [bujarbaruahAdapFIR].

## Appendix

### Choosing Offline Regressors

Following [YiMaBook, Theorem 3.3.4], there exists a numerical constant such that if

is a random matrix with entries independent

random variables, with high probability, , provided

 q≥2~C¯klog(num2¯k)¯δ2. (20)

Offline regressor vectors are thus chosen satisfying (20), with each entry of as . The constant

can be exactly found using properties of phase transition in compressed sensing

[donoho2006compressed].

### Matrix Notations

We define

 H(t)=[H1,H2,…,Hny]⊤∈Rnynum×1, Φ(t)=diag(Φ⊤(t),Φ⊤(t),…,Φ⊤(t))∈Rny×nynum,

where denotes the row of . Moreover,

 ¯q =⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣00⋯0010⋯0001⋯00⋮⋮⋱⋮⋮00⋯10⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦∈Rm×m

Based on this matrix we get:

 W=diag(¯q,¯q,…,¯q)∈Rnum×num, z=[1,0,…,0]⊤∈Rm,

which gives .

### Chance Constraint to Conic Optimization

For this part we repeat the derivation of [bujarbaruahAdapFIR]. To ensure satisfaction of (3b) with unknown true system , we must satisfy them for all . Thus is treated as a deterministic variable (mean and variance ) while reformulating (3b) to equivalent convex constraints. Therefore, for all , we have:

 P{EH(t)Φ(k|t)+Ew(k)≤p} ≥1−ϵ ⟺P{[EH(t)Ew(k)][Φ⊤(k|t)1]⊤≤p} ≥1−ϵ.

We denote,

 a⊤1(t) =[EH(t)Ew(k)],⟹^a⊤1(t)=[EH(t)0],

where denotes the mean of a quantity . Moreover,

 ¯Φ(k|t) =[Φ⊤(k|t)11]⊤, and,
 d1(t) =[a⊤1(t)−p]⊤,⟹^d1(t)=[^a⊤1(t)−p]⊤.

From these we can derive the variance of as

 Γ =σ2(d1(t))=σ2(⎡⎢⎣¯EH(t)Ew(k)−p⎤⎥⎦), where ¯E=diag(E,E,…,E), =diag(0,Eσ2wE⊤,0)⪰0.

As in [bujarbaruahAdapFIR], we assume no correlation between the disturbance and the impulse response distribution. Now, (3b) is reformulated into second-order cone constraints following [calafiore2006distributionally]

 κϵ√{¯Φ⊤(k|t)Γ¯Φ(k|t)}+^d⊤1(t)¯Φ(k|t)≤0, (21)

for all , where for any bounded disturbance distributions

with known moments. After simplifications, (

21) can be written for all as

 κϵ√{¯Φ⊤(k|t)Γ¯Φ(k|t)}+Φ⊤(k|t)¯EH(t)−p≤0.