# Network-Realised Model Predictive Control
## Part II: Distributed Constraint Management

Andrei Sperilă<sup>a</sup>, Alessio Iovine<sup>b</sup>, Sorin Olaru<sup>b</sup>, Patrick Panciatici<sup>b</sup>

### Abstract

A two-layer control architecture is proposed, which promotes scalable implementations for model predictive controllers. The top layer acts as both reference governor for the bottom layer, and as a feedback controller for the regulated network. By employing set-based methods, global theoretical guarantees are obtained by enforcing local constraints upon the network's variables and upon those of the first layer's implementation. The proposed technique offers recursive feasibility guarantees as one of its central features, and the expressions of the resulting predictive strategies bear a striking resemblance to classical formulations from model predictive control literature, allowing for flexible and easily customisable implementations.

**Key words:** Distributed control; convex optimisation; model predictive control; networked systems; scalable implementations.

## 1 Introduction

### 1.1 Setting and Context

Model Predictive Control (MPC, see [19] and also [5]) is one of the most popular control-based frameworks in both the industrial setting and in system-theoretical literature. One of the key features that has generated its current popularity is its affinity for enforcing set-based constraints in the centralised setting [4]. On the other hand, as technical applications are becoming increasingly interconnected and widely distributed, the synthesis of distributed control laws has evolved into one of the foci of contemporary research for the control community.

Among the myriad approaches developed for the MPC paradigm (see [15] for a comprehensive survey on the topic), some of the most promising solutions are those inspired by robust MPC (see [16] for a popular centralised setting), among which we highlight the techniques proposed in [8,26,25,11] for the distributed setting. A key aspect of this particular framework is given by the use of an ancillary control law, based upon static state feedback, which operates in tandem with its MPC policy.

Distinct from the latter approach, the procedure formalised in [1,2] aims to leverage the powerful design framework from [29] in order to tackle the intricacies of the distributed setting. One such intricacy, which arises chiefly in formation control, is the problem of attenuating the propagation of disturbance through said formation (see [9]). This problem has received a great deal of attention in MPC-oriented literature, with multiple solutions having been proposed in both the distributed setting [31,28,6,10], and also in the centralised one [27].

### 1.2 Motivation

Despite the valuable contributions made in the previously cited works, certain limitations which are inherent in these approaches serve as the chief motivation for the developments presented in this paper. More precisely, the control laws proposed in [8,26,11] notably require a set of *static and stabilising* state-feedback strategies, which are constrained to having block-diagonal sparsity structures. Furthermore, the procedure described in [25] revolves around intricate set-theoretical notions, such as *joint set invariance*, while implicitly assuming that the controlled network is bipartite and that the state matrix of its (linear) dynamics satisfies a restrictive spectral condition. Finally, the design framework presented in [1,2] tackles both dynamical disturbance decoupling and distributed constraint satisfaction jointly, in a *single design phase*, whereas the formation control-oriented approaches in [31,28,6,10,27] require *restrictive or computationally demanding* procedures to reach their goals.

### 1.3 Scope of Work

In contrast to the techniques from [8,26,25,11], the MPC-based control laws discussed in this paper do not focus on network regulation in the presence of constraints, but rather on *closed-loop supervision* of both the plant and its control laws, which based upon the Network Realisation Function (NRF, see the companion paper [23], along with [21]) representation. In our framework, disturbance rejection and reference tracking are handled predominantly by the distributed (and dynamical) controllers discussed in [23]. By exploiting the closed-loop properties ensured by these NRF-based subcontrollers [20,21,24], the MPC-inspired procedures stated in the sequel are more akin to distributed versions of the so-called reference governors (see [7] for a broad discussion on these policies). Indeed, a careful inspection of the survey presented in [22] reveals that our framework does not fit neatly into any of the architectures addressed therein.

<footnote><sup>a</sup>Department of Automatic Control and Systems Engineering, "Politehnica" University of Bucharest, Bucharest, Romania</footnote>

<footnote><sup>b</sup>Laboratoire des Signaux et Systèmes, CentraleSupélec, Paris-Saclay University, Gif-sur-Yvette, France</footnote>

<footnote>Email addresses: andrei.sperila@upb.ro (Andrei Sperilă), alessio.iovine@centralesupelec.fr (Alessio Iovine), sorin.olaru@centralesupelec.fr (Sorin Olaru), patrick.panciatici@ieee.org (Patrick Panciatici).</footnote>
## 1.4 Contribution

From a system-theoretical perspective, the chief contributions of the work presented in this paper can be listed via the following three features:

a) ensuring command constraint satisfaction, in an MPC-based framework, for a plant in closed-loop configuration with a *dynamical controller*;

b) coordinating the control authority of our MPC-based policies with the NRF-based controllers from [21,24], while splitting the MPC's actuation efforts between the latter and the controlled network;

c) enabling constraint-compliant and NRF-based control, in which the MPC-based subcontrollers are designed to act only as a safeguard mechanism against constraint violation.

A further point to emphasise is the fact that all of the above features are achieved in a fully *distributed setting*, and we point out that the mechanism by which a dynamical controller's output is reliably constrained (for the purpose of command limitation) has not, to the best of our knowledge, been adequately addressed in system-theoretical literature, even in the *centralised setting*.

## 1.5 Paper Structure

The rest of the paper is structured as follows. Section 2 contains a set of preliminary notions, while Section 3 discusses our two-layer architecture and its overarching objectives, focusing on the MPC layer. Section 4 presents the means of obtaining tractable prediction models for an area-partitioned control scheme, and Section 5 handles the technicalities of constraint management in the distributed setting. Section 6 holds the main theoretical results and implementation procedures, while Section 7 showcases a numerical example based upon the vehicle platooning application discussed in [20] and Section 8 offers a set of concluding remarks. In addition to this, we include a pair of appendices: one focusing on the key concepts of NRF-based control laws (Appendix A), as discussed in our companion paper [23], and one which holds the proofs of this paper's main results (Appendix B).

We now move on to introducing the various preliminary notions which are necessary to state our main results.

# 2 Preliminaries

## 2.1 Notation and Definitions

Let $\mathbb{N}$, $\mathbb{R}$ and $\mathbb{C}$ denote, respectively, the set of natural, real and complex numbers. Let $\mathcal{M}^{p \times m}$ be the set of all $p \times m$ matrices whose entries are scalar elements that are part of a set denoted $\mathcal{M}$. Similarly, $\mathcal{M}^p$ denotes the set of vectors with dimension $p$ and entries in $\mathcal{M}$. For any matrix $M$, row$_i(M)$ is its $i^{th}$ row, col$_j(M)$ stands for its $j^{th}$ column and $\text{elm}_{ij}(M) := \text{row}_i(\text{col}_j(M))$. Moreover, the transpose of an arbitrary matrix $M$ is denoted $M^\top$. We denote by $\mathcal{R}(z)$ the set of real-rational functions of indeterminate $z \in \mathbb{C}$. A matrix with entries in $\mathcal{R}(z)$ is called a Transfer Function Matrix (TFM), and all such TFM's will be denoted using boldface letters.

The vector $e_i$ stands for the $i^{th}$ column belonging to the identity matrix of compatible dimensions. For any matrix $M \in \mathbb{R}^{p \times m}$ and any set $\mathcal{X} \subseteq \mathbb{R}^m$, we define $M\mathcal{X} := \{Mx | x \in \mathcal{X}\}$ and $(-\mathcal{X}) := \{-x | x \in \mathcal{X}\}$. In addition to this, for any two sets $\mathcal{X}_1 \subseteq \mathbb{R}^m$ and $\mathcal{X}_2 \subseteq \mathbb{R}^m$, we will henceforth denote their Minkowski sum as follows $\mathcal{X}_1 \oplus \mathcal{X}_2 := \{x_1 + x_2 | x_1 \in \mathcal{X}_1, x_2 \in \mathcal{X}_2\}$, while the operation $\mathcal{X}_1 \ominus \mathcal{X}_2 := \{y \in \mathbb{R}^m | \{y\} \oplus \mathcal{X}_2 \subseteq \mathcal{X}_1\}$ will represent the Pontryagin difference between the two sets.

## 2.2 System Representations

The class of systems considered in this paper are represented in the time domain by the set of equations

$$x[k+1] = Ax[k] + B_u u[k] + B_d d[k], \quad (1a)$$

$$y[k] = Cx[k] + D_u u[k] + D_d d[k], \quad (1b)$$

for any $k \in \mathbb{N}$ with $k \ge k_0 \in \mathbb{N}$, where $k_0$ represents the initial time instant. For the representations of type (1a)-(1b), which are referred to as *state-space realisations*, $u$ is the realisation's *controlled input vector*, $d$ its *disturbance vector*, $y$ its *output vector* and $x$ its *state vector*.

The type of network that our procedure aims to control is described by a more particular class of realisations than (1a)-(1b). These systems also satisfy the fact that

$$C = I, \quad D_u = O, \quad D_d = O, \quad (1c)$$

implying that the states are available<fnref n="1" /> for measurement.

**Remark 2.1** The disturbance vector $d[k]$ along with the noise affecting the state measurements (1b)-(1c) are key components of our framework. One of the main goals of the distributed MPC-based procedure proposed in this paper is to ensure constraint satisfaction for $x[k]$ and $u[k]$ in the presence of the aforementioned exogenous signals (along with additional ones, introduced in the sequel).

For systems of type (1a)-(1c), we denote by $x_c \in \mathbb{R}^{n_x}$ the initial condition of the state vector and we also consider $A \in \mathbb{R}^{n_x \times n_x}$, $B_u \in \mathbb{R}^{n_x \times n_u}$, $B_d \in \mathbb{R}^{n_x \times n_d}$, $C \in \mathbb{R}^{n_y \times n_x}$, while, more generically, we consider $D_u \in \mathbb{R}^{n_y \times n_u}$ and $D_d \in \mathbb{R}^{n_y \times n_d}$. The positive integer $n_x$ will be referred to as the order of the realisation, and the link between a system described by (1a)-(1b) and its TFM is given by

$$\mathbf{G}(z) = \begin{bmatrix} A - zI_{n_x} & B_u & B_d \\ C & D_u & D_d \end{bmatrix} := \\ := C(zI_{n_x} - A)^{-1} [B_u & B_d] + [D_u & D_d]. \quad (2)$$

Finally, given any $\mathbf{H} \in \mathcal{R}(z)^{n_y \times n_u}$ which can be described by (1a)-(1b) and (2), the shorthand notation $y[k] = \mathbf{H}(z) \star u[k]$ represents the time-response of $\mathbf{H}(z)$ to an input signal vector $u[k]$, which can be computed as

$$y[k] = \mathbf{H}(z) \star u[k] = D_H u[k] + \sum_{i=1}^{\infty} C_H A_H^{i-1} B_H u[k-i], \quad (3)$$

for any realisation $(A_H, B_H, C_H, D_H)$ of $\mathbf{H}(z)$.

<footnote n="1">¹ While being subject to measurement noise, as will be discussed (and suitably compensated for) in the sequel.</footnote>
<figure><img src="image_1.png" /><figcaption>Figure 1. High-level implementation scheme depicting the proposed distributed control strategy</figcaption></figure>

## 2.3 Distributed Networks

Although we consider that the vectors $x$ and $u$ are completely accessible, in terms of measurement and actuation, respectively, we do not assume that manipulating all of these variables is possible from any one location in our physical network. This naturally separates the network into a number of $N \in \mathbb{N}$ distinct areas, with $N > 1$.

To each of the network's $N$ areas, we now assign two sets of indices, which denote the entries of $x$ and $u$ that are (uniquely) accessible for measurement or actuation in that area. Without loss of generality, we assume that these indices are assigned to the sets in ascending order, since the rows and columns of the matrices expressed in (1a)-(1c) may be permuted at will. Therefore, each area will be represented by the pair

$$ \mathcal{A}_i := (\mathcal{A}_{xi}, \mathcal{A}_{ui}), \forall i \in \{1:N\}, \quad (4) $$

which collects the two aforementioned index sets.

To streamline the definition of the sets given in (4), we will provide generic expressions via the placeholder subscript $\bullet$, which stands for one of the subscripts $\{x, u\}$. Let each set contain a number of $n_{\bullet i} \in \mathbb{N}$ indices, where $n_{\bullet i} > 0, \forall i \in \{1:N\}$, such that we define

$$ \left\{ \begin{array}{l} \mathcal{A}_{\bullet i} := \{\alpha_{\bullet i} + j \mid j \in 1:n_{\bullet i}\}, \forall i \in \{1:N\}, \\ \alpha_{\bullet 1} := 0, \alpha_{\bullet \ell} := \alpha_{\bullet (\ell-1)} + n_{\bullet (\ell-1)}, \forall \ell \in \{2:N\}. \end{array} \right. \quad (5) $$

To each partition $^2$ in (5), we attach the selection matrix

$$ S_{\bullet i} := [e_{\alpha_{\bullet i+1}} \dots e_{\alpha_{\bullet i+n_{\bullet i}}}], \forall i \in \{1:N\}, \quad (6) $$

with which we define $x_i[k] := S_{xi}^T x[k]$, $x_{ci} := S_{xi}^T x_c$ and $u_i[k] := S_{ui}^T u[k]$. Once the network is partitioned, the design of our control scheme can begin.

# 3 The Control Architecture

## 3.1 Global Perspective and Functioning

The aim of this manuscript and of its companion paper [23] is to formalise the design framework of the two-layer control architecture for the interconnected dynamical system depicted in Figure 1. By identifying the subset of areas $\mathcal{N}_i \subseteq \{1:N\}$ which are able to transmit information to the network's $i^{th}$ area (and which trivially

satisfy $i \in \mathcal{N}_i$), for all $i \in \{1:N\}$, the communication graph is clearly designated and the fundamental operating principle of our control scheme is the following:

a) the state-space-implemented and NRF-based subcontrollers receive network state feedback and the command signals of other first-layer subcontrollers from their areas' neighbourhoods, using these values to compute the signals denoted $u_{fi}$ in Figure 1;

b) the MPC subcontrollers also receive network state information along with the state variables of the NRF subcontroller implementations, denoted $w_i$ in Figure 1, from their areas' neighbours to compute:

b1) the command signals denoted by $u_{s1i}$ in Figure 1, which are broadcast and combined additively with the state vectors $x_i$, before the later are fed to the NRF subcontrollers;

b2) the command signals denoted by $u_{s2i}$ in Figure 1, which remain local and are combined aditively with the NRF outputs $u_{fi}$, to obtain the control signals being sent to the actuators.

**Remark 3.1** In terms of control architecture design, it should be clear that both the first and the second layers employ feedback and feedforward components. However, their goals and their communication-based mechanisms are completely different. Where the first layer focuses on dynamical decoupling and disturbance rejection, the second layer concerns itself with constraint satisfaction and feasibility, in the sense of robust control invariance.

As previously mentioned in Remark 2.1, disturbance plays a key role in our chosen framework. Thus, we recall the disturbance vector $d[k]$ from (1a)-(1c) and we denote the measurement noise mentioned in said remark as $\zeta[k] \in \mathbb{R}^{n_x}$. Moreover, we also introduce $\beta_f[k] \in \mathbb{R}^{n_u}$, for the first layer's communication disturbance, along with $\beta_{s1}[k] \in \mathbb{R}^{n_x}$ and $\beta_{s2}[k] \in \mathbb{R}^{n_u}$, the communication disturbance associated with the second layer. The latter three signal vectors act as additive disturbance on the command signals depicted in Figure 1, namely $u_f[k] := [u_{f1}^T[k] \dots u_{fN}^T[k]]^\top$, $u_{s1}[k] := [u_{s11}^T[k] \dots u_{s1N}^T[k]]^\top$ and $u_{s2}[k] := [u_{s21}^T[k] \dots u_{s2N}^T[k]]^\top$, respectively. All of the exogenous signals mentioned in this paragraph may be combined into the global disturbance vector

$$ d_s[k] := [(\zeta[k] + \beta_{s1}[k])^\top \ \beta_{s2}^T[k] \ \beta_f^T[k] \ d^T[k]]^\top, \quad (7) $$

with respect to which we now proceed to formally define the objectives of our MPC layer.

<footnote n="2">2 As an example, a network with 12 states and 7 inputs could be partitioned as in (4)-(5) via the following collection of sets $\{(1:2), \{1\}), (3:5), \{2:3\}), (6:11), \{4:5\}), (12), \{6:7\})$</footnote>
$$
\begin{equation}
\begin{array}{@{}l@{}}
\displaystyle \min_{u_{s1i}[k], \dots, u_{s1i}[k+T_i-1], u_{s2i}[k], \dots, u_{s2i}[k+T_i+\bar{T}_i-1]} \sum_{t=1}^{T_i+\bar{T}_i} g_t(\xi_i[k+t], u_{s1i}[k+t-1], u_{s2i}[k+t-1]), \\
\text{subject to} \quad \left\{
\begin{array}{@{}l@{}}
\xi_i[k+t] = A_{si}\xi_i[k+t-1] + B_{s1i}u_{s1i}[k+t-1] + B_{s2i}u_{s2i}[k+t-1], \quad \forall t \in \{1: T_i + \bar{T}_i\}, \\
\xi_i[k+t] \in \Xi_{it}, \quad \forall t \in \{1: T_i\}, \quad \xi_i[k] = O, \\
u_{s1i}[k+t-1] \in \mathcal{U}_{s1i}, \quad u_{s2i}[k+t-1] \in \mathcal{U}_{s2i}, \quad \forall t \in \{1: T_i\}.
\end{array}
\right.
\end{array}
\tag{8}
\end{equation}
$$

3.2 Objectives of the MPC Layer

The second layer of the control architecture described in
Section 3.1 must be formed by $N$ separate MPC-based
subcontrollers, which accomplish the following goals:

I) Given two sets $\mathcal{X} \subseteq \mathbb{R}^{n_x}$ and $\mathcal{U} \subseteq \mathbb{R}^{n_u}$, the second layer's global control signals $u_{s1}[k]$ and $u_{s2}[k]$ robustly satisfy the constraints $x[k] \in \mathcal{X}$ and $u[k] \in \mathcal{U}$ in the presence of $d_s[k]$ from (7);

