Control predictivo distribuido basado en la teoría de los juegos

Universidad de Sevilla Departmento de Sistemas y Autom´ atica Doctoral Thesis Distributed Model Predictive Control Based on Game Theory Author: Jos´e Mar´ıa Maestre Torreblanca Supervisors: Dr. Eduardo Fern´andez Camacho Dr. David Mu˜noz de la Pe˜na Sequedo October 22, 2010 Los problemas de los dem´as son a menudo la mitad de tus soluciones La Buena Suerte Alex Rovira Celma y Fernando Tras de Bes The problems of others often hide half of the solution to your problems The Good Luck Alex Rovira Celma y Fernando Tras de Bes Acknowledgements In first place, I would like to thank my family, particularly my parents. I hope they will be glad to know that there is going to be another doctor carrying their surnames. I would also like to thank my academic parents, the doctors Eduardo Fern´andez Camacho and David Mu˜noz de la Pe˜na. They have been very supportive during these years and have taught me how to become a true researcher. I appreciate all their confidence and the countless hours they have spent working with me. There have been many other researchers that have helped me with this work. In particular, I would like to recognize the role of professor Anders Rantzer, who granted me the possibility of working seven months with his team in Sweden. I also want to thank all my colleagues in Seville. Working in the Systems and Automation Department would not be such a privilege without people like you. At last, but not least, I would like to thank my girlfriend, Isabel, for all the time I have not been able to spend with her because of this thesis you have in your hands. Contents 1 Distributed Model Predictive Control and Game Theory 1 1.1 Sensors and actuators networks . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Distributed control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1 Types of control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.2 Types of neighborhood . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.3 Types of communication . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.4 Types of variable update . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.5 Types of agent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3 Model predictive control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Distributed model predictive control . . . . . . . . . . . . . . . . . . . . . . . 9 1.5 Game theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.5.1 Non-cooperative game theory . . . . . . . . . . . . . . . . . . . . . . . 11 1.5.2 Cooperative game theory . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.6 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 i ii CONTENTS 1.6.1 Decentralized schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.6.2 DMPC based on information broadcast . . . . . . . . . . . . . . . . . 17 1.6.3 DMPC based on agent collaboration . . . . . . . . . . . . . . . . . . . 18 1.7 Objectives of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 1.8 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.9 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2 Distributed Model Predictive Control Based on Game Theory for Two Agents 27 2.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.2 Stability properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.2.1 Design procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.3 Simulation examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.3.1 Two double integrators with coupled inputs . . . . . . . . . . . . . . . 42 2.3.2 Application to a supply chain problem . . . . . . . . . . . . . . . . . . 44 2.4 Robustness of the proposed approach against data losses . . . . . . . . . . . . 54 2.5 The four tank process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.5.1 The four tank process . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 2.5.2 Four tank plant model . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 2.5.3 The benchmark problem . . . . . . . . . . . . . . . . . . . . . . . . . . 65 2.5.4 Controllers under test . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 CONTENTS iii 2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 3 Distributed Model Predictive Control Based on Game Theory for Multiple Agents 75 3.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.2 Proposed DMPC controller . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 3.3 Stability....................................... 83 3.4 Controller design procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.5 Example....................................... 94 3.6 Application to a supply chain problem . . . . . . . . . . . . . . . . . . . . . . 98 3.7 Application to control of irrigation canals . . . . . . . . . . . . . . . . . . . . 103 3.8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 4 Distributed Receding Horizon Kalman Filter 109 4.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 4.2 Dual decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 4.2.1 Coordination alternatives for the price update . . . . . . . . . . . . . . 117 4.3 Examples ...................................... 119 4.3.1 Application to mobile robot localization . . . . . . . . . . . . . . . . . 120 4.3.2 Application to traffic and speed control . . . . . . . . . . . . . . . . . 122 4.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4DMPC and Game Theory the control system is hierarchical. In this case, upper layers manage the global objectives of the process and provide references for the lower layers, which control directly the plant. In this thesis we will focus mainly, but not exclusively, on distributed controllers. Note that, as we stated before, the performance of the closed-loop system depends on the decisions that all the agents take, so cooperation and communication policies become very important issues. At this point it is convenient to remark that the models used to design the control system have to be coherent with the type of control adopted. For this reason, models can also be classified in a similar way. 1.2.2 Types of neighborhood In a distributed control system, the state and control actions of a given subsystem may affect other subsystems. Different types of interaction between subsystems can be defined based on the relative degree of interaction between the inputs and states of each subsystem. We speak of an interaction of type zero when two different agents share a state variable or control directly the same manipulated variable. Interaction of type one is referred to the case in which the coupling between the variables controlled by two different agents is strong enough to require communication and coordination between them. In other words, it is not possible to consider the interaction as a mere disturbance. This is the most frequent type of neighborhood in distributed systems. This classification can be extended to a general case. For example, a type two would correspond to the case in which the interaction of a given agent iover an agent jis strong enough to induce a considerable disturbance on a third agent k, which in general will by a type 1 neighbor of agent j. Note that interactions and neighborhood may depend on the way the centralized control problem is distributed between the agents. It is beyond the scope of this work to provide a mechanism to break the centralized models into subsystems. Nevertheless there are works on the literature that deal with this problem, see for example [67, 32]. In this thesis we will deal with the types of neighborhood zero and one. We will assume that the disturbances due to other types of interactions are neglectful. 1.2. DISTRIBUTED CONTROL 5 1.2.3 Types of communication Communication plays an essential role in distributed control. It can be classified in several classes according to different factors related with the protocol that is followed during the communication. In first place, it is possible to distinguish between synchronous and asynchronous communication depending whether there is or not a strict timing in the communication process that determines when an agent can communicate. Second, we can classify the communication as serial or parallel depending on whether one or more than one agent is allowed to communicate at the same time. Finally, we can classify the communication taking into account the number of communications employed. This is an important factor for some systems. It is needless to say that decentralized algorithms have no communication steps. With one communication step it is possible to establish a distributed control scheme in which the agents can inform their neighbors about their plans so that this information considered in the optimization procedure. A consensus or agreement between the agents is only possible if two or more communication steps are used. The number of communications is finite only if the negotiation process converge. However note that it is possible that the number of iterations required to converge is greater than the maximum number of communications allowed. The algorithms that we have developed in this work are designed to be implemented in parallel. Naturally, these algorithms also admit a serial implementation. In general we require synchronous communication although asynchronous communication can be used in our multiple agent distributed control scheme. Finally, the number of communications that is used by the agents is low taking into account that all the decisions are cooperative. For two agents we have developed an scheme with two communication steps and for multiple agents we have proposed an algorithm that offers a good performance with an average number of five communications per agent. 1.2.4 Types of variable update Agents may update their manipulated variables in several ways. It is possible to establish a classification as a function of the way in which the agents update their variables from a social point of view. According to this, we can classify the variable update as exclusive, if each agent update its variables autonomously, shared, if an agent allows other agents to manipulate its variables, or democratical, if an agent takes into account the suggestions of other agents when he calculates the new values of its manipulated variables. The type of variable update that we have considered in this thesis is the democratical one. The decisions are taken in a social way whenever it is possible. 6DMPC and Game Theory 1.2.5 Types of agent The agents are the essence of distributed control systems. There are important features that allow us to classify them according to different parameters. In first place, we can classify them according to the amount of information each agent has about the whole system. This information has a great impact in the distributed control algorithm used. We say an agent is blind if the agent has only information about its corresponding local subsystem. In other words, the agent ignores the rest of subsystems and their influence is modeled as a disturbance. This is the case, for example, of decentralized control systems. An agent is said to be a standard agent if it knows his own dynamics and also how other agents affect him. This knowledge is the base of a possible negotiation procedure with other agents. An agent is samaritan if it is an standard agent that in addition knows how it affects the other agents as well. This information allows the agent to choose his control actions minimizing the possible negative effects on its neighbors. Finally, an agent is omniscient if it has full centralized information or at least it knows all the information it needs about the system. In second place we can classify the agents according with their attitude, which is another important factor. In this context, attitude is related with the will of collaboration between the agents. In this sense the agent’s attitude is noncooperative if the agent behaves selfishly, that is, it only tries to maximize its own utility function (i.e., to minimize its cost function). On the other hand, the agent’s attitude is cooperative when the agent tries to minimize not only its cost but the cost of its neighbors. Thus, the agent can make a sacrifice in terms of his own welfare to help the system reach a better global situation. One of the main assumptions in our control schemes is to consider standard agents from the informational point of view. In addition, the agents are considered to have a cooperative attitude. 1.3 Model predictive control Although there are numerous different control techniques, in this thesis we will focus on Model Predictive Control (MPC). MPC, also known as receding horizon control (RHC), is a popular control strategy for the design of high performance model-based process control systems because of its ability to handle multi-variable interactions, constraints on control (manipulated) inputs and system states, and optimization requirements in a systematic manner. MPC takes advantage of a system model to predict its future evolution starting from the current system state along a given prediction horizon. Nominal MPC controllers consider discrete models of the following form: 1.3. MODEL PREDICTIVE CONTROL 7 x(t+ 1) = f(x(t), u(t)) y(t) = g(x(t), u(t)) where trepresents the sample time and x(t)∈Rn,u(t)∈Rmand y(t)∈Rqare the state and input vectors of the system respectively. The model is used recursively at time tto predict the state or outputs in a finite horizon of length Nwhen a given input trajectory U(t) = {u(t|t), u(t+ 1 |t),...,u(t+N−1|t)}is applied. The resulting state at the end of the horizon, x(t+N|t), is called terminal state. The mathematical properties of the functions fand gdepend on the type of dynamics that are modeled (for example linear, hybrid, nonlinear...). The class of models considered by a given MPC defines the properties of the controller as well as the complexity of the resulting optimization problem. In this thesis we will focus on linear models, although some of the results can be extended to more general classes of systems. One of the most appealing features of MPC is its ability to handle the constraints on the values of the states and inputs in an explicit way. Real processes have limits in the values of all their variables. For example, a valve cannot be opened negatively. Mathematically, this is modeled by defining a set of admitted values for the state, output and input variables, for example: x(t)∈ X u(t)∈ U (1.1) where Xand Uare the sets that define the admissible values for the state and input variables. The objective of model predictive control is to minimize a given performance index that depends on the future predictions of the state, output and input variables. The performance index is a cost function which defines an optimization criterium that is used to determine which control action sequence offers a better performance. Mathematically it is function that expresses the cost associated to a certain evolution of the system in the horizon interval considered. A typical form of the cost function is J(x(t), U(t) = N−1 P j=0 L(x(t+j|t), u(t+j|t)) + V(x(t+N|t)) (1.2) where L(·,·) represents the stage cost of the system at time t+jand V(·) is the terminal cost. Note that the cost function depends on the current state x(t) of the system and the sequence of possible control actions U(t). Therefore, it can be used as a minimization criterium to choose 8DMPC and Game Theory the best possible control sequence for the system. Moreover, the minimization procedure has to include the constraints. This can be seen as if cost function punished those control sequences that result in a violation of any of the constraints. When using continuous time models MPC optimizes over a family of piecewise constant trajectories with a fixed sampling time and a finite prediction horizon so that a finite dimensional optimization problem is obtained. At each sampling time, the function (1.2) is minimized in order to find the optimal control sequence U∗(t). Once the optimization problem is solved, only the first manipulated input value is implemented and the rest of the trajectory is discarded; this optimization procedure is then repeated in the next sampling step. This is the so-called sliding or receding horizon scheme. Once the elements that characterize an MPC controller have been defined, it is possible to summarize the control algorithm followed by this technique: 1. Measure the current state of the system x(t). 2. Calculate which actions provide the best performance over the horizon by solving the following optimization problem: min U(t)J(x(t), U(t)) s.t. x(t+j|t)∈ X,∀j∈[1, N] u(t+j|t)∈ U,∀j∈[0, N −1] (1.3) 3. Apply the optimal inputs calculated for the first time sample of the prediction horizon, that is, u(t|t) and return to step 1. The success of MPC in industrial applications [12] has motivated an important amount of research on the stability, robustness and optimality of model predictive controllers. This success is logical given the advantages of MPC in comparison with its drawbacks. In the positive side we have that MPC is easy to tune if an appropriate model is available. MPC also works very well with control problems that are difficult to solve with other control techniques. For example, it deals naturally with multiple inputs and multiple outputs systems, constraints, delays and disturbances. On the other hand, there are disadvantages that make the implementation of MPC difficult on some systems. Possibly, the main drawback of MPC has to do with its strong computational requirements. Factor as nonlinearities or constraints may make MPC unsuitable for systems with fast dynamics. In particular, nonlinearities lead to non convex optimization problems, which in general are very hard to solve . Likewise, the 1.4. DISTRIBUTED MODEL PREDICTIVE CONTROL 9 performance of MPC heavily depends on the quality of the model used. The worse the model is in comparison with the real system, the poorer results are obtained. As it has been seen, the cost function and the system model have a great impact in the difficulty of the MPC optimization problem. In this thesis we use linear models and quadratic functionals, which allow to solve the optimization problem using quadratic programming (QP). The existing methods for QP problems are very efficient and allow us to solve the optimization problem on line. Moreover, in the case that there is no constraints, the optimization problem can be solved explicitly. 1.4 Distributed model predictive control The drawbacks of MPC hinder its application to large-scale systems. Typical examples of large scale systems are transportation systems such us traffic, water or power networks [78]. In these systems, the computational requirements or the impossibility of obtaining a centralized model are major problems that MPC cannot avoid. Besides large-scale systems, it is also difficult to apply centralized MPC to networked systems; the distributed nature of these systems require control schemes that do not depend on any centralized element. During the last years, we have assisted to the proliferation of networked control systems (NCS) which have emerged from the augmentation of the dedicated local control networks with additional networked (wired and/or wireless) actuator/sensor devices, which have become cheap and easy-to-install [101] [73]. In practice, most large scale and networked control systems are based on a decentralized architecture; that is, the system is divided into several subsystems, each controlled by a different agent which may or may not share information with the rest. Each of the agents implements a controller based on a reduced model of the system and on partial state information, which in general results in an optimization problem with a lower computational burden. As it was shown in the list of taxonomies of section 1.2, different possibilities arise at this point depending on several factors. In first place, the degree of centralized information shared by the agents is important. In general, the more information the agents have, the better control performance can be obtained. However, the amount of information shared by the different agents may have a cost that sometimes outweighs the improvement in the performance. Other issues such as scalability also have to be taken into account when defining access to the model and the state of the system. In second place, the way the centralized cost function is separated between the agents, that is, the way the local cost functions are formed, has relevance too. This factor is related too with the agent’s attitude, which is another important factor that defines the distributed control scheme. Agents can be cooperative or selfish. Local controllers may not be willing to cooperate if cooperation implies that they 10 DMPC and Game Theory have to sacrifice their local performance. Finally, the number of communications that the agents can make in order to take a decision is also a relevant factor. Communication requires time and energy, which may be design limiting factors. These factors generate a family of different distributed MPC problems which demand different solutions adapted to the particularities of each concrete problem. For these reason, different DMPC algorithms have emerged in the literature. Later in the chapter the most important ones will be surveyed. Now it is preferable to see the role that game theory plays in DMPC. 1.5 Game theory Game theory is a mathematical field that studies the process of interactive decision making, that is, situations in which there are several entities, namely players or agents, whose individual decisions determine jointly the final outcome. In other words, game theory translates into mathematical models situations of conflict and/or cooperation between rational and intelligent agents. Given that conflict is present everywhere in the world, its results have been applied to a great variety of fields such as biology, sociology, politics, engineering or economy, to cite just a few examples. The origins of this mathematical field date back to 1921, with the publication of several papers about “la th´eorie du jeu” by the French mathematician ´ Emile Borel. Some years later, in 1928, John Von Neumann published the article “Theory of Parlor Games”, with important contributions that gave the emerging field mathematical respectability [77]. Finally, it is worthy to mention that in 1944 Von Neumann published together with the economist O. Morgenstern the book Theory of Games and Economic Behavior [100], an influential and pioneering work that showed the potential contribution of game theory to economics. The need of game theory tools depends of the degree of interaction between the agents. Interaction is actually a key concept in this context. It is hard to find a situation in which there is no interaction with other entities that induce at least a degree of uncertainty in the decision making process. If this uncertainty is small enough then individual decision making techniques such as classical optimization can be used. For example, the problem of how many goods has to produce a firm in order to maximize its own profit may be seen as an individual decision making problem, but this is just a mere simplification because there are other factors that have more or less influence in the firm’s profit: price of competence products, substitutive goods, or just the general situation of the economy. In the context of game theory, some common words acquire a special meaning; they are still close to their ordinary meaning in most cases but there are some subtleties that cannot 1.5. GAME THEORY 11 be ignored. Next, some of these terms will be introduced and defined: •Game. Situation whose outcome depends on the individual decision of several entities or players. •Players. Each of the two or more entities whose decisions determine the outcome of the game. In this text players will be called indistinctively agents. A player will be said to be rational if he takes his decisions consistently to his objectives, that will be expressed as the optimization of a certain utility function. A player will be said to be intelligent if he knows everything that can be known about the game and can make inferences from that information. For example, in economics agents are assumed to be rational, but not intelligent because a complete understanding of the whole economic model from players is not expected. •Actions. Possible choices a player can make when is his turn to decide. •Strategies. Complete plan of actions of a player for the game, that is, a strategy determines which action must be chosen for a player at any time he has to decide. A profile of strategies is a set of possible strategies for a given player. •Outcome. One of the possible endings or results for the game. Each outcome implies consequences for the players. •Payoff. What is given to a player at the ending of the game, that is, the utility the player attributes to that particular ending. Games can be classified in different categories. The most meaningful one classifies games as cooperative or non-cooperative. Again, the term cooperative has a special meaning in this context. A game is said to be cooperative if agents can negotiate among them and commit themselves to follow common binding strategies. If the agents cannot negotiate, the game is said to be non-cooperative because the agents have no guarantees about the behavior that the rest of the agents will have in the present and the future. The most important and famous results of game theory emerged in non-cooperative game theory. Nevertheless, cooperative game theory provides mathematical tools that suit very well for distributed control problems. For this reason we will review briefly these two subfields. 1.5.1 Non-cooperative game theory Von Neumann and Morgenstern only studied two person non-cooperative games, which are trivial since it is impossible to talk about coalitions in this case. In the case of games with 12 DMPC and Game Theory more than two agents, they focused on the possible coalitions that players would establish. N-person games were not studied from a non-cooperative perspective until the 1940s and 1950s, when the mathematician John F. Nash extended game theory in that direction [68]. In this subsection we will introduce some of the results of noncooperative game theory that have been used in this thesis. The reader interested in a deeper insight on the topic is recommended to see some of the classical references in the literature such as [9, 66]. The most common way to represent a game is the strategic or normal form. This representation is static in the sense that the time instants in which the decisions are taken are not important an it can be assumed that all the agents make their decisions simultaneously. Mathematically, a game Γ in strategic form is given by: Γ = (M, (Ui)i∈M,(Ji)i∈M) (1.4) where M is a nonempty set of agents, Uiis a nonempty set of strategies for agent i and (Ji)i∈Mis a payoff function that maps the set of possible strategies chosen by the players in the set of reals, specifying the payoff for every agent. For simplicity we will use a 2-person game to introduce the most relevant concepts of noncooperative game theory. In this case, it is possible to pose the strategic form of the game as a table. A well known example that can be found in the literature example is the prisoner’s dilemma, which is presented next [77]: Two members of a criminal gang are arrested and imprisoned. Each prisoner is in solitary confinement with no means of speaking to or exchanging messages with the other. The police admit they don’t have enough evidence to convict the pair on the principal charge. They plan to sentence both to a year in prison on a lesser charge. Simultaneously, the police offer each prisoner a Faustian bargain. If he testifies against hist partner, he will go free while the partner will get three years in prison on the main charge. Oh, yes, there is a catch... If both prisoners testify against each other, both will be sentenced to two years in jail. B refuses deal B turns state’s evidence A refuses deal 1 year, 1 year 3 years, 0 years A turns state’s evidence 0 years, 3 years 2 years, 2 years As it can be seen in the table, there are four possible outcomes for this game. Each player has two possible strategies, accept or refuse the deal. However, given that the game is symmetric, if both players are intelligent and rational, only one of the two outcomes in the 1.5. GAME THEORY 13 diagonal will be the final outcome of the game. Each of this solutions has an special meaning in game theory. The first possible outcome that we will analyze is the case in which both player turn state’s evidence. This situation corresponds to the famous Nash equilibrium, which roughly speaking can be defined as the no-regrets outcome, that is, none of the players is unhappy with his decisions during the game once the game has ended. Mathematically a Nash equilibrium is defined as a strategy profile U∗in which no unilateral deviation in strategy by any single player iis profitable for that player, that is, ∀i, ui∈Ui, ui6=u∗ i:Ji(u∗ i, u∗ −i)≥Ji(ui, u∗ −i) where uiand u∗ iare different strategies carried out by the player iand u∗ −istands for the strategy implemented by the rest of agents in the Nash equilibrium. As it can be seen, in this outcome none of the players regrets his choice. For example, let us suppose that A would have choice to refuse the deal, then A would have got the worst possible outcome for his interests. The same analysis can be made by player B. For this reason, both are equally happy with the decisions they made, although, paradoxically, they would have been better if both of them have refused the deal. The other possible logical outcome of the game is the social optimum, which happens when the two players cooperate and refuse the deal. This is related with the concept of Pareto efficiency. Roughly speaking, Pareto optimality has to do with an outcome in which no agent can be made better off without making worse off at least other agent. Mathematically, a payoff vector v= (vi)i∈M) is weakly efficient in the sense of Pareto if there is no other payoff vector w= (wi)i∈M) such that wi≥vi∀i∈M. In addition, if the inequality holds strictly, that is wi> vi∀i∈M, then vis strongly efficient in the sense of Pareto. It can be seen that the payoffs that are Pareto efficient are maximals of the set of possible outcomes. It is important to remark one of the most important features of noncooperative games: the logical outcome of a prisoner’s dilemma played by selfish players is its Nash equilibrium, which paradoxically is worse than the Pareto efficient outcome. To avoid this paradox, it is necessary to introduce some kind of social concern in the payoffs of the players, so that they are not only interested in their own welfare. A distributed control problem can be usually posed as a game of the form 1.4. The process is not straight-forward, though. There is a set of agents whose strategies or control actions have side effects in the costs of the rest of the agents. Nevertheless, defining the cost functions of the agents may not be an easy task. This basic situation may be enhanced by, 20 DMPC and Game Theory L(S1, S1 12, S2, S2 12) = J1(S1, S1 12) + J2(S2, S2 12) + λT(S1 12 −S2 12),(1.11) which can be split in two different functions of the local variables of the agents. Dual decomposition algorithms can be summarized in these steps: 1. Each agent ifinds the values of (Si, Si 12) that minimizes its dual function with a fixed value of the lagrange multiplier λ. φ1(λ) = min S1,S1 12 J1(S1, S1 12) + λS1 12 φ2(λ) = min S2,S2 12 J2(S2, S2 12)−λS2 12 s.t. s.t. S1∈ S1S2∈ S2 (1.12) 2. The distributed problem is reduced to a problem of the complicating variable. The minimum of the primal problem (1.10) is attained at the maximum of the lagrangrian (1.11) with respect λ. The problem max λφ1(λ) + φ2(λ)(1.13) can be solved in a distributed fashion. 3. It can be proved that when λis maximum, the variables S1 12 and S2 12 have the same value. That means that the minimum of (1.10) will be attained when the gradient of λ is zero. In order to obtain the maximum on λ, we can use a distributed gradient search λ=λ−α(S2 12 −S1 12) where αis the step size. 4. Return to step 1 until the enough precision has been obtained, that is, when the update of the variable λis below a given threshold. Note that if the algorithm is stopped before the convergence, it is likely that the solution obtained is not feasible, that is, S1 12 6=S2 12. Nevertheless, this can be solved in part by taking S12 = (S1 12 +S2 12)/2. Dual decomposition has been used for DMPC in [82]. An augmented lagrangian formulation is proposed in [71] and applied to the control of irrigation canals in [70]. The problem of dual decomposition is, in general, the same that primal decomposition has: it requires a great number of iterations to obtain a solution. 1.6. LITERATURE REVIEW 21 Jacobi algorithm This algorithm is an iterative method for the parallel optimization of nonlinear problems and its description can be found in [11] (see pages 219-223). This algorithm is the core idea of one of Venkat’s feasible cooperation-based MPC [99, 93, 83], which is one of the most popular approaches to solve the DMPC problem. Cooperation-based DMPC assumes that the set of coupling variables S12 is empty. The coupling comes in the cost functions of both agents, which depend both on S1and S2. At each sample time, the agents begin an iterative procedure in order to reach a joint solution which can be described as follows: 1. Step 1. Begin the iteration procedure with p= 0. 2. Step 2. At iteration peach agent solves a centralized optimization problem assuming that the neighboring variables have the value of the previous iteration: Sp 1=arg min S1 J(S1, Sp 2)Sp 2=arg min S2 J(Sp 1, S2) s.t. s.t. S1∈ S1S2∈ S2 Sp 2=Sp−1 2Sp 1=Sp−1 1 (1.14) 3. If the algorithm converges, that is, Sp i=Sp+1 i, or a maximum number of iterations has been exceeded, then the iteration procedure stops. Otherwise, increase the iteration index pand go to step 2. Under certain assumptions it can be proved that this procedure is convergent. The sequence of iterates converges to an optimal limit point which coincides with the centralized MPC solution. This solution also holds the condition of Pareto optimality. Note that the algorithm is somehow similar to the primal decomposition one. The difference between them is in the amount of centralized information the agents have. In this case all the agents optimize the centralized cost function with respect their input variables instead of only their particular cost functions. This difference allows the agents not to have to exchange the subgradient of the complicating variables. The price to pay is the extra centralized information that all the agent in the system must have. Descent direction algorithm This algorithm is another iterative method and it is a distributed version of the method of the feasible directions given in [74]. Again it is assumed that the set of coupling variables 22 DMPC and Game Theory S12 is empty. The coupling comes in the cost functions of both agents, which depend both on S1and S2. In addition, it is assumed that the local cost functions include the effect of each agent control actions in the other agent objectives. Thus, the local cost functions are cooperative in their nature. In each time sample the agents begin an iterative procedure in order to reach a joint solution which can be described as follows: 1. Step1. Begin the iteration procedure with p= 0. 2. Step 2. At iteration p, agent one solves the following optimization problem Sp+1 1=arg min S1∇S1J1(Sp 1, Sp 2)T(S1−Sp 1) s.t. S1∈ S1 (1.15) while agent two solves Sp+1 2=arg min S2∇S2J2(Sp 1, S2)T(S2−Sp 2) s.t. S2∈ S2 (1.16) 3. If the algorithm converges, that is, Sp i=Sp+1 i, or a maximum number of iterations has been exceeded, then the iteration procedure stops. Otherwise, increase the iteration index pand go to step 2. As it can be seen, this algorithm follows a gradient search to find the solution for the problem. This method has been applied in [14] to an urban traffic network. A different gradient-based distributed dynamic optimization method is proposed in [89] and applied to an experimental four tanks plant in [5]. The method of [89] is based on the exchange of sensitivities. This information is used to modify the local cost function of each agent adding a linear term which partially allow to consider the other agents’ objectives. 1.7 Objectives of the thesis As it has been seen in the literature review, there is a good number of distributed MPC techniques. Some of them, provide the same performance as centralized MPC, at least theoretically [71, 99]. However, obtaining optimal results is not an easy task and requires in general an intensive use of the communication network. The main objective of this thesis is 1.8. THESIS OUTLINE 23 to provide a good control performance while minimizing the number of communications between the agents and guaranteeing properties such as the stability of the closed-loop system. Likewise, we have tried to minimize the amount of centralized information the agents need. Some sacrifices have to be done in order to reduce the communicational burden of the control schemes developed in this work. The first one is to reduce the number of possible actions that can be implemented by the agents. This reduction aims to limit the disturbances that the agents induce to their neighbors and simplifies the search of a joint solution for the control problem. It is important to stand out that the decisions are taken collectively. Collective decision making proves to be a good strategy even in scenarios in which players are selfish [8]. In addition, we provide design methods that guarantee the stability of the closedloop system for the distributed control strategies that are proposed in this thesis. Moreover, we have developed algorithms to obtain a particular type of invariant set for distributed systems which can be used to analyze decentralized and distributed control schemes. A second goal of the thesis has been to transpose results of distributed control to distributed estimation, which is natural since distributed estimation is needed in distributed applications and because the estimation problem is the dual of the control problem. Finally, our last objective has been to develope techniques based on the dynamical switching of the links of the network, so that links which provide small gains in the control performance are disabled. Based on this work, we have also proposed a method to determine what agents and links are more important in a distributed system. 1.8 Thesis outline The outline of the thesis is the following: •Chapter 2: Distributed Model Predictive Control Based on Game Theory for Two Agents. The interaction between two agents is simplified using an strategic game in which the number of possible control actions is reduced to three for each agent. The three possible choices are to cooperate with the other agent, to behave selfishly or to implement a stabilizing option. Agents choose the best decision of the possible nine combinations trying to optimize the global cost. The stability of the closed-loop system is guaranteed by the controller design method that we propose in the chapter. This control scheme is applied to a supply chain simulation example and a benchmark of the project HD-MPC, the four tank plant. In addition, the robustness of the algorithm is tested against data losses in a stirred tank reactor simulation example. •Chapter 3: Distributed Model Predictive Control Based on Game Theory for N Agents. The combinatorial explosion of the distributed problem hinder the 24 DMPC and Game Theory extension of the previous scheme to a general case. For this reason, a simplified scheme is developed in which each agent sends proposal for the control actions to his neighbors. Neighbors answer quantifying the positive or negative impact of the proposal. Only those proposals that reduce the global cost are accepted. In case that no proposal is accepted a stabilizing control action is applied. The stability of the closed-loop system is guaranteed by the controller design method that we propose in the chapter. Likewise, we define the concept of jointly invariant set and propose a method for its obtention. This control scheme is applied to a supply chain and an irrigation canal simulation examples. •Chapter 4: Distributed Receding Horizon Kalman Filter. Given that estimation is the dual problem of control, it is natural to apply distributed control techniques to the estimation problem In this chapter the estimation problem is reduced to a dynamic programming problem which is distributed between the agents by means of dual decomposition. •Chapter 5: Applications of Cooperative Game Theory to the Control and Estimation of Distributed Systems. Cooperative game theory provides mathematical tools very appropriated for the analysis of situations of conflict in which the agents may establish binding agreements. In this chapter of the thesis some results of coalitional game theory are transposed to distributed control. Based on these results, a control scheme to manage dynamically the links of a network is developed. In addition, a method to analyze the relative importance of links and agents in a distributed system is proposed. •Chapter 6: Conclusions. The thesis ends with a chapter that analyzes the most relevant contributions and, additionally, points out future research lines in the field of distributed systems. 1.9 Contributions As it has been seen in this section, this thesis can be located at the intersection between the fields of distributed control and game theory. In this subsection we will enumerate the major contributions made in this work and the results from the thesis that have been published or submitted to conferences and journals. Given the difficulty of the distributed control problem, it was natural to face a simplified case in the first years of work of the thesis. For this reason, we focused initially in the case of two agents. As a result of this we developed a distributed control scheme with low communications step. This scheme was put to test with a MIMO system described by its 1.9. CONTRIBUTIONS 25 transfer function in [51] and applied to systems described in state space in [53, 52]. The proposed control scheme provides a good performance and guarantees the stability of the closed-loop system under certain assumptions. Besides these three conference papers, this work has been published in an international journal in [54] is compared to other distributed schemes with a real benchmark in a paper submitted to a journal [5]. Once the problem of two agents was solved, we extended our methodology to the general case of systems with Nagents. In order to overcome the difficulties of the combinatorial explosion, a simplified version of the scheme was developed and as a result of this a negotiation algorithm was proposed. As our experiments show, the control scheme provides a good performance in a low number of communication steps. Again, a controller design method that guarantees the stability of the closed-loop system is provided. Additionally, a method to compute the jointly distributed invariant set of the system is presented. This work has been submitted for publication to an international journal [57, 102] and to an international conference [58]. Contributions have also been made to the distribution of the state estimation problem. Concretely, a technique based on the dual decomposition of the state estimation problem was developed together with professor Anders Rantzer, from the LTH in Lund (Sweden). This work has been accepted for publication in [50]. Besides these algorithms, other meaningful contributions to the state of the art have been made. In particular, it is remarkable the novel application of coalitional game theory to distributed control in order to dynamically decompose the agents in coalitions with low interaction. As a result of this we have proposed a distributed control scheme in which the decomposition of the centralized system is done dynamically. Additionally, the cooperative game theory tools provide an interesting method to analyze the relative relevance of agents and links in a network. This work has been presented in [48] and submitted to an international conference [55] and a journal [56]. 26 DMPC and Game Theory Chapter 2 Distributed Model Predictive Control Based on Game Theory for Two Agents As it was stated in the introduction, one of the goals of this thesis is to design distributed cooperative controllers with low communicational burden. In this chapter we focus on the particular case in which the system is divided in two subsystems. In particular, we propose a distributed model predictive control scheme based on a cooperative game in which two different agents communicate in order to find a solution to the problem of controlling two constrained linear systems coupled through the inputs. We assume that each agent only has partial information of the model and the state of the system. In the proposed scheme, the coordination problem between the agents is reduced to a team game where they have to choose one out of three options. To this end, the agents communicate twice each sampling time in order to share enough information to take a cooperative decision. Concretely, we provide sufficient conditions that guarantee practical stability of the closed-loop system as well as an optimization based procedure to design the controller so that these conditions are satisfied. In addition, we study the robustness of the distributed scheme against data losses due to failures in the communication channel. The theoretical results and the design procedure are illustrated using different simulation examples. Finally, we show the results that the control scheme has got in a benchmark of the european project “Hierarchical and distributed model predictive control” (HD-MPC)1. The outline of the chapter is as follows. In section 2.1 the proposed DMPC scheme for two agents is presented. Section 2.2 shows the stability properties of the proposed scheme. 1HD-MPC project, contract number INFSO-ICT-223854. 27 28 DMPC Based on Game Theory for Two Agents Centralized MPC Subsystem 1 Subsystem 2 X1 X2 U1 U2 Figure 2.1: Centralized MPC scheme. The next section shows different simulation examples to illustrate different features of the algorithm. A supply chain and two double integrators coupled through the inputs are used as examples. Section 2.4 tests the robustness of the scheme against communicational failures using a stirred tank reactor. In section 2.5 a real application of the control scheme to a four tank plant is shown. In addition, the controller is compared to other distributed algorithms in this section. Finally, conclusions are presented in section 2.6. This chapter is based on the results and ideas published in [51, 53, 52, 54, 64]. 2.1 Problem formulation Consider the following class of distributed linear systems in which two subsystems coupled with the neighbor subsystem through the inputs are defined x1(t+ 1) = A1x1(t) + B11u1(t) + B12u2(t) x2(t+ 1) = A2x2(t) + B21u1(t) + B22u2(t)(2.1) where xi∈Rni, i = 1,2 are the states of each subsystem and ui∈Rmi, i = 1,2 are the different inputs. This class of systems are of relevance when identifications techniques are used to obtain the transfer function of a process. We consider the following linear constraints in the state and the inputs xi∈ Xi, ui∈ Ui, i = 1,2 2.1. PROBLEM FORMULATION 29 Subsystem 1 Subsystem 2 X1 X2 U1 U2 Decentralized MPC 1 Decentralized MPC 2 Figure 2.2: Decentralized MPC scheme. where Xiand Uiwith i= 1,2 are defined by a set of linear inequalities. The control objective is to regulate the system to the origin while guaranteeing that the constraints are satisfied. Centralized MPC solves a single optimization problem to decide the optimal sequences of the inputs u1and u2with respect to a given performance index based on the full model of the system and on measurements from all the sensors, see figure 2.1. In distributed and decentralized schemes two independent controllers (hereby denoted agents) are defined. Agent 1 has access to the model of subsystem 1, its state x1and decides the value of u1. On the other hand, agent 2 has access to the model of subsystem 2, its state x2 and decides the value of u2. This implies that neither agent has access to the full model or state information and that in order to find a cooperative solution, they must communicate. A control system is decentralized if there is no communication among the agents, see figure 2.2. This is the worst scenario from the performance point of view because each agent has to cope alone with its control problem with the risk that the absence of coordination in the agents decisions may lead to the instability of the system. The control system is distributed if there is communication between agents, see figure 2.3. The degree of communication depends on the control problem and the communication constraints. In this section we present a distributed MPC controller based on a cooperative game scheme between two different agents. The objective of the proposed DMPC scheme is to minimize a performance index that 36 DMPC Based on Game Theory for Two Agents 2.2 Stability properties Controlling a system between two independent agents may lead to an unstable system. The resulting closed-loop system is a multiprocess system and studying the stability of this class of systems is in general a difficult task. Following a terminal region/terminal constraint approach [61, 26], in this section we provide sufficient conditions that guarantee practical stability of the closed-loop system as well as an optimization based procedure to design the controller so that these conditions are satisfied. This result is stated in the following theorem: Theorem 1 Assume that there exist linear feedbacks u1=K1x1and u2=K2x2such that the following conditions hold: F1((A1+B11K1)x1+B12K2x2)−F1(x1) + L1(x1, K1x1)−d1≤0,∀x2∈Ω2(2.9a) F2((A2+B22K2)x2+B21K1x1)−F1(x2) + L1(x2, K2x2)−d2≤0,∀x1∈Ω1(2.9b) x1∈Ω1→(A1+B11K1)x1+B12K2x2∈Ω1,∀x2∈Ω2(2.9c) x2∈Ω2→(A2+B22K2)x2+B21K1x1∈Ω2,∀x1∈Ω1(2.9d) K1x1∈ U1,∀x1∈Ω1(2.9e) K2x2∈ U2,∀x2∈Ω2(2.9f) Ω1∈ X1(2.9g) Ω2∈ X2(2.9h) Then, if at t= 0,Us 1(0), Us 2(0) are given such that Problems (2.2) and (2.3) are feasible for x1,0=x1(0), x2,0=x2(0), U1=Us 1(0) and U2=Us 2(0), then the proposed algorithm is feasible for all time steps t≥0and system (2.1) in closed-loop with the proposed distributed MPC controller is ultimately bounded in a region that contains the origin in its interior. Proof: The proof consists of two parts. We first prove recursive feasibility of Problems (2.2) and (2.3) if at time t,Us 1(t), Us 2(t) are given such that (2.2) and (2.3) are feasible for x1,0= x1(t), x2,0=x2(t), U1=Us 1(t) and U2=Us 2(t). Then we prove that, under the stated assumptions, J(t) = J1(x1(t), Ud 1(t), Ud 2(t)) + J2(x2(t), Ud 2(t), Ud 1(t)) is a decreasing sequence of values with a lower bound. This implies that system (2.1) in closed-loop with the proposed distributed MPC controller is ultimately bounded in a region that contains the origin in its interior. 2.2. STABILITY PROPERTIES 37 Part 1. We will prove this part by recursion. First, we prove that if the state and input trajectories obtained from x1(t−1) and x2(t−1) applying Ud 1(t−1) and Ud 2(t−1) satisfy all the constraints of Problems (2.2) and (2.3), then Ud 1(t) and Ud 2(t) also satisfy all the constraints. Recalling step 5 of the proposed algorithm, to prove this statement it is sufficient to prove that there exists at least a pair of input trajectories that satisfy all the constraints. To this end, we will prove that Us 1(t), Us 2(t) provide a feasible solution for x1(t) and x2(t). Note that in general, it is not possible to guarantee that any of the other options are feasible. Taking into account that by definition Ud 1(t−1) and Ud 1(t−1) satisfy the constraints of Problems (2.2) and (2.3), the following statements hold x1,k ∈ X1, k = 0,...N ud 1,k ∈ U1, k = 0,...N −1 x1,N ∈Ω1 x2,k ∈ X2, k = 0,...N ud 2,k ∈ U2, k = 0,...N −1 x2,N ∈Ω2 where x1,k, x2,k are the k-steps ahead predicted state obtained from x1(t−1), x2(t−1) respectively applying the input trajectories Ud 1(t−1), Ud 2(t−1) defined by ud 1,k, ud 2,k with k= 0,...,N −1. At time step t−1, the first input of the chosen trajectories Ud 1(t−1), Ud 2(t−1) are applied; that is, u1(t−1) = ud 1,0and u2(t−1) = ud 2,0. This implies that x1(t) = A1x1(t−1) + B11u1(t−1) + B12u2(t−1) = A1x1,0+B11ud 1,0+B12ud 2,0=x1,1. Taking into account the definitions of Us 1(t) and Us 2(t), it can be proved that the k-steps ahead predicted state obtained from x1(t), x2(t) respectively applying the input trajectories Us 1(t), Us 2(t) satisfy all the constraints from k= 0 to N−1. Moreover, as x1,N ∈Ω1, x2,N ∈Ω2 it holds that (A1+B11K1)x1,N +B12K2x2,N ∈Ω1 (A2+B22K2)x2,N +B21K1x1,N ∈Ω2 and hence all the constraints of Problems (2.2) and (2.3) are satisfied which implies that Us 1(t), Us 2(t) and hence Ud 1(t), Ud 2(t) provide a feasible solution for x1(t) and x2(t). Taking into account that by assumption, Us 1(0), Us 2(0) satisfy all the constraints for x1(0) and x2(0) and using the above result recursively, the statement of this part is proved. Part 2. In this part we will prove that J1(x1(t), Us 1(t), Us 2(t)) + J2(x2(t), Us 2(t), Us 1(t)) ≤J(t−1) + d1+d2 38 DMPC Based on Game Theory for Two Agents where J(t−1) = J1(x1(t−1), Ud 1(t−1), Ud 2(t−1)) + J2(x2(t−1), Ud 2(t−1), Ud 1(t−1)). Taking into account the definitions of Ud 1(t−1) and Us 1(t) it follows that J1(x1(t), Us 1(t), Us 2(t)) −J1(x1(t−1), Ud 1(t−1), Ud 2(t−1)) is equal to F1((A1+B11K1)x1,N +B12K2x2,N )−F1(x1,N ) + L1(x1,N , K1x1,N )−L1(x1,0, K1x1,0) As L1(x1,0, K1x1,0)≥0 and taking into account (2.9a) and (2.9b), that x1,N ∈Ω1and that x2,N ∈Ω2it follows that J1(x1(t), Us 1(t), Us 2(t)) −J1(x1(t−1), Ud 1(t−1), Ud 2(t−1)) −d1≤0 Following the same steps for J2we obtain that J2(x2(t), Us 2(t), Us 1(t)) −J2(x2(t−1), Ud 2(t−1), Ud 1(t−1)) −d2≤0 and hence J1(x1(t), Us 1(t), Us 2(t)) + J2(x2(t), Us 2(t), Us 1(t)) ≤J(t−1) + d1+d2 As the proposed algorithm chooses Ud 1(t), Ud 2(t) as the pair of input trajectories that yield the minimum cost, it is easy to see that J(t)≤J(t−1) + d1+d2 Following standard Lyapunov arguments and taking into account that recursive feasibility is guaranteed (see the first part of the proof), it is proved that system (2.1) in closed-loop with the MPC controller defined by the proposed controller is ultimately bounded in a region that contains the origin in its interior.  Remark. Theorem 1 guarantees that the closed-loop system is ultimately bounded in a closed region that contains the origin. However, it is possible to prove that the proposed controller provides asymptotic stability if the assumptions 2.9a and 2.9b are modified so that asymptotic stability of the centralized system is guaranteed. In the next chapter, we will study this topic in depth. 2.2. STABILITY PROPERTIES 39 2.2.1 Design procedure In the previous section, we have provided sufficient conditions to guarantee that the closedloop system with the proposed distributed MPC scheme is practically stable. In general, designing the controller parameters so that these conditions are satisfied is a hard problem because the design constraints are coupled; for example, the constraints that define the invariant sets Ω1depend on the set of Ω2and viceversa. For centralized MPC controllers, there are various methods described in the literature on how to design a stabilizing local controller, terminal cost function and terminal region [61, 26, 37] (for example, the local controller and the terminal cost can be obtained solving a LQR problem). These results however cannot be applied to the distributed case. In this section we present an optimization based procedure to find local controllers K1, K2, matrices P1, P2and regions Ω1,Ω2such that (2.9) holds for a given system. The procedure determines first matrices K1, K2, P1and P2such that (2.9a) and (2.9b) hold for any given sets Ω1and Ω2solving a linear matrix inequality (LMI) optimization problem. Once the local feedbacks K1and K2are fixed, the invariant sets Ω1and Ω2are obtained. Note that constants d1and d2are determined a posteriori, once the local feedbacks, terminal costs and terminal regions are fixed. From the point of view of each agent, its neighbor’s input can be viewed as a disturbance. This allows us to use well known tools from control of linear uncertain systems in order to determine a local controller such that a given degree of robustness is guaranteed. In [59, 1, 41] several methods to solve this class of problems are presented. In particular, constraint (2.9a) can be transformed into an LMI and solved using standard techniques, moreover, is equivalent to designing an H-infinity controller for subsystem 1 assuming that u2is the disturbance [41]. The same technique can be followed to design K2and P2. The following theorem defines an LMI constraint that only depends on the system matrices that guarantees that there exist K1, K2, P1and P2such that (2.9a) and (2.9b) hold. Theorem 2 Consider system (2.1). If there exist matrices Wi, Yiand a constant γisuch that the following inequality holds         γiI0BT i,ni 0 0 ∗WiWiAT ii +YT iBT ii WiQ 1 2 iYT iR 1 2 i ∗ ∗ Wi0 0 ∗ ∗ ∗ I0 ∗ ∗ ∗ ∗ I         >0.(2.10) then (2.9a) (or (2.9b), depending on the agent) is satisfied for Pi=W−1 i,Ki=YiW−1 i, and di=γimax x∈Ωni (Knix)TKnix 40 DMPC Based on Game Theory for Two Agents Proof: In [41] it is proved that if (2.10) holds, then the following constraint is satisfied Fi((Ai+BiiKi)xi+Bi,niv)−Fi(xi) + Li(xi, Kixi)−γivTv≤0,∀v(2.11) It follows that (2.9a) (or (2.9b), depending on the agent) holds for di=γimax x∈Ωni (Knix)TKnix  Once the local controllers and the terminal cost functions are fixed, in order to design a distributed MPC scheme that satisfies the assumptions of Theorem 1 one needs to find sets Ω1,Ω2such that (2.9c) to (2.9g) hold. In general this is a difficult problem because each of the sets depends on the other. The size of the terminal region for agent 1 is determined by the magnitude of the disturbances induced by its neighbor agent 2 and viceversa. We provide next an optimization based procedure to solve this problem. In order to present the algorithm we need the following definitions Definition 1 Given the following discrete-time linear system subject to bounded additive uncertainties x+=Ax +Bu +Dw with w∈ W, subject to constraints in the state and the input x∈ X, u ∈ U and a linear feedback u=Kx; a set Ωis said to be a robust positive invariant set for the system if the following constraints hold x∈Ω→(A+BK)x+BKx ∈Ω,∀w∈ W Kx ∈ U Ω∈ X Given system matrices A, B, D, K and the sets X,U,W, there exists several methods to find a set Ω that satisfies these constraints, see for example [39] for a procedure to find the maximal robust positive invariant and [80] for a procedure to find an approximation of the minimal robust positive invariant. We denote Ω(A, B, D, K, X,U,W) the corresponding maximal robust positive invariant set. Taking into account that the input of the neighbor agent can be considered as an unknown bounded disturbance, in order to decouple the computation of the sets Ω1and Ω2, we use the following result based on finding a robust positive invariant set for each subsystem: 2.3. SIMULATION EXAMPLES 41 Theorem 3 Given constants λ1∈(0,1] and λ2∈(0,1], if the sets defined as Ω1= Ω(A1, B11, B12,X1, K1, λ1U1, λ2U2) Ω2= Ω(A2, B22, B21,X2, K2, λ2U2, λ1U1) are not empty, then constraints (2.9c) to (2.9g) are satisfied. Proof: The theorem is proved taking into account the definition of the operator Ω and that λ1U1⊆ U1and λ2U2⊆ U2.  The main idea is that to determine the invariant sets both agents limit its inputs by a factor λi∈(0,1] with i= 1,2 so the other agent can find the maximal robust positive invariant set with respect to a known bounded disturbance. For example, agent 1 finds the maximal robust positive invariant with respect to a disturbance bounded in λ2U2assuming that its input is bounded in λ1U1. Agent 2 does the same. Note that these sets may be empty depending on the value of λ1and λ2. If both sets exists, then they satisfy the stability constraints. In general an infinite number of possible values of λ1and λ2such that both sets are non empty may exist. In order to choose one, we propose to solve the following optimization problem to maximize the feasibility region of the distributed MPC controller: max λ1∈(0,1],λ2∈(0,1] f(Ω1×Ω2) Ω1= Ω(A1, B11, B12, K1,X1, λ1U1, λ2U2) Ω2= Ω(A2, B22, B21, K2,X2, λ2U2, λ1U1) (2.12) where f(·) is a measure of the size of a polyhedra (for example, its Chebyshev radius). Once matrices K1, K2, P1, P2and the sets Ω1and Ω2are determined, constants d1and d2 can be calculated in order to obtain an estimation of the set in which the closed-loop system is ultimately bounded. 2.3 Simulation examples In this section the theoretical results and the design procedure are illustrated using two different examples. The first example is focused on the controller design procedure. The second controller shows the application of the proposed approach to a supply chain problem. The simulations presented in this chapter were performed using Matlab in a computer equipped with a 2.2GHz Core 2 duo processor and 3 GB of RAM memory. 42 DMPC Based on Game Theory for Two Agents 2.3.1 Two double integrators with coupled inputs The system considered is composed by two double integrators with coupled inputs. The first subsystem is defined by the following matrices A11 ="1 1 0 1#, B11 ="0 1#, B12 ="0 0.4# and the second subsystems is defined by A22 ="1 1 0 1#, B22 ="0 1#, B21 ="0 1# The state and the input must satisfy the following constraints: kx1k∞≤1,kx2k∞≤2,|u1| ≤ 1,|u2| ≤ 1 The stage cost functions of each agent are defined by Qi=Iand Ri= 1 for i= 1,2. In order to determine the local controllers Kiand the corresponding weight matrices Pi that define the terminal cost function, a LMI problem based on (2.10) that minimizes the constant γiis solved for each agent. The following matrices are obtained: Ki=h−0.2023 −0.9254i, i = 1,2 P1="32.6719 −17.5149 −17.5149 54.6366 #, P2="38.4509 −5.6447 −5.6447 50.1686 # The last step necessary to apply the proposed algorithm is to determine an invariant region for the two agents, Ω1and Ω2. Different approaches can be used to determine the values of λithat maximize the size of the terminal regions. In this example the terminal region was calculated for a grid with different values of λi. The criterion to select the maximum invariant region was the Chebyshev radius of the maximum ball inside the region. The results were λ∗ 1= 0.3 and λ∗ 2= 0.5. Figure 2.4 shows a 3D plot of the Chebyshev radius as a function of λ1and λ2. The constants λ1and λ2define a trade-off between the degree of freedom that the agents have in order to stabilize the system, and the size of the terminal region which determines the size of the disturbance. As λ2increases, the set defined by K2x∈λ2U2increases. This implies that the set Ω2= Ω(A2, B22, B21, K2,X2, λ2U2, λ1U1) becomes larger because the feasibility region of the input is larger, while the set Ω1= Ω(A1, B11, B12, K1,X1, λ1U1, λ2U2) has to take into account bigger disturbances and may even cease to be defined (i.e., is empty). This 2.3. SIMULATION EXAMPLES 43 0 0.1 0.2 0.3 0.4 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 λ2 λ1 Figure 2.4: Chebyshev radius of the set Ω1×Ω2for different values of λ1and λ2. happens when the minimum positive robust invariant set for an uncertainty bounded in λ2U2 is not included in the feasibility region defined by X1and K1x1∈λ1U1. In figure 2.5(a), inner approximations of the minimum positive invariant sets of subsystem 1 for different values of λ2and a fixed value of λ1are shown. It can be seen that for large values of λ2, the inner approximation is not contained in the feasibility region of agent 2 (shown in red dashed line), and hence, it is empty. In figure 2.5(b) the maximum positive invariant set for the same values of λ1and a fixed value of λ2are shown. It can be seen how the size of the set always increases with λ2. In the first time step, a feasible solution for the centralized problem is used as the shifted trajectories. Simulations results are shown in the next figures for the following initial conditions: x1(0) = "0.7 −1#, x2(0) = "1 0.8# Figure 2.6 shows the trajectories of the states of each agent, the inputs and the cost index. Figures 2.7(b) and 2.7(a) show the state trajectories of each agent along with its corresponding invariant set. 44 DMPC Based on Game Theory for Two Agents −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 x1 1 x2 1 (a) −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 −1.5 −1 −0.5 0 0.5 1 1.5 x1 2 x2 2 (b) Figure 2.5: (a) Minimum robust invariant set for agent 1 as a function of λ2with λ1fixed. (b) Maximum robust invariant set for agent 2 as a function of λ2with λ1fixed. 2.3.2 Application to a supply chain problem In this section, we apply the proposed controller to a reduced version of the MIT beer game and compare the performance with other control schemes. The MIT beer game was developed by Jay Forrester in the late 1950’s to show his managements students how oscillations arise in a supply chain, see for example [92]. A supply chain is the set of structures and processes used by an organization to provide a service or a good to a consumer. Typically three phenomena take place in supply chains flows: oscillation, amplification and phase lag. Due to material or informational delays in the flows of the supply chain, they are prone to oscillation; that is, production and inventories overshoot and undershoot the optimal levels. The magnitude of the fluctuations increase as they propagate from the customer to the factory, with each upstream stage tending to lag behind the downstream one in what is commonly known as the bullwhip effect. The original MIT beer game is composed of four agents: retailer, wholesaler, distributor and factory. Customers demand beer from the first one, who orders beer from the wholesaler, who orders and receives beer from the distributor, who finally orders and receives orders from the factory. There are shipping and processing delays at each stage. In [92], the original model and all the difficulties of the corresponding stock management problem are explained in detail. This problem has been widely used in the literature. In particular, in [23] it has been used as application example for a DMPC scheme. The main difference between the proposed scheme and the DMPC proposed in [23] is that in [23] the agents only communicated once and the only information shared was the future input trajectories (a strategy similar to Iter 1). In this example, a reduced version of the problem with two agents is considered: the 2.3. SIMULATION EXAMPLES 45 0 5 10 15 20 −1 0 1State trajectories of agent 1 x1 1 x2 1 0 5 10 15 20 −2 0 2State trajectories of agent 2 x1 2 x2 2 0 5 10 15 20 −5 0 5Input trajectories u1 u2 0 5 10 15 20 0 100 200 Global cost k J Figure 2.6: State, input and global cost J(t) trajectories of the double integrators in closedloop with the proposed controller. retailer and its supplier. There is no loss of generality since the structure of the game is regular: there is a cascade of firms, each maintaining and controlling its stock. The continuous time equations for the supplier are [92]: ˙sS(τ) = oS r(τ−τ2)−oR r(τ−τ1)−bS(τ)/tb ˙oS u(τ) = oS r(τ)−oS r(τ−τ2)) ˙ bS(τ) = oR r(τ)−oR r(τ−τ1)−bS(τ)/tb (2.13) The equations for the retailer are: ˙sR(τ) = oR r(τ−τ1−τ2) + bS(τ−τ2)/tb−dr(τ−τ1)−bR(τ)/tb ˙oR u(τ) = oR r(τ)−oR r(τ−τ1−τ2)−bS(τ−τ2)/tb ˙ bR(τ) = dR r(τ)−dR r(τ−τ1)−bR(τ)/tb (2.14) The super-scripts R, S denote variables from the retailer and the supplier respectively. Variable si(τ) is the stock level; that is, the number of items available in that stage for shipment downstream. The unfulfilled order of stock oi u(τ) stands for the number of ordered items that the agent is waiting to receive from the upstream stage. The backlog of unfulfilled orders bi(τ) accounts for the number of committed items that have to be shipped to the downstream stage. The parameter tbstands for the average backlog clearance time and introduces a first 52 DMPC Based on Game Theory for Two Agents 0 10 20 30 40 50 60 0 200 400 Retailer Stock Stock reference 0 10 20 30 40 50 60 0 200 400 Supplier Stock Stock reference 0 10 20 30 40 50 60 0 100 200 Retailer Items ordered Extern demand 0 10 20 30 40 50 60 0 100 200 Supplier Items ordered Supplier demand Figure 2.11: Iter1 closed-loop trajectories for scenario 1. 0 10 20 30 40 50 60 0 200 400 Retailer Stock Stock reference 0 10 20 30 40 50 60 100 200 300 Supplier Stock Stock reference 0 10 20 30 40 50 60 0 100 200 Retailer Items ordered Extern demand 0 10 20 30 40 50 60 0 100 200 Supplier Items ordered Supplier demand Figure 2.12: Iter2 closed-loop trajectories for scenario 1. 2.3. SIMULATION EXAMPLES 53 0 10 20 30 40 50 60 100 200 300 Retailer Stock Stock reference 0 10 20 30 40 50 60 200 250 Supplier Stock Stock reference 0 10 20 30 40 50 60 0 100 200 Retailer Items ordered Extern demand 0 10 20 30 40 50 60 0 100 200 Supplier Items ordered Supplier demand Figure 2.13: Iter5 closed-loop trajectories for scenario 1. 0 10 20 30 40 50 60 0 200 400 Retailer Stock Stock reference 0 10 20 30 40 50 60 100 200 300 Supplier Stock Stock reference 0 10 20 30 40 50 60 0 100 200 Retailer Items ordered Extern demand 0 10 20 30 40 50 60 0 100 200 Supplier Items ordered Supplier demand Figure 2.14: Iter10 closed-loop trajectories for scenario 1. 54 DMPC Based on Game Theory for Two Agents does not guarantee an improvement in the performance. It can be seen that Iter5 sometimes is worse than Iter2 from the performance point of view. The simulations demonstrate that the proposed distributed scheme provides a good performance with only two communication cycles because it obtains a cooperative solution; that is, the decision is taken in order to optimize a global performance index. The iterative controllers do not take a cooperative decision and this implies, that in general, the solutions provided are worst. This can be clearly seen in the figures of scenario 1. In this scenario, the centralized MPC is able to react in advance to the demand peak, maintaining the stocks close to the references. The trajectories of the proposed distributed MPC scheme show a larger deviation of the stocks from the references, however, these trajectories do not present oscillations as the trajectories corresponding to Iter1, Iter2 and Iter5. Oscillations are a common result of non-cooperative bargaining processes. In this scenario, Iter10 however provides a better response that the proposed DMPC, at the cost of a high computational burden and a large number of communication steps. 2.4 Robustness of the proposed approach against data losses In this section we carry out a set of simulations to study the robustness of the proposed approach when data losses occur in the control of a stirred tank reactor controller by two agents. The proposed algorithm assumes flawless communications between both agents. In a real distributed environment, errors in the communications and delays in the packets transmission should be expected. For simplicity, we will assume that an error in the communication link will affect the transmissions in both ways, so there is no possibility that only one of the agents is affected by the error. To test the effect of data losses in the closed-loop system performance,, we assume that the probability of flawless communications is given by the parameter reliability ∈[0,1]. This parameter characterizes the quality of the communication network. The higher it is, the better for the communication. In this section we show the results of simulations corresponding to different values of the parameter reliability. In the original algorithm the agents chose among three options for the control signal (U0 i, U∗ i,Uw i) with the goal of minimizing J. When data losses occur, the agents do not receive Uw ior the information needed to build the global cost table. In this case, each agent must decide whether to keep applying the last optimal input trajectory U0 i, or behave selfishly and try to minimize its local cost function choosing U∗ i. In order to test the robustness of the proposed approach on the worst possible case, we assume that when communication errors occur each controller operates in a decentralized way, that is, applying U∗ i. 2.4. ROBUSTNESS OF THE PROPOSED APPROACH AGAINST DATA LOSSES 55 Figure 2.15: Continuously stirred tank reactor (CSTR) Remark 1 Note that as the parameter reliability tends to zero, the amount of information shared by the agents decreases, and the controller tends to operate in a decentralized manner. To demonstrate the robustness of the proposed controller against communication errors, we use the linearized model of a continuously stirred tank reactor (CSTR) presented in [12]. The linearized process around a given equilibrium point is described in continuous time by the following transfer matrix2 "y1(s) y2(s)#="1 1+0.7s 5 1+0.3s 1 1+0.5s 2 1+0.5s#"u1(s) u2(s)#, where the manipulated variables u1and u2are respectively the flow rate and the flow of coolant in the jacket. The controlled variables y1and y2are respectively the effluent concentration and the reactor temperature, see figure 2.18. The sampling time is defined as Ts= 0.03s. The control objective is to track a given constant reference from a random initial state. We first design a centralized MPC scheme for comparison purposes. This controller decides both inputs simultaneously. For this reason, the MPC optimization problem that has been used for the simulations is based on minimizing the following cost function using the linearized discrete model of the process J= N−1 P k=0 (ref1(k)−y1(k))TWy,1(ref1(k)−y1(k)) +(ref2(k)−y2(k))TWy,2(ref2(k)−y2(k)) +∆u1(k)TW∆u,1∆u1(k) +∆u2(k)TW∆u,2∆u2(k) 2The notation x(s) refers to the Laplace transform of the signal x. 56 DMPC Based on Game Theory for Two Agents where refi(k) is the reference signal for the controlled variables y1and y2. Note that the cost function depends on the predicted values of the inputs and outputs. In particular, x(k) stands for the predicted value of the signal x k steps ahead in the future. The following values were used for the controller parameters: N= 5 ref1= 0 ref2= 0 Wy,1=Wy,2= 1 W∆u,1=W∆u,2= 0.05 (2.16) Note that we have not considered constraints on the input or the outputs in the simulations we have performed to test the robustness of the distributed scheme with respect to data losses. In order to analyze the performance of our distributed scheme, we will make a comparison with a decentralized and a centralized controller. Note that the centralized controller provides the optimal solution from point of view of the performance while the decentralized one allows to determine what to expect in the absence of communication. The comparison will be made based on the following parameters: •λ: Convergence rate of the global cost function. This parameter is computed as the smallest value λsuch that the following inequality holds J(kTs)≤J0·λk, λ > 0 (2.17) where J(kTs) is the value of the global cost function evaluated at time t=kTsfor the applied input trajectories and J(0) is its initial value, that is, at time t= 0. If λ > 1 the controlled system is unstable. •tr: Rise time. Number of sample times required in order to get a relative error below 5%, where the relative error is defined as Eri = refi−yi refi·100.(2.18) In first place, we will focus on the centralized controller. More than twenty simulations of the closed-loop system with the centralized controller were done, all of them beginning with different initial states. Half of the simulations were done for a number of kmax = 100 time 2.4. ROBUSTNESS OF THE PROPOSED APPROACH AGAINST DATA LOSSES 57 samples and the other half with kmax = 300. The average performance parameters obtained for the centralized controller were: λ= 0.77 Jv = 13. In order to apply decentralized and distributed MPC schemes, we considered that the CSTR is controlled by two different agents. Agent 1 controls the flow rate u1based on the measurements of the y1, while agent 2 controls u2based on the measurements of y2. Each agent has an incomplete model of the system; that is, they only know the first row of the system model (how their measured output is affected by each of the inputs). A decentralized MPC scheme is based on the idea that each agent tries to control its own subsystem without communicating with the other agent. Each agent tries to minimize a local cost function. For agent 1 the local cost function is J1= N−1 P k=0 (ref1(k)−y1(k))TWy,1(ref1(k)−y1(k)) + ∆u1(k)TW∆u,1∆u1(k) and for agent 2 the local cost function is: J2= N−1 P k=0 (ref2(k)−y2(k))TWy,2(ref2(k)−y2(k)) + ∆u2(k)TW∆u,2∆u2(k). At each time step, agent 1 receives y1and finds the optimal sequence of inputs such that J1 is minimized assuming that u2= 0. Agent 2 follows the same protocol. For this particular system the decentralized controller is not able to stabilize the system. These simulations demonstrate that even a simple system may become unstable when the control agents are not coordinated by a proper scheme. After carrying out the simulations of the centralized and decentralized controllers, we focus on the proposed DMPC scheme. As in the decentralized scheme, each agent tries to minimize its corresponding local cost function Ji. In this case, however, agents do communicate and try to minimize the sum of their optimization functions following the proposed DMPC scheme. In first place we consider the case in which there are no data losses or delays (reliability = 1). As it can be seen in figure 2.16 , the proposed controller scheme is able to stabilize the closedloop system. We carried out over twenty simulations with different initial states and constant references. For this set of simulations the performance parameters were tr= 41.1458 and λ= 0.8858. It can be seen that the performance of the distributed scheme is worst than the one of the centralized controller (although much better than the decentralized scheme which is not able to stabilize the system). As mentioned before, the centralized scheme is the best possible controller from the communication point of view. 58 DMPC Based on Game Theory for Two Agents 0 5 10 15 20 25 −20 0 20 40 60 k U u1 u2 0 5 10 15 20 25 −10 −5 0 5 k Y y1 y2 0 5 10 15 20 25 0 50 100 150 200 k J J1 J2 Figure 2.16: Trajectories of the system in closed-loop with the proposed DMPC and reliability = 1. Table 2.6: λand trfor reliability ∈[0.1,0.9]. reliability 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 λ, kmax = 100 0.89 0.89 0.90 0.93 0.99 1.17 1.39 1.60 1.93 tr, kmax = 100 44.63 48.02 64.13 76.57 84.15 98.86 100 100 100 λ, kmax = 300 0.95 0.96 0.96 0.97 0.98 1.16 1.38 1.60 1.92 tr, kmax = 300 45.37 52.81 71.31 90.93 129.62 291.03 300 300 300 In order to test the robustness of the proposed DMPC with respect to communications errors, a set of simulations with different reliability values were carried out. for each value reliability over 20 simulations were done. The results obtained are shown on table 2.4. Note that the value of kmax affects the value of the comparison parameters. An increment in the value of λis found when kmax increases if the system stays stable (λ < 1). This is logical given the definition of the parameter λ. For example, if the system has reached the reference after k1 sample times, the evolution during the rest of the time steps until the end of the simulation won’t be significant. Thus, it can be concluded that the time samples after the system has reached the reference only degrade quantitatively the value of λ. Note too that trgrows with kmax. Again, this is expected because when the controller is unable to regulate the system to the reference, the value tr=kmax is taken. Hence, it is obvious that the kmax has influence over tr. Finally, note how the value of the performance parameters improves as the parameter reliability grows. In particular, it can be seen how the value tends to the one obtained with the flawless communication simulations. The simulation results also show that depending on 2.5. THE FOUR TANK PROCESS 59 0 5 10 15 20 25 −500 0 500 1000 k U u1 u2 0 5 10 15 20 25 −150 −100 −50 0 50 k Y y1 y2 0 5 10 15 20 25 0 1 2 3 4 5x 104 k J J1 J2 Figure 2.17: Trajectories of the system in closed-loop with the proposed DMPC and reliability = 0.5. the value of reliability, the DMPC is able to stabilize the closed-loop system or not. For reliability ≥0.5, the performance parameter is λ < 1, which implies that the closed-loop system is stable. However, if more than 50% of the communications fails, then the DMPC is not able to stabilize the closed-loop system. These results show that when communication network becomes faulty, the proposed controller tends to operate in a decentralized manner, and hence, is not able to stabilize the system. 2.5 The four tank process In this section we show experimental results of the proposed controller in a four tank process, which is one of the benchmarks of the european project HD-MPC. The physical plant is situated in facilities of the University of Seville and was presented in [6]. Different universities are working in the project as Delf Institute of Technology (Netherlands), Aachen University (Denmark) and Universidad Nacional de Colombia. In this chapter we present the results of the benchmark for the different DMPC policies developed by these universities. The results presented in this section have been submitted for publication in a joint work with the rest of the benchmark participants [5]. The use of benchmarks is useful for evaluating the capabilities of different approaches to control systems for real problems. Benchmarks allow to test, evaluate, and compare different control solutions at real or simulated plants. The research and the industry community 60 DMPC Based on Game Theory for Two Agents 0 5 10 15 20 25 −5000 0 5000 10000 k U u1 u2 0 5 10 15 20 25 −1500 −1000 −500 0 500 k Y y1 y2 0 5 10 15 20 25 0 2 4 6x 106 k J J1 J2 Figure 2.18: Evolution of the system with reliability = 0.2 benefit from these activities since the design of a good simulation test bed is often time and resource consuming. Furthermore, good simulation test beds are often subjected to heavy criticism as they either cover only a narrow part of the problem or they are purposely designed to get biased rather than objective performance results. The benchmark examples would effectively overcome these problems by: a) allowing an objective evaluation of control, b) reducing resources and time spent on developing validation models, c) giving researchers the possibility to evaluate their proposals on a variety of cases, and d) opening up a forum to compare the performance of various solutions and to discuss the quality of the results. The four tank process has proven to be a very interesting system for control education and research [35]. The main characteristic of this process is that it is a simple multivariable system with highly coupled nonlinear dynamics that can exhibit transmission zeros dynamics. The four tank system has been used to illustrate advanced control strategies [19] such as internal model control and dynamic matrix control [25], multivariable robust control [95] and distributed MPC [62]. In addition, it has also been utilized as an educational tool in teaching advanced multivariable control techniques. 2.5.1 The four tank process The four tank process that we have in the university of Seville is an educational plant designed to test control techniques using industrial instrumentation and control systems. The plant is a hydraulic process of four interconnected tanks inspired by the educational quadruple tank 2.5. THE FOUR TANK PROCESS 61 process proposed by Johansson [35]. The main characteristic of this process is that it is a simple multivariable system with highly coupled nonlinear dynamics that can exhibit transmission zeros dynamics. The four tank plant retains the structure of Johansson’s process (see Figure 2.20(a)), but has been modified to enable different configurations and interconnections of the tanks. Figure 2.19: A photo of the four tank plant A photograph of the plant can be seen in figure 2.19 and a schematic plot of the plant is shown in Figure 2.20(b). The inlet flow of each tank of the plant is measured by an electro-magnetic flow-meter (Siemens Sitrans FM Flowsensor 711/S and transmitters Intermag/transmag) and regulated by a pneumatic valve (Siemens VC 101 with a positioner Sipart PS2 PA). This allows the plant to emulate the three-ways valve of Johansson’s quadruple tank process by providing suitable set-points to the flow controllers. The level of each tank is measured by means of a pressure sensor (Siemens Sitrans P 7MF4020 and 7MF4032). All the measurements and commands are 4-20 mA current signals and these are connected to a PLC (Siemens S7-200). The output flow of each tank can also be adjusted by means of a manual tuning valve which allows to adjust the speed of the dynamics of the plant. In order to achieve a safe operation of the plant and to prevent the overflow of tanks, each tank has a high level switching sensor used as alarm to switch off the pumps. 68 DMPC Based on Game Theory for Two Agents the value of Usto a feasible solution. This solution is obtained solving a feasibility problem, in particular an LP, based on the full model of the system. Note that for this particular benchmark, no terminal region has been taken into account. The weighting matrices were chosen to minimize the benchmark objective function, that is, Q1=Q2=I,R1=R2= 0.01. The local controller gains for each controller were K1= [0.17 0.21] and K2= [−0.16 −0.14]. These gains were designed with LMI techniques based on the full model of the system in order to stabilize both subsystems independently while assuring the stability of the centralized system. The role of these gains is important because the option in the game that allows to guarantee closed-loop stability is constructed shifting the last decided control action; that is, the first element is dropped after it is applied in the system and a term evaluated with these gains is added at the end of the horizon control vector, see [54] for more details. The proposed distributed MPC controller only needs three communication steps in order to obtain a cooperative solution to the centralized optimization problem, has low communication and computational burdens and provides a feasible solution to the centralized problem. The simulation and experimental results show that the distributed scheme is able to control the system. The designed controller has been sucessfully tested on the real plants and the trajectories are shown in Figure 2.22. The performance index of the test is J= 29.5787. The performance index is close to the performance index of the centralized MPC for regulation. Note however that the input trajectories are not smooth because the controllers switch between different modes during the simulation. Other controllers In the benchmark, the subsystems have been chosen according to the pairings derived from the relative gain array (RGA) analysis. Considering the values of the RGA the sensible pairing is to control the output h1with qb(y1with u2) and h2with qa(y2with u1). In first place, we mention the only decentralized controller that was tested. In particular a decentralized MPC for tracking was implemented, that is, a MPC for tracking [44] was designed for each subsystem. A communication based DMPC based on nonlinear dynamic optimization was tested as well. This controller is based on nonlinear dynamic optimization. The optimization is based on an iterative procedure of information broadcast in which the two local controllers exchange the value of the interaction variables. A more sophisticated cooperative version of the last algorithm was also tested. Concretely, the gradient-based distributed dynamic optimization (GBDDO) method [89] was put to test. Besides the communication of interac- 2.5. THE FOUR TANK PROCESS 69 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000 0.2 0.4 0.6 0.8 1 1.2 Level [m] y1 s1 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000 0.2 0.4 0.6 0.8 1 1.2 Level [m] y2 s2 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000 0 2 4 Flow [m3/h] qa qb 3000 4000 5000 6000 7000 8000 9000 10000 11000 12000 0 0.5 1 time [s] Level [m] h3 h4 Figure 2.22: Evaluation of the control test in the real plant of the DMPC based on a cooperative game tion variables, the GBDDO-based MPC requires a calculation and exchange of sensitivities. This information is used to modify the local cost function of each agent adding a linear term which partially allow to consider the other agents’ objectives. The scheme proposed in [70, 72], which consists on a serial DMPC based on dual decomposition, was tested as well. This scheme is derived from a serial decomposition of an augmented Lagrangian formulation of the centralized overall MPC problem. Finally, it is remarkable that another controller based on game theory concepts was also implemented. Concretely, a feasible-cooperation distributed model predictive controller based on bargaining game theory concepts was tested. This approach models the DMPC as a game. According to [69], the solution of the cooperative game associated with the DMPC problem is given by a control vector which maximizes the product of the difference between the corresponding cost for each system and the cost of the subsystem when it does not cooperate with the rest[29, 76]. This last cost is computed through a minmax problem and is called disagreement point. The solution for this problem is computed in a distributed fashion following the same algorithm used in [98]. Evaluation of the controllers Table 2.8 shows some qualitative properties of all the controllers that took part on the HDMPC benchmark. The entry Model Requirements shows whether the controllers need full or partial knowledge of the system and whether the model used is linear or nonlinear. The entry Control Objectives shows whether the controller is optimal from a centralized point of view (i.e., provides the same solution as the centralized MPC for regulation), guarantees 70 DMPC Based on Game Theory for Two Agents constraint satisfaction if a feasible solution is obtained and whether it can be designed to guarantee closed-loop stability in a regulation problem. The Auxiliary software entry shows which type of additional software is needed by each controller of the distributed scheme. The centralized controller is based on a linear model of the full plant and is included as a reference for the performance of the distributed MPC schemes. On the other hand, the decentralized controller provides a reference on what can be achieved with no communication among the controllers at all. The other controllers assume that each controller has access only to its local state and model. All the controllers but the DMPC based on dynamic optimization are based on a linear model. From the control objectives point of view, Table 2.8 shows whether the controller is optimal from a centralized point of view (i.e., provides an optimal solution to the corresponding centralized MPC), considers state and input constraints and whether it can be designed to guarantee closed-loop stability in a regulation problem. The decentralized controller considered cannot guarantee optimality, constraint satisfaction, nor stability. Note that in order to guarantee closed-loop stability, the the DMPC controller proposed in this thesis needs full model knowledge as we have mentioned in this chapter. The distributed controllers that guarantee optimality (provided sufficient evaluation time) are the Serial DMPC and the DMPC based on dynamic optimization with GBDDO. Note that this controllers are also the ones with a larger communication and computational burden. Another key issue in distributed schemes is the class of computational capabilities that each controller must have. In particular, for the schemes considered each controller must be able to solve either QP problems or general nonlinear optimization problems. In the experiments, the controllers used MATLAB’s optimization toolbox, in particular quadprog and fmincon. Evaluation of the experimental results The experimental results demonstrated how a centralized solution provides the best performance while the performance of a fully decentralized controller is worse. Distributed schemes in which the controllers communicate in general improve this performance, although the experimental results also demonstrated that a distributed MPC scheme is not necessarily better than a decentralized scheme and it depends on the formulation of the controller and its design. It is also clear how those controllers that incorporate offset-free techniques provide a better performance index. In order to obtain a measure of the performance without the effect of 2.5. THE FOUR TANK PROCESS 71 Qualitative Indexes Model Requirements Control Objectives Auxiliary Software Centralized Regulation Linear system Full model Optimal Constraints Stability QP Decentralized Linear system Local model Suboptimal QP DMPC Coop. game Linear system Local model (Full model) Suboptimal Constraints (Stability) QP DMPC D.O. (GBDDO On) Nonlinear system Local model Optimal Constraints NLP DMPC D.O. (GBDDO Off) Nonlinear system Local model Suboptimal Constraints NLP DMPC Bargaining game Linear system Local model Suboptimal Constraints NLP Serial DMPC Linear system Local model Optimal Constraints QP Table 2.8: Table of qualitative benchmark indexes of each tested controller. 72 DMPC Based on Game Theory for Two Agents the steady offset, the transient performance index Jthas been calculated. This index is evaluated computing the cumulated cost during the transient. The entry tsshows cumulated the settling time of the three reference changes. This shows that the offset-free controller (MPC for regulation) has a transient performance index similar to the total performance index while for the rest of the controllers, the transient index is better. The effect of the communication between controllers on the calculation of the control inputs can be seen comparing the decentralized MPC with the DMPC based on a cooperative game. The size of the QPs to be solved at each sampling time is similar in both cases, while the computational time is larger in the DMPC case. The DMPC based on dynamic optimization exhibits the largest computational time due to the NLP to be solved at each iteration and the number of iterations. All the controllers were implemented using a MATLAB function and were not designed to optimize the evaluation time. For this reason, the computation time has not been taken into account. In particular, both DMPCs based on dynamic optimization had a computation time lower than ten seconds, while the rest were of the order of one second. These computation times were lower than the sampling time chosen for each controller and moreover, they could be dramatically reduced using an appropriate implementation framework. Motivated by these issues, the computational burden is best measured on the number and size of the optimization problems solved at each sampling time. The centralized schemes solve a single QP problem with 2Noptimization variables while the decentralized controller solves 2 QP problems with Noptimization variables. The difference in the computational burden between these schemes grows with the prediction horizon and the number of controllers. Distributed schemes try to find a tradeoff between computational/communicational burden and optimality. The DMPC based on a cooperative game and the DMPC based on a bargaining game solve a fixed number of low complexity optimization problems. DMPC based on dynamic optimization and Serial DMPC provide optimality at the cost of a higher computational burden. In this particular benchmark, the best results are provided by the DMPC based on a cooperative game, however, there are several issues that must be taken into account. First of all, because the controller chooses among nine different modes of operation, the resulting input trajectories are not smooth. Figure 2.22 shows how the input seems to switch among at least two optimal trajectories. Depending on the application, this switching may not be acceptable. In addition, this control scheme is specially designed for only two controllers, because the number of possible modes grows in a combinatorial way with the number of controllers. 2.6. CONCLUSIONS 73 Quantitative Indexes J Jtts# floats # trans Centralized Regulation 25.46 23.78 2735 N.D N.D. Decentralized 39.54 21.2 1685 0 0 DMPC Coop. game 30.71 28.19 2410 20 3 DMPC D.O. (GBDDO On) 33.91 33.36 2555 150 5 DMPC D.O. (GBDDO Off) 35.65 34.63 1700 75 5 DMPC Bargaining game 46.32 39.52 3715 6 2 Serial DMPC 44.59 41.94 3130 Table 2.9: Table of the quantitative benchmark indexes of each tested controller 2.6 Conclusions In this chapter we have proposed a novel distributed MPC algorithm based on game theory for a class of systems controlled by two agents. The proposed controller only needs two communication steps in order to obtain a cooperative solution to the centralized optimization problem. Each agent solves an optimization problem that only depends on its local model and partial state information. After sharing information about the local cost, the agents choose the solution that yields the best global performance among a set of suboptimal possibilities. The options are suboptimal because each agent has an incomplete view of the system and they propose the best solutions from their point of view. The proposed algorithm has low communication and computational burdens and provides a feasible solution to the centralized problem. In addition, we provide sufficient conditions that guarantee practical stability of the closed-loop system as well as an optimization based procedure to design the controller so that these conditions are satisfied. Examples and real experiments have shown the properties the good performance of the controller, specially taking into account its low communicational and informational requirements. The robustness of the proposed scheme against failures in the communication channel has been proved as well, at least when the probability of failure is lower than a fifty percent. Finally, it is worthwhile to mention that the scheme has been tested and compared with other distributed algorithms in a benchmark of the european project HD-MPC. 74 DMPC Based on Game Theory for Two Agents Chapter 3 Distributed Model Predictive Control Based on Game Theory for Multiple Agents In the previous chapter we presented a distributed scheme based on a cooperative game for the particular case in which the system is controlled by two agents. Unfortunately, the complexity of that scheme grows exponentially with the number of agents. In this chapter we propose a distributed model predictive control scheme based on agent negotiation suitable for problems in which the number of control agents is greater than two. Once more, we consider the control of several subsystems coupled through the inputs by a set of independent agents that are able to communicate and we assume that each agent has access only to the model and the state of one of the subsystems. This implies that in order to take a decision which is cooperative from a global point of view, i.e. for the whole system, the agents must negotiate. At each sampling time, following a given protocol, agents make proposals to improve an initial feasible solution on behalf of their local cost function, state and model. These proposals are accepted if the global cost improves the cost corresponding to the current solution. In addition, we study the stability properties of the proposed distributed controller and provide precise conditions based on a new concept of invariance for distributed and decentralized systems that guarantee that the closed-loop system is practically stable along with an optimization based controller and invariant design procedure. The theoretical results and the design procedure are illustrated using different simulation examples. In particular, we use an academical example to show the main theoretical contributions of this chapter. Next, we test the scalability of the distributed scheme with supply chains composed by an increasing number of nodes. The last example consists on the application of the proposed scheme to the control of irrigation canals. 75 76 DMPC Based on Game Theory for N Agents The outline of the chapter is as follows. In section 3.1 the problem is formulated. Section 3.2 defines the proposed DMPC controller. Stability is studied in section 3.3 and a design method is given in section 3.4. An academical example to show the theoretical properties of the controller is presented in section 3.5. Section 3.6 deals with the application of the controller to a supply chain problem and section 3.7 studies its application to an irrigation canal problem. Finally, conclusions and future work are presented in section 3.8. This chapter is based on the results and ideas submitted for publication in [57, 102, 103, 58]. 3.1 Problem formulation We consider the following class of distributed linear systems in which there are Mxsubsystems coupled with their neighbors through Muinputs xi(t+ 1) = Aixi(t) + Pj∈niBijuj(t) (3.1) where xi∈Rqiwith i= 1,...,Mxare the states of each subsystem, and uj∈Rrjwith j= 1,...,Muare the different inputs1. The set of indices niindicates the set of inputs uj which affect the state xiand the set of indices mjindicates the set of states xiaffected by the input uj. We define mathematically the concept of neighborhood of agent ias Ni:= [ j∈ni mj.(3.2) Therefore, any agent jincluded in Niis a neighbor of agent i. Note that this does not imply that i∈Nj, that is, the neighborhood is not a symmetrical property in this context. We consider the following linear constraints in the states and the inputs xi∈ Xi, i = 1,...,Mx uj∈ Uj, j = 1,...,Mu (3.3) where Xiand Ujare closed polyhedra that contain the origin in their interior defined by the following set of linear inequalities xi∈ Xi↔Hxixi≤bxi, i = 1,...,Mx uj∈ Uj↔Hujuj≤buj, j = 1,...,Mu (3.4) Note that, as these polyhedra contain the origin in their interior, then bxi>0 and buj>0. 1Throughout this chapter the time dependence is omitted when possible for notational convenience. 3.1. PROBLEM FORMULATION 77 x 1 x Mx u 1 u Mu System Controller MPC Figure 3.1: Centralized MPC. There are many different physical systems that can be modeled under this formulation. For example, in [14] this model is used to represent the dynamics of a traffic network. In [70] the dynamics of an irrigation canal system are described with a similar formulation. In [54] the beer game, a typical supply chain problem, is described likewise. This class of systems can be represented by a graph in which to each node either the state of one of the subsystems or one of the inputs available is assigned, and the arcs connect the inputs to the states they affect. The control objective is to regulate the states of all the subsystems to the origin while satisfying the state and input constraints. To this end, centralized MPC follows a receding horizon approach and at each sampling time obtains the current states and solves a single finite horizon optimal control problem based on a performance index that depends on all the states and inputs. See figure 3.1 for a scheme of a centralized MPC controller. In distributed MPC schemes there are several agents that decide all the control inputs. It can be seen that although the states are not dynamically coupled, the agents need to negotiate in order to decide the value of the shared inputs. There are many possible distributed schemes depending on the available information and communication constraints. Figure 3.2 shows a scheme of a distributed controller in which each agent has access to partial state information and can communicate with the rest of the agents. This is the class of distributed control scheme considered in this work that is presented in the next section. Remark 2 One of the differences between the proposed approach and other cooperative MPC schemes is that the agents do not have a global model of the system. This may be important in some applications in which the centralized model is not available or the agents do not want to 84 DMPC Based on Game Theory for N Agents The requirements of Assumption 1 are twofold, first, the local feedbacks must satisfy constraint (3.15a) which implies that the system in closed-loop with these set of local controllers is stable. Second, sets Ωisuch that (3.15b) to (3.15d) are satisfied must exist. We denote these sets as jointly positive invariants for system (3.1) in closed-loop with the controllers defined by matrices Kij. It is important to note that although the cartesian product of these sets is a positive invariant of system (3.1), in general it is not possible to obtain the jointly positive invariant sets from an invariant set of system (3.1) obtained following standard procedures because each Ωiis defined only in a subspace of the whole state space; that is, in the space corresponding to the state xi. This property is necessary in order to define for each agent a set of constraints that depend only on its state, and hence, only on its model. See the constraints of problem (3.9). Theorem 4 If Assumption 1 holds and at time step t= 0,Us(0) is given such that each of the Mxoptimization problems (3.9)2are feasible for xi,0=xi(0) and Uj=Us j(0) with i= 1, .., Mxand j∈ni, then the proposed algorithm is feasible for all time steps t≥0and system (3.1) in closed-loop with the proposed distributed MPC controller is asymptotically stable. Proof The proof consists of two parts. We first prove that there is always a proposal which satisfies all the constraints (3.9) and then we prove that, under the stated assumptions, J(t) = Mx X i=1 Ji(xi(t),{Ud j(t)}j∈ni) (3.16) is decreasing sequence lower-bounded by zero. Part 1. Taking into account that {Ud j(t−1)}j∈nisatisfies all the constraints of (3.9) and Assumption 1, it is easy to prove that {Us j(t)}j∈niprovides a feasible solution for xi(t). It follows, that Ud(t) provides a feasible solution for the optimization problem of agent ibecause it is chosen among a set of proposals which are required to be feasible in order to be accepted. Note that a proposal which is unfeasible for any of the agents cannot be chosen because the corresponding local cost is infinite. Taking into account that by assumption, Us(0) satisfies all the constraints for all the agents at time step t= 0 and using the above result recursively, the statement of this part is proved. Part 2. Taking into account the definitions of Ud i(t−1) and Us i(t) it follows that Ji(xi(t),{Us j(t)}j∈ni)−Ji(xi(t−1),{Ud j(t−1)}j∈ni) (3.17) 2Although we used the index pin the definition of the optimization problems solved to obtain each proposal, in the proof of Theorem 4 we will use the index i. 3.3. STABILITY 85 is equal to Fi(Aixi,N +Pj∈niBij Pp∈mjKjpxp,N )−Fi(xi,N ) +Li(xi,N ,{Pp∈mjKjpxp,0}j∈ni)−Li(xi,0,{Pp∈mjKjpxp,0}j∈ni)(3.18) Taking into account (3.15a), this implies that Mx P i=1 Ji(xi(t),{Us j(t)}j∈ni)−J(t−1) ≤ −Mx P i=1 Li(xi,0,{Pp∈mjKjpxp,0}j∈ni)(3.19) As the proposed algorithm chooses Ud(t) as an input trajectory that improves the cost, it is easy to see that J(t)≤J(t−1) − Mx X i=1 Li(xi,0,{X p∈mj Kjpxp,0}j∈ni) (3.20) Taking into account that recursive feasibility is guaranteed (see the first part of the proof) and the definitions of Fiand Liand following the same lines of though as in [61] or [42], attractiveness and stability can also be proved. This implies that system (3.1) in closed-loop with the proposed distributed MPC controller is asymptotically stable.  The proof of Theorem 4 follows the standard terminal region/terminal constraint approach, see [61]. Stability is inherited from the set of local controllers defined by matrices Kij which by (3.15a) are known to stabilize the system. In fact this result is based on the well known idea “Feasibility implies stability”, see [90]. Remark 7 The stability properties of the proposed scheme rely heavily on the fact that Us satisfies all the constraints of the optimization problem. This implies, that in the start-up and when the controller looses feasibility due to disturbances, Ushas to be calculated either by a centralized supervisor or in a distributed manner by the agents. Remark 8 When applied to a real system in the presence of disturbances and/or possible model errors, if the controller operates close to the state constraints in practice the shifted input trajectory may become unfeasible and in would have to be evaluated again (in a centralized manner or using an appropriate distributed approach). This issue must be taken into account in the implementation procedure of this control strategy. Remark 9 Although in order to implement the proposed controller, the agents don’t need information about the state or the dynamics of the rest of the subsystems, a centralized model of the full system is needed to design the controller so that closed-loop stability is guaranteed. This issue will be shown in the next section. 86 DMPC Based on Game Theory for N Agents 3.4 Controller design procedure The local controllers Kij must satisfy two necessary conditions. First, the centralized system composed by the Mxsubsystems (3.1) in closed-loop with the local controllers must be stable. Second, the jointly invariant sets must exist. The local controllers that depend on each agent; that is, matrices Kji such that i∈mj, must be designed in a way such that (3.15a) holds. To take this condition into account, we will use the following centralized model of the system x(t+ 1) = Ax(t) + Bu(t) (3.21) where x= [xT 1,...,xMx]T, u = [uT 1,...,uMu]T(3.22) and matrices Aand Bare appropriate matrices that depend of the model (3.1) of each subsystem. In addition, stability of each subsystem in closed-loop with its corresponding local feedback must be guaranteed. A sufficient condition to guarantee stability of each of the subsystems is to require that the cost function defined by the matrices Piis a Lyapunov function for the subsystem in closed-loop with its corresponding local feedback. To take into account this condition, we will use the following uncertain model of each of the Mxsubsystems xi(t+ 1) = Aixi(t) + Bivi(t) + Eiwi(t) (3.23) where viis made of the part of the inputs that depend on xiand wiis the part of the inputs that depend on the rest of the states when the local controllers are applied; that is, Bivi(t) = Pj∈niBijKjixi Eiwi(t) = Pj∈niBij Pp∈mj−{i}Kjpxp (3.24) In this case, the objective is to design a controller Ki={Kji}j∈nithat stabilizes the subsystem considering wian unknown disturbance. Matrices Biand Diare appropriate matrices that depend of the model (3.1) of each subsystem. We provide next a set of linear matrix inequalities (LMI) that guarantees that (3.15a) holds and that Kistabilizes the subsystem i. These LMI constraints are obtained following standard procedures, see for example [40, 1, 41]. Theorem 5 Consider system (3.1). If there exist matrices Wi, Yiwith i= 1,...,Mxsuch 3.4. CONTROLLER DESIGN PROCEDURE 87 that the following inequalities hold3      Υ Φ Ψ Ξ ∗Υ 0 0 ∗ ∗ I0 ∗∗∗I     ≥0 (3.25) with Ri=Pj∈niRij,R=diag(R1,...,RMx),K= [K1,...,KMx],KT i= [K1i,...,KMui] and Φ =      W1AT 1+YT 1BT 1YT 1BT 2··· YT 1BT Mx YT 2BT 1W2AT 2+YT 2BT 2··· YT 2BT Mx . . .. . ..... . . YT MxBT 1YT MxBT 2··· WMxAT Mx+YT MxBT Mx      (3.26) Υ =      W10··· 0 ∗W2··· 0 ∗ ∗ .... . . ∗ ∗ ∗ WMx      ,Ξ =       YT 1R1 2 YT 2R1 2 . . . YT MxR1 2       ,Ψ =        W1Q 1 2 10··· 0 ∗W2Q 1 2 2··· 0 ∗ ∗ .... . . ∗ ∗ ∗ WMxQ 1 2 Mx        (3.27) and      WiWiAT i−YT iBT iWiQ 1 2 iYT iR 1 2 i ∗Wi0 0 ∗ ∗ I0 ∗ ∗ ∗ I     ≥0 (3.28) for i= 1,...,Mxthen (3.15a) is satisfied for the matrices Pi=W−1 i,Ki={Kji}j∈ni= YiW−1 iand systems (3.23) are stable in closed-loop with vi=Kixi. Proof We will prove the theorem in two parts. In the first part we will prove that if (3.25) holds, then (3.15a) is satisfied for the matrices Pi=W−1 i,Ki={Kji}j∈ni=YiW−1 i. In the second part, we will prove that if (3.28) holds then system (3.23) is stable in closed-loop with vi=Kixi. Part 1: In this part, we will prove that (3.25) is equivalent to (3.15a). Taking into account the definition of the centralized system (3.21), (3.15a) can be posed as follows (A+BK)TP(A+BK)−P+Q+KTRK ≤0 (3.29) 3The symbol “∗” stands for the symmetric part of a matrix. 88 DMPC Based on Game Theory for N Agents with R= diag(R1,...,RMx) Q= diag(Q1,...,QMx) P= diag(P1,...,PMx) (3.30) with Ri=Pj∈niRij. Taking into account that the Pand P−1are positive defined matrices, if we multiply (3.29) by minus one and apply the Schur’s complement we can recast (3.29) as the following constraint "P−Q−KTRK (A+BK)T (A+BK)P−1#≥0 (3.31) This LMI can be transformed into an equivalent one by pre and post multiplying it by a positive definite matrix "P−10 0I#"P−Q−KTRK (A+BK)T (A+BK)P−1#"P−10 0I#≥0 (3.32) The resulting equivalent matrix inequality is given by "P−1−P−1QP−1−P−1KTRKP −1P−1(A+BK)T (A+BK)P−1P−1#≥0 (3.33) In order to obtain a LMI inequality let Υ = P−1=diag(W1, W2,...,WMx) with Wi=P−1 i for i= 1,2, .., Mxand Y=KΥ = [Y1Y2. . . YMx]. It follows that "Υ−ΥQΥ−YTRY ΥAT+YTBT AΥ + BY Υ#≥0 (3.34) Using the decomposition Q=Q1/2Q1/2and applying Schur’s complement we obtain "Υ−YTRY ΥAT+YTBT AΥ + BY Υ#−"ΥQ1/2 0#IhQ1/2Υ 0 i≥0 (3.35)    Υ−YTRY ΥAT+YTBTΥQ1/2 AΥ + BY Υ 0 Q1/2Υ 0 I  ≥0 (3.36) The same procedure is repeated for R=R1/2R1/2obtaining    Υ ΥAT+YTBTΥQ1/2 AΥ + BY Υ 0 Q1/2Υ 0 I  −   YTR1/2 0 0  IhR1/2Y0 0 i≥0 (3.37)      Υ ΥAT+YTBTΥQ1/2YTR1/2 AΥ + BY Υ 0 0 Q1/2Υ 0 I0 R1/2Y0 0 I     ≥0 (3.38) 3.4. CONTROLLER DESIGN PROCEDURE 89 Defining Φ = ΥAT+YTBT, Ψ = ΥQ1/2and Ξ = YTR1/2, the following LMI constraint is obtained and hence the proof is completed      Υ Φ Ψ Ξ ∗Υ 0 0 ∗ ∗ I0 ∗∗∗I     ≥0 (3.39) Part 2: In this part, we will prove that if (3.28) holds then system (3.23) is stable in closed-loop with vi=Kixi. To this end, we will prove that (3.28) is equivalent to the following constraint (Ai+BiKi)TPi(Ai+BiKi)−Pi+Qi+KT iRiKi≤0 (3.40) which implies that Vi(x) = xT iPixiis a Lyapunov function of the closed-loop system and hence is stable. To prove this part of the theorem the constraint (3.40) is transformed in its equivalent LMI constraint (3.28) following the same procedure used in the first part.  Remark 10 Additional constraints can be added to the design procedure so that there is no need to know the state xiin order to calculate the input uj. This is relevant because in order to evaluate the shifted input trajectory, all the subsystems whose state affects a given input must communicate, so in certain cases, it may be desirable to limit these communications. Once the local controllers and the terminal cost functions are fixed, in order to design a distributed MPC scheme that satisfies the assumptions of Theorem 4 one needs to find sets Ωisuch that (3.15b) to (3.15d) hold. In general this is a difficult problem because each of the sets depends on the others. The size of the terminal region for agent iis determined by the magnitude of the disturbances induced by its neighbor agents and viceversa. A similar class of invariant systems was studied in [81] within the polytopic games framework. We provide next an optimization based procedure to solve this problem. In order to present the algorithm we need the following definitions. Definition 2 Given the following discrete-time linear system subject to bounded additive uncertainties x+=ˆ Ax +ˆ Bu +ˆ Ew (3.41) with w∈ˆ W, subject to constraints in the state and the input x∈ˆ X, u ∈ˆ Uand a linear feedback u=ˆ Kx; a set Ωis said to be a robust positive invariant set for the system if the 90 DMPC Based on Game Theory for N Agents following constraints hold x∈Ω→(ˆ A+ˆ Bˆ K)x+ˆ Ew ∈Ω,∀w∈ˆ W ˆ Kx ∈ˆ U Ω⊆ˆ X (3.42) Given system matrices ˆ A, ˆ B, ˆ E, ˆ Kand the sets ˆ X,ˆ U,ˆ W, there exists several methods to find a set Ω that satisfies these constraints, see for example [39] for a procedure to find the maximal robust positive invariant and [80] for a procedure to find an approximation of the minimal robust positive invariant. We denote Ω( ˆ A, ˆ B, ˆ E, ˆ X,ˆ K, ˆ U,ˆ W) the corresponding maximal robust positive invariant set. In order to obtain sets Ωisuch that Assumption 1 is satisfied, we will use the uncertain model (3.23) of each agent; that is, each agent assumes that the contribution of its neighbors to the inputs that affect its dynamics are an unknown bounded disturbance. The size of the set in which these disturbances are bounded depend on the size of the sets Ωi. This implies that finding these sets is in general a complex problem. In order to decouple the design of each set, each agent ilimits its contribution to each input jby a factor λji ∈(0,1] with Pi∈mjλji ≤1; that is, Kjixi∈λjiUj,∀i, j (3.43) Using the same notation introduced in (3.23), this implies that vi∈ Vi(Λ), wi∈ Wi(Λ) (3.44) with Vi(Λ) = λ1iU1×λ2iU2×...×λMuiUMu Wi(Λ) = (Pp∈m1−{i}λ1p)U1×(Pp∈m2−{i}λ2p)U2×...×(Pp∈mMu−{i}λMupUMu)(3.45) where Λ = {λij}∀i,j is a vector made of all the parameters λij. Note that the maximum contribution of a given agent inside Ωi, is the maximum contribution to the disturbance for the rest of the agents. In order to decouple the computation of the jointly invariant sets Ωi, we use the following result based on finding a robust positive invariant set for each subsystem: Lemma 1 Given constants λji ∈(0,1) with Pi∈mjλji ≤1, if the sets defined as Ωi= Ω(Ai, Bi, Ei,Xi, Ki,Vi(Λ),Wi(Λ)) (3.46) are not empty, they satisfy the constraints (3.15b) to (3.15d). 3.4. CONTROLLER DESIGN PROCEDURE 91 The lemma stems from the definition of the operator Ω. If all the sets exists, then they satisfy the stability constraints. Note that there exists an infinite number of possible values of λji such that these sets exist. In order to chose one, we propose to solve the following optimization problem which maximizes the feasibility region of the distributed MPC controller: max λji f(Ω1×Ω2.. ×ΩMx) Ωi= Ω(Ai, Bi, Ei,Xi, Ki,Vi(Λ),Wi(Λ)) λji ∈(0,1),∀j, i Pi∈mjλji ≤1,∀i (3.47) where function f(·) is a measure of the size of a polyhedron (for example, its Chebyshev radius). Solving problem (3.47) may be difficult in general, however, under certain assumptions it can be posed as a convex problem. In [81] it was proved that the feasibility region of this problem is convex. In the next lemma we prove that the jointly invariant sets Ωiare polyhedra defined by a set of inequalities whose right hand side can be expressed as an affine combination of the constants λij. This implies, that if an appropriate function f(·) is chosen, problem (3.47) can be cast into a convex optimization problem. Lemma 2 If Ai+PjBijKji is stable, then the set Ωi= Ω(Ai, Bi, Ei,Xi, Ki,Vi(Λ),Wi(Λ)) (3.48) is a polyhedron that can be defined as a set of inequalities whose independent term can be expressed as an affine combination of the constants λij, that is, Ωi={xi:Mixi≤bi+X j∈niX p∈mj λjpbij}(3.49) Proof: The calculation of the robust invariant for a linear system is a well known problem and several procedures can be found in the literature, for instance in [37] or [39]. In order to prove the lemma, we will follow the procedure presented in [39]. The main idea is to find the set of states such that the trajectories of the closed-loop system starting from these states fulfill all the state and input constraints for all possible disturbances. This is done in an iterative manner. The set of states that fulfill the constraints after ksteps is determined for increasing values of k. This process is repeated until convergence is obtained, that is, the same set of states is obtained for kand k+ 1. The resulting set is the maximum invariant set. Note that each value of kadds new constraints that the invariant set must fulfill, so the number of restrictions grows with each iteration. 92 DMPC Based on Game Theory for N Agents First of all, we will define the state constraints that the closed-loop system has to satisfy taking into account the constraints in uji, the contribution of the state xito the different inputs uj. By definition of Ωi,uji has to verify uji =Kjixi∈λjiUj, j ∈ni(3.50) Hence, the input constraint condition for the input j(3.4) can be transformed into the following set of inequalities: HujKjixi≤λjibuj, j ∈ni(3.51) Note that as λji ∈(0,1), the set of inequalities is equal or more restrictive than the original input constraints. These inequalities have to be taken into account in the state constraints of the closed-loop system. The new set of state constraints can be written as ˆ Hxixi≤ˆ bxi(3.52) For example, if ni={1,2,...,Mu}, that is, subsystem iis affected by all the inputs, then ˆ Hxi=     Hxi Hu1K1i . . . HuMuKMui      ,ˆ bxi=     bxi λ1ibu1 . . . λMuibuMu      (3.53) Note that the right hand side of the inequalities can be expressed as an affine combination of the constants λij with j∈ni. Let ACLi= (Ai+PjBijKji). If ACLiis stable, then for each value of Λ, the robust invariant set Ωican be determined in a finite number of steps backward k(Λ). Let k∗= maxΛk(Λ). We can compute the robust invariant set for all possible values of Λ as the set of states such that its k-steps ahead predictions satisfy all the constraints for all possible future disturbances; that is, ˆ Hxi(Ak CLixi+ k−1 X g=0 Ag CLiX j∈niX p∈mj−{i} Bijujp)≤ˆ bxi, k = 1,...,k∗(3.54) for all ujp ∈λjpUjwith j∈niand p∈mj−{i}. Taking into account that for all ujp ∈λjpUj with p∈mj−{i}there exists zj∈ Ujsuch that zjX p∈mj−{i} λjp =X p∈mj−{i} ujp (3.55) constraint (3.54) is equivalent to ˆ Hxi(Ak CLixi+ k−1 X g=0 Ag CLiX j∈ni BijzjX p∈mj−{i} λjp)≤ˆ bxi, k = 1,...,k∗(3.56) 3.4. CONTROLLER DESIGN PROCEDURE 93 for all zj∈ Ujwith j∈ni. In order to eliminate the disturbance from the constraints and obtain a deterministic set, let us focus on each of the nrrows of ˆ Hxi(which define a different constraint for each time step ktaken into account in the definition of the invariant set). To denote the r-th row of a matrix Awe will use the [A]r. Using this notation, constraint (3.56) is equivalent to [ˆ HxiAk CLi]rxi≤[ˆ bxi]r−[ˆ Hxi k−1 P g=0 Ag CLiP j∈ni BijzjP p∈mj−{i} λjp)]r, k= 1,...,k∗, r = 1,...,nr (3.57) Let us define σgr ij =maxzj∈Uj([( ˆ HxiAg CLiBij)]rzj) (3.58) Note that σgr ij is a scalar that can be calculated from the system model and constraints. This definition allows us to rewrite constraint (3.57) as: [ˆ HxiAk CLi]rxi≤[ˆ bxi]r− k−1 X g=0 X j∈ni σgr ij (X p∈mj−{i} λjp), k = 1,...,k∗, r = 1,...,nr(3.59) Taking into account that the second term of each of the constraints of (3.59) is an affine combination of the constants {λip}it is possible to find matrix Miand vectors biand bij with j∈nisuch that Ωi={xi:Mixi≤bi+X j∈niX p∈mj λjpbij}(3.60)  Using this result, the problem of finding a matrix Λ that maximizes a measure of the distance can be cast into a convex optimization problem. For instance, let us suppose that our criterium to compare the invariant sets is the radium of a Chebyshev ball inside the invariant region. In this case we are interested in obtaining the maximum xTxas function of Λ that verifies all the constraints, which is a convex problem. 100 DMPC Based on Game Theory for N Agents Procuder Firm 1 Firm 2 FirmMx Demand u0(t) u1(t) u2(t) uMx(t) d0,1 d1,2 d2,3 dMx-1,Mx s1(t) s2(t) sMx(t) Figure 3.8: Linear supply chain. assumed to have no time delays and the material flows have a delay modeled by di,j which corresponds to the time taken by the shipments from node ito node j. The only information shared by the agents is their inputs. In particular, the model of a node needs to keep track of the shipments made by its upstream node. This implies that a model of the form (3.1) can be obtained assigning a different subsystem to each firm iwith the following state vector xi xi(t) =         si(t) ui−1(t−1) ui−1(t−2) . . . ui−1(t−di−1,i)         Note that this model takes into account the different delays by augmenting the state vectors. The inputs are defined by the different shipments variables uj. In this model the first firm, with state x1(t) is the supplier which demands items directly to the factory by u0(t) which is modeled as a pure delay of value d0,1. The last firm is the retailer which must satisfy the external demand uMx(t) which is an external signal not controlled by the system. The control objective is to regulate the stock levels to a desired value ri(t). In addition, the last node of supply chain, the retailer, has to satisfy the external demand. To this end, we consider the following local cost function for each firm Ji= N P k=1 2i(ri k−si k− di−1 P l=k−1 ui−1 k−l)2 where Nis the prediction horizon, the subindex kdenotes the k-steps predicted value of a signal. No terminal cost function is considered. The cost penalizes the deviation of the sum of current stock and the items traveling from the upstream node from the desired reference. Note that if the controller ignores those units that have to arrive in the future, it would ask for more units than needed. The weights of the local cost grow with 2i, that is, the closer a node is to the retailer the more important is. This way of weighting the error is natural since the most important goal of a supply chain is to satisfy the external demand. 3.6. APPLICATION TO A SUPPLY CHAIN PROBLEM 101 The class of linear supply chains considered in this example is defined by the number of firms Mxand the delay parameters di,j. In the following tables we show the results of a set of simulations with three different supply chains of 5, 10 and 20 firms. We denote these scenarios as SUPPLY5, SUPPLY10 and SUPPLY20 respectively. The delay parameters of each supply chain have been randomly chosen with values between 2 and 5. The initial stock si(0) was chosen randomly between 100 and 300. In all these simulations, we assume that the external demand uMx(t) is null and that the objective of the controller is to regulate the stocks to their references. The references ri(t) were supposed to be constant and were chosen randomly between 180 and 280. The simulation times Tfwere set respectively to 50, 100 and 200 sample times. In order to study the effect of the number of proposals considered at each sampling time in the performance of the proposed DMPC scheme, we have applied several controllers which consider a different number of proposals Nprop. Given that the proposals are made randomly, each simulation was repeated 10 times. In addition, a centralized MPC controller has also been applied to the three scenarios as a reference of the performance that can be obtained with a centralized approach. The tables show the cumulated cost mean ¯ Jcum and the corresponding standard deviation σJof each controller. The cumulated cost of each simulation was computed as: Jcum = Tf P t=0 Mx P i=1 2i(ri(t)−si(t))2 In addition the tables show the mean number of sample times ¯ tss that an agent needs to have less than a 5% of error with respect to its reference, as well as the average number of sample times tss that the slowest agent needs to have less than a 5% of error with respect to its reference. These two entries provide additional information on the performance of each controller. In general, the simulations show that increasing the number of proposals Nprop improves the performance of the proposed DMPC scheme. It can be seen that ¯ Jcum and σJare decreasing functions of the parameter Nprop. However, communications can be a scarce resource for some systems and it is important to find a trade-off between the number of communications and the performance. In our example it can be seen that a good trade-off happens when Nprop is around 5Mxcommunications, where Mxis the number of agents. This implies that each agent makes an average of 5 proposals to its neighbors. 102 DMPC Based on Game Theory for N Agents Controller ¯ Jcum σJ¯ tss tss Nprop DMPC 2.36e+6 1.80e+6 25.57 30.70 1 DMPC 9.39e+5 2.99e+6 15.20 18.30 3 DMPC 6.64e+5 3.05e+5 11.25 13.80 5 DMPC 5.53e+5 1.40e+5 9.85 12.20 7 DMPC 5.39e+5 1.69e+5 9.05 11.30 10 DMPC 4.34e+5 8.30e+4 8.25 10.50 15 DMPC 3.88e+5 1.25e+4 7.70 9.50 20 DMPC 3.86e+5 1.24e+4 7.65 9.30 30 MPC 3.71e+5 - 7.50 9.00 - Table 3.1: Simulation results for SUPPLY5 Controller ¯ Jcum σJ¯ tss tss Nprop DMPC 8.50e+7 1.95e+7 58.32 93.50 1 DMPC 4.57e+7 1.20e+7 29.48 46.30 3 DMPC 2.61e+7 3.28e+6 21.78 34.30 5 DMPC 2.62e+7 4.23e+6 20.34 29.50 7 DMPC 2.06e+7 2.98e+6 16.18 24.20 10 DMPC 1.71e+7 1.98e+6 13.81 21.30 15 DMPC 1.70e+7 2.00e+6 13.68 21.50 20 DMPC 1.63e+7 1.25e+6 12.82 20.50 30 DMPC 1.53e+7 9.56e+5 12.91 20.50 50 DMPC 1.52e+7 7.07e+5 12.36 20.10 70 DMPC 1.52e+7 6.51e+5 12.07 20.10 100 MPC 1.45e+7 - 13.00 20.00 - Table 3.2: Simulation results for SUPPLY10 3.7. APPLICATION TO CONTROL OF IRRIGATION CANALS 103 Controller ¯ Jcum σJ¯ tss tss Nprop DMPC 5.84e+11 1.24e+11 162.49 199.80 1 DMPC 2.66e+11 6.02e+10 132.96 187.00 3 DMPC 1.46e+11 2.83e+10 93.34 137.50 5 DMPC 1.30e+11 1.48e+10 76.09 112.90 7 DMPC 9.85e+10 1.73e+10 60.21 88.90 10 DMPC 8.23e+10 1.11e+10 49.48 73.40 15 DMPC 6.19e+10 7.09e+9 42.24 61.70 20 DMPC 5.55e+10 4.52e+9 39.02 56.40 30 DMPC 5.24e+10 3.01e+9 32.38 46.70 50 DMPC 5.07e+10 1.38e+9 31.15 44.30 70 DMPC 5.01e+10 9.34e+8 30.36 42.90 100 DMPC 4.98e+10 8.90e+8 29.91 42.10 125 DMPC 4.97e+10 5.79e+8 29.65 41.80 150 DMPC 4.98e+10 4.03e+8 29.23 41.90 175 DMPC 4.95e+10 6.09e+8 28.31 41.10 200 MPC 3.84e+10 - 19.53 26.00 - Table 3.3: Simulation results for SUPPLY20 3.7 Application to control of irrigation canals In this chapter we apply the proposed controller to a model of a section of the “postrasvase Tajo-Segura” in the south-east of Spain. The ‘postrasvase Tajo-Segura’ is a set of canals which distribute water coming from the Tajo river in the basin of the Segura river. This water is mainly used for irrigation (78%), although a 22% of it is drinking water. The selected section is a Y-shape canals (see figure 3.9), a main canal that splits into two canals with a gate placed at the input of each one of them: Canal de la Pedrera, with a total length of 6,680 kilometres, and Canal de Cartagena, with a lenght of 17,444 kilometres. The total length of the canals is approximately of 24 kilometres. At the end of Canal de Cartagena there is a reservoir with limited capacity. The main elements in the canals are the main gates, which regulate the level of water along the canals, and the off-take gates, where the farmers take water from the canals for irrigation. There are 7 main gates and 17 off-take gates in the section studied. Figure 3.9 and table 3.4 show a scheme of the location of the gates, the off-take gates and the milestones where they are located. 104 DMPC Based on Game Theory for N Agents Code Type P/G Description Km Canal del Campo de Cartagena Starting Campo de Cartagena canal 0,000 CCMICAR-01 Gate G Initial Gate 0,200 MICAR-01 Off-take G Off-take 5 - Fuensanta and Estafeta 1,170 MICAR-02 Off-take G Off-take 5’ - Palacete 2,540 MICAR-03 Off-take P Off-take 6 - Santo Domingo 2,840 CCMICAR-04 Gate Gate Canal Pedrera 4,485 MICAR-04 Off-take P Off-take 7 - Campo Salinas 5,970 MICAR-05 Off-take G Off-take 8 - San Miguel 6,550 MICAR-06 Off-take G Off-take 9 - Las Caadas 8,050 MICAR-07 Off-take G Off-take 10 - San Miguel 9,390 MICAR-08 Off-take P Off-take 11 - Campo Salinas 9,590 CCMICAR-05 Gate Gate Tunel San Miguel 10,480 MICAR-09 Off-take G Off-take 12 - San Miguel 12,630 MICAR-10 Off-take P Off-take 13 - Campo Salinas 12,780 CCMICAR-06 Gate Gate La Rambla La Fayona (start) 14,433 CCMICAR-07 Gate Gate La Rambla La Fayona (end) 14,579 MICAR-11 Off-take P Off take 14 - Villamartin 16,540 CCMICAR-08 Gate Gate Caada La Estacada 17,444 Canal de la Pedrera CCMIPED-01 Gate Starting La Pedrera canal 0,000 MIPED-01 Off-take G Off-take 1P - S. Domingo 0,770 MIPED-02 Off-take G Off-take 2P - S. Domingo y Mengoloma 3,740 MIPED-03 Off-take P Off-take 3P - S. Domingo 4,260 MIPED-04 Off-take G Off-take Riegos Levante 1 5,260 MIPED-05 Off-take G Off-take 4P - Santo Domingo 6,440 MIPED-06 Off-take G Off-take Riegos Levante 2 y 3 6,680 Table 3.4: Data of irrigation canal Cartagena-La Pedrera 3.7. APPLICATION TO CONTROL OF IRRIGATION CANALS 105 Figure 3.9: Scheme of the canal. The dynamics of water flowing in irrigation open canals can be obtained by applying the Saint Venant equations [16, 17], which are nonlinear partial differential equations. The irrigation canals considered are divided into several sections separated by gates; the controlled variables are the downstream water levels, hi(t)∈R+(m) and the manipulated variables are the check point to gates, ui(t)∈R+(m). Each canal reach has an inflow from an upstream canal reach, Qin,i ∈R+(m3/s), and an outflow to a downstream canal reach, Qo,i ∈R+(m3/s). Also, other flows are considered as perturbation variables. In particular, qin,i ∈R+(m3/s) models the flows due to rainfall, failures in upstream gate and other unknown disturbances and qo,i ∈R+(m3/s) models the known offtake outflows from farmers, considered as measurable perturbations. The discrete model that has been considered using the previous variables is: Ai(hi(k+ 1) −hi(k)) = Td(Qin,i(k−td) + qin,i(k)−Qo,i(k)−qo,i(k)) (3.68) where Td(s) is the sampling time, Aiis the surface of the reach and tdthe delay of the input Qin (the level is measured downstream). 106 DMPC Based on Game Theory for N Agents The discharge through a submerged flow gate Qo(t) is usually determined as[16]: Qo(t) = CdLp2gu(t)qhup(t)−hdn(t),(3.69) where Cdis the gate discharge coefficient, Lis the gate width, u(t) the gate opening and hup(t), hdn(t) are the upstream and downstream water levels, respectively. The main target is to manage the water in the canals in order to guarantee the flows requested by users. Another objective to be considered is the minimization of the leaks and evaporation (which can be obtained as function of the levels) and also to minimize maintenance costs (the maintenance of concrete blocks and junctions is better if they are submerged, so high levels are preferred for that purpose). For these purposes, it is necessary to maintain the level of the canal over the off-take gate when flow is requested. The controlled variables are the upstream levels beside the gates, which have to satisfy maximum and minimum level constraints. The minimum level is determined by the demand of irrigation of surrounded lands; it must be guaranteed that off-take points are submerged. The maximum level is determined to provide regulation capabilities to the system in order to avoid floods. The manipulated variables are the flow at the head of the canal and the position of the gates, which are are also constrained. The flow at the head is limited by the total amount of available water and the gates have maximum and minimum openings. The actual reference in levels are sent to the DMPC in the low level. For this controller, the sample time has been considered 1 minute and the horizon, N, has been set to 5. The cost function has been designed to include the water traveling from the upstream gates, that is, the states which correspond to the traveling water are also weighted. It is necessary to consider the water coming from the upstream gate as a part of the level under consideration to be able to compensate the effect of the delay in the controller. The weights of the local costs in the canals grow with 2i, that is, the farther a node is from the beginning , the more important is. This way of weighting the error facilitates a faster flow of water towards the last canals. Finally, the matrix that weights the control effort Rihas been set to zero for simplicity. In figure 3.10 it can be seen a simulation in which all the reaches begin with a water level of 3 meters. At sampling time k= 0, the reference is set to 3.15mto all the reaches. After a whole day (k= 1440) there is another change of reference for all the reaches to 3.6m. These changes are originated in the upper control level as a function of the risk mitigation policy. It can be seen how the level in the reaches follow the reference even in the presence of disturbances (farmers take water, rains...). It is important to remark that these results have been obtained with an average number of 5 communications per agent and sampling time. 3.8. CONCLUSIONS 107 In other words, each agent makes up to 5 proposals in a minute in order to get a cooperative solution with the rest of the agents. 0 500 1000 1500 2000 2500 3000 2.9 3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 levels k s1 s2 s3 s4 s5 s6 s7 Figure 3.10: Water level evolution. 3.8 Conclusions In this chapter we have presented a novel distributed MPC algorithm based on negotiation for a class distributed linear systems coupled through the inputs. We assume that each agent has access only to the model and the state of one of the subsystems and that the agents must negotiate in order to reach a cooperative solution. The proposed algorithm has low communication and computational burdens and provides a feasible solution to the centralized problem. In addition, we provide sufficient conditions that guarantee practical stability of the closed-loop system as well as an optimization based procedure to design the controller so that these conditions are satisfied. 108 DMPC Based on Game Theory for N Agents Chapter 4 Distributed Receding Horizon Kalman Filter The most common approach to estimate the state of stochastic systems is the Kalman filter [36], developed in 1960 and named after his discoverer. The Kalman filter is the optimal state estimator for unconstrained linear systems subject to gaussian state and output noise. It is not possible to apply directly the centralized Kalman filter to a multiagent problem unless there is a node in the network that receives all the information. For example in [87], it can be seen how a central agent gathers the information from the moving devices and then distributes the position estimation back to them. An alternative is to calculate a decentralized version of the Kalman filter that takes into account the communications restrictions [31, 87]. In this thesis we follow a different approach to solve the estimation problem in a distributed manner. First, the Kalman filter is posed as a dynamic programming problem [18]. Then, the resulting optimization problem is distributed among the agents by means of dual decomposition. This idea has been successfully applied to distributed control in [82] and [28]. Given that the observation problem is the dual of the control problem, it is natural to apply and enhance the techniques presented in these papers to deal with the state estimation problem. Moreover, some of the results that will be shown in this chapter have direct application to distributed control based on dual decomposition. In this context, the application of state estimation schemes to problems in which the state represents the position of an object is very attractive [4], [31]. The localization of moving entities, such as robots or people, is important for many applications. Military applications in which the goal is to track a target that moves in a distributed sensor environment are typical examples. Other examples in which these techniques play an important role would 109 116 Distributed Receding Horizon Kalman Filter Then, the extended local cost function can be posed as VN,p i(ˆ Xi(τ−N:τ), P(τ−N:τ)) = τ P k=τ−N [li(ˆxi(k), di(k)) + pT i(k)di(k) + ˆxi(k)TPi6=jCT jipj(k)] The local stage cost can also be extended to include the terms due to the lagrangian prices lp i(ˆxi(k), di(k), P(k)) = li(ˆxi(k), di(k)) + pT i(k)di(k) + ˆxi(k)TPi6=jCT jipj(k), which allows to write the extended local cost function as: VN,p i(τ) = τ P k=τ−N lp i(ˆxi(k), di(k), P(k)) (4.15) Remark 15 After the introduction of dual variables, and assuming that the prices of the neighbors are given, it is possible to interpret the distributed optimization procedure in economic terms. Each agent behavior can be represented as a two player game. The first player objective is to minimize the price-extended stage cost τ X k=τ−N lp i(ˆxi(k), di(k), P(k)), which is composed of three elements that are interpretable as lp i(ˆxi(k), di(k), P(k)) = li(ˆxi(k), di(k)) | {z } local cost + neighbor help cost z }| { pT i(k)di(k) + ˆxi(k)TX i6=j CT jipj(k) | {z } incomes due to required help . The second player chooses the prices pi(τ−N),...,pi(τ)to maximize pT i(k)(di(k) + X i6=j Cij ˆxj(k)). This game is repeated iteratively. First, an estimate is calculated according to the given prices. Then, the prices are updated and the cycle starts again. As a result of the repeated interaction 4.2. DUAL DECOMPOSITION 117 of both players in each node the prices evolve until a maximum is reached. The consequence of this standard Lagrangian optimization procedure is that the minimum for the cost function (4.12) is attained and the constraints are satisfied when the price gradient is zero. The algorithm that is followed by the agents in the system can be summarized as: •Step 1: Each agent iestimates his own current state trajectory {ˆxi(τ−N),ˆxi(τ−N+ 1), .., ˆxi(τ)}solving the optimization problem given in (4.15) for a set of given prices pii= 0,...,J. •Step 2: Once the state trajectory has been calculated then the prices of agent iare updated by a gradient step as follows. pk+1 i(τ) = pk i(τ) + γk i[di(τ) + X i6=j Cij ˆxj(k)] (4.16) Convergence of such gradient algorithms has been proved under different type of assumptions on the step size sequence γk i. See for example [91]. Note that in order to update the prices the agents must communicate. •Step 3: If the precision obtained with the estimation is enough then there is no need to continue iterating. In the next section precise conditions are given. If enough precision is not attained and the number of iterations Kexceeds a given threshold maxiter, then the algorithm also stops. In other case then the process is repeated from step 1 for K=K+ 1. 4.2.1 Coordination alternatives for the price update It can be seen that the calculation of the estimate ˆxi(t) for t=τ−N,...,τ is completely decentralized once that prices are given. Therefore it is mandatory for an agent to keep the track of its neighbor prices. Nevertheless, in order to update the prices, coordination among the agents is necessary. The agents send their estimates ˆxi(τ) to their neighbors so that equation (4.16) can be applied. For some systems it could be desirable not to share the state information with their neighbors. To avoid the exchange of the state estimates we propose two alternatives: •Decentralized approach: The need for the shared information comes from term Pi6=jCij ˆxj(k) in equation equation (4.16). According to the dynamics of the subsystems P j6=i Cijxi(τ) = 118 Distributed Receding Horizon Kalman Filter yi(τ)−Ciixi(τ)−vi(τ), and thus it could be approximated by P j6=i Cij ˆxi(τ)≈yi(τ)− Cii ˆxi(τ). •Market approach: This alternative consists on changing the way in which prices are updated. To understand better this approach it is convenient to use the behavior model that represents each agent as a two player game. Then, it is possible to think on the centralized problem as a game with 2Jplayers. The objective of the first player in each node is to minimize his own cost according to the given prices. However, the second player in each node bargains with the the rest of the second players to maximize (4.15) with respect to the prices. The second players can be seen as market makers that fix the prices of the help services that the agents provide each other according to the offer and demand of such services. To do so, a gradient optimization of the cost function (4.15) is implemented. Each update is based on the addition of contributions of the different agents. The contribution of agent iis ∇pk i(τ) =              ˆxi(k)TCT 1ip1(k) . . . ˆxi(k)TCT i−1,ipi−1(k) di(k) ˆxi(k)TCT i+1,ipi+1(k) . . . ˆxi(k)TCT JipJ(k)              Theorem 6 The price update mechanism defined in the market approach provides the same results than the one presented in equation (4.16). Proof: It is straight forward to check that both methods provide the same centralized price vector. It is enough to sum the contribution ∇pk i(τ) for all i pk+1(τ) = pk(τ) + γX i∇pk i(τ). Then it can be seen that the price for agent iis just pk+1 i(τ) = pk i(τ) + γk i[di(τ) + X i6=j Cij ˆxj(k)]  If we move back to the agents and forget the game theory interpretation, it can be seen that under the market approach agents update their prices and also the prices of their neighbors and therefore there is no need to exchange the state estimate. All 4.3. EXAMPLES 119 the public information needed are the prices and their updates. The estimation of the agents through the different iterations bring increments or decrements in the prices until equilibrium prices are reached. However, there is a price to pay in terms on the amount of model information that agents have. With this price mechanism it is needed that agent ihas knowledge of the terms Cji. In other words, agents have knowledge of the collateral effects they induce in their neighbors. Remark 16 From an economic point of view, the situation can be interpreted as a market of help services. The price pi(τ)is the unit price that agent ihas to pay to his neighbors for them to change their current contribution to his output. The fact that neighbors of an agent ichange their estimates affects to his price in such a way that it reflects how costy is for his neighborhood to help him after the estimate update. On the other hand helping his neighbors is rewarded in (4.15). Taking all of this into account, agents are both service offerers and demanders. All of them behave selfishly according to the prices fixed by the market, that is, the distributed price mechanism proposed in this approach. Remark 17 In welfare economics, under certain assumptions such as the absence of externalities in transactions, it is proved that market prices guarantee that, despite of agents selfish behavior, a Pareto optimum is achieved [97]. In the optimization problem that we have, unfortunately we have to deal with the presence of externalities, taking this term in a wide sense. That is, decisions taken by agent ialso affect other agents. In order to overcome this problem and still reach a Pareto optimum while keeping selfish, i.e. decentralized, behavior in the agents, some modifications have to be introduced in the market: first, all the agents behave as price takers as they were in a competitive market when they really have power to modify the prices and, second, prices are updated globally according to the proposed mechanism. 4.3 Examples The problem of estimating the position of a moving object can be faced using different approaches. For outdoor applications in which the precision requirements are low GPS estimation is the most used choice. Radar measurements help to improve the quality of the estimation. When it comes to indoor applications the problem of localization is normally solved by means of a sensor network. In cases in which low precision is needed some it may be enough with a network of presence detectors. If more precision is required then more sophisticated techniques have to be used. In the case that the application is executed in a very controlled scenario, it is possible to use cameras to estimate the position. Infrared or ultrasonic sensors also provide a greater accuracy than the presence detectors. In the last 120 Distributed Receding Horizon Kalman Filter years the use of the link quality between wireless transceivers has been used too for this kind of applications [49]. 4.3.1 Application to mobile robot localization This subsection is based on the simulation scenario proposed in [31]. Let us consider a system consisting a set of µ={1, .., M}reference nodes or beacons and a set η={1, .., J}of mobile devices. In this example we will consider M= 6 beacons and J= 8 mobile devices, which are located in the positions depicted in figure 4.1. The goal is to estimate the position of the moving devices. If the sample time is assumed to be low enough, it is possible to simplify the dynamics considering that the devices move at every sample time a bit with respect their position. The equations for each device are: xi(τ+ 1) = xi(τ) + ∆xi(τ)∀i∈η={1, .., J} with xi(0) = x0 i. The beacon position is fixed so that xi(K+ 1) = xi(0) ∀i∈µ={1, .., M}. The distance between the nodes and the mobile devices can be calculated using d2 ij = (xi−xj)T(xi−xj)∀i, j ∈η, µ. The distance can be linearized around the steady state positions xiusing a first order Taylor approximation, which leads to d2 ij =d2 ij + 2(xi−xj)T(xi−xj) + 2(xi−xj)T(xi−xj) with d2 ij =d2 ij(xi, xj). Now, system variables can be introduced for all the mobile devices such that: xi(τ) = xi(τ)−xi∀i∈η yji(τ) = d2 ij −d2 ij ∀i∈η, ∀j∈η, µ Cji = 2(xi−xj)∀i∈η, ∀j∈η, µ. 4.3. EXAMPLES 121 Figure 4.1: Initial situation of the devices. Each moving device’s output provides information about the distance with respect the other moving devices and the beacons. If white gaussian additive noise is assumed in the state and output then each device can be modeled according to equation (4.1). In order to make the situation more realistic it is assumed that only devices and beacons within a range can communicate. Thus a communication radius ρis defined. In general two devices iand jcan communicate if dij < ρ. A communication graph can be defined to reflect which devices can communicate at each sample time. The communication graph at initial time is given by the following matrices: Aη 0=               1 1 0 0 1 1 0 0 1 1 1 0 0 1 1 0 0 1 1 0 0 1 1 1 0 0 0 1 1 1 1 0 1 0 0 1 1 1 0 0 1 1 1 1 1 1 1 0 0 1 1 1 0 1 1 1 0 0 1 0 0 0 1 1               , Aµ 0=               1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 1               where Aη 0(i, j) = 1 if the mobile device iis able to communicate with the mobile device jand Aµ 0(i, j) = 1 if the mobile robot iis able to communicate with the beacon j. 122 Distributed Receding Horizon Kalman Filter The simulations have been done considering a dynamic graph, that is, a situation where the movement of the devices is big enough to guarantee that the communication graph changes with the relinearization of the system. At each time in which the system is relinearized it is necessary not only to update the equations but the information about the last samples that is kept in the agents. Let us assume that in time τthere is a change of linearization point of the system. Then, xi(t) = ¯xi(t) + xi(t)−¯xi(τ)∀t∈[τ−N, τ] yi(t) = Ci(τ)ˆxi(t)∀t∈[τ−N, τ] This change of coordinates in the state estimates and the outputs allow to compute the distributed problem without suffering estimation disturbances after the change of linearization point. The system has been simulated for 40 time samples with a state noise stronger than the original one. The first 10 samples a centralized Kalman filter is working and the second 20 the distributed strategy. In t= 20 and t= 30 the system is relinearized. The window size used for the estimation was 4. In blue it is depicted the real trajectory and in red the estimation. The results for the estimation of the position of the mobile devices can be seen in figure 4.2. The overall picture is shown in figure 4.3. The quality of the estimation depends on several parameters. For example, the more iterations are made the better the estimation gets. In this figure it can be seen that the estimation is very precise for most agents. 4.3.2 Application to traffic and speed control Imagine a scenario where there are reference nodes in roads and streets, and cars are equipped with wireless transceivers. Each car could estimate its position using the link quality indicator from the packets in the reference nodes and other cars as an indicator of the distance. In such set up it is possible to imagine many useful applications that benefit from the algorithm presented in this chapter. Some ideas to take advantage of the position estimation would be to use the information as a way to monitor and control the traffic, or to check if the cars are moving without exceeding the speed limits, or, for example, for cars with free parking places around there to communicate it to the network, reducing search time for neighbors. The price to pay to fully enjoy this possibilities would be a loss of privacy because the network would be aware of the position of all the cars in this hypothetic scenario. An alternative could be to hide the identity of the cars so that applications such as traffic control could be maintained. 4.4. CONCLUSIONS 123 0 10 20 30 40 1 2 3 t x11, x12 0 10 20 30 40 0 5 t x21, x22 0 10 20 30 40 0 5 10 t x31, x32 0 10 20 30 40 0 5 t x41, x42 0 10 20 30 40 0 2 4 t x51, x52 0 10 20 30 40 0 2 4 t x61, x62 0 10 20 30 40 3 3.5 4 t x71, x72 0 10 20 30 40 2 4 6 t x81, x82 Figure 4.2: Robots’ state evolution with noise model mismatch. In figures 4.5 and 4.4 it can be seen a very simple simulation of what would happen if the mobile robots of the previous subsection scenario would behave as cars going left or right as in a normal road. It is assumed that once a robot has been located it is only necessary to estimate its xcomponent. For this simulation the window size was fixed to 3 and the communication radius to 5. The average number of iterations to reach a 95% of accuracy was 6.27. 4.4 Conclusions A distributed version of the Kalman filter based on dynamic programming has been developed. The use of dual decomposition allowed the problem distribution. In the simulations presented promising results of the future applications of these techniques are shown. The different coordination alternatives for the price update that have been presented are also remarkable. In particular, the market approach allows to use dual decomposition without revealing the state of the subsystems. This feature may be interesting in control applications in which dual decomposition is used as well. It will be important for future work some kind of suboptimality bounds to determine the precision obtained in the estimation after a number of iterations. Practical experiments will be developed too to see how the distributed estimation works in real application. 124 Distributed Receding Horizon Kalman Filter 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 1 1.5 2 2.5 3 3.5 4 4.5 x1 x2 x est. x Figure 4.3: Robots’ trajectories with noise model mismatch. 0 10 20 30 40 0 2 4 t x11, x12 0 10 20 30 40 0 5 10 t x21, x22 0 10 20 30 40 0 5 10 t x31, x32 0 10 20 30 40 3 4 5 t x41, x42 0 10 20 30 40 0 2 4 t x51, x52 0 10 20 30 40 0 5 10 t x61, x62 0 10 20 30 40 2 4 6 t x71, x72 0 10 20 30 40 0 5 10 t x81, x82 Figure 4.4: Robots’ state evolution with noise model mismatch in straight movement. 4.4. CONCLUSIONS 125 1 2 3 4 5 6 7 8 0.5 1 1.5 2 2.5 3 3.5 4 x1 x2 x est. x Figure 4.5: Robots’ straight trajectories with noise model mismatch. 132 Applications of Cooperative Game Theory the cost that would suppose for the agents in Sachieve the goal without the cooperation of the agents out of Sfor a given communication architecture L. 5.2 Distributed control problem formulation In this chapter an application of coalitional games to distributed control is proposed. The main objective of the proposed approach is to provide an a priori analysis of the best possible use of the communication network at each sampling time assuming that the use of each link has a cost. In addition, a qualitative interpretation of the results is given so that it can be determined which agents need most to communicate and which links are more important for a given distributed scheme. We consider a linear system divided in i= 1,...,N subsystems defined by the following model xi(t+ 1) = Aiixi(t) + Biiui(t) + di(t), di(t) = Pj6=iAijxj(t) + Pj6=iBijuj(t),(5.3) where xi∈Rqiand ui∈Rriwith i= 1,...,N1are the states and inputs of each subsystem respectively. The variable di(t) is the influence of the neighbors’ states and inputs in the update of xi. Each agent ihas access only to its state xiand decides at each sample time the value of its corresponding input ui. We assume that there exists a network Lwhich allows the agents to exchange information. The type of information exchanged depends on the distributed control algorithm that is being used. Any two agents that are not indirectly connected by the network will not be able to exchange any type of information. The control objective is to regulate the state of all the subsystems to the origin while minimizing a cost which depends on the state and input trajectories and on the communications. This cost will be used to define the characteristic function of a cooperative game. The stage cost of each agent is defined as follows ℓi(t) = xT i(t)Qixi(t) + uT i(t)Riui(t). The objective is to minimize the total cumulated cost taking into account the communication costs defined in the previous section. In some applications the stage cost can be interpreted in economic terms. 1In this chapter, the letter Nstands for the number of players in the game, not the horizon used in the MPC schemes. 5.2. DISTRIBUTED CONTROL PROBLEM FORMULATION 133 The communication costs depend on the number of links that are being used. Note that given an initial network Lnot all the links have to be used all the time. The term network mode will be used to denote each of the possible subset of links Ain L. The change of performance for the control system will be analyzed for the case in which some of the links were not present. The analysis of the cost trade-off between control performance and communicational costs will determine which is the best network mode Aat state x. Remark: The value of the cost of the use of a communication link during a sample time has to be determined ad hoc. For example, in a wireless network this value can be function of the inverse of the remaining battery of the nodes that support the link. In case that the characteristic function of the games can be interpreted in economic terms then it could be calculated the actual cost of transmissions through the link. In general, a simple way to provide a value is to assign is to average the cost impact between enabling or disabling the link in the system or just to set a bound on cost improvement for the link to be enabled. Remark: In systems where the number of agents is too big the calculations can be simplified assuming that w(USV) = w(U)+w(V), that is, the value of the coalition of the players in the set USVis equal to the sum values of Uand V. This approximation is much better when Uand Vare in different communication components. Remark: Note that initially the characteristic function for a game defined from a distributed control problem should be subadditive, that is w(S) + w(T)≤w(SST), that is, the control performance gets better as there are more agents involved. However there is a hidden implication in the last property: communication is costless. Theoretically the last statement can be a good starting point to develop some results, but in practice this results to be a fallacy. 5.2.1 Distributed control algorithm In this section we present a distributed control scheme that at each sampling time, implements a certain communication strategy defined by a network mode A. The communication strategy Ais chosen every Dsample times. To this end, the agents must broadcast their state and take a decision about the communication strategy that will be used in the next Dtime steps. This leads to a double sample rate control system. As a result of this policy, the agents are separated into separated groups Cthat are able to communicate defined as communication components. We will denote by N/A the set of all communication components in which the set Nis partitioned. Note that S ∀i∈N/A Ci=Nand CiTCj= 0 for all i6=j. We assume that for each network mode A, a different controller that stabilizes the whole 134 Applications of Cooperative Game Theory system which takes into account which agents can communicate is defined. The details about the calculation of the controller will be presented later in this section. In particular, we assume that each communication component C⊆N/A implements a linear controller uC=KA CxC, where uC∈RPi∈Criis the input of a given communication component defined as uC= {ui}i∈C,xC∈RPi∈Cqiis the the state of a given communication component defined as xC= {xi}i∈Cand KA Cis the matrix which defines the controller implemented by the communication component for the network mode A. From the set of the matrices KA C, the following centralized linear controller for the overall system, characterized by the absence of communication for agents in different communication components can be obtained u=KAx where u∈RPi∈Nriis the input of the centralized system defined as u={ui}i∈N,x∈RPi∈Nqi is the the state of the centralized system defined as x={xi}i∈Nand KAis the matrix which defines the centralized controller implemented for the network mode A. Note that matrix KAtakes into account the communications constraints in A. For example, if the i-th element of uand the j-th element of xbelong to different communication components, then KA(i, j) = 0; that is, the i-th input does not depend on the j-th state. For the particular case in which each communication component is composed by systems with consecutive numeration in the set N, then KA=diag(KA C1, KA C2,...). Note that all the matrices KA Chave to be designed so that KAguarantees closed-loop stability for the centralized system. If a given mode is not able to stabilize the system, then it is not taken into account. In order to decide which communication strategy must be implemented, we assume that there exists a quadratic function that satisfies xTPAx≥X j∈NX k=0,..,∞ ℓj(k) (5.4) that is, PAis a weight matrix that provides an upper bound of the cost to infinity of the centralized system in closed-loop with the controller u(k) = KAx(k) starting from the initial state x(0) = x. This quadratic function will be used so decide the optimal communication mode as well as to define the link problem. Summing up, the proposed distributed control scheme is implemented as follows: 1. If the sample time is a multiple of D, all the agents broadcast their state and calculate which is the network mode Athat minimizes the function rH(A, x). Otherwise, each agent sends his state only to those agents that belong to his communication component. 5.2. DISTRIBUTED CONTROL PROBLEM FORMULATION 135 2. Each agent uses the state information received in order to update its control action using its corresponding communication component feedback matrix KC A. Globally, this implies that the linear controller u=KAxis applied. Remark: It is possible to establish more sophisticated assumptions about the system. However, it has been preferred to simplify the system as much as possible in order to focus on the three mentioned problems. The definitions done are enough to identify both the agent and link games and to solve the decision problem of what agents should communicate at each sample time. Remark: The ideas presented in this chapter can also be applied to different distributed control strategies (such as distributed MPC schemes) or to more complex systems (such as non-linear systems). In that case, an appropriate definition of the utility function has to be provided. We propose to use a bound of the cost-to-go of the different modes, but other approaches are also possible. 5.2.2 Network modes In order to analyze which agents should communicate, the link game associated to the costextended communication situation Hof the current state xis studied. The characteristic function that assigns a value to each communication mode Ais based on the upper bound of the cost-to-go of all the communication components in the network and the corresponding communication costs and is defined as rH(A, x) = xTPAx+c|L(A)|,∀A⊆L. The best possible communication mode provides the winner network configuration choice. This mode is obtained for a given state xby minimizing rH(A, x) over A. The function rH(A, x) will be also used to define the link game. We introduce next the concept of dominance between modes. The mode Adominates the mode Bif rH(A, x)< rH(B, x). In general, we say that mode iis dominant if the last inequality holds ∀B6=A⊆L. Therefore, the set of points CRAfor which a network mode A is dominant is characterized by the states such that Aprovides the best (lower) cost, that is, CRA={x∈R|rH(A, x)≤rH(B, x),∀B⊆L}. The union of all these regions covers the whole state space. The distributed controller will switch between the different modes as the state moves from one region to another. In order to check online which is the optimal network mode for the particular distributed control 136 Applications of Cooperative Game Theory scheme presented, it is sufficient to evaluate the quadratic functions which define the cost of each mode. Note that in general, the regions associated to the network modes are non convex. For more complex cases it may not be possible to calculate a explicit characteristic function rH(A, x) for each network mode and the application of techniques might be necessary to find off-line the regions of each of the modes in order to check the optimal mode. Suboptimal approaches based on exhaustive simulation of all the possible network modes can be used to find these regions of dominance. For the distributed control scheme considered, the boundaries of the regions are defined by quadratics which depend on the different weight matrices PAand number of links |L(A)|which define each mode. In particular, it is possible to calculate explicit frontiers. To do so, let Ji=xTPix+ciand Jj=xTPjx+cjbe the costs for network modes iand jrespectively. Mode iis chosen if Ji< Jjand viceversa. Thus, the frontier between the regions of dominance of networks modes iand jis given by the following equation: Ji=Jj xTPix+ci=xTPjx+cj xT(Pi−Pj)x+ci−cj= 0 (5.5) The shape of the bound will depend on the matrix Pi−Pj, which in general does not have to be positive definite. 5.2.3 Link analysis The link game is constructed by a set of players composed by the links of the network L and a characteristic function that assigns to each communication mode Aan given utility. In this case, we propose to use rH(A, x) to define the link game. Once the link game is constructed, a qualitative and quantitative analysis of the relevance of the links in the game may be obtained from the corresponding Shapley value γ(L, rH, x). Each component of this vector represents the cost of a given link for the system. In other words, the lower value a link has, the higher utility it has for the system. 5.2.4 Agent analysis A qualitative and quantitative analysis of the relevance of the agents of the system may be obtained from the Shapley value γ(N, wH, x) of the corresponding agent game. To build such game it is necessary to define the characteristic function that assigns a value to each coalition Sof agents for a given network L. To this end, it is not possible to use the controllers defined for each communication component for a given network mode Abecause those controllers take into account the particular communication constraints of A. For this reason, for each communication component of L, we define a controller KCand a weight 5.2. DISTRIBUTED CONTROL PROBLEM FORMULATION 137 matrix PCsuch that xT CPCxC≥X j∈CX k=0,..,∞ ℓj(k) (5.6) where xcis obtained from the current state xand is composed of the states of all the subsystems that belong to C; PCis a weight matrix that provides an upper bound for the cost to infinity of the systems that belong to Cin closed-loop with the controller uC(k) = KCxC(k) starting from the initial state xC(0) = xCassuming that all the inputs and states that belong to agents outside the coalition are zero. Then, the characteristic function of the agent game that defines the utility of a coalition Sis defined as wH(S, x) = X C∈S/L xT CPCxC+c|L(S)|,∀S⊆N. The Shapley value of the game (N, vH, x) provides concise information about the relevance of all the agents in the game. The lower the value the is, the more relevant role the agent has in the game. It is important to stand out that the Shapley value of the agent game, as it is defined, does not have any physical meaning. Remark: This assumption allows only to calculate a simple approximation of the Shapley value of the agent game. More conservative choices could have been made, for example the agents outside of the coalition could have been considered as disturbances and then a minmax approach used. Nevertheless, we must not forget that the goal of all the agents is to regulate the system to the origin. Moreover, with the design method presented in the next subsection it is possible to calculate feedback gains such that the coalitions Cand N−C are able to stabilize the overall system. For this reason, we consider the approximation made appropriate. 5.2.5 Coalitional game design method In this section we present a method to design for a given system (5.3) all the matrices that define the controllers and the weights of the upper bound functions for each of the networks modes (KAand PAfor all A⊆L) as well as each of the possible communication components of the agent problem (KCand PCfor all C⊆N/L). To this end, we present two different theorems. First, we deal with the problem of finding the matrices that define the controllers and the weights of the upper bound functions for both each of the possible communication components of the agent problem (KCand PCfor all C⊆N/L). In this case, KCmust stabilize the states of Cassuming that the outputs and inputs that do not belong to that communication 138 Applications of Cooperative Game Theory component are zero. In addition, PCmust guarantee that (5.16) holds. The following theorem presents an LMI constraint that can be used to solve this design problem. Theorem 7 Let C∈S/L be a set of independent communication components for a given communication situation Hwhose dynamics are given by AC={Aij},∀i, j ∈Cand BC= {Bij},∀i, j ∈Cand its stage cost defined by QC=diag(Qi)and RC=diag(Ri),∀i∈C. If there exist matrices WCand YCsuch that the following constraint is satisfied      WCWCAT C+YT CBT CWCQ1/2 CYT CR1/2 C ACW+BCY WC0 0 Q1/2 CWC0I0 R1/2 CYC0 0 I      >0 (5.7) then matrices PC=W−1 Cand KC=YCW−1 Csatisfy (5.16) and stabilize the states of C assuming that the outputs that do not belong to that communication component are zero. Proof: Applying iteratively backwards Schur’s complement to equation (5.7) and taking into account the proposed variable change, it can be seen that if (5.7) holds then the following constraint can be obtained (AC+BCKC)TPC(AC+BCKC)−PC+QC+KT CRCKC≤0 (5.8) This constraint guarantees that the system defined by matrices AC, BCis stable in closedloop with a the linear controller defined by KC. In addition, pre and post multiplying (5.8) by xC(k) we obtain the following inequality xC(k+ 1)TPCxC(k+ 1) −xC(k)PCxC(k) + X j∈C Lj(k)≤0 Summing the previous inequality from k= 0 to k=∞and assuming that limk→∞ xC(k) = 0 (recall that the closed-loop system is stable) we obtain that xC(0)PCxC(0) ≥ ∞ X k=0 X j∈C Lj(k).  5.2. DISTRIBUTED CONTROL PROBLEM FORMULATION 139 In the agent game the approximation of the cost of the communication component Cis based on the cost to infinity given by the xT CPCxC. This upper bound is calculated assuming that the rest of the agents states and inputs are zero. The agent game only provides grounds for distributing the benefits or costs between the agents during the game and the resulting feedback gains are never used to control the system, so it is not necessary to impose centralized stability as a requirement. The link game is based on a different point of view of the communication situation. In this case the grand coalition controls the system taking into account all the possible network configurations A⊆L. Each network configuration Adivides the system in a set of communication components C∈N/A. In this case the stability of the centralized system has to be guaranteed because these linear feedback will be applied to control the system.. Following the same approach as in the agent game, the following theorem is presented. Theorem 8 Let A∈Lbe a set of active links for a given communication situation H. The dynamics of the whole system are given by AN={Aij},∀i, j ∈Nand BN={Bij},∀i, j ∈N and its stage cost defined by QN=diag(Qi)and RN=diag(Ri),∀i∈N. If there exist matrices WN={Wij},∀i, j ∈N, where Wi,j ∈Rqi×qj, and YN={Yij},∀i, j ∈N, where Yi,j ∈Rri×qj, such that the following constraints are satisfied      WNWNAT N+YT NBT NWNQ1/2 NYT NR1/2 N ANWN+BNYNWN0 0 Q1/2 NWN0I0 R1/2 NYN0 0 I      >0 (5.9a) s.t. Wij = 0, Yij = 0 ∀i, j such that xi∈C, xj/∈C(5.9b) then matrices PA=W−1 Nand KA=YNW−1 Nsatisfy (5.4), all the communication constraints imposed by the network mode Aand stabilize the whole system. Proof: The proof follows the same reasoning as the proof of Theorem 9. In this case the LMI constraint (5.9b) and the proposed variable change guarantee that the following inequality holds (AN+BNKA)TP(AN+BNKA)−PA+QN+KT ARNKA≤0.(5.10) Stability and (5.4) follow. The constraints imposed to the LMI guarantee that KAand PAsatisfy the communication restrictions of network mode A. Let PAand KAbe descomposed in blocks analogously to 140 Applications of Cooperative Game Theory WNand YN, that is PA={PA ij }, where PA ij ∈Rqi×qj, and KA={KA ij }, where KA i,j ∈Rri×qj. Without loss of generality let us assume that all the states of the subsystems that belong to a same communication component have been grouped, that is, PAcan be written as a block diagonal matrix PA=diag(PC)∀C∈N/A, where PC={PA ij }such that i, j ∈C∈N/A. Given that the inverse of a block diagonal matrix another block diagonal matrix in which the original blocks are inverted, that is WN=P−1 A=diag(P−1 C), it can be concluded that Wij = 0 ∀i, j such that xi∈C∈N/A, xj/∈C∈N/A implies that PA ij = 0 ∀i, j such that xi∈C∈N/A, xj/∈C∈N/A. Reordering the states in communication components also allows to write KAas a block diagonal matrix KA=diag(KC)∀C∈N/A, where KC={KA ij }such that i, j ∈C∈N/A. As KCand PCare dimensioned for the same state xCthen YN=KAWN=diag(KCP−1 C)∀C∈N/A. Thus, the fact that Yij = 0∀i, j such that ui∈C∈N/A, xj/∈C∈N/A is equivalent to make Kij = 0∀i, j such that ui∈C∈N/A, xj/∈C∈N/A. Remark: This theorem can be used also to provide a cost approximation for the agent game following a coherent criterium for all the agents that are not in communication component C. For example, the matrices PCcould be calculated assuming that all agents outside Cwork in a decentralized manner. 5.3 Distributed control simulation results In this section we show an academic example that illustrates the concepts and techniques presented in the chapter. The distributed system considered in the example is shown in figure 5.1. It consists of four agents, represented by boxes, which are coupled by pairs (the coupling interactions are represented by arrows). For example, agent 1 disturbs agents 2 and 3 and is 5.3. DISTRIBUTED CONTROL SIMULATION RESULTS 141 U 1 2 3 4 U U U 1 2 3 4 U1 U3 U2 U4 I IV II III Figure 5.1: Four systems coupled through the inputs. also disturbed by these two agents. The matrices that define the system are the following: A11 ="1 0.8 0 0.7#B11 ="0 1#B12 ="0 0.15#B13 ="0 0.15# A22 ="1 0.6 0 0.7#B21 ="0 0.15#B22 ="0 1#B24 ="0 0.15# A33 ="1 0.9 0 0.8#B31 ="0 0.15#B33 ="0 1#B34 ="0 0.15# A44 ="1 0.8 0 0.5#B42 ="0 0.15#B43 ="0 0.15#B44 ="0 1# Aij ="0 0 0 0#∀i6=j (5.11) where xi∈R2with i∈ {1, .., 4}are the states of each subsystem and ui∈Rwith i∈ {1, .., 4} are the corresponding inputs. The stage costs ℓiof all the subsystems are defined by matrices Qi=diag(1,1), Ri= 1 with i={1,2,3,4}. In order to implement the proposed control scheme, matrices KAand PAhave to be designed for each of the possible modes. Figure 5.2 shows the set of network modes for which each link is enabled. The number of possible modes is 16 and they have been numbered from 0 to 15. For example, in mode 0 no link is enabled while in mode 5 links I and II are enabled. Modes 11 to 15 have been omitted in the figure and only appear in the legend. These modes constitute all the cases where the grand coalition is formed, that is, there are at least 3