Scalable Forward Reachability Analysis of Multi-Agent Systems with
  Neural Network Controllers
Neural networks (NNs) have been shown to learn complex control laws
successfully, often with performance advantages or decreased computational cost
compared to alternative methods. Neural network controllers (NNCs) are,
however, highly sensitive to disturbances and uncertainty, meaning that it can
be challenging to make satisfactory robustness guarantees for systems with
these controllers. This problem is exacerbated when considering multi-agent
NN-controlled systems, as existing reachability methods often scale poorly for
large systems. This paper addresses the problem of finding overapproximations
of forward reachable sets for discrete-time uncertain multi-agent systems with
distributed NNC architectures. We first reformulate the dynamics, making the
system more amenable to reachablility analysis. Next, we take advantage of the
distributed architecture to split the overall reachability problem into smaller
problems, significantly reducing computation time. We use quadratic
constraints, along with a convex representation of uncertainty in each agent's
model, to form semidefinite programs, the solutions of which give
overapproximations of forward reachable sets for each agent. Finally, the
methodology is tested on two realistic examples: a platoon of vehicles and a
power network system.
Oliver Gates, Matthew Newton, Konstantinos Gatsis
# Scalable Forward Reachability Analysis of Multi-Agent Systems with Neural Network Controllers

###### Abstract

Neural networks (NNs) have been shown to learn complex control laws successfully, often with performance advantages or decreased computational cost compared to alternative methods. Neural network controllers (NNCs) are, however, highly sensitive to disturbances and uncertainty, meaning that it can be challenging to make satisfactory robustness guarantees for systems with these controllers. This problem is exacerbated when considering multi-agent NN-controlled systems, as existing reachability methods often scale poorly for large systems. This paper addresses the problem of finding overapproximations of forward reachable sets for discrete-time uncertain multi-agent systems with distributed NNC architectures. We first reformulate the dynamics, making the system more amenable to reachability analysis. Next, we take advantage of the distributed architecture to split the overall reachability problem into smaller problems, significantly reducing computation time. We use quadratic constraints, along with a convex representation of uncertainty in each agent's model, to form semidefinite programs, the solutions of which give overapproximations of forward reachable sets for each agent. Finally, the methodology is tested on two realistic examples: a platoon of vehicles and a power network system.

## I Introduction

There has been recent interest in the use of neural networks (NNs) for control in closed-loop feedback systems. Neural network controllers (NNCs) can be used to imitate traditional control policies, such as model predictive control (MPC), with reduced computational cost [1], or to implement deep reinforcement learning (RL) policies [2]. Even for simple linear systems, NNCs can be used to implement complex non-linear control laws (which may not be easy to achieve with existing methods). NNs are, however, highly sensitive to input perturbations, so disturbances in the closed-loop system can have adverse effects [3]. This is problematic when NNCs are used in safety-critical systems, and recent work has focused on reachability analysis for systems with NNCs; if we can overapproximate the forward reachable sets, then it can be verified that certain regions of the state space are avoided over a given horizon.

The problem of computing forward reachable sets becomes more challenging when considering a multi-agent system in which each agent is controlled by an NN (or a series of NNs); the effects of small perturbations to the input of one NNC are propagated through the system. Control of multi-agent systems is well-studied, and common goals include consensus, formation control and flocking/swarming [4]. Multi-agent control architectures can be categorised according to the dependence of each agent's control input on other agents' states: (a) centralised control, in which each agent's control input is a function of all agents' states; (b) distributed control, in which each agent's control input is a function of a subset of the other agents' states; (c) decentralised control, in which each agent's control input is a function of only its own state [5]. In distributed control, common approaches include state feedback [6] and distributed model predictive control (DMPC) [7].

A number of methods have been proposed for the reachability analysis of systems with NNCs [8, 9, 10, 11, 12, 13, 14, 15, 16]. In general, the problem of reachability for discrete-time LTI systems is undecidable [17], so methods for safety verification of closed-loop systems often aim to find tight overapproximations of the forward reachable sets. Alternatively, additional restrictions can be placed on the problem to allow the exact reachable sets to be computed. Generally, there is a tradeoff between scalability and tightness of the bounds [3].

In [8], semidefinite programming (SDP) is used to compute overapproximations of the forward reachable sets, and [9] builds on this work by considering parameter-varying systems. In [10] and [11], the input set is partitioned into smaller sets, and a linear programming (LP) approach is used to overapproximate the reachable sets; in [10], the input set partitioning approach is also applied to the method in [8], and a comparison is made between the LP and SDP approaches, demonstrating that the former is faster but the latter results in tighter bounds. The work in [12] restricts the input sets to be constrained zonotopes, allowing for the exact computation of reachable sets, and the work in [13] represents the input sets as hybrid zonotopes, allowing for a class of non-convex input sets. In [14], polynomial zonotopes are used to abstract the closed-loop dynamics, providing tight overapproximations. Other approaches include the use of polynomial optimisation [15] and Bernstein polynomials [16].

Existing reachability methods for NN-controlled systems do not explicitly consider multi-agent systems with distributed control architectures. Similarly to the single-agent case, in which NNCs can be used to learn complex control laws, a distributed NNC architecture can be used to learn complex distributed control laws [18]. An example is DMPC, in which each agent's controller solves an optimisation problem based on its own state and those of its neighbours. An overview of DMPC is given in [7], and its applications
include vehicle platooning [19], frequency regulation in power systems [20] and formation control of UAVs [21].

In this paper, we present a scalable method to compute overapproximations of the forward reachable sets for uncertain multi-agent systems with distributed NNC architectures. To the best of our knowledge, this is the first work which explicitly deals with the multi-agent NNC reachability problem. The main contributions are as follows:

* we reformulate the dynamics, making the system more amenable to reachability analysis;
* we take advantage of the distributed architecture to split the overall reachability problem into smaller problems, using an existing SDP-based approach to overapproximate the forward reachable sets for each agent, and further extend this result to incorporate model uncertainty in the agents' dynamics;
* we demonstrate the effectiveness of this methodology on two realistic multi-agent systems with different structures: a platoon of vehicles and a power network;
* we compare our approach to the approach of treating the multi-agent system as one overall system, and we demonstrate that our approach outperforms the alternative approach.

In Section II, we describe the dynamics, give some remarks about the form of the controller and describe the forward reachability problem. In Section III, we provide a simplification of the dynamics and control input, and in Section IV, we present the reachability method. In Section V, we introduce model uncertainty. In Section VI, we present experiments to demonstrate the method.

### _Notation_

The set of real \(n\times m\) matrices is denoted by \(\mathbb{R}^{n\times m}\), the set of real \(n\)-length vectors by \(\mathbb{R}^{n}\) and the set of real numbers by \(\mathbb{R}\). The set of symmetric \(n\times n\) matrices is denoted by \(\mathbb{S}^{n}\). The set of diagonal \(n\times n\) matrices is denoted by \(\mathbb{D}^{n}\). The set of positive integers is denoted by \(\mathbb{Z}^{+}\). The cardinality of a set \(\mathcal{S}\) is denoted by \(|\mathcal{S}|\). The \(n\times n\) identity matrix is denoted by \(I_{n}\). The symbols \(\geq\) and \(\leq\) apply elementwise to vectors and matrices. \(A\preceq 0\) implies that matrix \(A\) is negative semidefinite. The number \(0\) is used to represent the scalar, vector or matrix of appropriate size; the size will be clear from the context.

For clarity, whenever the letters \(i\) and \(j\) are used in this paper, they refer to the index of an agent.

## II Problem Statement