II) The area-based components of $u_{s1}[k]$ and $u_{s2}[k]$ are computed using *receding horizon-based optimization* (see, for example, Chapter 1 in [19]), with each area admitting tractable formulations of type (8), located at the top of this page, where:

a) $g_t(\cdot)$ is the stage cost at time instant $k+t$,

b) $\xi_i$ denotes the local state vector of the $i^{th}$ area's closed-loop system, to be addressed (along with its particular initialisation) in the sequel,

c) the scalar values $T_i \in \mathbb{N} \setminus \{0\}$ and $\bar{T}_i \in \mathbb{N}$ are, respectively, the (set-)constrained and the unconstrained prediction horizons of each area,

d) the sets denoted $\mathcal{U}_{s1i}$ along with $\mathcal{U}_{s2i}$ constrain the local second-layer commands and $\Xi_{it}$ is the set linking the local constraints with the global ones, to be constructed in the sequel;

III) For all $i \in \{1:N\}$, the $i^{th}$ area's MPC-based subcontroller from (8) computes $u_{s1i}[k]$ and $u_{s2i}[k]$ using *only information available in* $\mathcal{N}_i$;

IV) Ensuring the feasibility of each optimisation problem from (8) at some $k = k_0$ ensures all of them remain feasible for all $k > k_0$.

**Remark 3.2** As previously indicated in Remark 3.1, the goals stated above show that our second layer is only charged with enforcing constraints in a distributed manner. Thus, all objectives pertaining to reference tracking or robust stabilisation are to be handled by the NRF-based first layer, in order to take full advantage of its remarkable closed-loop guarantees (see, for example, the main results of [20,21,24]). In our framework, the MPC layer acts primarily as a safeguard for the NRF-based control laws, further highlighting the particular nature of our approach. As shown by the numerical example in the sequel, we aim to design our second layer so as to avoid acting on the first-layer-based closed-loop system as much as possible, only generating command signals of non-negligible amplitude when constraint violation becomes imminent.

In order to achieve all of the above-stated objectives,
however, tractable and area-based models will have to
be derived for the first-layer-based closed-loop system.

4 Obtaining an MPC-Oriented Formulation

This section is dedicated to discussing the derivation and
the manipulation of the models which form the equality
constraints of the problems stated in (8). Readers famil-
iar with the robust MPC framework (see Chapter 3 in
[19] for a comprehensive treatment of the topic) will un-
doubtedly recognise many of the formulations and con-
ventions adopted in the sequel. Nevertheless, the usage of
dynamical feedback instead of the classical, static gains
in the so-called tube-based MPC framework (see, for ex-
ample, Section 3.5 in [19]) renders our problem distinct
from other approaches in system-theoretical literature.
Therefore, given this distinct element of novelty, we pro-
ceed to lay out our expressions in an explicit manner.

4.1 Interfacing with the First Layer

We begin by briefly summarising several key concepts
from our companion paper [23]. Our architecture’s first
layer is implemented as depicted in Figure 2, located on
the next page, in which the distributed control laws are
represented by the TFM

$$
\mathbf{K}_{\mathrm{D}}(z) := \left[ \mathbf{K}_{\mathrm{D}1}^{\top}(z) \dots \mathbf{K}_{\mathrm{D}N}^{\top}(z) \right]^{\top}, \quad (9)
$$

where each $\mathbf{K}_{\mathrm{D}i}(z) = S_{ui}^{\top}\mathbf{K}_{\mathrm{D}}(z)$ represents the TFM of
the $i^{th}$ area's local subcontroller.

**Remark 4.1** We operate under the assumption that the communication infrastructure of our control architecture is endowed with error-detecting/correcting transmission mechanisms. Thus, all communication errors affecting transmitted information are modelled as disturbance which originates at the source of transmission, not at the receiving ends. Consequently, the resulting signal diagram of the first layer's closed-loop system will be precisely the one shown in Fig. 2, even when all $\mathbf{K}_{\mathrm{D}i}(z)$ are implemented separately, and the closed-loop mappings will have the beneficial properties discussed in [21,23,24].

Crucially, these subcontrollers are implemented via the
state-space realisations given in (10), located at the top
of the next page, whose component matrices (see Sec-
tion 4.1 in [23] for more details) are given, $\forall \ell \in \mathcal{A}_{ui}$, by

$$
\left\{
\begin{array}{l}
\tilde{A}_{r\ell} := [-a_{1\ell} \ \dots \ -a_{(n_{r\ell}-1)\ell}]^\top, \quad A_{r\ell} := \begin{bmatrix} \tilde{A}_{r\ell} & I_{n_{r\ell}-1} \\ -a_{n_{r\ell}\ell} & O \end{bmatrix}, \\
B_{r\ell} := [K_{1i}^\top \ \dots \ K_{n_{r\ell}i}^\top]^\top, \quad C_{r\ell} := [1 \ 0 \ \dots \ 0]^\top.
\end{array}
\right.
\tag{11}
$$

Henceforth, we define the orders $n_{wi} := \sum_{\ell=1}^{n_{ui}} n_{r\ell}$ and $n_w := \sum_{i=1}^{N} n_{wi}$. We also denote the state vectors of the
$$
\mathbf{K}_{\mathrm{Di}}(z) = \left[ \frac{A_{wi} - zI_{n_{wi}}}{C_{wi}} \middle| \frac{B_{wi}}{D_{wi}} \right] := \left[ \begin{array}{c|c} \begin{array}{c} A_{r(\alpha_{ui}+1)} - zI_{n_r(\alpha_{ui}+1)} \\ \vdots \\ A_{r(\alpha_{ui}+n_{ui})} - zI_{n_r(\alpha_{ui}+n_{ui})} \end{array} & \begin{array}{c} B_{r(\alpha_{ui}+1)} \\ \vdots \\ B_{r(\alpha_{ui}+n_{ui})} \end{array} \\ \hline \begin{array}{c} C_{r(\alpha_{ui}+1)} \\ \vdots \\ C_{r(\alpha_{ui}+n_{ui})} \end{array} & \begin{array}{c} O \\ \vdots \\ O \end{array} \end{array} \right] . \quad (10)
$$

<figure><img src="image_2.png" /><figcaption>Figure 2. Feedback loop of a network’s model G(z) with the
NRF-based distributed implementation K<sub>D</sub>(z)</figcaption></figure>

implementations from (10) as $w_i[k]$ and their initial conditions as $w_{ci} \in \mathbb{R}^{n_{wi}}$. Moreover, we concatenate the latter into $w_c := [w_{c1}^\top \dots w_{cN}^\top]^\top \in \mathbb{R}^{n_w}$, which are the initial conditions of the first layer’s global implementation.

One of the main features of this control scheme is that
it is able to *independently* compute the control signals
of the first layer as follows

$$
u_{fi}[k] := S_{ui}^{\top} u_f[k] = \mathbf{K}_{\mathrm{Di}}(z) \star \left[ \begin{matrix} u_f[k] + \beta_f[k] \\ x[k] + \zeta[k] + u_{s1}[k] + \beta_{s1}[k] \end{matrix} \right], \quad (12)
$$

for all $i \in \{1:N\}$, while employing only information
from each area's neighbourhood. Algebraically, this re-
duces to the fact that, for all $j \in \{1:N\} \setminus N_i$, the
implementations from (10)-(11) satisfy

$$
\mathrm{col}_q(B_{r\ell}) = O, \forall \ell \in \mathcal{A}_{ui}, q \in \mathcal{A}_{xj} \cup (\{n_x\} \oplus \mathcal{A}_{uj}). \quad (13)
$$

**Remark 4.2** Although we consider the first layer’s control laws to be obtained through an NRF-based framework, as in our companion paper [23], note that any set of distributed subcontrollers implemented akin to (9)-(13) may be employed for the procedures described in the sequel.

By using the signal vector introduced in (7), it is possible
to express the first layer’s closed-loop dynamics as

$$
\begin{bmatrix} x[k] \\ u_f[k] \end{bmatrix} = \left( \mathbf{F}_{\mathbf{Q}}(z) \begin{bmatrix} I_{(n_x+n_u)} \\ O \end{bmatrix} \right) \star \begin{bmatrix} u_{s1}[k] \\ u_{s2}[k] \end{bmatrix} + \mathbf{F}_{\mathbf{Q}}(z) \star d_s[k] + \mathcal{I}_{\mathbf{Q}}[k] \begin{bmatrix} x_c \\ w_c \end{bmatrix}, \quad (14)
$$

where $\mathbf{F}_{\mathbf{Q}}(z)$ is a strictly proper TFM, which is to say
that $\lim_{z \to \infty} \mathbf{F}_{\mathbf{Q}}(z) = O$. Notably, $\mathbf{F}_{\mathbf{Q}}(z)$ models the
closed-loop's response to exogenous inputs, while $\mathcal{I}_{\mathbf{Q}}[k]$ is a time-varying matrix which governs the closed-loop
system's evolution with respect to its initial conditions.

**Remark 4.3** While the signal vectors $d[k]$ and $\zeta[k]$ in (7) and (implicitly) in (14) possess classical interpretations in system-theoretical literature, the three vectors denoted generically as $\beta_{\bullet}[k]$ in (7) are caused primarily by encoding errors in the employed communication infrastructure. Due to finite-precision representation, the

absolute values of the signals in $\beta_{\bullet}[k]$ may be consid-
ered to be upper bounded and, unless stated otherwise,
the samples making up $\beta_{\bullet}[k]$ are implicitly assumed to be
uniformly distributed within their respective bounds. In
addition to this, we point out that boundedness (not nec-
essarily in terms of amplitude) also extends to both $d[k]$ and $\zeta[k]$, in most practical applications. Indeed, due to the stabilising properties of NRF-based control, the effect of all these signals upon the closed-loop system can be made bounded in some sense, by an appropriate choice of systemic norm optimisation (see Chapter 4 of [32] and Section 5.2 of [23]). Moreover, the results from [12] can be neatly incorporated into our considered approach, if one wishes to delve further into the stochastic setting.

The main result of our companion paper [23] shows that $\mathbf{F}_{\mathbf{Q}}(z)$ and $\mathcal{I}_{\mathbf{Q}}[k]$ are easily tunable via a freely chosen $\mathbf{Q}$-parameter (hence, the employed notation), thus enabling the *complete design* of the first layer’s closed-loop response. The full details pertaining to the synthesis and tuning of these closed-loop maps are given in our companion paper [23], yet, for the benefit of the reader, we have included a summary of the key properties ensured by NRF-based control in Appendix A.

4.2 Nominal Area-based Prediction Models

In order to formulate tractable prediction models for our
architecture's second layer, we need to account for the
state variables of the implementations from (10). Thus,
aside from the *x* and *u* subscripts appearing in (5)-(6),
we introduce the additional subscript *w* to refer to the
state vector of each area's NRF subcontroller, and we
extend the index collection defined originally in (4) to
$\mathcal{A}_i := (\mathcal{A}_{xi}, \mathcal{A}_{ui}, \mathcal{A}_{wi})$, $\forall i \in \{1:N\}$.

We begin our reformulation of the first layer’s closed-loop dynamics by pointing out the fact that, since the TFM $\mathbf{F}_{\mathbf{Q}}(z)$ from (14) is strictly proper, it is always possible to obtain state-space realisations of type

$$
\left[ \begin{array}{c} S_{\bar{x}_i}^{\top} \\ O \\ \bar{O}_{\bar{u}_i}^{\top} \end{array} \right] \mathbf{F}_{\mathbf{Q}}(z) \left[ \begin{array}{c} S_{\bar{x}_i} \\ O \\ \bar{O}_{\bar{u}_i} \end{array} \right] = \left[ \begin{array}{c} A_{\bar{s}_i} - zI_{n_{\bar{s}_i}} \\ B_{\bar{s}_i} \\ \bar{O}_{\bar{u}_i} \end{array} \right] \left[ \begin{array}{c} B_{\bar{s}_i} \\ O \\ \bar{O}_{\bar{u}_i} \end{array} \right], \quad (15)
$$

for all $i \in \{1:N\}$. By defining now $Z_i := \operatorname{diag}(S_{xi}, S_{ui})$ along with $Z_{ci} := \operatorname{diag}(S_{xi}, S_{wi})$, we employ the realisations from (15) to describe the MPC-oriented area model of the $i^{\text{th}}$ area, which can be expressed as

$$
\xi_i[k+1] = A_{si}\xi_i[k] + B_{s1i}u_{s1i}[k] + B_{s2i}u_{s2i}[k], \quad (16a)
$$

$$
x_i[k] = C_{xi} \xi_i[k] + (\psi_{xi}[k] + \theta_{xi}[k] + \delta_{xi}[k]), \quad (16b)
$$

$$
u_{fi}[k] = C_{ui}\xi_i[k] + (\psi_{ui}[k] + \theta_{ui}[k] + \delta_{ui}[k]), \quad (16c)
$$

$$
k \in \mathbb{N}, k \ge k_0, \xi_i[k_0] = O, \forall i \in \{1:N\}, (16d)
$$
where we have that

a) $\xi_i[k]$ is the state vector of the realisation obtained in (15) and which always has zero initial conditions at initial time $k_0$ (recall the identity from (3));

b) the contribution of exogenous disturbance to the area dynamics is given by

$$ egin{cases} \psi_i[k] := Z_i^\top \mathbf{F}_Q(z) \star d_s[k], \\ \psi_{xi}[k] := [I_{n_{xi}} \circ \psi_i[k], \psi_{ui}[k] := [\circ I_{n_{ui}}] \psi_i[k], \end{cases} \quad (17) $$

c) the contribution of the neighbouring areas' initial conditions to the area dynamics is given by

$$ egin{cases} \theta_i[k] := \sum_{j \in \mathcal{N}_i} Z_i^\top \mathcal{I}_Q[k] Z_{cj} \begin{bmatrix} x_{cj} \\ w_{cj} \end{bmatrix}, \\ \theta_{xi}[k] := [I_{n_{xi}} \circ \theta_i[k], \theta_{ui}[k] := [\circ I_{n_{ui}}] \theta_i[k], \end{cases} \quad (18) $$

d) the residual contributions of initial conditions from outside of the $i^{th}$ area's neighbourhood and of cross-coupling with other areas' MPC subcontrollers is

$$ egin{cases} \delta_i[k] := \sum_{j \in \{1:N\} \setminus \{i\}} Z_i^\top \mathbf{F}_Q(z) \begin{bmatrix} Z_j \\ O \end{bmatrix} \star \begin{bmatrix} u_{s1j}[k] \\ u_{s2j}[k] \end{bmatrix} + \\ \qquad + \sum_{j \in \{1:N\} \setminus \mathcal{N}_i} Z_i^\top \mathcal{I}_Q[k] Z_{cj} \begin{bmatrix} x_{cj} \\ w_{cj} \end{bmatrix}, \\ \delta_{xi}[k] := [I_{n_{xi}} \circ \delta_i[k], \delta_{ui}[k] := [\circ I_{n_{ui}}] \delta_i[k], \end{cases} \quad (19) $$

The signals from (17)-(19) separate naturally into two classes, in the form of information that is:

i) available to a second-layer subcontroller, represented by $\theta_{xi}[k]$ and $\theta_{ui}[k]$;

ii) unavailable to a second-layer subcontroller, represented by $(\psi_{xi}[k] + \delta_{xi}[k])$ and $(\psi_{ui}[k] + \delta_{ui}[k])$.

**Remark 4.4** The quantities described in (19) are called residual because it is the primary goal of the first layer to attenuate the effect that these signals have on the area dynamics from (16a)-(16d). Since these terms, along with those in (17), represent the contribution of variables which are not available (in terms of measurement or prediction) to the $i^{th}$ area, the only sensible recourse is to minimise their influence (see Section 5 in [23]).

Thus, each subcontroller must employ the information that is encoded in the signals from (18) to compute a pair of commands $u_{s1i}[k]$ and $u_{s2i}[k]$ which robustly account for the effects of the signals from (17) and (19). Yet, before translating this objective into a set-theoretical formulation, we investigate the effect of unreliable measurements for $x_c$ and $w_c$ in the distributed setting.

### 4.3 Accounting for Inaccurate Initial Conditions

