Robust coalitional model predictive control with predicted topology transitions

This is a repository copy of Robust coalitional model predictive control with predicted topology transitions in the Depósito de Investigación de la Universidad de Sevilla. Version: Author Accepted Version. Citation: Masero, E., Maestre, J. M., Ferramosca, A., Francisco, M., Camacho, E.F. (2021). Robust coalitional model predictive control with predicted topology transitions. IEEE Transactions on Control of Network Systems, DOI: https://doi.org/10.1109/TCNS.2021.3088806 To cite this publication, please use the final published version (if applicable). Please check the document version above. Copyright: © 2021 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. Takedown policy: Please contact us ([email protected]) and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim. 1 Robust coalitional model predictive control with predicted topology transitions Eva Masero, Jos´ e M. Maestre, Antonio Ferramosca, Mario Francisco, and Eduardo F. Camacho Abstract—This paper presents a novel clustering model predictive control technique where transitions to the best cooperation topology are planned over the prediction horizon. A new variable, the so-called transition horizon, is added to the optimization problem to calculate the optimal instant to introduce the next topology. Accordingly, agents can predict topology transitions to adapt their trajectories while optimizing their goals. Moreover, conditions to guarantee recursive feasibility and robust stability of the system are provided. Finally, the proposed control method is tested via a simulated eight-coupled tanks plant. Index Terms—Model predictive control; Control by clustering; Distributed control; Coalitional control; Networked control. I. INTRODUCTION MODEL predictive control (MPC) is an optimizationcontrol method that uses a model to predict the system behavior through a given time horizon. At each time step, the control sequence that minimizes a cost function based on the system evolution is computed subject to a set of constraints. Only the first control signal of the obtained sequence is applied, and the rest is discarded. In the next instant, the time window of the problem is displaced one step ahead following a receding horizon strategy, and the same procedure is repeated. Compared with other methods, MPC presents major advantages, such as the ability to deal explicitly with delays and dead times. It also handles the multi-variable, nonminimum phase, and unstable systems [1]. For these reasons, this method is widely applied in the process industry [2]. Nevertheless, the application of MPC in large-scale processes is sometimes impossible due to the incapacity to find a centralized model of the overall system and the high computational requirements, to name two typical problems in this context. In these cases, a distributed approach may be required. Indeed, the improvements in computational and communication technologies made distributed MPC (DMPC) a popular framework with applications in large-scale systems such as road traffic networks [3], power grids [4], gas networks [5], irrigation canals [6], water distribution systems [7], supply chains [8], and energy management in buildings [9], [10], among others. E. Masero, J. M. Maestre and E. F. Camacho are with the Dept. Ingenier´ ıa de Sistemas y Autom´ atica, Escuela T´ ecnica Superior de Ingenier´ ıa, Universidad de Sevilla, C/ Camino de los Descubrimientos, s/n., Seville 41092, Spain. J. M. Maestre is also currently with the Graduate School of Science and Technology, Keio University, Yokohama, Japan. E-mails: {evamasero, pepemaestre, efcamacho}@us.es Antonio Ferramosca is with Dept. of Management, Information and Production Engineering, University of Bergamo, Via Marconi 5, 24044 Dalmine (BG), Italy. E-mail: [email protected] Mario Francisco is with the Dept. Inform´ atica y Autom´ atica, Escuela T´ ecnica Superior de Ingenier´ ıa Industrial, Universidad de Salamanca, Av. Fernando Ballesteros, B´ ejar 37700, Salamanca, Spain. E-mail: [email protected] The main idea of DMPC is to divide the overall system into subsystems governed by local controllers or agents with individual control goals that may cooperate to improve performance; see [11], [12] for further details. Recently, it has been proposed to adapt the cooperation burden in real-time depending on the coupling between control tasks by promoting the formation of coalitions or clusters of cooperating controllers. This approach is useful for autonomous systems that perform cooperative tasks and deal with saving energy, such as UAVs and robots. For example, [13] proposed an adaptive cluster formation scheme based on active coupling constraints, and [14] presented hierarchical controllers that are dynamically adapted to the operational conditions. Furthermore, the controller topology can be selected by a supervisory controller (top-down approach) [15], [16], or be obtained disabling communication links with a poor contribution to the local performance (bottom-up architecture) [17], [18]. Time 𝑘 𝑁t 𝑁p 𝑘 − 1 𝑘 + 2 𝑘 + 𝑁t𝑘 + 4 𝑘 + 5 𝑘 + 6 𝑘 + 7 𝑘 + 8 𝑘 + 9 𝑁p 𝑁p 𝑁p 𝑘 + 1 Current Topology Next Topology Fig. 1. Example scheme of the prediction horizon Npand the transition horizon Nt. This work presents a novel hierarchical coalitional MPC method, where the information is exchanged between two control layers. In the lower layer, local controllers optimize their control and state sequences, and the upper layer evaluates the information received from agents and chooses the best cooperation topology to attain a global objective, including performance and cooperation costs. Once the topology is selected, it is sent to the lower layer, where local controllers must follow it. The instant when they adopt the new topology affects the results. Hitherto, the new topology was chosen and abruptly implemented by local controllers at the same instant, as presented in [17], [15]. In this work, the moment to switch topology becomes a variable of the optimization problem, the so-called transition horizon. Unlike the prediction horizon, which is displaced toward the future at each time instant, the transition horizon is not moved ahead but actualized at each time instant as depicted in Fig. 1. In other words, the optimization problem solved by the upper control layer considers the cost over the prediction horizon in a receding horizon fashion with the moment to switch topology being also optimized. Hence, our approach 2 provides more degrees of freedom for optimization, and the resulting performance will always be no worse regarding the MPC receding horizon objective. It also results in the controller’s capacity to gradually prepare for a switch of topology within the prediction horizon (see Fig. 1) and, hence, adapt its control sequence accordingly. Regarding robustness, several techniques can be found in the literature to deal with constrained systems’ stability in the presence of bounded additive disturbances. Some of them are the input-to-state (ISS) property [19], widely used to analyze non-linear systems; the min-max strategy, where the optimization problem is solved for the worst-case disturbances [20]; and the tube approach [21], [22], which is suitable for linear systems, and based on computing a region around the nominal trajectory that contains the system state under any uncertainties. Our approach relies on the inherent robustness of the predictive controllers due to the feedback nature of MPC [23], [24], [25]. To illustrate the benefits of our strategy, an eight-coupled tanks plant, which is an extension of the quadruple-tank process proposed by [26], has been designed. The quadrupletank process has popularly been used as a benchmark, e.g, to study the effects of multi-variable plants with dead times [27]. Finally, a much earlier and abridged version of the proposed method was accepted for presentation at a conference [28]. There are significant differences between the current article and the conference paper: i) the conference version does not include proofs for its stability claims; ii) robustness is now considered and guarantees are given and proved in this regard; iii) we present here new improved methods for the selection of topology; iv) several refinements in the main algorithm have been introduced in the current version; v) the plant configuration has been modified to couple subsystems by the states; and vi) the experiments performed are new for the article version and include the previously mentioned novelties. Index of contents: Section II introduces the problem formulation. Section III defines the control objective and the coalitional MPC algorithm. Controllers are designed in Section IV. Section V provides recursive feasibility and stability guarantees. Section VI details the eight-coupled tanks used as a benchmark, the control techniques assessed, and results via simulation. Section VII provides conclusions and discussion. Notation: N0`and N`are the sets of non-negative and positive integers. Rnrefers to an n-dimension set of real numbers. For sets X,YĎRn, the Cartesian product is XˆYfitpx, yq:xPX, y PYu. If tXiuiPNis a finite family of sets indexed by N, then the Cartesian product ŚiPNXiis X1ˆX2ˆ¨¨¨ˆXN“ tpx1, x2, . . . , xNq:x1P Xi, . . . , xNPXNu. The set subtraction operator is z. The image of set XĎRnunder a linear mapping A:RnÞÑ Rm is AXfitAx :xPXu. For sets X,YĎRn, the Minkowski sum is X‘Yfitx`y:xPX, y PYu; and for YĎX, the Pontryagin difference is XaYfitzPRn:Y‘tzu Ď Xu. Matrix Iis the identity matrix. The lz-norm defined in Rnis represented as ||¨||zwith zPN`. II. PROBLEM SETTING In this section, the problem formulation is provided. The system dynamics are firstly introduced, and, subsequently, it is explained how information is exchanged between agents. A. System dynamics and constraints The system is divided into N“ t1,2, . . . , Nucoupled subsystems whose discrete-time dynamics are modeled as xipk`1q “ Aiixipkq`Biiuipkq`wipkq, wipkq “ ř jPMi`Aijxjpkq`Bijujpkq˘`we ipkq,(1) where kPN`is the discrete-time index, xiPRqiand uiPRriare, respectively, the state and input vectors of each subsystem iPN, which are assumed to be constrained in convex sets containing the origin XifitxiPRqi:Ax ixiď bx iuand UifituiPRri:Au iuiďbu iu. Moreover, Aii P Rqiˆqi, Bii PRqiˆriare matrices of the proper dimension for each subsystem. The variable wiPRqicomprises neighbor coupling through state vector xjPXjand input vector ujPUj with jPMifitjPNztiu:Aij ‰0_Bij ‰0u, and the external disturbance we i, which is assumed to be bounded by the convex set We i. Therefore, wiis bounded by WifiÀjPMi`AijXj‘Bij Uj˘‘We i, where WiĂRqi. Aggregating all the subsystems states and inputs as xN“ pxiqiPNand uN“ puiqiPN, the global system evolution is xNpk`1q “ ANxNpkq`BNuNpkq`wNpkq,(2) where AN“ rAijsi,jPN, and BN“ rBij si,jPN. Note that only external uncertainties affect the system, i.e. wN“ pwe iqiPNbecause disturbances by coupling are implicitly included in the overall system dynamics. Thus, global additive disturbances are constrained in set WN“ŚiPNWe i. Additionally, global states and inputs are constrained in sets XN“ŚiPNXiand UN“ŚiPNUi, respectively. B. Information exchange and network controllability Each subsystem iPNis managed by a local controller or agent with partial information of the overall system. Thus, agents must communicate and cooperate to attain their control goals. The cooperation network is described by the undirected graph pN,Lq, where the set of subsystems is N, and the set of links is LĎNˆN. Each link lij “ ti, ju “ tj, iu “ lji PL, which connects agents iand j, provides a bidirectional information flow that can be enabled or disabled by the control scheme. Each enabled link involves a fixed cooperation cost clink PRą0. As a result of the trade-off between performance and cooperation costs, the set of active links at instant kdefines the network topology ΛĎL. Considering L“ |L|the number of links, there are 2Ldifferent topologies, which are grouped in a set TfitΛ1,Λ2,...,Λ2Lu. For convenience, Λ1denotes decentralized topology, which corresponds to all cooperation links being disabled, and Λ2Lcorresponds to the centralized topology, i.e. when full network cooperation is established. The disjoint set of cooperation clusters of subsystems, referred to as coalitions, resulting from a given topology 3 ΛPTis denoted by N{ΛfitC1,C2,...,C|N{Λ|u. Coalitions Care non-empty disjoint sets, and satisfy ŤCPN{ΛC“N. Each coalition CPN{Λis described by the discrete-time dynamics xCpk`1q “ ACxCpkq`BCuCpkq`wCpkq, wCpkq “ ř jPMC`ACjxjpkq`BCjujpkq˘`we Cpkq,(3) where xC“ pxiqiPCand uC“ puiqiPCare, respectively, the aggregate state and input vectors of the subsystems in C, which are constrained in the sets XC“ŚiPCXiand UC“ŚiPCUi. The corresponding state matrix is AC“ rAijsi,jPC, and the input matrix is BC“ rBijsi,jPC. The vector wC“ pwiqiPC, which is contained in the set WC“ŚiPCWi, represents the disturbances caused by the coupling with neighbors plus the external disturbances we C. The set of neighbors that affects coalition Cis MCfitjPNzC:ACj‰0_BCj‰0u. From a centralized viewpoint, coalitions dynamics can be aggregated as (2) defining matrices AN“diagpACqCPN{Λ, BN“diagpBCqCPN{Λand the overall state and input vectors as xN“ pxCqCPN{Λand uN“ puCqCPN{Λ, respectively. Assumption 1. Each coalition pAC, BCq, for CPN{Λ, is controllable for any topology ΛPT. Remark 1. The coalitions’ controllability matrices need to be full rank for any topology as a necessary and sufficient network controllability condition [29]. Although it is a conservative assumption, the objective of this work is to evaluate the performance when the controller can optimize the time instant of topology changes. III. CONTROL ALGORITHM Coalitional control balances control performance and cooperation costs adapting the topology in real-time, enabling and disabling links. The price per active cooperation link clink is a tuning parameter in the design of the coalitional controller and plays a significant role in deciding the new topology. Its value needs to be selected properly to represent the cooperation effort. The effect of this parameter was studied in [30]. A. Control goal At time instant k, the goal of each coalition CPN{Λis the minimization of a function JΛ Cp¨q, which sums a stage cost lCp¨q and a cooperation cost gΛ Cp¨q over a prediction horizon Np, and also includes a terminal cost function fΛ Cp¨q, i.e., JΛ C`xCpkq, uCp0 : Np´1|kq,ΛCp0 : Np|kq˘ “fΛ C`xCpNp|kq,ΛCpNp|kq˘` Np´1 ř t“0 lC`xCpt|kq, uCpt|kq˘` Np ř t“0 gΛ C`ΛCpt|kq˘, (4) with fΛ C`xCpNp|kq,ΛCpNp|kq˘“xCpNp|kqJPC`ΛCpNp|kq˘xCpNp|kq, lC`xCpt|kq, uCpt|kq˘“xCpt|kqJQCxCpt|kq`uCpt|kqJRCuCpt|kq, gΛ C`ΛCpt|kq˘“clink|ΛCpt|kq|, where uCp0 : Np´1|kqfiruCp0|kq, . . . , uCpNp´1|kqs and ΛCp0 : Np|kqfirΛCp0|kq,ΛCp1|kq,...,ΛCpNp|kqs are, respectively, the sequences of control actions and topologies in the prediction horizon1;pt|kqindicates the prediction step 1The topology is also relevant for the terminal cost of the controller. For this reason, its sequence is considered until the end of the prediction horizon. tcomputed at instant k;QCě0and RC, PCą0are weighting matrices; the cost per active link clink PRą0is a tuning parameter; |ΛCpt|kq| represents the number of links in coalition Cat prediction step t; and superscript Λdenotes a topology dependency. Note that the size of coalitions can range from an agent C“ tiuto even the global system C“N. Thus, the overall cost can be calculated aggregating (4) as JNp¨q “ řCPN{ΛJΛ Cp¨q. The optimization of (4) has to be performed subject to the coalition dynamics and disturbances (3), constraints XCand UC, and topology ΛĎL. Furthermore, the final state of coalition Cmust be constrained in a terminal region, which is a polyhedron defined as XΛ CffitxCPRqC:AXf CxCďbXf Cu, and whose calculation will be detailed later. Finally, note that since JΛ Cp¨q is convex by definition for any possible topology sequence, we conclude the subsection with some properties of interest of JΛ Cp¨q, which will become useful later. Definition 1 (Np-step stabilizable set [23]).Let XΛ Cfbe a terminal region for the nominal system xCpk`1q “ ACxCpkq` BCuCpkqof coalition CPN{Λ, subject to XCand UC. The Np-step stabilizable set SΛ NppXC,XΛ Cfqis the set of states in XCthat evolves to XΛ Cfin a trajectory of Npsteps by feasible control sequences. Definition 2 (Robust positively invariant (RPI) set).Given topology Λ, the set XΛ Cfis an RPI set for xCpk`1q “ ACxCpkq ` BCuCpkq ` wCpkqif and only if there exists a linear control law uC“KΛ CxCsuch that its evolution satisfies: xCpkq P XΛ CfÝÑ xCpk`1q“pAC`BCKΛ CqxCpkq`wCpkq P XΛ Cf,@wCpkq P WΛ Cand @kPN`. Following [23], it can be proved that the perturbed coalition Cevolves to an RPI set w.r.t. disturbances wC. This fact is summarized in the following lemma. Lemma 1 (Robust Stability [23]). 1) Let the nominal model of coalition be xCpk`1q “ F`xCpkq˘, such that Fp¨q is continuous, and the origin is a fixed point. Let JΛ CpxCqbe a Lyapunov function of the closed-loop system such that aΛ C||xCpkq||σďJΛ C`xCpkq˘ďbΛ C||xCpkq||σ, JΛ C`xCpk`1q˘´JΛ C`xCpkq˘ď ´cΛ C||xC||σ,(5) where aΛ C, bΛ C, cΛ Care positive constants, and σě1. Then, there exists a positive constant ρΛ C“1´ pcΛ C{bΛ Cq ă 1 that satisfies JΛ C`xCpk`1q˘ďρΛ CJΛ C`xCpkq˘. 2) Let JΛ Cp¨q be also Lipschitz in a neighborhood of the origin Ωr CfitxCPXC:JΛ Cp¨q ď rΛ Cuand λΛ Cbe the Lipschitz constant of JΛ Cp¨q in Ωr C, i.e, JΛ CpxCpkq`wCq ď JΛ CpxCpkqq ` λΛ C||wC||. Then, there exists a constant µΛ Cą0that satisfies µΛ Cď`p1´ρΛ Cq{λΛ C˘rΛ C,(6) such that if wCpkq P Bµ CfitwCPRqC:||wCpkq|| ă µΛ Cu for all k, then the perturbed coalition xCpk`1q “ F`xCpkq˘`wCis asymptotically ultimately bounded @xCp0q P Ωr C. 4 As shown in Section V, our goal is to assure that the closedloop system stays asymptotically ultimately bounded when tÑ 8 despite the disturbances by the topology switchings. B. Coalitional MPC algorithm The coalitional method can be prohibitive for large-scale systems since it solves a convex optimization problem for every possible topology. Nevertheless, it can be employed a subset of topologies Tnew ĂTto relieve this issue. Let kPN`denote the time instant, Tup PN`indicate how often the upper control layer is executed, Λcur PTnewpkq be the current topology, and Tnewpkq Ď Tdenote the set of potential successor topologies at instant kaccording to a given criterion, e.g., considering only the topologies whose links differ in one with regard to the current topology. At each time instant k, the proposed hierarchical coalitional MPC algorithm, which is divided into upper and lower control layers, is implemented as detailed in Algorithm 1. Every Tup time instants, the upper control layer is executed, where a centralized controller optimizes Problem 1obtaining the optimal control sequence, the new topology, and the transition horizon to attain the global objective. Algorithm 1 (Upper control layer) Inputs: k, Tup, xNpkq,Λcur Outputs: uNp0 : Np´1|kq,Λnew, Nt 1: The upper control layer measures xC,@CPN{Λ. 2: The upper layer solves the following mixed-integer optimization problem: Problem 1: min uNp0:Np´1|kq,Λnew, Nt JN`xNpkq, uNp0 : Np´1|kq,Λp0 : Np|kq˘, subject to xCpt`1|kq “ ACxCpt|kq`BCuCpt|kq, t “0, . . . , Np´1, xCpkq “ ˜xCpkq,@CPN{Λ, xCpt|kq P XCaWC`ΛCpt|kq˘, t “1, . . . , Np´1, uCpt|kq P UC, t “0, . . . , Np´1, xCpNp|kq P XCf`ΛCpt|kq˘, for all CPN{Λ, and Λpt|kq “ "Λcur if tăNt Λnew if těNt , t “0, . . . , Np, Λnew PTnewpkq, where ˜xCpkqis the coalition state measurement; WCrefers to the disturbances received by coalition C;XCfis the robust invariant terminal region computed offline and imposed as a terminal state constraint (see Section IV for details of its calculation); Λcur and Λnew are the current and the next topology, respectively. 3: The resulting control sequence, the new topology, and the transition horizon are sent to the lower layer. In the lower control layer, each coalition solves Problem 2 to compute its control action in accordance with the last information (i.e., the new topology and transition horizon) received by the upper layer: Algorithm 1 (Lower control layer) Inputs: k, xCpkq,Λcur,Λnew, Nt Outputs: k, Nt, xCpk`1q 1: The state xCis measured for each CPN{Λ. 2: Each coalition Ccalculates its control actions aware of the topology will switch after Nttime instants, as commanded by the upper layer. To this end, it is minimized (4) as: Problem 2: min uCp0:Np´1|kqJΛ C`xCpkq, uCp0 : Np´1|kq,ΛCp0 : Np|kq˘, subject to xCpt`1|kq “ ACxCpt|kq`BCuCpt|kq, t “0, . . . , Np´1, xCpkq “ ˜xCpkq, xCpt|kq P XCaWC`ΛCpt|kq˘, t “1, . . . , Np´1, uCpt|kq P UC, t “0, . . . , Np´1, xCpNp|kq P XCf`ΛCpt|kq˘, and ΛCpt|kq “ "Λcur Cif tăNt Λnew Cif těNt , t “0, . . . , Np, where Λcur Cand Λnew Cdenote coalition Cin the current and next topology, respectively. 3: Each coalition applies the first element of its optimal control sequence uCp0 : Np´1|kqto obtain xCpk`1q @CPN{Λ. 4: Set k“k`1, and update the transition horizon as Nt“Nt´1for the next time instant. Remark 2. Since disturbances are non-manipulable inputs, the feasibility of the state can be endangered. For that reason, a tightened constraint set XCaWCis considered to avoid violating the state constraints. Conversely, the input uCis an optimization variable and, thus, the controller holds it in UC. Remark 3. The transition horizon Ntdoes not appear explicitly in the cost but in the topology constraints of the optimization problems. The problem solved by the upper layer has an extra degree of freedom on choosing the moment to switch topology. Hence, the resulting performance will always be optimal for the MPC receding horizon objective. Remark 4. This basic algorithm may be varied in the way that the transition horizon is determined, e.g., Ntcan be a fixed value, computed in the upper control layer, or obtained in the lower control layer. Note that, in the latter case, feasibility guarantees may be endangered, as detailed later. These alternatives will be assessed in the simulation section. IV. CONTROLLER DESIGN PROCEDURE This section provides requirements to ensure that the closedloop system is feasible and stable with the proposed coalitional MPC scheme. Details regarding the controller design procedure, based on linear matrix inequalities (LMI), are given. A. Feedback design Global sparse matrices KΛand PΛadapted to topology Λ are designed to solve an LMI problem following [17]. These 5 matrices are sufficient to guarantee stability and to find the invariant terminal region, as shown later. Lemma 2 (Coalitional feedback controller design [17]).Let the global discrete system matrices be AN“ rAij si,jPN and BN“ rBijsi,jPN, with the stage cost matrices QN“ diagpQiqiPNand RN“diagpRiqiPN. If there are matrices HΛ“HJ Λ“ rHijsi,jPN, where Hij PRqiˆqj, and YΛ“ rYijsi,jPN, where Yij PRriˆqjwith Hij “0and Yij “0if the link lij PΛis disabled, i.e., if iPCand jRC, such that the following constraint is satisfied » — — — – HΛHΛAJ N`YJ ΛBJ NHΛQ1{2 NYJ ΛR1{2 N ANHΛ`BNYΛHΛ0 0 Q1{2 NHΛ0I0 R1{2 NYΛ0 0 I fi ffi ffi ffi flě0, (7) then there exist feedback control law KΛ“ rKΛ CsCPN{Λ“ YΛH´1 Λthat stabilizes the closed-loop system, and a Lyapunov function fN`xNpkq˘“xNpkqJPΛxNpkq, with PΛ“ rPΛ CsCPN{Λ“H´1 Λthat satisfies fN`xNpkq˘ě 8 ÿ t“0 lN`xNpk`tq, uNpk`tq˘.(8) As proved in [17], if (7) remains true, then the nominal system model xNpk`1q “ ANxNpkq ` BNuNpkqis exponentially stable with the linear control law uNpkq “ KΛxNpkq, with PΛsatisfying xNpkqJPΛxNpkq´xNpk`1qJPΛxNpk`1q ě lN`xNpkq˘, (9) i.e., xNpkqJPΛxNpkqis a a Lyapunov function. Consequently, (8) holds, and the exponential stability of the nominal system is guaranteed. Assumption 2. There exist matrices HΛand YΛsuch that (7) holds for the decentralized topology, which presents the highest disturbances due to couplings, and also provides a feasible solution for any other topology. In order to design the controller, we maximize the trace of HΛsubject to (7) to indirectly minimize PΛ“H´1 Λand, therefore, the control cost-to-go. Remark 5. Feedbacks KΛare defined for all cases ranging from decentralized to centralized topologies. Nevertheless, this method can be implemented with just a subset of topologies Tnew ĂTto avoid combinatorial explosion issues in largescale systems. B. Invariant set design From the coalitions’ viewpoint, the effect of their neighboring agents can be considered as unknown bounded disturbances when calculating invariant set XΛ Cfto decouple the terminal regions and solve a more manageable problem. Assumption 3. For each Cin N{Λ, there exists an RPI set XΛ Cf(whose size depends on the disturbance set WΛ Cthat relies on the coupling with its neighbors in topology Λ) that satisfies pAC`BCKΛ CqXΛ Cf‘WΛ CĎXΛ Cf,XΛ CfĎXC, and KΛ CXΛ CfĎUC under the linear control law uC“KΛ CxC. Hence, the robust positively invariant set for the overall system is computed as XNf“ą CPN{Λ XΛ Cf,(10) which satisfies pAN`BNKΛqXNf‘WNĎXNf,XNfĎXN, and KΛXNfĎUNunder the linear control law uN“KΛxN. Theorem 1. The RPI set XNf“ŚCPN{Λ1pXΛ1 Cfqcorresponding to the decentralized topology Λ1can be used for any other topology. Proof. For topology Λ1, each coalition C“ tiu,@iPN presents the highest disturbances set due to the highest coupling with neighbors, i.e., WΛ CĎWΛ1 C,@ΛPT. Since greater WΛ Cinvolves greater XΛ Cf, it is satisfied that XΛ CfĎXΛ1 Cf. There are numerous procedures to find RPI sets given unknown bounded disturbances [31]. The Multi-Parametric Toolbox (MPT) of MATLAB® [32] can be employed to compute the maximal RPI set for each coalition of topology ΛPTtaking into account its constraints and disturbances. For a conservative approach, the RPI set for decentralized topology can be generally used for any topology because it is computed for the worst-case disturbances. Note that, since an RPI set can be computed for each agent, the advantage of this approach is to avoid combinatorial explosion issues in large-scale systems. V. CONTROLLER FEASIBILITY AND STABILITY Hitherto, it has been proved the stability of the nominal system for a fixed topology and its corresponding feedback in the unconstrained case. In this section, recursive feasibility of the system is ensured when the coalitional MPC is implemented. Due to the persistent bounded disturbances, asymptotic stability cannot be assured. Consequently, we desire the system to converge to a bounded set around the origin [25]. A. Recursive feasibility Assumption 4. At time instant k“0, Problem 1has a feasible control sequence for the centralized topology, i.e., Λcur “Λ2L. Theorem 2. Let Assumption 2, which implies the existence of the most conservative KΛand PΛ(feasible solution for any topology), and Assumption 4 hold. Given topology ΛPTand matrices KΛ,PΛcomputed by maximizing the trace of HΛ subject to the LMI problem (7), then recursive feasibility of Algorithm 1is always ensured. Proof. First, the proof deals with the recursive feasibility for a fixed topology ΛPTfrom the global and coalitions viewpoints. Secondly, it is proved that the system is feasible despite topology switching. Part 1.Consider the state xNpkq P XN@k, the predicted nominal state is xNp1|kq “ ANxNpkq`BNuNpkq. We can prove that xNp1|kqfulfills constraints despite disturbances. Since xNpk`1q “ xNp1|kq`wNpkq, and due to the problem constraints: xNp1|kq P XNaWNand wNpkq P WN, then 6 xNpk`1q P XN. Thus, the recursive feasibility of the state is guaranteed. Under Assumption 4, at least a topology Λhas an overall feasible control solution at instant k“0: uNp0 : Np´1|kq“ruNp0|kq, . . . , uNpNp´1|kqs,(11) which fulfills constraints xNpkq P XN,uNpkq P UN, and xNpNp|kq P XNffor all disturbance sequences. At next instant k`1, a shifted feasible sequence us Ncan be formed as the tail of (11) as us Np0 : Np´1|k`1q “ rus Np0|k`1q, . . . , us NpNp´1|k`1qs “ “ ruNp1|kq, . . . , uNpNp´1|kq, KΛxNpNp|kqs. Since the overall feasible control solution uNp0 : Np´1|kq fulfills the terminal state constraint, the term xNpNp|kq P XNf, and, at k`1it holds pAN`BNKΛqxNpNp|kq`wNPXNf,@wNPWN. It follows from applying this procedure recursively that the feasibility of Problem 1is guaranteed for topology Λand @k. Note that KΛand XNfcan be respectively dis-aggregated in KΛ Cand XΛ Cffor each CPN{Λ. Thus, there are also feasible sequences for any CPN{Λat instant k“Tup, i.e., uCp0 : Np´1|Tupq“ruCp0|Tupq, uCp1|Tupq, . . . , uCpNp´1|Tupqs, subject to xCpTupq P XC,uCpTupq P UC, and xCpNp|Tupq P XΛ Cffor all disturbance sequences. At the following instant, a shifted feasible sequence us Ccan be composed with the feedback gain KΛ Cas us Cp0 : Np´1|Tup `1q “ ruCp1|Tupq, . . . , uCpNp´1|Tupq, KΛ CxCpNp|Tupqs. Given that xCpNp|Tupq P XΛ Cf, at the next instant, it holds pAC`BCKΛ CqxCpNp|Tupq`wCPXΛ Cf,@wCPWΛ C. Therefore, there are feasible solutions available for Problem 2 at time instant Tup `1. Recall that if the recursive feasibility from the centralized viewpoint is assured, it can also be guaranteed from the coalitional viewpoint because KΛis aggregated from the KΛ Cfor any CPN{Λand XNfis computed by (10). Applying this procedure recursively, it can be proved that coalitional feasibility holds for all k. Part 2.Switching between different topologies requires that the overall feasibility is not endangered. First, note that if a change of topology occurs at instant kc, a feasible global input sequence must exist according to the constraints of Problem 1. Otherwise, the same topology Λcur could be used because it is a feasible control sequence, as proved in Part 1. Secondly, if global feasibility is not endangered when switching to topology Λnew, neither is the coalitional one endangered. In particular, global feasibility implicitly considers the dynamics and constraints of all coalitions CPN{Λnew. Therefore, feasibility is also guaranteed despite topology changes. B. Robust stability Here it is analyzed the asymptotic ultimate boundedness of the constrained linear system in the presence of bounded additive uncertainties for our proposed coalitional MPC scheme. In particular, we consider the worst-case coupling between subsystems. By using this conservative approach, we can guarantee ultimate boundedness for any topology. Assumption 5. The most disturbed case is bounded by the constant ¯µNď p1´¯ρNq{¯ λN¯rN, such that wipkq P Bµ ifi twiPRqi:||wipkq|| ă ¯µNu@k, for any iPNand ΛP T. Here, ¯ρN,¯ λN,¯rNare calculated according the Lemma 1 based on a common upper bound of the cost function ¯ JNp¨q for all system topologies, which is defined as ¯ JNp¨q “ lNp¨q`¯gNpΛ2Lq` ¯ fNp¨q,(12) where lNp¨q is the stage cost, ¯gNpΛ2Lqis the cooperation cost of the centralized topology (all cooperation links enabled), and ¯ fNp¨q is the terminal cost, whose weighting matrix ¯ Pą PΛ@Λis computed such that (7) is satisfied for any KΛ.2 Note that, although Assumption 5 is conservative, it is needed to ensure robust stability and feasibility. Theorem 3 (Asymptotic ultimate boundedness).Let xNpk` 1q “ ANxNpkq ` BNuNpkqbe the global model of the nominal system, subject to xNpkq P XN,uNpkq P UN, and let the linear control law uNpkq “ KΛxNpkqstabilize the nominal system exponentially in SNppXN,XNfq. Considering that its cost function JNp¨qis Lipschitz in SNppXN,XNfq, and Assumption 5 holds, then system (2) with the linear control law uNpkq “ KΛxNpkqis asymptotically ultimately bounded @xNp0q P ¯ Ωr N, with ¯ Ωr NfitxNPXN:¯ JNp¨q ď ¯rNu Ď SNppXN,XNfq Ď XN. Proof. Firstly, given a fixed topology ΛPT, it is proved exponential stability and asymptotic ultimate boundedness for the nominal and the uncertain system, respectively. Secondly, we prove that the system is asymptotically ultimately bounded despite the switching of topologies. Part 1.With respect to Lemma 1, the overall nominal system is continuous, and we have to prove that ¯ JNp¨q is Lyapunov function for the nominal system, i.e., a non-increasing MPC cost. This upper-bound cost function decreases in time if ¯ JN`xNpk`1q, us Np0 : Np´1|k`1q,Λp0 : Np|k`1q˘ď ¯ JN`xNpkq, uNp0 : Np´1|kq,Λp0 : Np|kq˘ holds, which can be rewritten as Np´1 ř t“0 lN`xNpt|k`1q, uNpt|k`1q˘` Np ř t“0 ¯gN`|Λ2Lpt|k`1q|˘`¯ fN`xNpNp|k`1q˘ď Np´1 ř t“0 lN`xNpt|kq, uNpt|kq˘` Np ř t“0 ¯gN`|Λ2Lpt|kq|˘`¯ fN`xNpNp|kq˘. Since in the nominal case xNpt`1|kq “ xNpt|k`1q, deleting terms in common and rearranging terms, we have lN`xNpNp|kq, uNpNp|kq˘`¯ fN`xNpNp|k`1q˘´¯ fN`xNpNp|kq˘ ďlN`xNp0|kq, uNp0|kq˘, 2It is straightforward to prove that such weighting matrix exists. A very simple and conservative choice would be to take ¯ P“řΛPTPΛ. 7 where the left-hand side is lower than or equal to zero according to (9). Likewise, the right-hand side is positive definite by construction. Since ¯ JNp¨q is positive definite by (12), it is proved that ¯ JNp¨q is a Lyapunov function and fulfill the condition of Lemma 1. Moreover, since ¯ JNp¨q is a quadratic function with QNě0and RN,¯ Pą0, and the system is linear, we claim that Problem 1ensures exponential stability in SNppXN,XNfqfor the nominal model system (see Theorem 2.24 (c) in [33]). Nevertheless, when additive disturbances are considered, exponential stability may be endangered. Regarding Lemma 1, asymptotic ultimate boundedness can be assured for system (2), which has bounded additive uncertainties, if ¯ JNp¨q is also Lipschitz in SNppXN,XNfq. Since the system is linear and continuous, there exists a Lipschitz constant ¯ λNof ¯ JNp¨q in SNppXN,XNfqsuch that ¯ JNpxNpkq`wNq ď ¯ JN`xNpkq˘` ¯ λN||wN||. Since previous conditions and Assumption 5 hold, the overall perturbed system (2) is asymptotically ultimately bounded for any xNp0q P ¯ Ωr NĎSNppXN,XNfq Ď XN. Likewise, coalitions states are also ultimately bounded while satisfying constraints of Problem 2. Part 2.When a change from Λcur to Λnew is performed, the disturbances that affect the subsystems change. Since the asymptotic ultimate boundedness of the overall system and coalitions are guaranteed for the most disturbed case and a common Lyapunov function, we conclude that the closed-loop system remains asymptotically ultimately bounded despite the topology switchings. VI. SIMULATION This section presents the eight-coupled tanks plant and describes the techniques used to assess the proposed algorithm. ?3·qb 1 2 43 5 6 7 8 h1h2 h4 h3 h5 h7h8 h6 qaqbqcqd (1- ?1)·qa (1-?2-?3)·qb (1-?5)·qd (1-?4-?6)·qc ?6·qc ?1·qa?2·qb?4·qc?5·qd Subsystem 1 Subsystem 2 Subsystem 3 Subsystem 4 Fig. 2. Schematic diagram of the eight-coupled tanks plant. A. Plant description The plant is formed by eight interconnected tanks, as shown in Fig. 2. The four upper tanks (#3,#4,#7, and #8) discharge flows to the lower ones (#1,#2,#5, and #6), and these, in turn into sinking tanks. Moreover, the third and the fourth tanks are connected through a pipe, and also the fifth and the sixth tanks. The plant is controlled by four pumps that keep the water circulation between tanks. Additionally, six three-way valves, whose task is to divide into two ways the flow on arrival, are provided. The overall system is divided into N“4subsystems: the first one is composed of tanks #1 and #3; the second one by tanks #2 and #4; the third one by tanks #5 and #7; and the last one by tanks #6 and #8. The control objective is that the lower tanks reach target levels taking into account control and cooperation costs subject to operational constraints. Consequently, it is considered a multi-variable control problem with four inputs (qa,qb,qc,qd) and four outputs (h1,h2,h5,h6). B. Plant model The model used in the simulation is non-linear and can be obtained applying mass balance and Bernoulli’s law. It is given by the following equations: S1dh1 dt “a3?2gh3´a1?2gh1`γ1qa 3600 , S2dh2 dt “a4?2gh4´a2?2gh2`γ2qb 3600 , S3dh3 dt “ ´a3?2gh3´a34a2gph3´h4q`p1´γ2´γ3qqb 3600 , S4dh4 dt “ ´a4?2gh4`a34a2gph3´h4q`p1´γ1qqa 3600 `γ6qc 3600 , S5dh5 dt “a7?2gh7´a5?2gh5`a56a2gph6´h5q`γ4qc 3600 , S6dh6 dt “a8?2gh8´a6?2gh6´a56a2gph6´h5q`γ5qd 3600 , S7dh7 dt “ ´a7?2gh7`p1´γ5qqd 3600 `γ3qb 3600 , S8dh8 dt “ ´a8?2gh8`p1´γ4´γ6qqc 3600 , (13) where hnis the water level of tank nP t1,2,...,8u;Sn“ 13.89 ¨10´3m2is its corresponding cross-section, which is assumed to be equal for all tanks; an“50.265 ¨10´6m2 stands for the cross-section of the outlet pipes; and a34, a56 “ 50.265 ¨10´7m2refers to the cross-section of pipes which are connecting tanks 3-4and 5-6. The parameter γiP r0,1s, with iP t1,2,...,6u, refers to the opening of the six threeway valves (γ1, γ4“0.3;γ2, γ5“0.4and γ3, γ6“0.1), g“9.81 m{s2corresponds to the gravity, and qmrefers to the caudal flow provided by pump mP ta,b,c,du. In order to obtain the control model, let us define an operating point of each tank, measured in meters, by h0 1“2.7831; h0 2“4.8066;h0 3“1.3388;h0 4“1.7346;h0 5“2.5740; h0 6“2.6479;h0 7“1.0317;h0 8“1.2613, and the operating point of each pump, measured in cubic meters per hour, by q0 1“1.2975;q0 2“1.6663;q0 3“1.4253;q0 4“1.0113. Hence, the discrete linear state-space model can be expressed as ¯xNpk`1q “ AN¯xNpkq`BN¯uNpkq` ˆ dNpkq,(14) where ¯xN“ rh1pkq´h0 1, . . . , h8pkq´h0 8sJbecomes the state vector, and ¯uN“ rqapkq´q0 a, . . . , qdpkq´q0 dsJis the input vector. Note that the system is linearized around the operating point. We have striven to minimize the linearization error using offset free concept proposed in [34], where it is estimated a disturbance term ˆ dNpkq “ xNpkq ´ ¯xNpkqat each time instant k. Note that the system-modeling error ˆ dNcan also be considered as an external disturbance wN“ˆ dN. 8 Likewise, the state evolution of each subsystem iPN“ t1,2,3,4uis given by ¯xipk`1q “ Aii ¯xipkq`Bii ¯uipkq` ¯wipkq` ˆ dipkq, ¯wi“řjPMi`Aij ¯xjpkq`Bij ¯ujpkq˘,(15) where ¯wirefers to the coupling through states and inputs with its neighbors jPMi, which is constrained in the convex set WifiÀjPMi`AijXj‘Bij Uj˘, and the subsystem-modeling error is defined as ˆ dipkq “ xipkq´ ¯xipkq,@k. The subsystem matrices are given by: A11 “„0.9763 0.0336 0 0.9659, A12 “„0 0 0 0.0062, B11 “„0.0296 0.002 , B12 “„0.0008 0.0491, A22 “„0.9819 0.0297 0 0.9700, A21 “„0 0 0´0.0062, B21 “„0.0010 0.0689, B22 “„0.0396 0, B23 “„0.0001 0.0098, A33 “„0.9752 0.0382 0 0.9613, A34 “„0.0144 0 0 0, B32 “„0.0002 0.0098, B33 “„0.0296 0, B34 “„0.0014 0.0588, A44 “„0.9756 0.0346 0 0.9642, A43 “„´0.0144 0 0 0, B43 “„0.0008 0.0589, B44 “„0.0395 0, Additionally, states and inputs are constrained by 0.02 m ăhnď5 m,@nP t1,2,...,8u, 0 m3{hăqmď5 m3{h,@mP ta,b,c,du. C. Control methods Each agent ican cooperate through a network formed by the edges L“ t1,2u,t2,3u,t3,4u(, which can be enabled or disabled. Since the number of links is L“3, there are eight cooperation topologies T“ tΛ1,Λ2,...,Λ8u, which are determined by the set of active links as displayed in Fig. 3. For example, all disabled links L“ t0,0,0udefine the decentralized topology Λ1, and all active links L“ t1,1,1u refer to centralized topology Λ8. Ʌ8 {1, 1, 1} Ʌ7 {1, 1, 0} Ʌ6 {1, 0, 1} Ʌ4 {0, 1, 1} Ʌ5 {1, 0, 0} Ʌ3 {0, 1, 0} Ʌ2 {0, 0, 1} Ʌ5 {1, 0, 0} Ʌ2 {0, 0, 1} Ʌ3 {0, 1, 0} Ʌ1 {0, 0, 0} 3 links2 links1 link 0 links Fig. 3. The partially ordered set containing the eight cooperation topologies. In this work, the set of potential successor topologies Tnew ĎTis formed by those topologies that only differ in zero or one link with respect to the current topology Λcur, i.e., those ΛPTwhose Hamming distance is 0or 1, as follows Tnew fitΛĂΛcur :|Λ|“|Λcur|´1u YtΛcur ĂΛ : |Λ|“|Λcur|`1uYtΛcuru. The following five MPC control methods are assessed: 1) CEN consists of a centralized MPC scheme characterized by full system information. 2) ABRUPT presents a coalitional MPC algorithm where the topology Λcur is instantly switched to Λnew PTnew, i.e., Nt“0, which corresponds to schemes in [15], [35]. In this method, the topology Λ˚ new is chosen every Tup “ 20 samples after evaluating by exhaustive search |Tnew| quadratic problems. 3) PRED10 implements a coalitional MPC scheme where the transition to Λnew is made throughout a fixed transition horizon Nt“10. Every Tup “20 time steps from k“10, the upper control layer assesses |Tnew|problems and selects the optimal topology Λ˚ new to make transitions at the same time that ABRUPT. 4) PREDNtup is a coalitional MPC algorithm that every Tup “20 samples solves |Tnew|¨|Nt|convex problems and chooses that with the lowest cost to find the pair tΛ˚ new, N˚ tu. The variable NtPN`takes values between 0-10, and the considered number of topologies is |Tnew|. 5) PREDNtlow presents a coalitional MPC scheme where the transition horizon Ntis defined as an optimization variable of the lower control layer. The variable Ntcan take values between 0and 10. Note that the new topology Λ˚ new is still chosen by the upper layer every Tup “20 time steps. Table I displays a comparative summary of these methods. Remark 6. Feasibility guarantees may be endangered in PREDNtlow since the lower layer can delay the instant in which the new topology, chosen by the upper layer, is implemented. However, we did not have any feasibility problem in the implementation of the PREDNtlow method. The weighting Qiand Riused for the controller design are Qi“„1 0 0 0, Ri“0.1,@iPN“ t1,2,3,4u. The global feedback gain KΛ“ rKisiPNand global weighting matrix PΛ“ rPisiPN, which defines the terminal cost, are obtained according to Lemma 2. As explained before, to compute the RPI set XNf“ŚiPNXif, standard methods [36] have been implemented with the Multi-Parametric Toolbox (MPT) of MATLAB® [32]. To reduce complexity, the RPI set XNfemployed is the one of the decentralized topology Λ1 because it is the most conservative (see Theorem 1). Fig. 4 shows its corresponding invariant sets Xif,@iP t1,2,3,4u. D. Results The sampling time used in the simulation is Tm“5s, with a simulation length of Nsim “240 time steps. The initial topology is the centralized one Λ8. Recall that the topology will gradually change according Fig. 3. The MPC methods consider a prediction horizon Np“20 samples, and a cost per active cooperation link clink “2¨10´3. N.B.: This value is chosen large enough to represent cooperation costs, considering that a small value and a huge value will always lead to a centralized and a decentralized network, respectively. Further details about the effect of the tuning