Adaptive Learning-based Model Predictive Control for Uncertain Interconnected Systems: A Set Membership Identification Approach
Ahmed Aboudonia John Lygeros





 

[footnoteinfo1]Funded by NCCR Automation, the Swiss National Science Foundation (Corresponding Author: Ahmed Aboudonia)


Zurich]Ahmed Aboudoniaahmedab@ethz.ch,    
Zurich]John Lygeroslygeros@ethz.ch               


[Zurich]Automatic Control Laboratory, ETH Zurich  

          
                           
Model predictive Control; Optimization-based Control; Learning-based Control.   
  



  
We propose a novel adaptive learning-based model predictive control (MPC) scheme for interconnected systems which can be decomposed into several smaller dynamically coupled subsystems with uncertain coupling. The proposed scheme is mainly divided into two main online phases; a learning phase and an adaptation phase. Set membership identification is used in the learning phase to learn an uncertainty set that contains the coupling strength using online data. In the adaptation phase, rigid tube-based robust MPC is used to compute the optimal predicted states and inputs. Besides computing the optimal trajectories, the MPC ingredients are adapted in the adaptation phase taking the learnt uncertainty set into account. These MPC ingredients include the prestabilizing controller, the rigid tube, the tightened constraints and the terminal ingredients. The recursive feasibility of the proposed scheme as well as the stability of the corresponding closed-loop system are discussed. The developed scheme is compared in simulations to existing schemes including robust, adaptive and learning-based MPC.






§ INTRODUCTION


Model Predictive Control (MPC) is an advanced control scheme that optimally controls a given system by repeatedly solving an online optimal control problem (OCP) (<cit.>). MPC was successfully used in various applications such as interconnected systems (e.g. see <cit.> and the references therein).
MPC, however, heavily relies on the availability of an accurate model which is, in many cases, difficult to find and/or complex to use. 
Various efforts were devoted to finding MPC variants that do not require an accurate model. 
These variants include robust, stochastic and learning-based MPC (<cit.>). 
One drawback of these schemes is that the model is usually not updated online leading to potentially conservative performance. 
Recently, various adaptive and safe learning-based MPC schemes (<cit.>) have been proposed to ensure safety while enhancing closed-loop performance by improving the model using data collected online. This goal is achieved by combining physics-based knowledge with machine learning and system identification.

In this paper, we propose a novel adaptive learning-based MPC scheme for interconnected systems which can be decomposed into several smaller subsystems. 
Each subsystem is dynamically coupled with a set of other subsystems called the set of neighbours. 
The source of uncertainty in each subsystem is the strength of the coupling with its neighbours.
The proposed scheme learns the uncertainty set online by applying distributed set membership identification on collected data and adapts the MPC ingredients accordingly. 
In particular, the developed control scheme adapts the rigid tube, the prestabilizing controller, the tightened constraint sets and the terminal ingredients by considering them as decision variables in the OCP. 
The developed scheme is compared in simulation to existing schemes including robust MPC, adaptive MPC and learning-based MPC using a network of dynamically coupled double integrators and found to lead to a better trade-off between closed-loop performance and computation cost.

The main contributions and challenges of this work can be listed as follows:



    
  * Unlike many MPC schemes which usually only compute optimal state and input trajectories online, we formulate an MPC problem which can also compute various MPC ingredients in the adaptation phase. The main challenge here is to maintain the convexity of the OCP since these ingredients should generally satisfy non-convex infinite dimensional conditions.
    
  * Although set membership identification is extensively used for monolithic systems, we extend this approach here to interconnected systems with uncertain couplings where parameter estimation is not handled by a centralized unit.
    
  * We show that the closed-loop system under the developed controller is input-to-state stable. The main challenges here are that the closed-loop system is switched due to updating the MPC ingredients online and that artificial equilibrium points are used. Hence, standard MPC stability arguments do not apply.
    
  * We compare the proposed scheme to other off-the-shelf MPC schemes which prove to be successful in many applications.












In Section 2, we formulate the considered control problem. In Section 3, we develop the novel adaptive learning-based MPC scheme. In Section 4, we prove the recursive feasibility of the proposed scheme and the input-to-state stability of the corresponding closed-loop system. In Section 5, we validate the efficacy of the developed scheme by comparing it to other existing schemes in the literature. Finally, we draw some conclusions in Section 6. 
In the sequel, we refer to symmetric, symmetric positive semidefinite and symmetric positive definite matrices of size n as 𝕊^n, 𝕊^n_+ and 𝕊^n_++, respectively.
Similarly, we denote diagonal, diagonal positive semidefinite and diagonal positive definite matrices of size n by 𝔻^n, 𝔻^n_+ and 𝔻^n_++,  respectively. Additionally, we use ⊕ and  to refer to the Minkowski set addition and the Pontryagin set difference, respectively.



§ PROBLEM FORMULATION

 
We consider discrete-time linear time-invariant interconnected systems with a sparse structure.
We assume that such systems can be decomposed into M dynamically-coupled subsystems and that the i-th subsystem has a set of in-neighbours 𝒩_i^- including all subsystems which affect its dynamics and a set of out-neighbours 𝒩_i^+ including all subsystems whose dynamics are affected by the i-th subsystem.
We assume that the dynamics of the i-th subsystem is given by,



    z_i(t+1)  
        = 
         A_iz_i(t)+B_iv_i(t) + E_i ∑_j ∈𝒩_i^- a_ij C_j z_j(t) + d_i(t),

where z_i ∈ℝ^n_i, v_i ∈ℝ^p_i and d_i ∈ℝ^n_i are the state, input and noise vectors of the i-th subsystem, A_i∈ℝ^n_i × n_i, B_i ∈ℝ^n_i × p_i, C_i∈ℝ^m_0 × n_i and E_i ∈ℝ^n_i × m_0 are known matrices and a_ij∈ℝ are unknown fixed parameters that represent the coupling strengths. 
We assume that the scalars a_ij satisfy 0 ≤ a_ij^min≤ a_ij≤ a_ij^max where the bounds a_ij^min and a_ij^max are known. 
We also assume that each subsystem is subject to the state and input constraints z_i(t) ∈𝒵_i = {z_i : G_i z_i ≤ 1} and v_i(t) ∈𝒱_i = { v_i : H_i v_i ≤ 1 } where G_i ∈ℝ^q_i × n_i, H_i ∈ℝ^r_i × p_i, 𝒵_i and 𝒱_i are compact and include the origin in their interior. The noise vector d_i(t) is assumed to lie in a known compact polytopic set 𝒟_i = {d_i : Ξ_i d_i ≤ 1 } that includes the origin in its interior.
Following <cit.>, we express (<ref>) as,



    z_i(t+1) = A_iz_i(t)+B_iv_i(t) + E_iw_i(t) + d_i(t),




    w_i(t) = ∑_j ∈𝒩_i^- a_ijC_jz_j(t) ∈𝒲_i = ⊕_j ∈𝒩_i^- a_ij^max C_j Z_j.







By denoting the nominal state and input vectors of the i-th subsystem by x_i ∈ℝ^n_i and u_i ∈ℝ^p_i, we also define the nominal dynamics of the i-th subsystem to be, x_i(t+1) = A_ix_i(t)+B_iu_i(t).
Assuming that v_i(t)=u_i(t)+K_i(z_i(t)-x_i(t)) where K_i ∈ℝ^p_i × n_i is a stabilizing control gain and defining e_i(t)=z_i(t)-x_i(t), we introduce the error dynamics,



    e_i(t+1) = (A_i+B_iK_i)e_i(t)+E_iw_i(t)+d_i(t).

We consider here a tracking problem where the state is required to follow the target point x_r_i which is an equilibrium point of the nominal system. The target point x_r_i and the corresponding input u_r_i are assumed to satisfy the state and input constraints. For this purpose, we define the artificial equilibrium point (x_e_i,u_e_i) satisfying x_e_i = A_i x_e_i + B_i u_e_i 
and the cost function
J_i = ∑_k=0^N-1( x_i(k|t)-x_e_i_Q_i^2 + u_i(k|t)-u_e_i_R_i^2 ) + x_e_i-x_r_i_S_i^2
where x_i(k|t) and u_i(k|t) are the k-step ahead predicted nominal state and input vectors of the i-th subsystem at the t-th time instant. 
Note that (x_e_i,u_e_i) is used to enlarge the feasible region in tracking problems (<cit.>) and is considered as an artificial reference that converges asymptotically to the target point.