Given some pre-specified set $V \subseteq \mathbb{R}^{n_x+n_w}$, we consider the additively disturbed initial conditions defined as

$$ \begin{bmatrix} \tilde{x}_c \\ \tilde{w}_c \end{bmatrix} := \begin{bmatrix} x_c + \nu_x \\ w_c + \nu_w \end{bmatrix}, \quad \begin{bmatrix} \tilde{x}_{ci} \\ \tilde{w}_{ci} \end{bmatrix} := \begin{bmatrix} x_{ci} + S_{x_i}^\top \nu_x \\ w_{ci} + S_{w_i}^\top \nu_w \end{bmatrix}, \quad (20) $$

for any index $i \in \{1:N\}$, along with any $[\nu_x^\top \nu_w^\top]^\top \in V$.

**Remark 4.5** In light of the arguments made in Remark 4.1, notice that the type of disturbance discussed in this section can always be modelled as in (20). Practical causes for these unwanted terms include noisy measurement of the network's state variables and faulty reading of the NRF-based controller's state vector.

The main consequence of taking into account these terms is that the prediction models from (16b)-(16c) become

$$ \tilde{x}_i[k] = C_{xi}\xi_i[k] + (\psi_{xi}[k] + \tilde{\theta}_{xi}[k] + \tilde{\delta}_{xi}[k]), \quad (21a) $$

$$ \tilde{u}_{fi}[k] = C_{ui}\xi_i[k] + (\psi_{ui}[k] + \tilde{\theta}_{ui}[k] + \tilde{\delta}_{ui}[k]), \quad (21b) $$

when considering the noise-affected initial conditions introduced in (20) instead of the real ones, where

$$ \left\{ \begin{aligned} \tilde{\theta}_i[k] &:= \theta_i[k] + \sum_{j \in \mathcal{N}_i} Z_i^\top \mathcal{I}_Q[k] Z_{cj} \begin{bmatrix} S_{x_j}^\top \nu_x \\ S_{w_j}^\top \nu_w \end{bmatrix}, \\ \tilde{\theta}_{xi}[k] &:= [I_{n_{xi}} \circ \tilde{\theta}_i[k], \tilde{\theta}_{ui}[k] := [\circ I_{n_{ui}}] \tilde{\theta}_i[k], \end{aligned} \right. \quad (22) $$

along with

$$ \left\{ \begin{aligned} \tilde{\delta}_i[k] &:= \delta_i[k] + \sum_{j \in \{1:N\} \setminus \mathcal{N}_i} Z_i^\top \mathcal{I}_Q[k] Z_{cj} \begin{bmatrix} S_{x_j}^\top \nu_x \\ S_{w_j}^\top \nu_w \end{bmatrix}, \\ \tilde{\delta}_{xi}[k] &:= [I_{n_{xi}} \circ \tilde{\delta}_i[k], \tilde{\delta}_{ui}[k] := [\circ I_{n_{ui}}] \tilde{\delta}_i[k], \end{aligned} \right. \quad (23) $$

represent the noise-affected versions of the signals from (18) and (19), respectively, for all $i \in \{1:N\}$.

### 4.4 Connecting the Augmented Prediction Models to the Nominal Ones

Since the vector $[\nu_x^\top \nu_w^\top]^\top$ cannot be determined a priori, our MPC-based subcontrollers will have to make use of the initial conditions in (20) instead of the real ones. Consequently, our control policies will employ the prediction models affected by additive disturbance (21a)-(21b), instead of the ones given in (16b)-(16c). However, a straightforward computation yields the fact that the original area dynamics, given by the identities obtained in (16b)-(16c), can be retrieved from (21a)-(21b) via

$$ x_i[k] = \tilde{x}_i[k] - \eta_{xi}[k], \quad u_{fi}[k] = \tilde{u}_{fi}[k] - \eta_{ui}[k], \quad (24) $$

for all $i \in \{1:N\}$, where $\eta_{xi}[k]$ and $\eta_{ui}[k]$ are defined as

$$ \left\{ \begin{aligned} \eta_i[k] &:= \sum_{j=1}^{N} Z_i^\top \mathcal{I}_Q[k] Z_{cj} \begin{bmatrix} S_{x_j}^\top \nu_x \\ S_{w_j}^\top \nu_w \end{bmatrix}, \\ \eta_{xi}[k] &:= [I_{n_{xi}} \circ \eta_i[k], \eta_{ui}[k] := [\circ I_{n_{ui}}] \eta_i[k]. \end{aligned} \right. \quad (25) $$

Having obtained the area-based prediction models in (16a)-(16d) and having augmented them as in (21a)-(21b), we now proceed to discuss constraint formulation for the optimisation problems described in (8).

## 5 Constraints in the Distributed Setting

The most important step in breaking down the second layer's distributed synthesis problem into locally solvable subproblems is to assign independent control objectives to each of the network's $N$ areas. Recalling the sets $\mathcal{X}$ and $\mathcal{U}$, introduced in Section 3.2, we begin by defining the pairs of sets $(\mathcal{X}_i, \mathcal{U}_i)$, with $i \in \{1:N\}$, which satisfy

$$ \left\{ \begin{array}{l} x[k] \in \mathcal{X} \iff x_i[k] \in \mathcal{X}_i, \forall i \in \{1:N\}, \\ u[k] \in \mathcal{U} \iff u_i[k] \in \mathcal{U}_i, \forall i \in \{1:N\}. \end{array} \right. \quad (26) $$
**Remark 5.1** The assumption that the pair $(\mathcal{X}, \mathcal{U})$ can be decomposed as in (26) is not only required, when attempting to break down the global problem into any configuration of local subproblems, but it is also a mild one (see its equivalent counterpart in [1,2]). Indeed, the conditions located on the left-hand side of (26) are often represented (in practice) by box constraints, for which said decomposition may always be performed.

The chief aim of the subcontrollers which form the second layer (recall (8) and Figure 1) is to ensure the global, left-hand side inclusions from (26) by independently satisfying the local ones, on the right-hand side of (26). With respect to the closed-loop dynamics from (14), this problem can be broken down into five components:

i) accounting (in a set-theoretical framework) for the effect of the exogenous disturbance from (7);

ii) translating the command-related objectives from (26) into first-layer command constraints;

iii) linking command constraints to the state variables of the implementations from (10);

iv) accounting for the inaccurate measurements modelled in (20);

v) accounting for area cross-coupling.

We now proceed to tackle each of these subproblems, in the order in which they are indicated above.

## 5.1 Placing Constraints on Exogenous Disturbance

We begin by addressing the components of the signal vector from (7). Recalling the measurement noise set from (20) and the arguments made in Remark 4.5, we begin by considering that $\zeta[k] \in [I_n \circ \mathcal{V}]$. We contain the rest of the exogenous signals from (7) by introducing

$$ \begin{cases} \mathcal{D}_d \subseteq \mathbb{R}^{n_d} \text{ s.t. } d[k] \in \mathcal{D}_d, & \mathcal{D}_{\beta_f} \subseteq \mathbb{R}^{n_u} \text{ s.t. } \beta_f[k] \in \mathcal{D}_{\beta_f}, \\ \mathcal{D}_{\beta_{s1}} \subseteq \mathbb{R}^{n_x} \text{ s.t. } \beta_{s1}[k] \in \mathcal{D}_{\beta_{s1}}, \mathcal{D}_{\beta_{s2}} \subseteq \mathbb{R}^{n_x} \text{ s.t. } \beta_{s2}[k] \in \mathcal{D}_{\beta_{s2}}. \end{cases} \tag{27} $$

By doing so, it becomes possible to construct the set $\mathcal{D}_s \subseteq \mathbb{R}^{2n_u+n_x+n_d}$ which satisfies $d_s[k] \in \mathcal{D}_s$ and which, in tandem with the expressions from (17), enables us to compute the sets containing the effects of $d_s[k]$, with $k \ge k_0$, on each area's dynamics, henceforth denoted as

$$ \left\{ \begin{array}{l} \Psi_{xit} \subseteq \mathbb{R}^{n_{xi}} \text{ s.t. } \psi_{xi}[k_0 + t] \in \Psi_{xit}, \forall d_s[k] \in \mathcal{D}_s, \\ \Psi_{uit} \subseteq \mathbb{R}^{n_{ui}} \text{ s.t. } \psi_{ui}[k_0 + t] \in \Psi_{uit}, \forall d_s[k] \in \mathcal{D}_s. \end{array} \right. \tag{28} $$

**Remark 5.2** One of the key benefits of our two-layer approach is the fact that, from the point of view of the NRF-based layer, the signals $u_{s1}[k]$ and $u_{s2}[k]$ are perceived in the same manner as the all the other components the vector from (7), namely as additive and exogenous disturbance. Consequently, the special disturbance attenuation properties (such as string stability, see [20] and [9], along with the references therein) of our control system are preserved when adding the second layer on top of the first. Note that such special properties are particularly challenging to ensure using exclusively MPC-based techniques, a topic which has generated a great deal of research in system-theoretical literature [31,28,6,10,27].

One of the main contributions of the present work is the fact that it combines the native benefits of two separate techniques, disturbance attenuation for NRF and constraint management for MPC, into a single framework.

Moving onward, we address the variables associated with the first-layer subcontrollers and their implementations.

## 5.2 Placing Constraints on First Layer Commands

A cursory inspection of Figure 2 reveals the identity

$$ u_i[k] = u_{fi}[k] + u_{s2i}[k] + S_{ui}^{\top} \beta_{s2}[k], \forall i \in \{1:N\}. \tag{29} $$

For the state dynamics obtained in Section 4, the problem of constraining the command signals sent to the network reduces to judiciously constraining the control effort of both the first and the second layer. Thus, for each area, we choose the pair of sets introduced in (8), namely

$$ \mathcal{U}_{s1i} \subseteq \mathbb{R}^{n_{xi}} \text{ s.t. } u_{s1i}[k] \in \mathcal{U}_{s1i}, \forall i \in \{1:N\}, \tag{30a} $$

$$ \mathcal{U}_{s2i} \subseteq \mathbb{R}^{n_{ui}} \text{ s.t. } u_{s2i}[k] \in \mathcal{U}_{s2i}, \forall i \in \{1:N\}. \tag{30b} $$

and we employ the ones from (30b) in order to compute

$$ \mathcal{U}_{fi} := \mathcal{U}_i \ominus (\mathcal{U}_{s2i} \oplus (S_{ui}^{\top} \mathcal{D}_{\beta_{s2}})), \forall i \in \{1:N\}. \tag{31} $$

Notice that we have $u_{fi}[k] \in \mathcal{U}_{fi} \implies u_i[k] \in \mathcal{U}_i$, which allows us to translate the objectives given in (26) to the output vectors of our area dynamics from (16a)-(16d).

**Remark 5.3** The freely chosen sets given in (30a)-(30b) are the main design parameters of the second layer. They can be assimilated with the extent to which the second layer is permitted to interfere in the dynamics of the first layer's closed-loop system, so as to ensure that (26) is satisfied. The manner in which these sets may be chosen, and the resulting trade-off, will be addressed in the sequel.

## 5.3 Placing Constraints on First-Layer States

Since the initial conditions of the first layer's state-space implementation explicitly impact the signals located on the left-hand side of (14), it follows that the state vectors of the NRF-based subcontrollers from (10) should also possess a state-based description, which will be used to ensure the satisfaction of the inclusions in (26).

**Remark 5.4** The results presented in this section represent one of the key theoretical novelties introduced in this paper. In the classical tube MPC setting (see, for example, [8,26,11]), the command constraints on the signals analogous to $u_{fi}[k]$ are straightforward to compute by means of constant matrices. In our framework, however, it becomes necessary to pass through the dynamics in (10), which significantly complicates the resulting expressions. Given the delicate nature of this undertaking, said expressions are presented, within the current section, in as explicit a manner as possible, regardless of the apparent similarity with standard tube MPC (recall the first paragraph of Section 4). Nevertheless, the benefits of NRF-based control (as stated in Remark 5.2) more than make up for the increase in complexity when formulating the synthesis procedure (and not, we point out, when implementing the control laws). Indeed, a further advantage of our approach is the fact that, instead of limiting ourselves to block-diagonal and constant stabilising feedback
$$
\begin{equation}
\begin{aligned}
& \mathcal{D}_{wi} \subseteq \mathbb{R}^{n_{ui}+n_{xi}} \quad \text{s.t.} \quad \left[ \begin{matrix} u_{fi}[k] + S_{ui}^{\top} \beta_{f}[k] \\ x_i[k] + S_{xi}^{\top} \zeta[k] + u_{s1i}[k] + S_{xi}^{\top} \beta_{s1}[k] \end{matrix} \right] \in \mathcal{D}_{wi}, \\
& \forall u_{fi}[k] \in \mathcal{U}_{fi}, (S_{ui}^{\top} \beta_{f}[k]) \in S_{ui}^{\top} \mathcal{D}_{\beta_{f}}, x_i[k] \in \mathcal{X}_i, (S_{xi}^{\top} \zeta[k]) \in [s_{xi}^{\top} o] \mathcal{V}, u_{s1i}[k] \in \mathcal{U}_{s1i}, (S_{xi}^{\top} \beta_{s1}[k]) \in S_{ui}^{\top} \mathcal{D}_{\beta_{s1}}.
\end{aligned}
\tag{33}
\end{equation}
$$

$$
\mathcal{W}_{\ell n_{r\ell}} := (-a_{n_{r\ell}\ell}) \mathcal{W}_{\ell 1} \oplus \left( \bigoplus_{p \in \mathcal{N}_i} (K_{n_{r\ell}\ell} Z_{rp} Z_{rp}^{\top}) \mathcal{D}_{wp} \right), \forall \ell \in \mathcal{A}_{ui}, \quad (34a)
$$

$$
\mathcal{W}_{\ell j} := (-a_{j\ell}) \mathcal{W}_{\ell 1} \oplus \mathcal{W}_{\ell(j+1)} \oplus \left( \bigoplus_{p \in \mathcal{N}_i} (K_{j\ell} Z_{rp} Z_{rp}^{\top}) \mathcal{D}_{wp} \right), \forall \ell \in \mathcal{A}_{ui}, j \in \{2:n_{r\ell}-1\}. \quad (34b)
$$

(see Assumption 1 in [8,26,11]), we leverage a greater degree of freedom in considering dynamical stabilising feedback (10) endowed with generic sparsity patterns (11).

We begin by inspecting the distinctive structure of the matrices which make up the realisations from (10). The latter enables us to naturally partition each of the state vectors belonging to the first-layer subcontrollers into $w_i[k] = [w_{r(\alpha_{ui}+1)}^\top[k] \dots w_{r(\alpha_{ui}+n_{ui})}^\top[k]]$, in which the vectors $w_{r\ell}[k]$ have $n_{r\ell}$ rows for all $\ell \in \{1:n_u\}$, i.e., the same as the matrices $A_{r\ell}$ and $B_{r\ell}$ from (11). Recall- ing (11), the key observation which underpins the entire mathematical mechanism presented in this section is the fact that $\text{row}_\ell(u_f[k]) = C_{r\ell}w_{r\ell}[k] = \text{row}_1(w_{r\ell}[k])$, for all $\ell \in \{1:n_u\}$. Thus, our current objective becomes apparent: we aim to employ the bounds placed upon the partitions of $u_f[k]$ to satisfy further set inclusions with respect to $w_{r\ell}[k]$, in a distributed manner.

In order to state the main result of this section, we introduce the sets $\mathcal{W}_{\ell 1} \subseteq \mathbb{R}$, with $\ell \in \{1:n_u\}$, which satisfy

$$
\text{row}_\ell(u_f[k]) \in \mathcal{W}_{\ell 1}, \quad \forall \ell \in \mathcal{A}_{ui} \implies u_{fi}[k] \in \mathcal{U}_{fi}, \quad (32)
$$

for all $i \in \{1:N\}$. The following result employs the newly introduced sets with the aim of constraining the state variables of the first-layer subcontrollers.

**Proposition 5.1** Let the first-layer subcontrollers be implemented as in (9)-(13). For each $i \in \{1:N\}$, define $Z_{ri} := \text{diag}(S_{ui}, S_{xi})$ and the sets from (33) and (34a)-(34b), located at the top of this page. Assume that:

A1) the inclusions given in (27), (30a)-(30b) and (33) hold for $k = k_0$,

A2) $\text{row}_j(w_{r\ell}[k_0]) \in \mathcal{W}_{\ell j}, \forall \ell \in \mathcal{A}_{ui}, j \in \{1:n_{r\ell}\}$,

A3) $\text{row}_\ell(u_f[k_0+1]) \in \mathcal{W}_{\ell 1}, \forall \ell \in \mathcal{A}_{ui}, i \in \{1:N\}$,

given some $k_0 \in \mathbb{N}$. Then, $\forall i \in \{1:N\}$, we have that:

i) $\text{row}_j(w_{r\ell}[k_0+1]) \in \mathcal{W}_{\ell 1}, \forall \ell \in \mathcal{A}_{ui}, j \in \{1:n_{r\ell}\}$;

ii) $u_i[k_0] \in \mathcal{U}_i$.

**PROOF.** See Appendix B.