### _Multi-agent dynamics_

We consider a discrete-time multi-agent system of \(M\) agents, where each agent \(i\) has linear time-invariant dynamics

\[x_{k+1}^{[i]}=A_{ii}x_{k}^{[i]}+\sum_{j\in\mathcal{N}_{i}}A_{ij}x_{k}^{[j]}+B_ {i}u_{k}^{[i]}+w_{k}^{[i]}, \tag{1}\]

where for each agent \(i\in\{1,\ldots,M\}=\mathcal{I}\), \(x_{k}^{[i]}\in\mathbb{R}^{n_{x}}\) is the local state, \(x_{k}^{[j]}\in\mathbb{R}^{n_{x}}\) is the \(j^{\text{th}}\) neighbouring state, \(u_{k}^{[i]}\in\mathbb{R}^{n_{x}}\) is the control input, \(w_{k}^{[i]}\in\mathbb{R}^{n_{x}}\) is a known external input, \(A_{ii}\in\mathbb{R}^{n_{x}\times n_{x}}\) is the state matrix, \(A_{ij}\in\mathbb{R}^{n_{x}\times n_{x}}\) is the matrix describing the effect of state \(x_{k}^{[j]}\) on agent \(i\), \(B_{i}\in\mathbb{R}^{n_{x}\times n_{x}}\) is the input matrix, and \(\mathcal{N}_{i}\) is the set of neighbours. Note that for simplicity, we have assumed that all agents have the same dimensions (although this assumption could be relaxed). We also assume that \(|\mathcal{N}_{i}|>0\)\(\forall i\).

In this paper, without loss of generality, we focus on distributed control architectures, in which each agent's control input is a function of a subset \(\mathcal{N}_{i}\) of all other agents' states \(\mathcal{I}\). The methods presented in this paper can be easily extended to the other two cases by setting \(\mathcal{N}_{i}=\mathcal{I}\) (centralised) or \(\mathcal{N}_{i}=\{i\}\) (decentralised).

### _Control input_

In a traditional distributed control scheme with proportional feedback, the control input \(u_{\mathrm{trad},k}^{[i]}\) might be given by

\[u_{\mathrm{trad},k}^{[i]}=\sum_{j\in\mathcal{N}_{i}}K_{ij}\left(x_{k}^{[j]}-x_ {k}^{[i]}\right),\]

where \(K_{ij}\in\mathbb{R}^{n_{u}\times n_{x}}\) is some gain matrix (which could represent multiple gain matrices) [22]. Similarly, in a DMPC scheme, the \(i^{\text{th}}\) control input is generated by solving an optimisation problem based on the agent's state \(x_{k}^{[i]}\) and the neighbours' states \(x_{k}^{[j]}\)\(\forall j\in\mathcal{N}_{i}\). In this paper, we consider the extension of the traditional distributed control schemes to NN-based control, in which the \(i^{\text{th}}\) control input is generated by some non-linear function of the agent's state \(x_{k}^{[i]}\) and the neighbours' states \(x_{k}^{[j]}\)\(\forall j\in\mathcal{N}_{i}\). We also consider the possibility of controller saturation. Hence, the \(i^{\text{th}}\) control input \(u_{k}^{[i]}\) is given by

\[u_{k}^{[i]}=\mathrm{sat}_{\mathcal{U}_{i}}\left[\sum_{j\in\mathcal{N}_{i}} \pi_{ij}\left(\begin{bmatrix}x_{k}^{[i]}\\ x_{k}^{[j]}\end{bmatrix}\right)\right], \tag{2}\]

where \(\mathrm{sat}_{\mathcal{U}_{i}}\) is a projection onto the set \(\mathcal{U}_{i}=\{u\in\mathbb{R}^{n_{u}}\mid\underline{u}_{i}\leq u\leq \overline{u}_{i}\}\), where \(\underline{u}_{i}\) and \(\overline{u}_{i}\) are the lower and upper limits, respectively, for the \(i^{\text{th}}\) controller, and \(\pi_{ij}:\mathbb{R}^{2n_{x}}\rightarrow\mathbb{R}^{n_{u}}\) is a function representing the mapping of the input through a multi-layer perceptron (MLP).

**Remark 1**: _In (2), we consider a separate NN \(\pi_{ij}\) for each neighbouring agent, as this preserves the ability to consider the effects of each neighbour individually, by isolating the effect of a particular neighbour's contribution to the control input. However, the control input could also be generated by feeding the state of the agent and the states of its neighbours into a single MLP \(\Pi_{i}:\mathbb{R}^{(1+|\mathcal{N}_{i}|)n_{x}}\rightarrow\mathbb{R}^{n_{u}}\). A transformation between the two architectures is given in Section III-B._

### _Multi-layer perceptron_

The mapping \(s\mapsto\pi_{ij}(s)\) for an \(L\)-layer MLP is

\[z_{ij}^{0} =s, \tag{3a}\] \[z_{ij}^{\ell+1} =\sigma^{\ell}\left(W_{ij}^{\ell}z_{ij}^{\ell}+b_{ij}^{\ell} \right),\quad\ell=0,\ldots,L-1,\] (3b) \[\pi_{ij}(s) =W_{ij}^{L}z_{ij}^{\ell}+b_{ij}^{L}, \tag{3c}\]
where \(z_{ij}^{\ell}\in\mathbb{R}^{n_{\ell}}\) is the \(\ell^{\text{th}}\) vector of activation values (note that \(n_{0}=2n_{x}\)), \(W_{ij}^{\ell}\in\mathbb{R}^{n_{\ell+1}\times n_{\ell}}\) is the \(\ell^{\text{th}}\) weight matrix, \(b_{ij}^{\ell}\in\mathbb{R}^{n_{\ell+1}}\) is the \(\ell^{\text{th}}\) bias vector, and \(\sigma^{\ell}:\mathbb{R}^{n_{\ell+1}}\rightarrow\mathbb{R}^{n_{\ell+1}}\) is the \(\ell^{\text{th}}\) ReLU activation function, i.e. \(\sigma^{\ell}(y)=\max(y,0)\), applied elementwise, such that \(\sigma^{\ell}(y)=\begin{bmatrix}\max(y_{1},0)&\cdots&\max(y_{n_{\ell}+1},0) \end{bmatrix}^{\top}\), where \(y=\begin{bmatrix}y_{1}&\ldots&y_{n_{\ell}+1}\end{bmatrix}^{\top}\) is the vector of pre-activation values. For simplicity, we assume that all MLPs have the same structure (size and number of hidden layers) for all \(i,j\).

### _Overapproximation of forward reachable sets_

We denote the set of all possible states of the \(i^{\text{th}}\) agent at time \(k\) as \(\mathcal{X}_{k}^{[i]}\), such that \(x_{k}^{[i]}\in\mathcal{X}_{k}^{[i]}\). Given \(\mathcal{X}_{k}^{[i]}\) and the sets of neighbouring states at time \(k\), i.e. \(\mathcal{X}_{k}^{[j]}\)\(\forall j\in\mathcal{N}_{i}\), we aim to find an overapproximation \(\widehat{\mathcal{X}}_{k+1}^{[i]}\) of the reachable set \(\mathcal{X}_{k+1}^{[i]}\) at the next time step. Specifically, we aim to find the tightest polytopic overapproximation of the reachable set

\[\min \quad\operatorname{vol}\left(\widehat{\mathcal{X}}_{k+1}^{[i]}\right)\] (4a) subject to \[\widehat{\mathcal{X}}_{k+1}^{[i]}\supseteq\mathcal{X}_{k+1}^{[i]}, \tag{4b}\] \[\widehat{\mathcal{X}}_{k+1}^{[i]}\text{ is a polytope}, \tag{4c}\]

