Robust coalitional model predictive control with negotiation of mutual interactions

Journal of Process Control 123 (2023) 64–75 Contents lists available at ScienceDirect Journal of Process Control journal homepage: www.elsevier.com/locate/jprocont Robust coalitional model predictive control with negotiation of mutual interactions A. Sánchez-Amores ∗, P. Chanfreut, J.M. Maestre, E.F. Camacho Department of Systems and Automation Engineering, University of Seville, Camino de los Descubrimientos, no number E-41092, Seville, Spain article info Article history: Received 28 October 2022 Received in revised form 24 January 2023 Accepted 30 January 2023 Available online xxxx Keywords: Model predictive control Coalitional control Distributed control Robust control Constrained systems abstract This article presents a robust coalitional model predictive control (MPC) approach where neighboring agents negotiate the bounds of their coupling variables. Also, the control variables of each agent are divided into a private part that is locally optimized, and a public part that is controlled by the corresponding neighbors. Under certain conditions, the agents communicate to update the bounds that determine the constraints on these private and public variables. Moreover, the mutual disturbances induced by coupling are considered using a tube-based approach, guaranteeing recursive feasibility and stability of the closed-loop system. The proposed method is tested on a simulated eight-input coupled tank benchmark to illustrate its benefits. ©2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 1. Introduction Model predictive control (MPC) is a computer-based control method that uses a model to predict the evolution of a system as a function of the sequence of inputs provided along a given horizon. Then, the minimization of a cost function allows us to determine an optimal sequence to steer the system according to the control designer’s goals. To this end, the MPC controller solves an optimization problem at each time step, which may include explicit constraints on the system variables, among other complex issues such as delays, information about expected disturbances, and uncertainties affecting the system. Once a minimizer has been obtained for the optimal control problem, the element of the sequence corresponding to the current time step is implemented and the problem is solved again at the next time instant following a receding horizon approach. One of the principal drawbacks of MPC is that it cannot be applied to large-scale systems in a centralized way due to the number of decision variables in the optimization problem. Naturally distributed systems such as smart grids and water networks require distributed MPC (DMPC) strategies [1–4]. The main idea is to decompose the overall system into subsystems, which are assigned to local interacting controllers, also known as agents. As a consequence, the global problem is partitioned and solved in a distributed manner, preserving the essence of MPC while providing increased scalability and flexibility [5,6]. In the DMPC ∗Corresponding author. E-mail addresses: [email protected] (A. Sánchez-Amores), [email protected] (P. Chanfreut), [email protected] (J.M. Maestre), [email protected] (E.F. Camacho). framework, the communication and cooperation between the set of agents in their distributed computations play a key role to optimize the global performance and provide theoretical guarantees. While complete decentralized strategies minimize communication demands, fully cooperative DMPC algorithms can attain centralized-like behavior. Nonetheless, it may be interesting to restrict the information exchange and minimize the cooperation effort when the latter does not compromise the overall performance. A relevant issue in this regard is that of the coupling that often exists between the subsystems’ dynamics. In this regard, authors in [7] apply a distributed tube-based MPC approach where local controllers share information about the size of their currently used constraint space to reduce conservatism by exploiting the varying degree of coupling. Within the DMPC framework, the so-called coalitional control analyzes the dynamic couplings between the different parts of the system so that only strongly coupled elements communicate with each other. In this respect, the communication topology can be chosen based on a function that penalizes the communication cost [8,9] or the coupling strength [10,11]. Consequently, the global control structure is adjusted in real-time, optimizing the existing computational and communication resources. The underlying idea is to dynamically cluster the local agents into cooperative groups, i.e., coalitions, whenever this improves the overall system’s performance [12], thus leading to timevarying partitions. Accordingly, there may be coalitions of agents exchanging information, while others may be operating in a decentralized fashion. This idea of dynamic clusters of controllers has drawn attention within the context of multi-agent networked systems [13,14]. In particular, coalitional control is presented as an alternative that is halfway between fully cooperative schemes https://doi.org/10.1016/j.jprocont.2023.01.017 0959-1524/©2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). A. Sánchez-Amores, P. Chanfreut, J.M. Maestre et al. Journal of Process Control 123 (2023) 64–75 and decentralized control. Note that in distributed systems, the lack of coordination may compromise the global performance [5]. Likewise, coordinating all control decisions increases the communication and computation demands [15]. Considering this, coalitional schemes promotes data sharing and coordination under certain conditions. See also [16], which reviews the literature on control by clustering strategies, discussing different criteria for partitioning the system and providing examples of application. In this article, we present a novel coalitional MPC approach in which input variables are decomposed into private and public versions sharing a common constraint space. The proposed control scheme shares the nature of coalitional MPC strategies introduced in [17–19], which permit partial modes of cooperation, ranging from the fully distributed MPC algorithm, where all agents share information, to agents working in a decentralized manner without any communication. In this way, there may be a group of agents operating cooperatively for the sake of overall performance, while others work independently without exchanging data. Finally, note that this article is closely related to [20], which is used as the starting point for the problem and the algorithm formulation. However, the previously referred work lacks robustness and formal stability guarantees because it is based on a stochastic approach. To remedy this, we employ here a tube-based framework, which leads to significant changes in the control strategy. In particular, we follow a tube-based MPC approach applied in a distributed manner, where disturbances induced by the coupling between subsystems are modeled as bounded disturbances. Finally, we proceed to summarize the main contributions of this paper: •First, this work strengthens the control algorithm presented in [20] by using a robust approach to guarantee constraint satisfaction, recursive feasibility, and stability. The main idea of the algorithm is to decompose coupling so that agents can cede a portion of their inputs to neighboring agents. •Second, we propose a coalitional scheme where the agents’ communication is event-based, and they negotiate the constraint space of the shared inputs by performing an iterative procedure based on dual decomposition. Unlike traditional distributed algorithms, where state or input sequences are communicated, here the agents only broadcast information about the scale factors that bound their inputs’ constraint space. In this way, we minimize the amount of data to be shared. •Third, our approach promotes communication between agents only when it generates a significant benefit in terms of overall performance. In this way, it is possible to reduce the communication and computational burden without significantly compromising the system performance with respect to the centralized behavior. The outline of the rest of the article is organized as follows. In Section 2, the model of the system is presented, introducing the decomposition of input variables into private and public parts, and defining the setting for coalitional MPC. In Section 3, we describe the application of tube-based MPC to deal with uncertainties. Furthermore, we distribute the problem using the dual decomposition algorithm to negotiate the bounds of coupling variables. Section 4presents the proposed control scheme. In Section 5, conditions for recursive feasibility and stability are given. Section 6shows simulation results on an eight input-coupled tank system. Finally, concluding remarks are provided in Section 7. 2. System description Consider a discrete-time, linear time-invariant system that can be divided into a set N= {1,2,...,N}of input-coupled subsystems modeled as x+ i=Aiixi+Biiui+di, with di=∑ j∈Ni Bijuj,(1) where xi∈Rnxiand u∈Rnuiare the state and control input of subsystem i∈N, and x+ iis the state at the next time instant. Accordingly, di∈Rnxirepresents the input coupling, with the set Niof the neighboring agents defined as Ni≜{j∈N:Bij = 0,j= i}. Moreover, each subsystem i∈Nis subject to local constraints in its state, i.e., xi∈Xi, and input, i.e., ui∈Ui, where Xiand Ui are convex sets containing the origin. By aggregating all subsystems’ states x=[xi]i∈N∈Rnxand inputs u=[ui]i∈N∈Rnu, we can describe the global behavior of the system as x+=Ax +Bu,(2) where matrices A=diag(Aii)i∈Nand B=[Bij]i,j∈Nrepresent the global model and are defined as the aggregation of (1) for all subsystems. Note that the centralized model (2) does not include uncertainties because interactions are already present in the global matrices. 2.1. Decomposition of coupling variables In this work, we apply the decomposition method proposed in [20], where input variables are partitioned into public and private parts. Consequently, we can decompose a local variable uias ui=upr i+∑ j∈Mi upu ij ,(3) where (i) upr iis the private part of the variable, which is controlled exclusively by the agent that owns it, i.e., i, and it must verify upr i∈αiUiwith αi∈ [0,1]. (ii) upu ij is the public part of uicontrolled by agent jand must verify upu ij ∈αijUi. The set of agents jthat can manipulate the public part of uiare defined as affected subsystems Mi= {j∈N:Bji = 0,j= i}. Remark 1. In general, the sets Niand Miare different and depend on the dynamics of the system, i.e., while Nicontains the set of neighbors that affect i,Midefines the set of agents affected by the input of subsystem i. Remark 2. Without loss of generality, in this article we consider input-coupled subsystems as defined in (1), i.e., Aij =0 for all i= j. However, this decomposition could similarly be applied to state-coupled subsystems. In particular, each agent iwill locally control the private part of its input variable, upr iand the public part of its neighboring inputs {upu ji }∀j∈Ni. Accordingly, since the dynamics of the subsystem iis affected by the inputs ujof neighboring subsystems, agents j∈Ni will decide the value of the private input upr j. Moreover, the set j∈Miof agents affected by subsystem’s iinput, will decide the public part of ui, i.e., {upu ij }∀j∈Ni. Therefore, variables that cannot be locally controlled such as {upr j}∀j∈Niand {upu ij }∀j∈Miare treated as bounded disturbances from subsystem’s iviewpoint. Fig. 1 shows a diagram with three agents illustrating the decomposition and negotiation that occurs in our method. 65 A. Sánchez-Amores, P. Chanfreut, J.M. Maestre et al. Journal of Process Control 123 (2023) 64–75 Fig. 1. Scheme of the proposed variable decomposition for three agents. Each agent imanipulates its private variable upr iand the public input upu ji ceded by its neighbor j∈Ni. Agents can negotiate the value of input scale factors αi, αji. Focusing on agent 1, we identify N1= {3}. Therefore, agent 1 can manipulate a public part of u3. Also, M1= {2}, meaning that subsystem 2 is affected by u1and will manipulate a public part of it. Hereafter, it is considered that agents should determine the values of upr iand upu ij ,∀i∈Nand j∈Mi, so that the next inequality is satisfied: αi+∑ j∈Mi αij ≤1,∀i∈N.(4) Note that if uiis calculated according to (3) and (4) holds, then the input constraint ui∈Uiwill be satisfied for all i∈N. 2.2. Control architecture and strategy Let us assume that a local controller or agent governs each subsystem i∈N. Also, consider that local agents are interconnected through a configurable data network that allows communication among them. We can describe this network using the graph G=(N,L), being Nthe set of nodes that represent the agents, and Lthe set of links, with L⊆LN= {{i,j}|i,j∈N}. It is considered that the state of the links can change dynamically to suit the changing requirements of the control scheme, so that the set of agents is partitioned into disjoint communication elements or coalitions, i.e., only agents inside a coalition can exchange data. In particular, let ˆαbe a threshold to enable or disable the agents’ communication. Concerning this, the state of the links Lwill change depending on the values of αij with respect ˆα. Let P(k)= {C1,C2,...,C|P(k)|}define the partition of the system at time instant k, where Ci⊆Ndescribes the ith coalition within the partition, with Ci∩Cj= ∅ for any i= j. The size of Ci can range between the following two extremes: (i) If all agents work in a decentralized manner, i.e., there is no communication among them, there will be |N|singletons, i.e., P(k)= {{1},{2},...,{N}} and |Ci| = 1,∀i. This happens when αij <ˆα, for all upu ij ∈αijUi, with i∈Nand j∈Mi. (ii) If all agents work in a centralized manner, they will form the grand coalition, so there will be a single coalition that will group all agents, i.e. P(k)= {N}. This happens when αij ≥ ˆα, for all upu ij ∈αijUi, with i∈Nand j∈Mi. As coefficients αij take values over and below ˆαfor different combinations of {i,j}, there will be clusters of agents that communicate, i.e., coalitions, whereas others may be working in a decentralized manner. Therefore, agents operate in a flexible fashion that allows partial modes of cooperation. The agents’ goal is to maximize global performance with a minimum exchange of information. In particular, agents aim to minimize the cost function ∞ ∑ t=0∑ i∈N(x⊤ i(t+1) Qixi(t+1)+ upr i ⊤(t)Rpr iupr i(t))+∑ j∈Ni upu ji ⊤(t)Rpu iupu ji (t), (5) with Qi,Rpr i,Rpu ibeing positive definite weighting matrices. Any agent iwill optimize the private part of its input variable upr i, the public part of neighboring input variables upu ji ,∀j∈Ni, and will be able to negotiate scale factors αiand αji. In contrast to standard distributed schemes, which negotiate the values of coupling variables, we propose a negotiation of the scale factors that bound these variables. Also, the use of public variables will be penalized with a higher cost to discourage unnecessary cooperation efforts, i.e., Rpu i≫Rpr i. 3. Tube-based MPC for regulation In contrast to [20], where uncertainties were considered using a scenario-based MPC, we propose a tube-based MPC approach to deal with neighboring uncertainties. In this way, we are able to satisfy the system’s constraints regardless of the realization of the disturbance. 3.1. Nominal control problem Considering the definition of the local state (1), and the decomposition of private and public inputs (3), we can rewrite the dynamics of subsystem ias: x+ i=Aiixi+Biiupr i+∑ j∈Ni Bijupu ji +wi, with wi≜∑ j∈Mi Biiupu ij +∑ j∈Ni Bijupr j. (6) As a result of the interaction between subsystems, the disturbance wigroups the uncertainties generated by the agents in Mi, which control the public part of ui, i.e., {upu ij }∀j∈Mi, and those 66 A. Sánchez-Amores, P. Chanfreut, J.M. Maestre et al. Journal of Process Control 123 (2023) 64–75 Table 1 Notation summary. uiLocal variable of subsystem i, decomposed as (3). Xi,UiState and input set of constraints of xiand ui. upr iPrivate part of uicontrolled exclusively by agent i. upu ij Public part of uiceded to agents j∈Mi. upr jPrivate part of ujcontrolled by agent j∈Ni. upu ji Public part of ujwith j∈Ni, controlled by agent i. αiTightening factor for set Ui, with upr i:∈ αiUi. αji Tightening factor for set Uj, with upu ji :∈ αjiUj. NiNeighboring subsystems of i:{j∈N:Bij =0,j= i}. MiSubsystems affected by ui:{j∈N:Bji =0,j= i}. WiBound of agent i’s disturbances: wi∈Wi(7). RiRobust positively invariant set of subsystem i. xi,uiNominal state and input local variables of i(8). Xi,UiState and input set of nominal constraints. in Ni, which control the private part of their input variables, i.e., {upr j}∀j∈Ni. Therefore, uncertainties are bounded by: wi∈Wi≜⨁ j∈Mi BiiWpu ij ⊕⨁ j∈Ni BijWpr j, with Wpu ij =αijUi,and Wpr j=αjUj, (7) where Wpu ij and Wpr jare convex closed sets containing the origin, satisfying upu ij ∈Wpu ij and upr j∈Wpr j. To take into account these uncertainties, we will follow a robust approach. The tube-based MPC approach [21,22], is characterized by solving an MPC problem for the nominal system, that is, without considering disturbances, and by adding an auxiliary control law to keep the evolution of the real system within a tube around the nominal trajectory. Based on (6), we can derive the nominal model for each subsystem iby ignoring interactions wi, i.e., x+ i=Aiixi+Biiupr i+∑ j∈Ni Bijupu ji .(8) In the nominal control problem, subsystem iignores private variables of neighboring agents Ni, and the public part of its input variable ui, which is controlled by agents belonging to Mi. For the sake of clarity, Table 1 has been included to summarize the notation used in this article. Let us aggregate local input-to-state matrices as Bi≜[Bii Bij] for all j∈Ni. In what follows, we consider the following assumptions: Assumption 1. For every subsystem i: •There exists a local feedback gain Kithat guarantees that AKi≜(Aii +BiKi) is stable. •It is possible to find a robust positively invariant (RPI) set Ri that satisfies: AKiRi⊕Wi⊆Ri,Ri⊆Xi, and KiRi⊆Ui. •Nominal state and input constraints are non-empty sets: Xi= ∅ and Ui= ∅, with Xi≜Xi⊖Riand Ui≜Ui⊖KiRi. Remark 3. Assumption 1 is generally considered in the tubebased MPC framework [21,22], as it must be satisfied to have a non-empty solution space for the nominal problem. For the nominal problem, the goal of each local controller iis to minimize the objective function Jci(xi,upr i,[upu ji ]∀j∈Ni,αi)= k+Np−1 ∑ t=k ℓi(xi(t),upr i(t),[upu ji (t)]j∈Ni) +fi(αi)+Vf i(xi(k+Np)), (9) where tdenotes the time step along the prediction horizon Np. Let us describe upr iand upu ji as the sequence of inputs upr i(·) and upu ji (·) from t=kto k+Np−1: upr i=[upr i(0)⊤,upr i(1)⊤,..., upr i(Np−1)⊤]⊤, upu ji =[upu ji (0)⊤,upu ji (1)⊤,..., upu ji (Np−1)⊤]⊤. (10) Let αigroup the value of the scale factors controlled by i: αi=[αi,[αji]∀j∈Ni]⊤ .(11) Note that scale factors αiand αji are recalculated every time instant k, but are kept constant over the prediction horizon. Furthermore, the first term of (9) corresponds to the stage cost ℓi(·), described by the quadratic function: ℓi(xi(t),upr i(t),[upu ji (t)]j∈Ni)= x+ i ⊤Qix+ i+upr i ⊤Rpr iupr i+∑ j∈Ni upu ji ⊤Rpu iupu ji .(12) Moreover, fi(·) is introduced as a penalization for the scale factors αiand αji for all j∈Ni, i.e., fi(αi)=ρprαi+∑ j∈Ni ρpuαji,(13) being ρpr and ρpu positive weighting factors. The last term of Jci represents the terminal cost function, with Pi>0: Vf i(xi(k+Np))=xi(k+Np)⊤Pixi(k+Np).(14) Considering (9), the centralized MPC problem for the nominal system at each time instant kis defined as: min [upr i,[upu ji ]∀j∈Ni ,αi]i∈N ∑ i∈N Jci(xi,upr i,[upu ji ]∀j∈Ni,αi) s.t. xi(k)∈xi⊕(−Ri), xi(t+1) =Aiixi(t)+Biiupr i(t)+∑ j∈Ni Bijupu ji (t), upr i(t)∈αiUi, upu ji (t)∈αjiUj, xi(t)∈Xi, xi(Np)∈Xf i= {0}, ∀i∈N,∀j∈Ni, ∀t=k,...,k+Np−1. (15) where Xf irepresents the terminal region of the nominal model, which has been chosen as the origin. In this way, at the end of the prediction horizon, the state of the plant subject to uncertainties will stay in a neighborhood of the origin given by Ri. Moreover, let upr,∗ iand upu,∗ ji be the optimal sequences of private and public inputs, which are solution for subsystem ito problem (15). Note that there will be |Ni|different sequences upu,∗ ji , one from each neighboring subsystem. Since the MPC follows a receding horizon strategy, only the first component of the optimized sequences upr,∗ iand upu,∗ ji for all j∈Niis applied at time instant k, yielding the nominal control law κi(xi)=[upr,∗ i(0),[upu,∗ ji (0)]∀j∈Ni]⊤ .(16) To deal with the mutual disturbances ignored by the nominal model, an auxiliary control law is used so that the trajectory of 67 A. Sánchez-Amores, P. Chanfreut, J.M. Maestre et al. Journal of Process Control 123 (2023) 64–75 the real system follows closely the nominal one, trying to cancel the error between the nominal and real states, i.e., xiand xi. Consequently, the private and public inputs of each subsystem iare calculated as: [upr i [upu ji ]∀j∈Ni]=κi(xi)+Ki(xi−xi).(17) Remark 4. For all i∈N, the implemented input of subsystem i, expressed according to (1), is the sum of the private and public parts of the manipulated variable uidefined in (3), whereas (17) is the vector of input variables computed by agent i. In this regard, note that the private part of uiis computed by agent i, but the public part is decided by agents in Mi. Solving (15) provides for each subsystem i∈Nthe value of the private and public inputs used to obtain the real implemented input (17), i.e., upr,∗ i(0) and upu,∗ ji (0), and the corresponding optimal scale factors α∗ i, i.e., αi=α∗ iand αji =α∗ ji , for all neighboring agents j∈Ni. In what follows, we will describe how the centralized problem is solved in a distributed manner among the set of agents by using dual decomposition. Consequently, (15) is not intended to be solved directly. 3.2. Distributed MPC based on dual decomposition Let us consider the dual decomposition algorithm described in [23], which allows us to compute the solution of (15) in a distributed fashion. In this context, convergence is attained throughout an iterative negotiation procedure where Lagrange multipliers λiare used to coordinate coupling variables. Consider the problem in (15) and note that if the dynamics of different agents is affected by a shared input ui, they will need to negotiate the value of the scale factors that define the input constraints, i.e., αi,[αij]j∈Ni. Moreover, the solution needs to satisfy condition (4). The latter will be enforced by the introduction of Lagrange multipliers, which become new parameters of the local objective functions, and thus influence the agents’ solutions. These multipliers are assumed to remain constant during the prediction horizon, however, note that they may vary at each time instant kand throughout the dual decomposition iterations. In particular, to comply with constraints (3) and (4), agents that carry out the negotiation are required to satisfy: λi⎛ ⎝αi+∑ j∈Mi αij −1⎞ ⎠≤0 with λi≥0.(18) Let Sibe the set containing subsystem iand its neighboring agents, i.e., Si≜{{i}∪Ni}. Accordingly, there will be |Si|Lagrange multipliers in the local objective function of subsystem i. In other words, there will be as many Lagrange multipliers as private and public input local variables. These auxiliary variables are introduced in the local objective function of each subsystem i∈ Nas Λi(αi,[λm]m∈Si)=λiαi+∑ j∈Ni λjαji.(19) Remark 5. Note that (18) needs to be fulfilled for every agent i∈N. Consequently, (19) is the result of keeping all the terms of the |N|expressions of (18) that multiply the variables controlled by agent i. As a result, to distribute the global problem (15), we can rewrite (9) taking into account the agents’ negotiation. The local objective function is formulated as the sum of the stage cost ℓi(·)(12), the penalization of scale factors fi(·)(13), the term introducing the Lagrange multipliers for the dual decomposition algorithm Λi(·)(19), and the terminal cost Vf i(·)(14): Ji(xi,upr i,[upu ji ]∀j∈Ni,αi,[λp m]m∈Si)= k+Np−1 ∑ t=k ℓi(xi(t),upr i(t),[upu ji (t)]j∈Ni) +fi(αi)+Λi(αi,[λp m]m∈Si)+Vf i(xi(k+Np)). (20) Consequently, at each time instant kand iteration step p, each agent i∈Nsolves: min [upr i,[upu ji ]∀j∈Ni ,αi]i∈N Ji(xi,upr i,[upu ji ]∀j∈Ni,αi,[λp m]m∈Si) s.t. xi(k)∈xi⊕(−Ri), xi(t+1) =Aiixi(t)+Biiupr i(t)+∑ j∈Ni Bijupu ji (t), upr i(t)∈αp iUi,upu ji (t)∈αp jiUj,∀j∈Ni, xi(t)∈Xi, xi(Np)∈Xf i= {0}, λp m≥0,∀m∈Si, ∀t=k,...,k+Np−1. (21) where [λp m]m∈Siare introduced as parameters. The iterative negotiation takes the values of the variables involved in (19) obtained at any iteration pand compares them with those resulting from the previous iteration p−1. Practical convergence is attained when ∆≜|αp−αp−1|is below a threshold ϵ > 0. In this way, while ∆> ϵ agents will negotiate the bounds on shared variables. On the other hand, when ∆≤ϵ, it is considered that consecutive values are similar enough to finish the negotiation for time step k. Meanwhile, Lagrange multipliers are updated in each step of the negotiation according to: λp+1 i=λp i+γ⎛ ⎝αp i+∑ j∈Mi αp ij −1⎞ ⎠,(22) where γ > 0 is the step size. Finally, the scale factors αiand αij involved in (18) come from the solution to (21) at time instant k and iteration step p. 4. Control scheme The proposed coalitional control approach aims to reduce communication among agents and is summarized in Algorithm 1. For convenience, in what follows we use superscript −to refer to a variable at the previous time instant, e.g., x−denotes xat time instant k−1. We consider event-based coordination, meaning that communication links will be enabled or disabled according to the condition introduced in Step 9 of Algorithm 1. The communication link between two agents iand jwill be disabled if the public variable upu ji is small enough during two consecutive time steps. It is considered that there is no need for coordination depending on the value of αji because it scales the constraint space of the public variable. When αji takes small values, the region of the public variable is negligible compared to the values that the private variable can take. Nevertheless, communication will be restored if at a single time step kagent irequires greater use of the public variable upu ji , being the value of its scale factor αji above the threshold ˆα. 68 A. Sánchez-Amores, P. Chanfreut, J.M. Maestre et al. Journal of Process Control 123 (2023) 64–75 Agents that are communicating follow an iterative procedure under the distributed dual decomposition algorithm. They share the values of the scale factors αiand αji until the negotiation converges, updating the value of the Lagrange multipliers at every iteration step p. It is worth mentioning that, unlike common distributed approaches, agents share the variables that bound input constraints and not input variables directly. In addition, agents will keep calculating the value of their public variables even when they are not communicating since the negotiation can be suddenly resumed. Algorithm 1 Control Scheme Initialization: At the first time instant k=0 and iteration p=0, the values of the scaling factors are set to a positive non-zero value: αp=0 i(0) =αaux,αp=0 ji (0) =αaux ∀i,j∈N, with αaux >0. If k>0, we will set the value of the scale factor at p=0 as the one in the previous time instant, i.e., αp=0 i(k)=αi(k−1), αp=0 ji (k)=αji(k−1) ∀i,j∈N. Moreover, set [λm]m∈Si=0 at p=0∀k. At each sample time k, each agent i∈Nproceeds as follows: 1: while ∆> ϵ do 2: Update the disturbance set Wp i(7) according to the scale values from the previous iteration, i.e., [αp−1 ij ]∀j∈Mi and [αp−1 j]∀j∈Ni . 3: Compute the corresponding invariant set Rp i. 4: if xi−xi/∈Rp ithen 5: Set Rp i=R− i. 6: end if 7: Solve (21) to obtain the optimal nominal sequences upr,∗ i, upu,∗ ji and the optimal scale factors α∗ i. 8: Calculate the real inputs upr iand upu ji ∀j∈Niaccording to equation (17). 9: if αp ji <ˆαand α− ji <ˆαthen 10: Disable the communication between agents jand i. 11: Ignore the value of the public variable: upu ji =0. 12: The associated Lagrange multiplier to αji will remain constant in the next iteration: λp+1 j=λp j. 13: Set ∆j=0. 14: else 15: Enable the communication between agents jand i. 16: Update the neighboring Lagrange multiplier λp:+1 jaccording to (22). 17: Compute ∆j=max{|αp j−αp−1 j|,|αp ji −αp−1 ji |}. 18: end if 19: Set p←p+1 and compute ∆=max{∆j}∀j∈Ni. 20: end while The set Rp ihas to be recalculated for each iteration pand for every time instant k. Note that, in the literature, there are methods such as [24,25] that allow us to perform this calculation. In particular, the invariant sets for Algorithm 1are designed according to [26], where a one-step procedure is applied to compute a polytopic minimal robust positively invariant (mRPI) set. In particular, this method solves a single LP to compute the so-called (P,r)-mRPI, where rdenotes the number of inequalities defining the set and Pis a predefined matrix. Therefore, the update of Rp i entails a light computational load. 4.1. Verify whether xi−xi∈Rp i Since input, state, and uncertainty constraint sets vary over time, the RPI set must be recomputed at every instant kto suit the new constraints. Following [7], a checking step is introduced to verify if the difference between the predicted and nominal states in k+1 belongs to the new RPI calculated at that time instant in order to maintain constraint satisfaction and recursive feasibility. According to [22], if the state belongs to the set xi⊕Rifor a certain time instant k, we can guarantee that for the successive time instant k+1 the state will remain within the set x+ i⊕Ri. At the beginning of the algorithm, k=0, R0 iis calculated, and we guarantee that x0 i−x0 i∈R0 i, since the initial nominal state is optimized according to xi(0) ∈xi⊕(−Ri). The problem arises when we have to recalculate for each time step kthe RPI set because the constraints that affect Rimay have changed with respect to the previous instant. According to xi∈xi⊕Ri, at time instant k, successor states x+ iand x+ iare calculated using the invariant set Ri. As condition x+ i∈x+ i⊕Riholds, we have to guarantee that the new RPI calculated at k+1 contains the difference between the predicted real and nominal current states: x+ i−x+ i∈R+ i.(23) To pursue this goal, Step 3 of Algorithm 1is introduced with the same procedure as the checking step introduced in [7] (Subsection 4.1). If (23) is met, the following is satisfied: AKi(x+ i−x+ i)⊕W+ i⊆R+ i.(24) Consequently, the predicted state and input trajectory will verify the original constraints Xiand Ui. Conversely, if x+ i−x+ i/∈R+ iwe will discard the new RPI set and use the previous one, i.e., Ri:= R− i, since x+ i∈x+ i⊕Riis always fulfilled. 5. Recursive feasibility and stability In what follows, we proceed as in [7] to develop the conditions that ensure recursive feasibility and stability since the applied schemes are somewhat similar except for the novel distribution of inputs among agents, and the exchanges of communication that are necessary to adjust the bounds on shared variables. 5.1. Decreasing trend of the disturbance set W Hereafter, we will introduce several assumptions that allow us to preserve the guarantees of recursive feasibility and stability when using the previous RPI if the new set does not satisfy condition (23). Assumption 2. The set of disturbances are non-increasing from one time instant kto another k+1, i.e., W+ i⊆Wi. Remark 6. Since we are working with a regulation control problem, the state will approach the origin as time goes by, and hence the inputs will also tend to zero. In this regard, any nonincreasing evolution of factors αj(j∈Ni) and αij (j∈Mi) will lead to coupling disturbances sets satisfying Assumption 2 (recall (6) and (7)). In this way, the bounds of the uncertainties evolve into smaller sets over time. This implies that if the disturbance set changes in two consecutive instants, the successor set will be a subset of the previous one. As the size of the RPI set depends on the disturbance set, a smaller disturbance set entails a smaller RPI set, leading to the following assumption: Assumption 3. The invariant sets are non-increasing from one time instant kto another k+1, i.e., R+ i⊆Ri. 69 A. Sánchez-Amores, P. Chanfreut, J.M. Maestre et al. Journal of Process Control 123 (2023) 64–75 Let us consider that at time instant k, subsystem isatisfies: AKiRi⊕Wi⊆Ri. Considering the above, for the successor instant, it holds that W+ i⊆Wi, and R+ i⊆Ri. Hence, for k+1 is satisfied: AKiR+ i⊕W+ i⊆R+ i. Let us define Wias any uncertainty set contained in Wi, i.e., Wi⊆Wi. As disturbances decrease over time, AKiRi⊕Wi⊆Riis fulfilled for any uncertainty set contained in Wi. A fail-safe option for k+1 if x+ i−x+ i/∈R+ i, is to replace the new RPI with the preceding one, i.e., set R+ i=Ri. This is why Step 4 of the control algorithm is proposed as an alternative when the necessary conditions are not met, as we are being more conservative with a larger RPI set that satisfies the problem’s constraints. 5.2. Conditions for recursive feasibility Theorem 1 (Recursive Feasibility).Consider that Assumptions 2and 3hold. Then, if at the initial time instant k =0it is possible to find a feasible solution of problem (21) for all agents i ∈N, it is guaranteed that it exists a feasible solution for any instant k ∈ {0,1,2, . . .}.□ Proof. Let x∗ i=[x∗ i(0)⊤,x∗ i(1)⊤,..., x∗ i(Np)⊤]⊤and upr,∗ i(k)=⎡ ⎢ ⎢ ⎢ ⎣ upr,∗ i(0) upr,∗ i(1) . . . upr,∗ i(Np−1) ⎤ ⎥ ⎥ ⎥ ⎦ ,upu,∗ ji =⎡ ⎢ ⎢ ⎢ ⎣ upu,∗ ji (0) upu,∗ ji (1) . . . upu,∗ ji (Np−1) ⎤ ⎥ ⎥ ⎥ ⎦ ∀j∈Ni be the optimal state and input sequences derived from the solution of (21) at instant k. Additionally, note that Assumptions 2 and 3imply that W+ i⊆Wiand R+ i⊆Ri, with W+ iand R+ ibeing respectively the disturbances and RPI sets at instant k+1, and Wi and Rithose at k. Finally, let us build up the following candidate solution for instant k+1: ˜ x+ i= ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ˜ x+ i(0) ˜ x+ i(1) . . . ˜ x+ i(Np−1) ˜ x+ i(Np) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ x∗ i(1) x∗ i(2) . . . x∗ i(Np) 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ,(25) ˜ upr,+ i= ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ upr,∗ i(1) upr,∗ i(2) . . . upr,∗ i(Np−1) 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ ,˜ upu,+ ji = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ upu,∗ ji (1) upu,∗ ji (2) . . . upu,∗ ji (Np−1) 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ∀j∈Ni.(26) From the constraints of (21), we have that xi−x∗ i(0) ∈Riand x∗ i(Np)= {0}. Additionally, if we apply control law (17), it holds that x+ i−x∗ i(1) ∈Ri. However, as the RPI is recalculated at each time step, condition x+ i−x∗+ i(1) ∈R+ imust be satisfied at k+1. Since ˜ x+ i(0) =x∗ i(1), then x+ i−˜ x+ i(0) ∈Ri. Moreover, according to the checking step described in Section 4.1, if x+ i−˜ x+ i(0) /∈R+ i, then we set R+ i:= Rias a fail-safe option, i.e., x∗ i(1) will always satisfy the constraint on the initial nominal state of problem (21). Additionally, given that R+ i⊆Ri, nominal constraints sets Ui and Xican only enlarge in time, i.e., U+ i⊇Ui,X+ i⊇Xi. Therefore, states x∗ i(n) and inputs upr,∗ i(n),upu,∗ ji (n)∀j∈Niin (25) and (26) are admissible for all n=1,...,Np. Finally, since x∗ i(Np)= {0}, we can always apply zero as the terminal input to stay at the origin. Therefore, (25) and (26) constitutes a feasible solution of problem (21) for time step k+1. By induction, the theorem is proven. ■ 5.3. Conditions for stability Theorem 2 (Stability).For all agents i ∈N, the origin is asymptotically stable. □ Proof. To establish stability, we have to prove that the nominal cost function Ji(20) decreases during the evolution of the system. Let the optimal sequences in 1x∗ i,upr,∗ iand upu,∗ ji at klead to the optimal cost J∗ i. Additionally, the candidate state and input sequence for k+1, i.e., ˜ x+ i,˜ upr,+ iand ˜ upu,+ ji , lead to the cost ˜ J+ i, which provides a value for the upper bound of the optimal cost: ˜ J+ i≥J∗+ i. Moreover, we can calculate the difference ˜ J+ i−J∗ ias: ℓi(˜ x+ i(Np−1),˜ upr,+ i(Np−1),[˜ upu,+ ji (Np−1)]∀j∈Ni) −ℓi(x∗ i(0),upr,∗ i(0),[upu,∗ ji (0)]∀j∈Ni) +˜ f+ i(αi)−f∗ i(αi) +˜ Λ+ i(αi,[λm]m∈S)−Λ∗ i(αi,[λm]m∈S) +Vf i(˜ x+ i(Np))−Vf i(x∗ i(Np)). (27) Unlike private and public inputs upr iand upu ji , optimal scale factors in α∗ iremain constant along the prediction horizon. Moreover, the values of the Lagrange multipliers for the candidate solution do not vary along Np, so ˜ λ+ m=λ∗ m. We have assumed that the values of optimization variables that remain constant over the prediction horizon are also kept constant from the optimal solution in kto the candidate solution in k+1. Therefore, ˜ f+ i=f∗ i and ˜ Λ+ i=Λ∗ i, so they cancel themselves in (27). As a result, in the subtraction ˜ J+ i−J∗ ithere are only left terms referring to the stage and terminal costs. Furthermore, the terminal cost is a continuous Lyapunov function at the origin, meaning that Vf i(k+1) −Vf i(k)≤0, which entails: ℓi(˜ x+ i(Np−1),˜ upr,+ i(Np−1),[˜ upu,+ ji (Np−1)]∀j∈Ni) +Vf i(˜ x+ i(Np))−Vf i(x∗ i(Np)). (28) Therefore, we can rewrite (27) as: ˜ J+ i−J∗ i≤ −ℓi(x∗ i(0),upr,∗ i(0),[upu,∗ ji (0)]∀j∈Ni),(29) which also implies J∗+ i−J∗ i≤ −ℓi(x∗ i(0),upr,∗ i(0),[upu,∗ ji (0)]∀j∈Ni).(30) We can extend the (30) for all subsystems i∈N: J=∑ ∀i∈N Ji→J∗+ −J∗≤ −ℓ(x(0),u(0)).(31) As the stage cost ℓ(·) is strictly positive, it is proven that the cost function Jidefined in (20) decreases over time, ensuring the stability of the system. ■ 6. Simulation results In this section, we apply the proposed coalitional control algorithm to the academic example shown in Fig. 2. The eight input-coupled tanks plant consists of four top tanks (5, 6, 7, 8) that discharge flow into four bottom tanks (1, 2, 3, 4), and these, in turn, discharge into a shared storage tank. Four pumps (Qa,Qb,Qc,Qd) are used to fill the tanks, carrying water from the storage tank to the tanks indicated in Fig. 2. As can be seen, flow regulation is done through three-way valves, which divide the pumped flow into two ways to fill the tanks. The global system can be divided into N=4 subsystems that consist of a top and bottom tank. Thus, the first subsystem 70 A. Sánchez-Amores, P. Chanfreut, J.M. Maestre et al. Journal of Process Control 123 (2023) 64–75 is formed by tanks #1 and #5; tanks #2 and #6 describe the second one; the third subsystem is composed of tanks #3 and #7; and the fourth one is formed by tanks #4 and #8. As shown in Fig. 4, subsystems are physically coupled through the colored pipes that connect their tanks. On the other hand, we refer to the data connections between their corresponding local controllers as communication links. In particular, link (i,j) represents the bidirectional data connection between agents iand j. Whenever a given link is enabled, the agents connected through it will form a coalition and will be able to share information. In this regard, we consider the following links: (4,1), (2,3), (1,2), and (3,4). The target is to regulate the lower tanks towards their operating point in terms of water level, i.e., h◦ i. To this end, the state of each subsystem is defined as the water level measured from the operating point. Taking the fourth subsystem as an example, its state is given by: x4=[h4−h◦ 4,h8−h◦ 8]T. Additionally, the inputs are given by the difference of the pump flow and its value at the operating point ui=Qk−Q◦ k, with i=1,...,Nand k∈ {a,b,c,d}. Let us describe the operating point of the plant as: h◦ 1=0.6487,h◦ 2=0.6639,h◦ 3=0.6534,h◦ 4=0.6521, h◦ 5=0.6498,h◦ 6=0.6592,h◦ 7=0.6594,h◦ 8=0.6587, Q◦ a=1.63,Q◦ b=2,Q◦ c=1.8,Q◦ d=2, (32) where water level is measured in meters and flows in cubic meters per hour. Moreover, we can characterize each subsystem with the following matrices: A11 =[0.8257 0.1178 0 0.8703],B11 =[0.0379 0],B14 =[0.0056 0.0843], A22 =[0.8163 0.1023 0 0.8867],B22 =[0.0503 0],B21 =[0.0053 0.0916], A33 =[0.8232 0.1077 0 0.8813],B33 =[0.0442 0],B32 =[0.0047 0.0783], A44 =[0.8194 0.1050 0 0.8840],B44 =[0.0441 0],B43 =[0.0050 0.0849]. Accordingly, the disturbance vectors (recall (7)) for the simulation example are such that: w1∈W1≜B11Wpu 12 ⊕B14Wpr 4, with Wpu 12 =α12U1,and Wpr 4=α4U4, w2∈W2≜B22Wpu 23 ⊕B21Wpr 1, with Wpu 23 =α23U2,and Wpr 1=α1U1, w3∈W3≜B33Wpu 34 ⊕B32Wpr 2, with Wpu 34 =α34U3,and Wpr 2=α2U2, w4∈W4≜B44Wpu 14 ⊕B43Wpr 3, with Wpu 14 =α14U4,and Wpr 3=α3U3. Lastly, the system is subject to the following constraints: 0.2≤ hn≤1.3,∀n=1,...,8, 0 ≤Qa≤3.26, 0 ≤Qb≤4, 0 ≤Qc≤ 3.6, and 0 ≤Qd≤4. Note that these conditions must be adapted to the state and input constraints by subtracting the operating point (32) from the previous limit values. The simulation has been done using as weighting matrices: Qi=I2,Rpr i=0.20, Rpu i=2Rpr i, and ρpr =Rpr,ρpu =Rpu penalize scale factors in the objective function (13), calculated by trial and error. The prediction horizon has been set to Np=10, and the simulation lasts 25 s. Finally, the threshold ˆαis a tunning parameter, and its value has been set as ˆα=0.015 by trial and error. In what follows, the proposed coalitional method is compared with a decentralized MPC, where agents optimize their local control objectives regardless of the control goals of the rest of Fig. 2. Diagram of the eight tanks plant, where the coupling relations are represented by the colored pipes. the system, and with a centralized MPC, representing a fully cooperative scheme. Note that the centralized MPC solution can also be found in a distributed fashion, e.g., using dual decomposition DMPC after convergence is attained. To that end, the agents need to share information regarding their input sequences at every time instant, whereas communication in the proposed coalitional method is less demanding and event-based (recall Algorithm 1). Fig. 3(a) represents the evolution of the water level of lower tanks towards their equilibrium point. As they all start with a water level that is above their operating point, they must be slightly emptied and, therefore, the flow rate of the four pumps must decrease. This is shown in Fig. 3(b), which illustrates the performance of Algorithm 1in terms of the evolution of flow rates. As can be seen, when using the coalitional scheme with ˆα=0.015, water level and pump flow trajectories stay close to the centralized solution. This fact highlights the benefits provided by the proposed method, as a very similar behavior to the centralized one is achieved without permanent communication between the agents. On the other hand, the water level evolution towards their operating point h◦ iis slower with the decentralized MPC due to the lack of coordination, which makes the agents act independently of each other’s control objectives. As a consequence, the pump operation differs substantially from that implemented with the centralized and coalitional approaches. Let us define Pas the index that evaluates the system’s performance: P= T ∑ k=1∑ i∈N u⊤ i(k)Rpr iui(k)+x⊤ i(k+1)Qixi(k+1),(33) being Tthe simulation length. Note that the definition of the performance index allows one to compare the three assessed methods. In the coalitional approach, the input variable ui(k) in (33) is defined as the sum of the private and public parts of the variable according to (3). In this regard, Table 2 provides a comparison of the performance index when using centralized, coalitional, and decentralized MPC strategies. Applying the proposed coalitional approach results in a 1.28% decrease in the overall performance regarding the centralized solution. This was expected since coalitional control leads to a slight loss of performance in exchange of savings in cooperation and communication efforts. In turn, the decentralized MPC strategy implies a performance loss of a 17.06% compared to the centralized solution, as there is no communication between the different agents. Thus, it can be seen that, in terms of performance, it is beneficial to have some coordination between agents, as there is a significant 71 A. Sánchez-Amores, P. Chanfreut, J.M. Maestre et al. Journal of Process Control 123 (2023) 64–75 Fig. 3. Evolution of the water level of the lower tanks and flow rate of the four pumps. Solids lines represent the result using the proposed coalitional algorithm using ˆα=0.015, while dashed and dotted lines represent the evolution using a centralized and decentralized MPC respectively. The fine dashed–dotted line represents the operating point of the plant (32). Table 2 Comparison of the overall performance for different control approaches. Performance index P Centralized MPC 4.8189 Coalitional algorithm with ˆα=0.015 4.8805 Decentralized MPC 5.6411 loss in the performance when agents work independently in the decentralized scheme. Our coalitional approach aims to reduce communication between agents while minimizing performance losses with respect to the centralized approach. Therefore, performance is slightly decreased in order to save on cooperation efforts. To check this, we have calculated a value for the communication cost as the total number of links activated along the simulation. While in Fig. 4. Communication topology by means of the state of the communication links, i.e., enabled/disabled, using ˆα=0.015. 72