**Remark 5.5** Note that the $i^{th}$ area's sets given in (34a)-(34b) are computed efficiently, using only for those indices belonging to $\mathcal{N}_i$. This is owed to the sparsity structure displayed in (13), which serves to incentivise the synthesis of sparse control laws for the first layer (see Appendix A, along with [23]). Note that the bounds associated with these sets could also be tightened by adapting the ideas presented in [26], through the sequential reduction of those $\mathcal{X}_j, \mathcal{U}_j, \mathcal{U}_{s1j}$ and $\mathcal{U}_{s2j}$ for which $j \in \mathcal{N}_i$.

The result formalised via Proposition 5.1 will prove instrumental in guaranteeing the recursive feasibility of the second layer and, in retrospect, the choice to impose constraints upon $u_{fi}[k]$, instead of the $u_i[k]$ vectors in (26), becomes justified. Not only do the specialised implementations presented in (10) allow us to construct set-based descriptions (34a)-(34b) for the states belonging to the first-layer implementations, but they also enable us to explicitly link the sets mentioned above to the second layer's design objectives. By eliminating the need to individually select each of these sets (which is an intricate and highly empirical problem, in the distributed setting), the design procedure becomes greatly simplified.

**Remark 5.6** By inspecting Figure 1 along with the various sets defined in this section, the intended use of $u_{s1i}[k]$ and $u_{s2i}[k]$ becomes clear. Indeed, just as $u_{s2i}[k]$ is primarily meant to ensure that $x_i[k] \in \mathcal{X}_i$, so too is $u_{s1i}[k]$ chiefly responsible for ensuring that $\text{row}_j(w_{r\ell}[k]) \in \mathcal{W}_{\ell j}$ and, implicitly, that $u_i[k] \in \mathcal{U}_i$ (recall Proposition 5.1).

## 5.4 Addressing Uncertainty in the Initial Conditions

Having obtained the means by which we will satisfy the inclusions given in (26), we now turn our attention to the disturbance discussed in Sections 4.3 and 4.4. In order to account for the inaccuracies modelled in (20), we begin by exploiting the expressions from (25) to compute

$$
\left\{
\begin{array}{l}
\mathcal{H}_{xit} \subseteq \mathbb{R}^{n_{xi}} \text{ s.t. } \eta_{xi}[k_0 + t] \in \mathcal{H}_{xit}, \forall [\nu_x^\top \nu_w^\top]^\top \in \mathcal{V}, \\
\mathcal{H}_{uit} \subseteq \mathbb{R}^{n_{ui}} \text{ s.t. } \eta_{ui}[k_0 + t] \in \mathcal{H}_{uit}, \forall [\nu_x^\top \nu_w^\top]^\top \in \mathcal{V},
\end{array}
\right.
\tag{35}
$$

which allows us to define, for all $i \in \{1:N\}$, the sets

$$
\tilde{\mathcal{X}}_{it} := \mathcal{X}_i \ominus (-\mathcal{H}_{xit}), \quad (36)
$$

and

$$
\tilde{\mathcal{W}}_{\ell 1t} := \mathcal{W}_{\ell 1} \ominus (-e_{(\ell-\alpha_{ui})}^{\top} \mathcal{H}_{uit}), \quad \forall \ell \in \mathcal{A}_{ui}. \quad (37)
$$

Thus, our problem reduces to employing the dynamics from (21a)-(21b) with the aim of constraining $\tilde{x}_i[k]$ and $\tilde{u}_{fi}[k]$ to the sets obtained in (36) and (37), respectively. Recalling the identities from (24) directly yields the desired inclusions from (26). All that remains is to compute bounding sets for the signals from (22)-(23).

## 5.5 Placing Constraints on Area Cross-coupling

In order to provide compact formulations in the sequel, we first construct the sets $\mathcal{W}_i \subseteq \mathbb{R}^{n_{wi}}$ which satisfy

$$
\text{row}_j(w_{r\ell}[k]) \in \mathcal{W}_{\ell j}, \forall j \in \{1:n_{r\ell}\}, \ell \in \mathcal{A}_{ui} \Rightarrow w_i[k] \in \mathcal{W}_i, \quad (38)
$$

for all $i \in \{1:N\}$, by simply taking the Cartesian product of the $\mathcal{W}_{\ell j}$ sets appearing in (38). Since the aim of
$$
\Xi_{it} := \left\{ \phi \in \mathbb{R}^{n_{si}} \left| \begin{bmatrix} H_{1it}C_{xi} \\ H_{2(\alpha_{ui}+1)t}\mathrm{row}_1(C_{ui}) \\ \vdots \\ H_{2(\alpha_{ui}+n_{ui})t}\mathrm{row}_{n_{ui}}(C_{ui}) \end{bmatrix} \right. \phi \le \left. \begin{bmatrix} h_{1it} - H_{1it}\tilde{\theta}_{xi}[k_0 + t] \\ h_{2(\alpha_{ui}+1)t} - H_{2(\alpha_{ui}+1)t}\mathrm{row}_1(\tilde{\theta}_{ui}[k_0 + t]) \\ \vdots \\ h_{2(\alpha_{ui}+n_{ui})t} - H_{2(\alpha_{ui}+n_{ui})t}\mathrm{row}_{n_{ui}}(\tilde{\theta}_{ui}[k_0 + t]) \end{bmatrix} \right\}, \forall t \in \{1: T_i\}. \tag{44}
$$

the second layer is to constrain the state variables of the
network and those of the implementations given in (10)
to the sets $\mathcal{X}_i$ and $\mathcal{W}_i$, respectively, we consider that the
initial conditions of the first layer's closed-loop system
satisfy $x_{ci} \in \mathcal{X}_i$ along with $w_{ci} \in \mathcal{W}_i$, for all $i \in \{1:N\}$.
In this case, the definitions introduced in (20) yield

$$
\tilde{x}_{ci} = x_{ci} + S_{xi}^{\top} \nu_x \in \tilde{\mathcal{X}}_{ci} := \mathcal{X}_i \oplus [s_{xi}^{\top} o] \mathcal{V}, \quad (39a)
$$

$$
\tilde{w}_{ci} = w_{ci} + S_{wi}^{\top} \nu_w \in \tilde{W}_{ci} := W_i \oplus [o S_{wi}^{\top}] V, \quad (39b)
$$

for all $i \in \{1:N\}$, which allows us to construct

$$
\left\{
\begin{array}{l}
\tilde{\Theta}_{xit} \subseteq \mathbb{R}^{n_{xi}} \text{ s.t. } \tilde{\theta}_{xi}[k_0 + t] \in \tilde{\Theta}_{xit}, \\
\quad \forall \tilde{x}_{cj} \in \tilde{\mathcal{X}}_{cj}, \tilde{w}_{cj} \in \tilde{W}_{cj}, j \in \mathcal{N}_i, \\
\tilde{\Theta}_{uit} \subseteq \mathbb{R}^{n_{ui}} \text{ s.t. } \tilde{\theta}_{ui}[k_0 + t] \in \tilde{\Theta}_{uit}, \\
\quad \forall \tilde{x}_{cj} \in \tilde{\mathcal{X}}_{cj}, \tilde{w}_{cj} \in \tilde{W}_{cj}, j \in \mathcal{N}_i, \\
\end{array}
\right.
\tag{40}
$$

along with

$$
\left\{
\begin{array}{l}
\tilde{\Delta}_{xit} \subseteq \mathbb{R}^{n_{xi}} \text{ s.t. } \tilde{\delta}_{xi}[k_0 + t] \in \tilde{\Delta}_{xit}, \\
\quad \forall u_{s1j}[k] \in \mathcal{U}_{s1j}, u_{s2j}[k] \in \mathcal{U}_{s2j}, j \in \{1:N\} \setminus \{i\}, \\
\quad \forall \tilde{x}_{cl} \in \tilde{\mathcal{X}}_{cl}, \tilde{w}_{cl} \in \tilde{\mathcal{W}}_{cl}, j \in \{1:N\} \setminus \mathcal{N}_i, \\
\tilde{\Delta}_{uit} \subseteq \mathbb{R}^{n_{ui}} \text{ s.t. } \tilde{\delta}_{ui}[k_0 + t] \in \tilde{\Delta}_{uit}, \\
\quad \forall u_{s1j}[k] \in \mathcal{U}_{s1j}, u_{s2j}[k] \in \mathcal{U}_{s2j}, j \in \{1:N\} \setminus \{i\}, \\
\quad \forall \tilde{x}_{cl} \in \tilde{\mathcal{X}}_{cl}, \tilde{w}_{cl} \in \tilde{\mathcal{W}}_{cl}, l \in \{1:N\} \setminus \mathcal{N}_i.
\end{array}
\right.
\tag{41}
$$

**Remark 5.7** Note that, in choosing the sets from (30b), one negotiates a command-budgeting balance (as limited by $\mathcal{U}_i$) between the first and second layers. This induces a trade-off associated with selecting the sets from (30a)-(30b), since the more freedom to act $u_{s1i}[k]$ is given, via the set from (30a), its effect will expand throughout the sets from (33) and (34a)-(34b). Similarly, if the network's areas are not all perfectly decoupled (recall Remark 4.4), then too much leeway given to $u_{s1i}[k]$ and $u_{s2i}[k]$ may result in significant area cross-coupling. In this case, one may derive an extended version of the technique from [25] (which manifests its procedure's effectiveness on only two collaborating areas) to mitigate these particular aspects. Section 4 in [25] also provides valuable insight into negotiating this particular type of trade-off.

With all of these sets now defined according to their cor-
responding areas, we now move on to the following sec-
tion, in which we present guarantees for the (recursive)
satisfaction of the second layer's constraints.

**6 Distributed Implementation of the MPC Layer**

In this section, we proceed to show how to ensure global
theoretical guarantees for the first layer's closed-loop

system, by *independently* implementing an MPC con-
troller for each area’s dynamics described in (16a)-(16d).
The results presented here are inspired by the core idea
investigated in [17], and our proposed technique repre-
sents the MPC-oriented analogue of the indicated work.

**Remark 6.1** As discussed in this manuscript’s companion paper [23], our approach showcases a series of notable similarities with the SLS-based procedure from [1,2]. Indeed, in both instances, the focus rests upon designing the closed-loop response of the network in accordance with a number of structural restrictions. Additionally, the d-locality constraints from [1,2] closely match our own area-based sparsity patterns discussed in Remark 4.4, while the feasibilities of these two approaches are also closely intertwined, since both the SLS [29] and the Youla Parametrisation (upon which the first layer is based, see [21]) yield the set of all achievable and stable closed-loop maps (see [30]). On the other hand, our procedure leverages the fact that the NRF-based layer can always be designed and implemented offline, thus relegating all online computation to the significantly less resource-demanding MPC layer, as shown in the sequel. This stands in stark contrast to the approach from [1,2], where sparsity constraints must be factored into the online computation, alongside all the MPC-like constraints.

6.1 Distributed Enforcement of Constraint Satisfaction

We begin by addressing the final and most important
component of the formulations introduced in (8). Recall
that, in stating the objectives of our control architec-
ture's second layer, the optimisation problems from (8)
make use of the sets Ξ<sub>it</sub> to ensure the satisfaction of the
constraints pertaining to x[k] and u[k]. The aforemen-
tioned sets, whose precise expressions have been deferred
up to this point, may now finally be stated in a formal
manner. This is done by first computing the sets

$$
\left\{
\begin{array}{l}
\mathcal{P}_{1it} := \{v \in \mathbb{R}^{n_{xi}} \mid H_{1it}v \le h_{1it}\}, \\
\mathcal{P}_{1it} \subseteq \tilde{\mathcal{X}}_{it} \ominus (\Psi_{xit} \oplus \tilde{\Delta}_{xit}),
\end{array}
\right.
\tag{42}
$$

and

$$
\left\{
\begin{array}{l}
\mathcal{P}_{2\ell t} := \{v \in \mathbb{R} \mid H_{2\ell t}v \le h_{2\ell t}\}, \\
\mathcal{P}_{2\ell t} \subseteq \tilde{W}_{\ell t} \ominus ((e_{(\ell-\alpha_{ui})}^{\top} \Psi_{uit}) \oplus (e_{(\ell-\alpha_{ui})}^{\top} \tilde{\Delta}_{uit})),
\end{array}
\right.
\quad (43)
$$

for every index $\ell \in \mathcal{A}_{ui}$ and any set of real-valued tuples $(H_{1it}, h_{1it}) \in \mathbb{R}^{q_{1it} \times n_{xi}} \times \mathbb{R}^{q_{1it}}$ and $(H_{2\ell t}, h_{2\ell t}) \in \mathbb{R}^{q_{2\ell t}} \times \mathbb{R}$ which satisfy (42)-(43). These tuples allow us to define $\Xi_{it}$ as in (44), located at the top of this page.

**Remark 6.2** We implicitly assume that the first-layer closed-loop system has been designed in such a way as to
ensure that the sets from (42)-(44) are non-empty. If this is not the case, then the reason is the fact that the first layer has not been appropriately tuned in order to tackle cross-coupling and the exogenous disturbance. Provided that the sets $\mathcal{X}$ and $\mathcal{U}$ (along with their area-based compositions) are available during the first layer's design phase, a judicious choice of systemic norm optimisation in Section 5.2 of [23] (tied to the nature of these sets) can be used to ensure non-emptiness. Indeed, given that our first layer is specifically designed to improve upon the open-loop response characteristics of the network described by (1a)-(1c), our assumption is by no means restrictive. Moreover, it highlights one of the key benefits of our approach: the fact that constraint feasibility is certified after first-layer synthesis, with any potential infeasibility being detected early on during the design procedure.

With all of their components now clarified, and assuming that every problem described in (8) is feasible<fnref n="3" /> for each time instant $k \ge k_0$, each MPC-based subcontroller in our second layer can be implemented as per Algorithm 1.

**Algorithm 1 Implementation of the i-th subcontroller**

**Data:** The i-th network's area model and the design parameters of the optimisation problem from (8)

**Result:** The values of $u_{si1}[k]$ and $u_{si2}[k]$

**Step 1a:** Measure the noise-affected initial condition vector $[\tilde{x}_{ci}^\top \tilde{w}_{ci}^\top]^\top$ and transmit it to the subcontrollers of all network areas of index $j \neq i$, for which $i \in \mathcal{N}_j$;

**Step 1b:** Receive the noise-affected initial condition vectors $[\tilde{x}_{cj}^\top \tilde{w}_{cj}^\top]^\top$ from the subcontrollers of all network areas of index $j \in \mathcal{N}_i \setminus \{i\}$;

**Step 2:** Form the vectors $\tilde{\theta}_{xi}[k+t]$ and $\tilde{\theta}_{ui}[k+t]$ for all $t \in \{1: T_i\}$ and use them to construct the sets in (44);

**Step 3:** Solve (8) for two sequences of optimal command signal vectors $u_{s1i}^*[k+t-1]$ and $u_{s2i}^*[k+t-1]$, in which $t \in \{1: T_i + \bar{T}_i\}$;

**Step 4a:** Set $u_{s1i}[k] \leftarrow u_{s1i}^*[k]$ and $u_{s2i}[k] \leftarrow u_{s2i}^*[k]$, then transmit $u_{s1i}[k]$ to the subcontrollers of all network areas of index $j \neq i$, for which $i \in \mathcal{N}_j$;

**Step 4b:** Receive $u_{s1j}[k]$ from the subcontrollers of all network areas of index $j \in \mathcal{N}_i \setminus \{i\}$, then pass them along to the local first-layer subcontroller;

**Step 5:** Wait for the start of the following sampling period and return to Step 1;

**Remark 6.3** It is straightforward to notice that the tuples formed by $(H_{1it}, h_{1it})$ and $(H_{2\ell t}, h_{2\ell t})$ characterise polyhedral inner-approximations of the sets on the right-hand sides of (42)-(43). Crucially, these pairs can always be determined offline, which greatly simplifies the online construction of the sets from (44) in **Step 3** of Algorithm 1. In addition to this computational benefit,

note that no explicit construction of any invariant sets occurs in our setting. For especially large or highly interconnected networks, the type of invariant sets discussed in [8] or [25] may be highly complex, or even numerically intractable for implementation. Our framework allows us to bypass explicit invariance-related concerns, and it allows us to control the complexity of the sets in (44), by means of the inner approximations given in (42)-(43).

Although the implementations proposed in Algorithm 1 show that the MPC-based subcontrollers interact with their counterparts from neighbouring areas, with the aim of exchanging local state measurements, we highlight the fact that all of the individual command signals are computed separately, in a distributed manner.

**Remark 6.4** Regarding the computation of said command signals, the MPC formulations provided in (8) are particularly flexible, in the context of (distributed) constraint satisfaction. The two prediction horizons, introduced in point c) of Section 3.2 and associated with these problems, permit the careful management of the computational burdens which come with solving (8). Since we are chiefly interested in satisfying the latter's constraints only for the sake of recursive feasibility (see the numerical example presented in the sequel), it is possible to improve closed-loop performance of the MPC-based layer by increasing the unconstrained prediction horizon without dramatically impacting the overall computational cost.

The following result shows that the local efforts of these subcontrollers ensure global theoretical guarantees.

**Proposition 6.1** For each of the area dynamics described in (16a)-(16d), let an MPC subcontroller be implemented as in Algorithm 1 and let the communicated information between these subcontrollers be affected as in (20). Assume that:

A1) The first layer's initial condition at time $k = k_0$ satisfies $x_c \in \mathcal{X}$ along with $\text{row}_j(w_{cr\ell}) \in W_{\ell j}$, for all $\ell \in \{1: n_u\}$ and all $j \in \{1: n_{r\ell}\}$;

A2) All of the optimisation problems described in (8) are feasible at time $k = k_0$.

Then, by computing all partitions of $u_{s1}[k_0]$ and $u_{s2}[k_0]$ as in Algorithm 1 and by applying these two vectors as shown in Figure 1, it follows that $x[k_0+1] \in \mathcal{X}$, $u[k_0] \in \mathcal{U}$ and $\text{row}_j(w_{r\ell}[k_0+1]) \in W_{\ell j}$, $\forall \ell \in \{1: n_u\}$, $j \in \{1: n_{r\ell}\}$.

**PROOF.** See Appendix B.

**Remark 6.5** Note that Assumption A1) from Proposition 6.1 hinges upon the initial condition of the controlled network and upon the sets $\mathcal{X}_i$ in (26), since the initial conditions of the implementations in (10) can be imposed at will. Should the system from (1a)-(1c) be initialised in some $x_{nf} \notin \mathcal{X}$, the procedure presented in this paper and in [23] can also be used as a recovery mechanism. The key is to relax the state constraints given by $\mathcal{X}$ into some feasible and connected set $\mathcal{X}_r \supseteq (\mathcal{X} \cup \{x_{nf}\})$, while keeping the same $\mathcal{U}$ and resynthesising the second layer in accordance with the relaxed constraints. Since, in our