for each agent \(i\), where \(\operatorname{vol}\) is the \(n_{x}\)-dimensional volume. This could then be applied recursively to overapproximate the next forward reachable sets \(\mathcal{X}_{k+2}^{[i]},\ldots,\mathcal{X}_{k+N}^{[i]}\) over a finite horizon \(N\), and used to verify that certain unsafe regions of the state space are avoided over this horizon.

## III Reformulation of Dynamics

An obvious approach to the forward reachability problem is to augment the agents' states into one 'overall' state

\[x_{k}=\begin{bmatrix}x_{k}^{[1]}{}^{\top}&\cdots&x_{k}^{[M]}{}^{\top}\end{bmatrix}^ {\top},\]

then form a recursion for the overall dynamics and use existing methods to perform reachability analysis on this system. The main issue with this approach is that it ignores the distributed architecture of the system - the dynamics of each agent depend only on its neighbours, not on all other agents. As a result, we reformulate the dynamics given by (1) and (2) to allow reachability analysis to be performed for each agent. In Section VI, we show that there is a significant computational advantage to solving \(M\) smaller reachability problems over one large reachability problem for the SDP-based method. We note that the approach of decomposing the reachability problem into smaller problems has been used more generally in work on reachability analysis [23, 24].

### _Simplification of dynamics_

Let \(g(i,j)\) be a function which returns the \(j^{\text{th}}\) neighbour of the \(i^{\text{th}}\) agent. For example, if agent \(3\) has \(\mathcal{N}_{3}=\{2,5,6\}\), then \(g(3,1)=2\), \(g(3,2)=5\) and \(g(3,3)=6\). Then we can write \(\mathcal{N}_{i}=\{g(i,1),\ldots,g(i,q_{i})\}\), where \(q_{i}=|\mathcal{N}_{i}|\), and let

\[\tilde{A}_{i} =\begin{bmatrix}A_{ii}&A_{ig(i,1)}&\cdots&A_{ig(i,q_{i})}\end{bmatrix}, \tag{5}\] \[\tilde{x}_{k}^{[i]} =\begin{bmatrix}x_{k}^{[i]}{}^{\top}&x_{k}^{[g(i,1)]}{}^{\top}& \cdots&x_{k}^{[g(i,q_{i})]}{}^{\top}\end{bmatrix}^{\top}, \tag{6}\]

such that \(\tilde{x}_{k}^{[i]}\) is the concatenation of the state of agent \(i\) and the states of the neighbours of agent \(i\) at time \(k\), then (1) can be written as

\[x_{k+1}^{[i]}=\tilde{A}_{i}\tilde{x}_{k}^{[i]}+B_{i}u_{k}^{[i]}+w_{k}^{[i]}.\]

### _Simplification of control input_

To simplify (2), we aim to represent it as the mapping of \(\tilde{x}_{k}^{[i]}\) through a single MLP \(\Pi_{i}:\mathbb{R}^{(1+q_{i})n_{x}}\rightarrow\mathbb{R}^{n_{u}}\). Note that the layer sizes can differ from agent to agent, depending on \(|\mathcal{N}_{i}|\). The mapping from \(\tilde{x}_{k}^{[i]}\) to \(\Pi_{i}(\tilde{x}_{k}^{[i]})\) is then

\[z_{i,k}^{0} =\tilde{x}_{k}^{[i]}, \tag{7a}\] \[z_{i,k}^{\ell+1} =\sigma_{i}^{\ell}\left(W_{i}^{\ell}z_{i,k}^{\ell}+b_{i}^{\ell} \right),\quad\ell=0,\ldots,L-1,\] (7b) \[\Pi_{i}\left(\tilde{x}_{k}^{[i]}\right) =W_{i}^{L}z_{i,k}^{L}+b_{i}^{L}, \tag{7c}\]

where \(z_{i,k}^{\ell}\in\mathbb{R}^{\tilde{n}_{i}^{i}}\) is the \(\ell^{\text{th}}\) vector of activation values (note that \(\tilde{n}_{0}^{i}=(1+q_{i})n_{x}\)), \(\sigma_{i}^{\ell}:\mathbb{R}^{\tilde{n}_{i+1}^{i}}\rightarrow\mathbb{R}^{\tilde {n}_{i+1}^{i}}\) is the \(\ell^{\text{th}}\) ReLU activation function, and the weight matrices and bias vectors are given by \(W_{i}^{0}=\mathfrak{W}_{i}^{0}\Lambda_{i}\), \(b_{i}^{0}=\mathfrak{W}_{i}^{0}\), \(W_{i}^{\ell}=\mathfrak{W}_{i}^{\ell}\), \(b_{i}^{\ell}=\mathfrak{b}_{i}^{\ell}\) for \(\ell=1,\ldots,L-1\), and \(W_{i}^{L}=\Omega_{i}\mathfrak{W}_{i}^{\ell}\), \(b_{i}^{L}=\Omega_{i}\mathfrak{b}_{i}^{\ell}\), where

\[\mathfrak{W}_{i}^{\ell} =\operatorname{blkdiag}\left(W_{ig(i,1)}^{\ell},\ldots,W_{ig(i,q _{i})}^{\ell}\right),\] \[\mathfrak{b}_{i}^{\ell} =\begin{bmatrix}b_{ig(i,1)}^{\ell}&\cdots&b_{ig(i,q_{i})}^{\ell} \end{bmatrix}^{\top},\]

for \(\ell=0,\ldots,L\), where \(\Lambda_{i}\in\mathbb{R}^{2q_{i}n_{x}\times\tilde{n}_{0}^{i}}\) and \(\Omega_{i}\in\mathbb{R}^{n_{u}\times q_{i}n_{u}}\) are given by

\[\Lambda_{i}=\begin{bmatrix}I_{n_{x}}&0&0&0&\cdots&0\\ 0&I_{n_{x}}&0&0&\cdots&0\\ I_{n_{x}}&0&0&0&\cdots&0\\ 0&0&I_{n_{x}}&0&\cdots&0\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots\\ I_{n_{x}}&0&0&0&\cdots&0\\ 0&0&0&0&\cdots&I_{n_{x}}\end{bmatrix},\quad\Omega_{i}=\begin{bmatrix}I_{n_{u}} \\ \vdots\\ I_{n_{u}}\end{bmatrix}^{\top}.\]

Finally, we can write

\[\sum_{j\in\mathcal{N}_{i}}\pi_{ij}\left(\begin{bmatrix}x_{k}^{[i]}\\ x_{k}^{[j]}\end{bmatrix}\right)\equiv\Pi_{i}\left(\tilde{x}_{k}^{[i]}\right),\]

which is illustrated in Figure 1.

Figure 1: Reformulation of individual MLPs into one larger MLP
### _Summary of reformulation_

**Lemma 1**: _The dynamics given by (1) and (2) can be written as_

\[x_{k+1}^{[i]}=\tilde{A}_{i}\tilde{x}_{k}^{[i]}+B_{i}\mathrm{sat}_{\mathcal{U}_{i }}\left[\Pi_{i}\left(\tilde{x}_{k}^{[i]}\right)\right]+w_{k}^{[i]}, \tag{8}\]

_where \(\tilde{A}_{i}\) and \(\tilde{x}_{k}^{[i]}\) are defined in (5) and (6), respectively, and \(\Pi_{i}(\tilde{x}_{k}^{[i]})\) is given by (7a)-(7c)._