Although the main aim is to drive the actual state to the target point, we use the nominal state in J_i since we can only predict the behavior of the nominal dynamics.
Note, however, that the closed-loop stability guarantees derived in Section 4 apply to the actual states.
In summary, an MPC problem can be formulated for the i-th subsystem with the local OCP,



    min_ 
            J_i   
            
            
            s.t. {    for all  k ∈{0,...,N-1}, 
        z_i(t) - x_i(0|t) ∈𝒵_K_i , 
        x_i(k+1|t)=A_ix_i(k|t)+B_iu_i(k|t), 
        x_i(k|t) ∈𝒳_i,    u_i(k|t) ∈𝒰_i, 
        x_i(N|t) = x_e_i=A_ix_e_i+B_iu_e_i, 
        x_e_i∈ξ𝒳_i,    u_e_i∈ξ𝒰_i,
            .

where x_i(k|t)|_k ∈{0,,N}, u_i(k|t)|_k ∈{0,,N-1}, x_e_i and u_e_i are the decision variables, z_i(t) is the measured state of the i-th subsystem, 𝒵_K_i is a robust positive invariant (RPI) set for its error dynamics (<ref>) under the disturbances w_i(t) and d_i(t), ξ∈ (0,1), 𝒳_i = {x_i : G_i x_i ≤ g_i} and 𝒰_i = { u_i : H_i x_i ≤ h_i} are its tightened constraint sets with g_i and h_i positive scalars. The last three constraints in (<ref>) indicate that x_i(N|t) reaches an invariant terminal set satisfying the state and input constraints. 
The scalar ξ is used to ensure that the artificial equilibrium lies in the interior of the constraint sets. This requirement is needed later in the stability proof.
Note that the global prestabilzing control gain is denoted by K=diag(K_1,...,K_M). Similarly, the global system matrices A and B can be obtained in the obvious way. Following <cit.>, the control input v_i(t)=u_i(0|t)+K_i(z_i(t)-x_i(0|t)) is applied to (<ref>) to ensure robust constraint satisfaction.

Interconnected systems with uncertain coupling (<ref>) can be exploited in various scenarios with varying-topology networks and large-scale systems. For example, varying-topology networks comprise subsystems which can join and leave the network. While accurate models can be identified offline for the new subsystems, the coupling strength between the network and these subsystems can only be identified online after the plug-in operation. 
Practical applications comprising such systems include microgrids, power systems and vehicle platoons.


Note that considering the whole coupling term as a disturbance might be conservative. This, however, leads to decentralized OCPs rather than distributed ones. Decentralized OCPs are usually more robust against failures and requires less communication and computation. This is more appealing for safety-critical systems with fast dynamics.



Moreover, the inclusion of parametric uncertainties within an additive term might be also conservative. An alternative consideration involves employing adaptive MPC tailored for parametric uncertainties, as applied in monolithic systems (<cit.>). However, the complexity arises in interconnected systems when uncertainties are embedded within coupling terms, necessitating the formulation of a distributed OCP, which imposes substantial communication and computational burdens.
Deceralized OCPs have been widely used with interconnected systems (<cit.>).
Finally, note that coupling constraints among neighbours could in principle be incorporated in decentralized OCPs by finding a decentralized inner-approximation of the constraint sets. This, however, adds to the conservatism of decentralized OCPs.

In this paper, we make use of the OCP (<ref>) to develop the proposed scheme which, at each time instant, performs two phases: a learning phase to learn the model uncertainty and an adaptation phase to compute the optimal control action. We assume that every subsystem can communicate with its neighbours at each time instant.  
While this assumption is not required for the adaptation phase which relies on the decentralized OCP (<ref>) that requires no communication, it is still required for the learning phase. This is because the model uncertainty in (<ref>) relies on the states of the neighbours. Hence, communication is required among neighbours only once at each time instant in the learning phase. Since communication is allowed by assumption, the decentralized OCP (<ref>) can in principle be replaced by a distributed OCP. While this may alleviate the conservatism of the developed scheme, the computation and communication complexity will increase significantly since distributed optimization is required.




In the learning phase, we aim to reduce the volume of the uncertainty set 𝒲_i in (<ref>) by learning the bounds a_ij^min and a_ij^max of the uncertain parameter a_ij as more data is collected online. Hence, we develop a distributed set membership identification approach which learns these bounds by sharing information among the neighbours only once at each time instant. 
While the set 𝒲_i is a function of a_ij^max only, a_ij^min is also updated since it can speed up the shrinkage of  𝒲_i by finding smaller upper bounds.
In the adaptation phase, we aim to adapt all MPC ingredients depending on the learnt uncertainty set and compute the local optimal control inputs taking the updated uncertainty set into consideration.

The adapted MPC ingredients include the rigid tube represented by 𝒵_K_i, the prestabilizing controller K_i, the tightened state constraint set 𝒳_i and the tightened input constraint set 𝒰_i. 
Note that these ingredients are updated at each time instant t based on the learnt uncertainty set 𝒲_i which is a function of the updated parameters a_ij^max and a_ij^min at this time instant. Hence, these MPC ingredients together with the uncertainty set are now functions of time.


In the sequel, we use the ellipsoidal rigid tube 𝒵_K_i(t)={e_i : e_i^⊤ Z_i e_i ≤α_i^2(t)} and the linear prestabilizing controller κ_i(x_i)=K_i(t)x_i. We assume that K_i(t)=K_o_iT_i(t) where K_o_i is a full-rank stabilizing control gain and T_i(t) ∈𝕊^n_i_++. Recall also that the tightened polytopic state and control constraint sets are given, respectively, by 𝒳_i(t) = {x_i : G_i x_i ≤ g_i(t)} and 𝒰_i(t) = { u_i : H_i x_i ≤ h_i(t)}. 






In addition to the decision variables in (<ref>), α_i(t), T_i(t), g_i(t) and h_i(t) are also considered as decision variables in the developed scheme in order to update the MPC ingredients online; K_o_i, on the other hand, is assumed to be known in advance, computed, for example, with methods as those in <cit.>. 
When updating these ingredients at each time instant, one should ensure the following conditions,

    

    
  I. The set 𝒵_K_i(t) is an RPI set for the dynamics (<ref>).
    
  II. The local matrices A_i+B_iK_i(t) are Schur. 
    
  III. The global matrix A+BK(t) is Schur.
    
  IV. The condition 𝒳_i(t) ⊕𝒵_K_i(t) ⊆𝒵_i holds. 
    
  V. The condition 𝒰_i(t) ⊕ K_i(t)𝒵_k_i(t) ⊆𝒱_i holds. 
    


The aforementioned conditions are ensured by adding extra constraints in the OCP of the adaptation phase while preserving its convexity. These additional constraints are derived in the next section. 
Although the constraint sets are polytopic, the rigid-tube is chosen to be ellipsoidal. As shown in the next section, this leads in the adaptation phase to an OCP which considers the prestabilizing controller as a decision variable while maintaining the convexity of the OCP by making use of the S-lemma (<cit.>). Recall also the scalability properties of ellipsoidal sets in previous related works when controlling networked dynamical systems (<cit.>).






The system (<ref>) can be described by a graph 𝒢(ℳ,ℰ,𝒜) where ℳ = {1, ..., M} represents the set of nodes, ℰ⊆ (ℳ×ℳ) the set of edges and 𝒜 = {a_ij∈ℝ ,(i, j) ∈ℰ} the set of weights. Each subsystem is represented by a node. An edge exists from the j-th node to the i-th node if the j-th subsystem belongs to the set of in-neighbours of the i-th subsystem. We define the adjacency matrix A_g ∈ℝ^M × M of the graph 𝒢 such that the element in its i-th row and j-th column is given by A_g_ij=a_ij if i ∈𝒩_j^+ and zero otherwise. 









The degree matrix D_g ∈𝔻_++^M of the graph 𝒢 is defined such that the diagonal element in the i-th row is given by ∑_j ∈𝒩_i^- a_ij and the Laplacian matrix is given by L_g = D_g - A_g  ∈𝕊_+^M. 

These graph theory concepts are required to locally ensure Condition III that the global matrix A+BK(t) is Schur. The global dynamics is then given by,

    z(t+1)=Az(t)+Bv(t)+d(t),

where 
z(t)=[z_1^⊤(t),,z_M^⊤(t)]^⊤, 
v(t)=[v_1^⊤(t),, v_M^⊤(t)]^⊤, 
d(t)=[d_1^⊤(t),,d_M^⊤(t)]^⊤, 
A_d=diag(A_1, , A_M), 
Ã_g = A_g ⊗ I_m_0,
A=A_d+E Ã_g C, 
B=diag(B_1, ,B_M), 
C=diag(C_1,,C_M) and 
E=diag(E_1,, E_M).
Note that the proposed method is mainly developed for sparse graphs where only neighbours can communicate. While this method can be used with fully connected graphs where all subsystems can share information with each other, other efficient methods (e.g. <cit.>) can be used where all subsystems share information with a central entity that does all the computations instead.
In the sequel, we make use of the following definitions.



    
    A set Z_K⊂ℝ^n is robustly positively invariant under z^+=Az+Bv, v=Kz  and z ∈𝒵, v ∈𝒱, if and only if z^+ ∈𝒵_k, z ∈𝒵, v ∈𝒱 for all z ∈𝒵_k.




    
    The system z^+=Az+Bv, y=Cz+Dv is strictly passive if there exist functions V_p:𝒵→ℝ and γ_p:𝒵→ℝ_+ so that V_p(0)=0 and V_p(z^+)-V_p(z) ≤ y^⊤ v - γ_p(z) for all v ∈𝒱.











§ ADAPTIVE LEARNING-BASED MPC


In this section, we start by discussing the distributed set membership identification approach used to learn the uncertainty set online in the learning phase. We then derive the constraints added to the local online OCP of each subsystem in the adaptation phase.
First, we discuss the learning phase where we aim to learn the uncertainty set 𝒲_i(t) in real time, as more data is collected online. According to (<ref>), this set is a function of the bounds of the uncertain parameters. Hence, a better estimation of these bounds leads to a better description of this set. Thus, we develop here a distributed set-membership identification technique to improve the estimates of these bounds. Although such techniques are widely used with monolithic systems (e.g. see <cit.>), we extend this technique here to interconnected systems with uncertain couplings where the identification is not carried out by a central entity.

Assume that the bounds of the uncertain parameter a_ij at the (t-1)-th time instant are a_ij^min(t-1) and a_ij^max(t-1). Our goal is to improve the estimates of these bounds at the t-th time instant. For this purpose, we construct for every a_ij the feasible parameter set defined by the inequality F_ij a_ij≤ f_ij(t-1) where F_ij = [1  -1]^⊤ and f_ij(t-1) = [a_ij^max(t-1)  a_ij^min(t-1)]^⊤. We then define the vector θ_i which comprises all the uncertain parameters to be identified by the i-th subsystem (i.e. a_ij, j ∈𝒩_i^-). In this case, the feasible parameter set for the i-th subsystem is given by F_i θ_i ≤ f_i(t-1)



where F_i is a block diagonal matrix with the matrices F_ij for all j ∈𝒩_i^- on the diagonal and f_i is a vector vertically concatenating the subvectors f_ij for all j ∈𝒩_i^-.

Recall that each subsystem can communicate with its neighbours in the learning phase. Hence, the states of the neighbours of the i-th subsystem at the (t-1)-th time instant are available for the i-th subsystem at the t-th time instant. Furthermore, the states and inputs of the i-th subsystem at the (t-1)-th time instant are available. Finally, the states of the i-th subsystem at the t-th time instant can be measured. Recall also that 𝒟_i = {d_i : Ξ_i d_i ≤ 1}. Defining χ_i(t-1) to be a matrix horizontally concatenating the vectors C_j z_j(t-1) for all j ∈𝒩_i^-, Ω_i(t-1)=-Ξ_i E_i χ_i(t-1) and ω_i(t-1)=1 + Ξ_i A_i z_i(t-1) + Ξ_i B_i v_i(t-1) - Ξ_i z_i(t), we can construct for the i-th subsystem the non-falsified parameter set Δ_i = {θ_i: Ω_i(t-1) θ_i ≤ω_i(t-1) }.




























In order to update the bounds of the uncertain parameters, the i-th subsystem should solve 2|𝒩_i^-| optimization problems at the t-th time instant where |·| denotes set cardinality. Each optimization problem updates one of these bounds. For each j ∈𝒩_i^-, the lower and upper bounds of the parameter a_ij are given respectively by,





    a_ij^min(t) 
        = min_a_ij e_ij^⊤θ_i 
         s.t.[      F_i; Ω_i(t-1) ]θ_i ≤[ f_i(t-1); ω_i(t-1) ],






    a_ij^max(t) 
        = max_a_ij 
        e_ij^⊤θ_i 
         s.t.[      F_i; Ω_i(t-1) ]θ_i ≤[ f_i(t-1); ω_i(t-1) ],


where e_ij is a unit vector which selects the parameter a_ij in the vector θ_i. 

Note that computing the new lower (upper) bound is a minimization (maximization) problem because we aim to find the smallest (largest) value of the corresponding parameter that lies inside the feasible parameter set and the non-falsified set.
In conclusion, the i-th subsystem learns the uncertainty set 𝒲_i(t) in (<ref>) by solving at the t-th time instant the optimization problems (<ref>) which identify the bounds of the uncertain parameter. These bounds are communicated among the in-neighbours and out-neighbours to be used in the adaptation phase. 
We also highlight the dependency of the upper bound a_ij^max(t) on the lower bound a_ij^min(t-1) which is part of the vector f_i(t-1) in (<ref>) to support the fact that updating the lower bound might improve the updates of the upper bound.
Note that using asymmetric duality of linear programs, (<ref>) can be written compactly as a single optimization problem. We guarantee that the true parameters always lie within the updated bounds in the following theorem.

    Assume that a_ij^min(0) ≤ a_ij≤ a_ij^max(0). Then, a_ij^min(t) ≤ a_ij≤ a_ij^max(t) for all t ≥ 0 where a_ij^min(t) and a_ij^max(t) are computed using (<ref>).
    



    Note that z_i(1) = A_iz_i(0)+B_iv_i(0) + E_i ∑_j ∈𝒩_i^- a_ij C_j z_j(0) + d_i(0) and recall that d_i(0) ∈𝒟_i for all i ∈ℳ. 
    Thus, the true parameters satisfy Ξ z_i(1) - Ξ A_i z_i(0) + Ξ B_i v_i(0) + Ξ E_i ∑_j ∈𝒩_i^- a_ij C_j z_j(0)) ≤ 1 and hence, Ω_i(0) θ_i ≤ω(0) for all i ∈ℳ.
    Since a_ij^min(0) ≤ a_ij≤ a_ij^max(0), the true parameters also satisfy F_i θ_i ≤ f_i(0) for all i ∈ℳ by assumption.
    Since a_ij^min(1) and a_ij^max(1) are computed using (<ref>), then these are the smallest and largest values respectively satisfying Ω_i(0) θ_i ≤ω(0) and F_i θ_i ≤ f_i(0) that the true parameters also satisfy and thus, a_ij^min(1) ≤ a_ij≤ a_ij^max(1). 
    By induction, the same results hold for all t.
    □
    