<footnote n="3">³ This property will be rigorously addressed in the sequel.</footnote>
framework, the first layer is the one responsible for dis-
turbance rejection and reference tracking, the NRF-based
control laws can then be used to drive the network's state
within X, all while the second layer handles command
limitation. The first layer's subcontrollers can then be
reinitialised to satisfy Assumption A1), and nominal op-
eration may commence using the original second layer.

When starting from a point inside the state-constraining
set X, Proposition 6.1 shows that our control scheme is
able to maintain the first layer's state variables inside X
while also satisfying the command constraints encoded
by U, provided that all of the MPC-based problems from
(8) are feasible. Given the critical importance of satisfy-
ing Assumption A2) from Proposition 6.1, we now pro-
ceed to investigate a set of conditions which offer strong
theoretical guarantees regarding the feasibility of the op-
timisation problems from (8).

6.2 Guaranteeing Recursive Feasibility

In order to extend the theoretical guarantees obtained in
Proposition 6.1 beyond the initial time instant, we must
ensure that each of the MPC-based control laws from (8)
is feasible for every time instant $k \in \mathbb{N}$ with $k \ge k_0$. This
property, known as *recursive feasibility* (see Section 12.3
in [5]), plays a key role in our chosen approach and, con-
sequently, the following theorem (which directly builds
upon the results obtained in Proposition 6.1) represents
the paper's main theoretical result.

**Theorem 6.1** For each of the area dynamics described in (16a)-(16d), let an MPC subcontroller be implemented as in Algorithm 1 and let the communicated information between these subcontrollers be affected as in (20). For all $i \in \{1:N\}$, let $\rho_i \in \mathbb{N}$ with $\rho_i \ge 1$ and define

$$
\tilde{\Theta}_{it} := \left\{ \begin{array}{l} \left[ \theta_x^\top \theta_{u1}^\top \dots \theta_{un_{ui}}^\top \right]^\top \quad | \quad \theta_x \in \tilde{\Theta}_{xit}, \\ \theta_{uj} \in (e_j^\top \tilde{\Theta}_{uit}), \quad j \in \{1:n_{ui}\} \end{array} \right\}, \quad (45a)
$$

$$
\tilde{P}_{it} := \left\{ \begin{array}{l} \left[ p_x^\top \ p_{u1}^\top \ \dots \ p_{un_{ui}}^\top \right]^\top \quad | \quad p_x \in \mathcal{P}_{1it}, \\ p_{uj} \in \mathcal{P}_{2(j+\alpha_{ui})t}, \quad j \in \{1:n_{ui}\} \end{array} \right\}. \quad (45b)
$$

Assume that:

A1) The first layer's initial condition at time $k = k_0$ satisfies $x_c \in \mathcal{X}$ along with $\text{row}_j(w_{cr\ell}) \in W_{\ell j}$, for all $\ell \in \{1:n_u\}$ and all $j \in \{1:n_{r\ell}\}$;

A2) For each $i \in \{1:N\}$, there exists $\rho_i \in \mathbb{N}$ with $\rho_i \ge 1$ which satisfies

$$
\tilde{\Theta}_{it} \subseteq \left( \left[ \begin{array}{c} -C_{xi} \\ -C_{ui} \end{array} \right] \bigoplus_{j=1}^{t} A_{si}^{j-1} \left( (B_{s1i} \mathcal{U}_{s1i}) \oplus (B_{s2i} \mathcal{U}_{s2i}) \right) \right) \oplus \tilde{P}_{it}, (46)
$$
for all $t \in \{1 : \rho_i\}$.

Then, by setting $T_i = \rho_i, \forall i \in \{1:N\}$, we have that:

i) Each of the optimisation problems described in (8) becomes recursively feasible;

ii) By successively computing all partitions of $u_{s1}[k]$ and $u_{s2}[k]$ as in Algorithm 1 and by applying these two vectors as shown in Figure 1, it follows that $x[k+1] \in \mathcal{X}$, $u[k] \in \mathcal{U}$ and $\text{row}_j(w_{r\ell}[k+1]) \in W_{\ell j}$, for all $\ell \in \{1:n_u\}$, $j \in \{1:n_{r\ell}\}$ and $k \ge k_0$.

PROOF. See Appendix B.

**Remark 6.6** Notice the fact that it is possible to retrieve the above-mentioned guarantees by taking $\rho_i = 1$, for every $i \in \{1:N\}$, resulting in especially efficient implementations (see the numerical example from the sequel) from a computational standpoint. Indeed, as the values of $\rho_i$ increase, not only does (46) become more difficult to satisfy, but the additional constraints associated with taking $T_i = \rho_i$ may dramatically increase the cost of solving (8), depending on the complexity of the sets in (44). As discussed in Remark 6.4, taking larger values for $\bar{T}_i$ may yield similar degrees of performance as increasing $T_i$, without greatly impacting overall computational costs.

To summarise the design workflow of our architecture's
second layer, we provide a high-level algorithm (located
at the top of the next page), which condenses the various
topics discussed in this paper into a sequence of steps.

**Remark 6.7** Performing Steps 3a-3c in Algorithm 2 makes the paramount importance of our architecture's first layer readily apparent. It is precisely due to the proper tuning of the NRF-based control laws that area cross-coupling may be diminished (and potentially rendered negligible, for all practical purposes) to such a degree that the MPC-based subcontrollers may be tuned individually via (30a)-(30b). While this effectively bypasses the need to carefully construct the sets from (30a)-(30b) and it crucially prevents the design procedure of the MPC-based subcontrollers from being co-dependent in any significant way, we point out that said procedure does indeed require information pertaining to the first-layer response of other areas. However, this information is fully available once the NRF-based control laws have been obtained, and it is precisely the role of these distributed controllers to ensure the feasibility (recall Remark 6.2) and the relaxed co-dependency of the second layer's design.

6.3 Computational Considerations

The following discussion is centred around the three core
concepts which underpin the computational burden of
our procedure, namely: synthesis costs, design scalabil-
ity, and implementation complexity.

6.3.1 Synthesis Costs

One of the main features associated with Algorithm 2 is
the fact that all of its constituent steps (along with the
design phase of the first layer, see our companion paper
[23]) can be performed offline, and that any reconfigu-
ration of the second layer's constraints or cost functions
does not impact the first layer's control laws (in contrast
to the technique from [1,2], recalling Remark 6.1). The
chief computational burden associated with our solution
is represented by the cost of computing the sets from
**Algorithm 2** MPC layer design procedure

**Data:** The first layer's closed-loop dynamics from (14), the objective sets $\mathcal{X}$ and $\mathcal{U}$, along with the exogenous signal sets $\mathcal{V}$, $\mathcal{D}_s$ and the ones from (27)

**Result:** MPC policies obtained by solving (8)

**Initialisation:** The state-space dynamics (16a)-(16d);

**Step 1:** Split the sets $\mathcal{X}$ and $\mathcal{U}$ into their area-wise components, as shown in (26);

**Step 2a:** Form the noise propagation sets given in (35) and compute the sets from (36) and (39a);

**Step 2b:** Form the disturbance propagation sets in (28);

**Step 3a:** Choose the sets from (30a)-(30b);

**Step 3b:** Form the sets in (31)-(33), (34a)-(34b), those from (37)-(38) and the ones from (39b)-(41);

**Step 3c:** Compute the polyhedral inner-approximations given in (42)-(43) and the sets from (45a)-(45b);

**Step 3d:** Set $\rho_i = 1$ and check (46) for all $i \in \{1:N\}$;

if (46) holds for all $i \in \{1:N\}$ then

go to **Step 4a**;

else

go to **Step 3a** and retune the sets from (30a)-(30b) in accordance with Remark 5.6;

end

**Step 4a:** Try increasingly larger values for $\rho_i \ge 2$ and then check (46) for all $i \in \{1:N\}$;

**Step 4b:** For all $i \in \{1:N\}$, set $T_i$ as the largest found value of $\rho_i$ which satisfies (46);

**Step 4c:** For all $i \in \{1:N\}$, select an appropriate $\bar{T}_i$ and the weighting matrices of the cost function in (8);

(34a)-(34b) along with (42)-(43) for large-scale applications. Although the measures described in Remarks 5.5 and 6.3 can alleviate some of the complexity associated with these operations, the latter remain heavily dependent on the shape and the structure of the involved sets, making it challenging to come up with a one-size-fits-all computational cost. In spite of this, switching to a homothetic [18] or even zonotope-based [13] framework (while employing various approximations, either in terms of the operation or of the operands themselves) is guaranteed to have a significant, positive impact on the complexity of the synthesis procedure. For a practical and specific example of the costs associated with the design of the second layer, see the numerical example presented in the sequel, and the source files from its GitHub repository.

### 6.3.2 Design Scalability

While all of the steps which form Algorithm 2 can, indeed, be iterated in a fully offline setting, yet another appealing feature is the fact that once the second layer's command budgeting has been negotiated via the selection of the sets from (30a)-(30b), the subsequent set-based computations from Algorithm 2 can be performed in parallel. As indicated in Remark 6.7, the design procedure of each area can proceed independently of the

others, once the aforementioned choice has been made. This fact is a direct consequence of not employing any notions of joint set invariance (see, for example, [25]) and of relying, instead, upon the decoupling-focused control laws discussed in our companion paper [23].

### 6.3.3 Implementation Complexity

Due to the separate implementations in (10) and in Algorithm 1, a crucial observation is the fact that the runtime cost of our architecture is divided among the network's areas. Said burden is appealingly modest, since it reduces to several inexpensive operations, namely: four matrix-vector multiplications and two vector-vector additions for an area's first-layer subcontroller (10), along with the solving of (8) for the second layer. We point out that the latter problem boils down to a standard Quadratic Programming (QP) optimisation problem, whenever the cost in (8) has a quadratic (and symmetrical) expression and the sets from (30a)-(30b) are chosen to be polyhedral in nature (which is always possible, as per Remark 5.3). Despite the fact that the overall complexity of this problem is highly dependent upon the number of constraints it includes (recall Remark 6.3), note that numerous solvers and optimisation toolboxes (such as [3] and [14], respectively) are widely available for the efficient solving of such problems.

## 7 Numerical Example

### 7.1 Network Model and Control Objectives

#### 7.1.1 The Dynamics of the Application

In order to showcase the effectiveness of the proposed design framework, we now focus on a numerical example centred around the vehicle platoon application discussed in [20]. We consider a formation of $N = 10$ identical vehicles, in which the continuous-time models from (22) in [20] are configured with the following parameters: $m_i = 1$, $\tau_i = 0.1$ and $\sigma_i = 1$, $\forall i \in \{1:N\}$. The resulting zero-order hold discretised dynamics of each car, for a sampling period of $T_s = 0.1$, are given by

$$ \begin{bmatrix} p_i[k+1] \\ v_i[k+1] \\ \mu_i[k+1] \end{bmatrix} = \underbrace{\begin{bmatrix} 1 & 0.1 & -0.0331 \\ 0 & 1 & -0.5689 \\ 0 & 0 & 0.3679 \end{bmatrix}}_{A_{car}} \begin{bmatrix} p_i[k] \\ v_i[k] \\ \mu_i[k] \end{bmatrix} + \underbrace{\begin{bmatrix} 0.0381 \\ 0.6689 \\ 0.6321 \end{bmatrix}}_{B_{car}} u_i[k], $$

where $u_i$ is the control signal of the $i^{th}$ car, $p_i[k]$ denotes its (absolute) position, $v_i[k]$ is its speed and $\mu_i[k]$ is the state of the $i^{th}$ car's actuator (the $\Phi$ terms in [20]).

In addition to these variables, we also consider the "position" of the $0^{th}$ car (which acts as a reference signal to the lead car in the platoon) denoted $p_0[k]$, the latter's increment in one time instant, denoted $\Delta p_0[k]$, and the length of the vehicle located in front of the platoon's $i^{th}$ car, denoted $\ell_i[k]$. Forming $\ell_i^T[k] := [\ell_1[k] \dots \ell_N[k]]$, the dynamics of these variables are captured by

$$ \begin{bmatrix} p_0[k+1] \\ \ell_v[k+1] \end{bmatrix} = I_{N+1} \begin{bmatrix} p_0[k] \\ \ell_v[k] \end{bmatrix} + e_1 \Delta p_0[k]. $$
For the homogeneous platoon considered in this application, we consider $\ell_i[k] = 5$ for all $i \in \{2:N\}$ along with $\ell_1[k] = 0$, since the 0^th$ (reference) car is a virtual one.

### 7.1.2 First-Layer Objectives and Platoon Model

In keeping with the control laws proposed in [20], we wish to designate each car as an individual area and to obtain a first layer in which the control signals are computed (in the absence of exogenous or second-layer signals) as

$$ \text{row}_i(u_f[k]) = \begin{cases} \mathbf{\Gamma}_1(z) \star (y_1[k] + h v_1[k]), & i = 1, \\ \mathbf{\Gamma}_i(z) \star (y_i[k] + h v_i[k]) + \mathbf{\Phi}_i(z) \star \text{row}_{i-1}(u_f[k]), & i \neq 1 \end{cases} \quad (47) $$

for all $i \in \{2:N\}$, where $y_i[k] := p_i[k] + \ell_i[k] - p_{i-1}[k]$ stands for the $i^{th}$ car's length-based interdistance and $h > 0$ is the platoon's time headway (see [20] and the references therein for its theoretical justification), which has been chosen as $h = 5$ for the considered application.

Given the control laws of type (47) we wish to obtain, and the fact that we will want to impose constraints on the $y_i$ signals (to prevent traffic collisions), we will model our platoon in such a way as to have the headway-based interdistances as state variables. Thus, we form

$$ T = \begin{bmatrix} I_{N+1} & O & O & \cdots & O & O \\ T_{\ell 1} & I_3 & O & O & O & O \\ T_{\ell 2} & T_{\ell d} & I_3 & O & \cdots & O & O \\ T_{\ell 3} & O & T_{\ell d} & I_3 & \cdots & O & O \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ T_{\ell N} & O & O & O & \cdots & T_{\ell d} & I_3 \end{bmatrix} \in \mathbb{R}^{(4N+1) \times (4N+1)}, $$

where $T_{\ell d} = [-e_1 \ O \ O] \in \mathbb{R}^{3 \times 3}$ and

$$ T_{\ell i} = \begin{cases} [e_2 - e_1 \ O \ O]^\top \in \mathbb{R}^{3 \times (N+1)}, & i = 1, \\ [e_{i+1} \ O \ O]^\top \in \mathbb{R}^{3 \times (N+1)}, & \forall i \in \{2:N\}. \end{cases} $$

to obtain the platoon's model of type (1a)-(1c), given by

$$ x[k+1] = Ax[k] + Bu[k] + e_1 \Delta p_0[k], \quad (48) $$

where we have that:

$$ \text{a) } A := T \cdot \operatorname{diag}(I_{N+1}, \underbrace{A_{car}, \dots, A_{car}}_{N \text{ times}}) \cdot T^{-1}; $$

$$ \text{b) } B := T \cdot [O \ I_{3N}]^\top \cdot \operatorname{diag}(\underbrace{B_{car}, \dots, B_{car}}_{N \text{ times}}); $$

c) the platoon's state vector is expressed as follows

$$ x[k] := [p_0[k] \ \ell_v^\top[k] \ x_1^\top[k] \ \dots \ x_N^\top[k]]^\top, \text{ where} $$

$$ x_i[k] := [y_i[k] \ v_i[k] \ \mu_i[k]]^\top \text{ denotes the vector containing the platoon states of the } i^{\text{th}} \text{ area;} $$