## IV Forward Reachability Analysis

To overapproximate the forward reachable sets, we extend the reachability method outlined in [3] and [8] for the multi-agent case. The method in [8], called Reach-SDP, uses quadratic constraints (QCs) and semidefinite programming; the input set, NN and reachable set are abstracted using QCs. This involves forming quadratic inequalities for sets and functions by pre- and post-multiplying a matrix by a 'basis vector'. Using 'change-of-basis' matrices allows the same basis vector to be used for all inequalities, which allows linear matrix inequalities (LMIs) to be formed for forward reachability analysis of the closed-loop system.

### _Incorporation of control limits into MLP_

We first add two layers to the MLP in (7a)-(7c) to account for the control limits [8]. Recall that \(\tilde{x}_{k}^{[i]}\) is the concatenation of the state of agent \(i\) and the states of the neighbours of agent \(i\) at time \(k\). The mapping from \(\tilde{x}_{k}^{[i]}\) to \(u_{k}^{[i]}\) is now defined by (7a), (7b) and

\[z_{i,k}^{L+1} =\sigma_{i}^{L}\left(W_{i}^{L}z_{i,k}^{L}+b_{i}^{L}-\underline{u} _{i}\right), \tag{9a}\] \[z_{i,k}^{L+2} =\sigma_{i}^{L+1}\left(-z_{i,k}^{L+1}+\overline{u}_{i}-\underline{ u}_{i}\right),\] (9b) \[u_{i}^{[i]} =-z_{i,k}^{L+2}+\overline{u}_{i}. \tag{9c}\]

We then define the overall basis vector for the QCs as \(\left[\underline{\mathbf{z}}_{i,k}\hskip 1.0pt^{\top}\hskip 1.0pt1\right]^{\top}\), where \(\underline{\mathbf{z}}_{i,k}^{\top}=\left[{z_{i,k}^{0}}\hskip 1.0pt^{\top} \hskip 1.0ptz_{i,k}^{1}\hskip 1.0pt^{\top}\hskip 1.0pt\cdots\hskip 1.0ptz_{i,k}^{L+2 ^{\top}}\right]\).

### _Input set_

We describe the input set \(\widetilde{\mathcal{X}}_{k}^{[i]}\), such that \(\tilde{x}_{k}^{[i]}\in\widetilde{\mathcal{X}}_{k}^{[i]}\), by a hyper-rectangle, i.e.

\[\widetilde{\mathcal{X}}_{k}^{[i]}=\left\{x\in\mathbb{R}^{\tilde{n}_{0}^{i}} \mid\underline{\tilde{x}}_{k}^{[i]}\leq x\leq\overline{\tilde{x}}_{k}^{[i]} \right\}, \tag{10}\]

where \(\underline{\tilde{x}}_{k}^{[i]},\overline{\tilde{x}}_{k}^{[i]}\in\mathbb{R}^{ \tilde{n}_{0}^{i}}\) are known bounds on the state of the agent and the states of its neighbours at time \(k\). Using [3, Definition 1] and [3, Proposition 1], and noting that the first activation value \(z_{i,k}^{0}\) is equal to \(\tilde{x}_{k}^{[i]}\), we can write (10) as

\[\begin{bmatrix}z_{i,k}^{0}\\ 1\end{bmatrix}^{\top}P_{k}^{i}(\Gamma)\begin{bmatrix}\tilde{z}_{i,k}^{0}\\ 1\end{bmatrix}\geq 0, \tag{11}\]

where

\[P_{k}^{i}(\Gamma)=\begin{bmatrix}-2\Gamma&\Gamma\left(\underline{\tilde{x}}_ {k}^{[i]}+\overline{\tilde{x}}_{k}^{[i]}\right)\\ \left(\underline{\tilde{x}}_{k}^{[i]}+\overline{\tilde{x}}_{k}^{[i]}\right)^{ \top}\Gamma&-2\underline{\tilde{x}}_{k}^{[i]}{}^{\top}\Gamma\overline{ \tilde{x}}_{k}^{[i]}\end{bmatrix},\]

where

\[\Gamma\in\mathbb{D}^{\tilde{n}_{0}^{i}}\]

and \(\Gamma\geq 0\). By using a change-of-basis matrix [8]

\[E_{\mathrm{in}}^{i}=\begin{bmatrix}I_{\tilde{n}_{0}^{i}}&0&\cdots&0&0\\ 0&0&\cdots&0&1\end{bmatrix},\]

we can write (11) as

\[\begin{bmatrix}\underline{\mathbf{z}}_{i,k}\\ 1\end{bmatrix}^{\top}\Delta_{k}^{i}(\Gamma)\begin{bmatrix}\underline{\mathbf{ z}}_{i,k}\\ 1\end{bmatrix}\geq 0, \tag{12}\]

where \(\Delta_{k}^{i}(\Gamma)={E_{\mathrm{in}}^{i}}^{\top}P_{k}^{i}(\Gamma)E_{\mathrm{ in}}^{i}\). This is summarised in the following lemma.

**Lemma 2**: _Consider the state of agent \(i\) and the states of the neighbours of agent \(i\) at time \(k\). If these are are bounded by a hyper-rectangular input set, i.e. \(\underline{\tilde{x}}_{k}^{[i]}\leq\tilde{x}_{k}^{[i]}\leq\overline{\tilde{x} }_{k}^{[i]}\), as in (10), then (12) holds._

### _ReLU activation functions_

First, we define \(\hat{z}_{i,k}^{\ell+1}=W_{i}^{\ell}z_{i,k}^{\ell}+b_{i}^{\ell}\) for \(\ell=0,\ldots,L-1\), \(\hat{z}_{i,k}^{L+1}=W_{i}^{L}z_{i,k}^{L}+b_{i}^{L}-\underline{u}_{i}\), and \(\hat{z}_{i,k}^{L+2}=-z_{i,k}^{L+1}+\overline{u}_{i}-\underline{u}_{i}\). Then, we can write \(z_{i,k}^{\ell}=\sigma_{i}^{\ell-1}(\hat{z}_{i,k}^{\ell})\) for \(\ell=1,\ldots,L+2\), and write the concatenation of activation functions as

\[\boldsymbol{\nu}_{i,k}=\sigma_{i}\left(\hat{\boldsymbol{\nu}}_{i,k}\right), \tag{13}\]

where \(\boldsymbol{\nu}_{i,k}^{\top}=\begin{bmatrix}z_{i,k}^{1\top}&\cdots&z_{i,k}^{L +2^{\top}}\end{bmatrix}\), \(\hat{\boldsymbol{\nu}}_{i,k}^{\top}=\begin{bmatrix}\hat{z}_{i,k}^{1\top}& \cdots&\hat{z}_{i,k}^{L+2^{\top}}\end{bmatrix}\), and \(\sigma_{i}:\mathbb{R}^{n_{i}^{i}+2n_{u}}\rightarrow\mathbb{R}^{n_{i}^{i}+2n_{u}}\) is applied elementwise, where \(n_{i}^{i}=\sum_{\ell=1}^{L}\tilde{n}_{\ell}^{i}\). Using [3, Definition 2] and [3, Lemma 3], we can relax (13) as

\[\begin{bmatrix}\hat{\boldsymbol{\nu}}_{i,k}\\ \boldsymbol{\nu}_{i,k}\\ 1\end{bmatrix}^{\top}Q^{i}(\lambda,\nu,\eta)\begin{bmatrix}\hat{\boldsymbol{ \nu}}_{i,k}\\ \boldsymbol{\nu}_{i,k}\\ 1\end{bmatrix}\geq 0, \tag{14}\]