As a consequence of the above theorem, it is guaranteed that a_ij^min(t) ≥ a_ij^min(t-1) because a_ij^min(t) is computed by minimizing over all possible values satisfying a_ij≥ a_ij^min(t-1). The same arguments apply to the upper bound and hence, 𝒲_i(t) ⊆𝒲_i(t-1) for all t (i.e. the volume of the uncertainty set can never increase).
Note that the proposed scheme is not necessarily restricted to set membership identification. While this technique is used here to update uncertain parameters in linear coupling terms, the proposed scheme paves the way to develop other schemes which can learn more complex uncertainties with nonlinear couplings and other disturbances. In this case, one can use other machine learning methods which have shown success in learning uncertainty sets.

We now move to discussing the adaptation phase.
As mentioned in Section <ref>, the online OCP is solved in the adaptation phase to adapt the MPC ingredients based on the learnt uncertainty set and compute the optimal control inputs. To adapt the MPC ingredients online while ensuring robust constraint satisfaction, the five conditions mentioned in Section <ref> should be satisfied. We derive here additional constraints which ensure the satisfaction of these conditions while maintaining the convexity of the OCP (<ref>).

We start with the condition which ensures the robust positive invariance of the set 𝒵_K_i(t). First, we define for the i-th subsystem the output vector y_i(t) = C_i z_i(t) and the constraint set 𝒴_i, obtained from 𝒵_i in the obvious way, such that y_i ∈𝒴_i. We then compute a minimum volume ellipsoidal set 𝒴̂_i = {y_i : y_i-c_y_i^2_Y_i≤ 1} that contains the set 𝒴_i. For this purpose, the vertex description of 𝒴_i is needed. Note that this is in general combinatorial in the linear inequalities needed to define 𝒴_i. However, the sets 𝒴_i of all subsystems do not change over time and hence, their minimum volume ellipsoidal sets 𝒴̂_i can be computed once offline by solving a convex optimization problem for each. Finally, we construct the outer approximation 𝒲̂_i(t) = ⊕_j ∈𝒩_i^- a_ij^max(t) 𝒴̂_̂ĵ of the uncertainty set 𝒲_i(t) = ⊕_j ∈𝒩_i^- a_ij^max(t) 𝒴_j. Note that  𝒲̂_i(t) is not necessarily an ellipsoidal set, but it is the Minkowski sum of the ellipsoidal sets 𝒲_ij(t) = a_ij^max(t)𝒴̂_j = {w_ij : w_ij-a_ij^max(t)c_y_j_Y_j^2 ≤ a_ij^max^2(t)}. 
Similarly, we also compute a minimum volume ellipsoidal set 𝒟̂_i = {d_i : d_i_Ξ̂_i^2 ≤ 1} offline for the set 𝒟_i of the i-th subsystem for all i ∈{1,,M}. 