d) the platoon's input vector is given explicitly by

$$ u[k] := [u_1^\top[k] \ \dots \ u_N^\top[k]]^\top; $$

$$ \text{e) } B_d := e_1 \in \mathbb{R}^{4N+1} \text{ and } d[k] := \Delta p_0[k], \text{ by direct correspondence with (1a)-(1c).} $$

By choosing $C_g = \{z \in \mathbb{C} \mid |z| < 1\}$, the first layer will aim to stabilise platoon's closed-loop dynamics, in the sense that, by implementing the control laws of type (47)

as detailed in our companion paper [23], the closed-loop dynamics will be described by

$$ \begin{bmatrix} p_0[k+1] \\ \ell_v[k+1] \\ x_f[k+1] \end{bmatrix} = \begin{bmatrix} 1 & O & O \\ O & I_N & O \\ O & O & A_f \end{bmatrix} \begin{bmatrix} p_0[k] \\ \ell_v[k] \\ x_f[k] \end{bmatrix} + \begin{bmatrix} B_0 \\ B_l \\ B_f \end{bmatrix} d_s[k], \quad (49) $$

in which:

a) the eigenvalues of $A_f$ are all in $C_g$;

b) $x_f[k] = [x_1^\top[k] \ w_1^\top[k] \ \dots \ x_N^\top[k] \ w_N^\top[k]]^\top$, with $w_i[k]$ being the state vectors belonging to the NRF implementations of type (10), for all $i \in \{1:N\}$;

c) $d_s[k]$ is the signal vector from (7).

**Remark 7.1** Since $p_0[k]$ and $\ell_i[k]$, for any $i \in \{1:N\}$, cannot be influenced by the $u_i[k]$ signals, the first layer will be unable to stabilise their dynamics. However, given the nature of our application and the physical interpretation of the variables, the aforementioned stabilisation is of no practical interest. The second layer will focus only on the identity $x_f[k+1] = A_f x_f[k] + B_f d_s[k]$, which will represent our first layer's closed-loop dynamics.

### 7.1.3 Second-Layer Objectives

The main goal of the second layer is to ensure constraint satisfaction, with respect to the variables that make up $x_i[k]$ and $u_i[k]$, in the presence of recursive feasibility. Thus, for all $i \in \{1:N\}$ and all $k \ge k_0$, we wish to have:

$$ y_i[k] \in [-360, 0], \quad v_i[k] \in [0, 36], \quad u_i[k] \in [-10, 10]. \quad (50) $$

These constraints must be satisfied while employing the same communication scheme as the one considered in (47), i.e., the $i^{th}$ car can employ only information that is accessed locally or from the car directly in front of it.

Additionally, due to the fact that the prototype NRF controller from [20] (if properly tuned) ensures a high level of performance, we want the second layer to intervene as little as possible in the closed-loop dynamics of the first one. In other words, most of the control action should be undertaken by the $u_{fi}[k]$ signals, whereas the entries of $u_{s1i}[k]$ and $u_{s2i}[k]$ should have negligible absolute values (compared to those of $u_{fi}[k]$) whenever the constraints from (50) are not in danger of being violated.

## 7.2 Control Architecture Design

### 7.2.1 The First Layer

Since the application considered in this numerical example already has a dedicated NRF-based solution, we defer to [20] for the procedural details regarding the design of our first layer. In so doing, we briefly discuss the numerical implementation of these control laws, before focusing our attention on the design of the second layer.

For the homogeneous platoon whose physical parameters were given in the previous subsection, the first-layer command signals will be computed as stated in (12), via the state-space implementations of type (10) given by

$$ u_{f1}[k] = \begin{bmatrix} \frac{a_1 - z}{1} & B \Gamma_1 \\ \frac{1}{O} & O \end{bmatrix} \star (x_1[k] + \zeta_1[k] + u_{s11}[k] + \beta_{s11}[k]), \quad (51) $$
$$ \begin{bmatrix} x_1[k+1] \\ w_1[k+1] \end{bmatrix} = \begin{bmatrix} A_{car} & B_{car} \\ B_{\Gamma 1} & a_1 \end{bmatrix} \begin{bmatrix} x_1[k] \\ w_1[k] \end{bmatrix} + \begin{bmatrix} B_{car} & O \\ O & B_{\Gamma 1} \end{bmatrix} \begin{bmatrix} u_{s21}[k] \\ u_{s11}[k] \end{bmatrix} + \begin{bmatrix} -e_1 & B_{car} & O \\ 0 & O & B_{\Gamma 1} \end{bmatrix} \begin{bmatrix} \Delta p_0[k] \\ \beta_{s21}[k] \\ \zeta_{1}[k] + \beta_{s11}[k] \end{bmatrix}, \quad (53a) $$

$$ \begin{bmatrix} x_i[k+1] \\ w_i[k+1] \end{bmatrix} = \begin{bmatrix} A_{car} & B_{car} \\ B_{\Gamma i} & a_i \end{bmatrix} \begin{bmatrix} x_i[k] \\ w_i[k] \end{bmatrix} + \begin{bmatrix} B_{car} & O \\ O & B_{\Gamma i} \end{bmatrix} \begin{bmatrix} u_{s2i}[k] \\ u_{s1i}[k] \end{bmatrix} + \begin{bmatrix} A_w & B_w \\ O & b_{\Phi i} \end{bmatrix} \begin{bmatrix} x_{i-1}[k] \\ w_{i-1}[k] \end{bmatrix} + \begin{bmatrix} B_w & B_{car} & O \\ O & O & b_{\Phi i} & B_{\Gamma i} \end{bmatrix} \begin{bmatrix} u_{s2(i-1)[k]} + \beta_{s2(i-1)[k]} \\ \beta_{s2i}[k] \\ \beta_{f(i-1)[k]} \\ \zeta_i[k] + \beta_{s1i}[k] \end{bmatrix}, \quad \forall i \in \{2:N\}. \quad (53b) $$

$$ \hat{\chi} := \{ \hat{x} \in \mathbb{R}^4 \mid -360 \le e_1^\top \hat{x} \le 0, 0 \le e_2^\top \hat{x} \le 36, |e_3^\top \hat{x}| \le 10, |e_4^\top \hat{x}| \le 4.99, 5.01 |e_1^\top B_{car}| \le W_x \hat{x} \le 36T_s - 5.01 |e_1^\top B_{car}| \}. \quad (54) $$

and, for all $i \in \{2:N\}$, by

$$ u_{fi}[k] = \begin{bmatrix} a_i - z & b_{\Phi i} & B_{\Gamma i} \\ 1 & 0 & O \end{bmatrix} \star \begin{bmatrix} u_{f(i-1)[k]} + \beta_{f(i-1)[k]} \\ x_i[k] + \zeta_i[k] + u_{s1i}[k] + \beta_{s1i}[k] \end{bmatrix}, \quad (52) $$

with $a_i$, $b_{\Phi i}$ and $b_{\Gamma i}$ being given in Table 1 below.

<table><caption>Table 1
Coefficients of the first-layer controller implementations
By using the state-space systems in (48) and (51)-(52) to form the first layer's closed-loop system, we obtain a realisation of type (49) in which the spectral radius of $A_f$ is 0.9936, thus validating all of the design objectives associated with the first layer. In addition to this, the resulting closed-loop dynamics of each car in the platoon will be given by the realisations from (53a)-(53b), which are located at the top of this page and in which</caption><thead><tr><th>i</th><th>a<sub>i</sub></th><th>b<sub>Φi</sub></th><th>B<sub>Γi</sub></th></tr></thead><tbody><tr><td>1</td><td>0.9690</td><td>0</td><td>-0.0038 -0.0192 0</td></tr><tr><td>2</td><td>0.9799</td><td>0.0199</td><td>-0.0030 -0.0152 0</td></tr><tr><td>3</td><td>0.9799</td><td>0.0200</td><td>-0.0032 -0.0161 0</td></tr><tr><td>4</td><td>0.9798</td><td>0.0200</td><td>-0.0034 -0.0171 0</td></tr><tr><td>5</td><td>0.9797</td><td>0.0200</td><td>-0.0036 -0.0182 0</td></tr><tr><td>6</td><td>0.9796</td><td>0.0201</td><td>-0.0039 -0.0195 0</td></tr><tr><td>7</td><td>0.9795</td><td>0.0201</td><td>-0.0042 -0.0209 0</td></tr><tr><td>8</td><td>0.9794</td><td>0.0202</td><td>-0.0045 -0.0224 0</td></tr><tr><td>9</td><td>0.9793</td><td>0.0202</td><td>-0.0049 -0.0243 0</td></tr><tr><td>10</td><td>0.9792</td><td>0.0203</td><td>-0.0053 -0.0265 0</td></tr></tbody></table>

$$ A_w = \begin{bmatrix} 0 & -0.1 & 0.0331 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad \text{and} \quad B_w = \begin{bmatrix} -0.0381 \\ 0 \\ 0 \end{bmatrix}. $$

It is precisely the dynamics from (53a)-(53b) which will be tackled by our distributed MPC-based policies.

### 7.2.2 The Second Layer

The first step in designing our control architecture's second layer, with the aim of ensuring its recursive feasibility, consists of defining all of the sets discussed in Section 5. We begin by addressing the command signals

generated by the MPC layer, for which we impose:

$$ \left\{ \begin{array}{l} U_{s1i} = \{ u_{s1i} \in \mathbb{R}^3 \mid |e_1^\top u_{s1i}| \le 720, |e_2^\top u_{s1i}| \le 72, e_3^\top u_{s1i} = 0 \}, \forall i \in \{1:N\}, \\ U_{s2i} = \{ u_{s2i} \in \mathbb{R} \mid |u_{s2i}| \le 5 \}, \forall i \in \{1:N\}. \end{array} \right. $$

Next, we consider the following box-type constraints for the exogenous signals:

$$ \left\{ \begin{array}{l} V = \{ \nu \in \mathbb{R}^{n+n_w} \mid |e_j^\top \nu| \le 0.02, \forall j \in \{1:n+n_w\} \}, \\ D_{\beta_f} = \{ \beta_f \in \mathbb{R}^m \mid |e_j^\top \beta_f| \le 0.02, \forall j \in \{1:m\} \}, \\ D_{\beta_{s1}} = \{ \beta_{s1} \in \mathbb{R}^n \mid |e_j^\top \beta_{s1}| \le 0.01, \forall j \in \{1:n\} \}, \\ D_{\beta_{s2}} = \{ \beta_{s2} \in \mathbb{R}^m \mid |e_j^\top \beta_{s2}| \le 0.01, \forall j \in \{1:m\} \}, \\ D_d = \{ \Delta p_0 \in \mathbb{R} \mid 0 \le \Delta p_0 \le 36T_s \}, \end{array} \right. $$

and, finally, we address the constraints related to the platoon's states by first introducing the matrix denoted $W_x := \text{row}_1([A_{car} \ B_{car}]) \text{diag}(0, I_3)$ and then constructing the set $\hat{\chi} \subseteq \mathbb{R}^4$ given explicitly in (54), which is located at the top of this page.

By choosing $T_i = 1$ for all $i \in \{1:N\}$ and by leveraging $\hat{\chi}$, while also accounting for the disturbance discussed in Section 4.3, we proceed to implement our design procedure. The resulting numerical routines and all of the data sets from this numerical example are available at https://github.com/AndreiSperila/NRMPC. Moreover, Table 2 (located at the top of the next page) shows the computation times associated with the complete design of every car's second-layer subcontroller.

**Remark 7.2** As explained in the comments that accompany the source code from the above-indicated link, our approach entails the implicit design of the robust control invariant sets $C_i \subseteq \hat{\chi}$, which (by their very construction; see, for example, Section 10.4 in [5]) guarantee the recursive feasibility of our second layer. We stress the fact that (as previously indicated in Section 6) the sets $C_i$ are not computed explicitly, with their existence being a desirable and natural consequence of our proposed approach. Finally, we choose $\bar{T}_i = 0$ and we set the costs in (8) as

$$ 10^{-9} y_i^2[k+1] + u_{s2i}^2[k] + u_{s1i}^\top[k] u_{s1i}[k], \quad (55) $$

thereby concluding the design phase of the second layer.
<table><caption>Table 2
Computation times associated with second-layer synthesis</caption><thead><tr><th>car number</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th></tr></thead><tbody><tr><td>time (seconds)</td><td>211.740</td><td>757.767</td><td>390.860</td><td>396.871</td><td>295.800</td><td>235.603</td><td>155.952</td><td>204.719</td><td>74.838</td><td>71.048</td></tr></tbody></table>

### 7.3 Simulation Scenario and Results

In order to test the obtained control laws, we simulate a scenario in which $\Delta p_0[k] = T_s v_0[k]$, where

$$v_0[k] = 10 \cdot 1[k] - 7 \cdot 1[k-400] + 30 \cdot 1[k-1200] - 30 \cdot 1[k-1300]$$

is the speed of the virtual 0th car and 1[k] denotes the discrete-time Heaviside step function. As can be seen in Figure 3 below, in which $v_0[k]$ is depicted in orange, the speed of each car in the platoon will try to match that of the 0th car in order to achieve steady-state evolution.

<figure><img src="image_3.png" /><figcaption>Figure 3. 0th car speed and speeds of the platoon's cars</figcaption></figure>

Moreover, as can be seen in Figure 4, our two-layer control architecture successfully maintains all of the inter-distance variables within their prescribed bounds, thus avoiding collisions and preventing the cars from drifting too far apart (leading to traffic congestion or the break-down of the communication between consecutive cars).

<figure><img src="image_4.png" /><figcaption>Figure 4. Length-based interdistances of the platoon's cars</figcaption></figure>

With respect to the limitations placed upon the signals sent to the car actuators, our control scheme is once again successful. As can be seen from Figures 5 and 6 below, the sum of the signals from the two aforementioned figures is well within the bounds from (50), even taking into account the amplitudes of the entries in $\beta_{s2}[k]$.

<figure><img src="image_5.png" /><figcaption>Figure 5. Outputs of the first-layer subcontrollers</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 6. Second-layer outputs to the platoon</figcaption></figure>

Finally, due to the choice of cost functions in (55), the entries of $u_{s1}[k]$ and $u_{s2}[k]$ have absolute values of at most $10^{-6}$ for most of the simulation scenario's duration, as can be seen in Figures 6 and 7.

<figure><img src="image_7.png" /><figcaption>Figure 7. Norm values for second-layer outputs to the first-layer subcontrollers</figcaption></figure>

The brief intervals in which the two signal vectors are of the same order of magnitude as the first-layer control signals are those in which the second layer intervenes in order to saturate either the output of the first car's NRF-based subcontroller (see Figure 8 below), or the speed of the lead car in the platoon (see Figure 9, also located below) with respect to these variables' prescribed bounds.
<table><caption>Table 3
NRF-based subcontroller runtime statistics in milliseconds</caption><thead><tr><th>car number</th><th>min</th><th>max</th><th>mean</th><th>median</th><th>mode</th></tr></thead><tbody><tr><td>1</td><td>0.003</td><td>0.213</td><td>0.005</td><td>0.004</td><td>0.004</td></tr><tr><td>2</td><td>0.001</td><td>0.057</td><td>0.003</td><td>0.002</td><td>0.002</td></tr><tr><td>3</td><td>0.001</td><td>0.023</td><td>0.002</td><td>0.002</td><td>0.002</td></tr><tr><td>4</td><td>0.001</td><td>0.022</td><td>0.002</td><td>0.002</td><td>0.002</td></tr><tr><td>5</td><td>0.001</td><td>0.022</td><td>0.002</td><td>0.002</td><td>0.002</td></tr><tr><td>6</td><td>0.001</td><td>0.233</td><td>0.002</td><td>0.002</td><td>0.002</td></tr><tr><td>7</td><td>0.001</td><td>0.008</td><td>0.002</td><td>0.002</td><td>0.002</td></tr><tr><td>8</td><td>0.001</td><td>0.010</td><td>0.002</td><td>0.002</td><td>0.002</td></tr><tr><td>9</td><td>0.001</td><td>0.013</td><td>0.002</td><td>0.002</td><td>0.002</td></tr><tr><td>10</td><td>0.001</td><td>0.011</td><td>0.002</td><td>0.002</td><td>0.002</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 8. Constraint satisfaction for the first-layer command of the platoon's first car (second-layer subcontroller active in the green shaded area)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 9. Constraint satisfaction for the speed of the platoon's first car (second-layer subcontroller active in the green shaded area)</figcaption></figure>

Finally, as anticipated in Section 6.3.3 and as showcased by the runtime statistics showcased in Tables 3 and 4 above, the implementation costs associated with our proposed solution are sufficiently modest such that command computation is well below the sampling period of our control system, namely $T_s = 0.1$ seconds. Note that all simulations associated with these runtimes were performed via MATLAB R2024b, on a mid-range laptop equipped with an i9-100885H Intel Core processor.

## 7.4 Comparison with Alternative Solutions

In order to conclude this numerical example, we present a comparison with several prominent alternative tech-

<table><caption>Table 4
MPC-based subcontroller runtime statistics in milliseconds</caption><thead><tr><th>car number</th><th>min</th><th>max</th><th>mean</th><th>median</th><th>mode</th></tr></thead><tbody><tr><td>1</td><td>1.678</td><td>4.375</td><td>2.202</td><td>2.127</td><td>1.872</td></tr><tr><td>2</td><td>1.688</td><td>5.028</td><td>2.028</td><td>1.955</td><td>1.877</td></tr><tr><td>3</td><td>1.659</td><td>4.347</td><td>1.949</td><td>1.884</td><td>1.881</td></tr><tr><td>4</td><td>1.661</td><td>3.882</td><td>1.930</td><td>1.868</td><td>1.742</td></tr><tr><td>5</td><td>1.648</td><td>3.929</td><td>1.921</td><td>1.861</td><td>1.853</td></tr><tr><td>6</td><td>1.630</td><td>4.936</td><td>1.931</td><td>1.864</td><td>1.826</td></tr><tr><td>7</td><td>1.629</td><td>4.895</td><td>1.913</td><td>1.848</td><td>1.807</td></tr><tr><td>8</td><td>1.624</td><td>4.829</td><td>1.924</td><td>1.852</td><td>1.852</td></tr><tr><td>9</td><td>1.639</td><td>4.233</td><td>1.956</td><td>1.838</td><td>1.847</td></tr><tr><td>10</td><td>1.614</td><td>4.680</td><td>1.989</td><td>1.838</td><td>1.786</td></tr></tbody></table>

niques from literature, which tackle the platooning problem in an MPC-based framework. We refer here to the procedures proposed in [31,28,6,10,27] and we point out that a major limitation for all of them is the fact that, unlike our proposed solution, no source of disturbance (either from exogenous sources or from model uncertainty) is considered in any of these approaches. Moreover, the control laws presented in [31,28,6] all share in common the fact that they do not explicitly support constraints for the state variables of the platoon's cars, except for a restrictive use of terminal constraints, which are meant to ensure a number of theoretical guarantees.

Another common feature of the cited works is the fact that the lead cars in [28,6,10] must broadcast information to all other vehicles in the formation, which places an impractical amount of strain on the communication infrastructure of the platoon. In addition to this fact, the property of asymptotic closed-loop stability is retrieved in [31,6,10] by carefully tuning the stage costs of the MPC-based subcontrollers, rather than it being an inbuilt feature, as ensured by the NRF-based first layer (see [23] and [20]). This being said, we now address the specific aspects of the five aforementioned techniques:

a) One of the key features of the procedure from [31] is the fact that the above-mentioned terminal constraint requires the computation of a *regional average value*, evaluated with respect to a number of neighbouring cars. In addition to the demands this places on the communication infrastructure of the platoon, the feasibility of the resulting control problem is noticeably dependent on choosing a large enough prediction horizon. While the overall computational complexity of the optimisation problem to be solved is similar to that of a QP-based formulation in (8), increasing the prediction horizon may have a significant computational impact which can be bypassed in our framework (recall Remark 6.6).

b) The main limitation of the technique in [28] is the fact that the optimisation problem associated with the MPC-based subcontrollers has to be run twice during a car's command computation routine, with said problem having a comparable cost to a QP-based formulation of our problem in (8). Furthermore, we point out that both the asymptotic and
string stability of the closed-loop system are only guaranteed when the platoon's initial condition satisfies a series of restrictive assumptions.