where \(Q^{i}(\lambda,\nu,\eta)\in\mathbb{S}^{2(n_{u}^{i}+2n_{u})+1}\) is defined in [3, Lemma 3] as \(Q\). By using a change-of-base matrix [8]

\[E_{\mathrm{mid}}^{i}=\begin{bmatrix}W_{i}^{0}&0&\cdots&0&0&0\\ 0&W_{i}^{1}&\cdots&0&0&0\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots&b_{i}^{L-1}\\ 0&0&\cdots&W_{i}^{L}&0&0&b_{i}^{L}-\underline{u}_{i}\\ 0&0&\cdots&0&-I_{n_{u}}&0&\overline{\underline{u}}_{i}-\underline{u}_{i}\\ \hline 0&I_{\tilde{n}_{1}^{i}}&\cdots&0&0&0\\ \vdots&\vdots&\ddots&\vdots&\vdots&0\\ 0&0&\cdots&I_{\tilde{n}_{1}^{i}}&0&0&\boldsymbol{0}\\ 0&0&\cdots&0&I_{n_{u}}&0\\ 0&0&\cdots&0&0&I_{n_{u}}&0\\ \hline\end{bmatrix},\]

we can write (14) as

\[\begin{bmatrix}\underline{\mathbf{z}}_{i,k}\\ 1\end{bmatrix}^{\top}\Theta^{i}(\lambda,\nu,\eta)\begin{bmatrix}\underline{ \mathbf{z}}_{i,k}\\ 1\end{bmatrix}\geq 0, \tag{15}\]

where \(\Theta^{i}(\lambda,\nu,\eta)={E_{\mathrm{mid}}^{i}}^{\top}Q^{i}(\lambda,\nu, \eta)E_{\mathrm{mid}}^{i}\). This is summarised in the following lemma.

**Lemma 3**: _Given that the activation values for the extended NN satisfy \(z_{i,k}^{\ell}=\sigma_{i}^{\ell-1}(\hat{z}_{i,k}^{\ell})\) for \(\ell=1,\ldots,L+2\), as in (13), then (15) holds._
### _Reachability algorithm_

We parameterise the overapproximation \(\widehat{\mathcal{X}}_{k+1}^{[i]}\) of the reachable set \(\mathcal{X}_{k+1}^{[i]}\) using a polytope, i.e. the intersection of \(m\) halfspaces

\[\widehat{\mathcal{X}}_{k+1}^{[i]}=\left\{x\in\mathbb{R}^{n_{x}}\ |\ H_{1}^{\top}x\leq h_{1},\ldots,H_{m}^{\top}x\leq h_{m}\right\}, \tag{16}\]

where \(H_{1},\ldots,H_{m}\in\mathbb{R}^{n_{x}}\) and \(h_{1},\ldots,h_{m}\in\mathbb{R}\), which can be written as [8]

\[\widehat{\mathcal{X}}_{k+1}^{[i]}=\bigcap_{p=1}^{m}\left\{x\in\mathbb{R}^{n_{ x}}\ \left[\begin{array}{c}x\\ 1\end{array}\right]^{\top}S_{p}(h_{p})\begin{bmatrix}x\\ 1\end{bmatrix}\leq 0\right\},\]

where

\[S_{p}(h_{p})=\begin{bmatrix}0&H_{p}\\ H_{p}^{\top}&-2h_{p}\end{bmatrix},\]

so each halfspace can be equivalently written as

\[\begin{bmatrix}x_{k+1}^{[i]}\\ 1\end{bmatrix}^{\top}S_{p}(h_{p})\begin{bmatrix}x_{k+1}^{[i]}\\ 1\end{bmatrix}\leq 0. \tag{17}\]

By using a change-of-basis matrix

\[E_{\mathrm{out},k}^{i}=\begin{bmatrix}\tilde{A}_{i}&0&\cdots&0&-B_{i}&B_{i} \overline{u}_{i}+w_{k}^{[i]}\\ 0&0&\cdots&0&0&1\end{bmatrix},\]

we can write (17) as

\[\begin{bmatrix}\mathbf{z}_{i,k}\\ 1\end{bmatrix}^{\top}\Psi_{k}^{i}(h_{p})\begin{bmatrix}\mathbf{z}_{i,k}\\ 1\end{bmatrix}\leq 0, \tag{18}\]

where \(\Psi_{k}^{i}(h_{p})=E_{\mathrm{out},k}^{i}{}^{\top}S_{p}(h_{p})E_{\mathrm{out },k}^{i}\). Note that the structure of the change-of-basis matrix differs to that in [8], as \(\tilde{A}_{i}\) is not square, to account for the multi-agent dynamics. The result is summarised in the following lemma.

**Lemma 4**: _If (18) holds for \(p=1,\ldots,m\), then \(\widehat{\mathcal{X}}_{k+1}^{[i]}\) is a polytope defined by \(H_{1},\ldots,H_{m}\) and \(h_{1},\ldots,h_{m}\), as in (16)._

### _Reachability algorithm_

By combining the results in Lemmas 1-4, we arrive at the following result for the overapproximation of the forward reachable set of agent \(i\).

**Theorem 1**: _Consider a discrete-time LTI system of the form in (1) and (2), where the structure of each MLP is given by (3a)-(3c). At time \(k\), let the state of agent \(i\) and the states of its neighbours be bounded by a hyper-rectangular input set, as in (10). If there exists a solution to_

\[\min_{\Gamma,\lambda,\nu,\eta,h_{p}} \quad h_{p}\] subject to \[\quad\Delta_{k}^{i}(\Gamma)+\Theta^{i}(\lambda,\nu,\eta)+\Psi_{k}^{i}(h_{p}) \preceq 0,\]

_for \(p=1,\ldots,m\), where \(H_{1},\ldots,H_{m}\) are specified by the user, then the resulting polytope is the solution to (4a)-(4c)._

_Proof:_ Note that we are performing reachability analysis on the system given in (8). However, from Lemma 1, the original dynamics of the form given in (1) and (2) are equivalent to those given in (8), so performing reachability analysis on the system given in (8) is identical to performing reachability analysis on the system given in (1) and (2). The remainder of the proof is similar to that in [8, Theorem 1]. First, we pre-multiply each of the three terms in the LMI by \(\left[\mathbf{z}_{i,k}{}^{\top}\quad 1\right]\) and post-multiply each term by \(\left[\mathbf{z}_{i,k}{}^{\top}\quad 1\right]^{\top}\), resulting the scalar inequality

\[\begin{bmatrix}\mathbf{z}_{i,k}\\ 1\end{bmatrix}^{\top}\Delta_{k}^{i}(\Gamma)\begin{bmatrix}\mathbf{z}_{i,k}\\ 1\end{bmatrix}+\begin{bmatrix}\mathbf{z}_{i,k}\\ 1\end{bmatrix}^{\top}\Theta^{i}(\lambda,\nu,\eta)\begin{bmatrix}\mathbf{z}_{i, k}\\ 1\end{bmatrix}\\ +\begin{bmatrix}\mathbf{z}_{i,k}\\ 1\end{bmatrix}^{\top}\Psi_{k}^{i}(h_{p})\begin{bmatrix}\mathbf{z}_{i,k}\\ 1\end{bmatrix}\leq 0.\]

From Lemmas 2 and 3 and our definitions of the input set and MLP, we know that the first two terms are non-negative, and as the overall expression is non-positive, then the last term is non-positive, i.e.