In Proposition <ref>, we consider Condition I in Section <ref> and derive constraints which ensure the robust positive invariance of the set 𝒵_K_i(t) for the error dynamics (<ref>) under the disturbances w_i ∈𝒲̂_i(t) and d_i ∈𝒟̂_i. Although this might be more conservative than considering the disturbances w_i ∈𝒲_i(t) and d_i ∈𝒟_i, the convexity of the resulting OCP is preserved in this case. Another source of conservatism in Proposition <ref> is using the approximate S-lemma to reach local convex constraints. Recall here that K_i(t) = K_o_i T_i(t). The dependence on time in (<ref>) is omitted in the interest in space




    
    The set Z_K_i(t) is robustly positively invariant for the error dynamics (<ref>) with the uncertainty sets 𝒲̂_i(t) and 𝒟̂_̂î if there exist non-negative scalars λ_i(t), λ_d_i(t)  and λ_ij(t) for all j ∈𝒩_i^- such that (<ref>) holds.
    
    






    According to Definition <ref>, the set Z_K_i(t) is robustly positively invariant under the considered dynamics and uncertainty sets if and only if e_i^+_Z_i^2 ≤α_i^2(t) for all e_i satisfying e_i_Z_i^2 ≤α_i^2(t), w_ij satisfying w_ij-a_ij^max(t) c_y_j_Y_j^2 ≤ a_ij^max^2(t) for all j ∈𝒩_i^- and d_i satisfying d_i_Ξ̂_i^2 ≤ 1. 
    By defining s_i such that e_i = α_i(t) s_i, the invariance condition becomes (A_iα_i(t) + B_i K_o_i T_i(t) α_i(t)) s_i + E_i w_i + d_i_Z_i^2 ≤α_i^2(t) for all s_i satisfying s_i_Z_i^2 ≤ 1, w_ij satisfying w_ij-a_ij^max(t) c_y_j_Y_j^2 ≤ a_ij^max^2(t) for all j ∈𝒩_i^- and d_i satisfying d_i_Ξ̂_i^2 ≤ 1. 
    We now define the block diagonal matrix Y_λ_i(t) with the submatrices λ_ij(t) Y_j for all j ∈𝒩_i^- on the diagonal and the vector τ_i(t) which vertically concatenates the subvectors a_ij^max(t) c_y_j. 
    Using the S-lemma and applying the Schur complement yields (<ref>) where I_n_i is an identity matrix of size n_i and ℰ_i is a matrix horizontally concatenating the matrix E_i for |𝒩_i^-| times.
    
    □



























Note that the dimension of the constraint in (<ref>) increases linearly with respect to the number of neighbours. 
Rigid tubes are known for their little computation burden in comparison to homothetic tubes (<cit.>). Since we use time-varying rigid tubes, one might think that rigid tubes lose their main computational advantage and could be better replaced by homothetic tubes. This is, however, not the case because when using homothetic tubes, the nominal dynamics would be replaced by tube dynamics. Since the prestabilizing controller is a decision variable and the tube is ellipsoidal, the tube dynamics is represented by one matrix inequality for each prediction step. In this case, the problem becomes computationally challenging as the prediction horizon gets longer.

Next, we consider Conditions II and III in Section 2 by deriving constraints ensuring that A_i+B_iK_i(t) for all i ∈{1,,M} and A+BK(t) are Schur. 

For this purpose, we rely on passivity arguments for discrete-time interconnected systems. 
First, we derive a constraint ensuring that the noise-free dynamics of the i-th subsystem (i.e. (<ref>) with d_i(t)=0 for all t) under the controller u_i(t)=K_o_iT_i(t)x_i(t) is strictly passive. We show that this implies that A_i+B_iK_i(t) is Schur. We then derive other constraints ensuring that the strict passivity of all subsystems implies the asymptotic stability of the whole system ensuring that A+BK(t) is Schur. 
To ensure strict passivity of the i-th subsystem, we define the virtual output ỹ_i(t)=C_iz_i(t)+D_iw_i(t) where D_i ∈𝔻_++^m_0 is a decision variable. We also assume that a matrix P_i ∈𝕊_++^n_i is computed offline together with K_o_i, for example as in <cit.>, so that V_p_i(z_i)=z_i_P_i^2 is a storage function over the set 𝒵_i. We finally define a minimum volume ellipsoidal set 𝒵̂_i={z_i:z_i-c_z_i_Ẑ_i≤ 1} for the state constraint set 𝒵_i. The dependence on time in (<ref>) is omitted in the interest in space.


    
    The noise-free dynamics of the i-th subsystem z_i^+=(A_i+B_iK_o_iT_i(t))z_i + E_iw_i is strictly passive with respect to w_i ∈𝒲̂_i(t) and ỹ_i ∈ C_i𝒵̂_i ⊕ D_i 𝒲̂_i(t) if there exist Γ_i(t)  ∈𝔻_++^n_i, D_i(t)  ∈𝔻_++^n_i, T_i(t) ∈𝕊_++^n_i, α_i(t) > 0 such that,
    
    [
            [            (P_i - Γ_i) α_i               0.5 C_i^⊤α_i (A_iα_i+B_i K_o_iT_iα_i)^⊤;    [2pt/2pt]
                *                    D_i α_i                   E_i^⊤α_i;    [2pt/2pt]
                *                          *                  P_i^-1α_i ]]
            ≥ 0.

    

 
    To ensure strict passivity, we have to satisfy the condition in Definition <ref> for all w_i ∈𝒲̂_i(t). Recall also that the storage function V_p_i(z_i) is defined over the set 𝒵̂_i.
    Hence, the noise-free dynamics of the i-th subsystem under the controller u_i(t)=K_o_iT_i(t)x_i(t) is strictly passive according to Definition <ref> if z_i^+_P_i^2 - z_i_P_i^2 ≤ỹ_i^⊤ w_i - z_i_Γ_i^2 for all z_i satisfying z_i-c_z_i_Ẑ_i^2 ≤ 1 and w_ij satisfying w_ij-a_ij^max(t) c_y_j_Y_j^2 ≤ a_ij^max^2(t) for all j ∈𝒩_i^-. 
    By substituting for the dynamics, the strict passivity condition is given by z_i_L_1_i^2 + 2 w_i^⊤ L_2_i z_i + w_i_L_3_i^2 ≥ 0 for all z_i satisfying z_i-c_z_i_Ẑ_i^2 ≤ 1 and w_ij satisfying w_ij-a_ij^max(t) c_y_j_Y_j^2 ≤ a_ij^max^2(t) for all j ∈𝒩_i^- where L_1_i=P_i-Γ_i-(A_i+B_i K_o_i T_i(t))^⊤ P_i (A_i+B_i K_o_i T_i(t)), L_2_i=0.5C_i-E_i^⊤ P_i (A_i+B_i K_o_i T_i(t)) and L_3_i=D_i-E_i^⊤ P_i E_i.
    
    We now define the block diagonal matrix Y_σ_i(t) with the submatrices σ_ij(t) Y_j for all j ∈𝒩_i^- on the diagonal and the vector τ_i(t) which vertically concatenates the subvectors a_ij^max(t) c_y_j for all j ∈𝒩_i^-. 
    We also recall that the ellipsoidal set 𝒴̂_i contains the origin in its interior and hence, c_y_i_Y_i^2 ≤ 1.
    Using the S-lemma, applying the Schur complement and multiplying the resulting matrix inequality by α_i(t) thus leads to (<ref>).
    

    
    
    
    
    



















































    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    













    □

