A Coalitional Distributed Model Predictive Control Perspective for a Cyber-Physical Multi-Agent Application

Abstract

Following the current technological development and informational advancement, more and more physical systems have become interconnected and linked via communication networks. The objective of this work is the development of a Coalitional Distributed Model Predictive Control (C- DMPC) strategy suitable for controlling cyber-physical, multi-agent systems. The motivation behind this endeavour is to design a novel algorithm with a flexible control architecture by combining the advantages of classical DMPC with Coalitional MPC. The simulation results were achieved using a test scenario composed of four dynamically coupled sub-systems, connected through an unidirectional communication topology. The obtained results illustrate that, when the feasibility of the local optimization problem is lost, forming a coalition between neighbouring agents solves this shortcoming and maintains the functionality of the entire system. These findings successfully prove the efficiency and performance of the proposed coalitional DMPC method.

1. Introduction

Presently, manifold systems are modular, interconnected and have a cyber-physical setup, meaning they can be viewed as coupled physical sub-systems, which are connected via communication networks [1,2,3,4,5]. For such processes, Distributed Model Predictive Control (DMPC) is a reliable control solution that uses local controllers that compute the control action using both (i) the local information derived from specific sensors and (ii) coupling data received/transmitted using the communication network [6].

As recent studies attest, the DMPC strategy was successfully applied on multi-agent systems in varying applications, such as formation control of autonomous surface and aerial vehicles [7], leader–follower platoons [8,9], traffic signal control [10], temperature regulation systems [11], battery energy storage systems [12] and microgrids [13,14]. In [15], a DMPC strategy for multi-agent systems based on error upper bounds is provided. This criterion is used in a min–max optimization of the cost function to minimize the communication between neighbouring agents. An event-triggered synchronous DMPC for multi-agent systems is introduced in [16]. The method is tailored for dynamically decoupled sub-systems, coupled through a cost function. An event-triggered mechanism designed using the forward difference of the cost function is deployed to activate the local optimization problem at each sampling time; otherwise the agents use the solutions computed in the previous sampling period. In [17], a DMPC to reach consensus for time-varying, multi-agent systems is proposed. The consensus DMPC algorithm is designed for heterogeneous, time-varying decoupled sub-systems, connected uni-directionally with a coupled cost function.

In all the research mentioned above, regardless of the application or the methodology details, one key feature is noticeable, namely that the architecture of both sub-systems and agents (i.e., local controllers) is fixed. The latter is predefined in the initialization phase of the control design, based on the sensors placements and interconnection between the local sub-systems. Therefore, the configuration of the DMPC neighbourhoods (i.e., groups of local sub-systems that are interconnected either dynamically or through cost functions or constraints) is established and predefined [18].

To overcome this shortcoming, a new approach emerged from cooperative game theory framework named Coalitional Control was introduced, with the following characteristics [19,20]: (i) the topology of the communication links between agents is flexible (i.e., links can be enabled or disabled when necessary), (ii) the control strategy encourages the agents to group in cooperative clusters called coalitions (to reduce the communication burden), and (iii) the communication links between agents that ensure their cooperation are weighted and introduce supplementary costs in the cost functions when activated.

Using this foundation, a Coalitional Model Predictive Control (C-MPC) strategy was developed and applied on different applications such as cellular networks in [21,22] or an eight-coupled tank process [23]. Thus, the agents charged with controlling the local sub-systems can form coalitions depending on the activation of the communication links between them. Several topologies can be derived, starting from the default one, which is the decentralized MPC strategy (i.e., with no communication between agents), to the most complex one. described as centralized MPC (i.e., in which all the communication links are active). In between, there are the coalitions between several agents in a neighbourhood, while the remaining ones work independently. In [24], details regarding the feasible regions for tube-based MPC controllers are discussed. The coalitions are associated with different partitions of a large-scale system (i.e., several sub-systems can be joined in a single entity), and their feasible region is analysed.

The main contribution of this work is the development of a novel perspective of the DMPC algorithm, called Coalitional DMPC (C-DMPC), which combines both the advantages and features of classical DMPC strategy with the characteristics of Coalitional MPC. Hence, the envisioned solution is to dynamically reshape the controller network by merging some of the agents within a neighbourhood into coalitions when needed. The advantages of this approach are shown when, due to various reasons (e.g., reference changes, unknown disturbances, etc.), one or more local optimal solutions become infeasible. In this case, to maintain the feasibility and functionality of the interconnected cyber-physical multi-agent system, the agents will decide to form a coalition. This means that inside a coalition, the sub-systems become a single entity and the controllers aggregate and solve a cooperative optimization problem (i.e., a global cost function is minimized) [25].

The main difference of our approach with respect to the cited coalitional literature is that the default topology is a non-cooperative DMPC (i.e., each agent minimizes a local cost function, using received information from its neighbours) [26]. This means that when a coalition occurs, the remaining agents outside the coalitions are not independent but retain their previous status and solve a non-cooperative optimization problem. Thus, depending on the topology, it is possible that the coalition must exchange information with its neighbours (if not all the agents inside a neighbourhood are merged into the coalition). Since all the agents start as non-cooperative players, they use the communication network to share relevant data, according to their coupling within neighbourhoods, and all pertinent communication links are activated and not weighted.

Another key difference is our proposed merging procedure, which selects the agents that will form a coalition. This is done at each agent level, without using a hierarchical supervisory layer. Moreover, when the local optimization problem becomes unfeasible (due to the coupling information), the coalition is activated. Furthermore, two simplified versions of this method with different agent merging procedures were published in [27,28]. In [27], each agent considers that the coupling information received from the neighboring agents is an uncertainty in the local nominal model. When a predefined threshold value for the local uncertainty level is crossed, than a coalition between the agents is formed. The further development of this idea is given in [28], in which the coalition between the agents is formed, when the local optimization problems become infeasible due to the received uncertainty level.