\[\begin{bmatrix}\mathbf{z}_{i,k}\\ 1\end{bmatrix}^{\top}\Psi_{k}^{i}(h_{p})\begin{bmatrix}\mathbf{z}_{i,k}\\ 1\end{bmatrix}\leq 0.\]

Hence, from Lemma 4, we can conclude that \(\widehat{\mathcal{X}}_{k+1}^{[i]}\) is a polytope defined by \(H_{1},\ldots,H_{m}\) and \(h_{1},\ldots,h_{m}\), where \(\widehat{\mathcal{X}}_{k+1}^{[i]}\supseteq\mathcal{X}_{k+1}^{[i]}\).

This result is implemented in Algorithm 1, which presents a method for computing hyper-rectangular output sets (a special case of the polytopic sets), given hyper-rectangular input sets. \(\mathrm{Reach}(\widehat{\mathcal{X}}_{k}^{[i]})\) represents the solution to the semidefinite program in Theorem 1, given \(\widehat{\mathcal{X}}_{k}^{[i]}\), and \(\delta_{\cdot,\cdot}\), is the Kronecker delta. This algorithm can be used iteratively to find approximate reachable sets for \(k+2\), \(k+3\), etc. Note also that this approach is parallelisable across agents.

```
1:input sets \(\mathcal{X}_{k}^{[1]},\ldots,\mathcal{X}_{k}^{[M]}\)
2:for\(i=1,\ldots,M\)do
3:\(\widehat{\mathcal{X}}_{k}^{[i]}\leftarrow\mathcal{X}_{k}^{[i]}\times\mathcal{X}_ {k}^{[g(i,1)]}\times\ldots\times\mathcal{X}_{k}^{[g(i,q_{i})]}\)
4:for\(p=1,\ldots,2n_{x}\)do
5:if\(p\leq n_{x}\)then
6:\(H_{p}\leftarrow\left[\delta_{1,p}\quad\cdots\quad\delta_{n_{x},p}\right]^{\top}\)
7:else
8:\(H_{p}\leftarrow-\left[\delta_{n_{x}+1,p}\quad\cdots\quad\delta_{2n_{x},p} \right]^{\top}\)
9:endif
10:\(h_{p,k+1}^{[i]}\leftarrow\mathrm{Reach}\left(\widehat{\mathcal{X}}_{k}^{[i]}\right)\)
11:endfor
12:\(\widehat{\mathcal{X}}_{k+1}^{[i]}\leftarrow\left[h_{1,k+1}^{[i]}\quad\cdots \quad h_{n_{x},k+1}^{[i]}\right]^{\top}\)
13:\(\hat{\mathcal{X}}_{k+1}^{[i]}\leftarrow-\left[h_{n_{x}+1,k+1}^{[i]}\quad\cdots \quad h_{2n_{x},k+1}^{[i]}\right]^{\top}\)
14:\(\widehat{\mathcal{X}}_{k+1}^{[i]}\leftarrow\left\{x\in\mathbb{R}^{n_{x}}\ |\ \hat{\mathcal{X}}_{k+1}^{[i]}\leq x\leq\overline{\mathcal{X}}_{k+1}^{[i]}\right\}\)
15:endfor
16:approximate reachable sets \(\widehat{\mathcal{X}}_{k+1}^{[1]},\ldots,\widehat{\mathcal{X}}_{k+1}^{[M]}\)
```

**Algorithm 1** One-step forward reachability analysis with hyper-rectangular constraints

## V Model Uncertainty

In this section, we introduce model uncertainty into the dynamics and extend Theorem 1 to account for this. Consider the dynamics in (1). Instead of assuming direct knowledge of the state and input matrices, we now assume that they lie
in the convex hull of a finite number of matrices, i.e.

\[A_{ii}\in\mathrm{co}\left\{A_{ii}^{1},\ldots,A_{ii}^{C_{i}}\right\},\quad B_{i} \in\mathrm{co}\left\{B_{i}^{1},\ldots,B_{i}^{D_{i}}\right\}, \tag{19}\]

where \(C_{i},D_{i}\in\mathbb{Z}^{+}\), for \(i=1,\ldots,M\). Note that \(\Psi_{k}^{i}\) in (18) depends on \(A_{ii}\) and \(B_{i}\), so we can write \(\Psi_{k}^{i}(h_{p};A_{ii},B_{i})\).

**Theorem 2**: _Consider the formulation in Theorem 1, but with the addition of model uncertainty in the state and input matrices described by (19). If there exists a solution to_

\[\min_{\Gamma,\lambda,\nu,\eta,h_{p}} h_{p}\] \[\text{subject to} \Delta_{k}^{i}(\Gamma)+\Theta^{i}(\lambda,\nu,\eta)+\Psi_{k}^{i}(h _{p};A_{ii}^{c},B_{i}^{d})\preceq 0,\] \[\forall\ c\in\left\{1,\ldots,C_{i}\right\},\ d\in\left\{1,\ldots, D_{i}\right\},\]

_for \(p=1,\ldots,m\), where \(H_{1},\ldots,H_{m}\) are specified by the user, then the resulting polytope is the solution to (4a)-(4c)._

_Proof:_ The proof relies on the fact that \(\Psi_{k}^{i}(h_{p};A_{ii}^{c},B_{i}^{d})\) depends affinely on \(A_{ii}^{c}\) and \(B_{i}^{d}\). The proof has some similarities to that in [9], but the forms of the matrices and dynamics differ, so we include the proof for completeness. First, note that we can write (19) as

\[A_{ii}=\sum_{c=1}^{C_{i}}\alpha_{c}A_{ii}^{c},\quad B_{i}=\sum_{d=1}^{D_{i}} \beta_{d}B_{i}^{d}, \tag{20}\]

where \(\sum_{c=1}^{C_{i}}\alpha_{c}=1,\ \alpha_{c}\geq 0\ \forall c\in\left\{1, \ldots,C_{i}\right\}\) and \(\sum_{d=1}^{D_{i}}\beta_{d}=1,\ \beta_{d}\geq 0\ \forall d\in\left\{1,\ldots,D_{i}\right\}\). From the definitions of \(\Psi_{k}^{i}\) and \(\bar{A}_{i}\), we can write

\[\Psi_{k}^{i}(h_{p};A_{ii}^{c},B_{i}^{d})=\begin{bmatrix}0&\Phi_{k}^{i}(A_{ii} ^{c},B_{i}^{d})\\ \Phi_{k}^{i}(A_{ii}^{c},B_{i}^{d})^{\top}&\Xi_{k}^{i}(h_{p};B_{i}^{d})\end{bmatrix},\]

where

\[\Phi_{k}^{i}(A_{ii}^{c},B_{i}^{d})=\begin{bmatrix}\left[A_{ii}^{c}&A_{ig(i,1)} &\cdots&A_{ig(i,q_{i})}\right]^{\top}H_{p}\\ &0&\left[-B_{i}^{d}\right]^{\top}H_{p}\end{bmatrix},\]

and

\[\Xi_{k}^{i}(h_{p};B_{i}^{d})=2H_{p}^{\top}\left(B_{i}^{d}\overline{u}_{i}+w_{ k}^{[i]}\right)-2h_{p},\]

so it can be seen that \(\Psi_{k}^{i}(h_{p};A_{ii}^{c},B_{i}^{d})\) depends affinely on \(A_{ii}^{c}\) and \(B_{i}^{d}\). Hence, multiplying by \(\alpha_{c}\) and \(\beta_{d}\), summing over \(c=1,\ldots,C_{i}\) and \(d=1,\ldots,D_{i}\), and using \(\sum_{c=1}^{C_{i}}\alpha_{c}A_{ii}^{c}=A_{ii}\), \(\sum_{c=1}^{C_{i}}\alpha_{c}=1\), \(\sum_{d=1}^{D_{i}}\beta_{d}B_{i}^{d}=B_{i}\) and \(\sum_{d=1}^{D_{i}}\beta_{d}=1\) from (20) gives