Note that (<ref>) implies that P_i-(A_i+B_i K_o_i T_i(t))^⊤ P_i (A_i + B_i K_o_i T_i(t)) ≥Γ_i(t) where Γ_i(t) ∈𝔻_++^n_i and hence, A_i+B_iK_i(t) is Schur. 

One source of conservatism in Proposition <ref> is that the local controller is designed to not only stabilize the local dynamics, but also to passivize it. The benefit of passivizing the local dynamics will become obvious in the next proposition.
We now consider Condition III in Section 2 by deriving constraints which ensure that the strict passivity of all noise-free subsystems implies the asymptotic stability of the overall noise-free system. For this purpose, we consider the noise-free global dynamics
z(t+1) = A_cl(t) z(t) + E w(t)
where A_cl(t)=diag(A_1+B_1K_1(t),,A_M+B_MK_M(t)), E=diag(E_1,,E_M) and w(t)=[w_1^⊤(t),,w_M^⊤(t)]^⊤.
We also define the Lyapunov function V(z) = z_P^2 where P=diag(P_1,...,P_M). Furthermore, we define the matrices Γ(t) = diag(Γ_1(t),...,Γ_M(t)), C = diag(C_1,...,C_M) and D(t) = diag(D_1(t),...,D_M(t)). 
Recall that A_g and D_g are the adjacency and degree matrices, respectively. The global virtual output vector is given by ỹ(t) = C z(t) + D w(t) and the global coupling vector is given by w(t) = Ã_g C z(t). 

Finally, we construct A_g^max(t) and D_g^max(t) by replacing all uncertain parameters in A_g and D_g with their upper bounds. Consequently, we define U=diag(U_1,...,U_M)=C^⊤ D_g C, V=[V_1^⊤,...,V_M^⊤]^⊤=C^⊤Ã_g^⊤ and W=[W_1^⊤,...,W_M^⊤]^⊤=Ã_g C where U_i ∈ℝ^n_i × n_i, V_i ∈ℝ^n_i × m_0 and W_i ∈ℝ^m_0 × n_i. We also construct U^max(t), V^max(t) and W^max(t) by replacing all uncertain parameters in U, V and W with their upper bounds. In the sequel, we denote the diagonal element in the j-th row of a matrix by [·]_j and the summation of the absolute values of the entries in the j-th row of a matrix by |·|_j. The dependence on time in (<ref>) is omitted in the interest in space.


    
    The strict passivity of the noise-free dynamics of all subsystems under the corresponding controller u_i(t)=K_o_iT_i(t)x_i(t) (i.e.  z_i(t+1)=(A_i+B_iK_o_iT_i(t))z_i(t) + E_iw_i(t)) implies the asymptotic stability of the noise-free dynamics of the overall system (i.e.  z(t+1) = A_cl(t) z(t) + E w(t)) if,
    
    ([Γ_i]_j-ϵ) α_i ≥ |U_i^max|_j α_i + |V_i^max|_j α_i, ∀ j ∈{1,...,n_i},

    

    
    [D_i]_j α_i ≤α_i/|W_i^max|_j,  ∀ j ∈{1,...,m_0} :  |W_i^max|_j ≠ 0.

    




    The passivity of the noise-free dynamics of each subsystem implies that z_i^+_P_i^2- z_i_P_i^2 ≤ỹ_i^⊤ w_i - z_i_Γ_i(t)^2. Summing up these inequalities for all subsystems yields V(z^+)-V(z) ≤ỹ^⊤ w - z_Γ(t). To ensure asymptotic stability, it suffices to guarantee that ỹ^⊤ w - z_Γ(t) < 0. Using the expressions of ỹ and w, this condition can be written as z^2_Γ(t) - C^⊤Ã_g C - C^⊤Ã_g^⊤ D(t) Ã_g C≥ϵz^2 where ϵ > 0. Recall that L_g = D_g - A_g and hence, z^2_Γ(t) + C^⊤ L_g C - C^⊤ D_g C - C^⊤Ã_g^⊤ D(t) Ã_g C≥ϵz^2. Note also that L_g is positive semidefinite by definition and hence, it is sufficient to ensure that z^2_Γ(t) - C^⊤ D_g C - C^⊤Ã_g^⊤ D(t) Ã_g C≥ϵz^2. This condition can be reformulated using the Schur Complement as,
    
    [
            [ Γ(t) - C^⊤ D_g C -ϵ I_n                C^⊤Ã_g^⊤; [2pt/2pt]
                *                 D(T)^-1;                         ]]
            ≥ 0.

    Using diagonal dominance, the constraint (<ref>) is approximated for each subsystem by [Γ_i(t)]_j -ϵ≥ |U_i|_j + |V_i|_j for all j ∈{1,...,n_i} and [D_i(t)^-1]_j ≥ |W_i|_j for all j ∈{1,...,m_0}. Multiplying the former by α_i(t) and the latter by α_i(t)^-1 and inverting, these inequalities become ([Γ_i(t)]_j-ϵ) α_i(t) ≥ |U_i|_j α_i(t) + |V_i|_j α_i(t) for all j ∈{1,...,n_i} and
    [D_i(t)]_j α_i(t) ≤α_i(t)/|W_i|_j for all j ∈{1,...,m_0} such that |W_i|_j ≠ 0. Note that the resulting inequalities are functions of the uncertain parameters. More conservative inequalities can be obtained by replacing the uncertain parameters with their upper bounds, leading to (<ref>). 
    □


In Proposition <ref>, (<ref>) ensures that the Lyapunov condition is satisfied for the noise-free global dynamics z(t+1)=A_cl(t)z(t)+Ew(t) or equivalently, z(t+1)=(A+BK(t))z(t) and hence, A+BK(t) is Schur. Sources of conservatism in this proposition include implying ỹ^⊤ w - z_Γ(t) < 0 instead of V(z^+)-V(z) < 0, ignoring the term C^⊤ L_g C, using diagonal dominance and replacing the true values of the parameters with their upper bounds. These sources are utilized to eventually reach local convex constraints.
We now move to Conditions IV and V in Section 2.