With respect to our previous papers, the method proposed in the current work has significant improvements, such as the following: (i) the network topology is tailored for in-chain coupled sub-systems, with unidirectional communication links; (ii) a more realistic academic example is used for simulation tests, with four heterogeneous sub-systems dynamically coupled through the inputs; (iii) each sub-system model is augmented with an additional state defined as the integral of control error to ensure a non-zero reference tracking; and (iv) multiple coalitions between agents can be simultaneously active at each sampling time.

The remaining of this paper is structured as follows: Section 2 presents the problem formulation and details the proposed method, whereas the simulation configuration, results, and discussions are provided in Section 3. The conclusions of this work and future work plans are addressed in Section 4.

2. Problem Formulation

A cyber-physical multi-agent system (CP-MAS), as depicted in Figure 1, is composed of N interconnected cyber-physical sub-systems (CPsS). Each CPsS is defined by the pair ($S_i$, $A_i$), $\forall i\in \mathcal{N}$, where $\mathcal{N}$ denotes the set $\{1,\dots ,N\}\subseteq \mathbb{N}$, with $N\in \mathbb{N}$ the number of sub-systems and $\mathbb{N}$ the set of natural numbers. The physical part of the CPsS is denoted with $S_i$, whereas the cyber part of the CPsS is denoted with $A_i$ and represents the corresponding local controller or agent. All the interconnected sub-systems $S_i$ form the physical layer (depicted with grey colour), while the cyber layer (depicted with blue colour) is composed of all the agents and the communication networks.

Let each sub-system $S_i$ be defined by the following model:

$$\begin{array}{ccc}\hfill x_p^i(k+1)& =& A_p^{i,i}{x}_p^i\left(k\right)+B_p^{i,i}{u}_p^i\left(k\right)+w_p^i\left(k\right)\hfill \\ \hfill w_p^i\left(k\right)& =& B_p^{i,i-1}{u}_p^{i-1}\left(k\right)\hfill \\ \hfill y_p^i\left(k\right)& =& C_p^i{x}_p^i\left(k\right),\forall i\in \mathcal{N}\hfill \end{array}$$

(1)

with the notations for the state $x_p^i\in {\mathbb{R}}^{n_i}$, input $u_p^i\in {\mathbb{R}}^{m_i}$, input uncertainty $w_p^i\in {\mathbb{R}}^{p_i}$ and output $y_p^i\in {\mathbb{R}}^{q_i}$. $A_p^{i,i},B_p^{i,i},B_p^{i,i-1}$ and $C_p^i$ are matrices with adequate dimensions. $n_i$, $m_i$, $p_i$ and $q_i$ are the number of states, inputs, input uncertainties and outputs, respectively. $\mathbb{R}$ denotes the set of real numbers. Note that, $u_p^{i-1}\in {\mathbb{R}}^{m_i}$ denotes the input signal received from the predecessor sub-system with index $i-1$.

Note that (1) defines a model in which the input-coupling information $w_p^i$ is considered an uncertainty in the nominal model. Moreover, all the sub-systems $S_i$, $\forall i\in \mathcal{N}$, are in a chain architecture, and for sub-system indexed i, the information is received through an unidirectional link from its predecessor and neighbour, defined as the sub-system with index $i-1$.

To ensure that the reference tracking control problem has a zero error in stationary regime, the state vector $x_p^i$, $\forall i\in \mathcal{N}$, from (1) is extended with an additional state ${\overline{x}}_p^i$ defined as integral of the control error, using the following definition [29,30]:

$${\overline{x}}_p^i(k+1)={\overline{x}}_p^i\left(k\right)-C_p^i{x}_p^i\left(k\right)+r_i\left(k\right),$$

(2)

obtaining

$$\begin{array}{ccc}\hfill \underset{x_i(k+1)}{\underbrace{\left[\begin{array}{c}{x}_p^i(k+1)\\ {\overline{x}}_p^i(k+1)\end{array}\right]}}& =& \underset{A_{i,i}}{\underbrace{\left[\begin{array}{cc}{A}_p^{i,i}& O\\ -C_p^i& I\end{array}\right]}}\underset{x_i\left(k\right)}{\underbrace{\left[\begin{array}{c}{x}_p^i\left(k\right)\\ {\overline{x}}_p^i\left(k\right)\end{array}\right]}}+\underset{B_{i,i}}{\underbrace{\left[\begin{array}{c}{B}_p^{i,i}\\ O\end{array}\right]}}{u}_i\left(k\right)\hfill \\ & +& \underset{B_{i,i-1}}{\underbrace{\left[\begin{array}{c}{B}_p^{i,i-1}\\ O\end{array}\right]}}{u}_{i-1}\left(k\right)+\underset{R_{{\mathrm{sp}}_i}}{\underbrace{\left[\begin{array}{c}O\\ I\end{array}\right]}}{r}_i\left(k\right)\hfill \\ \hfill y_i\left(k\right)& =& \underset{C_i}{\underbrace{\left[\begin{array}{c}{C}_p^{i}O\end{array}\right]}}\underset{x_i\left(k\right)}{\underbrace{\left[\begin{array}{c}{x}_p^i\left(k\right)\\ {\overline{x}}_p^i\left(k\right)\end{array}\right]}}\hfill \end{array}$$

(3)

where $r_i\left(k\right)$ is the imposed reference value at time k. $x_i\left(k\right)$, $u_i\left(k\right)$ and $y_i\left(k\right)$ are the extended state, input and output vectors, respectively. Note that the input uncertainty $w_i\left(k\right)$ is defined based on the input vector received from the predecessor $u_{i-1}\left(k\right)$. I and O are the identity and zero matrix, respectively, each with appropriate dimensions.

Hereafter, each sub-system $S_i$, $\forall i\in \mathcal{N}$, will be represented by the compact extended model:

$$\begin{array}{ccc}\hfill x_i(k+1)& =& A_{i,i}{x}_i\left(k\right)+B_{i,i}{u}_i\left(k\right)+w_i\left(k\right)+R_{{\mathrm{sp}}_i}{r}_i\left(k\right)\hfill \\ \hfill w_i\left(k\right)& =& B_{i,i-1}{u}_{i-1}\left(k\right)\hfill \\ \hfill y_i\left(k\right)& =& C_i{x}_i\left(k\right)\hfill \end{array}$$

(4)

where $A_{i,i}$, $B_{i,i}$, $B_{i,i-1}$, $C_i$ and $R_{{\mathrm{sp}}_i}$ are matrices with adequate dimensions.

Consider linear inequality constraints for the outputs, inputs and uncertainties defined with:

$$y_i\in {\mathcal{Y}}_i,u_i\in {\mathcal{U}}_i,w_i\in {\mathcal{W}}_i,\forall i\in \mathcal{N}$$

(5)

where ${\mathcal{Y}}_i$, ${\mathcal{U}}_i$ and ${\mathcal{W}}_i$ are sets defined by linear inequalities.

At every sampling time, each agent $A_i$, $\forall i\in \mathcal{N}$, solves a min–max optimization problem, which aims to obtain the minimum optimal input with respect to the maximum level of uncertainty received from its neighbour.

$$\begin{array}{cc}\hfill J_i\left(x_i^0\right)& =\underset{\begin{array}{c}{u}_i^0,\dots ,u_i^{N_p-1}\\ {\widehat{u}}_i^{\max }\end{array}}{\min }\underset{w_i^0,\dots ,w_i^{N_p-1}}{\max }\hfill \\ \hfill & \Vert R_{w_i}{\widehat{u}}_i^{\max }{\Vert }_1+\sum _{l=0}^{N_p-1}\Vert r_i^l-y_i^l{\Vert }_1+\Vert R_{u_i}{u}_i^l{\Vert }_1+\Vert r_i^{N_p}-y_i^{N_p}{\Vert }_1\hfill \\ \hfill \mathrm{s}.\mathrm{t}.\left(4\right)\\ \hfill & y_i^l\in {\mathcal{Y}}_i,l=1,\dots ,N_p-1\hfill \\ \hfill & y_i^{N_p}\in {\mathsf{\Omega }}_i\hfill \\ \hfill & u_i^l\le {\widehat{u}}_i^{\max }\le u_i^{\max },l=0,\dots ,N_p-1\hfill \\ \hfill & w_i^l\le w_i^{\max }\hfill \end{array}$$

(6)