\[\sum_{c=1}^{C_{i}}\sum_{d=1}^{D_{i}}\alpha_{c}\beta_{d}\Psi_{k}^{i}(h_{p};A_{ii }^{c},B_{i}^{d})=\Psi_{k}^{i}(h_{p};A_{ii},B_{i}). \tag{21}\]

Finally, multiplying the LMIs in the semidefinite program by \(\alpha_{c}\) and \(\beta_{d}\) and summing over \(c=1,\ldots,C_{i}\) and \(d=1,\ldots,D_{i}\), such that

\[\sum_{c=1}^{C_{i}}\sum_{d=1}^{D_{i}}\alpha_{c}\beta_{d}\left[\Delta _{k}^{i}(\Gamma)+\Theta^{i}(\lambda,\nu,\eta)\right]\\ +\sum_{c=1}^{C_{i}}\sum_{d=1}^{D_{i}}\alpha_{c}\beta_{d}\Psi_{k}^ {i}(h_{p};A_{ii}^{c},B_{i}^{d})\preceq 0,\]

results in

\[\Delta_{k}^{i}(\Gamma)+\Theta^{i}(\lambda,\nu,\eta)+\Psi_{k}^{i}(h_{p};A_{ii},B_ {i})\preceq 0,\]

which follows from the fact that \(\sum_{c=1}^{C_{i}}\sum_{d=1}^{D_{i}}\alpha_{c}\beta_{d}=1\) and from (21). \(\blacksquare\)

This result is useful, as we only have to solve a finite number of semidefinite programs to find an overapproximation of the forward reachable set for all convex combinations of \(A_{ii}^{1},\ldots,A_{ii}^{C_{i}}\) and \(B_{i}^{1},\ldots,B_{i}^{D_{i}}\).

## VI Experiments

In this section, we use two realistic examples of multi-agent systems to demonstrate our results. We also give a comparison to the approach proposed at the start of Section III, in which the states of each agent are augmented into one overall state and existing reachability methods are applied on the overall system. We also introduce model uncertainty into one of the systems and analyse this case. Simulations were performed on MATLAB, and CVX with MOSEK was used to solve the semidefinite programs. The NNs were trained to approximate distributed MPC schemes with a given horizon; the systems were simulated with MPC, and the resulting input-output data pairs were used to train the NNs.

### _Vehicle platooning_

In the first example, we consider the example of control of a platoon of vehicles. There are several forms of this problem, and control of a vehicular platoon has a number of benefits, including improved safety, higher road capacity, lower emissions, and/or reduced congestion [25, 26, 27, 28]. In this example, we consider the adaptive cruise control (ACC) problem, in which each vehicle aims to maintain a fixed distance from the vehicle in front, whilst travelling at a given velocity. The continuous-time longitudinal dynamics of each vehicle are given by [29]

\[\dot{x}^{[i]}(t)=\bar{A}_{ii}x^{[i]}(t)+\bar{A}_{ii-1}x^{[i-1]}(t)+\bar{B}_{i}u^ {[i]}(t),\]

(note that \(\mathcal{N}_{i}=\{i-1\}\)) where the state vector is

\[x^{[i]}(t)=\begin{bmatrix}e^{[i]}(t)&v^{[i]}(t)&a^{[i]}(t)\end{bmatrix}^{\top},\]

where \(e^{[i]}(t)\) is the distance error between vehicle \(i\) and vehicle \(i-1\) (i.e. if the desired distance between vehicle \(i\) and the vehicle in front, \(i-1\), is \(\bar{d}^{[i]}\) and the actual distance is \(d^{[i]}(t)\), then \(e^{[i]}(t)=d^{[i]}(t)-\bar{d}_{i}\)), \(v^{[i]}(t)\) is the velocity of the \(i^{\text{th}}\) vehicle, \(a^{[i]}(t)\) is the acceleration of the \(i^{\text{th}}\) vehicle, and

\[\bar{A}_{ii} =\begin{bmatrix}0&-1&0\\ 0&0&1\\ 0&0&-\frac{1}{\tau}\end{bmatrix},\quad\bar{A}_{ii-1}=\begin{bmatrix}0&1&0\\ 0&0&0\\ 0&0&0\end{bmatrix},\] \[\bar{B}_{i} =\begin{bmatrix}0&0&\frac{1}{\tau}\end{bmatrix}^{\top},\]

where \(\tau\) is the engine time constant [30] and \(u^{[i]}(t)\) is the \(i^{\text{th}}\) acceleration input. Note that the lead vehicle (\(i=1\)) has no physical neighbour, but this can be resolved by imagining a virtual vehicle [29] with state \(x^{[0]}=\begin{bmatrix}0&\bar{v}&0\end{bmatrix}^{\top}\), where \(\bar{v}\) is the reference velocity (to be maintained by the platoon). This
is shown in Figure 2. The dynamics are discretised assuming zero-order hold (ZOH) with a sample period \(T=0.1\) s, treating \(x^{[i-1]}(t)\) and \(u^{[i]}(t)\) as exogenous inputs, such that the discrete-time dynamics are in the form in (1), where \(w_{k}^{[i]}=0\)\(\forall i\in\mathcal{I}\). The NNs have \(2\) hidden layers, both with \(15\) neurons. The \(i^{\text{th}}\) control input is

\[u_{k}^{[i]}=\operatorname{sat}_{\mathcal{U}_{i}}\left[\pi_{ii-1}\left(\begin{bmatrix} e_{k}^{[i]}\\ v_{k}^{[i-1]}-v_{k}^{[i]}\\ a_{k}^{[i-1]}-a_{k}^{[i]}\end{bmatrix}\right)\right],\]

where \(e_{k}^{[i]}\equiv e^{[i]}(kT)\), \(v_{k}^{[i]}\equiv v^{[i]}(kT)\) and \(a_{k}^{[i]}\equiv a^{[i]}(kT)\). The forward reachable sets were computed for five time steps, \(M=9\) agents, initial conditions given by \(\mathcal{X}_{\mathbf{l}}^{[i]}=\{x\in\mathbb{R}^{3}\mid\underline{x}\leq x \leq\underline{\overline{x}}\}\)\(\forall i\in\mathcal{I}\), where \(\underline{x}^{[i]}=\begin{bmatrix}-0.1&19.95&-0.01\end{bmatrix}\) and \(\overline{x}^{\top}=\begin{bmatrix}0.1&20.05&0.01\end{bmatrix}\), and controller limits given by \(\overline{u}_{i}=-\underline{u}_{i}=5\)\(\forall i\in\mathcal{I}\). A step change of \(-2\) was applied to the reference velocity at \(k=0\) (from \(\bar{v}=20\) to \(\bar{v}=18\)). The results are shown in Figure 3 for the first three vehicles.

Note that the step change in \(\bar{v}\) takes some time to propagate down the platoon, hence the difference in range between the plots. A comparison between the computation time for this approach (Reach-SDP-MA) and the existing method (Reach-SDP) is shown in Table 1 for different values of \(M\).