c) In addition to the list of common restrictions detailed above, the solution proposed in [6] is affected by another limitation; namely, that it can only be initialised when the platoon is already in equilibrium. In spite of these disadvantages, we point out that the procedure's computational cost is appealingly modest, being roughly analogous to that of our own procedure's QP-based formulation, since it does not require the double iteration from [28].

d) The approach from [10] possesses some structural similarity with the one presented in this paper, in the sense that it includes a rudimentary, PID-based lower control layer, while also enabling the use of limited, box-based constraints for the platoon's states. However, we point out that the lower-level PIDs are not only unable to ensure closed-loop string stability (in contrast to our first layer), but the NRF-based subcontrollers offer many more degrees of freedom, in terms of both sparsity patterns and closed-loop performance (see [23]). Since [10] is heavily inspired by the approach presented in [6], the same computational benefits apply to it as well.

e) Finally, the procedure laid out in [27] proposes a centralised take on the problem and, unsurprisingly, it showcases a high degree of closed-loop performance in the presence of state constraints. In addition to being intractable in the distributed setting, its critical aspect is given by the fact that closed-loop string stability is incompatible with the aforementioned constraints. In contrast to our solution, when safety-critical constraints become active, the technique given in [27] sacrifices the platoon's string stability for guaranteed constraint satisfaction.

## 8 Conclusions

As a contribution to the synthesis of distributed constraint management techniques, we have shown that it is possible to elegantly decouple structural design constraints (which represent the focal point of the distributed control paradigm) from the time-based constraints which are ubiquitous in MPC-oriented literature. By leveraging a promising class of distributed control strategies [21,23,24] and by adapting a series of set-theoretical notions [4] to the distributed setting, we obtain a computationally efficient control architecture which provides *powerful theoretical guarantees*, and whose components can be *designed fully offline*.

## Acknowledgements

The authors would like to thank both Şerban Sabău and Cristian Oară for the many fruitful discussions which ultimately led to the elaboration of this manuscript.

## References

[1] Carmen Amo Alonso, Jing Shuang Li, James Anderson, and Nikolai Matni. Distributed and Localized Model-Predictive

Control—Part I: Synthesis and Implementation. *IEEE Transactions on Control of Network Systems*, 10(2):1058–1068, 2023.

[2] Carmen Amo Alonso, Jing Shuang Li, Nikolai Matni, and James Anderson. Distributed and Localized Model Predictive Control—Part II: Theoretical Guarantees. *IEEE Transactions on Control of Network Systems*, 10(3):1113–1123, 2023.

[3] MOSEK ApS. MOSEK Optimization Toolbox for MATLAB 11.0.25, 2025.

[4] Franco Blanchini and Stefano Miani. Set-Theoretic Methods in Control, Second Edition. Birkhäuser, 2015.

[5] Francesco Borrelli, Alberto Bemporad, and Manfred Morari. Predictive Control for Linear and Hybrid Systems. Cambridge University Press, 2017.

[6] William B. Dunbar and Derek S. Caveney. Distributed Receding Horizon Control of Vehicle Platoons: Stability and String Stability. *IEEE Transactions on Automatic Control*, 57(3):620–633, 2012.

[7] Emanuele Garone, Stefano Di Cairano and Ilya Kolmanovsky. Reference and command governors for systems with constraints: A survey on theory and applications. *Automatica*, 75:306–328, 2017.

[8] Marcello Farina and Riccardo Scattolini. Distributed predictive control: A non-cooperative algorithm with neighbor-to-neighbor communication for linear systems. *Automatica*, 48(6):1088–1096, 2012.

[9] Shuo Feng, Yi Zhang, Shengbo Eben Li, Zhong Cao, Henry X. Liu, and Li Li. String stability for vehicular platoon control: Definitions and analysis methods. *Annual Reviews in Control*, 47:81–97, 2019.

[10] Yangyang Feng, Shuyou Yu, Hao Chen, Yongfu Li, Shuming Shi, Jianhua Yu, and Hong Chen. Distributed mpc of vehicle platoons with guaranteed consensus and string stability. *Scientific Reports*, 13:1–21, 06 2023.

[11] Bernardo Hernández and Paul Trodden. Distributed model predictive control using a chain of tubes. In *2016 UKACC 11th International Conference on Control*, pages 1–6, 2016.

[12] Ernesto Kofman, José A. De Doná, and María M. Seron. Probabilistic set invariance and ultimate boundedness. *Automatica*, 48(10):2670–2676, 2012.

[13] Vu Tuan Hieu Le, Cristina Stoica, Teodoro Alamo, Eduardo Camacho, and Didier Dumur. *Zonotopes: From Guaranteed State-estimation to Control*. John Wiley & Sons, 2013.

[14] Jan Löfberg. YALMIP: A Toolbox for Modeling and Optimization in MATLAB. In *In Proceedings of the CACSD Conference*, Taipei, Taiwan, 2004.

[15] José Maestre and Rudy Negenborn. *Distributed Model Predictive Control Made Easy*. Intelligent Systems, Control and Automation: Science and Engineering. Springer, 2014.

[16] D.Q. Mayne, M.M. Seron, and S.V. Raković. Robust model predictive control of constrained linear systems with bounded disturbances. *Automatica*, 41(2):219–224, 2005.

[17] Marco Mirabilio, Alessio Iovine, Elena De Santis, Maria Domenica Di Benedetto, and Giordano Pola. Scalable Mesh Stability of Nonlinear Interconnected Systems. *IEEE Control Systems Letters*, 6:968–973, 2022.

[18] Saša V. Raković, Basil Kouvaritakis, Rolf Findeisen, and Mark Cannon. Homothetic tube model predictive control. *Automatica*, 48(8):1631–1638, 2012.

[19] James Rawlings, David Mayne, and Moritz Diehl. *Model Predictive Control: Theory, Computation, and Design, Second Edition*. Nob Hill Publishing, 2020.
[20] Şerban Sabău, Cristian Oară, Sean Warnick, and Ali Jadbabaie. Optimal Distributed Control for Platooning via Sparse Coprime Factorizations. *IEEE Transactions on Automatic Control*, 62(1):305–320, 2017.

[21] Şerban Sabău, Andrei Sperilă, Cristian Oară, and Ali Jadbabaie. Network Realization Functions for Optimal Distributed Control. *IEEE Transactions on Automatic Control*, 68(12):8059–8066, 2023.

[22] Riccardo Scattolini. Architectures for distributed and hierarchical Model Predictive Control – A review. *Journal of Process Control*, 19(5):723–731, 2009.

[23] Andrei Sperilă, Alessio Iovine, Sorin Olaru, and Patrick Panciatici. Network-Realized Model Predictive Control – Part I: NRF-Enabled Closed-loop Decomposition, 2025. Available Online: https://arxiv.org/abs/2502.13042.

[24] Andrei Sperilă, Cristian Oară, Bogdan Ciubotaru, and Şerban Sabău. Distributed Control of Descriptor Networks: A Convex Procedure for Augmented Sparsity. *IEEE Transactions on Automatic Control*, 68(12):8067–8074, 2023.

[25] A. Sánchez-Amores, J. M. Maestre, P. A. Trodden, and E. F. Camacho. A Bound on the Existence of the Maximum Jointly Invariant Set of Input-Coupled Systems. *IEEE Control Systems Letters*, 7:2293–2298, 2023.

[26] Paul A. Trodden and J.M. Maestre. Distributed predictive control with minimization of mutual disturbances. *Automatica*, 77:31–43, 2017.

[27] Ellen van Nunen, Joey Reinders, Elham Semsar-Kazerooni, and Nathan van de Wouw. String Stable Model Predictive Cooperative Adaptive Cruise Control for Heterogeneous Platoons. *IEEE Transactions on Intelligent Vehicles*, 4(2):186–196, 2019.

[28] Jiange Wang, Xiaolei Li, Ju H. Park, and Ge Guo. Distributed MPC-Based String Stable Platoon Control of Networked Vehicle Systems. *IEEE Transactions on Intelligent Transportation Systems*, 24(3):3078–3090, 2023.

[29] Yuh-Shyang Wang, Nikolai Matni, and John Doyle. A System-Level Approach to Controller Synthesis. *IEEE Transactions on Automatic Control*, 64(10):4079–4093, 2019.

[30] Yang Zheng, Luca Furieri, Antonis Papachristodoulou, Na Li, and Maryam Kamgarpour. On the Equivalence of Youla, System-Level, and Input-Output Parameterizations. *IEEE Transactions on Automatic Control*, 66(1):413–420, 2021.

[31] Yang Zheng, Shengbo Eben Li, Keqiang Li, Francesco Borrelli, and J. Karl Hedrick. Distributed Model Predictive Control for Heterogeneous Vehicle Platoons Under Unidirectional Topologies. *IEEE Transactions on Control Systems Technology*, 25(3):899–910, 2017.

[32] Kemin Zhou, John Doyle, and Keith Glover. *Robust and Optimal Control*. Prentice-Hall, 1996.

## A Summary of NRF-based Control

We include here several key definitions and properties associated with the class of control laws from [21].

Fundamentally, the concept of an NRF-based representation arises when attempting to represent the TFM of a centralised controller, denoted henceforth as $K(z)$ and assumed to be of dimension $n_u \times n_x$, via the expression

$$K(z) = (I - \Phi(z))^{-1} \Gamma(z),$$

for which $\det(I - \Phi(z)) \neq 0$ and $\mathrm{elm}_{ii}(\Phi(z)) = 0$, for all $i \in \{1: n_u\}$. Thus, the TFM pair $(\Phi(z), \Gamma(z))$ allows for

the computation of the signal $u[k]$ in (1a)-(1c) as follows

$$u[k] = \Phi(z) \star u[k] + \Gamma(z) \star x[k]. \quad (A.1)$$

Since the diagonal elements of $\Phi(z)$ are all identically zero this enables a distinct feedback-feedforward expression of the resulting control laws, in which an arbitrary command signal $u_i[k]$ may be computed via the use of a feedback component $x[k]$ alongside the feedforward of the command signals $u_j[k]$ belonging to other nodes in the network ($j \neq i$). It is precisely this feature which, when expressed via (A.1), yields the identities from (12).

In the NRF framework, the problem of imposing sparsity on $K(z)$ for a distributed implementation (a highly restrictive objective, see [21] and [24]) is relaxed into ensuring sparsity for its NRF pair $(\Phi(z), \Gamma(z))$. The latter goal is far more tractable, due to the manner in which the pair may be computed. To this end, we now introduce the notion of *Doubly Coprime Factorisation* (DCF) over a set $C_g$, which can be employed as a design parameter.

**Remark A.1** The choice of $C_g$ underlines several key properties of the resulting NRF-based control laws. For example, choosing $C_g \subseteq \{z \in \mathbb{C} \mid |z| < 1\}$ ensures that the effect of bounded exogenous disturbance on the signals arising in closed-loop configuration remains bounded.

With respect to the realisation from (1a)-(1c), we define $G_u(z) := (zI_{n_x} - A)^{-1}B_u$ and we say that any collection of eight TFM's ($N(z), M(z), X(z), Y(z), \tilde{N}(z), \tilde{M}(z), \tilde{X}(z), \tilde{Y}(z)$) is a DCF over $C_g$ of $G_u(z)$ if:

a) $M(z)$ along with $\tilde{M}(z)$ are invertible and they satisfy $G_u(z) = N(z)M^{-1}(z) = \tilde{M}^{-1}(z)\tilde{N}(z)$;

b) all eight TFM's are proper (taking $|z| \to \infty$ results in a matrix having only finite entries) and have no poles outside of $C_g$;

c) the eight TFM's satisfy the (Bézout-like) identity

$$\begin{bmatrix} \tilde{Y}(z) - \tilde{X}(z) \\ -\tilde{N}(z) - \tilde{M}(z) \end{bmatrix} \begin{bmatrix} M(z) & X(z) \\ N(z) & Y(z) \end{bmatrix} = \begin{bmatrix} I_{n_u} & O \\ O & I_{n_x} \end{bmatrix}.$$

Using just such a DCF over $C_g$, it is possible to obtain a proper (and, thus, implementable) NRF pair for a controller of the system in (1a)-(1c) as follows

