Coalitional Model Predictive Control based on Partial Cooperation

Tesis Doctoral Coalitional Model Predictive Control based on Partial Cooperation Autor: Ana Sánchez Amores Directores: José María Maestre Torreblanca Eduardo Fernández Camacho Dpto. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2024 Tesis Doctoral Coalitional Model Predictive Control based on Partial Cooperation Autor: Ana Sánchez Amores Directores: José María Maestre Torreblanca Catedrático de Universidad Eduardo Fernández Camacho Catedrático de Universidad Dpto. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla 2024 Agradecimientos A gradecer, en primer lugar, a mis directores de tesis la ayuda, implicación y todas las enseñanzas que me han inculcado durante estos años de doctorado. En especial, agradecer a Pepe la confianza depositada en mí desde aquel día en 2018 que entré por su despacho para solicitarle un trabajo de fin de grado. Gracias por ayudarme a desarrollar una parte de mí que no conocía, por alentarme cuando más lo necesitaba, y por hacer que el mundo de la investigación y los conceptos abstractos me apasionen tanto. Agradecer a Eduardo la oportunidad que me brindó de trabajar dentro del proyecto OCONTSOLAR durante estos años, así como su disposición para ayudar y orientarme en el camino cuando este se volvía un poco más complicado. También quiero agradecer a Antonio Ferramosca por haberme acogido durante mi estancia en la Universidad de Bérgamo como uno más; ha sido un placer trabajar contigo. Y, por supuesto, quiero agradecer a Paula todas sus enseñanzas durante los primeros años; sin ti, el inicio de esta etapa habría sido mucho más difícil. Agradecer a mi hermana y a mis padres su cariño y apoyo incondicional durante estos años. Gracias por confiar en mí, por el orgullo que me expresáis que sentís, y por ayudarme a ver que puedo lograr todo lo que me proponga. Por último, agradecer a mis amigas, y en especial a Asun, por no soltarme la mano en ningún momento durante este camino, y por motivarme siempre a seguir avanzando. Ana Sánchez Amores Sevilla, 2024 I Resumen E sta tesis doctoral investiga el desarrollo controladores predictivos coalicionales basados en modelo. Nuestro estudio hace hincapié en la necesidad de gestionar sistemas complejos y de gran escala manteniendo la eficiencia y la escalabilidad de las estrategias de control. Con este fin, se exploran los límites de la descentralización dentro del control coalicional, desarrollando algoritmos que minimizan la cantidad de datos que se comunican en la red para preservar la privacidad y mejorar la autonomía, logrando un equilibrio entre los esfuerzos computacionales y el rendimiento. Aunque este último pueda degradarse en comparación con otros métodos coalicionales, se disfruta a cambio de una mayor autonomía y un menor esfuerzo de comunicación. En primer lugar, se ha realizado una revisión exhaustiva de la literatura existente, resaltando la necesidad de implementar estrategias de control descentralizadas para gestionar sistemas de gran escala. También se introducen los principales objetivos de esta tesis doctoral, cuyas contribuciones se pueden dividir principalmente en dos partes. En la primera parte de la tesis, se presenta una nueva descomposición de variables diseñada para sistemas parcialmente cooperativos. Esta metodología implica la distribución de las variables de acoplamiento entre los agentes locales. Se desarrollan también estrategias de control predictivo coalicional estocásticas y robustas, permitiendo a los agentes compartir la autoridad de control sobre ciertas variables de entrada. Además, se han investigado los límites de esta delegación de autoridad, estableciendo la porción máxima que un agente puede ceder de su entrada a otro sin comprometer la factibilidad del problema de optimización global. Este enfoque asegura que las operaciones descentralizadas respetan las restricciones del sistema con una comunicación mínima de datos. La segunda parte de la tesis se enfoca en sistemas con recursos compartidos y se desarrollan métodos coalicionales para desacoplar restricciones conjuntas, permitiendo que los agentes operen con mayor autonomía. Dado que esta tesis se ha desarrollado en el marco del proyecto OCONTSOLAR, los algoritmos de control presentados en esta segunda parte se aplican al caso de estudio de plantas solares a gran escala. Finalmente, se ofrecen conclusiones y se sugieren direcciones para futuras investigaciones basadas en las principales contribuciones de esta tesis doctoral. III Abstract T his doctoral thesis investigates the development of coalitional strategies based on model predictive control (MPC). Our study emphasizes the necessity of managing increasingly complex systems while maintaining efficiency and scalability. To this end, we explore the limits of decentralization in coalitional control, developing algorithms that minimize the amount of data communicated in the network to preserve privacy and enhance autonomy, while still maintaining acceptable levels of performance. The proposed strategies achieve a trade-off between computational efforts and performance. While performance may degrade compared to other coalitional methods, this trade-off allows for increased autonomy and reduced communication requirements. First, we survey the existing literature, motivating the need for non-centralized control strategies to manage large-scale systems, and we introduce the main objectives of this PhD dissertation. The contributions of this thesis can be broadly classified into two main parts. In the first one, we introduce a novel variable decomposition for partially cooperative systems. In this regard, we distribute coupling variables among local agents, developing stochastic and robust coalitional MPC approaches that allow agents to share control authority over certain input portions while maintaining privacy over others. Moreover, we study the limits of control authority delegation by establishing the maximum portion of an agent’s input space that can be ceded to another agent while ensuring the feasibility of the global optimization problem. This ensures that decentralized operations respect system constraints with minimal data communication. On the other hand, the second part focuses on systems with shared resources, developing coalitional methods to decouple joint constraints, allowing agents to operate with greater autonomy. Since this thesis has been developed in the framework of OCONTSOLAR project, the control algorithms presented in this second part are applied to the case study of large-scale solar plants. Finally, we provide conclusions and future research directions of our main contributions. V XII Notation αji Scale factor of the public input such that upu ji ∈αjiUj. NiNeighboring agents of i. MiAgents affected by ui:{j∈N:Bji =0,j=i}. wiDisturbances induced by coupling of subsystem i. WiDisturbance constraint set of subsystem isuch that wi∈Wi. Tube-based setting RiRobust positively invariant set of subsystem i. ¯xiNominal state of subsystem i. ¯uiNominal input variable of subsystem i. ¯ Xi Shrinked constraint set for the nominal state of subsystem i such that ¯xi∈¯ Xi. ¯ Xi Shrinked constraint set for the nominal input of subsystem i such that ¯ui∈¯ Ui. Tracking setting ¯xs iArtificial state for tracking. ¯upr,s iArtificial private input for tracking. ¯upu,s ji Artificial public input for tracking. Sets AΩImage of a set Ω⊂Rnunder the linear mapping A:Rn→Rp. A−1ΩPreimage of a set Ω⊂Rn. Ω⊕ΦMinkowski set addition for Ω,Φ⊂Rn:{x+y:x∈Ω,y∈Φ}. Ω⊖Φ Minkowski set difference for Ω⊂Rn , Φ⊂Ω : {x∈Rn:Φ+x⊂ Ω}. λΩScaling of a set Ω⊂Rnby a factor λ∈R:λΩ≜{λx|x∈Ω}. hΩ(η) Support function of a non-empty set Ω⊂Rn in the direction of the vector η: sup{η⊤x:x∈Ω}. C-set Compact and convex set containing the origin. PC-set Compact and convex set containing the origin in its interior. Glossary DMPC Distributed Model Predictive Control. DNI Direct Normal Irradiance. HTF Heat Transfer Fluid. LP Linear Programming. LSS Large Scale Systems. LTI Linear Time Invariant. MPC Model Predictive Control. PTC Parabolic Trough Colector. RCI Robust Control Invariant. RPI Robust Positively Invariant. XIII Part I Introduction 1 1 Introduction T his thesis has been performed within the OCONTSOLAR (Optimal Control of Thermal Solar Energy Systems) project [1], which aims to develop new algorithms to enhance the efficiency and operation of large scale solar fields. Being large scale systems (LSS), solar fields lead to a high computational burden due to their high number of optimization variables. Thus, solving the optimization problem in a centralized manner can be computationally demanding considering the available processing times. While centralized control provides a comprehensive approach to manage complex systems and takes advantage of the concentration of decision-making and computational tasks in a single entity, it becomes increasingly burdensome as systems grow in size and complexity [2], generating bottlenecks, elevated communication costs, and susceptibility to a single point of failure. The inability to find a centralized model for the global system and the significant computational load of the underlying optimization are other issues that explain the lack of scalability of centralized approaches. Additionally, in scenarios where the system is spatially distributed, e.g., micro-grids [3, 4] and traffic networks [5, 6], or exhibits decentralized dynamics such as swarm robotics [7,8], the demands for information exchange and coordination in centralized control may become impractical. In such cases, it is preferable to rely on non-centralized control methods, where the overall system is divided and several independent controllers work in parallel. These methods are also suitable for LSS, which can be modeled as a set of interconnected geographically-distributed subsystems with coupled dynamics [9]. From the different control frameworks available, this thesis focuses on model predictive control (MPC) [10 – 12], which is a continuous replanning method that recalculates control actions at each time step in a receding horizon fashion. In particular, it 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 and the problem constraints. The MPC optimization problem can also deal with other complex issues such as delays, information about expected disturbances, and the 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 3 4Chapter 1. Introduction step is implemented and the problem is solved again at the next time instant following a receding horizon approach. This well-known strategy is the most general method to work with control problems in the time domain and has been applied in multiple industrial and academic applications for decades [13–15]. As happens with other control methods, MPC has difficulties when dealing with LSS in a centralized way, e.g., due to the number of decision variables in the optimization problem. In response to these challenges, non-centralized control strategies partition the overall system into a set of dynamically coupled smaller subsystems, each of them assigned to a local controller or agent [16, 17], as depicted in Figure 1.1. During the last decade, we have witnessed a significant growth of distributed MPC (DMPC) approaches due to its capacity to retain the benefits of the MPC formulation, while overcoming scalability issues, enhancing robustness, and improving overall adaptability. The basic idea behind DMPC is to tackle the overall control problem by distributing the decision-making authority among multiple local MPC controllers that, in turn, may exchange information. Within the DMPC 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. Depending on the information-sharing among the different control entities, we can broadly differentiate distributed approaches, where agents share data through a communication network [18 – 21], from decentralized strategies, where no data exchange occurs between local controllers [22 – 24]. Simultaneously, within distributed strategies we can broadly distinguish non-cooperative algorithms from cooperative schemes [25,26]. In non-cooperative approaches, such as communication-based MPC [27], each agent optimizes a local objective function regardless the control objectives of other components of the system, potentially leading to a Nash equilibrium. Conversely, cooper- Figure 1.1 Scheme of a distributed control system composed of 7 subsystems. Each subsystem is assigned to a local agent i that calculates its control input ui based on the local state xi and coupling data. Dashed lines represent the communication network among the local agents, and solid lines represent the dynamical coupling between the different subsystems. 1.1 Coalitional control 5 ative algorithms, such as feasible cooperation-based MPC [28], or dual decomposition methods [29], can preserve the optimal centralized performance, allowing for a Pareto efficient solution. However, this is achieved at the expense of a greater communication and computation burden. Conversely, while the absence of coordination in decentralized schemes minimizes communication demands, it compromises overall performance [30,31]. A halfway approach between decentralized and fully-cooperative control schemes is provided by distributed control approaches that promote the formation of dynamic clusters of local controllers [32 – 34]. Indeed, much of the recent research within the OCONTSOLAR project has focused on developing clustering-based algorithms to reduce computation times and increase controller scalability for the type of solar plants considered [35 – 39]. 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 [40 – 42]. More precisely, local controllers merge into cooperative clusters or coalitions only when it improves global performance, thus reducing the overall cooperation burden (see an example of these coalitions in Figure 1.2). In particular, agents within a coalition operate independently from other local controllers, fostering a scenario of partial cooperation. Accordingly, there may be coalitions of agents exchanging information, while others may be operating in a decentralized fashion. This flexible approach aims to efficiently manage computational and communication resources by dynamically adapting the controller structure [43,44] . This thesis contributes to this partially cooperative framework by developing different coalitional MPC schemes [45 – 51] by providing local agents with limits regarding the extent of their partial cooperation [52,53]. In particular, we explore different coalitional approaches that minimize drastically the amount of data to be shared between cooperative controllers with respect to other coalitional approaches. This way, we propose various coalitional control schemes that are more alike to decentralized strategies than to distributed ones considering the amount of shared data within the communication network. 1.1 Coalitional control Considerable effort has been dedicated within the DMPC framework towards restricting information exchange and minimizing cooperation, given the substantial volume of data transferred between controllers. For example, flexibility regarding the extent of coupling Figure 1.2 Different coalition topologies for a network with 7 agents. 6Chapter 1. Introduction is exploited in [54], where local controllers broadcast information regarding how much of their constraint space they are using to improve performance in a distributed tube-based MPC. In this regard, coalitional control promotes data exchange only when it leads to a significant increment of performance [40, 55,56]. As a consequence, the control network topology changes dynamically, leading to dynamic partitioning of agents into disjoint clusters known as coalitions [43,57 – 59]. This idea of dynamic clusters of controllers has drawn attention within the context of multi-agent networked systems [32,60]. The introduction of network topology as an additional variable in optimization problems comes with added complexity. For instance, top-down coalitional strategies often involve solving mixed-integer optimization problems, treating network link states as discrete variables [40, 55]. Alternatively, some approaches consider randomly activating links, retaining only those that enhance performance [41, 56]. However, both methods have drawbacks, since the first type increases computational burden and typically requires a lower rate of network topology changes by a supervisory control layer to maintain feasibility. On the other hand, bottom-up methods are more difficult to analyze and introduce a level of suboptimality due to the lack of centralized information in link switching operations. 1.1.1 Dynamics of multi-agent systems Coalitional control also starts from a system partitioned into loosely coupled local controllers that pursue a common goal and coalesce whenever this enhances overall performance, resulting in time-varying clusters. Throughout this thesis, we consider a networked system defined by the following discrete, linear time-invariant (LTI) dynamics: x+=Ax +Bu +de,(1.1) where x∈Rnx and u∈Rnu are respectively the state and input vector of the system, variable d∈Rnx represents the external disturbances, and x+ is the state at the next time instant. From a decentralized viewpoint, the system can be divided into a set N={1,2,...,N} of coupled subsystems with LTI dynamics: x+ i=Aiixi+Biiui+di, with di=∑ j∈NiAi jxj+Bi juj+de i,(1.2) where xi∈Rnxi and ui∈Rnui are respectively the state and input vector of subsystem i∈N , and variable di∈Rnxi denotes the coupling dynamics with neighboring agents, i.e., Ni={j∈N:[Ai j,Bi j]=0,j=i} , plus the measurable external uncertainties de i . Also, we assume that the subsystems’ states and inputs are subject to the following constraints: xi∈Xiand ui∈Ui , where Xi and Ui are convex sets containing the origin, for all agents i∈N. 1.1.2 Control architecture and strategy The main goal of coalitional control is to leverage the strengths of collaboration while addressing challenges such as communication overhead and scalability, achieving a tradeoff between overall performance and coordination burden. To this end, we assume 1.2 Objectives of the thesis 7 that a local controller or agent governs each subsystem i∈N . Also, local agents can communicate through a configurable data network described by an undirected graph G= (N,L) , being N the set of nodes that represent the agents, and L the set of links, with L⊆LN={(i,j)|(i,j)∈N,i=j} . In coalitional control, these links can be switched dynamically to suit the changing requirements of the control scheme, so that the set of agents can be partitioned into disjoint communication elements or coalitions. Inside each coalition, agents exchange data and coordinate their actions, but there is no inter-cluster communication. At time instant k , the partition of the system, say P(k) = {C1,C2,...,C|P(k)|} , is defined as the set of |P(k)| non-empty clusters of agents. The size of Ci can range from a singleton |Ci|=1 , to the grand coalition |Ci|=N . Note that the grand coalition, corresponds to the centralized case, that is Ci={N} . Depending on the configuration of the control scheme, agents within a given coalition Ci operate either as a single entity or as distinct partially cooperative entities. This thesis addresses both scenarios. 1.2 Objectives of the thesis This PhD dissertation deals with partially cooperative model-based control methods. Specifically, we explore decentralized and distributed optimization techniques for input and constraint-coupled systems, aiming to develop flexible approaches that enable the distribution of interconnected optimization problems while minimizing the data shared in the communication network, thus enhancing the autonomy of local agents. More precisely, we address the following objectives: • Objective 1: Development of a variable decomposition scheme to distribute coupling variables among local agents. This setup grants agents the ability to share a public portion of their input, ceding control authority to neighboring agents while ensuring the confidentiality of the remaining private portion. Each part of the coupling variable is associated with a different mode of interaction, and the bounds on the latter can be negotiated, as the need arises, in a distributed fashion. This approach enhances autonomy and privacy, minimizing the necessity for extensive cooperation while preserving overall system performance. Within this objective we can distinguish two subgoals: – Objective 1.1: Development of stochastic and robust coalitional MPC approaches that permit partial cooperation. In this subgoal, we pursue the development of coalitional schemes that include the latter variable decomposition method, enabling agents to negotiate the bounds on the shared portion rather than the complete sequence of common inputs. Consequently, the negotiation characterizes the portion each agent receives from a shared input, while local agents will independently determine how much of the shared input to utilize within these negotiated limits. To assess variables that are not locally controlled, introducing coupling uncertainty, we leverage stochastic and robust methods. 14 Chapter 2. Ceding control authority in input-coupled systems disregarded, giving rise to the following expression x+ i=Aiixi+Biiui+di,with di=∑ j∈Ni Bi juj.(2.1) Note that equation (2.1) is naturally satisfied in various scenarios, including multiple industrial processes [62], electrical networks [33], traffic flow models [63], and irrigation canals [64]. As mentioned before, we propose the decomposition of a local input variables ui into private and public parts as ui=upr i+∑ j∈Mi upu i j ,(2.2) where (i) upr i is the private part of the variable, which is controlled exclusively by the agent that owns it, that is, agent i, and it must verify upr i∈αiUiwith αi∈[0,1]. (ii) upu i j is the public part of ui controlled by agent j and must verify upu i j ∈αi jUi with αi j ∈[0,1] . The set of agents j that can manipulate the public part of ui are defined as affected subsystems Mi={j∈N:Bji =0,j=i} . In general, the sets Ni and Mi are different and depend on the dynamics of the system. While Ni contains the set of neighbors that affect i , Mi defines the set of agents affected by the input of agent i. In particular, each agent i will locally control the private part of its input variable, upr i , and the public part of its neighboring inputs upu ji for all j∈Ni, satisfying upr i∈αiUi,upu ji ∈αjiUj,∀j∈Ni.(2.3) Consequently, as the dynamics of subsystem i are influenced by neighboring inputs uj , neighboring agents j∈Ni will determine the value of their private input upr j . Additionally, the set Mi of agents affected by subsystem i ’s input will determine the public part of ui , denoted as upu i j for all j∈Ni . Therefore, variables that cannot be locally controlled, such as {upr j}j∈Ni and {upu i j }j∈Mi , are treated as bounded disturbances from the viewpoint of subsystem i . Figure 2.1 illustrates this decomposition and negotiation process with three agents. To satisfy the input constraints ui∈Ui for all i∈N , we consider that, at each time instant k, agents decide the values of upr iand upu i j , for all i∈Nand j∈Mi, so that αi+∑ j∈Mi αi j ≤1,∀i∈N(2.4) holds. 2.2 Communication network 15 Figure 2.1 The proposed variable decomposition scheme for three agents is illustrated. Focusing on agent 1 , we can identify N1=3 and M1=2 . Consequently, agent 1 can manipulate a public portion of u3 , while agent 2 will manipulate a public portion of u1. 2.2 Communication network In contrast to other coalitional schemes, agents forming a coalition do not operate as a single controller. Instead of forming a unique control entity sharing the entire sequence of coupling variables, communicating agents share and negotiate the bounds of their coupled variables. This enables different types of cooperation and minimizes the amount of data to be communicated, reducing the necessity for extensive cooperation while preserving overall system’s performance. Consequently, agents negotiation characterizes the portion each agent receives from a shared input, while local agents will independently determine how much of the shared input to utilize within these negotiated limits. The sharing of control authority in this partially cooperative framework fosters collaboration among agents while maintaining a higher degree of autonomy and privacy in their operations. As described in Chapter 1.1.2, the state of communication links in coalitional control is dynamically adjusted according to the control system’s requirements, resulting in the creation of disjoint communication components or coalitions. For the control approaches proposed in Part II of this thesis (Chapters 3 to 5), communication links change based on the evolution of scaling factors αi j , which determine the size of the constraint sets of public variables upu i j . Specifically, in Chapters 3 and 4, we disable communication link (j,i) if the corresponding value of αi j falls below a threshold ¯ α , defined to enable or disable agents’ communication. Conversely, in Chapter 5, a communication link (j,i) is deactivated as the derivative of αi j, which scales the constraint set of upu i j , approaches zero. 16 Chapter 2. Ceding control authority in input-coupled systems 2.3 Local dynamics of each subsystem Considering the model (2.1) and the variable decomposition proposed in (2.2) , the dynamics of each subsystem i∈Ncan be modeled as: x+ i=Aiixi+Biiupr i+∑ j∈Ni Bi jupu ji +wi, with wi≜∑ j∈Mi Biiupu i j +∑ j∈Ni Bi jupr j.(2.5) In particular, wi groups the uncertainties generated as a result of the interaction between subsystems, and can be decomposed into two terms. The first term contains the public part of ui , which is managed by agents in Mi (those affected by it), satisfying upu i j ∈αi jUi . On the other hand, the second term includes the private variables of neighboring agents Ni, which satisfy upr j∈αjUj. This way, disturbances wiare bounded by: wi∈Wi≜M j∈Mi Biiαi jUi⊕M j∈Ni Bi jαjUj,(2.6) where Wi is a compact and convex polyhedral set containing the origin in its interior, for all i∈N . In this regard, throughout Part II of this thesis, we will address disturbances wi in various ways. Specifically, Chapter 3 explores a stochastic approach through scenariobased MPC, while Chapters 4 and 5 will adopt a robust approach employing tube-based MPC, expanding the addressed scope with each chapter. 2.4 Distribution of the coupled constraint Since agents remain coupled through a common constraint regarding their constraint space, and the optimization problem is not meant to be directly solved in a centralized manner, we employ dual decomposition [29] to decouple local optimization problems. Agents affected by inputs coupled through shared constraint (2.4) need to negotiate the scale factors that limit the decision space of their public and private variables. Convergence to the centralized solution is attained via an iterative process where Lagrange multipliers λi ensure the satisfaction of the coupled constraint, yielding coordinated scale factors αi and αi j for all j∈Mi . To accomplish this, Lagrange multipliers are incorporated into the local objective functions of the agents involved in the negotiation such that λi αi+∑ j∈Mi αi j −1!≤0with λi≥0.(2.7) Thus, Lagrange multipliers are introduced into the local objective functions as Λiα α αi,[λm]m∈Si=λiαi+∑ j∈Ni λjαji,(2.8) 2.5 Underlying optimization 17 where Si≜{i}∪Ni defines the set comprising agent i and its neighbors. Consequently, each local objective function will include |Si| Lagrange multipliers, corresponding to the number of inequalities (2.4) that a particular agent must satisfy. As noted, dual decomposition involves a distributed negotiation of Lagrange multipliers to ensure that condition (2.4) holds for each subsystem. At each iteration p , the values of variables involved in (2.8) , specifically αi and αi jj∈Mi , are compared with those from the previous iteration p−1. Practical convergence is attained when ∆i=maxn|αp i−αp−1 i|,|αp i j −αp−1 i j |o∀j∈Mi,(2.9) falls below the threshold value ε>0 , concluding the negotiation for time step k . In this way, while ∆i>ε agents will negotiate the bounds on shared variables. On the other hand, when ∆i≤ε , it is considered that consecutive values are similar enough to finish the negotiation for time step k . Meanwhile, Lagrange multipliers are adjusted at each iteration step pusing γ>0as step size: λp+1 i=λp i+γ αp i+∑ j∈Mi αp i j −1!.(2.10) 2.5 Underlying optimization In this section, we present the basic objective function to be employed in configurations that apply this variable decomposition. It is important to emphasize that this expression serves as the foundation upon which more complex expressions will be constructed. Thus, the primary objective function minimized by each local controller i∈N over the prediction horizon Npis: Vixi,upr i,hupu ji ij∈Ni ,α α αi,[λp m]m∈Si= Np−1 ∑ n=0 ℓi(n)+ fi(α α αi)+Λiα α αi,[λp m]m∈Si.(2.11) Here, ¯ upr i and ¯ upu ji describe the private and public input sequences for Np time steps, that is, ¯ upr i=¯upr i(n)Np−1 n=0 and ¯ upu ji =h¯upu ji (n)iNp−1 n=0 , for all j∈Ni . Let us also define α α αi as the aggregation of the scale factors controlled by i: α α αi=hαi,αjij∈Nii⊤.(2.12) Regarding the objective function (2.11), the stage cost ℓi(·)measures the performance of the MPC controller along the prediction horizon, which is described by a quadratic function that will be specified in each chapter of Part II, since it varies according to the type of configuration applied (scenario-based, tube-based, regulation, tracking...). Moreover, 18 Chapter 2. Ceding control authority in input-coupled systems the term fi(·)penalizes the values of the scale factors within α α αisuch that fi(α α αi) = ρprαi+∑ j∈Ni ρpuαji,(2.13) where ρpr and ρpu are positive weighting factors. Finally, the term Λi(·) in (2.11) incorporates the elements that define the dual decomposition algorithm, facilitating the decentralized optimization across the local MPC agents as detailed in Section 2.4. Note that scale factors within α α α and the Lagrange multipliers λm remain constant over the prediction horizon Npbut are recalculated at every time instant k. Accordingly, at each time step k and iteration p , the local MPC optimization problem of each controller i∈Ncan be described as min upr i,hupu ji ij∈Ni ,α α αi Vi(·) s.t. Initial conditions, Local dynamics (2.5), Constraints on the local state, Constraints on the private and public inputs (2.3). (2.14) As will be seen throughout Part II of this thesis, other issues and variables can be taken into account within this setup , possibly modifying some of the basic expressions presented. 3 Stochastic coalitional MPC with negotiation of mutual interactions The results of this chapter, published in [45], introduce a novel coalitional MPC approach designed for systems interconnected via system variables. Here, we implement the decomposition of coupled variables into public and private components, as described in Chapter 2. The idea is to distribute coupling variables sharing a common constraint space between agents. The amount of constraint space for each of these modes is to be decided by agents in a decentralized manner. Additionally, a threshold is defined to activate different interaction modes, including classical decentralized and distributed approaches. The outline of the remainder of this chapter is structured as follows. Section 3.1 introduces the control strategy, describing how flexible modes of cooperation are achieved. Section 3.2 outlines the basic framework for the scenario-based MPC and presents the objective function to be solved by each MPC agent. In Section 3.3, the proposed control scheme is introduced, leveraging the dual decomposition algorithm to negotiate shared variables. Section 3.4 presents simulation results conducted on a simulated eight-input- coupled tanks benchmark. Finally, conclusions are provided in Section 3.5. 3.1 Control strategy This chapter introduces a dual decomposition-based algorithm that exploits the concept of private and public variables described in Chapter 2. Each agent i∈N will optimize the private part of its input variable upr i and the variable αi that bounds its constraints, as well as the public part of neighboring input variables upu ji , and the variable αji that determines its constraints. To avoid a fully cooperative operation, the use of public variables will be penalized with a higher cost than the private ones. Also, to protect the agents against the neighboring uncertainty, we consider scenario-based constraints to reduce the degree of conservativeness. Regarding agents’ communication, we define the threshold ˆ α to trigger between the different interaction modes to permit partial cooperation. In this regard, depending on the 19 20 Chapter 3. Stochastic coalitional MPC with negotiation of mutual interactions values of αi j , agents can switch between different modes of operation ranging between these two extreme cases: • All agents work in a decentralized manner, i.e., there is no communication, which happens when αi j <ˆ αfor all upu i j ∈αi jUi. • All agents work in a distributed manner, i.e., all of them share data if αi j ≥ˆ α for all upu i j ∈αi jUi. Unlike standard distributed schemes, here the boundaries (and not the values) of coupling variables are negotiated. The rationale of this choice is to reduce the cooperation burden. Hence, agents work in a flexible manner that admits partial modes of cooperation. That is, there may be clusters of agents that communicate, i.e., coalitions, whereas others may be operating in a decentralized manner. 3.2 Scenario-based MPC As detailed in Chapter 2, variables that are not locally controlled, that is, the private part of neighboring agents upr j and the public parts {upu i j }j∈Ni , will be considered as bounded disturbances. To take into account these uncertainties, it is possible to follow a conservative robust approach, e.g., tube-based MPC [65, 66]. However, in order to reduce conservativeness, we propose a stochastic approach that generates a fixed number ns of equiprobable random scenarios considering different possible realizations of the disturbance. Let ¯xi be the part of the state that does not depend on the disturbances, i.e., nominal state: ¯xi=Aiixi+Biiupr i+∑ j∈Ni Bi jupu ji .(3.1) Considering ns possible realizations of the uncertainty, the prediction model (2.5) of Chapter 2 can be expressed as x+ i=¯xi+           ∑j∈MiBiiw1 i j +∑j∈NiBi jw1 j ∑j∈MiBiiw2 i j +∑j∈NiBi jw2 j . . .. . . ∑j∈MiBiiwns i j +∑j∈NiBi jwns j. (3.2) As a result, xi aggregates subsystems’ i state at time instant k for all possible scenarios, i.e., xi= [xd i]d∈[1,...,ns] . The control action computed in scenario-based MPC is optimized to account for possible scenarios [67, 68], which can be used to obtain robustness regarding closed-loop constraint satisfaction, e.g., following [69]. 3.2.1 Local objective function According to Chapter 2, we can formulate the objective function of each local controller i as the sum of the stage cost ℓi(·) , penalization of scale factors fi(·) and the corresponding 3.2 Scenario-based MPC 21 terms for the dual decomposition algorithm Λi(·)as Vixi,upr i,hupu ji ij∈Ni ,α α αi,[λp m]m∈Si= Np−1 ∑ n=0 ℓi(n)+ fi(α α αi)+Λiα α αi,[λp m]m∈Si.(3.3) Note that the terms fi(·) and Λi(·) of (3.3) were defined in Chapter 2. Moreover, the performance of each scenario is defined as the following quadratic function ℓd i(n) = ∥xd i(n+1)−xt i(n+1)∥2 Qi+∥upr i(n)∥2 Rpr i +∑ j∈Ni ∥upu ji (n)∥2 Rpu i .(3.4) Weighting matrices Qi , Rpr i , Rpu i are positive definite and satisfy Rpu i>> Rpr i to promote the use of local resources first, and xt i represents the desired reference for subsystem i in time instant k. Accordingly, the stage cost ℓi(·) for each subsystem i is calculated as the weighted sum of the performance cost of each of the scenarios ℓi(n) = ns ∑ d=1 piℓd i(n).(3.5) The weights pi are the probabilities of the occurrence of each scenario, and we assume that all scenarios have the same probability pi=1/ns, thus, ∑∀nspi=1. Therefore, at each time step kand iteration p, agent isolves min upr i,hupu ji ij∈Ni ,α α αi Vi(·) s.t. (3.2), xd i(0) = xi, upr i(n)∈αiUi, upu ji (n)∈αjiUj,∀j∈Ni, xd i(n)∈Xi, ∀d∈[1,2,...,ns], λm≥0,∀m∈Si, ∀n=0,...,Np−1, (3.6) where n is the time index along the prediction horizon Np . As noted in Chapter 2, variables upr i(·) and upu ji (·) vary along the prediction horizon while αi and αji are kept constant in the optimization. The solution of (3.6) provides for each i∈N the value of the private and public inputs to be implemented at instant k , i.e., upr i=upr∗ i(0) and upu ji =upu∗ ji (0) , and the corresponding optimal scale factors, i.e., αi=α∗ iand αji =α∗ ji, for all j∈Ni. 22 Chapter 3. Stochastic coalitional MPC with negotiation of mutual interactions 3.3 Control scheme Considering the above, Algorithm 1 presents a pseudo-code with the steps that have to be followed by the different agents to implement the proposed coalitional MPC scheme. The idea behind it is to reduce communication in comparison with the fully distributed MPC approach. Coordination between agents will be event-based, as communication links will be disabled if Step 2 of the algorithm is verified. We consider that if the public part that an agent can manipulate is small in two consecutive time steps, there is no need for communication. To avoid recomputing an optimization, the implementation assumes that upu ji =0 , but note that an extra step could be included here to re-optimize imposing that the public variable is zero. However, if at a single time step an agent demands a greater use of the public variable, being its value above the threshold ˆ α, communication is retaken. When communication occurs, agents follow a distributed approach according to the dual decomposition algorithm. They share the values of their local variables, updating the corresponding Lagrange multiplier until convergence is attained. It is worth mentioning that agents keep calculating public variables and bounds αji even when they are not communicating with their neighbors in case negotiation should be resumed. Algorithm 1: Control Scheme At each sample time k, each agent i∈Nproceeds as follows: 1: Solve the optimization problem (3.6) for nsequiprobable scenarios and find optimal sequence upr i,upu ji and optimal bounds of the variables αi,αji,∀j∈Ni. 2: if αji(k)<ˆ αand αji(k−1)<ˆ αthen 3: Disable communication flag: fji =0. 4: Ignore the value of the public variable: upu ji =0. 5: else 6: Enable communication flag: fji =1. 7: end if 8: if fji =1then 9: while ∆>εdo 10: Update local and neighboring Lagrange multipliers {λi,λj:∀j∈Ni} according to (2.10). 11: Compute ∆=max{∥αp j−αp−1 j∥,∥αp ji −αp−1 ji ∥} ∀ j∈Ni. 12: Set p←p+1. 13: end while 14: end if 3.4 Simulation Results An eight input-coupled tanks plant has been used as an example to illustrate the behavior of the proposed technique. The system is composed by four top tanks (5, 6, 7, and 8) that discharge into four bottom tanks (1, 2, 3, and 4). In turn, tanks in the lower level 3.4 Simulation Results 23 Figure 3.1 Scheme of the eight tanks system. The colored pipes represent the coupling connections. discharge into a shared storage tank. The eight tanks are filled thanks to four pumps (Qa,Qb,Qcand Qd) that carry water from the storage tank. Moreover, four three-way valves are employed to regulate the pumped flow, dividing it into two ways to fill the upper and lower tanks. We can distinguish N=4 subsystems formed by a top and bottom tank. In this way, tanks 1 and 5 define the first subsystem; the second is composed of tanks 2 and 6 ; tanks 3 and 7 describe the third subsystem; finally, tanks 4 and 8 form the fourth one. These subsystems are physically coupled through the colored pipes that interconnect the tanks in Figure 3.1. Also, the following bidirectional links, i.e., data connections between the corresponding local controllers, are considered: (4,1) , (1,2) , (2,3) , and (3,4) . Here, link (i,j)denotes that local controllers iand jcan exchange information. The objective is to regulate the four lower tanks towards their desired reference in the water level. For this purpose, we define the state of each subsystem as the water level in meters of the two tanks that represent it, e.g., for subsystem one: x1= [h1h5]T . Furthermore, inputs are given by the value of the pump flow in cubic meters per hour ui=Qk , with i=1,...,N and k∈ {a,b,c,d} . Besides, the system is subject to the following constraints: 0.2≤xi≤1.3,[0,0,0,0]≤ui≤[3.26,4,3.6,4],with i=1,...,N . Lastly, each subsystem is characterized by 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#, 30 Chapter 4. Robust coalitional MPC with negotiation of mutual interactions Depending on the values of αi j relative to ˆ α for different combinations of {i,j} , there may be clusters of agents that communicate for the sake of overall performance, while others may operate in a decentralized manner without any communication. Thus, agents operate flexibly, allowing for partial modes of cooperation. This approach enables the reduction of communication and computational burden without significantly compromising the system performance compared to centralized behavior. 4.2 Tube-based MPC for regulation According to the decomposition of coupling variables introduced in Chapter 2, any agent i will optimize the private part of its input variable upr i and the public part of neighboring input variables upu ji for all j∈Ni . Conversely, variables that are not locally controlled, namely, the private variables of the set of neighboring agents Ni , i.e., {upr j}j∈Ni , and the public part of its input variable ui , which is controlled by the set of agents within Mi , i.e., {upu i j }j∈Mi , will be treated as bounded disturbances. In contrast to Chapter 3, where uncertainties were addressed through a stochastic approach, in this chapter, mutual disturbances induced by coupling are considered using a tube-based approach. This way, we ensure robust satisfaction of the system’s constraints, regardless of the disturbance realization. The tube-based MPC approach [65,66], 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. Accordingly, we can derive the nominal model for each subsystem i by ignoring interactions wi, i.e., ¯x+ i=Aii ¯xi+Bii ¯upr i+∑ j∈Ni Bi j ¯upu ji .(4.1) In what follows, we consider the following assumptions: Assumption 4.2.1 For every subsystem i: i. There exists a local feedback gain Ki that guarantees that AKi≜(Aii +BiKi) is stable, where Bi is defined as the aggregation of the local input-to-state matrices, that is, Bi≜Bii {Bi j}j∈Ni. ii. A robust positively invariant (RPI) set Ri can be found such that AKiRi⊕Wi⊆Ri , Ri⊆int(Xi), and KiRi⊆int(Ui). iii. The nominal state and input constraints are non-empty sets: ¯ Xi=/0 and ¯ Ui=/0 , with ¯ Xi≜Xi⊖Riand ¯ Ui≜Ui⊖KiRi. 4.2.1 Local optimization problem The objective function that characterizes the nominal control problem is defined starting from the basic function introduced in Chapter 2. The goal of each local controller i is to 4.2 Tube-based MPC for regulation 31 minimize the objective function Vi¯xi,¯ upr i,h¯ upu ji ij∈Ni ,α α αi,[λp m]m∈Si= Np−1 ∑ n=0 ℓi(n)+ fi(α α αi)+Λiα α αi,[λp m]m∈Si+∥¯xi(Np)∥2 Pi. (4.2) The first term of (4.2) corresponds to the stage cost ℓi(·) , described by the quadratic function: ℓi(n) = ∥¯xi(n)∥2 Qi+∥¯upr i(n)∥2 Rpr i +∑ j∈Ni ∥¯upu ji (n)∥2 Rpu i (4.3) where Qi , Rpr i , Rpu i are positive definite weighting matrices. The use of public variables will be penalized with a higher cost to discourage unnecessary cooperation efforts, that is, Rpu i>> Rpr i . Moreover, the terms fi(·) and Λi(·) were defined in Chapter 2, and the last term of (4.2) represents the terminal cost function, with Pi>0. Consequently, at each time instant kand iteration step p, each agent i∈Nsolves: min ¯xi(0),¯ upr i,h¯ upu ji ij∈Ni ,α α αi Vi(·) s.t. ¯xi(0)∈xi⊕(−Ri), ¯xi(n+1) = Aii ¯xi(n)+Bii ¯upr i(n)+ ∑ j∈Ni Bi j ¯upu ji (n), ¯upr i(n)∈αp i¯ Ui, ¯upu ji (n)∈αp ji ¯ Uj,∀j∈Ni, ¯xi(n)∈¯ Xi, ¯xi(Np)∈¯ Xf i={0}, λp m≥0,∀m∈Si, ∀n=0,...,Np−1. (4.4) where n denotes the time step along the prediction horizon Np . The terminal region of the nominal model is represented by ¯ Xf i , and 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. Solving (4.4) provides for each subsystem i∈N the value of the optimal nominal private and public input sequences, i.e., ¯ upr∗ i and ¯ upu∗ ji , and the corresponding optimal scale factors α α α∗ i , i.e., αi=α∗ i and αji =α∗ ji , for all neighboring agents j∈Ni . Note that there will be |Ni|different sequences ¯ upu∗ ji , one from each neighboring subsystem. 4.2.2 Tube-based control law To deal with the mutual disturbances ignored by the nominal model, an auxiliary control law is used so that the trajectory of the real system follows closely the nominal one, trying 32 Chapter 4. Robust coalitional MPC with negotiation of mutual interactions to cancel the error between the nominal and real states, i.e., ¯xi and xi . Moreover, since the MPC follows a receding horizon strategy, only the first component of the optimized sequences ¯ upr∗ i and ¯ upu∗ ji for all j∈Ni is applied at time instant k for the nominal control law. Consequently, the private and public inputs of each subsystem iare calculated as:   upr i hupu ji ij∈Ni  =  ¯upr∗ i(0) h¯upu∗ ji (0)ij∈Ni  +Ki(xi−¯xi).(4.5) 4.3 Control scheme The proposed coalitional control approach aims to reduce communication among agents and is summarized in Algorithm 2. For notational convenience, in what follows we use the superscript − to refer to a variable at the previous time instant, e.g., x− denotes x at 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 2. The communication link between two agents i and j will 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 k agent i requires greater use of the public variable upu ji , being the value of its scale factor αji above the threshold ˆ α. Agents that are communicating follow an iterative procedure under the distributed dual decomposition algorithm. They share the values of the scale factors αi and α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. The set Rp i has to be recalculated for each iteration p and for every time instant k . Note that, in the literature, there are methods such as [70, 71] that allow us to perform this calculation. In particular, the invariant sets for Algorithm 2 are designed according to [72], 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 r denotes the number of inequalities defining the set and Pis a predefined matrix. Therefore, the update of Rp ientails a light computational load. 4.3.1 Checking step Since input, state, and uncertainty constraint sets vary over time, the RPI set must be recomputed at every instant k to suit the new constraints. Following [54], a checking step is introduced to verify if the difference between the predicted and nominal states in 4.3 Control scheme 33 Algorithm 2: Control Scheme Initialization: At time instant k=0and 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 for all {i,j} ∈ N , with αaux >0 . If k>0 , the value of the scale factor at p=0 is set according to k−1:αp=0 i(k) = αi(k−1),αp=0 ji (k) = αji(k−1)for all i,j∈N. Moreover, set [λm]m∈Si=0at p=0for all k. At each sample time k, each agent i∈Nproceeds as follows: 1: while ∆>εdo 2: Update the disturbance set Wp i(2.6) according to the scale values from the previous iteration, i.e., hαp−1 i j i∀j∈Mi and hαp−1 ji∀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 (4.4) 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 (4.5). 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 (2.10). 17: Compute ∆j=max{|αp j−αp−1 j|,|αp ji −αp−1 ji |}. 18: end if 19: Set p←p+1and compute ∆=max{∆j}∀j∈Ni. 20: end while k+1 belongs to the new RPI calculated at that time instant in order to maintain constraint satisfaction and recursive feasibility. According to [65], if the state belongs to the set ¯xi⊕Ri for 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 i is 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 k the RPI set because the constraints that affect Ri may have changed with respect to the previous instant. According to xi∈¯xi⊕Ri , at time instant k , successor states x+ i and ¯x+ i are calculated using the invariant set Ri . As condition x+ i∈¯x+ i⊕Ri holds, we have to guarantee that the new RPI calculated at k+1 contains the difference between the predicted real and nominal current 34 Chapter 4. Robust coalitional MPC with negotiation of mutual interactions states: x+ i−¯x+ i∈R+ i.(4.6) To pursue this goal, Step 3 of Algorithm 2 is introduced with the same procedure as the checking step introduced in [54] (Subsection 4.1). If (4.6) is met, the following is satisfied: AKix+ i−¯x+ i⊕W+ i⊆R+ i.(4.7) Consequently, the predicted state and input trajectory will verify the original constraints Xi and Ui . Conversely, if x+ i−¯x+ i/∈R+ i we will discard the new RPI set and use the previous one, i.e., Ri:=R− i, since x+ i∈¯x+ i⊕Riis always fulfilled. 4.4 Recursive feasibility and stability In what follows, we proceed as in [54] 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. 4.4.1 Decreasing trend of the disturbance set 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 (4.6). Assumption 4.4.1 The set of disturbances are non-increasing from one time instant k to another k+1, i.e., W+ i⊆Wi. Remark 4.4.1 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 non-increasing evolution of factors αj ( j∈Ni ) and αi j ( j∈Mi ) will lead to coupling disturbances sets satisfying Assumption 4.4.1 (recall (2.5) and (2.6)). 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 4.4.2 The invariant sets are non-increasing from one time instant k to another k+1, i.e., R+ i⊆Ri. Let us consider that at time instant k , subsystem i satisfies: 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 Wi as any uncertainty set contained in Wi , i.e., Wi⊆Wi . As disturbances decrease over time, AKiRi⊕Wi⊆Ri is 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 4.4 Recursive feasibility and stability 35 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. 4.4.2 Conditions for recursive feasibility Theorem 4.4.1 (Recursive feasibility) Consider that Assumptions 4.4.1 and 4.4.2 hold. Then, if at the initial time instant k=0 it is possible to find a feasible solution of problem (4.4) 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 (4.4) at instant k . Additionally, note that Assumptions 4.4.1 and 4.4.2 imply that W+ i⊆Wi and R+ i⊆Ri , with W+ i and R+ i being respectively the disturbances and RPI sets at instant k+1 , and Wi and Ri those 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          ,(4.8) ˜ 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.(4.9) From the constraints of (4.4) , we have that xi−¯x∗ i(0)∈Ri and ¯x∗ i(Np) = {0} . Additionally, if we apply control law (4.5) , 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+ i must 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.3.1, if x+ i−˜x+ i(0)/∈R+ i , then we set R+ i:=Ri as a fail-safe option, i.e., ¯x∗ i(1)will always satisfy the constraint on the initial nominal state of problem (4.4). 36 Chapter 4. Robust coalitional MPC with negotiation of mutual interactions Additionally, given that R+ i⊆Ri , nominal constraints sets ¯ Ui and ¯ Xi can 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∈Ni in (4.8) and (4.9) 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, (4.8) and (4.9) constitutes a feasible solution of problem (4.4) for time step k+1 . By induction, the theorem is proven. ■ 4.4.3 Conditions for stability Theorem 4.4.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 Vi(4.2) decreases during the evolution of the system. Let the optimal sequences in (4.8) ¯ x∗ i , ¯ upr∗ i and ¯ upu∗ ji at k lead to the optimal cost V∗ i . Additionally, the candidate state and input sequence for k+1 , i.e., ˜ x+ i , ˜ upr+ i and ˜ upu+ ji , lead to the cost ˜ V+ i , which provides a value for the upper bound of the optimal cost: ˜ V+ i≥V∗+ i . Moreover, we can calculate the difference ˜ V+ i=V∗ i−∥¯x∗ i(0)∥2 Qi−∥ ¯upr∗ i(n)∥2 Rpr i −∑ j∈Ni ∥¯upu∗ ji (n)∥2 Rpu i ∥˜x+ i(Np−1)∥2 Qi+∥˜upr+ i(Np−1)∥2 Rpr i +∑ j∈Ni ∥˜upu+ ji (Np−1)∥2 Rpu i +˜ f+ i(α α αi)−f∗ i(α α αi)+∥˜x+ i(Np)∥2 Pi−∥¯x∗ i(Np)∥2 Pi +˜ Λ+ iα α αi,[λm]m∈Si−Λ∗ iα α αi,[λm]m∈Si. (4.10) Unlike private and public inputs ¯upr i and ¯upu ji , optimal scale factors in α α α∗ i remain 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 k to the candidate solution in k+1 . Therefore, ˜ f+ i=f∗ i and ˜ Λ+ i=Λ∗ i , so they cancel themselves in (4.10) . As a result, in the subtraction ˜ V+ i−V∗ i there are only left terms referring to the stage and terminal costs. Furthermore, the terminal cost is a continuous Lyapunov function at the origin, which entails: ∥˜x+ i(Np−1)∥2 Qi+∥˜upr+ i(Np−1)∥2 Rpr i +∑ j∈Ni ∥˜upu+ ji (Np−1)∥2 Rpu i +∥˜x+ i(Np)∥2 Pi−∥¯x∗ i(Np)∥2 Pi≤0. (4.11) Therefore, we can rewrite (4.10) as: ˜ V+ i−V∗ i≤ −∥¯x∗ i(0)∥2 Qi−∥ ¯upr∗ i(n)∥2 Rpr i −∑ j∈Ni ∥¯upu∗ ji (n)∥2 Rpu i .(4.12) Identifying the right-hand side of the inequality as ℓi¯x∗ i(0),¯upr∗ i(0),h¯upu∗ ji (0)ij∈Ni,(4.13) 4.5 Simulation results 37 and considering ˜ V+ i≥V∗+ i, (4.12) implies V∗+ i−V∗ i≤ −ℓi¯x∗ i(0),¯upr∗ i(0),h¯upu∗ ji (0)ij∈Ni.(4.14) We can extend the (4.14) for all subsystems i∈N: V=∑ ∀i∈N Vi→V∗+−V∗≤ −ℓ(¯x(0),¯u(0)).(4.15) The stage cost ℓ(·) is strictly positive, therefore, it is proven that the cost function Vi defined in (4.2) decreases over time, ensuring the stability of the system. ■ 4.5 Simulation results In this section, we apply the proposed coalitional control algorithm to the academic example shown in Figure 4.1. 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 Figure 4.1. 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 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 Figure 4.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 i and 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 example is similar to that of Chapter 3, but in this case, we consider the operating point of the plant, that is, in this case 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, (4.16) 38 Chapter 4. Robust coalitional MPC with negotiation of mutual interactions Figure 4.1 Diagram of the eight tanks plant, where the coupling relations are represented by the colored pipes. 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#. (4.17) Accordingly, the disturbance vectors (recall (2.6) ) for the simulation example are such that: w1∈W1≜B11α12U1⊕B14α4U4, w2∈W2≜B22α23U2⊕B21α1U1, w3∈W3≜B33α34U3⊕B32α2U2, w4∈W4≜B44α14U4⊕B43α3U3. (4.18) 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 4.5 Simulation results 39 must be adapted to the state and input constraints by subtracting the operating point (4.16) 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 (4.2) , calculated by trial and error. The prediction horizon has been set to Np=10 , and the simulation lasts 25s. 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 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. In Figures 4.2 and 4.3, solid 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 (4.16) . Figure 4.2 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 Figure 4.3, which illustrates the performance of Algorithm 2 in 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 0 5 10 15 20 25 0.6 0.8 1 1.2 0 5 10 15 20 25 0.6 0.8 1 1.2 0 5 10 15 20 25 0.6 0.8 1 1.2 0 5 10 15 20 25 0.6 0.8 1 1.2 Figure 4.2 Evolution of the water level of the lower tanks. 5 Partially cooperative robust coalitional MPC for tracking This chapter extends our previously proposed robust approach for regulation, as detailed in Chapter 4, to specifically tackle the tracking scenario. Accordingly, we investigate a coalitional MPC setting for input-coupled multi-agent systems to robustly track changing set points, which can compromise the controller’s feasibility and result in tracking failures. To address this, we adopt the tracking formulation presented in [73, 74], and adapt it to the previously mentioned partially cooperative framework. This way, we manage to guide multi-agent systems toward the desired target points while ensuring recursive feasibility and stability through advanced tracking techniques. The content of this chapter is based on [47,48]. The remainder of the chapter is organized as follows. Section 5.1 outlines the control strategy and describes the communication topology among the different MPC agents. In Section 5.2, we describe the problem formulation for tube-based MPC for tracking incorporating the public and private variable decomposition, and defines the local optimization problem of the different agents. In Section 5.3, the proposed control algorithm is presented, and Section 5.4 provides theoretical guarantees of our approach. Finally, Section 5.5 presents simulation results, and concluding remarks are given in Section 5.6. 5.1 Control strategy In this chapter, we introduce a coalitional model predictive control scheme that builds upon the decomposition of coupling variables into public and private components presented in Chapter 2. This decomposition was initially applied in Chapter 3 within the context of a stochastic framework and further extended in Chapter 4 to a robust approach. In contrast to the preceding chapters, the current one delves into a coalitional control setting for input-coupled multi-agent systems, aiming to robustly track changing set points. These techniques have been applied in a robust and stochastic manner [75 – 77] and also in a distributed setting [78–80]. 47 48 Chapter 5. Partially cooperative robust coalitional MPC for tracking In line with the previous chapters, communication between agents is event-based. In Chapters 3 and 4, we defined a threshold ˆ α on the scaling factor αi j , which characterizes the ceded portion of control action that triggers communication. Conversely, in this chapter, the state of the communication links will change based on the evolution of αi j , which scales the constraint set of the public variable upu i j . Specifically, as the derivative of the variable αi j approaches zero, the corresponding communication link is deactivated. 5.2 Tube-based MPC for tracking To robustly track different setpoints (xt i,ut i) , we employ tube-based MPC for tracking [75]. This approach introduces an artificial setpoint as an optimization variable in the nominal MPC problem, extending the domain of attraction of the MPC controller. Additionally, it addresses the deviation between the artificial and desired setpoints by introducing a penalty term in the cost function. As seen in Chapter 4, the tube-based approach aims to manage the bounded disturbances arising from the interaction between subsystems. Recall that the control law is characterized by an optimal input for the nominal system and an auxiliary feedback control law that ensures the evolution of the real system remains in the vicinity of the nominal predictions, defining the nominal model by disregarding coupling interactions: ¯x+ i=Aii ¯xi+Bii ¯upr i+∑ j∈Ni Bi j ¯upu ji .(5.1) We consider the following assumptions, which implicitly includes Assumption 4.2.1: Assumption 5.2.1 For every subsystem i: i. The local state xi is measurable, and the pair (Aii,Bi) is controllable, where Bi is defined as Bi≜Bii {Bi j}j∈Ni. ii. There exists a local feedback gain Ki guaranteeing that AKi≜(Aii +BiKi) is stable. iii. A robust positively invariant (RPI) set Ri can be found such that AKiRi⊕Wi⊆Ri , Ri⊆int(Xi), and KiRi⊆int(Ui). iv. The nominal state and input constraints are non-empty sets: ¯ Xi=/0 and ¯ Ui=/0 , with ¯ Xi≜Xi⊖Riand ¯ Ui≜Ui⊖KiRi. 5.2.1 Local optimization problem Building upon the basic objective function defined in Chapter 2, we incorporate additional terms required for the tracking setting. This way, each local controller i∈N minimizes Vi¯xi,¯xs i,¯ upr i,¯upr,s i,h¯ upu ji ,¯upu,s ji ij∈Ni ,α α αi,[λp m]m∈Si= Np−1 ∑ n=0 ℓi(n)+ fi(α α αi)+Λiα α αi,[λp m]m∈Si+∥¯xs i−xt i∥2 T+∥¯xi(Np)−¯xs i∥2 Pi. (5.2) 5.2 Tube-based MPC for tracking 49 The terms fi(·) and Λi(·) are included according to Chapter 2. Also, ℓi(·) corresponds to the stage cost, which penalizes the deviation between the predicted nominal state of the plant and the artificial steady conditions along the prediction horizon Np , and is described as a quadratic function ℓi(n) = ∥¯xi(n)−¯xs i∥2 Qi+∥¯upr i(n)−¯upr,s i∥2 Rpr i +∑ j∈Ni ∥¯upu ji (n)−¯upu,s ji ∥2 Rpu i .(5.3) Weighting matrices Qi , Rpr i , Rpu i are positive definite, and (¯xs i,¯upr,s i,¯upu,s ji ) correspond to the artificial variables. Since we are applying the variable decomposition proposed in Chapter 2, we need to introduce as many artificial inputs as the number of private and public variables we have, satisfying ¯xs i=Aii ¯xs i+Bii ¯upr,s i+∑j∈NiBi j ¯upu,s ji . These artificial state and inputs are considered as decision variables that are kept constant over the prediction horizon, extending the domain of attraction of the MPC controller and mitigating the potential feasibility issues arising from changes in the setpoint. The considered cost function minimizes both the predicted tracking error with respect to this artificial reference and the offset deviation from the desired setpoint. In this regard, the term ∥¯xs i−xt i∥2 T corresponds to the so-called offset cost, which ensures convergence to a vicinity of the target setpoint penalizing the distance between the artificial reference and the setpoint, with T>0 . Moreover, the terminal cost ∥¯xi(Np)−¯xs i∥2 Pi penalizes the deviation of the predicted trajectory with respect to the artificial reference, with Pi>0 . Consequently, the nominal system converges to an admissible target steady state. If such a state is not admissible, the system is steered towards the nearest admissible steady state. This way, the controller maintains the feasibility of the closed-loop system under any change of the setpoint, ensuring convergence to the admissible setpoints. Considering (5.2) , at every time instant k and iteration step p , each agent i∈N solves min ¯xi(0),¯xs i,¯ upr i,¯upr,s i,h¯ upu ji ,¯upu,s ji ij∈Ni ,α α αi Vi(·) s.t. ¯xi(0)∈xi⊕(−Ri), ¯xi(n+1) = Aii ¯xi(n)+Bii ¯upr i(n)+ ∑ j∈Ni Bi j ¯upu ji (n), ¯upr i(n)∈αp i¯ Ui,¯upu ji (n)∈αp ji ¯ Uj, ¯xi(n)∈¯ Xi, ¯xi(Np) = ¯xs i, ¯upr,s i∈αp i¯ Ui,¯upu,s ji ∈αp ji ¯ Uj, ¯xs i∈¯ Xi, ¯xs i=Aii ¯xs i+Bii ¯upr,s i+∑ j∈Ni Bi j ¯upu,s ji , λp m≥0,∀m∈Si, ∀j∈Ni, ∀n=0,...,Np−1, (5.4) 50 Chapter 5. Partially cooperative robust coalitional MPC for tracking The index n is the time step along the prediction horizon Np . To ensure asymptotic stability, besides the terminal cost function ∥¯xi(Np)−¯xs i∥2 Pi included in (5.2) , a terminal constraint is added, requiring the nominal state at the end of the prediction horizon, i.e, ¯xi(Np) , to match the artificial state for tracking. Furthermore, the artificial state is penalized through the offset cost to coincide with the target setpoint. Consequently, the real state of the system, subject to uncertainties, will remain in the vicinity of the target set defined by Ri . In model predictive control strategies with terminal equality constraint it is generally assumed that, the prediction horizon Np is such that the controllability matrix Ci≜hBiAiiBiA2 iiBi... ANp−1 ii Biiis full rank. 5.2.2 Tube-based control law Following the same considerations as in Chapter 4.2.2, let ¯ upr∗ i and ¯ upu∗ ji denote the optimal sequences of private and public inputs, respectively, which are solution of (5.4) for a given subsystem i . In this regard, we can define the control law as the sum of the optimal nominal input and a linear feedback law as   upr i hupu ji ij∈Ni  =  upr i(0) hupu∗ ji (0)ij∈Ni  +Ki(xi−¯xi).(5.5) 5.3 Control scheme The proposed coalitional control scheme, detailed in Algorithm 3, is designed to minimize communication among agents while ensuring constraint satisfaction and recursive feasibility in tracking different setpoints. For notational convenience, we adopt the convention of using superscripts “ + " and “ − " to denote a variable at the subsequent and preceding time instants, respectively (i.e., k+1and k−1). At the beginning of the simulation, all links (j,i) , where i∈N , j∈Ni , and i=j , are enabled. Then, the communication topology of the system is selected at regular intervals of Ttop time instants. Agents coordinating at a given time instant are engaged in an iterative negotiation, exchanging the values of the scale factors αi and αji until convergence is achieved. In particular, communication between agents follows an eventbased approach, where the communication link (j,i) becomes inactive if the derivative of the scale factor αi j —defining the input constraint set of upu i j —approaches zero, and becomes active otherwise. 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 the time derivative of αji approaches zero, it implies that the system is converging towards the desired reference. Consequently, while maintaining this operating point, the system is not expected to require additional input space. Hence, there is no requirement for agents to engage in negotiation regarding this value, and it will be locally determined by agent i and incorporated into agent j as part of its disturbance, saving communication resources. 5.3 Control scheme 51 Algorithm 3: Control scheme Start: At k=0and p=0, set αp=0 i=α,αp=0 ji =αfor all j∈N, with α>0, and λp=0 m=0for all m∈Si. At every time instant k, each agent i∈Nproceeds as: 1: while ∆>εdo 2: Compute Wp iaccording to (2.6), with the values of αjand αi j at p−1. 3: Calculate the RPI set Rp i. 4: if k>0and xi−xi/∈Rp ithen 5: Set Rp i=R− i. 6: end if 7: Solve (5.4) to obtain upr,∗ i,upu,∗ ji and α α α∗ i. 8: Apply (5.5) to calculate upr iand upu ji for all j∈Ni. 9: if communication link (j,i)is enabled then 10: Update the Lagrange multiplier λp+1 j(2.10). 11: Compute ∆j=max{|αp j−αp−1 j|,|αp ji −αp−1 ji |}. 12: else 13: The Lagrange multiplier will remain constant in p+1, i.e., λp+1 j=λp j. 14: Set ∆j=0. 15: end if 16: Update p←p+1and compute ∆=max{∆j}∀j∈Ni. 17: end while 18: Initialize p=0at k+1with αp=0 i(k+1) = αi(k),αp=0 ji (k+1) = αji(k)∀j∈N, and λp=0 m(k+1) = λm(k)for all m∈Si. 5.3.1 Checking step If, at a given time instant k , the state satisfies xi∈¯xi⊕Ri , then, according to [65], we can ensure that at the subsequent time instant k+1 , the state will satisfy x+ i∈¯x+ i⊕Ri . However, in our proposed approach the RPI set can change from one instant k to another k+1 (i.e., R+ i=Ri ) due to the optimization of scale factors α at each time instant k and iteration p. The latter may result in the non-satisfaction of constraint x+ i−¯x+ i∈R+ i.(5.6) To avoid this and ensure the satisfaction of (5.6) , Step 4 of Algorithm 3 follows the same considerations as in Chapter 4.3.1. We introduce a checking step to determine whether the discrepancy between the predicted and nominal states at k+1 falls within the updated RPI set at that particular time instant, that is, we verify if (5.6) holds. If it does, then we can ensure AKix+ i−¯x+ i⊕W+ i⊆R+ i,(5.7) implying that the predicted state and input trajectory remain within the original constraints Xi and Ui . Conversely, if (5.6) is not met, i.e., x+ i−¯x+ i/∈R+ i , the new RPI set, R+ i , is 52 Chapter 5. Partially cooperative robust coalitional MPC for tracking discarded, setting R+ i=Ri , since x+ i∈¯x+ i⊕Ri always holds because successor states at k+1(x+ iand ¯x+ i) are calculated using the invariant set at k, i.e., Ri. In addition, the initial nominal state at k=0 is optimized with respect to ¯xi(0)∈ xi⊕(−Ri) . Consequently, we ensure x0 i−¯x0 i∈R0 i , where R0 i represents the RPI set at the initial time step. Thus, we initiate Algorithm 3 with a feasible solution. With this checking step, we can ensure constraint satisfaction and preserve recursive feasibility, as further examined in the subsequent section. Remark 5.3.1 Given the potential computational demands associated with recomputing the set Rp i , we propose employing the one-step procedure introduced in [72]. This method efficiently computes a polytopic minimal robust positively invariant (mRPI) set by solving a single LP problem to determine the (P,r) -mRPI, where r represents the number of inequalities defining the set, and P is a predetermined matrix. Consequently, updating Rp i incurs a relatively low computational burden using this method. 5.4 Recursive feasibility and stability This section establishes the conditions for recursive feasibility and stability by extending the proofs outlined in Chapter 4 to the tracking scheme. 5.4.1 Decreasing trend of the disturbance set We introduce a set of assumptions, previously outlined in Chapter 4.4, which reinforce the maintenance of recursive feasibility and stability guarantees when using Step 4 of Algorithm 3, in which the previous RPI set is used if the newly computed set fails to satisfy condition (5.6). Assumption 5.4.1 The set of disturbances are non-increasing from one time instant k to another k+1, i.e., W+ i⊆Wi. Remark 5.4.1 For a constant setpoint, the system’s state and input trajectories will evolve towards the desired reference over time. Consequently, the set of disturbances will become smaller as the system approaches a given setpoint, satisfying Assumption 5.4.1. Assumption 5.4.2 The invariant sets are non-increasing from one time instant k to another k+1, i.e., R+ i⊆Ri. Remark 5.4.2 The construction of the RPI set depends on the definition of the uncertainty set, hence, a reduction in the disturbance set results in a smaller RPI set, giving rise to Assumption 5.4.2. If, at a given time instant k , AKiRi⊕Wi⊆Ri is satisfied by subsystem i , for k+1 condition AKiR+ i⊕W+ i⊆R+ i is also satisfied, since W+ i⊆Wi and R+ i⊆Ri . Let Wi denote any uncertainty set such that Wi⊆Wi . Given the contracting nature of the uncertainty sets, the condition AKiRi⊕Wi⊆Ri is satisfied for any disturbance set Wi contained in Wi . Therefore, if at k+1 : x+ i−¯x+ i/∈R+ i , a fail-safe option is to replace the RPI at k+1 (i.e., R+ i ) with the preceding one, setting R+ i=Ri . Thus, Step 4 of Algorithm 3 constitutes a conservative alternative that complies with the system’s constraints using a larger RPI set if the necessary conditions are not satisfied. 5.4 Recursive feasibility and stability 53 5.4.2 Conditions for recursive feasibility Theorem 5.4.1 (Recursive feasibility) Consider that Assumptions 5.4.1 and 5.4.2 hold. Then, if at the initial time instant k=0 it is possible to find a feasible solution of problem (5.4) 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(0),¯xs∗ i,¯ upr∗ i,¯upr,s∗ i,¯ upu∗ ji ,¯upu,s∗ ji ,α α α∗ i) be the optimal solution of (5.4) at instant k for all j∈Ni . Accordingly, let us define ¯ x∗ i= [¯x∗ i(n)]Np n=0 as the optimal sequence of predicted states. Moreover, let us consider ˜ x+ i=h˜x+⊤ i(0)˜x+⊤ i(1)··· ˜x+⊤ i(Np−1)˜x+⊤ i(Np)i⊤ =h¯x∗⊤ i(1)¯x∗⊤ i(2)··· ¯x∗⊤ i(Np)¯xs∗⊤ ii⊤ ˜ upr+ i=h˜upr+⊤ i(0)˜upr+⊤ i(1)··· ˜upr+⊤ i(Np−2)˜upr+⊤ i(Np−1)i⊤ =h¯upr∗⊤ i(1)¯upr∗⊤ i(2)··· ¯upr∗⊤ i(Np−1)¯upr,s∗⊤ ii⊤, ˜ upu+ ji =h˜upu+⊤ ji (0)˜upu+⊤ ji (1)··· ˜upu+⊤ ji (Np−2)˜upu+⊤ ji (Np−1)i⊤ =h¯upu∗⊤ ji (1)¯upu∗⊤ ji (2)··· ¯upu∗⊤ ji (Np−1)¯upu,s∗⊤ ji i⊤∀j∈Ni, (5.8) and ˜xs+ i=¯xs∗ i,˜upr,s+ i=¯upr,s∗ i,and ˜upu,s+ ji =¯upu,s∗ ji ∀j∈Ni,(5.9) as a candidate solution for instant k+1 . These candidate state and input sequences correspond to the shifted optimal sequences ¯ x∗ i,¯ upr∗ i,[¯ upu∗ ji ]j∈Ni . Admissibility of the optimal sequences entail the admissibility of ˜x+ i(0 : Np−1) , ˜upr+ i(0 : Np−2) , and ˜upu+ ji (0 : Np−2) . In addition, since ¯x∗ i(Np) = ¯xs∗ i , with ¯x∗ i being a steady state for subsystem i characterized by (¯xs∗ i,¯upr,s∗ i,[¯upu,s∗ ji ]j∈Ni) , we can always apply these private and public inputs, ˜upr∗ i(Np−1) = ¯upr∗ i and ˜upu+ ji (Np−1) = ¯upu∗ ji for all j∈Ni , as the terminal input to stay at the artificial state ˜x+ i(Np) = ¯xs∗ i. Furthermore, constraint xi−¯x∗ i(0)∈Ri is imposed for the initial state at instant k on the local control problem (5.4) . Once the control law (5.5) has been applied, it holds that x+ i−¯x∗ i(1)∈Ri . However, the RPI set needs to be recalculated at each time step, since the state, input and disturbance constraint sets vary over time. Therefore, condition x+ i−¯x∗ i(1)∈R+ i must be satisfied at k+1 . Since we have defined ¯x∗ i(1) as the initial state for the candidate solution, i.e., ˜x+ i(0) = ¯x∗ i(1) , we have x+ i−˜x+ i(0)∈Ri . According to the checking step in the RPI sets, if x+ i−˜x+ i(0)/∈R+ i we can establish R+ i=Ri as a fail-safe option since ¯x∗ i(1) will always satisfy the constraint on the initial nominal state of the optimization problem (5.4). As stated in Assumption 5.2.1, the nominal constraints are constructed as ¯ Xi=Xi⊖ Ri and ¯ Ui=Ui⊖KiRi . Moreover, according to Assumption 5.4.2, the RPI set in two 54 Chapter 5. Partially cooperative robust coalitional MPC for tracking consecutive time steps satisfies R+ i⊆Ri . Thus, the nominal constraints sets ¯ Ui and ¯ Xi can only enlarge in two consecutive time steps. Hence, the candidate sequences (5.8) are admissible, and constitute a feasible solution of problem (5.4) for time step k+1 . By induction, recursive feasibility is guaranteed. ■ 5.4.3 Conditions for stability Theorem 5.4.2 (Stability) For all agents i∈N , we can guarantee asymptotic stability. □ Proof. The cost Vi(5.2) is a positive definite function. To prove the decreasing trend of Vi over time we define V∗ i as the optimal nominal cost at time instant k and for a given subsystem i∈N, that is, Vi¯x∗ i,¯xs∗ i,¯ upr∗ i,¯upr,s∗ i,[¯ upu∗ ji ,¯upu,s∗ ji ]j∈Ni,α α α∗ i,[λ∗ m]m∈Si. Also, the candidate solution for instant k+1 , which has been proven to be feasible, to ˜ V+ i , which is defined as Vi˜x+ i,˜xs+ i,˜ upr+ i,˜upr,s+ i,[˜ upu+ ji ,˜upu,s+ ji ]j∈Ni,˜ α α α+ i,[˜ λ+ m]m∈Si. The candidate solution ˜ V+ i provides an upper bound of the optimal cost V∗+ i at k+1 , that is, ˜ V+ i≥V∗+ i. Consequently, we calculate the following subtraction ˜ V+ i=V∗ i−∥¯x∗ i−¯xs∗ i∥2 Qi−∥ ¯upr∗ i(0)−¯upr,s∗ i∥2 Rpr i −∑ j∈Ni ∥¯upu∗ ji (0)−¯upu,s∗ ji ∥2 Rpu i +∥˜x+ i(Np−1)−˜xs+ i∥2 Qi+∥˜upr+ i(Np−1)−˜upr,s+ i∥2 Rpr i +∑ j∈Ni ∥˜upu+ ji (Np−1)−˜upu,s+ ji ∥2 Rpu i +˜ f+ i(˜ α α α+ i)−f∗ i(α α α∗ i)+∥˜xs+ i−xt∗ i∥2 T−∥¯xs i−xt i∥2 T +∥˜x+ i(Np)−˜xs+ i∥2 Pi−∥¯x∗ i(Np)−¯xs∗ i∥2 Pi +˜ Λ+ i˜ α α α+ i,[˜ λ+ m]m∈Si−Λ∗ iα α α∗ i,[λ∗ m]m∈Si. (5.10) We assume variables that do not vary over the prediction horizon are kept constant from the optimal solution at kto the candidate solution at k+1, that is, ˜xs+ i=¯xs∗ i,˜upr,s+ i=¯upr,s∗ i,˜upu,s+ ji =¯upu,s∗ ji ,˜ α α α+ i=α α α∗ i,˜ λ+ m=λ∗ m,(5.11) for all j∈Ni and m∈Si . In this way, the terms corresponding to the offset cost, penalization of scale factors (˜ f+ i,f∗ i) and Lagrange multipliers (˜ Λ+ i,Λ∗ i) cancel themselves in 5.5 Simulation results 55 (5.10). Moreover, the terminal cost satisfies ∥˜x+ i(Np−1)−˜xs+ i∥2 Qi+∥˜upr+ i(Np−1)−˜upr,s+ i∥2 Rpr i +∑ j∈Ni ∥˜upu+ ji (Np−1)−˜upu,s+ ji ∥2 Rpu i +∥˜x+ i(Np)−˜xs+ i∥2 Pi−∥¯x∗ i(Np)−¯xs∗ i∥2 Pi≤0. (5.12) As a result, we can express the subtraction (5.10) as ˜ V+ i−V∗ i≤ −∥¯x∗ i−¯xs∗ i∥2 Qi−∥ ¯upr∗ i(0)−¯upr,s∗ i∥2 Rpr i −∑ j∈Ni ∥¯upu∗ ji (0)−¯upu,s∗ ji ∥2 Rpu i ,(5.13) and, therefore, V∗+ i−V∗ i≤ −∥¯x∗ i−¯xs∗ i∥2 Qi.(5.14) Since (5.14) holds for every i∈N , we can express the subtraction for the general system by aggregating V=∑∀i∈NVi, that is, V∗+−V∗≤ −∥¯x∗−¯xs∗∥2 Q.(5.15) Consequently, it is proven that the cost function Vi(5.2) decreases over time, since the value of ∥¯x∗−¯xs∗∥2 Q is strictly positive. According to (5.15) , as time tends to infinity, we can infer that the nominal state of the system tends towards the artificial steady state and the cost function V∗converges to a constant. Considering previous studies on MPC for tracking such as [75, Lemma 1] and [81, Lemma 4 and 5], it can be demonstrated that the artificial steady state converges to a fixed point. For every subsystem i∈N , the nominal system state at infinity, denoted as ¯x∞ i , converges to a steady state ¯xs∗,∞ i . Additionally, based on the concept of tubes, it is ensured that the real state xi converges to a neighborhood of the nominal system defined by Ri , expressed as x∞ i∈¯xs∗,∞ i⊕Ri . In particular in [75, 76] it is established that this final operating point xs∗,∞ i is determined as the minimizer of the offset cost function, i.e., ∥¯xs i−xt i∥2 T, subject to the constraints specified in Equation (5.4). ■ 5.5 Simulation results In this section, we illustrate our coalitional approach through simulation results on an eight input-coupled tanks benchmark, which was the same benchmark used in Chapter 4. In summary, the system consists of four upper tanks feeding into four lower tanks that discharge into a common storage tank. When decoupling the system, we can identify N=4 input-coupled subsystems, each composed of a top and bottom tank. Figure 5.4 illustrates the system’s layout, depicting the physical couplings between subsystems with colored pipes connecting the corresponding tanks. On the other hand, data connections between the local controllers of physically coupled subsystems are denoted as communication links. Specifically, the bidirectional data connections between agents j and i are referred to as link (j,i) , and we can distinguish links (4,1) , (2,3) , (1,2) , and (3,4) . When a particular link is active, the connected agents share information. The state of each subsystem is defined in terms of the water level measured from the corresponding operating point. For instance, considering the first subsystem, its state is 6 Existence of the maximum jointly invariant set of input-coupled systems As seen in preceding chapters, the variable decomposition introduced in Chapter 2 allows agents to partition coupling variables into public and private components. Control authority is ceded in terms of the public portion that a given agent shares of its input with neighboring agents, which is particularly advantageous to enable partial cooperation. Within this framework, we study the maximal portion of each agent’s input space that must remain private to ensure the reachability of the state constraint set. This maximum private portion depends on both the portion ceded by the agent and the portion received from its neighbors. To achieve this, we introduce a set-theoretical characterization that establishes sufficient conditions for the non-existence of the maximum robust control invariant set. The research of this chapter is also related to [82], which considers a constrained autonomous system subject to additive uncertainties bounded through a scaled disturbance set. In particular, the authors characterize the bounds of the scale of the disturbance set that renders impossible to find a robust positive invariant set, and therefore, to guarantee constraint satisfaction. Moreover, [83] proposes a set-theoretical analysis in the context of actuation attacks, setting a condition on the size of the attacker’s input set so that the system cannot be robustly defended. The framework of the latter is taken as a reference point for developing the conditions presented here. In contrast to these previous results, here we study a multi-agent input-coupled system operating under the previously mentioned partially cooperative framework of Chapter 2, and provide limits of the private portion of the input space for which the non-existence of an invariant set is guaranteed. This way, agents can be provided with limits regarding the extent of their partial cooperation. A similar problem is studied in [84], where the concept of a distributed maximum robust positively invariant set is introduced, which may be of interest for works such as [85], where authors propose a clustering method based on the volume of the robust positively invariant set that a coalition of agents can attain. Finally, the content of this chapter has 63 64 Chapter 6. Existence of the maximum jointly invariant set of input-coupled systems led to two articles. In [52], we presented an analytical expression for an upper bound of the portion that an agent can cede of its input variable to another agent within a two-agent scenario. In [53], we expand the scope of analysis to encompass a generic multi-agent system. The remainder of the chapter is organized as follows. Section 6.1 introduces the problem setting for a generic input-coupled multi-agent system. In Section 6.2, definitions and previous results of reachability and invariant sets are given. Section 6.3 develops sufficient conditions on the portion that agents must keep private of their input in order to guarantee the feasibility of the global problem. Section 6.4 illustrates the main results of this chapter through different examples. Finally, conclusions are provided in Section 6.5. 6.1 Problem setting Consider a discrete-time, linear time-invariant system that can be partitioned into a set N={1,2,...,N}of input-coupled subsystems, whose dynamics can be modeled as x+ i=Aiixi+Biiui+di,with di=∑ j∈Ni Bi juj,(6.1) where xi∈Rnxi and ui∈Rnui are the state and control input of subsystem i∈N , and di∈Rnxi describes the input coupling with its set of neighboring agents, defined as Ni= {j∈N:Bi j =0,j=i} . Moreover, each subsystem i∈N is subject to state and input constraints, that is, xi∈Xiand ui∈Ui, for all k≥0. Assumption 6.1.1 The constraint sets Xi and Ui are PC-sets, and the input constraint set is assumed to be symmetrical about the origin, i.e., Ui=−Ui , for every subsystem i∈N . Remark 6.1.1 The symmetry of Assumption 6.1.1 can be easily satisfied in situations where norm bounds are imposed performing a change of variable to reset the origin. 6.1.1 Sharing coupling variables Let us consider the variable decomposition proposed in Chapter 2, where a local variable ui is decomposed into a private part upr i , which controlled exclusively by agent i , and a set of public parts upu i j where each is controlled by agent j∈Mi , where Mi denotes the set of agents affected by subsystem’s i input. That is, agent i will cede a portion of ui to those agents belonging to Mi. Accordingly, we can express the local input uias ui=upr i+∑ j∈Mi upu i j ,∀i∈N.(6.2) These portions are bounded by different scale factors that satisfy: upr i∈αiUi , with αi∈[0,1] , and upu i j ∈αi jUi , with αi j ∈[0,1] . In this regard, to fulfill the input constraint ui∈Ui , the following inequality must be met for all i∈N αi+∑ j∈Mi αi j ≤1.(6.3) 6.2 Definitions and previous results 65 Considering the decomposition in (6.2) we can rewrite the model (6.1) of every subsystem i∈Nas x+ i=Aiixi+Biiupr i+∑ j∈Mi upu i j +∑ j∈Ni Bi jupr j+∑ h∈Mj upu jh.(6.4) Let us define Mi j=Mj\{i} as the set of agents affected by subsystem’s j input, except agent i, so that (6.4) can be rewritten as x+ i=Aiixi+Biiupr i+∑ j∈Mi upu i j +∑ j∈Ni Bi jupr j+upu ji +∑ l∈Mi j upu jl .(6.5) We can rewrite the previous expression by separating the variables locally controlled by agent ifrom the disturbances induced by coupling as x+ i=Aiixi+Biiupr i+∑ j∈Ni Bi jupu ji +∑ j∈Mi Biiupu i j +∑ j∈Ni Bi jupr j+∑ l∈Mi j upu jl  | {z } wi .(6.6) The term wi groups the uncertainties induced by the interaction between subsystems, and is bounded by wi∈Wias follows Wi≜M j∈Mi Biiαi jUi⊕M j∈Ni Bi jUjαj+∑ l∈Mi j αjl,(6.7) where Wi is a C-set for all i∈N , as it is a mapped set derived from the input constraint sets (see Assumption 6.1.1). Assumption 6.1.2 The pair (Aii,Bii)is controllable for every i∈N. The objective of this chapter is to characterize a bound on the maximal portion αi that agent i can keep private of its input variable, considering that it is going to gain control authority of a portion of neighboring agents’ input Ni , and cede part of its input to agents within the set Mi . This characterization is performed by examining the non-existence of a jointly invariant set for which agents cannot determine control laws for partial cooperation capable of robustly ensure the infinite reachability of their constraint sets Xi for all i∈N . 6.2 Definitions and previous results This section introduces some relevant definitions and previous findings regarding the reachability of a given set, which is a key concept for the development of the results of this chapter. For this purpose, consider a general linear system: x+=Ax +Bv +Ed,(6.8) 66 Chapter 6. Existence of the maximum jointly invariant set of input-coupled systems subject to the constraints x∈X , v∈V and d∈D . Our goal is to maintain the system within its state constraint set X , applying an admissible input v∈V , regardless of the realization of the uncertainty d∈D. 6.2.1 Set invariance theory In this section, we delve into the study of the reachability of the set X [86 – 88]. A set Ω is considered infinitely reachable when there exists a control law v(·)and some initial state x(0)∈Ω such that, regardless of the uncertainty’s realization d∈D , the state satisfies x∈Ω , with v∈V , for all k>0 . Likewise, a set Ω is deemed strongly reachable if a control law v(·) exists that ensures, for all initial states x(0)∈Ω and under any possible uncertainty realization d∈D , the state satisfies x∈Ω , with v∈V , for all k>0 . In this case, Ω is a robust control invariant (RCI) set, ensuring that the system’s evolution will remain inside it and guaranteeing constraint satisfaction. Therefore, the evolution of the system is admissible if there exists an RCI set Ω⊆X and an initial state x(0)∈Ω , so that x∈Ω⊆X∀k>0 . Notice that both definitions are closely related, as a set is considered infinitely reachable if it contains a strongly reachable (or RCI) set. The largest strongly reachable set of a given set Ω is denoted as the maximal robust control invariant set. It is said to be maximal because it contains all the robust control invariant sets contained in Ω . The computation of this set stems from the k -step robust constraint admissible set Ck(Ω) , which is the set of all initial states that can be kept within the set Ωfor at least ksteps, fulfilling the input constraints v∈V, and for any realization of the disturbance d∈D. The set Ck(Ω)satisfies i. Ck+1(Ω)⊆Ck(Ω), with C0(Ω) = Ω, ii. Ck(Ω) = ∩k n=0Cn(Ω), iii. C∞(Ω):=limk→∞Ck(Ω) = ∩∞ n=0Cn(Ω), where C∞(Ω) corresponds to the maximal RCI set, which is finitely determined if, for a certain time step k , Ck+1(Ω) = Ck(Ω)=/0 , and therefore C∞(Ω) = Ck(Ω) . However, if, for a given k , Ck(Ω) = /0 , then C∞(Ω) = /0 . Since we are studying the reachability of the state constraint set X , let us focus on the existence of the set Ck(X) , which can be computed as Ck+1(X) = Q(Ck(X))∩X,with C0(X) = X,(6.9) where Q(·)is the robust one-step set [89]: Q(Φ)≜{x:∃v∈V|Ax +Bv ⊕ED∈Φ}= (A)−1((Φ⊖ED)⊕(−BV)).(6.10) Note that the term (A)−1(·) represents the preimage of the linear transformation A(·) , and exists even if the matrix Ais not invertible. Therefore, we can rewrite (6.9) as Ck+1(X) = (A)−1((Ck(X)⊖ED)⊕(−BV))∩X.(6.11) 6.2 Definitions and previous results 67 6.2.2 Closely-related known results In this section, we summarize the principal findings of [83], which is the basis for the results presented in this chapter. This previous work characterizes the existence of the maximal RCI set by analyzing the sequence of sets {Ck(X)} , since, as mentioned in the properties of the maximal RCI set, if the set Ck(X) becomes empty for a given k , then C∞(X) = /0. In particular, we consider the following propositions: •Proposition 9. The set Ck(X)is bounded as Ck(X)⊆ k \ n=0 (An)−1Tn(X),(6.12) where Tk+1(X) = (Tk(X)⊖AkED)⊕Ak(−BV),with T0(X) = X.(6.13) Proof of Proposition 9. We can compute the recursion, starting from k=0: C0=X= (A0)−1T0, C1=A−1((X⊖ED)⊕(−BV))∩X =A−1((T0⊖ED)⊕(−BV))∩T0 =A−1T1∩A−0T0. For simplicity, consider D≜EDand V≜(−BV): C2=A−1((C1⊖D)⊕V)∩X =A−1[(A−1T1∩A−0T0)⊖D])⊕V∩X =A−1[(A−1T1⊖D)∩(A−0T0⊖D)]⊕V∩X ⊆A−1(A−1T1⊖D)⊕V∩(A−0T0⊖D)⊕V∩X =A−2(T1⊖AD⊕AV)∩A−1(T0⊖D⊕V)∩T0 C2⊆A−2T2∩A−1T1∩A−0T0. Accordingly, the same recursion holds for k>2.■ • Proposition 12. If Tk∗(X) = /0 for some k∗>0 , then Ck(X) = /0 for all k≥k∗ . Proof of Proposition 12. Based on Proposition 9, if, for some k∗:Tk∗(X) = /0 , then Ck(X) = /0 for all k≥k∗, and therefore C∞(X) = /0.■ • Proposition 13. If Sk∗(X) = /0 for some k∗>0 , then Tk(X) = /0 for all k≥k∗ , where Sk(X)≜X⊕"k−2 M n=0 An(−BV)#⊖"k−1 M n=0 AnED#.(6.14) 68 Chapter 6. Existence of the maximum jointly invariant set of input-coupled systems Proof of Proposition 13. Applying (6.13) recursively from k=0we obtain: Tk=X⊖A0ED⊕A0B(−V) ⊖A1ED⊕A1B(−V) . . . ⊖Ak−1ED⊕Ak−1B(−V). Using the property (A⊖B)⊕C⊆(A⊕C)⊖Bfor non-empty sets, we have, for all k≥0, Ti−1⊆X⊕"k−2 M n=0 AnB(−V)#⊖"k−2 M n=0 AnED#. To determine Tk from Tk−1 we first have to subtract Ak−1ED . Then, Tk∗←→ Tk∗−1⊖ Ak∗−1ED=/0 for some k∗≥0. This is true if X⊕"k−2 M n=0 An(−BV)#⊖"k−2 M n=0 AnED#⊖Ak−1ED=/0. ■ Given D≜αV , of the system (6.8) , the main objective of [83] was to determine the smallest value of the scale factor α , denoted as αk∗ , for any step number k∗ , such that Ck∗(X) becomes empty: αk∗≜inf{α:Ck∗(X) = /0,α∈[0,1]} . This was achieved by identifying a critical value of α , denoted as ¯ αk∗ , from which, for any k≥k∗ , Sk(X) becomes empty, thus guaranteeing the emptiness of the set Ck(X) . This result establishes sufficient conditions under which the state constraint set is not infinitely reachable. Remark 6.2.1 The property (A∩B)⊕C⊆(A⊕C)∩(B⊕C) is used when deriving (6.12) using (6.11) and (6.9) , which means that the Minkowski sum is not distributive under set intersections (see [83, Remark 10] for details). As a consequence, (6.12) does not hold with strict equality, and the theoretical bound will satisfy ¯ αk∗≥αk∗ . Therefore, the admissible set Ck(X) may be empty for a value of α below the value determined by the theoretical bound ¯ αk∗, since the inclusion is not tight. 6.3 Limits of the maximum private portion of the input set This section presents the main results of this chapter, where we extend [83] to a distributed setting with multiple agents, using the decomposition in (6.6) . In particular, we study N input-coupled subsystems and define the limits of the input portion that agents must keep private in a partially cooperative framework. That is, agents must keep their states within their constraint sets Xi for any possible realization of neighboring inputs, considering that they can use the entire portion of the input that has been ceded to them, satisfying the input constraints. To this end, we can adapt the definitions of Section 6.2 for our particular 6.3 Limits of the maximum private portion of the input set 69 distributed setting, obtaining an expression of Ck , Tk , and Sk for each subsystem i∈N by identifying EiDi≡Wi, with Widefined in (6.7), and BiVi≡BiiαiUi⊕M j∈Ni Bi jαjiUj.(6.15) In what follows, we include the superindex ‘ αi ’ to identify the particularized sets of equations (6.11) , (6.13) and (6.14) . Hereafter, we proceed as in [83] to develop the conditions on the scale factor αi that bounds the private input constraint set so that Cαi k(Xi) becomes empty for some k , which implies Cαi ∞(Xi) = /0 . Note that, within our partially cooperative framework, we aim to achieve constraint satisfaction at the global system level. Consequently, if any of the sets Cαi ∞(Xi) were to become empty for a given subsystem, we would no longer be able to ensure constraint satisfaction for every agent i∈N. Our objective is to find, for every k∗∈N , the smallest scale factor αk∗ i≜inf{αi: Cαi k∗(Xi) = /0,αi∈[0,1]} . Consequently, we can define, for every k∗∈N , a set Fk∗ containing all the scale factors that have a value below the actual bounds αk∗ i . However, the values of these actual bounds need to be calculated by exhaustive search so that the set Cαi k(Xi) is empty for a given value of k . Although this may not seem critical in small systems, this chapter focuses on multi-agent systems, where the computation of these actual bounds can result in a significant computational burden. To overcome this, we define a theoretical bound ¯ αk∗ i for every subsystem i∈N that can be calculated by simply evaluating an expression. Accordingly, let ¯ Fk∗ be the set containing all the scale factors that have a value below the theoretical bounds {¯ αk∗ i}∀i∈N , for every k∗∈N . For this purpose, we will study the conditions on αi that guarantee the emptiness of the set Sαi k(Xi) , with Xi the local state constraint set for each subsystem i∈N. The sets Fk∗ and ¯ Fk∗ are defined according to the values of αi that imply the emptiness of the sets Cαi k(Xi) and Sαi k(Xi) , respectively. However, since the inclusion of (6.12) is not tight, for {αi}i∈N/∈¯ Fk∗ we can guarantee the non-existence of a RCI set, while for every {αi}i∈N∈Fk∗ we can guarantee the existence of a RCI set. To compute the theoretical bounds ¯ αk∗ i , a series of assumptions and definitions are needed, which are described below. Assumption 6.3.1 The dominant eigenvalue of Aii , say ρi , is real and positive for every subsystem i∈N. Remark 6.3.1 Assumption 6.3.1 is fulfilled by several types of systems, such as positive systems, systems with triangular A-matrices having real values along their diagonal, and first-order systems. Let us define the set Vi≜{v1 i,...,vnv i,−v1 i,...,−vnv i} as the nv linearly independent eigenvectors related to the dominant eigenvalue ρiand their opposites. Also, let ¯vi≜arg min vi∈V hXi(vi) hBiiUi(vi).(6.16) Assumption 6.3.2 We can find at least a non-zero support of the mapped set BiiUi in one of the directions of Vi . We also assume the support of Bi jUj to be non-zero in the direction ¯vi . 70 Chapter 6. Existence of the maximum jointly invariant set of input-coupled systems Considering the previous expressions, we can obtain, for any k∗∈N , and for every i∈Na theoretical bound ¯ αk∗ iso that Sαi k(Xi) = /0, if k≥k∗. Theorem 6.3.1 If, for some k∗∈N , αi>¯ αk∗ i , with ¯ αk∗ i<1 , then Cαi k(Xi) = /0 for all k≥k∗ . The value of ¯ αk∗ idepends on the dominant eigenvalue as follows •If ρi=1 ¯ αk∗ i≜hXi(¯vi)(1−ρi)+∑j∈NihBi jUj(¯vi)Ψi j +hBiiUi(¯vi)ξi hBiiUi(¯vi)ρk∗−1 i−1, with Ψi j =αji1−ρk∗−1 i+αj+∑ l∈Mi j αjlρk∗ i−1, ξi=ρk∗ i−1∑ j∈Mi αi j. (6.17) •If ρi=1 ¯ αk∗ i≜hXi(¯vi)+∑j∈NihBi jUj(¯vi)Ψi j −hBiiUi(¯vi)ξi hBiiUi(¯vi)(1−k∗), with Ψi j =αji(k∗−1)−k∗αj+∑ l∈Mi j αjl,ξi=k∗∑ j∈Mi αi j. (6.18) Proof. Let us consider (6.14) for a general subsystem (6.6) , with EiDi and BiVi defined in (6.7) and (6.15) respectively Sαi k(Xi) =Xi⊕"k−2 M n=0 An iiαiBiiUi#⊕"k−2 M n=0 M j∈Ni An iiαjiBi jUj!# ⊖"k−1 M n=0 M j∈Mi An iiαi jBiiUi!#⊖  k−1 M n=0 M j∈Ni An iiαj+∑ l∈Mi j αjlBi jUj  . Rearranging the above terms, and taking support functions in the direction of ¯vi , we can express the previous expression as hSαi k (¯vi)≤hXi(¯vi)+αi k−2 ∑ n=0 hAn iiBiiUi(¯vi) + ∑ j∈Ni αji k−2 ∑ n=0 hAn iiBi j Uj(¯vi)! −∑ j∈Mi αi j k−1 ∑ n=0 hAn iiBiiUi(¯vi)!−∑ j∈Niαj+∑ l∈Mi j αjl k−1 ∑ n=0 hAn iiBi j Uj(¯vi)!. By definition, hMΩ(η) = hΩ(M⊤η) , for a set Ω and matrix M . Since ¯vi is the eigenvector of Aii corresponding to the dominant eigenvalue ρi , it holds that Aj ii ¯vi=ρj i¯vi . Also, as 6.3 Limits of the maximum private portion of the input set 71 ρi>0, we can express hΩ(ρiη) = ρihΩ(η). hSαi k (¯vi)≤hXi(¯vi)+ k−2 ∑ n=0 ρn ihBiiUi(¯vi)αi−∑ j∈Mi αi j−ρk−1 ihBiiUi(¯vi)∑ j∈Mi αi j + k−2 ∑ n=0 ρn i∑ j∈Ni hBi jUj(¯vi)αji −αj−∑ l∈Mi j αjl−ρk−1 i∑ j∈Ni hBi jUj(¯vi)αj+∑ l∈Mi j αjl. For ρi=1, the geometric series ∑k−2 n=0ρn i=1−ρk−1 i 1−ρi hSαi k (¯vi)≤hXi(¯vi)+ 1−ρk−1 i 1−ρi hBiiUi(¯vi)αi−∑ j∈Mi αi j−ρk−1 ihBiiUi(¯vi)∑ j∈Mi αi j +1−ρk−1 i 1−ρi ∑ j∈Ni hBi jUj(¯vi)αji −αj−∑ l∈Mi j αjl−ρk−1 i∑ j∈Ni hBi jUj(¯vi)αj+∑ l∈Mi j αjl. If hSαi k∗(¯vi)<0then Sαi k∗(Xi) = /0. Hence, solving for αi: αi>hXi(¯vi)(1−ρi)+∑j∈NihBi jUj(¯vi)Ψi j +hBiiUi(¯vi)ξi hBiiUi(¯vi)ρk∗−1 i−1, with Ψi j =αji1−ρk∗−1 i+αj+∑ l∈Mi j αjlρk∗ i−1,ξi=ρk∗ i−1∑ j∈Mi αi j. Conversely, for ρi=1:∑k−2 n=0ρn i=k−1 hSαi k (¯vi)≤hXi(¯vi)+(k−1)hBiiUi(¯vi)αi−∑ j∈Mi αi j−hBiiUi(¯vi)∑ j∈Mi αi j +(k−1)∑ j∈Ni hBi jUj(¯vi)αji −αj−∑ l∈Mi j αjl−∑ j∈Ni hBi jUj(¯vi)αj+∑ l∈Mi j αjl. If hSαi k∗(¯vi)<0then Sαi k∗(Xi) = /0. Rearranging for αi: αi>hXi(¯vi)+∑j∈NihBi jUj(¯vi)Ψi j −hBiiUi(¯vi)ξi hBiiUi(¯vi)(1−k∗), with Ψi j =αji(k∗−1)−k∗αj+∑ l∈Mi j αjl,ξi=k∗∑ j∈Mi αi j. ■ In addition, it is possible to determine the theoretical bound ¯ α∞ i for which the emptiness of the maximal RCI set is guaranteed, i.e., Cαi ∞(Xi) = /0 , by taking limits when k∗→∞ of (6.17) and (6.18): 78 Chapter 6. Existence of the maximum jointly invariant set of input-coupled systems -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 Figure 6.5 Representation of the k -step constraint admissible sets Cαi k(Xi) , for the three subsystems of E1with k<10.. Focusing on the second example, E2 , and the third subsystem, which cedes control authority to the first subsystem through α31 , and receives a portion of subsystem’s two input through α23 , Figure 6.6 shows a comparison of the maximal RCI set of the third subsystem, Cαi ∞(X3) , for different pairs of {α23,α31} . The purpose of this figure is to emphasize the flexibility introduced through the decomposition of (6.6) concerning (6.1) . Note that the setting of (6.1) can be regarded as a special case of our decomposition when agents keep their entire input private without sharing a portion with those agents affected by them. That is, for the third subsystem, having {α23 =0,α31 =0} corresponds to the decentralized scenario. The resulting maximal RCI set for this case has been represented in Figure 6.6 with a white set delimited by solid lines. Our proposed decomposition allows for a dynamic adjustment of the size of the local RCI sets according to the system’s needs. See, for example, when {α23 =0,α31 = 0.17} (depicted with a green set in Figure 6.6), which represents an unfavorable case from subsystem’s 3 perspective. Here, control authority is ceded without receiving any neighboring input portion in return. The value of α31 =0.17 has been selected based on Figure 6.4b, representing a value close to the limit beyond which the existence of an invariant set for subsystem 3 cannot be guaranteed when it receives nothing from its neighbors. In this case, the cost of ceding a portion of its input is a decrease in the size of the local maximal RCI set, resulting in a 28.65% decrease of the Chebyshev radius compared to the decentralized case (white set). However, if {α23 =0,α31 =0} , the maximal RCI set of subsystem 3 will increase concerning the decentralized scenario, as it receives a portion of neighboring inputs while keeping its entire input private. This is represented by the blue and gray sets in Figure 6.6, resulting in increases of 13% and 17.86% in the Chebyshev radius compared to the the decentralized case, respectively. In particular, the blue set shows an intermediate setting where α23 =0.26 and α31 =0 , while the gray set corresponds to the extreme case where the neighboring agent cedes their whole input ( α23 =1 ) while α31 =0 . This hypothetical scenario, where an agent detaches from the entirety of its control action, shows how much could be potentially gained locally with our decomposition. 6.5 Conclusions 79 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Figure 6.6 Representation of the maximal RCI set for the third subsystem, Cαi ∞(X3) , for different values of {α23,α31}. 6.5 Conclusions In this chapter, we establish an upper bound on the portion an agent can keep private of its input variable, assuming that it is sharing a portion of its input to other agents, and that it will receive a portion of neighboring inputs. This bound has been characterized analytically so that we are able to find a jointly robust control invariant set for each subsystem. Our proposed results can also be of interest in some applications such as cybersecurity [90], where it is interesting to find the set of states for which it is possible to find a control law that keeps the system’s state within its states constraint set for a given number of time steps. The calculation of this type of sets is also interesting in coalitional control [41], since the partition of the system is periodically recalculated based on the dynamic interaction between subsystems. However, our results are based on several assumptions that might limit their applicability. Firstly, the system is required to be linear, although the extension to a nonlinear case might be achieved leveraging the robust nature of the approach to account for modeling errors. Also, the largest eigenvalue of the system is required to be real, but this condition is fulfilled naturally by some type of systems. Likewise, constraints must be expressed symmetrically with respect to the origin –which can be attained for most practical situations with a change of variable– and subsystem dynamics must be expressed using only input coupling. Part III Dynamic resource-sharing mechanisms based on coalitional MPC 81 7 Decentralized optimization in constraint-coupled systems In Part III of this PhD dissertation, we focus on developing coalitional resource-allocation mechanisms to dynamically adjust coupling constraints in systems such as the PTC solar fields within the OCONTSOLAR project. Specifically, PTC solar fields are a subtype of concentrating solar power systems designed to concentrate solar radiation [91,92]. This concentrated radiation heats a fluid, typically producing steam, which in turn drives a turbine to generate electrical power [93, 94]. In contrast to other power generation processes where the primary energy source can be handled as a control variable, solar energy cannot be manipulated and thus acts as a disturbance from a control perspective. Consequently, MPC stands out as one of the most widely used techniques for controlling PTC plants [95,96], as it effectively manages system disturbances and constraints in a receding-horizon fashion. However, solar fields are LSS, resulting in a considerable computational burden due to their numerous optimization variables. For this reason, solving the optimization problem in a centralized manner may require significant computational resources for real-time control. To address this challenge, recent research has aimed to reduce the computational burden associated with implementing centralized MPC in PTC solar fields [97 – 99]. In this context, various coalitional control schemes have been proposed to decompose the global control problem into multiple local subproblems [35 – 39], which are then solved in a decentralized manner by agents within a coalition. However, the collectors in PTC solar plants share a common resource, namely the total amount of heat transfer fluid that can be provided to the field. This common coupling constraint presents a challenge in partitioning the global control problem, particularly if we aim to satisfy the constraint in a decentralized manner with partial cooperation. The crucial challenge of how to decouple each coalition’s local problem is addressed in Part III of this thesis, where we develop various coalitional control strategies to distribute the coupled constraint. Such a decentralized framework seeks to reduce the computational burden of the overall approach by taking advantage of parallel executions. Unlike classical decentralized schemes, this is achieved by a dynamic 83 84 Chapter 7. Decentralized optimization in constraint-coupled systems distribution of the constraint across all local probems so as to minimize performance losses. 7.1 Description of the ACUREX solar field Throughout Part III of this thesis, the ACUREX parabolic-trough collector solar field located in Plataforma Solar de Almería (PSA), Spain, is used as case study. The solar field consists of a set of parabolic mirrors, known as collectors, which concentrate solar irradiance in a pipe located in its focal line, where a heat transfer fluid (HTF) heats up as it circulates. The HTF carries the thermal energy to produce steam, and the electrical energy is generated by a steam turbine. For this plant, Therminol 55 is used as the HTF, whose density (ρf) and specific heat capacity (cf) are temperature dependent and are given by ρfTf=903 −0.672Tf,cfTf=1820 +3.478Tf,(7.1) where Tfis the temperature of the HTF [100]. The distributed collector field can be modeled as a set N={1,..., N} of N parallel loops; in particular, the ACUREX solar collector field consists of 480 east-west aligned single-axis tracking collectors that form N=10 loops. Each loop measures 174 m in length, and consists of four 12-module collectors connected in series. For each loop i∈N , it is possible to distinguish an active part that receives solar irradiance ( 144 m length) from a passive part where solar radiation does not reach ( 30 m length). Moreover, it is assumed that the HTF inflow of each loop i∈N is controllable. In general, loops have different dynamics, because their mirrors might have different levels of cleanliness and thermal losses. However, there are some homogeneous parameters for all loops i∈N , such as the length L of each loop, the cross-sectional area af of the HTF, and the reflective surface S of each loop. As an illustration, Figure 7.1 shows a diagram of the ACUREX plant. 7.2 Concentrated parameter model of a loop This section presents the dynamical model considered to describe the behavior of the ACUREX solar plant [101]. For the sake of clarity, the model variables and parameters are summarized in Table 7.1. In what follows, each loop i∈N will be described using the concentrated parameter model, which provides a global description of a loop by modeling the variation of the internal energy of the HTF. More formally, each loop i∈N can be described by the continuous-time model given by Ci dTi dt =αiηSI −βiSHi(¯ Ti−Ta)−qiPiTi−Tin i,(7.2) 7.3 Characteristics and operational constraints of the system 85 BUFFER STORAGE TANK POWER CONVERSION SYSTEM qTPUMP . . . . . . . . . . . . Loop 1 Loop N . . . Loop j . . . . . . . . . . . . . . . . . . . . . . . . SOLAR COLLECTORS FIELD Figure 7.1 Schematic of the ACUREX parabolic-trough collector solar plant. with Ci=ρf(Ti)cf(Ti)afL, Pi=ρf(Ti)cf(Ti), ¯ Ti=Ti+Tin i 2. (7.3) Here, the global coefficient of thermal losses Hi depends on the outlet ( Ti ), inlet ( Tin i ), and ambient ( Ta ) temperatures, and its expression can be found in [102]. Moreover, η is the overall efficiency of the collectors, considering the optical and geometric efficiencies. In addition, the difference in loop dynamics is modeled through two scaling factors αi and βi that characterize the mirror cleanliness and thermal losses of each loop, respectively. Note that, as (7.2) describes the dynamics of the outlet temperature of loop i , the properties in (7.1) of the HTF in loop ibecome dependent on the corresponding outlet temperature Ti. 7.3 Characteristics and operational constraints of the system Let us assume that the inlet temperature of each loop is equal to the inlet temperature of the fluid to the field, i.e., Tin i=Tin , for all i∈N . Additionally, the inlet temperature, direct normal irradiance (DNI), and ambient temperature are considered disturbances that can be measured or estimated. The outlet temperature of the solar field is given by T=∑N i=1Tiqi qT,(7.4) 86 Chapter 7. Decentralized optimization in constraint-coupled systems Table 7.1 Model’s variables (top) and parameters (bottom). Symbol Description Units kDiscrete-time index - T,Tin Outlet and inlet temperature of the field ºC Ti,Tin iOutlet and inlet temperature of loop iºC TaAmbient temperature ºC ¯ TMean inlet-outlet temperature ºC qiHTF flow in loop il/s qTTotal HTF flow l/s IDirect solar irradiance W/m2 HiCoef. of thermal losses of loop iW/(m2ºC) δDiscretization time s NNumber of loops - LLength of each loop m SReflective surface of each loop m2 afCross-sectional area of the fluid m2 ρfDensity of the HTF kg/m3 cfSpecific heat capacity of the HTF J/(kg ºC) CiThermal capacity of a loop J/ºC ηEfficiency of the collectors - αiCleanliness scale factor of loop iβiLoss scale factor of loop iwhere the total HTF flow of the field is defined as qT= N ∑ i=1 qi.(7.5) On the other hand, several operating restrictions must be considered due to the minimum (Tmin) and maximum (Tmax) temperatures supported by the HTF, and also the operating limits of the pump and valves. Namely, it is required that Tmin ≤Ti≤Tmax,∀i∈N(7.6a) qmin ≤qi≤qmax,∀i∈N(7.6b) qT≤qT,max,(7.6c) where the minimum flow value in a loop, qmin , is such that there is a turbulent regime 7.4 Linear discrete-time model of the system 87 within the loop, and the maximum flow in a loop, qmax , is related to the maximum pressure drop allowed. Moreover, the pump limits the total HTF flow of the field within a maximum value qT,max . Note that although the dynamics of the different loops are independent, they are coupled by the maximum flow constraint (7.6c). Remark 7.3.1 According to the values of the operational constraints for the ACUREX solar field in [101], the inequality Nqmin ≤qT,max ≤Nqmax holds. 7.4 Linear discrete-time model of the system The main objective of the control system is to keep the outlet temperature and the HTF flow of the field close to a given set point despite disturbances. For this purpose, the concentrated parameter model of each loop (7.2) is linearized around a desired operating point (T◦ i,q◦ i)1 . This way, the temperature and flow of each loop i∈N are redefined as the sum of its value at the operating point plus a small increment represented by the deviation variables (xi,ui). That is, Ti=T◦ i+xi,qi=q◦ i+ui.(7.7) Assumption 7.4.1 Since only small variations around the operating point are considered, it is assumed that temperature-dependent parameters, e.g. Pi , Ci and Hi , are fixed at a temperature value equal to T◦ i . Consequently, these time-varying parameters are now assumed to be constant. For the sake of clarity, these parameters are identified with the superscript ◦in the linear model to emphasize their dependence on the operating point. Considering (7.7) and (7.2), we linearize the lumped model of each loop i∈Nas C◦ i dxi dt =−P◦ iq◦ i+βiSH◦ i 2xi−P◦ iT◦ i−Tinui +αiηSI −P◦ iq◦ iT◦ i−Tin−βiSH◦ i 2T◦ i+Tin −2Ta, (7.8) Moreover, we can group and identify Ac i=−P◦ iq◦ i+βiSH◦ i 2, Bc i=−P◦ iT◦ i−Tin, wc i=αiηSI −P◦ iq◦ iT◦ i−Tin−βiSH◦ i 2T◦ i+Tin −2Ta. (7.9) This way, the model in (7.8) can be rewritten as C◦ i dxi dt =Ac ixi+Bc iui+wc i,∀i∈N,(7.10) 1 Apart from considering a coalition approach, we propose the linearization of the lumped parameter model of the solar plant to relieve the computational burden. In this way, the use of a nonlinear distributed parameter models is avoided as it would increase the overall computational burden. 94 Chapter 8. Coalitional MPC of PTC solar fields with population-dynamics assistance Now, note that according to the dynamics in (7.14) , the state of the system at time instant k+nunder a constant input u(k)applied during times k,k+1,...,k+n−1is given by x(k+n) = Anx(k)+Gnu(k)+d(k+n−1),∀(k+n)≥1, where Gn= n−1 ∑ ℓ=0 AℓB d(k+n−1) = n−1 ∑ ℓ=0 Aℓw(k+n−1−ℓ). Consequently, the cost function in (8.4) can be equivalently written as J(u(k)) = Np ∑ n=1x(k)⊤An⊤QAnx(k)+2x(k)⊤An⊤QGnu(k) +2x(k)⊤An⊤Qd(k+n−1) + u(k)⊤G⊤ nQGnu(k) +2u(k)⊤G⊤ nQd(k+n−1) +d(k+n−1)⊤Qd(k+n−1)+Npu(k)⊤Ru(k), and the gradient of J(·)with respect to u(k), which is denoted as g(·), is given by g(u(k)) = Np ∑ n=12G⊤ nQAnx(k)+2G⊤ nQGnu(k)+2G⊤ nQd[k+n−1]+2NpRu(k). Based on these formulations and motivated by the results in [111], to solve the problem in Eq. (8.4) the following iterative dynamics (with state ¯v(κ) = col(¯v1(κ),¯v2(κ),...,¯vN(κ)) and output v(κ) = col(v1(κ),v2(κ),...,vN(κ))) are considered2 ϑj(κ) = maxqmax −qmin −¯vj(κ),0(8.5a) ρi j(κ) = minmaxgi(v(κ))−gj(v(κ)),0,γ(8.5b) µj(κ) = ∑ i∈N¯vi(κ)ϑj(κ)ρi j(κ)−¯vj(κ)ϑi(κ)ρji(κ)(8.5c) ¯vj[κ+1] = ¯vj(κ)+εµj(κ)(8.5d) vj(κ) = ¯vj(κ)+qmin −q◦ i,(8.5e) 2Given a vector a∈Rn, the notation airefers to the i-th element of a. 8.1 Proposed approach 95 for all i,j∈N , where γ,ε∈R>0 are strictly positive constants whose values are defined in Theorem 8.1.1, and the initial condition ¯v(0)satisfies that ¯v(0)∈(¯v∈RN ≥0:∑j∈N¯vj=qT,max −Nqmin ¯vj≤qmax −qmin,∀j∈N).(8.6) Moreover, it is worth to highlight that the updates in (8.5) occur at discrete-time instants κ , which are not necessarily the same time instants k of the discrete-time model in (7.14) . Finally, given the partition P(k) , the resource allocation for each coalition Ci∈P(k) at time κis calculated as qmax Ci(κ) = ∑ j∈Ci vj(κ)+ ∑ i∈Ci q◦ i.(8.7) The iterative dynamics in (8.5) correspond to a discretization of the so-called Smith population dynamics with carrying capacities [112], which are a particular case of the population dynamics studied in [111, 113], and whose evolutionary game theoretical interpretation is as follows. Remark 8.1.2 (Interpretation of the dynamics in (8.5) ) Population dynamics describe the temporal evolution of the strategic interaction of a large population of non-descriptive players under the mean-field approximation [107], i.e., assuming an infinite number of players modeled as continuum of a given mass. Under such a framework, ¯vj(t)∈R≥0 represents the mass of players choosing strategy i∈N at the continuous-time t∈R≥0 (notice that the set of loops N is interpreted as the set of strategies of the game), and the temporal evolution of ¯vj(t)is given by the (mean) evolutionary dynamics model ˙ ¯vj(t) = ∑ i∈N¯vi(t)ϑj(t)ρi j(t)−¯vj(t)ϑi(t)ρji(t). Here, ϑj(t)∈R≥0 denotes the carrying capacity of strategy j at time t (with ϑj(t) = 0 meaning that no more players can choose strategy j at time t ), and ρi j(t)∈R≥0 is the switch rate of players from strategy i to strategy j at time t . Hence, ¯vi(t)ϑj(t)ρi j(t) is the mass of players switching from strategy i to strategy j at time t , while ¯vj(t)ϑi(t)ρji(t) is the mass of players switching from strategy j to strategy i at time t . Adding these terms over all strategies i∈N gives the total variation µj(t)≜˙ ¯vj(t) of the mass of players playing strategy j at time t . Furthermore, to ease the practical application of the dynamics, the above continuous-time evolutionary dynamics model is discretized following a forward Euler discretization with step-size ε∈R>0 , as in (8.5d) . In addition, to relate the population dynamics to the underlying optimization problem, the carrying capacity ϑj(κ) is set according to the HTF capacity of loop i , as in (8.5a) ; and the strategic switch rate ρi j(κ) is defined based on the gradients of the cost function J(·) and the saturation parameter γ∈R>0 , as in (8.5b) . Finally, employing (8.5e) , the variable ¯vj(κ)∈[0,qmax −qmin] is shifted to vj(κ)∈[qmin −q◦ i,qmax −q◦ i] , so that (7.15b) is satisfied. Overall, the dynamics in (8.5) comprise a population-dynamics-inspired iterative algorithm to update the variables vj(κ) , for all i∈N , which are in turn related to the solution of (8.4) . For further details on how to apply population dynamics to solve convex 96 Chapter 8. Coalitional MPC of PTC solar fields with population-dynamics assistance optimization problems, as well as on the related analogies, we refer the interested reader to [108]. As mentioned above, our interest on the dynamics in (8.5) is to exploit certain invariance and asymptotic stability properties that these dynamics have. For the sake of completeness, such properties are formalized next. Theorem 8.1.1 Consider the iterative dynamics in (8.5) with3 γ=2qT,max −Nqminλmax Np ∑ n=1 G⊤ nQGn+NpR!, and 0<ε<γqmax −qmin(N−1)−1 . Moreover, let ¯v(0) satisfy the condition in (8.6) . The following facts hold. 1. For all κ≥0 , the vector v(κ) satisfies the constraints of the problem in Eq. (8.4) , subject to its corresponding constraints. 2. The set of solutions of the problem in Eq. (8.4) , subject to its corresponding constraints, is asymptotically stable under the considered dynamics. 3. For all κ≥0 , if v(κ) is not a solution of the problem in Eq. (8.4) , subject to its corresponding constraints, then J(v[κ+1]) <J(v(κ)). Proof. Fact 1) follows immediately from [111, Theorem 3]. Namely, following the same steps as in the proof of [111, Theorem 3], one can show that under the given conditions for γ , ε , and ¯v(0) , the set in (8.6) is positively invariant under the discrete-time dynamics in (8.5d) . That is, if ¯v(0) satisfies (8.6) , then ¯v(κ) belongs to the set in (8.6) , for all κ≥0 . Consequently, by (8.5e) it holds that qmin −q◦ i≤vj(κ)≤qmax −q◦ i,∀j∈N ∑ j∈N vj(κ) = qT,max −∑ j∈N q◦ i, for all κ≥0 , which means that the vector v(κ) indeed satisfies the constraints of the problem in Eq. (8.4), subject to its corresponding constraints, for all κ≥0. Similarly, Facts 2) and 3) also follow from [111]. More precisely, following the formulation in Remark 8.1.2 and [111, Section 3], it holds that the discrete-time dynamics in (8.5a) - (8.5d) can be rewritten compactly as ¯v(κ+1) = ¯v(κ) +εL(κ)f(κ) , where f(κ) = −g(v(κ)) and L(κ)∈RN×N is a symmetric positive semi-definite matrix, for all κ≥0 [111, Lemma 1]. Also, based on (8.5e) , the output vector v(κ) is a constant shift of the state vector ¯v(κ) , for all κ≥0 . Hence, consider the (convex quadratic) Lyapunov function candidate given by V(v(κ)) = J(v(κ))−J(v∗) , where v∗ is an arbitrary optimal 3Given a square matrix A∈Rn×n, the notation λmax (A)refers to the maximum eigenvalue of A. 8.2 Control scheme 97 solution of the problem in (8.4) . Following the same steps as in the proof of [111, Theorem 5] it can be shown that V(v(κ+1))−V(v(κ)) ≤ −ε(1−εθ)f(κ)⊤L(κ)f(κ), where θ=γqmax −qmin(N−1) . Moreover, from [111, Lemma 1] and [111, Theorem 1], one can conclude that f(κ)⊤L(κ)f(κ) = 0 if and only if v(κ) = v∗ , and that v(κ+1) = v(κ) if and only if v(κ) = v∗ . Thus, if 0<ε<1/θ , then V(v(κ+1))−V(v(κ)) <0 , for every v(κ+1)=v(κ)=v∗ . Consequently, v∗ is asymptotically stable under the dynamics in (8.5) and Fact 2) holds. Also, Fact 3) holds from the previous Lyapunov stability analysis by noting that V(v(κ+1))−V(v(κ)) = J(v(κ+1))−J(v(κ)).■ Based on Theorem 8.1.1, note that if the constraint in (8.6) is satisfied, then the feasible set of the problem in Eq. (8.4) , subject to its corresponding constraints, is positively invariant under the dynamics in (8.5) , which in turn implies that the resource allocation in (8.7) always satisfies ∑Ci∈P(k)qmax Ci=qT,max , as desired. In addition, such a property holds without relying on recursive projections onto the feasible set, which simplifies the execution of the proposed dynamics. Furthermore, the cost function in (8.4) strictly decreases under the trajectories of the dynamics in (8.5) , and its minimizer is achieved asymptotically. Thus, under the dynamics in (8.5) , the resource allocation v(κ2) at time κ2 is always better, in the sense of the problem in (8.4) , than the allocation v(κ1) at time κ1 , for all 0≤κ1<κ2 and granted that the allocation at time κ1 is not optimal (in this case v(κ2) = v(κ1) ). These facts show that, regardless of the number of update iterations applied to the dynamics in (8.5) , the corresponding resource allocation is always better than the initial one (again, here better is to be understood in the sense of the problem in Eq. (8.4) ). This property is attractive in the scenario where the updates in (8.5) can only be executed a limited number of times due to real-time constraints. Moreover, note that this property does not hold under other iterative approaches, such as the so-called primal-dual gradient dynamics [114], where the feasibility of the solution is guaranteed only in the asymptotic sense. 8.2 Control scheme The proposed coalitional model predictive control approach consists of two layers. At the top layer, the population-dynamics-assisted resource allocation method of Section 8.1.2 is employed to solve the unitary control horizon problem of Eq. (8.4) , and to allocate the resource qmax Ci to each coalition Ci∈P(k) . At the bottom layer, on the other hand, the controller of each coalition Ci∈P(k) solves the local model predictive control problem of Eq. (8.3) subject to the resource allocation provided by the top layer. Both layers can be executed at different time-scales, so that the bottom layer uses the most up-to-date value of qmax Ciwhen required. The proposed approach is summarized in Algorithm 4. Remark 8.2.1 The parameters Niter , ttop , and tbot in Algorithm 4 are set according to the computational capacities of the corresponding controllers, as well as to the associated time-scales of the solar plant. In particular, the time-scale tbot is set as a multiple of the 98 Chapter 8. Coalitional MPC of PTC solar fields with population-dynamics assistance discretization time δ of the model, while the time-scale ttop is set as a multiple of tbot . The rationale for setting ttop slower than tbot is to avoid recomputing the resource allocation and partitioning of the system too frequently, and thus to reduce the computational burden of the overall approach. In fact, to further reduce the computational burden, the upper control layer might be executed in an event-triggered fashion based on the DNI disturbances. Algorithm 4: Overall control scheme. The top and bottom control layers are repeatedly executed at possibly different time-scales. Top control layer (executed at a time-scale ttop) 1. Given an initial condition ¯v(0) satisfying (8.6) , the dynamics in (8.5) are updated a number of Niter iterations. The resulting vector v(Niter)provides a feasible approximation for the solution of the problem in Eq. (8.4) (with the error of the approximation converging to 0as Niter →∞) subject to its corresponding constraints. 2. The field partition P(k)is computed by forming coalitions using the selected criterion4. 3. The resource allocation for each coalition Ci∈P(k)is computed by means of (8.7). Namely, qmax Ci(Niter) = ∑j∈Civj(Niter)+∑i∈Ciq◦ i,∀Ci∈P(k). Bottom control layer (executed at a time-scale tbot) 1. For every coalition Ci∈P(k), the local optimization problem in (8.3) is solved. Such a local optimization problem is subject to the (most up-to-date) resource allocation provided by the top layer, and its solution corresponds to the profile u∗ Ci(k), where kdenotes the current time instant under the time-scale tbot. 2. Every coalition Ci∈P(k)applies the first component of the optimal sequence u∗ Ci(k)to its associated loops (the deviation variables are transformed back to the original inputs by adding their corresponding operating values). 8.3 Simulation results This section provides simulation results on a 10- and 100- loop parabolic-trough solar plant such as the one described in Chapter 7.2. The control objective is to operate the plant around a desired operating temperature of 250 ◦ C, while satisfying the coupled constraint in the total HTF flow, as well as the remaining operational constraints. Hereafter, the results obtained with the proposed coalitional approach, summarized in Algorithm 4, are 4 The partitioning of the field can be done following different criteria such as geographical proximity, e.g., grouping loops that are close to each other in the solar plant; difference in the DNI, e.g., associating loops that receive more solar irradiance with those that are dirty or shaded; or simply by creating random coalitions. In particular, this chapter does not focus on a specific partitioning method. Instead, it is assumed that the partition P(k)is given by an arbitrary previously selected criterion. 8.3 Simulation results 99 compared to those obtained with a centralized MPC. For this purpose, the index Jperf is introduced to measure the performance of a given control algorithm throughout the entire simulation. The selected performance index is formulated as Jperf = Tsim ∑ k=0 ∥x(k+1)∥2 Q+∥u(k)∥2 R,(8.8) where the simulation length is set to Tsim =60 min for both case studies. Moreover, a one-hour DNI profile is used to model the disturbances (the considered DNI profile is generated as in [39]), while the ambient temperature Ta and the inlet temperature Tin of the HTF to the field are considered constant for simplicity. As shown in the DNI profiles of Figure 8.1 and Figure 8.5, if the total length of a loop is discretized into one-meter segments, the passive parts of the loops correspond to the segments {37−42} , {79−96} , and {133−138} , respectively. Since these segments are not reached by concentrated solar radiation, their DNI is set to 0. For these numerical simulations, the partition P(k) of the field is calculated by grouping sets of |Ci| unbalanced loops according to the effective solar irradiance they receive. For instance, loops that are dirty or shaded due to a cloud are associated with loops that receive more DNI, either because they are not affected by any cloud or because they are cleaner. Moreover, let us assume that all coalitions have the same number of loops, i.e., |Ci|=|Cℓ| , for all i,ℓ ∈ {1,2,...,M}with i=ℓ. On the other hand, to compare the coalitional and centralized control strategies by means of the required computational time, τcomp is defined as the average computation time in the simulation, calculated as τcomp =∑Tsim k=1τ(k) Tsim ,with τ(k) = ∑Ci∈P(k)τCi(k) |P(k)|,(8.9) where τ(k) is the average computation time per coalition at the time instant k , and τCi(k) is the time it takes for the coalition Ci to solve the local optimization problem in Eq. 8.3 subject to its corresponding constraints. Seeking to be more rigorous with the results provided, Table 8.1 shows, in addition to the average computation time τcomp , the value of the standard deviation σcomp of the sample. The times included in Table 8.1 are measured in seconds. Finally, Table 8.2 includes the values of the homogeneous parameters for all loops, and the values of operational constraints, which are taken equal for both simulations. Table 8.1 Time comparison for the centralized and coalitional simulations. Strategy τcomp σcomp Units Cen. MPC N=10 0.2804 0.0466 s Coal. MPC N=10 0.1463 0.0154 s Cen. MPC N=100 2.7260 0.4522 s Coal. MPC N=100 0.2511 0.0045 s 100 Chapter 8. Coalitional MPC of PTC solar fields with population-dynamics assistance Table 8.2 Parameters (left) and constraints (right) of ACUREX field. Symbol Value Units Symbol Value Units L174 m Tmin 220 ºC S267.4 m2Tmax 300 ºC af7.55 ·10-4 m2qmin 0.2 l/s η0.64 - qmax 1.5 l/s αi[0.6, 1] - βi[1, 1.25] - 8.3.1 Results for the 10-loop PTC solar field The first case of study corresponds to the 10-loop ACUREX field (see Chapter 7.1). This scenario considers qT,max =5.7209 l/s as the maximum flow value, which corresponds to the total HTF flow at the desired operating point of the system. Moreover, without loss of generality, five coalitions of two loops are considered. Two moving clouds have been generated for the DNI profile, which cross the field at specific times. As shown in Figure 8.1, a small cloud moves through the plant between t=13 min and t=26 min, and a larger one appears between t=38 min and t=52 min. These clouds produce a decrease in the solar radiation received by the shaded loops, which in turn has an impact on the partition P(k) of the field. As an illustration, the coalitional topology of the entire simulation is depicted in Figure 8.2, where it can be seen that if the plant is not affected by the presence of any cloud, then the coalitions remain constant, but when a cloud enters the field, the partition P(k) of the system changes, resulting in new coalitions. Finally, it is worth to highlight that besides the two moving clouds, it has been considered that some loops’ collectors might be cleaner than others, and thus receive more solar radiation. This fact adds more heterogeneity to the components of the field. 20 40 60 80 100 120 140 160 1 2 3 4 5 6 7 8 9 10 100 200 300 400 500 600 700 800 900 1000 (a) Cloud passing at t=20min. 20 40 60 80 100 120 140 160 1 2 3 4 5 6 7 8 9 10 100 200 300 400 500 600 700 800 900 1000 (b) Cloud passing at t=45min . Figure 8.1 Effective DNI profile of the 10-loop collector field with clouds passing. 8.3 Simulation results 101 Figure 8.2 Coalition forming of the 10-loop collector field. Figure 8.3 Outlet temperature (top plot) and HTF flow (bottom plot) for each loop i∈ {1,2,...,10}in the 10-loop field, when using the coalitional approach. In Figure 8.3, the top plot shows the evolution over time of the outlet temperature of each loop, while the bottom plot depicts the evolution of the HTF flow circulating through each loop. Here, loops that have cleaner collectors and are therefore more efficient, have a higher flow value at the operating point q◦ i (see, for example, that flows in loops # 1 or # 6 converge to a value qi far above the corresponding values in loops # 3 or # 5 ). In addition, the passing of the larger cloud has a greater impact on both the temperature and the flows 102 Chapter 8. Coalitional MPC of PTC solar fields with population-dynamics assistance of the system, when compared to the passing of the smaller cloud (as expected). Moreover, such an impact is more notable in the less efficient loops (see, for example, the drop in the outlet temperature of loop # 3 , represented by a pink line in Figure 8.3, when the second cloud appears). A comparison of the proposed coalitional approach with respect to the centralized approach is presented in Figure 8.4 in terms of the outlet temperature and the total HTF flow of the field. We highlight that the temporal evolution of both T and qT under the coalitional approach follows up closely the ones of the centralized solution. In fact, the proposed coalitional algorithm results only in a 1.1771% drop in the overall performance Jperf with respect to the centralized solution. Likewise, as shown in Table 8.1, the computational time of the coalitional approach is lower than that of the centralized strategy, reducing the centralized time by half. This result is expected because the centralized approach solves the problem in Eq. (7.16) using a single controller, whereas the coalitional one employs N/|Ci|decentralized controllers to solve the local problems of Eq. (8.3) in parallel. 5 10 15 20 25 30 35 40 45 50 55 60 248 249 250 251 5 10 15 20 25 30 35 40 45 50 55 60 3 4 5 6 Figure 8.4 Outlet temperature (top plot) and total HTF flow (bottom plot) of the 10- loop field. Solid lines represent the evolution using the proposed coalitional algorithm, while dashed lines represent the result using a centralized MPC. Here, the coupled constraint qT,max is represented with a dotted line, and T and qTare given by (7.4) and (7.5), respectively. 8.3 Simulation results 103 8.3.2 Results for the 100-loop PTC solar field Let us now consider a 100-loop extension of the original 10-loop ACUREX field. As such, the values of the local constraints in (7.6) are also set according to Table 8.2, yet, since the number of loops is ten times greater, the value of qT,max is set an order of magnitude higher than the 10-loop plant. Likewise, the maximum flow value in the plant is considered equal to qT,max =53.8834 l/s, which corresponds to the total HTF flow at the desired operating point of the system. Besides, without loss of generality, in this scenario ten coalitions of ten loops are considered. For this simulation, the DNI profile generates two moving clouds that cross over the 100-loop field. A large cloud moves through the plant almost during the whole simulation, between t=4 min and t=52 min, while a smaller cloud appears between t=16 min and t=46 min. As an illustration, Figure 8.5 shows the two clouds over the plant halfway through the simulation, at t=30 min. In addition, as before, it is considered that some loop collectors may be cleaner than others and therefore receive more solar radiation. Figure 8.6 shows the evolution over time of the outlet temperature Ti and the HTF flow qi of ten random loops of the 100-loop plant. The small cloud crosses diagonally, affecting approximately the upper half of the field, while the large cloud crosses the entire plant, creating shadows on all of the loops. The consequences of passing clouds can be seen in Figure 8.6. As in the case of the 10-loop plant, the most efficient loops stabilize in a higher flow value. See, for example, in the upper part of the field, loop # 93 (light-purple line) reaches a higher value of the flow than loop # 82 (red line). Similarly, when clouds shade Figure 8.5 Effective DNI profile of the 100-loop collector field with clouds passing at t=30min. 110 Chapter 9. Resource allocation based on coalitional MPC for PTC solar fields Algorithm 5: Control Scheme Initialization. At k=0, we distribute qT,max as in (9.3). At each time instant k: Top layer: 1. The field partition P(k)is computed by grouping unbalanced loops regarding the effective solar irradiance they receive. To this end, we define a vector Ikthat sorts the effective solar irradiance received by loops i∈Nat time instant kin a descending order. Coalitions are made by grouping the loops at the beginning of the vector Ikwith those at the end. 2. If D(k)=/0 and S(k)=/0, the excess sT(9.8) is distributed among coalitions Ci∈D(k)following qmax,+ Ci=qmax Ci+∆Ci, where ∆Ci=|Ci|sT |ND(k)|,∀Ci∈D(k),(9.11) with |ND(k)| being the total number of loops in the subset D(k) . On the other hand, the maximum flow constraint for Ci∈S(k)is updated as (9.5). 3. If D(k) = /0, then S(k) = P(k). In this case, the total excess sT(9.8) is proportionally distributed among all coalitions using ∆Ci(k) = |Ci|sT |N|,∀Ci∈P(k).(9.12) The maximum flow constraint is updated as qmax,+ Ci=qmax Ci−sCi+∆Cifor all coalitions Ci∈P(k). Bottom layer: 4. For every coalition Ci∈P(k), the local optimization problem (9.9) is solved subject to the most recent maximum flow constraint provided by the top layer as a consequence of the flow-sharing between coalitions. 5. According to (9.6) and the optimal sequence uCi , loops i∈Ci will apply to the first component of their input vectors ui, for all Ci∈P(k). 6. For every coalition Ci∈P(k)the first component of sCiwill determine: a) If sCi(0) = 0, the coalition will send a flow request: Ci∈D(k). b) If sCi(0)>0, the coalition has an excess flow make available: Ci∈S(k). 9.3 Simulation results In this section, we apply the proposed coalitional scheme on a parabolic-trough solar collector field such as the one described in Chapter 7, with N=100 loops. The main objective is to control the plant around an operating temperature of 250◦C , meeting the 9.3 Simulation results 111 constraints described in the previous sections despite external disturbances. The 100-loop benchmark has been simulated using an hour and a half DNI profile to model the solar radiation, setting the radiation to zero in the passive segments of the loops, corresponding to the sections {37 −42} , {79 −96} , and {133 −138} , of the L=174 m loop. For simplicity, we consider the ambient temperature Taof the field and the inlet temperature Tin of the HTF as constant parameters along the simulation length Tsim =90min. The results obtained using the coalitional approach of Algorithm 5 will be compared with those obtained with I. Centralized MPC, which solves Problem (7.16). II. Decentralized MPC, where each loop i∈N is governed by a local controller with no information from the rest that solves min ui Np−1 ∑ n=0 ∥xi(n+1)∥2 Qi+∥ui(n)∥2 Ri, s.t. (7.13),(7.15a),(7.15b), ui(k)≤qmax Ci(0)−q◦ i, ∀n=0,...,Np−1. (9.13) where Qi and Ri are set according to the global matrices in (7.16) , and the vector ui is defined as the sequence of local inputs ui= [ui(n)]Np−1 n=0. Remark 9.3.1 Problem (9.13) is solved without the resource-sharing mechanism between loops proposed in Algorithm 5. Therefore, the maximum flow constraint is simply set according to (9.3) with |Ci|=1for all Ci∈P(k). Let us define the index Jperf = Tsim ∑ k=1 ∑ i∈N ∥xi(k+1)∥2 Qi+∥ui(k)∥2 Ri, which will be used to evaluate the performance of a given control strategy during the simulation length Tsim . In addition, we provide an average value for the computation time τcomp, which has been calculated as τcomp =∑Tsim k=1τ(k) Tsim ,with τ(k) = ∑Ci∈PτCi(k) |P(k)|,(9.14) where τ(k) corresponds to the average computation time per coalition at the time instant k , and τCi(k) represents the time it takes a given coalition to solve its optimization problem. For the centralized MPC, we work with the grand coalition Ci=N , i.e., |P(k)|=1 , and for the decentralized MPC we study |P(k)|=N singletons. In this regard, Table 9.1 shows the value of τcomp for the different simulations. Finally, we consider the parameters S=267.4m2 , af=7.55·10−4m2 , η=0.64 . The values of αi∈[0.6,1] and βi∈[1,1.25] have been generated following a uniform random distribution. Also, the operational 112 Chapter 9. Resource allocation based on coalitional MPC for PTC solar fields constraints have been set as Tmin =220◦C , Tmax =300◦C , qmin =0.2l/s and qmin =1.5l/s . Note that these values correspond to those described in Table 8.2 of Chapter 8. Table 9.1 Average computational time using different control strategies. Loops Strategy ¯ τ(s) N=10 Cen. MPC 0.3769 Coal. MPC |Ci|=20.1936 Dec. MPC 0.0359 N=100 Cen. MPC 3.5998 Coal. MPC |Ci|=20.1998 Coal. MPC |Ci|=10 0.4553 Dec. MPC 0.0363 9.3.1 Results for the 10-loop PTC solar field For this scenario, we consider coalitions of two-loops, i.e., |Ci|=2 , and a DNI profile with two moving clouds: one small that moves through the field during t={13−26} min, and a larger one that appears between t={38 −52} min. The top layer of Algorithm 5 calculates the partition P(k) by clustering loops that receive more DNI with other loops that are less efficient or shaded by a cloud, according to Step 1 of Algorithm 5. A shaded loop will receive lower DNI, becoming less efficient, so that it can provide its remaining flow to other loops in the same coalition to increase their performance. In this sense, Figure 9.1 shows the evolution of the field partition P(k) during the whole simulation, where we can see the formation of new coalitions when the plant is affected by passing clouds. Loops represented with the same color at a given time instant belong to the same coalition. Furthermore, the same colors have been used in Figure 9.2 to depict the evolution of the maximum flow qmax Ci of the different coalitions. As can be seen, coalitions start with the same maximum value at t=0 (recall (9.3) ; in this case, |Ci|=|Cℓ| for all i=ℓ ), represented as the dotted line in Figure 9.2. According to the dynamic sharing of the flow Figure 9.1 Coalition forming of the 10-loop collector field. 9.3 Simulation results 113 constraint described in Algorithm 5, coalitions change their maximum constraint value, guaranteeing the satisfaction of the common constraint (7.6c). Figure 9.3, represents the evolution over time of the outlet temperature of each loop 10 20 30 40 50 60 70 80 90 0.9 1 1.1 1.2 1.3 Figure 9.2 Evolution of the maximum flow constraint (9.5) for all coalitions. Figure 9.3 Outlet temperature (top plot) and HTF flow (bottom plot) for each loop i∈ {1,2,...,10}in the 10-loop field, when using the coalitional approach. 114 Chapter 9. Resource allocation based on coalitional MPC for PTC solar fields Ti , and the corresponding HTF flow qi circulating through each of them. As mentioned in Chapter 7.1, we consider loops to have different efficiency, which means that some of them have cleaner collectors and less thermal losses than others. Compared to those loops that have dirtier collectors and more thermal losses, efficient loops should stabilize at a higher flow value at the operating point q◦ i . By analyzing Figure 9.3, we can distinguish that loop # 5 (pink line) is less efficient than loop # 10 (green line). For example, see that the operating flow value of loop # 5 ( q◦ 5 ) is 0.4l/s , while q◦ 10 ≈0.7l/s . Also, the impact of cloud disturbances will be greater in less efficient loops, since, for example, loop # 3 (yelow line) suffers a temperature drop of approximately 15◦C when the large cloud shades it. Figure 9.4 compares the evolution of the outlet temperature and the total flow of the field when using the proposed coalition scheme (blue solid lines) with the results obtained with the centralized approach (orange dashed lines) and the decentralized MPC (yellow dashed-dotted lines). The value of the coupled constraint qT,max is represented in a black dashed-dotted line. The presence of clouds directly affects the solar radiation received by the collector, which decreases significantly, reducing the field outlet temperature and HTF flow. In the coalitional scheme, loops transfer and absorb flow according to the needs of the system, shrinking and enlarging their constraint set qmax Ci . As a result, regulation towards the operating point after the passing clouds is slower than with the centralized solution. However, the performance loss in index Jperf with respect to the centralized solution is 12.3% with the coalition approach, reducing the computation time by nearly half, as shown in Table 9.1. 10 20 30 40 50 60 70 80 90 250 252 254 256 10 20 30 40 50 60 70 80 90 3 4 5 6 Figure 9.4 Outlet temperature and total HTF flow of the 10-loop field. 9.3 Simulation results 115 Loops solve independent optimization problems using the local MPC (see Problem (9.13) ) with no information about other loops. For this reason, the maximum flow constraint of a given loop cannot change during the simulation. This entails that loops that need more flow to regulate toward their operating point cannot benefit from the flow excess of other loops, resulting in a huge performance loss regarding the centralized approach of 4.6·103% . In turn, those loops that have an extra flow cannot share it with other loops. As a consequence, when the local controllers are used, we do not manage to regulate the outlet temperature of the solar field, as depicted in Figure 9.4. 9.3.2 Results for the 100-loop PTC solar field Since this scenario considers a larger plant, we have simulated the proposed coalitional algorithm for two different cases: using two (|Ci|=2) and ten-loop (|Ci|=10) coalitions in the partition of the system. As in the previous case, we study the behavior of the field when two moving clouds shade the collectors. In this case, a large cloud moves through the plant during t={4−52} min, while a smaller cloud crosses the field diagonally between t={16−46} min. The top layer of Algorithm 5 calculates the partition P(k) by clustering loops that receive more DNI with other loops that are less efficient or shaded by a cloud. In this regard, Figure 9.5 shows the evolution of the field partition P(k) during the whole simulation when using |Ci|=10 , where we can appreciate the formation of new coalitions when the plant is affected by passing clouds. Figure 9.6 shows the evolution of the outlet temperature and HTF flow of 10 random loops of the 100-loop field. As in the case of the 10-loop field, the scale factors αi and βi will model the cleanliness and thermal losses of each loop, respectively. A small value of αi represents loops with dirtier collectors and therefore will receive less effective solar irradiance. In turn, a high value of βi represents loops with greater thermal losses. Both 1 10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90 100 1 2 3 4 5 6 7 8 9 10 Figure 9.5 Coalition forming of the 100-loop collector field, when using |Ci|=10. 116 Chapter 9. Resource allocation based on coalitional MPC for PTC solar fields 0 10 20 30 40 50 60 70 80 90 235 240 245 250 10 20 30 40 50 60 70 80 90 0.2 0.4 0.6 0.8 Figure 9.6 Outlet temperature (top plot) and HTF flow (bottom plot) for 10 random loops in the 100-loop field, when using the |Ci|=10 coalitional approach. of these characteristics make loops less efficient, stabilizing at low flow value q◦ i . See, for example, the differences in the evolution of loop # 5 (red line) and loop # 35 in Figure 9.6. Loop # 35 stabilizes at a high flow value at its operating point, and the passing clouds practically do not affect its outlet temperature. Therefore, loop # 35 is more efficient than # 5 , since the effect of the clouds in the latter is much more notable. Loops tend to decrease their flow value qi when a cloud shades them. In this regard, efficient loops that stabilize at a higher HTF flow value will have more margin to decrease their local flow when a cloud affects them, maintaining the outlet temperature in its operating value. Recall, for example, loop # 93 (dark-blue line), which manages to maintain its outlet temperature around 250◦C by decreasing its flow value. On the other hand, less efficient loops, such as loop # 82 , drop the value of its flow but saturates at its lower bound qmin. Figure 9.7 compares the evolution of the outlet temperature and the total flow of the field when using different control strategies. Blue solid lines represent the result using the coalition algorithm proposed with |Ci|=2 , while the green dotted lines represent the evolution using |Ci|=10 . The orange dashed lines represent the centralized MPC result, and the yellow dashed-dotted lines show the evolution when a decentralized MPC is applied. In the bottom plot, the black dashed-dotted line represents the maximum flow allowed in the field, which is equal to the operating point. The presence of clouds directly affects the solar radiation received by the collector, which decreases significantly, reducing the field outlet temperature and HTF flow. 9.3 Simulation results 117 10 20 30 40 50 60 70 80 90 250 252 254 256 10 20 30 40 50 60 70 80 90 45 50 55 10 20 30 40 50 60 70 80 90 249.8 250 250.2 Figure 9.7 Outlet temperature and total HTF flow of the 100-loop field using different control strategies. For the coalitional approach, different coalition sizes |Ci| have been considered. In the coalitional scheme, loops transfer and absorb flow according to the needs of the system, shrinking and enlarging their constraint set qmax Ci . As can be seen in Figure 9.7, the coalitional scheme takes longer to recover from the cloud event and regulate the value of the temperature, since conflicts with the shared resource appear when loops work close to the operating point. However, Algorithm 5 using |Ci|=10 and |Ci|=2 results in a decrease in the overall performance Jperf of 0.82% and 3.6% , respectively, with respect to the centralized solution. Considering a smaller size of coalitions implies that we have more agents working in parallel, but at the same time, these local controllers handle less amount of information. As shown in Table 9.1, we managed to significantly reduce computational time compared to the centralized strategy. The centralized approach solves Eq. (7.16) using a single controller, while in the coalitional approach, |P(k)| controllers work in parallel. Therefore, considering larger coalitions results in a longer computational time, but also in better performance than when considering smaller coalitions. Finally, applying a local MPC to each loop of the N=100 plant is not feasible even if the computing time is minimal for each loop, as shown in Table 9.1. In this case, loops solve independent optimization problems (see Eq. (9.13) ) without information about other loops. For this reason, the maximum flow constraint of a given loop cannot change during the simulation. We are not allowing the loops to dynamically distribute their excess flow, nor are we allowing them to demand more, so there are loops that need to increase their 118 Chapter 9. Resource allocation based on coalitional MPC for PTC solar fields flow and cannot. This, besides increasing the performance index Jperf in a 9.1·103% , may cause the violation of the operational constraints, in addition to not complying with the control objective. 9.4 Conclusions We have presented a two-layer coalitional model predictive control approach for the regulation of a parabolic-trough solar collector field that allows us to distribute the large-scale optimization problem. We have dealt with constraint-coupled subsystems, and distributed the coupled constraint using a resource-sharing mechanism to solve the optimization in a decentralized manner while guaranteeing the satisfaction of operational constraints. Furthermore, simulation results demonstrate that we achieve a small performance loss with respect to the centralized solution, as well as reducing the computational load, as can be seen in Table 9.1 through the average processing times. Part IV Conclusions and future research directions 119 126 List of Figures 4.5 Performance loss and the communication savings of the coalitional approach in comparison with the centralized control for different values of ˆ α42 4.6 Evolution of the public variables upu ji and scale factors αji controlled by each agent iwhen using ˆ α=0.015 42 4.7 Representation of the state constraint set Xi (blue), mRPI set Ri (gray) and (on the bottom plot) nominal constraint set ¯ Xi (light-blue) for every subsystem. The mRPI has been calculated using the methods proposed in [66] and [72], both leading to the same set 43 4.8 Evolution of the flow rates of the four pumps for the system with reduced coupling terms Bi j 44 5.1 Diagram of the eight tanks plant, where the coupling relations are represented by the colored pipes 56 5.2 State trajectories of the different control approaches for tracking along with the tube section evolution for all agents i∈Nfor the coalitional scheme 59 5.3 Evolution of the nominal public inputs ¯upu i j and its scale factors αi j for all agents 59 5.4 Communication topology depicted by the state of the communication links 61 6.1 Diagram of the input coupling between the three agents of the system under study, considering the private and public variable decomposition 73 6.2 Theoretical (dashed-lines) and actual (solid-lines) planes for system E1 and k∗=20 74 6.3 Theoretical ¯ F20 and actual F20 regions enclosed by the theoretical and actual planes of Figure 6.2 for system E1and k∗=20 75 6.4 Comparison of the projections of the theoretical ¯ Fk∗ and the actual Fk∗ regions for systems E1to E3for different values of k∗76 6.5 Representation of the k -step constraint admissible sets Cαi k(Xi) , for the three subsystems of E1with k<10. 78 6.6 Representation of the maximal RCI set for the third subsystem, Cαi ∞(X3) , for different values of {α23,α31}79 7.1 Schematic of the ACUREX parabolic-trough collector solar plant 85 8.1 Effective DNI profile of the 10-loop collector field with clouds passing 100 8.2 Coalition forming of the 10-loop collector field 101 8.3 Outlet temperature (top plot) and HTF flow (bottom plot) for each loop i∈ {1,2,...,10}in the 10-loop field, when using the coalitional approach 101 8.4 Outlet temperature (top plot) and total HTF flow (bottom plot) of the 10-loop field. Solid lines represent the evolution using the proposed coalitional algorithm, while dashed lines represent the result using a centralized MPC. Here, the coupled constraint qT,max is represented with a dotted line, and T and qT are given by (7.4) and (7.5), respectively 102 8.5 Effective DNI profile of the 100-loop collector field with clouds passing at t=30min 103 8.6 Outlet temperature (top plot) and HTF flow (bottom plot) for 10 random loops in the 100-loop field, when using the coalitional approach 104 List of Figures 127 8.7 Outlet temperature (top plot) and total HTF flow (bottom plot) of the 100-loop field. Solid lines represent the evolution using the proposed coalitional algorithm, while dashed lines represent the result using a centralized MPC. Here, the coupled constraint qT,max is represented with a dotted line, and T and qT are given by (7.4) and (7.5), respectively 105 9.1 Coalition forming of the 10-loop collector field 112 9.2 Evolution of the maximum flow constraint (9.5) for all coalitions 113 9.3 Outlet temperature (top plot) and HTF flow (bottom plot) for each loop i∈ {1,2,...,10}in the 10-loop field, when using the coalitional approach 113 9.4 Outlet temperature and total HTF flow of the 10-loop field 114 9.5 Coalition forming of the 100-loop collector field, when using |Ci|=10 115 9.6 Outlet temperature (top plot) and HTF flow (bottom plot) for 10 random loops in the 100-loop field, when using the |Ci|=10 coalitional approach 116 9.7 Outlet temperature and total HTF flow of the 100-loop field using different control strategies. For the coalitional approach, different coalition sizes |Ci| have been considered 117 List of Tables 3.1 Comparison of the performance costs for different modes of operation 25 4.1 Comparison of the overall performance for different control approaches 41 5.1 Comparative of the different MPC approaches for tracking 58 7.1 Model’s variables (top) and parameters (bottom) 86 8.1 Time comparison for the centralized and coalitional simulations 99 8.2 Parameters (left) and constraints (right) of ACUREX field 100 9.1 Average computational time using different control strategies 112 129 Bibliography [1] E. Camacho, “OCONTSOLAR project – optimal control of thermal solar energy systems,” 2018. [Online]. Available: https://www.europeandissemination.eu/ ocontsolar-project/2563 [2] R. Scattolini, “Architectures for distributed and hierarchical model predictive control – A review,” Journal of Process Control, vol. 19, no. 5, pp. 723–731, 2009. [3] F. Katiraei and R. Iravani, “Power management strategies for a microgrid with multiple distributed generation units,” IEEE Transactions on Power Systems, vol. 21, pp. 1821 – 1831, 12 2006. [4] M. Yazdanian and A. Mehrizi-Sani, “Distributed control techniques in microgrids,” IEEE Transactions on Smart Grid, vol. 5, no. 6, pp. 2901–2909, 2014. [5] G. Chaloulos, P. Hokayem, and J. Lygeros, “Distributed hierarchical MPC for conflict resolution in air traffic control,” in Proceedings of the 2010 American Control Conference, June 2010, pp. 3945–3950. [6] J. Hernández, S. Ossowski, and A. Garca-Serrano, “Multiagent architectures for intelligent traffic management systems,” Transportation Research Part C: Emerging Technologies, vol. 10, 10 2002. [7] E. Bonabeau, M. Dorigo, and G. Theraulaz, Swarm intelligence: From natural to artificial systems. Oxford University Press, 1999. [8] E. Şahin, “Swarm robotics: From sources of inspiration to domains of application,” in Swarm Robotics, E. Şahin and W. M. Spears, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2005, pp. 10–20. [9] M. Kordestani, A. A. Safavi, and M. Saif, “Recent survey of large-scale systems: Architectures, controller strategies, and industrial applications,” IEEE Systems Journal, vol. 15, no. 4, pp. 5440–5453, 2021. 131 132 Bibliography [10] E. F. Camacho and C. Bordons, Model Predictive Control, 2nd ed., ser. Advanced Textbooks in Control and Signal Processing. London: Springer, 2013. [11] J. Rawlings, “Tutorial overview of model predictive control,” IEEE Control Systems Magazine, vol. 20, no. 3, pp. 38–52, 2000. [12] D. Mayne, J. Rawlings, C. Rao, and P. Scokaert, “Constrained model predictive control: Stability and optimality,” Automatica, vol. 36, no. 6, pp. 789–814, 2000. [13] J. Richalet, A. Rault, J. L. Testud, and J. Papon, “Model predictive heuristic control: Applications to industrial processes,” Automatica, vol. 14, no. 2, pp. 413–428, 1978. [14] S. J. Qin and T. A. Badgwell, “A survey of industrial model predictive control technology,” Control Engineering Practice, vol. 11, no. 7, pp. 733–764, 2003. [15] P. Neumann, “Communication in industrial automation–What is going on?” Control Engineering Practice, vol. 15, no. 11, pp. 1332–1347, 2007. [16] J. M. Guerrero, M. Chandorkar, T.-L. Lee, and P. C. Loh, “Advanced control architectures for intelligent microgrids—part i: Decentralized and hierarchical control,” IEEE Transactions on Industrial Electronics, vol. 60, no. 4, pp. 1254– 1262, 2013. [17] D. Ding, Q.-L. Han, Z. Wang, and X. Ge, “A survey on model-based distributed control and filtering for industrial cyber-physical systems,” IEEE Transactions on Industrial Informatics, vol. 15, no. 5, pp. 2483–2499, 2019. [18] X. Ge, F. Yang, and Q.-L. Han, “Distributed networked control systems: A brief overview,” Information Sciences, vol. 380, pp. 117–131, 2017. [19] A. N. Venkat, I. A. Hiskens, J. B. Rawlings, and S. J. Wright, “Distributed MPC strategies with application to power system automatic generation control,” IEEE Transactions on Control Systems Technology, vol. 16, no. 6, pp. 1192–1206, 2008. [20] R. R. Negenborn and J. Maestre, “Distributed model predictive control: An overview and roadmap of future research opportunities,” IEEE Control Systems Magazine, vol. 34, no. 4, pp. 87–97, 2014. [21] J. M. Maestre and R. R. Negenborn, Distributed Model Predictive Control Made Easy, ser. Intelligent Systems, Control and Automation: Science and Engineering. Springer, 2014, no. 69. [22] J. Lunze, Feedback Control of Large Scale Systems. New York: Prentice Hall PTR, 1992. [23] D. Siljak, Decentralized Control of Complex Systems, ser. Mathematics in Science and Engineering. Academic Press, 1991. [24] H. Wu, “Decentralized adaptive robust control for a class of large-scale systems including delayed state perturbations in the interconnections,” IEEE Transactions on Automatic Control, vol. 47, no. 10, pp. 1745–1751, 2002. Bibliography 133 [25] P. D. Christofides, R. Scattolini, D. M. de la Pena, and J. Liu, “Distributed model predictive control: A tutorial review and future research directions,” Computers & Chemical Engineering, vol. 51, pp. 21–41, 2013. [26] M. Razzanelli, E. Crisostomi, L. Pallottino, and G. Pannocchia, “Distributed model predictive control for energy management in a network of microgrids using the dual decomposition method,” Optimal Control Applications and Methods, vol. 41, no. 1, pp. 25–41, 2019. [27] J. B. Rawlings and B. T. Stewart, “Coordinating multiple optimization-based controllers: New opportunities and challenges,” Journal of process control, vol. 18, no. 9, pp. 839–845, 2008. [28] A. N. Venkat, J. B. Rawlings, and S. J. Wright, “Stability and optimality of distributed model predictive control,” in Proceedings of the 44th IEEE Conference on Decision and Control, and the European Control Conference 2005, Seville, Spain, December 2005, pp. 6680–6685. [29] F. Farokhi, I. Shames, and K. H. Johansson, “Distributed MPC via dual decomposition and alternative direction method of multipliers,” in Distributed model predictive control made easy. Springer, 2014, pp. 115–131. [30] A. Bemporad and D. Barcelli, “Decentralized model predictive control,” in Networked Control Systems. Springer, 2010, pp. 149–178. [31] J. Maestre, M. Ridao, A. Kozma, C. Savorgnan, M. Diehl, M. Doan, A. Sadowska, T. Keviczky, B. De Schutter, H. Scheu et al., “A comparison of distributed MPC schemes on a hydro-power plant benchmark,” Optimal Control Applications and Methods, vol. 36, no. 3, pp. 306–332, 2015. [32] A. Jain, A. Chakrabortty, and E. Biyik, “Distributed wide-area control of power system oscillations under communication and actuation constraints,” Control Engineering Practice, vol. 74, pp. 132–143, 2018. [33] A. La Bella, P. Klaus, G. Ferrari-Trecate, and R. Scattolini, “Supervised model predictive control of large-scale electricity networks via clustering methods,” Optimal Control Applications and Methods, vol. 43, pp. 44–64, 03 2021. [34] R. Guicherd, P. A. Trodden, A. R. Mills, and V. Kadirkamanathan, “Superviseddistributed control with joint performance and communication optimisation,” International Journal of Control, pp. 1–12, 2020. [35] E. Masero, J. R. D. Frejo, J. Maestre, and E. Camacho, “A light clustering model predictive control approach to maximize thermal power in solar parabolic-trough plants,” Solar Energy, vol. 214, pp. 531–541, 12 2020. [36] E. Masero, J. M. Maestre, and E. F. Camacho, “Market-based clustering of model predictive controllers for maximizing collected energy by parabolic-trough solar collector fields,” Applied Energy, vol. 306, p. 117936, 2022. 134 Bibliography [37] P. Chanfreut, J. M. Maestre, D. Krishnamoorthy, and E. F. Camacho, “Aladin-based distributed model predictive control with dynamic partitioning: An application to solar parabolic trough plants,” in 2023 62nd IEEE Conference on Decision and Control (CDC), 2023, pp. 8376–8381. [38] P. Chanfreut, J. M. Maestre, A. J. Gallego, A. M. Annaswamy, and E. F. Camacho, “Clustering-based model predictive control of solar parabolic trough plants,” Renewable Energy, vol. 216, p. 118978, 2023. [39] J. R. D. Frejo and E. F. Camacho, “Centralized and distributed model predictive control for the maximization of the thermal power of solar parabolic-trough plants,” Solar Energy, vol. 204, pp. 190–199, 2020. [40] P. Chanfreut, J. M. Maestre, and E. F. Camacho, “Coalitional model predictive control on freeways traffic networks,” IEEE Transactions on Intelligent Transportation Systems, vol. 22, no. 11, pp. 6772–6783, 2021. [41] P. R. Baldivieso-Monasterios and P. A. Trodden, “Coalitional predictive control: Consensus-based coalition forming with robust regulation,” Automatica, vol. 125, p. 109380, 2021. [42] A. Maxim and C.-F. Caruntu, “A coalitional distributed model predictive control perspective for a cyber-physical multi-agent application,” Sensors, vol. 21, no. 12, 2021. [43] P. Chanfreut, J. M. Maestre, and E. F. Camacho, “A survey on clustering methods for distributed and networked control systems,” Annual Reviews in Control, vol. 52, pp. 75–90, 2021. [44] F. Fele, J. M. Maestre, and E. F. Camacho, “Coalitional control: Cooperative game theory and control,” IEEE Control Systems Magazine, vol. 37, no. 1, pp. 53–69, 2017. [45] A. Sánchez-Amores, P. Chanfreut, J. Maestre, and E. Camacho, “Coalitional model predictive control with different inter-agent interaction modes,” European Journal of Control, vol. 68, p. 100676, 2022, 2022 European Control Conference Special Issue. [46] ——, “Robust coalitional model predictive control with negotiation of mutual interactions,” Journal of Process Control, vol. 123, pp. 64–75, 2023. [47] A. Sánchez-Amores, A. Ferramosca, J. M. Maestre, and E. F. Camacho, “Tubebased coalitional model predictive control for tracking,” in Sent to 63nd IEEE Conference on Decision and Control (CDC), 2024. [48] ——, “Partially cooperative robust distributed model predictive control for tracking,” Sent to Automatica, 2024. Bibliography 135 [49] A. Sánchez-Amores, J. Martinez-Piazuelo, J. M. Maestre, C. Ocampo-Martinez, E. F. Camacho, and N. Quijano, “Population-dynamics-assisted coalitional model predictive control for parabolic-trough solar plants,” IFAC-PapersOnLine, vol. 56, no. 2, pp. 7710–7715, 2023, 22nd IFAC World Congress. [50] ——, “Coalitional model predictive control of parabolic-trough solar collector fields with population-dynamics assistance,” Applied Energy, vol. 334, p. 120740, 2023. [51] A. Sánchez-Amores, J. Maestre, and E. Camacho, “Two-layer coalitional model predictive control for parabolic-trough collector fields,” IFAC-PapersOnLine, vol. 56, no. 2, pp. 3146–3151, 2023, 22nd IFAC World Congress. [52] A. Sánchez-Amores, J. M. Maestre, P. A. Trodden, and E. F. Camacho, “A bound on the existence of the maximum jointly invariant set of input-coupled systems,” IEEE Control Systems Letters, vol. 7, pp. 2293–2298, 2023. [53] ——, “A set-theoretic characterization of the existence of the maximum jointly invariant set of input-coupled multi-agent systems,” Sent to Automatica, 2024. [54] P. A. Trodden and J. M. Maestre, “Distributed predictive control with minimization of mutual disturbances,” Automatica, vol. 77, pp. 31–43, 2017. [55] F. Fele, J. M. Maestre, S. M. Hashemy, D. M. de la Peña, and E. F. Camacho, “Coalitional model predictive control of an irrigation canal,” Journal of Process Control, vol. 24, no. 4, pp. 314–325, 2014. [56] F. Fele, E. Debada, J. M. Maestre, and E. F. Camacho, “Coalitional control for self-organizing agents,” IEEE Transactions on Automatic Control, vol. 63, no. 9, pp. 2883–2897, 2018. [57] J. M. Maestre, D. Muñoz de la Peña, A. Jiménez Losada, E. Algaba, and E. F. Camacho, “A coalitional control scheme with applications to cooperative game theory,” Optimal Control Applications and Methods, vol. 35, no. 5, pp. 592–608, 2014. [58] W. Ananduta and C. Ocampo-Martinez, “Decentralized energy management of power networks with distributed generation using periodical self-sufficient repartitioning approach,” in 2019 American Control Conference (ACC), 2019, pp. 3230– 3235. [59] A. Maxim, J. M. Maestre, C. F. Caruntu, and C. Lazar, “Min-max coalitional model predictive control algorithm,” in 22nd International Conference on Control Systems and Computer Science (CSCS), 2019, pp. 24–29. [60] Y. Zheng, Y. Wei, and S. Li, “Coupling degree clustering-based distributed model predictive control network design,” IEEE Transactions on Automation Science and Engineering, vol. 15, no. 4, pp. 1749–1758, 2018.