We then extend the vehicle platooning example to account for model uncertainty. Consider the case in which \(A_{ii}\in\operatorname{co}\left\{(1-\delta)A^{0},(1+\delta)A^{0}\right\}\) and \(B_{i}\in\operatorname{co}\left\{(1-\delta)B^{0},(1+\delta)B^{0}\right\}\), where \(A^{0}\in\mathbb{R}^{3\times 3}\) and \(B^{0}\in\mathbb{R}^{3}\) are the nominal values of \(A_{ii}\) and \(B_{i}\), respectively, and \(\delta=0.01\). The results are shown in Figure 4 for the first three vehicles.

Because of the uncertainty, the size of the exact reachable sets increases, so the overapproximations are larger (compared to Figure 3).

### _Power network system_

For the second example, we consider automatic generation control (AGC) of a power network system. Unlike the previous example, the dynamics are not identical across agents, and some agents have more than one neighbour. There is also an exogenous input term. This system consists of \(M\) generation areas, and the aim is to reduce the frequency deviation in each area, in spite of load changes. Common approaches to this control problem include decentralised and distributed MPC schemes [20, 31, 32]. The continuous-time dynamics of each area are given by [31, 33]

\[\dot{x}^{[i]}(t)=\bar{A}_{ii}x^{[i]}(t)+\sum_{j\in\mathcal{N}_{i}}\bar{A}_{ij }x^{[j]}(t)+\bar{B}_{i}u^{[i]}(t)+\bar{L}_{i}\Delta P_{L}^{[i]}(t),\]

where the state vector for area \(i\) is

\[x^{[i]}(t)=\left[\Delta\theta^{[i]}(t)\quad\Delta\omega^{[i]}(t)\quad\Delta P _{m}^{[i]}(t)\quad\Delta P_{v}^{[i]}(t)\right]^{\top},\]

where \(\Delta\theta^{[i]}(t)\), \(\Delta\omega^{[i]}(t)\), \(\Delta P_{m}^{[i]}(t)\) and \(\Delta P_{v}^{[i]}(t)\) are the deviations in rotor angle, frequency, mechanical power and steam valve position, respectively, from the nominal values [34], \(u^{[i]}(t)\) is the reference power, \(\Delta P_{L}^{[i]}(t)\) is the local power load, and

\[\bar{A}_{ii}=\begin{bmatrix}0&1&0&0\\ -\frac{\sum_{j\in\mathcal{N}_{i}}P_{ij}}{2H_{i}}&-\frac{D_{i}}{2H_{i}}&\frac{1 }{2H_{i}}&0\\ 0&0&-\frac{\bar{D}_{i}}{T_{i}}&\frac{1}{T_{i}}\\ 0&-\frac{1}{R_{i}T_{g_{i}}}&0&-\frac{1}{T_{g_{i}}}\end{bmatrix},\]

\[\bar{A}_{ij}=\begin{bmatrix}0&0&0&0\\ \frac{P_{ij}}{2H_{i}}&0&0&0\\ 0&0&0&0\end{bmatrix},\quad\bar{B}_{i}=\begin{bmatrix}0&0&0&T_{g_{i}}\end{bmatrix}^ {\top},\]

\[\bar{A}_{ij}=\begin{bmatrix}0&0&0&0\\ \frac{P_{ij}}{2H_{i}}&0&0&0\\ 0&0&0&0\end{bmatrix},\quad\bar{L}_{i}=\begin{bmatrix}0&-\frac{1}{2H_{i}}&0&0 \end{bmatrix}^{\top},\]

\begin{table}
\begin{tabular}{|c||c|c|c|c|c|} \hline \(M\) & \(1\) & \(2\) & \(3\) & \(4\) & \(5\) \\ \hline Reach-SDP-MA & \(3.85\) & \(7.04\) & \(11.60\) & \(14.17\) & \(20.82\) \\ \hline Reach-SDP [8] & \(4.01\) & \(75.39\) & \(562.94\) & \(3636.62\) & \(41025.05\) \\ \hline \end{tabular}
\end{table}
Table 1: Comparison of methods (times in s)

Figure 4: Plots of the reachable sets in red (solid) for distance error and velocity, and simulated trajectories (blue) for the vehicle platooning example with model uncertainty; the initial set is shown in red (dashed) – only agents \(1\) to \(3\) are shown

Figure 3: Plots of the reachable sets in red (solid) for distance error and velocity, and simulated trajectories (blue) for the vehicle platooning example; the initial set is shown in red (dashed) – only agents \(1\) to \(3\) are shown

Figure 2: Platoon of \(M\) vehicles, where each vehicle \(i\) only receives information from vehicle \(i-1\); the ‘virtual’ vehicle is shown with increased opacity
where \(P_{ij}\), \(H_{i}\), \(D_{i}\), \(T_{t_{i}}\), \(R_{i}\) and \(T_{g_{i}}\) are defined in [31]. In this example, we consider \(M=4\) generation areas (Scenario 1 in [34]), where \(\mathcal{N}_{1}=\{2\}\), \(\mathcal{N}_{2}=\{1,3\}\), \(\mathcal{N}_{3}=\{2,4\}\) and \(\mathcal{N}_{4}=\{3\}\), as shown in Figure 5.

The dynamics are discretised assuming ZOH with a sample period \(T=1\) s, treating \(x^{[j]}(t)\)\(\forall j\in\mathcal{N}_{i}\), \(u^{[i]}(t)\) and \(\Delta P_{L}^{[i]}(t)\) as exogenous inputs, such that the discrete-time dynamics are in the form in (1). The NNs have \(2\) hidden layers, both with \(10\) neurons. The \(i^{\text{th}}\) control input is

\[u_{k}^{[i]}=\operatorname{sat}_{\mathcal{U}_{i}}\left[\sum_{j\in\mathcal{N}_{ i}}\pi_{ij}\left(\begin{bmatrix}x_{k}^{[i]}\\ x_{k}^{[j]}\end{bmatrix}-\begin{bmatrix}x_{\text{ref},k}^{[i]}\\ x_{\text{ref},k}^{[j]}\end{bmatrix}\right)\right],\]

where \(x_{\text{ref},k}^{[i]}\) and \(x_{\text{ref},k}^{[j]}\) are the state and neighbour reference values, respectively. These are incorporated into the reachability analysis by modifying the first weight matrices and bias vectors. The forward reachable sets were computed for three time steps with initial conditions given by \(\mathcal{X}_{0}^{[i]}=\left\{x\in\mathbb{R}^{4}\mid\underline{x}\leq x\leq \overline{x}\right\}\)\(\forall i\in\mathcal{I}\), where \(\overline{x}^{\top}=-\underline{x}^{\top}=\left[10^{-4}\quad 10^{-7}\quad 10^{-3} \right]\), and controller limits given in [34]. A step change of \(-0.15\) was applied to \(\Delta P_{L,k}^{[i]}\equiv\Delta P_{L}^{[i]}(kT)\)\(\forall i\in\mathcal{I}\), and the results are shown in Figure 6.

## VII Conclusion and Future Work

In this paper, we presented a scalable method to overapproximate the forward reachable sets of multi-agent systems with distributed NNC architectures. After simplifying the dynamics, we presented a method to split the overall reachability problem into multiple smaller reachability problems. We then extended this approach to account for model uncertainty. The effectiveness of this method was demonstrated on realistic examples, and it was shown to be significantly faster than the overall reachability method whilst producing the same bounds. It should also be noted that the method presented in this paper can be applied to any system that can be decomposed into the form in (1); it is not necessarily specific to multi-agent systems.

Opportunities for future work include using the general framework presented in this paper to improve the efficiency of other reachability methods, such as LP-based methods. Also, further consideration could be given to synthesis of the NNCs, and reachability methods could be incorporated into the training process. Other types of activation functions and sets could also be considered.