In the sequel, we denote the j-th row of the matrices G_i and H_i by G_i^(j) and H_i^(j) respectively.
The dependence on time is omitted in (<ref>) and (<ref>) in the interest in space.


    
    The Minkowski sum of 𝒳_i(t) and 𝒵_K_i(t) is contained in 𝒵_i if there exist μ_ij(t) ≥ 0 for all j ∈{1,...,q_i} such that,
    
    [
            [              μ_ij(t)Z_i     -0.5G_i^(j)^⊤α_i(t); [2pt/2pt]
                *       1-μ_i_j(t)-g_i(t) ]]
            ≥ 0.




    The Minkowski sum of 𝒳_i(t) and 𝒵_K_i(t) is contained in 𝒵_i if and only if for all j ∈{1,...,q_i}, G_i^(j) (x_e_i+x_s_i) ≤ 1 for all x_e_i such that x_e_i_Z_i^2 ≤α_i^2(t) and for all x_s_i such that G_i x_s_i≤ g_i(t). 
    By defining s_i and t_i such that x_e_i=α_i(t) s_i and x_s_i=g_i(t) t_i, this condition becomes G_i^(j)α_i(t) s_i + G_i^(j) g_i(t) t_i ≤ 1 for all j ∈{1,...,q_i}, s_i satisfying s_i^⊤ Z_i s_i ≤ 1 and t_i satisfying G_i t_i ≤ 1. Using the S-lemma, this condition is given for all j ∈{1,...,q_i} by,
    
    [
            [                                             μ_ij(t) Z_i                                                       0                                  - 0.5 α_i(t) G_i^(j)^⊤;                                 [2pt/2pt]
                *                                                       0 ∑_l=1^q_i 0.5 μ_ijl(t) G_i^(l)^⊤ - 0.5 G_i^(j)^⊤ g_i(t);                                 [2pt/2pt]
                *                                                       *                             1-μ_ij(t)-∑_l=1^q_iμ_ijl(t);                                                         ]]
            ≥ 0.
    
    Note that (<ref>) holds if ∑_l=1^q_i 0.5 μ_ijl(t) G_i^(l) = 0.5 G_i^(j) g_i(t). By stacking this equality for all j ∈{1,...,q_i}, we reach g_i(t) G_i = μ_i(t) G_i where μ_ijl(t) is the entry of μ_i(t) in the j-th row and l-th column. To satisfy this equality, we choose the diagonal elements of μ_i(t) to be equal to g_i(t) and the off-diagonal elements to be equal to zero and hence, (<ref>) reduces to (<ref>).
    
    
    
    
    
    
    
    
    
    
    
    □



    
    The Minkowski sum of 𝒰_i(t) and K_o_iT_i(t)𝒵_k_i(t) is contained in 𝒱_i if there exist ν_ij(t) > 0 for all j ∈{1,...,r_i} and M_i(t) ∈𝕊^p_i_++ such that,
    
    [
            [                  M_i^-1                K_i^⊤α_i; [2pt/2pt]
                *                ν_ij Z_i;                         ]]
            ≥ 0,
            [
            [                  M_i^-1    0.5 M_i^-1 H_i^(j)^⊤; [2pt/2pt]
                *          1 - ν_ij - h_i;                         ]]
            ≥ 0.




    First, we define M_i(t) ∈𝕊_++^p_i and ν_ij(t) > 0 such that M_i(t) ≤ν_ij(K_o_i T_i(t) α_i(t) Z_i^-1α_i(t) T_i^⊤(t) K_o_i^⊤)^-1 for all j ∈{1,...,r_i}. Since T_i(t) ∈𝕊_++^n_i, α_i(t) > 0 and K_o_i is full-rank, this inequality is invertible and hence, M_i(t)^-1 - ν_ij(t)^-1 K_o_i T_i(t) α_i(t) Z_i^-1α_i(t) T_i^⊤(t) K_o_i^⊤α_i(t) ≥ 0  for all j ∈{1,...,r_i}.
    Applying the Schur Complement yields the first condition in (<ref>). 
    The Minkowski sum of 𝒰_i(t) and K_o_iT_i(t)𝒵_K_i(t) is contained in 𝒱_i if and only if for all j ∈{1,...,r_i}, H_i^(j) (u_e_i+u_s_i) ≤ 1 for all u_e_i satisfying u_e_i_(K_o_i T_i(t) Z_i^-1 T_i^⊤(t) K_o_i^⊤)^-1^2 ≤α_i^2(t) and u_s_i satisfying H_i u_s_i≤ h_i(t). 
    By defining s_i and t_i such that u_e_i= s_i and u_s_i=h_i(t) t_i, we reach the condition H_i^(j) s_i + H_i^(j) h_i(t) t_i ≤ 1 for all j ∈{1,...,r_i}, s_i and t_i satisfying s_i^⊤ (K_o_i T_i(t) α_i(t) Z_i^-1α_i(t) T_i^⊤(t) K_o_i^⊤)^-1 s_i ≤ 1 and H_i t_i ≤ 1. 
    Using the S-lemma, we reach  for all j ∈{1,...,r_i}, 
    
    [
            [                                           ν_ijψ_i^-1                                                    0                                      - 0.5 H_i^(j)^⊤;                              [2pt/2pt]
                *                                                    0 ∑_l=1^r_i 0.5 ν_ijl H_i^(l)^⊤ - 0.5 H_i^(j)^⊤ h_i(t);                              [2pt/2pt]
                *                                                    *                                1-ν_ij-∑_l=1^r_iν_ijl;                                                      ]]
            ≥ 0,

    where ψ_i(t) = K_o_i T_i(t) α_i(t) Z_i^-1 T_i^⊤ K_o_i^⊤(t) α_i(t).
    Note that (<ref>) holds if ∑_l=1^r_i 0.5 ν_ijl(t) H_i^(l) = 0.5 H_i^(j)^⊤ h_i(t). 
    By stacking this equality for all j ∈{1,...,r_i}, we reach h_i(t) H_i = ν_i(t) H_i where ν_ijl(t) is the entry of ν_i(t) in the j-th row and l-th column. 
    To satisfy this equality, we choose the diagonal elements of ν_i(t) to equal h_i(t) and the off-diagonal elements to equal zero. 
    Recalling that M_i(t) ≤ν_ij(t) (K_i(t) α_i(t) Z_i^-1 K_i^⊤(t) α_i(t))^-1, (<ref>) reduces to a matrix, that if pre and post-multiplied by diag(M_i(t)^-1,1), yields the second condition in (<ref>).
    
    
    
    
    
    
    
    
    
    
    
    
    □