$$\left\{ \begin{array}{l} \tilde{Y}_{\mathbf{Q}}(z) := \tilde{Y}(z) + \mathbf{Q}(z)\tilde{N}(z), \\ \tilde{X}_{\mathbf{Q}}(z) := \tilde{X}(z) + \mathbf{Q}(z)\tilde{M}(z), \\ \tilde{Y}_{\mathbf{Q}}^{\mathrm{diag}}(z) := \mathrm{diag}(\mathrm{elm}_{11}(\tilde{Y}_{\mathbf{Q}}(z)), \dots, \mathrm{elm}_{n_u n_u}(\tilde{Y}_{\mathbf{Q}}(z))), \\ \mathbf{\Phi}(z) := I_{n_x} - (\tilde{Y}_{\mathbf{Q}}^{\mathrm{diag}}(z))^{-1}\tilde{Y}_{\mathbf{Q}}(z), \\ \mathbf{\Gamma}(z) := (\tilde{Y}_{\mathbf{Q}}^{\mathrm{diag}}(z))^{-1}\tilde{X}_{\mathbf{Q}}(z). \end{array} \right.$$

for a TFM $\mathbf{Q}(z)$ that is proper, that has no poles outside of $C_g$, and which ensures that $\tilde{Y}_{\mathbf{Q}}(z)$ and $\tilde{Y}_{\mathbf{Q}}^{\mathrm{diag}}(z)$ have proper inverses. For systems of type (1a)-(1c), the DCF over $C_g$ can always be chosen so that such a $\mathbf{Q}(z)$ exists.

One of the main features of this representation is the fact that, by computing the NRF pair as above, its two terms *inherit the sparsity patterns of* $\tilde{Y}_{\mathbf{Q}}(z)$ and $\tilde{X}_{\mathbf{Q}}(z)$.
$$ \text{row}_j(w_{r\ell}[k+1]) = -a_{j\ell} \text{row}_1(w_{r\ell}[k]) + \text{row}_{j+1}(w_{r\ell}[k]) + K_{j\ell} \left[ \begin{array}{c} u_f[k] + \beta_f[k] \\ x[k] + \zeta[k] + u_{s1}[k] + \beta_{s1}[k] \end{array} \right], \forall j \in \{1:n_{r\ell}-1\}, \quad (B.1a) $$

$$ \text{row}_{n_{r\ell}}(w_{r\ell}[k+1]) = -a_{n_{r\ell}\ell} \text{row}_1(w_{r\ell}[k]) + K_{n_{r\ell}\ell} \left[ \begin{array}{c} u_f[k] + \beta_f[k] \\ x[k] + \zeta[k] + u_{s1}[k] + \beta_{s1}[k] \end{array} \right], \quad (B.1b) $$

$$ \text{row}_{\ell}(u_f[k]) = \text{row}_1(w_{r\ell}[k]), \forall k \ge k_0, \ell \in \mathcal{A}_{ui}. \quad (B.1c) $$

$$ x_i[k_0+1] = C_{xi}\xi_i[k_0+1] + (\psi_{xi}[k_0+1] + \theta_{xi}[k_0+1] + \delta_{xi}[k_0+1]), \quad \psi_{xi}[k_0+1] \in \Psi_{xi1}, \theta_{xi}[k_0+1] \in \Theta_{xi1}, \delta_{xi}[k_0+1] \in \Delta_{xi1}, \quad (B.2a) $$

$$ \begin{aligned} \text{row}_j(u_{fi}[k_0+1]) ={}& \text{row}_j(C_{ui})\xi_i[k_0+1] + (\text{row}_j(\psi_{ui}[k_0+1]) + \text{row}_j(\theta_{ui}[k_0+1]) + \text{row}_j(\delta_{ui}[k_0+1))), \quad \text{row}_j(\psi_{ui}[k_0+1]) \in (e_j^\top \Psi_{ui1}), \\ & \text{row}_j(\theta_{ui}[k_0+1]) \in (e_j^\top \Theta_{ui1}), \quad \text{row}_j(\delta_{ui}[k_0+1]) \in (e_j^\top \Delta_{ui1}), \quad \forall j \in \{1:n_{ui}\}. \end{aligned} \quad (B.2b) $$

Indeed, since the two TFMs are affine terms of $Q(z)$, ensuring desired sparsity patterns on these expressions is a particularly tractable problem (see [24], for example). A further benefit of employing this control formalism is the fact that all closed-loop maps from exogenous disturbance to the interconnection’s internal signals (recall Remark A.1) are similarly simple expressions of the $Q(z)$ parameter. For example, the TFM $\mathbf{F}_Q(z) [I_{(n_x+n_u)} O]^\top$ from (14) can be expressed as follows

$$ \mathbf{F}_Q(z) \begin{bmatrix} I_{(n_x+n_u)} \\ O \end{bmatrix} = \begin{bmatrix} \mathbf{N}(z)\tilde{\mathbf{X}}_Q(z) & \mathbf{N}(z)\tilde{\mathbf{Y}}_Q(z) \\ \mathbf{M}(z)\tilde{\mathbf{X}}_Q(z) & \mathbf{M}(z)\tilde{\mathbf{Y}}_Q(z) - I_{n_u} \end{bmatrix}, $$

with the other TFMs from (14) being given explicitly in the main result of our companion paper [23]. Consequently, systemic norm-based procedures (such as the ones proposed in [24]) can be used to manipulate these amenable expressions, with the aim of decoupling the areas of the network from (1a)-(1c), as per Remark 4.4.

## B Proofs and Auxiliary Results

We begin by proving the result, which enables us to constrain the state variables of the first-layer subcontrollers.

### Proof of Proposition 5.1

The result follows by straightforward algebraic manipulations of the state-space dynamics given in (10). For the sake of clarity, we write these equations out in an explicit manner and we obtain the identities from (B.1a)-(B.1c), which are located at the top of this page.

First, we tackle the inclusions showcased in point i). Assumption A2 implies the fact that $\text{row}_1(w_{r\ell}[k_0]) \in W_{\ell 1}$ and, due to Assumption A1), using (34a) and (B.1b) yields the inclusion $\text{row}_{n_{r\ell}}(w_{r\ell}[k_0+1]) \in W_{\ell n_{r\ell}}$. Moreover, by employing the same arguments for (34b) and by effecting successive substitutions in (B.1a), we retrieve the fact that $\text{row}_j(w_{r\ell}[k_0+1]) \in W_{\ell j}$, for all $j \in \{2:n_{r\ell}-1\}$. To get that $\text{row}_1(w_{r\ell}[k_0+1]) \in W_{\ell 1}$, we simply use Assumption A3) in conjunction with (B.1c). To retrieve point ii), recall Assumption A2) and then use (B.1c) along with (32) to get $u_{fi}[k_0] \in U_{fi}$. The desired inclusion follows from Assumption A1), (29) and (31).

Moving on, we now prove the result which represents the cornerstone of our constraint management technique.

### Proof of Proposition 6.1

The proof is constructive in nature and, to begin with, notice that Assumption A2) enables the local computation of the second-layer command pairs $u_{s1i}[k_0] \in U_{s1i}$ and $u_{s2i}[k_0] \in U_{s2i}$, for each $i \in \{1:N\}$, as in **Step 3** of Algorithm 1. Since the algebraic conditions from (13) are satisfied by each of the first layer’s subcontrollers, then there is no need to pass on the values of $u_{s1j}[k_0]$, for any $j \in \{1:N\} \setminus N_i$, to the $i^{th}$ subcontroller from the first layer. Therefore, it follows that by communicating and by applying these partial command signals which originate from the second layer, as indicated in **Steps 4a and 4b** of Algorithm 1, then the implemented command scheme is equivalent to the one depicted in Figure 2. Consequently, the implementation described in Algorithm 1 ensures that the closed-loop dynamics of our control scheme’s first layer are those in (16a)-(16d). By applying the above-mentioned commands, the $i^{th}$ area’s state and its first-layer subcontroller’s outputs at time $k = k_0+1$ will be given by the identities showcased in (B.2a)-(B.2b), located at the top of this page, where

$$ \xi_i[k_0+1] = B_{s1i}u_{s1i}[k_0] + B_{s2i}u_{s2i}[k_0], \quad (B.3) $$

for all $i \in \{1:N\}$. Allowing for the state measurements received by each second-layer subcontroller to be subject, at $k = k_0$, to a disturbance of the type described in (20), it follows due to Assumption A1) and to (38) that these measurements can be expressed as in (39a)-(39b).

By employing now the identities from (21a)-(21b) and (24), it follows that the variables described in (B.2a)-(B.2b) can be rewritten as in (B.4a)-(B.4b), located at the top of the next page. However, since $u_{s1i}[k_0]$ and $u_{s2i}[k_0]$ are obtained by solving (8), it follows that the vectors located on the left-hand side of (B.3) must satisfy

$$ \xi_i[k_0+1] \in \Xi_{i1}, \forall i \in \{1:N\}. $$

By considering the expressions of the sets defined in (42)-(44), it follows by simple algebraic manipulations that

$$ \left\{ \begin{array}{l} (C_{xi}\xi_i[k_0+1] + \tilde{\theta}_{xi}[k_0+1]) \in \mathcal{P}_{1i1}, \\ (\text{row}_j(C_{ui})\xi_i[k_0+1] + \text{row}_j(\tilde{\theta}_{ui}[k_0+1])) \in \mathcal{P}_{2(j+\alpha_{ui})1}, \end{array} \right. $$

for all $i \in \{1:N\}$, along with all $j \in \{1:n_{ui}\}$. Moreover, by using the identities which form (B.4a)-(B.4b) in tandem with the set inclusions from (42)-(43) and with the
$$x_i[k_0 + 1] = \left( (C_{xi}\xi_i[k_0 + 1] + \tilde{\theta}_{xi}[k_0 + 1]) + (\psi_{xi}[k_0 + 1] + \tilde{\delta}_{xi}[k_0 + 1]) \right) + (-\eta_{xi}[k_0 + 1]), \\ \psi_{xi}[k_0 + 1] \in \Psi_{xi1}, \tilde{\theta}_{xi}[k_0 + 1] \in \tilde{\Theta}_{xi1}, \tilde{\delta}_{xi}[k_0 + 1] \in \tilde{\Delta}_{xi1}, \eta_{xi}[k_0 + 1] \in \mathcal{H}_{xi1}, \quad (\text{B.4a})$$

$$\begin{aligned} \text{row}_j(u_{fi}[k_0 + 1]) = & \left( (\text{row}_j(C_{ui})\xi_i[k_0 + 1] + \text{row}_j(\tilde{\theta}_{ui}[k_0 + 1])) + (\text{row}_j(\psi_{ui}[k_0 + 1]) + \text{row}_j(\tilde{\delta}_{ui}[k_0 + 1])) \right) + (-\text{row}_j(\eta_{ui}[k_0 + 1))), \\ & \text{row}_j(\psi_{ui}[k_0 + 1]) \in (e_j^\top \Psi_{ui1}), \text{row}_j(\tilde{\theta}_{ui}[k_0 + 1]) \in (e_j^\top \tilde{\Theta}_{ui1}), \text{row}_j(\tilde{\delta}_{ui}[k_0 + 1]) \in (e_j^\top \tilde{\Delta}_{ui1}), \\ & \text{row}_j(\eta_{ui}[k_0 + 1]) \in (e_j^\top \mathcal{H}_{ui1}), \forall j \in \{1: n_{ui}\}. \end{aligned} \quad (\text{B.4b})$$

set definitions introduced in (36) and (37), we get that $x_i[k_0 + 1] \in \mathcal{X}_i$ and that $\text{row}_j(u_{fi}[k_0 + 1]) \in \mathcal{W}_{(j+\alpha_{ui})1}$ for all $i \in \{1: N\}$ and all $j \in \{1: n_{ui}\}$.

The first conclusion we are able to draw, at this point, is the fact that $x[k_0 + 1] \in \mathcal{X}$, by virtue of (26). Now, in order to retrieve the remaining inclusions given in our result's statement, notice that all of the assumptions associated with Proposition 5.1 hold. Thus, by applying said result, we get that $\text{row}_j(w_{r\ell}[k_0 + 1]) \in \mathcal{W}_{\ell j}$, for all $i \in \{1: N\}$, $\ell \in \mathcal{A}_{ui}$ and $j \in \{1: n_{r\ell}\}$, along with the fact that $u_i[k_0] \in \mathcal{U}_i, \forall i \in \{1: N\}$. Finally, the latter inclusions yield $u[k_0] \in \mathcal{U}$, by way of (26).

We now proceed to employ the arguments made in the proof of Proposition 6.1, in order to obtain guarantees of recursive feasibility for our architecture's second layer.

### Proof of Theorem 6.1

The central idea, which underpins the entire proof, is to show that all of the problems described in (8) are guaranteed to be feasible at an arbitrary time instant $k_0$, and considering some prediction horizons $T_i > 0$ and $\bar{T}_i \ge 0$. Once this is achieved, the proof boils down to a series of induction-based arguments in which Proposition 6.1 is repeatedly applied, with the aim of ensuring the sought-after inclusions from the result's statement.

We begin by choosing an arbitrary area index $i \in \{1: N\}$, and by investigating its associated area of the network. Due to Assumption A2), we have the fact that (46) holds for some $\rho_i \ge 1$. Thus, for any $\tilde{\theta}_{it} \in \tilde{\Theta}_{it}$ and any $t \in \{1: \rho_i\}$, there must exist $p_{it} \in \tilde{\mathcal{P}}_{it}$, $u_{1it} \in \mathcal{U}_{s1i}$ along with $u_{2it} \in \mathcal{U}_{s2i}$ such that the following identities hold

$$\tilde{\theta}_{it} = \left( \left[ \begin{matrix} -C_{xi} \\ -C_{ui} \end{matrix} \right] \sum_{j=1}^{t} A_{si}^{j-1} (B_{s1i} u_{1ij} + B_{s2i} u_{1ij}) \right) + p_{it}. \quad (\text{B.5})$$

Now, by exploiting (45a) and (45b), it becomes possible to rewrite (B.5) in a more explicit manner, which yields

$$\left\{ \begin{array}{l} \tilde{\theta}_{xit} + C_{xi} \sum_{j=1}^{t} A_{si}^{(j-1)} [B_{s1i} \quad B_{s2i}] \begin{bmatrix} u_{1ij} \\ u_{2ij} \end{bmatrix} = p_{xit}, \\ \\ \tilde{\theta}_{u\ell t} + \text{row}_{(\ell-\alpha_{ui})}(C_{ui}) \sum_{j=1}^{t} A_{si}^{(j-1)} [B_{s1i} \quad B_{s2i}] \begin{bmatrix} u_{1ij} \\ u_{2ij} \end{bmatrix} = p_{u\ell t}, \end{array} \right. \quad (\text{B.6})$$

for all $t \in \{1: \rho_i\}$ and all $\ell \in \mathcal{A}_{ui}$, where $\tilde{\theta}_{xit} \in \tilde{\Theta}_{xit}$, $\tilde{\theta}_{u\ell t} \in (e_{(\ell-\alpha_{ui})}^\top \tilde{\Theta}_{uit})$, $p_{xit} \in \mathcal{P}_{1it}$ and $p_{u\ell t} \in \mathcal{P}_{2\ell t}$. We now employ the expressions deduced in (B.6) in order to prove the recursive feasibility of the problems formulated in (8). To do so, choose any desired initial time

instant $k_0$ for the optimisation problem associated with the $i^{th}$ area index. For all $t \in \{1: \rho_i\}$, consider any $\tilde{\theta}_{xi}[k_0 + t] \in \tilde{\Theta}_{xit}$ and any $\tilde{\theta}_{ui}[k_0 + t] \in \tilde{\Theta}_{uit}$ which may arise in the context of the aforementioned optimisation problem, in order to define $\tilde{\theta}_{xit} := \tilde{\theta}_{xi}[k_0 + t]$ along with $\tilde{\theta}_{u\ell t} := \text{row}_{(\ell-\alpha_{ui})}(\tilde{\theta}_{ui}[k_0 + t])$, with $\ell \in \mathcal{A}_{ui}$. Due to the arguments made in the previous paragraph of this proof, it follows that there must exist some appropriately constrained vectors denoted $u_{1it}$ and $u_{2it}$ along with $p_{xit}$ and $p_{u\ell t}$, with $t \in \{1: \rho_i\}$ and $\ell \in \mathcal{A}_{ui}$, such that the identities from (B.6) hold. Therefore, by setting $T_i = \rho_i$ and by selecting an arbitrary $\bar{T}_i \in \mathbb{N}$, it is always possible to choose the command signal vectors which appear in (8) as $u_{s1i}[k_0 + t - 1] := u_{1ij} \in \mathcal{U}_{s1i}$ and as $u_{s2i}[k_0 + t - 1] := u_{2ij} \in \mathcal{U}_{s2i}$, for all $t \in \{1: T_i\}$. A direct consequence of all these choices is the fact that

$$\xi_i[k_0 + t] = \sum_{j=1}^{t} A_{si}^{(j-1)} (B_{s1i} u_{1ij} + B_{s2i} u_{2ij}), \quad (\text{B.7})$$

for all $t \in \{1: T_i\}$. When employed in tandem with (B.6), the identities obtained in (B.7) yield the fact that

$$\left\{ \begin{array}{l} C_{xi}\xi_i[k_0 + t] + \tilde{\theta}_{xi}[k_0 + t] \in \mathcal{P}_{1it}, \\ \text{row}_{(\ell-\alpha_{ui})}(C_{ui})\xi_i[k_0 + t] + \text{row}_{(\ell-\alpha_{ui})}(\tilde{\theta}_{xi}[k_0 + t]) \in \mathcal{P}_{2\ell t}, \end{array} \right. \quad (\text{B.8})$$

for all $t \in \{1: T_i\}$ and all $\ell \in \mathcal{A}_{ui}$. It follows by direct inspection of the inequalities used to define the sets introduced in (42)-(44) that the inclusions obtained in (B.8) are equivalent to $\xi_i[k_0 + t] \in \Xi_{it}$, for all $t \in \{1: T_i\}$.

At this point in the proof, we have shown that by setting $T_i = \rho_i$ and by choosing an arbitrary $\bar{T}_i \in \mathbb{N}$, there exists at least one sequence of command vectors which satisfies the constraints of the problems formulated in (8), thus rendering them feasible. Moreover, recall that the time instant $k_0$ and the area index $i \in \{1: N\}$ were chosen arbitrarily at the beginning of the proof, thus making the feasibility of all the aforementioned problems a generic property. By recalling now Assumption A1), we employ Proposition 6.1 to compute the second-layer command values at time instant $k = k_0$ and, by iterating the procedure detailed in Algorithm 1, we arrive at time instant $k = k_0 + 1$, in which Assumption A1) is still valid and in which all $N$ local optimisation problems remain feasible. It follows that this procedure may be successfully repeated at any subsequent instant $k \ge k_0$, thus guaranteeing recursive feasibility, i.e., point i) of our result, while also ensuring the satisfaction of the constraints which make up point ii) of our result.