where $y_i^l=y_i(k+l|k)$ denotes the output predictions for sub-system $S_i$ at time $k+l$, computed at time step k; this is calculated recursively starting from the initial state $x_i^0=x_i\left(k\right)$ measured at time $k$, using the model (4); the input sequence $u_i^l=\left\{u_i^0,\dots ,u_i^{N_p-1}\right\}$ computed over the prediction horizon $N_p$; and the uncertainty sequence ${\mathcal{W}}_i=B_{i,i-1}{\mathcal{U}}_{i-1}$ received from the neighbour (where ${\mathcal{W}}_i$ is the uncertainty polytope and ${\mathcal{U}}_{i-1}=\{u\in {\mathbb{R}}^{m_i}:Au\le b$ is a H-polytope); $r_i^l$ is the value at time $k+l$ for the output reference trajectory; $r_i^{N_p}$ and $y_i^{N_p}$ are the values for the reference and the output trajectories, at the end of the prediction horizon $k+N_p$, respectively; $u_i^{\max }$, $w_i^{\max }$ are the maximum limits for the input and the uncertainty sequences, respectively; ${\Vert .\Vert }_1$ denotes the 1-norm; $R_{u_i}\in {\mathbb{R}}^{m_i\times m_i}$ and $R_{w_i}\in {\mathbb{R}}^{m_i\times m_i}$ are the weight matrices for the input and self-imposed input limit ${\widehat{u}}_i^{\max }$. The latter is an additional optimization parameter introduced in the local cost function, and its value is communicated at each sampling period to the neighbour. This will guarantee that the uncertainty level received from the neighbour is smaller than this value, without actually transmitting the entire input sequence. The set ${\mathsf{\Omega }}_i$ is a robust positive invariant set used to ensure the closed-loop stability of the algorithm by means of the terminal invariant set.

Remark 1.

The uncertainty in each sub-system model refers to the coupling information that must be received from the neighbouring sub-system. Please note that the local optimization problem (6) minimizes the control input, for the worst-case scenario related to uncertainty level received from the predecessor agent. This means that, although unknown, this uncertainty must be bounded to a known value, which is shared between consecutive sub-system. Moreover, this ensures that each local sub-system is prepared for the disturbance signal, which is received via the coupling links.

Next, some details regarding the computation of the invariant set ${\mathsf{\Omega }}_i$, followed by the proposed coalitional DMPC method are given.

2.1. Robust Positive Invariant Set Computation

In this sub-section, the details regarding the computation of the robust positive invariant set ${\mathsf{\Omega }}_i$, $\forall i\in \mathcal{N}$, which acts as a constraint region for the terminal state $y_i^{N_p}\in {\mathsf{\Omega }}_i$ are presented. To this end, the procedure firstly introduced in [28] is briefly summarized below, tailored for the extended sub-system model.

For each sub-system $S_i$, $\forall i\in \mathcal{N}$, with the model defined in (4) and subject to constraints (5), only the nominal model (i.e., $w_i$ and $r_i$ are zero) is considered. Let us compute a local linear feedback $u_i=K_i{x}_i$, which ensures that the closed loop eigenvalues are in the unit circle. One suggestion to compute the state feedback matrix $K_i$ is to apply classical state-space feedback control designed for the nominal model using Ackermann’s formula (i.e., solving a pole allocation problem) [31], or to calculate it through the minimization of a linear-quadratic cost function, by solving a discrete-time Riccatti Equation [32].

The set ${\mathsf{\Omega }}_i$ is robust positive invariant for the nominal model from (4), if the following assumption holds [28,33]:

$$\begin{array}{cc}\hfill & x_i\in {\mathsf{\Omega }}_i\to (A_{i,i}+B_{i,i}{K}_i)x_i+w_i\in {\mathsf{\Omega }}_i,\\ & K_i{x}_i\in {\mathcal{U}}_i,C_i{\mathsf{\Omega }}_i\subseteq {\mathcal{Y}}_i,\forall w_i\in {\mathcal{W}}_i\end{array}$$

(7)

It is worth mentioning the following observations regarding the use of the invariant set in the C-DMPC context:

• the default working framework is non-cooperative DMPC, which implies that each agent $A_i$, $\forall i\in \mathcal{N}$, from the multi-agent application communicates with its neighbour, in order to compute the local solution;
• each sub-system model $S_i$, $\forall i\in \mathcal{N}$, is subject to input uncertainties received from the sub-system to whom it is connected (in our case its predecessor);
• to provide a simplified algorithm with minimal communication load in the network, only the self-imposed upper bound for the local input trajectory is broadcast in the network (i.e., the optimization variable ${\widehat{u}}^{\max }$ introduced in (6));
• a table with different predefined robust positive invariant sets ${\mathsf{\Omega }}_i$ is computed using the constraints limits from (5), in which each element is a particular combination of the variable bounds (see Algorithm 1);
• at each sampling period, after the uncertainty upper bound is received from the neighbour, each agent $A_i$ uses this information to compute the uncertainty polytope. Next, from the predefined terminal sets table, a set ${\mathsf{\Omega }}_i$ is searched for, which includes the received uncertainty polytope (i.e., which will ensure a local feasible solution in the terminal state framework).

Further on, the pseudo-code algorithm used to compute the invariant set table is provided (where for simplicity the sub-system indices are omitted):

Thus, each agent $A_i$, $\forall i\in \mathcal{N}$, uses Algorithm 1 in the initialization phase of the proposed method to compute a table of invariant sets ${\mathsf{\Omega }}_i$, for different input and uncertainty parametrizations (i.e., distinct combinations for the two parameters $\alpha$ and $\beta$). Note that the first set ${\mathsf{\Omega }}_i$ from the table corresponds to the largest value for the input constraint, denoted $u^{max}$, whereas the uncertainty has the smallest value. The latter is gradually increased with a step size denoted $ste{p}_{\beta }$, until it reaches its maximum admissible value $w^{max}$. In doing so, the size of the invariant set slowly reduces, as the input constraint limit value decreases with a step size denoted $ste{p}_{\alpha }$ and the uncertainty level rises.

In practice, a good start for $u^{max}$ and $w^{max}$ bounds are the values for the imposed constraints (5). The values of the step size $ste{p}_{\alpha }$, $ste{p}_{\beta }$ should be selected such that the table size remains reasonable, with various invariant sets. Moreover, the limits in the state constraints are considered fixed, according to the sub-systems dynamics and used to compute every set ${\mathsf{\Omega }}_i$ from the table.

Algorithm 1

For$\alpha =u^{max}:-ste{p}_{\alpha }:0.1$

For $\beta =0.1:ste{p}_{\beta }:w^{max}$

1. Compute the inequality constraints:

$$\begin{array}{c}\hfill A_{u}u\le \alpha b_u;A_{w}w\le \beta b_w;A_{x}x\le b_x\end{array}$$

2. Compute the robust positive-invariant set:

$$\begin{array}{c}\hfill \mathsf{\Omega }(A,B,K,A_u,b_u,A_w,b_w,A_x,b_x)\end{array}$$

3. Save the information $\alpha$, $\beta$, $\mathsf{\Omega }$

end

end

2.2. Coalitional Distributed Model Predictive Control (C-Dmpc) Methodology

As previously mentioned, what differentiates our proposed coalitional algorithm from the existing works is the flexible framework set for the cyber-physical multi-agent system with a chain architecture. Hence, at each step time, the agents architecture starts as non-cooperative DMPC and will switch to coalitional DMPC (C-DMPC)—when the local feasibility of the interconnected agents is lost. In the C-DMPC framework, the coalition procedure is initialized without a hierarchical level by the local agents with infeasible problems, because due to the coupling links between sub-systems, if not solved, this problem will propagate among neighbouring sub-systems. Using the communication links, these agents share their optimization status with their neighbour, and after that, one of them is randomly selected to start a coalition. Once the coalition procedure is activated, the agents framework changes.

To simplify the design and computational costs, the size of the coalition is increased gradually, if needed. That is, if a coalition of two agents, coupled with the remaining agents from the network, still does not provide feasible solutions for all involved actors, then more work needs to be done. The idea is to first activate all coalitions of two agents, if needed, then the coalitions of three agents, and so on, until in the end, in the extreme case, all the agents are involved in a single coalition. Note that this last case is equivalent to solving a centralized problem for the multi-agent system and will be used in the last resort, if nothing else solved the infeasibility problems that started the coalitional procedure. The reason for this is related with the coalition dynamics (i.e., when two or more agents form a coalition, their respective sub-system models are aggregated and become a single entity). Thus, the number of the optimization variables in a coalition increases with its size, and the local non-cooperative optimization problem becomes a cooperative one inside the coalition. The extreme case of a ‘grand’ coalition between all agents will aggregate all the sub-systems in a single entity (from the control point of view).

2.2.1. Coalition Dynamics

As described before, our C-DMPC algorithm is tailored specifically for cyber-physical multi-agent systems, linked in a unidirectional communication topology. Thus, the coupling information, which is treated as an uncertainty in the local nominal model of each sub-system $S_i$, is received from its predecessor sub-system $S_{i-1}$. To minimize the communication burden between consecutive agents, only the self imposed optimization variable ${\widehat{u}}^{\max }$ introduced in (6) is broadcast. This value is firstly used to search for an invariant set inside the predefined table, and secondly acts as the uncertainty limit constraint in the local optimization problem. Using this information, the local optimization problem is then solved, and if the solution is infeasible, then the coalition procedure must be started.

Inside a coalition between different consecutive agents, the aim is to solve a cooperative optimization problem; thus the uncertainty variable becomes fully known. Each agent $A_i$, $\forall i\in \mathcal{N}$, can form a coalition only with its predecessor, i.e., agent $A_{i-1}$, due to the particular dynamical coupling between their corresponding sub-systems (i.e., linked in a chain). When this occurs, the agents involved will form a compact set denoted generically $\mathcal{C}$. To simplify the notations, the coalition is described without sub-script indices with the following model:

$$\begin{array}{ccc}\hfill x_{\mathcal{C}}(k+1)& =& A_{\mathcal{C}}{x}_{\mathcal{C}}\left(k\right)+B_{\mathcal{C}}{u}_{\mathcal{C}}\left(k\right)+w_{\mathcal{C}}\left(k\right)\hfill \\ \hfill w_{\mathcal{C}}\left(k\right)& =& \sum _{j\in {\mathcal{N}}_{\mathcal{C}}}{B}_{\mathcal{C}}^j{u}_j\left(k\right)\hfill \\ \hfill y_{\mathcal{C}}\left(k\right)& =& C_{\mathcal{C}}{x}_{\mathcal{C}}\left(k\right)\hfill \end{array}$$

(8)

where $x_{\mathcal{C}}$ is the state vector of the coalition, $u_{\mathcal{C}}$ is the coalition’s input vector, $w_{\mathcal{C}}$ is the uncertainty vector of the coalition and $y_{\mathcal{C}}$ is output vector for the coalition. All these vectors are composed by aggregating the local vectors corresponding to each sub-system involved in the coalition (e.g., $x_{\mathcal{C}}={\left[x_i\right]}_{i\in \mathcal{C}}$). Moreover, the matrices $A_{\mathcal{C}}$, $B_{\mathcal{C}}$, $B_{\mathcal{C}}^j$ and $C_{\mathcal{C}}$ are computed according to the aggregation.

The set ${\mathcal{N}}_{\mathcal{C}}$ denotes the coalition’s $\mathcal{C}$ neighbour, defined as the predecessor sub-system for the sub-systems inside the coalition (e.g., if Agent 2 and 3 form a coalition, then ${\mathcal{N}}_{\mathcal{C}}=\left\{1\right\}$, because sub-system 2 is coupled to sub-system 1; thus the coalition in which Agent 2 is involved must receive relevant information from Agent 1, which is outside the coalition and solves a non-cooperative DMPC problem). Moreover, following this reasoning, a coalition involving Agent 1 does not have neighbours (i.e., ${\mathcal{N}}_{\mathcal{C}}=\varnothing$, because Agent 1 does not have predecessors).

2.2.2. Coalition Problem Definition

In this section, some details regarding the construction of the constraints sets imposed for the coalition and the optimization problem solved by the coalition are presented.

Hence, the constraint sets for the coalition $\mathcal{C}$ are computed as the union of the constraints sets (5) corresponding to each agent $A_i$, $i\in \mathcal{C}$:

$$\begin{array}{c}\hfill y_{\mathcal{C}}\in {\mathcal{Y}}_{\mathcal{C}}=\prod _{i\in \mathcal{C}}{\mathcal{Y}}_i,u_{\mathcal{C}}\in {\mathcal{U}}_{\mathcal{C}}=\prod _{i\in \mathcal{C}}{\mathcal{U}}_i,w_{\mathcal{C}}\in {\mathcal{W}}_{\mathcal{C}}=\prod _{i\in \mathcal{C}}{\mathcal{W}}_i,\end{array}$$

(9)

and the min–max optimization problem solved by the coalition is:

$$\begin{array}{ccc}\hfill J_{\mathcal{C}}\left(x_{\mathcal{C}}^0\right)& =& \underset{\begin{array}{c}{u}_{\mathcal{C}}^0,\dots ,u_{\mathcal{C}}^{N-1}\\ {\widehat{u}}_{\mathcal{C}}^{\max }\end{array}}{\min }\underset{w_{\mathcal{C}}^0,\dots ,w_{\mathcal{C}}^{N-1}}{\max }\hfill \\ & & \Vert R_{w_{\mathcal{C}}}{\widehat{u}}_{\mathcal{C}}^{\max }{\Vert }_1+\sum _{l=0}^{N_p-1}\Vert r_{\mathcal{C}}^l-y_{\mathcal{C}}^l{\Vert }_1+\Vert R_{u_{\mathcal{C}}}{u}_{\mathcal{C}}^l{\Vert }_1+\Vert r_{\mathcal{C}}^{N_p}-y_{\mathcal{C}}^{N_p}{\Vert }_1\hfill \\ \hfill \mathrm{s}.\mathrm{t}.(8)\\ \hfill y_{\mathcal{C}}^l& \in & {\mathcal{Y}}_{\mathcal{C}},l=1,\dots ,N_p-1\hfill \\ \hfill y_{\mathcal{C}}^{N_p}& \in & {\mathsf{\Omega }}_{\mathcal{C}}\hfill \\ \hfill u_{\mathcal{C}}^l& \le & {\widehat{u}}_{\mathcal{C}}^{\max }\le u_{\mathcal{C}}^{\max }\hfill \\ \hfill w_{\mathcal{C}}^l& \le & w_{\mathcal{C}}^{\max }\hfill \end{array}$$

(10)

The weighting matrices $R_{u_{\mathcal{C}}}$ and $R_{w_{\mathcal{C}}}$ are block diagonal, ${\mathsf{\Omega }}_{\mathcal{C}}$ is the aggregated terminal set and $r_{\mathcal{C}}^l$ and $x_{\mathcal{C}}^0$ are aggregated vectors containing the corresponding imposed references and initial state values, respectively. $r_{\mathcal{C}}^{N_p}$ and $y_{\mathcal{C}}^{N_p}$ are aggregated vectors containing the corresponding imposed references and output predictions values at time $k+N_p$, respectively. $u_{\mathcal{C}}^l$, $w_{\mathcal{C}}^l$ and $y_{\mathcal{C}}^l$ are are aggregated vectors containing the corresponding input, uncertainty and output sequences, respectively. ${\widehat{u}}_{\mathcal{C}}^{\max }$ is an aggregated vector containing the corresponding self-imposed input limits. $u_{\mathcal{C}}^{\max }$ and $w_{\mathcal{C}}^{\max }$ are aggregated vectors containing the corresponding input and uncertainty limits, respectively.

2.2.3. C-Dmpc Algorithm

To summarize the C-DMPC methodology, the following pseudo-code is provided: With regard to Algorithm 2, the following observations are in order:

• the default uncertainty value used in Step 1 is selected to ensure that optimization problems from Step 3 are feasible, thus ensuring that the proposed methodology is recursively stable (i.e., the terminal set for the coalition is obtained by aggregating the terminal sets of the involved individual agents).
• if the condition from Step 6 is satisfied, then at that sampling period, the working framework is non-cooperative DMPC; otherwise the framework changes to coalitional DMPC (since at least one coalition is activated).
• the priority value, which is used as a condition term to initialize a coalition, is defined by each agent as a random sub-unitary number. In this manner, there is no use of a hierarchical control level to assign these priorities.
• in the extreme, all the agents can be combined in a coalition ($\mathcal{C}=\mathcal{N}$), which corresponds to a centralized MPC working framework.
• one or more coalitions can be active simultaneously and are dissolved at the end of each sampling period.

Remark 2

([28]). In Algorithm 2, the coalitional control problem is feasible (i.e., Step 6. (c). ii.), because ${\mathcal{W}}_{\mathcal{C}}\subseteq {\prod }_{i\in \mathcal{C}}{\mathcal{W}}_i$, ${\mathcal{U}}_{\mathcal{C}}={\prod }_{i\in \mathcal{C}}{\mathcal{U}}_i$ and ${\mathsf{\Omega }}_{\mathcal{C}}={\prod }_{i\in \mathcal{C}}{\mathsf{\Omega }}_i$. The stability of the coalition is ensured by the terminal constraint set of the coalition, which is calculated as the Minkowski sum of the terminal sets polytopes defined for each individual agent from the coalition. The coalitional algorithm is recursive-feasible, contingent on Step 3, for which all the optimization problems are feasible, i.e., for which systems can work in a decentralized fashion.

Next, the C-DMPC methodology is validated in simulation, and the results are provided in Section 3.

3. Illustrative Example

In this section, the simulation results and discussion for the C-DMPC method are presented. The proposed simulation scenario for the cyber-physical multi-agent system described in Section 2, Figure 1, has the following characteristics:

• Four heterogeneous discrete-time sub-systems $S_i$, $\forall i\in \{1,\dots ,4\}$, coupled in a chain architecture were defined using (1), with the following numerical matrices:

  $$\begin{array}{c}\hfill S_1:A_p^{1,1}=\left[\begin{array}{cc}0.7913& 0.2020\\ 0.1010& 0.8417\end{array}\right]B_p^{1,1}=\left[\begin{array}{c}0.0271\\ 0.2291\end{array}\right]\\ \hfill B_p^{1,0}=\left[\begin{array}{c}0\\ 0\end{array}\right]C_p^1=\left[\begin{array}{cc}0& 1\end{array}\right]\end{array}$$

  (11)

  $$\begin{array}{c}\hfill S_2:A_p^{2,2}=\left[\begin{array}{cc}0.7936& 0.1996\\ 0.1198& 0.8236\end{array}\right]B_p^{2,2}=\left[\begin{array}{c}0.0269\\ 0.2265\end{array}\right]\\ \hfill B_p^{2,1}=\left[\begin{array}{c}0.0004\\ 0.0034\end{array}\right]C_p^2=\left[\begin{array}{cc}0& 1\end{array}\right]\end{array}$$

  (12)

  $$\begin{array}{c}\hfill S_3:A_p^{3,3}=\left[\begin{array}{cc}0.7888& 0.2043\\ 0.0817& 0.8604\end{array}\right]B_p^{3,3}=\left[\begin{array}{c}0.0273\\ 0.2316\end{array}\right]\\ \hfill B_p^{3,2}=\left[\begin{array}{c}0.0004\\ 0.0034\end{array}\right]C_p^3=\left[\begin{array}{cc}0& 1\end{array}\right]\end{array}$$

  (13)

  $$\begin{array}{c}\hfill S_4:A_p^{4,4}=\left[\begin{array}{cc}0.7912& 0.1994\\ 0.0997& 0.8211\end{array}\right]B_p^{4,4}=\left[\begin{array}{c}0.0269\\ 0.2263\end{array}\right]\\ \hfill B_p^{4,3}=\left[\begin{array}{c}0.0004\\ 0.0034\end{array}\right]C_p^4=\left[\begin{array}{cc}0& 1\end{array}\right]\end{array}$$

  (14)
• The limit constraints for the inputs, disturbances and outputs are the following:

  $$\begin{array}{c}\left\{\begin{array}{c}{u}_i^{\min }=-5\hfill \\ u_i^{\max }=5\hfill \end{array}\right.\hfill \end{array}\left\{\begin{array}{c}{w}_i^{\min }=-1\hfill \\ w_i^{\max }=1\hfill \end{array}\right.\left\{\begin{array}{c}{w}_i^{\min }=-8\hfill \\ w_i^{\max }=8\hfill \end{array}\begin{array}{}\forall i\in \{1,\dots ,4\}\end{array}\right.$$

  (15)
• For all sub-systems $S_i$, $\forall i\in \{1,\dots ,4\}$, the following optimization parameters are used: the prediction horizon $N_p=5$, the input weights $R_{u_i}=0.1$ and $R_{w_i}=0.01$.

Remark 3.

The optimization parameters were carefully selected after a thorough analysis from the point of view of achieved performances. Several tests were performed, with different values for the weights and the prediction horizon. The chosen values ensured the best performances.

• The feedback laws were computed using classical state-feedback control based on the Ackermann’s formula [31], applied for the extended model (4), obtaining:

  $$\begin{array}{c}{K}_1=[0.5494-2.60610.7488],\hfill \\ K_2=[0.3473-2.51990.7047],\hfill \\ K_3=[0.6249-2.65650.7401],\hfill \\ K_4=[0.4355-2.51070.7051].\hfill \end{array}$$

  (16)

Remark 4.

The Ackermann’s formula [31] was used to achieve specific closed-loop transient performances, chosen as an overshoot value of $5\%$ and settling time of 5 time units, for sub-systems $S_1$ and $S_3$, and an overshoot value of $4\%$, and the same settling time, corresponding to sub-systems $S_2$ and $S_4$. These performance values, were accordingly selected based on each sub-system dynamics.

• The reference tracking scenario was constructed for 12 time samples, using a sampling period $T_s=0.25s$, with the following imposed references:

  $$\begin{array}{c}{r}_1=\left[0.20.20.20.20.20.20.50.50.50.50.50.5\right],\hfill \\ r_2=\left[0.20.20.20.20.20.20.20.20.20.20.20.2\right],\hfill \\ r_3=\left[0.20.20.20.20.50.50.50.50.20.20.20.2\right],\hfill \\ r_4=\left[0.20.20.20.20.20.20.20.20.20.20.20.2\right].\hfill \end{array}$$

  (17)

Since our proposed scenario has four sub-systems in a chain architecture with unidirectional communication links between the agents, there are eight possible frameworks including coalitions of two, three or four agents defined as follows:

1.

default case—no coalitions between $A_1$, $A_2$, $A_3$, $A_4$;

2.

coalition ${\mathcal{C}}_{12}$ between $A_1$ and $A_2$, while $A_3$, $A_4$ remain outside the coalition but interconnected;

3.

coalition ${\mathcal{C}}_{123}$ between $A_1$, $A_2$, and $A_3$, while $A_4$ remains outside the coalition but interconnected;

4.

twosimultaneous active coalitions $C_{12}$ and ${\mathcal{C}}_{34}$ between $A_1$ and $A_2$ and $A_3$ and $A_4$, respectively, which are interconnected;

5.

coalition ${\mathcal{C}}_{23}$ between $A_2$ and $A_3$, while $A_1$, $A_4$ remain outside the coalition but interconnected;

6.

coalition ${\mathcal{C}}_{234}$ between $A_2$, $A_3$ and $A_4$, while $A_1$ remains outside the coalition but interconnected;

7.

coalition ${\mathcal{C}}_{34}$ between $A_3$ and $A_4$, while $A_1$, $A_2$ remain outside the coalition but interconnected;

8.

extreme case: coalition ${\mathcal{C}}_{1234}$ between all agents $A_1$, $A_2$, $A_3$, $A_4$.

Algorithm 2

Initialization: For each agent $A_i$, $\forall i\in \mathcal{N}$, compute a table ${\mathcal{T}}_i$, with potential terminal sets ${\mathsf{\Omega }}_i$.

At each sampling time k, each agent $A_i$, $\forall i\in \mathcal{N}$, receives the local state value and performs the following steps:

1. Computes the uncertainty polytope using default limit values for the constraints:

$$\begin{array}{ccc}\hfill {\mathcal{W}}_i& =& B_{i,i-1}{\mathcal{U}}_{i-1}^0\hfill \\ \hfill U_{i-1}& =& [u_{i-1}^{\max ,0};-u_{i-1}^{\max ,0}].\hfill \end{array}$$

2. Searches in the predefined table ${\mathcal{T}}_i$ for a terminal set ${\mathsf{\Omega }}_i^0$ that includes the default uncertainty ${\mathcal{W}}_i\subseteq {\mathsf{\Omega }}_i^0$.

3. Solves the local optimization problem (6) and obtains the optimal values $U_i^{*,0}$, ${\widehat{u}}_i^{\max ,0}$ using the default values ${\mathsf{\Omega }}_i={\mathsf{\Omega }}_i^0$ for the terminal set and the uncertainty constraint limit ($w_i^{\max }=u_{i-1}^{\max ,0}$).

4. Broadcasts to its successor the local optimal value ${\widehat{u}}_i^{\max ,0}$ and receives the corresponding value ${\widehat{u}}_{i-1}^{\max ,0}$ from its predecessor.

5. Repeats Steps 1–3 using the uncertainty constraint value received in Step 4.

6. Checks the feasibility of the local optimization problem:

If the optimization problem from Step 5 is feasible:

then: Coalitions between agents are not necessary.Each local agent $A_i$ sends to its sub-system $S_i$, the first value from the optimal trajectory $U_i^{*}$;

else: Coalitions between agents are necessary. In this case, in order to be included in a coalition, each agent $A_i$, $\forall i\in \mathcal{N}$, performs the following steps:

a. Receives, from its predecessor, a coalitional report containing the following information: the feasibility status (for the local optimization problem solved at Step 5) and priority value relating to all the predecessor agents from the chain architecture.

b. Sends to its successor, the updated coalitional report (i.e., all the relevant information received, together with its own local feasibility and priority data).

c. Initializes a coalition only if its local priority is the highest from the report. Within a coalition between two agents, the following steps are performed:

i. the coalition model is defined as (8);

ii. the optimization problem (10) subject to (9) is solved.

iii. the relevant information is broadcast to the coalition’s neighbour.

iv. a feasibility check for all the optimization problems is done.

If the all the optimization problems are feasible:

then: The existing coalition was successful and can be dissolved after every sub-system $S_i$ receives the first value from the optimal trajectory $U_i^{*}$;

else: The existing coalition was not successful. Another agent must be included in the existing coalition (if the coalition’s status is infeasible), or another coalition can be activated (if more agents outside the existing coalition have infeasible problems). At this stage, Step (c) is repeated as necessary.

7. End algorithm.

Remark 5.

Please note that our proposed Coalitional DMPC algorithm is tailored specifically for cyber-physical multi-agent systems. The key feature is its capability to switch between control architectures, whenever the feasibility of the multi-agent system is lost, due to uncertainties in the local sub-systems. One example of such multi-agent system is a vehicle platoon. In this case, it is clear that classical centralized MPC is not suitable for controlling this application. Moreover, decentralized MPC, in which the couplings between sub-systems are ignored, can render instability within the platoon. One compromise solution is distributed MPC, in which the interactions are taken into account when computing the local solutions. However, if the distributed MPC (i.e., in non-cooperative framework) fails at this task, then our proposed coalitional DMPC provides a backup plan, namely to merge different sub-systems into coalitions. Inside a coalition, all the information is known; thus, only the coupling signals with sub-systems outside the coalition must be accounted for.

The invariant sets obtained for sub-system $S_1$ using Algorithm 1 presented in Section 2.1 are depicted in Figure 2. For the computation, the following numerical values were used: $u^{max}=5$, $w^{max}=5$, $ste{p}_{\alpha }=ste{p}_{\beta }=0.5$. As expected, the larger invariant set (depicted with red colour) was obtained for $\alpha =u^{max}=5$ and $\beta =0.1$. Moreover, as the constraint limits become smaller, the set ${\mathsf{\Omega }}_i$ decreases in dimension and is included in the larger red set (Figure 3—the sets plotted with green, blue, magenta and black colours, respectively). Since, the state variable for the extended model (4) has three values, the computed invariant sets are three dimensional and can be plotted as convex hulls (ref. Figure 2 and Figure 3). This graphical representation of the invariant sets, which are predefined options for the terminal set constraint (6), are also useful when defining the reference target for the multi-agent system. Thus, one must take into account that the imposed trajectory for each sub-system $S_i$, $\forall i\in \{1,\dots ,4\}$ should be placed in the interior of the invariant set.

Remark 6.

It is worth mentioning that the step values $ste{p}_{\alpha }=ste{p}_{\beta }=0.5$ were selected taking into account the numerical values of the input and uncertainty constraints to ensure a sufficient number of invariant sets computed. If a smaller value, e.g., $ste{p}_{\alpha }=ste{p}_{\beta }=0.1$ is chosen, the result would be an increased size for the table containing the invariant sets. However, as depicted in Figure 2, these values also parametrize the dimensions of the invariant set polytopes. Thus, although we would have more available sets, their dimensions would be too similar, to justify the involved computational costs.

The reference tracking results and the formation of the coalitions during the simulation are presented in Figure 4.

As depicted in Figure 4, lower subplot, during the first seven time steps, the simulation runs in the default scenario, in which all the agents solve a non-cooperative DMPC algorithm without being involved in a coalition. This is marked with blue circles, at each time step, for each agent $A_i$, $\forall i\in \{1,\dots ,4\}$. At time step 8, due to the setpoint change of 0.3 units in $r_3$ for sub-system $S_3$ and the corresponding increase in the control effort $u_3$, the local feasibility for sub-system $S_4$ is lost. Hence, the coalition ${\mathcal{C}}_{34}$ is activated, which is plotted with a red star marker for $A_3$ and $A_4$. At the next time step, coalition ${\mathcal{C}}_{234}$ between agents $A_2$, $A_3$ and $A_4$ is active and is coupled with the remaining agent $A_1$, because sub-system $S_2$ is dynamically coupled through the input with $S_1$ and their corresponding agents share information. At time step 10, coalition ${\mathcal{C}}_{34}$ is activated, and for the remaining two time steps of the simulation, coalition ${\mathcal{C}}_{234}$ is active. Moreover, the reference tracking results show that all the imposed set-points are successfully reached in one sampling time, with zero offset error. This occurs for the first seven time steps, in which all the agents work outside a coalition, and also for the remaining simulation time, when coalitions of two or three agents are necessary to maintain the feasibility of the CP-MAS. The results clearly prove the efficiency of our proposed C-DMPC method in a reference tracking scenario.

4. Conclusions

In this work, a coalitional distributed model predictive (C-DMPC) methodology suitable for input coupled cyber-physical multi-agent systems was proposed. The algorithm was tailored for an in-chain system architecture with unidirectional communication links (i.e., the coupling information viewed as an uncertainty in the local nominal model was broadcast from a predecessor sub-system to a successor). The methodology was validated in simulation, using an academic cyber-physical multi-agent system as a proof of concept for the proposed algorithm. The simulation results show that if the uncertainty level received by the local agent is manageable, a non-cooperative DMPC algorithm could be locally solved. However, when the local feasibility of the optimization problem was lost, then forming coalitions between agents showed satisfactory performance and the usefulness of the C-DMPC algorithm was proven.

Future work will test the efficiency of the proposed algorithm on a vehicle platooning application.

5. Materials and Methods

The simulations from this work were performed using $\mathrm{MATLAB}$ R2020b on Windows 10, 64-bit Operating System with a laptop $\mathrm{Intel}$ $\mathrm{Core}$ i7-9850H CPU @ 2.60 GHz and 16 GB RAM.

The optimizations were implemented using the YALMIP toolbox [34].