One source of conservatism in Proposition <ref> is choosing the diagonal elements of μ_i(t) to equal g_i(t) and the off-digonal elements to equal zero. The same applies to ν_i(t) in Proposition <ref>. Additionally, introducing M_i(t) adds one more degree of conservatism in Proposition <ref>.
In conclusion, the OCP of the i-th subsystem becomes,

    min 
            
            
            
            J_i 
            
             s.t.{   ∀ k ∈{0,...,N-1}, 
        z_i(t) - x_i(0|t) ∈𝒵_K_i (t), 
        x_i(k+1|t)=A_ix_i(k|t)+B_iu_i(k|t), 
        x_i(k|t) ∈𝒳_i(t),    u_i(k|t) ∈𝒰_i(t), 
        x_i(N|t) = x_e_i(t)=A_ix_e_i(t)+B_iu_e_i(t), 
        x_e_i∈ξ𝒳_i(t),    u_e_i∈ξ𝒰_i(t), 
       (<ref>),  λ_i(t),  λ_d_i(t),  λ_ij(t) ≥ 0,  ∀ j ∈𝒩_i^-,
       (<ref>), (<ref>),  (<ref>),  μ_ij(t) ≥ 0  ∀ j ∈{1,...,q_i}, 
       (<ref>),  ν_ij(t) > 0  ∀ j ∈{1,...,r_i}, 
        M_i(t) ∈𝕊_++^n_i,  T_i(t) α_i(t) ∈𝕊_++^n_i, 
       Γ_i(t) α_i(t) ∈𝔻_++^n_i,  D_i(t) α_i(t) ∈𝔻_++^n_i.
            .

The decision variables are x_i(k|t) and u_i(k|t) for all k ∈{0,...,N}, x_e_i(t), u_e_i(t), g_i(t), h_i(t), α_i(t), T_i(t) α_i(t), Γ_i(t) α_i(t), D_i(t) α_i(t), M_i(t)^-1, λ_i(t), λ_d_i(t), λ_ij(t) for all j ∈𝒩_i^-, μ(t)_ij for all j ∈{1,...,q_i} and ν(t)_ij for all j ∈{1,...,r_i}. 
Although the OCP might seem nonconvex using these decision variables, one can turn it into an SDP by redefining the decision variables and, if desired, reversing the transformation to recover the original variables from the solution.
As the number of neighbours of each subsystem increases, the dimensions of all constraints remain unchanged except for (<ref>) which increases linearly and hence, the computation time is not significantly altered.
Similar to the OCP (<ref>), the online OCP (<ref>) is solved at each time instant and the optimal control input v_i(t)=u_i(0|t)+K_o_iT_i(t)(z_i(t)-x_i(0|t)) is applied to the i-th subsystem.























§ FEASIBILITY AND STABILITY


In this section, we prove that the online OCP (<ref>) is recursively feasible and the resulting closed-loop system is input-to-state stable. In the sequel, we use the superscript (·)^* to refer to the optimal solution. The proof is inspired from <cit.> and is divided into three main steps. First, we discuss the convergence of the nominal state x^*(t)=[x_1^*^⊤(t),,x_M^*^⊤(t)]^⊤ to the artificial equilibrium x_e^*(t)=[x_e_1^*^⊤(t),,x_e_M^*^⊤(t)]^⊤. Then, we discuss the convergence of the artificial equilibrium x_e^*(t) to the target point x_r=[x_r_1^⊤,,x_r_M^⊤]^⊤. Finally, we discuss the stability of the error dynamics between the actual state z^*(t)=[z_1^*^⊤(t),,z_M^*^⊤(t)]^⊤ and the nominal state x^*(t).
We start with the first step in Lemma <ref>.




    
    The nominal state trajectory converges asymptotically to the artificial equilibrium trajectory.

The proof uses standard MPC arguments (<cit.>) by comparing the optimal cost at time t, the tail sequence cost at time t and the optimal cost at time t+1.


    
    
    
    
    
    
    
    
    
    



    
Second, we discuss the convergence of x_e^*(t) to the target point x_r. For this purpose, we define (x̃_r,ũ_r) to be the projection of (x_r,u_r) onto the set of equilibrium points contained in the set (ξ𝒳_i(t),ξ𝒰_i(t)). We also define the matrix P̃ such that P̃-(A+BK(t))^⊤P̃ (A+BK(t)) - Q - K^⊤(t) R K(t) ≥ 0. Note that such P̃ exists since K(t) asymptotically stabilizes the actual global dynamics as shown in Proposition <ref>. Finally, we define E(c,s) = {x:x-c^2_P̃≤ s}. For brevity, we omit the dependence on t in the proof of Lemma <ref>. This proof is inspired from <cit.>, however, we prove that lim_t →∞ (x^*(t)-x_e^*(t))=0 implies lim_t →∞ (x_e^*(t)-x_r)=0 instead of x^*(t)-x_e^*(t)=0 implying x_e^*(t)-x_r=0.
    
    
        
        
        If x^*(t) and x_e^*(t) are such that lim_t →∞ (x^*(t)-x_e^*(t)) = 0, then lim_t →∞ (x_e^*(t)-x̃_r) = 0.
    
    
    
        According to the OCP (<ref>), (x_e^*,u_e^*) ∈int(𝒳×𝒰). 
        Hence, there always exists E(x_e^*,β_1) such that (x,K(x-x_e^*)+u_e^*) ∈γ𝒳×γ𝒰 for all x ∈ E(x_e^*,β_1) where γ∈ (0,1). 
        Note that E(x_e^*,β_1) is invariant since it is defined using the level sets of a Lyapunov function. 
        We now define an equilibrium point (x̅_e,u̅_e) ∈int(𝒳×𝒰) such that x̅_e-x̃_r_S^2 = β_2 x_e^*-x̃_r_S^2 and x̅_e-x_e^*^2_P̃=(1-λ)x̃_r-x_e^*^2_P̃ where β_2 ∈ (0,1) and λ∈ (0,1) is chosen sufficiently large such that (0,K(x̅_e-x_e^*)+u̅_e-u_e^*) ∈ (1-γ)𝒳× (1-γ) 𝒰. 
        Next, we define E(x̅_e,β_3) ⊂ E(x^*_e,β_1) such that x_e^* ∈int(E(x̅_e,β_3)). 
        This set is invariant since it is defined using the level sets of a Lyapunov function. Moreover, for all x ∈ E(x̅_e,β_3), (x,K(x-x̅_e)+u̅_e)=(x,K(x-x_e^*)+u^*_e)+(0,K(x̅_e-x_e^*)+u̅_e-u_e^*) ∈( γ𝒳×γ𝒰) ⊕((1-γ) 𝒳× (1-γ) 𝒰) = 𝒳×𝒰.
        
        We now define β_4 > 0 such that S>β_4P̃ and the set 𝒳̂ = {x:x-x_e^*_P̃^2 ≤β_5 x_e^*-x̃_r_P̃^2  & (x_e^*-x̅_e)^⊤P̃ |x-x_e^*| ≤β_5 x_e^*-x̃_r_P̃^2 } where 0<β_5<β_4(1-β_2)/3. 
        Hence, x-x̅_e_P̃^2 
        = 
        x-x^*_e_P̃^2 + 2 (x_e^*-x̅_e)^⊤P̃ (x-x_e^*) + (1-λ)x_e^*-x̃_r_P̃^2 
        ≤x_e^*-x̃_r_3β_5P̃+(1-λ)P̃^2 for all x ∈ E(x_e^*,β_6) ⊂𝒳̂. 
        Note that E(x_e^*,β_6) can be also chosen to lie inside E(x̅_e,β_3) ⊂ E(x_e^*,β_1) and hence is invariant. 
        Note also that x-x̅_e_P̃^2+x̅_e-x̃_r_S^2 ≤x_e^*-x̃_r_S^2 for all x ∈ E(x_e^*,β_6) if (1-β_2)S-(1-λ+3β_5) P̃≥ 0 or equivalently, β_4(1-β_2)-(1-λ+3β_5) ≥ 0. This can be achieved if λ≥ 1-(β_4(1-β_2)-3β_5). Since β_5<β_4(1-β_2)/3, this condition on λ does not contradict the aforementioned conditions on λ.
        
        We conclude the proof using contradiction. Assume that lim_t →∞ (x^*(t)-x_e^*(t)) = 0 and lim_t →∞ (x_e^*(t) - x̃_r) ≠ 0. Hence, it is always possible to pick an arbitrarily large t such that x^* ∈ E(x_e^*,β_6) ⊂ (x̅_e,β_3). In addition, there exists β_7 such that x_e^*-x̃_r≥β_7 for infinitely many t. Since E(x_e^*,β_6) and E(x̅_e,β_3) are both invariant and satisfy the state and input constraints, let J^* and J̅ be the cost corresponding to the optimal solution and the feasible solution aiming to converge to x̅_e and note that J^* < J̅. Note also that J̅≤x^*-x̅_e_P̃^2+x̅_e-x̃_r_S^2 due to the Lyapunov inequality that P̃ satisfies. Since x^* ∈ E(x_e^*,β_6), then x^*-x̅_e_P̃^2+x̅_e-x̃_r_S^2 ≤x_e^*-x̃_r_S^2. But, x_e^*-x̃_r_S^2 ≤x_e^*-x_r_S^2 ≤ J^* which contradicts the optimality of J^*.
        □
    
    
    
    
    Finally, we discuss the input-to-state stability of the global closed-loop system under the MPC scheme (<ref>). Recall from (<ref>) that the global system is given by z(t+1)=(A_d+EÃ_gC)z(t)+Bv(t)+d(t). Recall also that the optimal input is given by v^*(t)=u^*(t)+K^*(t)(z^*(t)-x^*(t)) where K^*(t)=diag(K_o_1T_1^*(t),,K_o_MT_M^*(t)) and u^*(t)=u^*(0|t)=[u_1^*^⊤(0|t),,u_M^*^⊤(0|t)]^⊤. Hence, the closed-loop system is given by z^*(t+1)=A_dz^*(t)+Bu^*(t)+BK^*(t)(z^*(t)-x^*(t))+EÃ_gCz^*(t)+d(t).
    By defining e^*(t)=z^*(t)-x^*(t) and η(t)=d(t)-x^*(t+1)+Ax^*(t)+Bu^*(t), the closed-loop system becomes,
    
    e^*(t+1) = (A+BK^*(t))e^*(t)+η(t).

    While e^*(t) is the state of this system, x^*(t) and u^*(t), together with d(t), are exogenous signals to this system. While d(t) is bounded by assumption, x^*(t) and u^*(t) are bounded since they are the outputs of the proposed MPC problem (15) which consider compact constraint sets. Hence, we prove input-to-state stability with respect to these bounded exogenous signals. Note that x^*(t+1) is not necessarily equal to Ax^*(t)+Bu^*(t) where A=diag(A_1,,A_M)+EA_gC. This is because x^*(·) and u^*(·) are the nominal states and inputs which satisfy x^*(t+1)=diag(A_1,,A_M)x^*(t)+Bu^*(t).   
    
    
    Moreover, the considered system is switched since K^*(t) is adapted at each time instant. Hence, the closed-loop system under the proposed MPC scheme can be considered as an uncertain switched system with bounded disturbances.







    
    
    
    
    



    
    
        The MPC scheme (<ref>) is recursively feasible and the global closed-loop system (<ref>) under this scheme is input-to-state stable with respect to η(t).
        
        
        
        
        
        
        
    
    
    
        Recursive feasibility is proven using standard MPC arguments (<cit.>) by making use of the tail sequence. Regarding stability, we first consider the disturbance-free system e^*(t+1)=(A+BK^*(t))e^*(t). The matrix A+BK^*(t) is always Schur for any K^*(t) since the constraints derived in Propositions <ref> and <ref> ensure that P-(A+BK^*(t))^⊤ P (A+BK^*(t)) ≥Γ(t) - C^⊤Ã_g C - C^⊤ L D(t) L C ≥ϵ I_n at the t-th time instant for all t. Since this condition is always satisfied using the same matrix P, then V(e^*)=e^*^2_P is a Lyapunov function for the disturbance-free system independent of K^*(t).
        Now, we consider the closed-loop system (<ref>). Using the Lyapunov function V(e^*) and defining η(t)=d(t)-x^*(t+1)+Ax^*(t)+Bu^*(t) , we can deduce that V((A+BK^*(t))e^*(t)+η(t))-V(e^*) ≤ -(ϵ-δ)e^*(t)^2 - δe^*(t)^2 + 2η(t)_P (A+BK^*(t))e^*(t)_P + η(t)_P^2 where 0 < δ < ϵ. To prove input-to-state stability, it suffices to ensure that δe^*(t)_2^2 ≥ 2η(t)_P(A+BK^*(t))e^*(t)_P + η(t)_P^2 for some ζ > 0 such that e^*(t)≥ζ. 
        We denote the maximum eigenvalue of P by λ_max(P) and the maximum norm of η(t) by η_max. Since (A+BK^*(t))e^*(t)_P < e^*(t)_P and (·)_P ≤√(λ_max(P))(·), it suffices instead to ensure that δe^*(t)^2 - 2 λ_max(P)η_maxe^*(t)-λ_max(P)η_max^2 ≥ 0. By computing the roots of the convex quadratic function on the left-hand side of the inequality, it is easy to verify that this inequality holds for e^*(t)≥ζ = (λ_max(P)+√(λ_max^2(P)+δλ_max(P)))η_max/δ. Hence, input-to-state stability is guaranteed.
        □
        
        
        
        
        
        
        
        
        
        
        
        
        
    
    
    



§ NUMERICAL RESULTS



In this section, we investigate the efficacy of the proposed scheme by comparing it to existing schemes in the literature. In particular, we compare the proposed scheme (ADP+LRN) to robust MPC (ROB) (<cit.>), adaptive MPC (ADP) (<cit.>) and learning-based MPC (LRN) (<cit.>). Note that ADP is an extension of the adaptive MPC scheme in <cit.> to the class of interconnected systems. We solve all optimization problems using MATLAB with YALMIP (<cit.>), Mosek and Gurobi on a computer equipped with a 1.9-GHz Intel core i7-8550U processor.

In the simulation example, we use a network of five double integrators connected in series as in <cit.>. 
We consider each double integrator as a subsystem with the dynamics in (<ref>) where A_i = [ 1 1; 0 1 ], B_i=[0.5  1]^⊤, C_i=[0  1] and E_i=[0.05  0.1]^⊤ for all i ∈{1,2,3,4,5}. 
All states and inputs are constrained between -2 and 2. 
The cost function matrices are Q_i=I_2, R_i=0.1 and S_i=10I_2 for all subsystems and the prediction horizon is 5. 
The noise d(t) ∈ℝ^2 is extracted uniformly with the first entry lying between -0.025 and 0.025 and the second between -0.05 and 0.05. 
The unknown true value of the parameters a_12 and a_54 is 1, whereas that of a_21, a_23, a_32, a_34, a_43 and a_45 is 0.5. The parameters a_12 and a_54 are assumed to lie initially between 0 and 4, whereas all other parameters are assumed to lie initially between 0 and 2.
The unknown true values of the parameters are not used within the proposed framework and, by design, communication is not used in the adaptation phase to locally solve the OCPs.
Each subsystem is required to follow a randomly generated piecewise constant reference changing every 100 seconds and constrained between -2 and 2. The simulation time is 1000 seconds and the initial condition is [0  0]^⊤ for all subsystems.







Fig.<ref>(Top) shows the closed-loop performance of the third subsystem for each scheme. As shown in this figure, ROB and LRN generally fail to track the reference trajectory since LRN updates the model only while keeping the initial tube and tightened constraint sets unchanged and ROB updates nothing. On the contrary, ADP follows the reference successfully since it updates the tube and the tightened constraint sets online. Similarly, ADP+LRN generally follows the reference (due to the online shrinkage of the uncertain parameter sets as shown in Fig.<ref>(Bottom) for the parameter a_31) expect for extreme references close to the actual state constraint bounds due to the conservative constraint tightening.


The resulting behavior of each scheme has a direct effect on the closed-loop cost shown in Fig.<ref>(Top). ROB has the highest cost, followed by LRN whose cost is slightly lower due to updating the model. The costs of ADP+LRN and ADP are considerably lower with ADP having the lowest cost. 


This is because ADP can quickly track the reference, even in extreme cases, since this scheme can directly handle parametric uncertainties unlike ADP+LRN which lumps all uncertainties in an additive term. 
ADP+LRN, however, results in a decentralized OCP which requires significantly reduced computations and communication among neighbours. On the other hand, ADP yields a distributed OCP that is solved using distributed optimization techniques (e.g. ADMM) leading to extensive communication and computations.



Fig.<ref>(Bottom) shows the mean and standard deviation of the computation time required by each scheme (including different versions of ADP which differ in the maximum ADMM iterations allowed) and the corresponding closed-loop cost. We use the ADP scheme whose computation time is around 1 second in the other figures. It is found, however, that if we run ADMM for around 0.1 second only, we still get a relatively satisfactory performance. Running ADMM for less time (e.g. same order of magnitude of the time required by ADP+LRN) result in severe closed-loop performance deterioration.

On the other hand, all other schemes have lower computation time due to their decentralized nature which do not require using distributed optimization techniques such as ADMM. Since the OCP of ADP+LRN is a semidefinite program, ADP+LRN requires slightly higher computation compared to LRN and ROB whose OCPs are quadratic programs. 





§ CONCLUSION


We propose an adaptive learning-based MPC scheme for uncertain interconnected systems. 

The developed scheme is found to yield a better trade-off between computation complexity and closed-loop performance compared to existing schemes. 

Future work include applying the developed scheme to real-world applications.


                               
The authors are grateful to Prof. Roy Smith, Prof. Andrea Iannelli, Dr. Goran Banjac and Dr. Anil Parsi for the fruitful discussions on the topic.


agsm 
 






















































