# On the Design of Integral Multiplex Control Protocols for Nonlinear Network
Systems with Delays

Shihao Xie  Giovanni Russo  School of Electrical and Electronic Engineering,
University College Dublin, Ireland (e-mail<EMAIL_ADDRESS>Department of Information and Electrical Engineering and Applied Mathematics,
University of Salerno, Italy (e-mail<EMAIL_ADDRESS>
###### Abstract

We consider the problem of designing control protocols for possibly nonlinear
networks with delays that not only allow the fulfillment of some desired
behaviour, but also simultaneously guarantee the rejection of polynomial
disturbances and the non-amplification of other classes of disturbances across
the network. To address this problem, we propose the systematic use of
multiplex architectures to deliver integral control protocols ensuring the
desired disturbance rejection and non-amplification properties. We then
present a set of sufficient conditions to assess these properties and hence to
design the multiplex architecture for both leaderless and leader-follower
networks with time-varying references consisting of possibly heterogeneous
nonlinearly coupled agents affected by communication delays. The effectiveness
of our conditions, which are also turned into an optimisation problem allowing
protocol design, is illustrated via both in-silico and experimental
validations with a real hardware set-up.

††thanks: An early version of this paper, without proofs and considering a
special case only, presented at NecSys $22$.

,

##  1 Introduction

Driven by the introduction of low-cost, high performance and connected
devices, network systems have considerably increased their size and
complexity. In this context, a key challenge is the design of control
protocols for the network that do not only guarantee the fulfillment of some
desired behaviour, but also satisfy the following requirements: (i) ensure
that certain classes of disturbances are rejected as they are compensated by
the protocol; (ii) guarantee that the disturbances that are not rejected are
not amplified across the network. Vehicle platoons (Stüdli et al., 2017; Feng
et al., 2019; Monteil et al., 2019), robot formations (Liu et al., 2018; Li et
al., 2019), power grids (Dörfler and Bullo, 2014; Silva et al., 2022), neural
(Jafarpour et al., 2021; Revay and Manchester, 2020; Xie et al., 2021) and
biochemical (Qian and Del Vecchio, 2018) networks are just few examples of
network systems for which these two requirements cannot be neglected in the
control design. Motivated by this, in the present paper we study the problem
of designing control protocols for possibly nonlinear networks affected by
delays that not only allow the network to achieve some desired behaviour, but
also, by means of integral actions, allow rejection of polynomial disturbances
and ensure the non-amplification of other classes of disturbances across the
network. We capture the fulfillment of these requirements via a scalability
property of the network (see Section 3 for the definitions).

#### Related works.

The study of multiplex networks has emerged within the physics and control
literature. We refer to (Mucha et al., 2010; De Domenico et al., 2013; Tran et
al., 2020; Gomez et al., 2013; Burbano Lombana and di Bernardo, 2016) and
references therein for a set of results, with literature reviews, on these
networks. We now briefly survey some related works on disturbance propagation
in networks. The study of how disturbances propagate within a network system
is a central topic for the platooning of autonomous vehicles. In particular,
the key idea behind several definitions of string stability (Swaroop and
Hedrick, 1996) in the literature is that of giving upper bounds on the
deviations induced by disturbances that are uniform with respect to the
platoon size, see e.g. (Ploeg et al., 2014; Besselink and Johansson, 2017;
Monteil et al., 2019). These works assume delay-free inter-vehicle
communications and an extension to delayed platoons with linear vehicle
dynamics can be found in e.g. (di Bernardo et al., 2015). Consensus stability
for (disturbance-free) robotic networks subject to delays is instead
considered in (Fan et al., 2021). The problem of designing a distributed
integral action (Freeman et al., 2006; Andreasson et al., 2014; Seyboth and
Allgöwer, 2015) has been recently considered in the context of string
stability to design protocols able to compensate constant disturbances for
delay-free platoons (Knorn et al., 2014b; Silva et al., 2021). For networks
with delay-free interconnections and no integral actions, we also recall
results on mesh stability (Seiler et al., 1999) where linear dynamics are
considered and its extension to nonlinear networks (Pant et al., 2002).
Leader-to-formation stability is instead considered in (Tanner et al., 2004)
and it characterizes network behaviour with respect to inputs from a leader.
For delay-free leaderless networks with regular topology, the scalability
property has been recently investigated in (Besselink and Knorn, 2018), where
Lyapunov-based conditions are given; for networks with arbitrary topologies
and delays, sufficient conditions for scalability can be found in (Xie et al.,
2021), which leverages contraction theory arguments for time-delayed systems.
Contracting systems exhibit transient and asymptotic behaviors (Lohmiller and
Slotine, 1998) that are desirable when designing network systems (Centorrino
et al., 2022); we refer to (Aminzare and Sontag, 2014; di Bernardo et al.,
2016; Tsukamoto et al., 2021) and references therein for further details. We
also recall (Wang and Slotine, 2006) which shows, using the Euclidean metric,
how contraction is preserved through certain time-delayed communications and
(Shiromoto et al., 2019) where conditions for the synthesis of distributed
controls are given by using separable metric structures. For delay-free
systems, based on the use of contraction, a sufficient condition for the
stability of a feedback loop consisting of an exponentially stable multi-input
multi-output nonlinear plant and an integral controller (to compensate
constant disturbances) has been obtained in (Simpson-Porco, 2021). The problem
of constant output regulation for a class of input-affine multi-input multi-
output nonlinear systems with constant disturbances, has also been recently
tackled via contraction in (Giaccagli et al., 2021). Other works have also
shown that contraction using non-Euclidean metrics can be useful to study a
wide range of biological (Russo et al., 2010a), neural (Davydov et al., 2021;
Ofir et al., 2022; Centorrino et al., 2022) and engineered (Monteil and Russo,
2017) networks. We also recall (Russo and Wirth, 2022) the recent extension of
contraction to dynamical systems on time scales, i.e. systems evolving on
arbitrary (potentially non-uniform) time domains.

#### Statement of Contributions.

We present an integral control design for nonlinear networks with delays that
not only allows the network to achieve some desired behaviour, but also: (i)
guarantees rejection of polynomial disturbances; and (ii) ensures the non-
amplification of other classes of disturbances. To the best of our knowledge,
these are the first results for the design of control protocols able to
simultaneously handle these requirements and, in the context of the above
literature, our main technical contributions are summarized as follows:

  * (i)

we propose a novel, systematic, use of multiplex architectures to deliver
protocols ensuring a scalability property for the network. This, besides
guaranteeing the fulfillment of a desired network behaviour, implies rejection
of the polynomial components of the disturbances and the non-amplification of
their, e.g. residual, non-polynomial components;

  * (ii)

our main technical results give sufficient conditions to assess these
properties and hence to design the multiplex architecture – we are not aware
of any other result that allows to fulfill these properties for networks
consisting of nonlinear (possibly, heterogeneous) agents affected by
communication delays. Moreover, the results allow to consider leaderless and
leader-follower networks, arbitrary topologies, nonlinear protocols and time-
varying references;

  * (iii)

we show that our conditions can be effectively used to design control
protocols. After discussing certain design implications of the results, we
show that the conditions can be recast into an optimisation problem that
allows to design the control protocol;

  * (iv)

finally, we validate our protocols both in-silico and via a real hardware
test-bed. All the experiments confirm the effectiveness of our approach111The
code, data and all the required information to replicate all our results are
available at https://tinyurl.com/4wyacf7z together with video recordings of
the experiments..

An early version of the results was presented in (Xie and Russo, 2022) where
only first-order polynomials rejection was considered and no proofs were
given.

##  2 Mathematical preliminaries

Let $A$ be a $m\times m$ real matrix, we denote by $\|A\|_{p}$ the matrix norm
induced by the $p$-vector norm $|\cdot|_{p}$. The matrix measure of $A$
induced by $|\cdot|_{p}$ is defined as $\mu_{p}(A):=\lim_{h\rightarrow
0^{+}}\frac{\|I+hA\|_{p}-1}{h}$. We write $A\succeq 0$ when $A$ is positive
semi-definite and $A\preceq 0$ when it is negative semi-definite. The
symmetric part of $A$ is denoted as $[A]_{s}:=\frac{A+A^{\mathsf{T}}}{2}$.
Given a piece-wise continuous signal $w_{i}(t)$, we let
$\|w_{i}(\cdot)\|_{\mathcal{L}_{\infty}^{p}}:=\sup_{t}|w_{i}(t)|_{p}$. We
denote by $I_{n}$ the $n\times n$ identity matrix and by $0_{m\times n}$ the
$m\times n$ zero matrix (if $m=n$ we simply write $0_{n}$). We let
$\textbf{diag}\\{a_{1},\ldots,a_{N}\\}$ be a diagonal matrix with diagonal
elements $a_{i}$. For a generic set $\mathcal{A}$, its cardinality is denoted
by $|\mathcal{A}|$. Let $f$ be a smooth function, then we denote by $f^{(n)}$
the $n$-th derivative of $f$. We recall that a continuous function
$\alpha:[0,a)\rightarrow[0,\infty)$ is said to belong to class $\mathcal{K}$
if it is strictly increasing and $\alpha(0)=0$. It is said to belong to class
$\mathcal{K}_{\infty}$ if $a=\infty$ and $\alpha(r)\rightarrow\infty$ as
$r\rightarrow\infty$. A continuous function
$\beta:[0,a)\times[0,\infty)\rightarrow[0,\infty)$ is said to belong to class
$\mathcal{KL}$ if, for each fixed $s$, the mapping $\beta(r,s)$ belongs to
class $\mathcal{K}$ with respect to $r$ and, for each fixed $r$, the mapping
$\beta(r,s)$ is decreasing with respect to $s$ and $\beta(r,s)\rightarrow 0$
as $s\rightarrow\infty$.

The following lemma is stated in (Xie et al., 2021) and follows directly from
(Russo et al., 2010b). We let $|\cdot|_{S}$ and $\mu_{S}(\cdot)$ be,
respectively, any $p$-vector norm and its induced matrix measure on
$\mathbb{R}^{N}$. In particular, the norm $|\cdot|_{S}$ is monotone, i.e. for
any non-negative $N$-dimensional vector $x,y\in\mathbb{R}_{\geq 0}^{N}$,
$x\leq y$ implies that $|x|_{S}\leq|y|_{S}$ where the inequality $x\leq y$ is
component-wise.

######  Lemma 1.

Consider the vector
$\eta:=[\eta_{1}^{\mathsf{T}},\ldots,\eta_{N}^{\mathsf{T}}]^{\mathsf{T}}$,
$\eta_{i}\in\mathbb{R}^{n}$. We let
$|\eta|_{G}:=|\left[|\eta_{1}|_{G_{1}},\ldots,|\eta_{N}|_{G_{N}}\right]|_{S}$,
with $|\cdot|_{G_{i}}$ being norms on $\mathbb{R}^{n}$, and denote by
$\|\cdot\|_{G},\mu_{G}(\cdot)$ ($\|\cdot\|_{G_{i}},\mu_{G_{i}}(\cdot)$) the
matrix norm and measure induced by $|\cdot|_{G}$ ($|\cdot|_{G_{i}}$). Finally,
let:

  1. 1.

$A:=({A_{ij}})_{i,j=1}^{N}\in\mathbb{R}^{nN\times nN}$,
$A_{ij}\in\mathbb{R}^{n\times n}$;

  2. 2.

$\hat{A}:=(\hat{A}_{ij})_{i,j=1}^{N}\in\mathbb{R}^{N\times{N}}$, with
$\hat{A}_{ii}:=\mu_{G_{i}}(A_{ii})$ and $\hat{A}_{ij}:=\|A_{ij}\|_{G_{i,j}}$,
$\|A_{ij}\|_{G_{i,j}}:=\sup_{|x|_{G_{i}}=1}|A_{ij}x|_{G_{j}}$;

  3. 3.

$\bar{A}:=(\bar{A}_{ij})_{i,j=1}^{N}\in\mathbb{R}^{N\times N}$, with
$\bar{A}_{ij}:=\|A_{ij}\|_{G_{i,j}}$.

Then:

  * (i)

$\mu_{G}(A)\leq\mu_{S}(\hat{A})$;

  * (ii)

$\|A\|_{G}\leq\|{\bar{A}}\|_{S}$.

The next lemma follows from (Wen et al., 2008, Theorem $2.4$).

######  Lemma 2.

Let $u:[t_{0}-\tau_{\max},+\infty)\rightarrow\mathbb{R}_{\geq 0}$ ,
$\tau_{\max}<+\infty$ and assume that

$D^{+}u(t)\leq au(t)+b\sup_{t-\tau(t)\leq s\leq t}u(s)+c,\ \ t\geq t_{0}$

with:

  * (1)

$\tau(t)$ being bounded and non-negative, i.e.
$0{\leq}\tau(t)\leq\tau_{\max}$, $\forall t$;

  * (2)

$u(t)=|\varphi(t)|$, $\forall t\in[t_{0}-\tau_{\max},t_{0}]$ where
$\varphi(t)$ is bounded in $[t_{0}-\tau_{\max},t_{0}]$;

  * (3)

$a<0$, $b\geq 0$ and $c\geq 0$ and that there exists some $\sigma>0$ such that
$a+b\leq-\sigma<0,\forall t\geq t_{0}$.

Then:

$u(t)\leq\sup_{t_{0}-\tau_{\max}\leq s\leq
t_{0}}u(s)e^{-\hat{\lambda}(t-t_{0})}+\frac{c}{\sigma},$

where $\hat{\lambda}:=\inf_{t\geq
t_{0}}\\{\lambda|\lambda(t)+a+be^{\lambda(t)\tau(t)}=0\\}>0$ is the
convergence rate of the system.

##  3 Statement of the control problem

We now describe the proposed control architecture (Section 3.1) and formalize
the control goal (Section 3.2).

###  3.1 Architecture description

We consider a network system of $N>1$ agents with the dynamics of the $i$-th
agent given by

$\displaystyle\begin{split}\dot{x}_{i}(t)&=f_{i}(x_{i},t)+u_{i}(t)+d_{i}(t),\
\ t\geq t_{0}\geq 0,\\\ y_{i}(t)&=g_{i}(x_{i}),\end{split}$ (1)

with initial conditions $x_{i}(t_{0})$, $i=1,\dots,N$, and where: (i)
$x_{i}(t)\in\mathbb{R}^{n}$ is the state of the $i$-th agent; (ii)
$u_{i}(t)\in\mathbb{R}^{n}$ is the control input; (iii)
$d_{i}(t)\in\mathbb{R}^{n}$ is an external disturbance signal on the agent;
(iv) $f_{i}:\mathbb{R}^{n}\times\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}^{n}$
is the intrinsic dynamics of the agent, which is assumed to be smooth; (v)
$g_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{q}$ is the output function for the
$i$-th agent. We consider disturbances of the form:

$\displaystyle
d_{i}(t)=w_{i}(t)+\bar{d}_{i}(t):=w_{i}(t)+\sum_{k=0}^{m-1}\bar{d}_{i,k}\cdot
t^{k},$ (2)

where $w_{i}(t)$ is a piece-wise continuous signal and $\bar{d}_{i,k}$’s are
constant vectors. Disturbances of the form of (2) can be thought of as the
superposition of a polynomial disturbance of order $m-1$ (denoted as
$\bar{d}_{i}(t)$) and the signal $w_{i}(t)$. This latter signal might capture
residual terms in the disturbance that are not modeled with the polynomial.
For example, in the special case where $m=1$ in (2) we have
$d_{i}(t)=w_{i}(t)+\bar{d}_{i,0}$. As noted in e.g. (Silva et al., 2021), in
the context of platooning, these types of disturbances model situations when a
platoon of vehicles encounters a slope: this gives rise to a constant
disturbance on the vehicle acceleration and the term $w_{i}(t)$ can then model
small bumps along the slope. See also Remark 2 below.

Figure 1: Schematic representation of the multiplex architecture. Only one
disturbance, i.e. $d_{1}(t)$ on agent $1$, is shown and the reference signal
is omitted. The layers can have different topologies, which can be both
directed and undirected.

Our goal is to design the control protocol $u_{i}(t)$ in (1) so that the
polynomial disturbance in (2) is rejected, while ensuring non-amplification of
the residual disturbance $w_{i}(t)$ – this is captured via a scalability
property, see Section 3.2 for the rigorous definition and statement of the
control goal. To do so, we propose the multiplex network architecture
schematically shown in Figure 1. In such a figure, the network system is in
layer $0$ and the multiplex layers (i.e. layer $1,\ldots,m$) concur to build
up the control protocol. This is of the form:

$\displaystyle\begin{split}u_{i}(t)&=h_{i,0}(x(t),x_{l}(t),t)\\\
&+h_{i,0}^{(\tau)}(x(t-\tau(t)),x_{l}(t-\tau(t)),t)+r_{i,1}(t),\\\
\dot{r}_{i,1}(t)&=h_{i,1}(x(t),x_{l}(t),t)\\\
&+h_{i,1}^{(\tau)}(x(t-\tau(t)),x_{l}(t-\tau(t)),t)+r_{i,2}(t),\\\ &\ \
\vdots\\\ \dot{r}_{i,m}(t)&=h_{i,m}(x(t),x_{l}(t),t)\\\
&+h_{i,m}^{(\tau)}(x(t-\tau(t)),x_{l}(t-\tau(t)),t),\end{split}$ (3)

where $r_{i,k}(t)$ is the output generated by the multiplex layer
$k\in\\{1,\ldots,m\\}$ and where $\tau(t)\leq\tau_{\max}$, $\forall t$. As
illustrated in Figure 1, the multiplex layer $k\in\\{1,\ldots,m\\}$ receives
information from the agents (on layer $0$) and outputs a signal to the layer
immediately below, i.e. layer $k-1$. In (3): (i)
$x(t):=[x_{1}^{\mathsf{T}}(t),\ldots,x_{N}^{\mathsf{T}}(t)]^{\mathsf{T}}$ is
the state of the network; (ii)
$x_{l}(t):=[x_{l_{1}}^{\mathsf{T}}(t),\ldots,x_{l_{M}}^{\mathsf{T}}(t)]^{\mathsf{T}}$
is the reference signal, possibly provided by a group of $M$ leaders; (iii)
the functions
$h_{i,k}:\mathbb{R}^{nN}\times\mathbb{R}^{nM}\times\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}^{n}$ and
$h_{i,k}^{(\tau)}:\mathbb{R}^{nN}\times\mathbb{R}^{nM}\times\mathbb{R}_{\geq
0}\rightarrow\mathbb{R}^{n},k\in\\{0,1,\ldots,m\\}$, include both (leader and
leaderless) delayed and delay-free couplings (see Remark 3 for an example).
These functions model, on layer $k$, the possible connections between the
agents and leaders. Note that not all the agents necessarily receive
information from leaders (if any). In what follows, we simply term the smooth
coupling functions $h_{i,k}(\cdot,\cdot)$ as delay-free coupling functions,
while the functions $h_{i,k}^{(\tau)}(\cdot,\cdot)$ are termed as delayed
coupling functions. As noted in e.g. (Xie et al., 2021) situations where there
is an overlap between delayed and non-delayed communications naturally arise
in the context of e.g. platooning, robotics formation control and neural
networks. Without loss of generality, in (3) we set, $\forall
s\in[t_{0}-\tau_{\max},t_{0}],\forall i=1,\ldots,N$, $\forall k=1,\ldots,m$,
$x_{i}(s)=\varphi_{i}(s)$ and $r_{i,k}(s)=\phi_{i,k}(s)$, with
$\varphi_{i}(s)$ and $\phi_{i,k}(s)$ being continuous and bounded functions in
$[t_{0}-\tau_{\max},t_{0}]$.

######  Remark 1.

The time varying delay in the multiplexed control protocol (3) is the same for
all the agents. Network systems affected by these homogeneous delays naturally
arise in the context of multi-agent systems. For example, in (Abolfazli et
al., 2021) string stability of platoon is studied when the agents are affected
by homogeneous constant delays. Homogeneous delays are also considered in the
context of consensus (Wang et al., 2013), synchronization of networks with
linear dynamics (Gao et al., 2008) and nonlinear dynamics (Xie and Liu, 2013).

######  Remark 2.

We consider disturbances consisting of a polynomial component and a piece-wise
continuous component. Polynomial disturbances are commonly considered in the
literature. See e.g. (Kim et al., 2010), where observers for these
disturbances are devised and (Park et al., 2012), where the problem of
rejecting these disturbances is considered. The signal $w_{i}(t)$ can be
physically interpreted as a (typically, small) discrepancy between the
polynomial disturbance model and the actual disturbance signal. For
platooning, rejection of constant disturbances (i.e. the disturbance in (2)
when $m=1$) has been considered in (Knorn et al., 2014a; Silva et al., 2021).

######  Remark 3.

Protocols of the form of (3) arise in a wide range of applications. For
example, in the context of formation control, typical choices for the coupling
functions, see e.g. (Xie et al., 2021; Lawton et al., 2003; Olfati-Saber and
Murray, 2004; Tanner et al., 2002), are
$h_{i,0}(x(t),x_{l}(t),t)=\sum_{j\in\mathcal{N}_{i}}h_{ij}(x_{i}(t),x_{j}(t),t)+\sum_{l\in\mathcal{L}_{i}}\underline{h}_{il}(x_{i}(t),x_{l}(t),t)$
and
$h_{i,0}^{(\tau)}(x(t-\tau(t)),x_{l}(t-\tau(t)),t)=\sum_{j\in\mathcal{N}_{i}}h_{ij}^{(\tau)}(x_{i}(t-\tau(t)),x_{j}(t-\tau(t)),t)+\sum_{l\in\mathcal{L}_{i}}\underline{h}_{il}^{(\tau)}(x_{i}(t-\tau(t)),x_{l}(t-\tau(t)),t)$,
where $\mathcal{N}_{i}$ and $\mathcal{L}_{i}$ denote, respectively, the set of
neighbours of the $i$-th robot and the set of leaders to which the $i$-th
robot is connected. Also, the coupling functions are typically of the
diffusive type and model delayed and delay-free communications.

$\begin{split}&\bar{A}_{ii}(t)=\left[\begin{matrix}\frac{\partial
f_{i}(x_{i},t)}{\partial x_{i}}+\frac{\partial h_{i,0}(x,x_{l},t)}{\partial
x_{i}}&I_{n}&0_{n}&\cdots&0_{n}\\\ \frac{\partial h_{i,1}(x,x_{l},t)}{\partial
x_{i}}&0_{n}&I_{n}&\cdots&0_{n}\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
\frac{\partial h_{i,m-1}(x,x_{l},t)}{\partial
x_{i}}&0_{n}&0_{n}&\cdots&I_{n}\\\ \frac{\partial h_{i,m}(x,x_{l},t)}{\partial
x_{i}}&0_{n}&0_{n}&\cdots&0_{n}\end{matrix}\right],\\\
\bar{A}_{ij}(t)=&\left[\begin{matrix}\frac{\partial
h_{i,0}(x,x_{l},t)}{\partial x_{j}}&0_{n}&\cdots&0_{n}\\\
\vdots&\vdots&\ddots&\vdots\\\ \frac{\partial h_{i,m}(x,x_{l},t)}{\partial
x_{j}}&0_{n}&\cdots&0_{n}\end{matrix}\right],\bar{B}_{ij}(t)=\left[\begin{matrix}\frac{\partial
h_{i,0}^{(\tau)}(x,x_{l},t)}{\partial x_{j}}&0_{n}&\cdots&0_{n}\\\
\vdots&\vdots&\ddots&\vdots\\\ \frac{\partial
h_{i,m}^{(\tau)}(x,x_{l},t)}{\partial
x_{j}}&0_{n}&\cdots&0_{n}\end{matrix}\right].\end{split}$ (4)

###  3.2 Control goal

We let
$u(t)=[u_{1}^{\mathsf{T}}(t),\ldots,u_{N}^{\mathsf{T}}(t)]^{\mathsf{T}}$ be
the stack of the control inputs,
$d(t)=[d_{1}^{\mathsf{T}}(t),\ldots,d_{N}^{\mathsf{T}}(t)]^{\mathsf{T}}$ be
the stack of the disturbances,
$w(t)=[w_{1}^{\mathsf{T}}(t),\ldots,w_{N}^{\mathsf{T}}(t)]^{\mathsf{T}}$ be
the stack of the residual disturbances and
$\bar{d}(t)=[\bar{d}_{1}^{\mathsf{T}}(t),\ldots,\bar{d}_{N}^{\mathsf{T}}(t)]^{\mathsf{T}}$
be the stack of the polynomial disturbances. Our control goal is expressed in
terms of the so-called desired solution of the disturbance-free (or
unperturbed in what follows) network system (Monteil et al., 2019). The
desired solution is the solution of the network system characterized by
having: (i) the state of the agents keeping some desired configuration; (ii)
the multiplex layers giving no contribution to the $u_{i}$’s. This is
formalized next:

######  Definition 1.

Consider the network system (1) controlled by (3) with $d_{i}(t)=0$, $\forall
i$ and $\forall t$. We say that

  * •

$x^{\ast}(t):=[x_{1}^{\ast\mathsf{T}}(t),\dots,x_{N}^{\ast\mathsf{T}}(t)]^{\mathsf{T}}$,
$\dot{x}_{i}^{\ast}(t)=f_{i}(x_{i}^{\ast}(t),t)$ and
$x_{i}^{\ast}(s)=x_{i}^{\ast}(t_{0}),s\in[t_{0}-\tau_{\max},t_{0}]$, where
$x_{i}^{\ast}(t)$ is the desired state of the $i$-th agent;

  * •

$r_{i,k}^{\ast}(t):=0$, $\forall i$, $\forall k$ and $\forall t$;

is the desired solution of the system.

It is intrinsic in this definition that when the desired solution is achieved
it must hold that $u_{i}(t)=0$; note that this property is satisfied by e.g.
any diffusive-type control protocol and implicitly used in e.g. (Besselink and
Knorn, 2018; Mirabilio et al., 2021)222Nevertheless, our technical results can
be adapted to tackle situations where the desired solution is defined
differently, e.g. allowing for $u_{i}(t)\neq 0$ on the desired solution.. In
what follows, for the sake of brevity, we simply say that $x^{\ast}(t)$ is the
desired solution. Analogously, the desired output of (1) is
$y^{\ast}(t):=[y_{1}^{*\mathsf{T}}(t),\dots,y_{N}^{\ast\mathsf{T}}(t)]^{\mathsf{T}}$,
with ${y}^{\ast}_{i}(t)=g_{i}(x_{i}^{\ast}(t)),\forall i$ (leaving it implicit
that $r_{i,k}^{\ast}(t)=0$, $\forall i$, $\forall k$ and $\forall t$). We aim
at designing the control protocol (3) so that the closed-loop system rejects
the polynomial disturbances $\bar{d}(t)$ while guaranteeing that the residual
disturbance $w(t)$ is not amplified within the network system. This
requirement is captured with the definition of scalability with respect to
$w(t)$ formalized next:

######  Definition 2.

Consider the closed-loop system (1) - (3) with disturbance
$d(t)=w(t)+\bar{d}(t)$. The system is

  * •

$\mathcal{L}_{\infty}^{p}$-Input-to-State Scalable with respect to $w(t)$: if
there exists class $\mathcal{KL}$ functions $\alpha(\cdot,\cdot)$,
$\beta(\cdot,\cdot)$, a class $\mathcal{K}$ function $\gamma(\cdot)$, such
that for any initial condition and $\forall t\geq t_{0}$,

$\displaystyle\begin{split}&\max_{i}|x_{i}(t)-x_{i}^{\ast}(t)|_{p}\leq\gamma\left(\max_{i}\|w_{i}(\cdot)\|_{\mathcal{L}_{\infty}^{p}}\right)\\\
&+\beta\left(\max_{i}\sup_{t_{0}-\tau_{\max}\leq s\leq
t_{0}}\sum_{k=1}^{m}|r_{i,k}(s)+\bar{d}_{i}^{(k-1)}(s)|_{p},t-t_{0}\right)\\\
&+\alpha\left(\max_{i}\sup_{t_{0}-\tau_{\max}\leq s\leq
t_{0}}|x_{i}(s)-x_{i}^{\ast}(s)|_{p},t-t_{0}\right),\forall N;\end{split}$

  * •

$\mathcal{L}_{\infty}^{p}$-Input-Output Scalable with respect to $w(t)$: if
there exists class $\mathcal{KL}$ functions $\alpha(\cdot,\cdot)$,
$\beta(\cdot,\cdot)$, a class $\mathcal{K}$ function $\gamma(\cdot)$, such
that for any initial condition and $\forall t\geq t_{0}$,

$\displaystyle\begin{split}&\max_{i}|y_{i}(t)-y_{i}^{\ast}(t)|_{p}\leq\gamma\left(\max_{i}\|w_{i}(\cdot)\|_{\mathcal{L}_{\infty}^{p}}\right)\\\
&+\beta\left(\max_{i}\sup_{t_{0}-\tau_{\max}\leq s\leq
t_{0}}\sum_{k=1}^{m}|r_{i,k}(s)+\bar{d}_{i}^{(k-1)}(s)|_{p},t-t_{0}\right)\\\
&+\alpha\left(\max_{i}\sup_{t_{0}-\tau_{\max}\leq s\leq
t_{0}}|x_{i}(s)-x_{i}^{\ast}(s)|_{p},t-t_{0}\right),\forall N.\end{split}$

In the special case when $\bar{d}(t)=0$ and there are no multiplex layers,
Definition 2 becomes the definition for scalability given in (Xie et al.,
2021). In this context we note that the bounds in Definition 2 are uniform in
$N$ and this in turn guarantees that the residual disturbances are not
amplified within the network system. Intuitively, scalability is a stronger
property than stability. In what follows, whenever it is clear from the
context, we simply say that the network system is
$\mathcal{L}_{\infty}^{p}$-Input-to-State Scalable
($\mathcal{L}_{\infty}^{p}$-Input-Output Scalable) if Definition 2 is
fulfilled. In the special case where $p=2$ we simply say that the the network
system is $\mathcal{L}_{\infty}$-Input-to-State Scalable
($\mathcal{L}_{\infty}$-Input-Output Scalable).

##  4 Technical results

For the network system (1) we give a sufficient condition on the control
protocol (3) guaranteeing that the closed-loop system affected by disturbances
of the form (2) is $\mathcal{L}_{\infty}^{p}$-Input-to-State Scalable (see
Definition 2). We also give a corollary for $\mathcal{L}_{\infty}^{p}$-Input-
Output Scalability of the closed-loop system. The results are stated in terms
of the block diagonal coordinate transformation matrix
$T:=\textbf{diag}\\{T_{1},\ldots,T_{N}\\}\in\mathbb{R}^{N\cdot
n\cdot(m+1)\times N\cdot n\cdot(m+1)}$ with

$T_{i}:=\left[\begin{matrix}I_{n}&\alpha_{i,1}\cdot I_{n}&&&\\\
&I_{n}&\ddots&&\\\ &&\ddots&&\alpha_{i,m}\cdot I_{n}\\\
&&&&I_{n}\end{matrix}\right]\in\mathbb{R}^{n\cdot(m+1)\times n\cdot(m+1)},$

where $\alpha_{i,k}\in\mathbb{R},\forall k$.

######  Proposition 1.

Consider the closed-loop network system (1) - (3) with $y_{i}(t)=x_{i}(t)$
affected by disturbances (2). Let the matrices
$\bar{A}_{ii}(t),\bar{A}_{ij}(t),\bar{B}_{ij}(t)$ be defined as in (4) and
assume that, $\forall t\geq t_{0}$, the following set of conditions are
satisfied for some $0\leq\underline{\sigma}<\bar{\sigma}<+\infty$:

  * C1

$h_{i,k}(x^{\ast},x_{l},t)=h_{i,k}^{(\tau)}(x^{\ast},x_{l},t)=0$, $\forall
i,\forall k$;

  * C2

$\mu_{p}(T_{i}\bar{A}_{ii}(t)T_{i}^{-1})+\sum_{j\neq
i}\|T_{i}\bar{A}_{ij}(t)T_{j}^{-1}\|_{p}\leq-\bar{\sigma}$, $\forall i$ and
$\forall x\in\mathbb{R}^{nN},\forall x_{l}\in\mathbb{R}^{nM}$;

  * C3

$\sum_{j}\|T_{i}\bar{B}_{ij}(t)T_{j}^{-1}\|_{p}\leq\underline{\sigma}$,
$\forall i$ and $\forall x\in\mathbb{R}^{nN},\forall x_{l}\in\mathbb{R}^{nM}$.

Then, the system is $\mathcal{L}_{\infty}^{p}$-Input-to-State Scalable. In
particular,

$\displaystyle\begin{split}&\max_{i}|x_{i}(t)-x_{i}^{\ast}(t)|_{p}\leq\frac{\kappa_{G}(T)}{\bar{\sigma}-\underline{\sigma}}\max_{i}\|w_{i}(\cdot)\|_{\mathcal{L}_{\infty}^{p}}\\\
&+\kappa_{G}(T)e^{-\hat{\lambda}(t-t_{0})}\max_{i}\sup_{t_{0}-\tau_{\max}\leq
s\leq t_{0}}|x_{i}(s)-x_{i}^{\ast}(s)|_{p}\\\
&+\kappa_{G}(T)e^{-\hat{\lambda}(t-t_{0})}\max_{i}\sup_{t_{0}-\tau_{\max}\leq
s\leq t_{0}}\sum_{k=1}^{m}|r_{i,k}(s)\\\
&+\sum_{b=0}^{m-k}\frac{(m-1-b)!}{(m-k-b)!}\cdot\bar{d}_{i,m-1-b}\cdot
s^{m-k-b}|_{p},\forall N,\end{split}$ (5)

where $\kappa_{G}(T):=\|T\|_{G}\|T^{-1}\|_{G}$ and

$\hat{\lambda}=\inf_{t\geq
t_{0}}\\{\lambda|\lambda(t)-\bar{\sigma}+\underline{\sigma}e^{\lambda(t)\tau(t)}=0\\}.$
(6)

The proof of the result is given in Appendix A and we now make the following
considerations on the conditions.

######  Remark 4.

Proposition 1 implies that, in order to reject polynomial disturbances of
order up to ${m-1}$, one needs to design a control protocol with $m$ multiplex
layers. See the last term of the upper bound in (5).

######  Remark 5.

Condition $C1$ implies that $u_{i}(t)=0$ at the desired solution. This rather
common condition (see e.g. (di Bernardo et al., 2015; Xie et al., 2021))
guarantees that $x^{\ast}(t)$ is a solution of the unperturbed dynamics. This
assumption is satisfied in e.g. all consensus/synchronization dynamics with
diffusive-type couplings. Condition $C2$, giving an upper bound on the matrix
measure of the Jacobian of the delay-free part of the closed-loop network
dynamics, is a diagonal dominance condition. Instead, condition $C3$ gives an
upper bound on the norm of the Jacobian of the dynamics containing delays.
Intuitively, Proposition 1 implies that the matrix measure should be negative
enough to compensate the presence of delays.

######  Remark 6.

Conditions $C2$ and $C3$ can be leveraged to shape the coupling functions
between the agents and to determine the maximum number of neighbours for each
agent in the network. Moreover, if $C2$ and $C3$ are satisfied, then the
network is also connective stable in the sense of (Siljak, 2011, Chapter 2.1).
Intuitively, a network is connective stable if the removal of couplings
preserves stability. This property and the related $\gamma$-scalability from
(Knorn and Besselink, 2020) have been investigated for delay-free networks. In
Section 5 we show that $C2$ and $C3$ can be effectively checked by recasting
the fulfillment of these conditions into an optimisation problem that allows
to design the control protocol for each agent independently on the other
agents.

######  Remark 7.

Following Proposition 1, the convergence rate of the closed-loop network is at
least $\hat{\lambda}$. We note that $\hat{\lambda}$ depends on $\tau(t)$ and
(6) highlights that the larger $\tau_{\max}$, the lower $\hat{\lambda}$. In
particular, $\hat{\lambda}=\bar{\sigma}$ when $\tau_{\max}=0$ (i.e. when there
are no delays) and decreases as the delay increases. From the design
viewpoint, the expression in (6) can be used to design the protocols so that
the network exhibits a given, desired, convergence rate.

######  Remark 8.

Proposition 1 generalises a number of results in the literature. Specifically,
our result yields the sufficient condition for string stability in (Silva et
al., 2021) when the network topology is a string, $\bar{d}_{i}(t)$ is constant
$\forall i$ and there are no delays. That is, we extend the results in (Silva
et al., 2021; Monteil et al., 2019) by considering dynamic compensation of
polynomial disturbances, general network topologies and delays. We also extend
our prior work (Xie et al., 2021) in which scalability is considered without
rejection of polynomial disturbances.

Before illustrating (Section 5) how our approach can be effectively used for
protocol design, we give here the next result that immediately follows from
Proposition 1.

######  Corollary 1.

Consider the closed-loop network system (1) - (3) affected by disturbances
(2). Assume that all the conditions in Proposition 1 are satisfied and that,
in addition, the output functions $g_{i}(\cdot)$ are Lipschitz. Then the
system is $\mathcal{L}_{\infty}^{p}$-Input-Output Scalable.

The proof, directly following from the Lipschitz hypothesis on $g_{i}(\cdot)$
and from the upper bound in (5), is omitted here for brevity. ∎

##  5 Using Proposition 1 to design multiplex protocols for formation
scalability

We consider the problem of designing a control protocol guaranteeing that a
network of $N$ unicycle robots is $\mathcal{L}_{\infty}$-Input-to-State
Scalable. Our goal is to design a protocol allowing the formation to: (i)
track a reference provided by a leader; (ii) reject certain polynomial
disturbances; (iii) ensure the non-amplification of residual disturbances.
Before presenting the validation results, we describe here the robot dynamics
and how Proposition 1 was used to synthesize the protocol. The results from
the in-silico and experimental validations are presented in Section 6. We made
use of the Robotarium (Wilson et al., 2020) for our experimental
validation333The Robotarium by Georgia Tech is an open research infrastructure
to perform experiments on real robots. Experiments can be launched remotely
via a web interface..

#### Agent dynamics.

We consider the following dynamics for unicycle robots (see e.g. (Wilson et
al., 2020) and references therein):

$\displaystyle\begin{split}\dot{p}_{i}^{x}(t)&=v_{i}(t)\cos\theta_{i}(t)+d_{i}^{x}(t),\\\
\dot{p}_{i}^{y}(t)&=v_{i}(t)\sin\theta_{i}(t)+d_{i}^{y}(t),\\\
\dot{\theta}_{i}(t)&=\Omega_{i}(t),\end{split}$ (7)

$\forall i$, where the state variables
$p_{i}(t)=[p_{i}^{x}(t),p_{i}^{y}(t)]^{\mathsf{T}}$ represent the inertial
position of the $i$-th robot and $\theta_{i}(t)$ is its heading angle. The
control input is denoted by $u_{i}(t)=[v_{i}(t),\Omega_{i}(t)]^{\mathsf{T}}$
with $v_{i}(t)$ being the linear velocity and $\Omega_{i}(t)$ being the
angular velocity. We consider the case where the disturbances affecting the
robots, i.e. $d_{i}(t)=[d_{i}^{x}(t),d_{i}^{y}(t)]^{\mathsf{T}}$, are of the
form $d_{i}^{x}(t):=\bar{d}_{i,0}^{x}+\bar{d}_{i,1}^{x}\cdot t+w_{i}^{x}(t)$
and $d_{i}^{y}(t):=\bar{d}_{i,0}^{y}+\bar{d}_{i,1}^{y}\cdot t+w_{i}^{y}(t)$.
The terms $\bar{d}_{i,0}^{x}$ and $\bar{d}_{i,0}^{y}$ model constant
disturbances on the velocity of the robots. These disturbances also naturally
arise in the context of unicycle-like marine robots, the dynamics of which is
also captured by (7). For these robots, the constant terms in $d_{i}(t)$ model
the disturbances due to the ocean current (Panagou and Kyriakopoulos, 2011)
and $w_{i}^{x}(t)$, $w_{i}^{y}(t)$ model e.g. transient variations of the
current. The ramp terms in the disturbance can instead model ramp attack
signals (Sridhar and Govindarasu, 2014). In what follows, we make use of the
compact notation $w_{i}(t):=[w_{i}^{x}(t),w_{i}^{y}(t)]^{\mathsf{T}}$,
$\bar{d}_{i}(t):=[\bar{d}_{i,0}^{x}+\bar{d}_{i,1}^{x}\cdot
t,\bar{d}_{i,0}^{y}+\bar{d}_{i,1}^{y}\cdot t]^{\mathsf{T}}$. Following (Lawton
et al., 2003), the dynamics for the robot hand position is given by

$\displaystyle\dot{\eta}_{i}(t)={\left[\begin{matrix}\cos\theta_{i}(t)&-l_{i}\sin\theta_{i}(t)\\\
\sin\theta_{i}(t)&l_{i}\cos\theta_{i}(t)\end{matrix}\right]}u_{i}(t)+d_{i}(t),$
(8)

where $l_{i}\in\mathbb{R}_{>0}$ is the distance of the hand position to the
wheel axis. The dynamics can be feedback linearised by

$\displaystyle
u_{i}(t)={\left[\begin{matrix}\cos\theta_{i}(t)&-l_{i}\sin\theta_{i}(t)\\\
\sin\theta_{i}(t)&l_{i}\cos\theta_{i}(t)\end{matrix}\right]^{-1}}\nu_{i}(t),$

which yields

$\displaystyle\dot{\eta}_{i}(t)=\nu_{i}(t)+d_{i}(t),\ \ \ \forall i$ (9)

Next we leverage Proposition 1 to design $\nu_{i}(t)$ so that network (9) is
$\mathcal{L}_{\infty}$-Input-to-State Scalable.

#### Protocol design.

We denote by $\eta_{l}(t)$ the hand position provided by a virtual leader.
Robots are required to keep a desired offset from the leader
($\delta_{li}^{\ast}$) and from neighbours ($\delta_{ji}^{\ast}$) while
tracking a reference velocity from the leader, say $v_{l}(t)$. Following these
requirements, we pick the desired solution (see Section 3.2) for the $i$-th
robot dynamics, say $\eta_{i}^{\ast}(t)$, so that: (i) the robot keeps the
desired offsets from the leader and from the neighbours, i.e.
$\eta_{l}^{\ast}(t)-\eta_{i}^{\ast}(t)=\delta_{li}^{\ast}$ and
$\eta_{j}^{\ast}(t)-\eta_{i}^{\ast}(t)=\delta_{ji}^{\ast}$; (ii) the reference
velocity is tracked, i.e. $\dot{\eta}_{i}^{\ast}(t)=v_{l}(t)$. Inspired by
(Liu et al., 2019; Li et al., 2013), we consider protocols of the form:

$\displaystyle\nu_{i}(t)=\bar{\nu}_{i}(t)+v_{l}(t)$ (10)

with

$\displaystyle\begin{split}\bar{\nu}_{i}(t)&=h_{i,0}(\eta(t),\eta_{l}(t),t)\\\
&+h_{i,0}^{(\tau)}(\eta(t-\tau(t)),\eta_{l}(t-\tau(t)),t)+r_{i,1}(t),\\\
\dot{r}_{i,1}(t)&=h_{i,1}(\eta(t),\eta_{l}(t),t)\\\
&+h_{i,1}^{(\tau)}(\eta(t-\tau(t)),\eta_{l}(t-\tau(t)),t)+r_{i,2}(t),\\\
\dot{r}_{i,2}(t)&=h_{i,2}(\eta(t),\eta_{l}(t),t)\\\
&+h_{i,2}^{(\tau)}(\eta(t-\tau(t)),\eta_{l}(t-\tau(t)),t).\end{split}$ (11)

In the above protocol, $r_{i,1}(t)$ and $r_{i,2}(t)$ are generated by two
multiplex layers (say, layer $1$ and layer $2$). With this protocol, the
closed-loop dynamics is given by

$\displaystyle\dot{\eta}_{i}(t)=v_{l}(t)+\bar{\nu}_{i}(t)+d_{i}(t),\ \ \
\forall i$ (12)

Hence, we recast the problem of designing the protocol $\nu_{i}(t)$ for (9)
into the problem of designing $\bar{\nu}_{i}(t)$ in (11) for the dynamics
(12). To this aim, we can apply Proposition 1. The coupling functions in (11)
are (time dependence inside these functions omitted in what follows):

$\displaystyle\begin{split}h_{i,0}(\eta,\eta_{l},t)&=k_{0}(\eta_{l}-\eta_{i}-\delta_{li}^{\ast}),\\\
h_{i,0}^{(\tau)}(\eta,\eta_{l},t)&=k_{0}^{(\tau)}\sum_{j\in\mathcal{N}_{i}}\psi(\eta_{j}-\eta_{i}-\delta_{ji}^{\ast}),\\\
h_{i,1}(\eta,\eta_{l},t)&=k_{1}(\eta_{l}-\eta_{i}-\delta_{li}^{\ast}),\\\
h_{i,1}^{(\tau)}(\eta,\eta_{l},t)&=k_{1}^{(\tau)}\sum_{j\in\mathcal{N}_{i}}\psi(\eta_{j}-\eta_{i}-\delta_{ji}^{\ast}),\\\
h_{i,2}(\eta,\eta_{l},t)&=k_{2}(\eta_{l}-\eta_{i}-\delta_{li}^{\ast}),\\\
h_{i,2}^{(\tau)}(\eta,\eta_{l},t)&=k_{2}^{(\tau)}\sum_{j\in\mathcal{N}_{i}}\psi(\eta_{j}-\eta_{i}-\delta_{ji}^{\ast}),\end{split}$
(13)

with $\psi(x):=\tanh(k^{\psi}x)$ being inspired from (Monteil et al., 2019).
In the above expression, $\mathcal{N}_{i}$ is the set of neighbours of robot
$i$. The cardinality of $\mathcal{N}_{i}$ is bounded, that is,
$|\mathcal{N}_{i}|\leq\bar{N},\forall i$. Also, the control gains $k_{0},$
$k_{1},$ $k_{2},$ $k_{0}^{(\tau)},k_{1}^{(\tau)},k_{2}^{(\tau)},k^{\psi}$ are
non-negative scalars designed next.

$\begin{gathered}\bar{A}_{ii}:=\left[\begin{matrix}-k_{0}I_{2}&I_{2}&0_{2}\\\
-k_{1}I_{2}&0_{2}&I_{2}\\\
-k_{2}I_{2}&0_{2}&0_{2}\end{matrix}\right],\tilde{B}_{ii}:=\left[\begin{matrix}-\bar{N}\bar{g}_{0}I_{2}&0_{2}&0_{2}\\\
-\bar{N}\bar{g}_{1}I_{2}&0_{2}&0_{2}\\\
-\bar{N}\bar{g}_{2}I_{2}&0_{2}&0_{2}\end{matrix}\right],\tilde{B}_{ij}:=\left[\begin{matrix}\bar{g}_{0}I_{2}&0_{2}&0_{2}\\\
\bar{g}_{1}I_{2}&0_{2}&0_{2}\\\
\bar{g}_{2}I_{2}&0_{2}&0_{2}\end{matrix}\right],\bar{T}:=\left[\begin{matrix}I_{2}&\alpha_{1}I_{2}&0_{2}\\\
0_{2}&I_{2}&\alpha_{2}I_{2}\\\
0_{2}&0_{2}&I_{2}\end{matrix}\right].\end{gathered}$ (14)

Specifically, we make use of Proposition 1 to select the control gains so that
the robotic network is $\mathcal{L}_{\infty}$-Input-to-State Scalable. In
particular, we note that the choice of the control protocol (10) with coupling
functions (13) guarantees the fulfillment of $C1$. Hence, we now only need to
find a set of control gains that fulfills $C2$ and $C3$. To this aim, in
Proposition 1, we pick all transformation matrices to be the same, say
$T_{i}=T_{j}=\bar{T}$ $\forall i,j$. Then, as shown in Appendix B, the problem
of finding control gains that fulfill these two conditions can be recast as
the following optimisation problem:

$\displaystyle\begin{split}&\min_{\xi}\ \mathcal{J}\\\ &s.t.\quad k_{0}\geq
0,k_{1}\geq 0,k_{2}\geq 0,\bar{g}_{0}\geq 0,\bar{g}_{1}\geq 0,\bar{g}_{2}\geq
0,\\\
&\phantom{s.t.\quad}k_{0}+\bar{g}_{0}>0,k_{1}+\bar{g}_{1}>0,k_{2}+\bar{g}_{2}>0,\\\
&\phantom{s.t.\quad}\bar{\sigma}>0,\underline{\sigma}\geq
0,\bar{\sigma}-\underline{\sigma}>0,[\bar{T}\bar{A}_{ii}\bar{T}^{-1}]_{s}\preceq-\bar{\sigma}I_{6},\\\
&\phantom{s.t.\quad}\left[\begin{matrix}\frac{\underline{\sigma}}{2\bar{N}}I_{6}&(\bar{T}\tilde{B}_{ij}\bar{T}^{-1})^{\mathsf{T}}\\\
\bar{T}\tilde{B}_{ij}\bar{T}^{-1}&\frac{\underline{\sigma}}{2\bar{N}}I_{6}\end{matrix}\right]\succeq
0,\\\
&\phantom{s.t.\quad}\left[\begin{matrix}\frac{\underline{\sigma}}{2}I_{6}&(\bar{T}\tilde{B}_{ii}\bar{T}^{-1})^{\mathsf{T}}\\\
\bar{T}\tilde{B}_{ii}\bar{T}^{-1}&\frac{\underline{\sigma}}{2}I_{6}\end{matrix}\right]\succeq
0.\end{split}$ (15)

The decision variables are
$\xi:=[k_{0},k_{1},k_{2},\bar{g}_{0},\bar{g}_{1},\bar{g}_{2},\allowbreak\underline{\sigma},\bar{\sigma}]$
with $\bar{g}_{0}=k^{\psi}k_{0}^{(\tau)},\bar{g}_{1}=k^{\psi}k_{1}^{(\tau)}$
and $\bar{g}_{2}=k^{\psi}k_{2}^{(\tau)}$. The matrices in the above expression
are defined in (14). In (15) we decided to use the cost
$\mathcal{J}:=-\bar{g}_{0}-\bar{g}_{1}-\bar{g}_{2}$, which was chosen in
accordance to (Monteil et al., 2019) with the aim of maximizing the upper
bound of the inter-robot coupling functions (other cost functions could be
chosen as the steps described in Appendix B formalizing the optimisation
problem are not dependent on $\mathcal{J}$).

##  6 Validation

We now validate the effectiveness of the protocol obtained in Section 5. We do
this by showing that: (i) robots are able to keep a given desired formation;
(ii) the protocol prohibits the amplification of residual disturbances; (iii)
the polynomial disturbances are compensated by the control protocol. Robots
need to keep a formation consisting of concentric circles (the $k$-th circle
consists of $4k$ robots) and their hand positions need to move following a
reference trajectory. Robots receive the velocity and position signals from
the virtual leader and each robot is connected to a maximum of $\bar{N}=3$
neighbours (i.e. the closest robots). Specifically, a given robot on the
$k$-th circle is connected to the robots immediately ahead and behind on the
same circle and with the closest robot on circle $k-1$ (if any). The set-up we
consider in our validation experiments is schematically illustrated in Figure
2, together with the reference trajectory from the virtual leader.

Figure 2: The set-up considered in the experiments of Section 6. The reference
trajectory of the hand position provided by the virtual leader is shown
together with an example of desired formation consisting of $3$ concentric
circles.

In the figure, for clarity, $3$ concentric circles are shown. We used a delay
of $\tau(t)=0.33$s in both our simulations and hardware experiments. We first
illustrate the results from the simulations and then the results obtained from
the experiments on the Robotarium. The code and data to replicate all the
experiments of Section 6 can be found at https://tinyurl.com/4wyacf7z.

#### In-silico validation.

We consider a formation of $30$ circles, with two robots on circle $1$ (say,
robot $1$ and robot $3$) affected by disturbances:

$\begin{split}d_{1}(t)&=\left[\begin{array}[]{*{20}c}0.04+0.4\sin(0.5t)e^{-0.1t}\\\
0.04+0.4\sin(0.5t)e^{-0.1t}\end{array}\right],\\\
d_{3}(t)&=\left[\begin{array}[]{*{20}c}-0.05t+0.4\sin(0.5t)e^{-0.1t}\\\
-0.05t+0.4\sin(0.5t)e^{-0.1t}\end{array}\right].\end{split}$ (16)

We computed the gains of the control protocol in (10) - (13) by solving the
optimisation problem in (15) for a grid of parameters $\alpha_{1}$ and
$\alpha_{2}$. We then selected the control parameters as the ones returning
the lowest cost for each fixed pair of $\alpha$’s. By doing so, we obtained
the gains $k_{0}=1.4155$, $k_{1}=1.5103$, $k_{2}=0.4803$,
$k_{0}^{(\tau)}=0.642$, $k_{1}^{(\tau)}=0.872$,
$k_{2}^{(\tau)}=0.425,k^{\psi}=0.1$ (corresponding to
$\alpha_{1}=-0.6,\alpha_{2}=-1.6$). In Figure 3 the maximum hand position
deviation is shown when the number of robots in the formation is increased,
starting with a formation of $1$ circle only to a formation with $30$ circles
(i.e. $1860$ robots).

Figure 3: Maximum hand position deviation (in meters) of each circle across as
number of circles increases from $1$ to $30$.

The figure was obtained by starting with a formation of $1$ circle and
increasing at each simulation the number of circles. We recorded at each
simulation the maximum hand position deviation for each robot on a given
circle and finally plotted the largest deviation on each circle across all the
simulations. The figure clearly shows that the polynomial components of the
disturbances in (16) are compensated by the integral control protocol and the
residual disturbances are not amplified through the formation. We also report
the behaviour of the full formation with $30$ circles in Figure 4.

Figure 4: Top panel: hand position deviations of all the robots (in meters).
Bottom panel: hand position deviations of the unperturbed robots only (in
meters). Robots on the same circle have the same color. Disturbances are the
ones in (16).

Both panels of the figure confirm that, in accordance with our theoretical
results, the protocol allows the robots to keep the desired formation, while
compensating the polynomial components of the disturbances and prohibiting the
amplification of the residual disturbances.

#### Experimental validation.

We further validate our results by carrying out experiments on the Robotarium,
which provides both hardware infrastructure and a high-fidelity simulator of
the hardware. In the experiments, a formation of $2$ concentric circles (hence
with $12$ robots) is considered and, for consistency with our previous set of
simulations, $2$ robots on circle $1$ are perturbed by the disturbances
$d_{1}(t)$ and $d_{3}(t)$ given in (16). The Robotarium documentation
(Robotarium, 2022) reports a nominal step size (for both the simulator and the
hardware infrastructure) of $0.033$s. Since the step size is used to implement
the multiplex layers of the integral control protocol in (10) - (13) as a
first step we measured the variability of the step size in the hardware
infrastructure. The result is given in Figure 5.

Figure 5: Measurement of the step size on the Robotarium hardware. Solid
markers represent the average step size from $10$ sets of experiments while
the shaded area represents the confidence interval corresponding to the
standard deviation.

Such a figure reports the average step size we measured using built-in timing
functions across $10$ experiments444See our code at
https://tinyurl.com/4wyacf7z for the details on how these measurements were
performed.. In the same figure, the shaded area represents the confidence
interval corresponding to the standard deviation. As illustrated in the
figure, while the average step size is indeed around $0.033$s and consistent
with the nominal value, it also introduces some variability in the
experiments. To mitigate this variability and take this implementation aspect
into account when deploying our control protocol on the hardware platform, we
imposed the control gains of the multiplex layers (i.e. $1$ and $2$) to be
smaller than the gains of layer $0$. This was done by solving the optimisation
problem in (15) this time with the following additional constraints:
$k_{0}\geq 2k_{1},k_{0}\geq 2k_{2},\bar{g}_{0}\geq
2\bar{g}_{1},\bar{g}_{0}\geq 2\bar{g}_{2}$. As an outcome of this process, we
obtained the following gains: $k_{0}=1.2674$, $k_{1}=0.6312$, $k_{2}=0.133$,
$k_{0}^{(\tau)}=0.325$, $k_{1}^{(\tau)}=0.162$,
$k_{2}^{(\tau)}=0.06,k^{\psi}=0.1$ (which correspond to
$\alpha_{1}=-1.1,\alpha_{2}=-2.6$). We then validated the control protocol
with this choice of parameters by first leveraging the Robotarium simulator
and the results, consistent with our theoretical findings, are shown in the
top panel of Figure 6.

Figure 6: Top panel: hand position deviations of the robots (in meters) from
the Robotarium simulator. Bottom panel: hand position deviations (in meters)
from the hardware experiments. Solid lines represent the average hand position
deviations across $10$ sets of experiments; the shaded area represents the
confidence interval corresponding to the standard deviation. Disturbances are
the ones in (16). See https://tinyurl.com/4wyacf7z for the animated version of
the plot in the bottom panel and for the recording of the experiments from
both the simulator and hardware set-up.

Next, we validated the control protocol on the Robotarium hardware
infrastructure and the outcome from these experiments are shown in the bottom
panel of Figure 6. In the figure, which was obtained from a set of $10$
experiments, the solid lines are the robots’ average hand position deviations
and the shaded area represents the confidence interval corresponding to the
standard deviation. The behaviour of the hardware experiments is in agreement
with the one obtained from the simulator. Both panels show that, in accordance
with Proposition 1, our control protocol allowed the robots to keep the
desired formation, while compensating the polynomial components of the
disturbances and ensuring non-amplification of the residual disturbances.

##  7 Conclusions and future work

We considered the problem of designing distributed control protocols for
possibly nonlinear networks affected by delays that not only allow the network
to achieve some desired behaviour, but also, by means of integral actions
delivered via multiplex architecture, guarantee a scalability property of the
network. This property implies the rejection of polynomial disturbances and
the non-amplification of the residual disturbances. To tackle this problem, we
presented a set of sufficient conditions for scalability enabling the design
of integral multiplex protocols for both leader-follower and leaderless
networks consisting of possibly heterogeneous agents affected by communication
delays. After turning the results into an optimisation problem, we
experimentally evaluated the effectiveness of our approach via in-silico and
hardware validations. Besides considering networks with heterogeneous delays,
we plan to build on the results of this paper to study disturbance rejection
in discrete-time networks with non-uniform sampling. We will also explore the
effectiveness of our results on the design of integral actions for biochemical
networks (Qian and Del Vecchio, 2018).

##  8 Acknowledgments

The authors would like to thank Prof. Francesco Bullo for the insightful
discussions on an early version of these results. The inputs provided by Prof.
Bullo helped us to improve the quality of our paper.

##  Appendix A Proof of Proposition 1

We start with augmenting the state of the original dynamics by defining

$z_{i}(t):=[x_{i}^{\mathsf{T}}(t),\zeta_{i,1}^{\mathsf{T}}(t),\zeta_{i,2}^{\mathsf{T}}(t),\cdots,\zeta_{i,m}^{\mathsf{T}}(t)]^{\mathsf{T}},$

and where

$\displaystyle\begin{split}\zeta_{i,k}(t)&=r_{i,k}(t)+\sum_{b=0}^{m-k}\frac{(m-1-b)!}{(m-k-b)!}\cdot\bar{d}_{i,m-1-b}\cdot
t^{m-k-b},\end{split}$

$k=1,\cdots,m$. In these new coordinates the dynamics of the network system
becomes

$\displaystyle\dot{z}_{i}(t)=\tilde{f}_{i}(z_{i},t)+\tilde{v}_{i}(z,t)+\tilde{w}_{i}(t),$
(17)

where $\tilde{f}_{i}(z_{i},t)=[f_{i}^{\mathsf{T}}(x_{i},t),0_{1\times
n},\ldots,0_{1\times n}]^{\mathsf{T}}$,
$\tilde{w}_{i}(t)=[w_{i}^{\mathsf{T}}(t),0_{1\times n},\ldots,0_{1\times
n}]^{\mathsf{T}}$, $\tilde{v}_{i}(z,t)=v_{i}(z,t)+v_{i}^{(\tau)}(z,t)$ with

$\displaystyle\begin{split}v_{i}(z,t)=\left[\begin{array}[]{*{20}c}h_{i,0}(x(t),x_{l}(t),t)+\zeta_{i,1}(t)\\\
h_{i,1}(x(t),x_{l}(t),t)+\zeta_{i,2}(t)\\\ \vdots\\\
h_{i,m-1}(x(t),x_{l}(t),t)+\zeta_{i,m}(t)\\\
h_{i,m}(x(t),x_{l}(t),t)\end{array}\right],\end{split}$

and

$\displaystyle\begin{split}v_{i}^{(\tau)}(z,t)=\left[\begin{array}[]{*{20}c}h_{i,0}^{(\tau)}(x(t-\tau(t)),x_{l}(t-\tau(t)),t)\\\
h_{i,1}^{(\tau)}(x(t-\tau(t)),x_{l}(t-\tau(t)),t)\\\ \vdots\\\
h_{i,m}^{(\tau)}(x(t-\tau(t)),x_{l}(t-\tau(t)),t)\end{array}\right].\end{split}$

Condition $C1$ implies that $x_{i}^{\ast}(t)$ is a solution of the unperturbed
network dynamics, i.e. $x_{i}^{\ast}(t)$ is a solution of (1) when there are
no disturbances. Moreover, when there are no disturbances, in the augmented
dynamics, the solution $z_{i}^{\ast}(t):=[x_{i}^{\ast\mathsf{T}}(t),0_{1\times
n},\ldots,0_{1\times n}]^{\mathsf{T}}$ satisfies
$\dot{z}^{\ast}_{i}(t)=\tilde{f}_{i}(z_{i}^{\ast},t),$ with
$\tilde{f}_{i}(z_{i}^{\ast},t):=[f_{i}^{\mathsf{T}}(x_{i}^{\ast},t),0_{1\times
n},\ldots,0_{1\times n}]^{\mathsf{T}}$. Hence, the dynamics of state deviation
(i.e. the error) $e_{i}(t)=z_{i}(t)-z_{i}^{\ast}(t)$ is given by

$\displaystyle\begin{split}\dot{e}_{i}(t)&=\tilde{f}_{i}(z_{i},t)-\tilde{f}_{i}(z_{i}^{\ast},t)+\tilde{v}_{i}(z,t)+\tilde{w}_{i}(t).\end{split}$
(18)

Following (Desoer and Haneda, 1972), we let $\eta_{i}(\rho)=\rho
z_{i}+(1-\rho)z_{i}^{\ast}$,
$\eta(\rho)=[\eta_{1}^{\mathsf{T}}(\rho),\ldots,\eta_{N}^{\mathsf{T}}(\rho)]^{\mathsf{T}}$
and then rewrite the error dynamics as

$\displaystyle\dot{e}(t)=A(t)e(t)+B(t)e(t-\tau(t))+\tilde{w}(t),$ (19)

where
$\tilde{w}=[\tilde{w}_{1}^{\mathsf{T}}(t),\ldots,\tilde{w}_{N}^{\mathsf{T}}(t)]^{\mathsf{T}}$
and $A(t)$ has entries: (i)
$A_{ii}(t)=\int_{0}^{1}\left(J_{\tilde{f}_{i}}(\eta_{i}(\rho),t)+J_{v_{i}}(\eta_{i}(\rho),t)\right)d\rho$;
(ii) $A_{ij}(t)=\int_{0}^{1}J_{v_{i}}(\eta_{j}(\rho),t)d\rho$. Similarly,
$B(t)$ has entries:
$B_{ij}(t)=\int_{0}^{1}J_{v_{i}^{(\tau)}}(\eta_{j}(\rho),t)d\rho$. In the
above expressions, the Jacobian matrices are defined as
$J_{\tilde{f}_{i}}(\eta_{i},t):=\frac{\partial\tilde{f}_{i}(\eta_{i},t)}{\partial\eta_{i}}$,
$J_{v_{i}}(\eta_{i},t):=\frac{\partial v_{i}(\eta,t)}{\partial\eta_{i}}$,
$J_{v_{i}^{(\tau)}}(\eta_{i},t):=\frac{\partial
v_{i}^{(\tau)}(\eta,t)}{\partial\eta_{i}}$. Now, consider the coordinate
transformation $\tilde{z}(t):=Tz(t)$ and $\tilde{e}(t):=Te(t)$, we have

$\displaystyle\dot{\tilde{e}}(t)=TA(t)T^{-1}\tilde{e}(t)+TB(t)T^{-1}\tilde{e}(t-\tau(t))+T\tilde{w}(t).$
(20)

Let $|x|_{G}:=\left||x_{1}|_{p},\ldots,|x_{N}|_{p}\right|_{\infty}$. Then, by
taking the Dini derivative of $|\tilde{e}(t)|_{G}$ we may continue as follows

$\displaystyle\begin{split}&D^{+}|\tilde{e}(t)|_{G}=\limsup_{h\rightarrow
0^{+}}\frac{1}{h}\left(|\tilde{e}(t+h)|_{G}-|\tilde{e}(t)|_{G}\right)\\\
&=\limsup_{h\rightarrow
0^{+}}\frac{1}{h}\bigl{(}|\tilde{e}(t)+hTA(t)T^{-1}\tilde{e}(t)+hTB(t)T^{-1}\tilde{e}(t-\tau(t))\\\
&\phantom{=}+hT\tilde{w}(t)|_{G}-|\tilde{e}(t)|_{G}\bigr{)}\\\
&\leq\limsup_{h\rightarrow
0^{+}}\frac{1}{h}\left(\|I+hTA(t)T^{-1}\|_{G}-1\right)|\tilde{e}(t)|_{G}+|T\tilde{w}(t)|_{G}\\\
&\phantom{=}+\|TB(t)T^{-1}\|_{G}|\tilde{e}(t-\tau(t))|_{G}\\\
&\leq\mu_{G}(TA(t)T^{-1})|\tilde{e}(t)|_{G}+\|TB(t)T^{-1}\|_{G}\sup_{t-\tau_{\max}\leq
s\leq t}|\tilde{e}(s)|_{G}\\\ &\phantom{=}+|T\tilde{w}(t)|_{G}\\\
&\leq\mu_{G}(TA(t)T^{-1})|\tilde{e}(t)|_{G}+\|TB(t)T^{-1}\|_{G}\sup_{t-\tau_{\max}\leq
s\leq t}|\tilde{e}(s)|_{G}\\\
&\phantom{=}+\|T\|_{G}\max_{i}\|\tilde{w}_{i}(\cdot)\|_{\mathcal{L}_{\infty}^{p}}.\end{split}$

Next, we find upper bounds for $\mu_{G}(TA(t)T^{-1})$ and
$\|TB(t)T^{-1}\|_{G}$ which allow us to apply Lemma 2. First, we give the
expression of the matrix $\bar{A}(t)$ which have entries: (i)
$\bar{A}_{ii}(t)=J_{\tilde{f}_{i}}(z_{i},t)+J_{v_{i}}(z_{i},t)$; (ii)
$\bar{A}_{ij}(t)=J_{v_{i}}(z_{j},t)$, and $\bar{B}(t)$ has entries:
$\bar{B}_{ij}(t)=J_{v_{i}^{(\tau)}}(z_{j},t)$. Then, by sub-additivity of
matrix measures and matrix norms, we get
$\mu_{G}(TA(t)T^{-1})\leq\int_{0}^{1}\mu_{G}(T\bar{A}(t)T^{-1})d\rho$ and
$\|TB(t)T^{-1}\|_{G}\leq\int_{0}^{1}\|T\bar{B}(t)T^{-1}\|_{G}d\rho$ (see also
Lemma $3.4$ in (Russo and Wirth, 2022)). Moreover, from Lemma 1 it then
follows that

$\displaystyle\begin{split}&\mu_{G}(T\bar{A}(t)T^{-1})\leq\\\
&\max_{i}\Bigl{\\{}\mu_{p}(T_{i}\bar{A}_{ii}(t)T_{i}^{-1})+\sum_{j\neq
i}\|T_{i}\bar{A}_{ij}(t)T_{j}^{-1}\|_{p}\Bigr{\\}},\end{split}$

and

$\displaystyle\|T\bar{B}(t)T^{-1}\|_{G}\leq\max_{i}\Bigl{\\{}\sum_{j}\|T_{i}\bar{B}_{ij}(t)T_{j}^{-1}\|_{p}\Bigr{\\}}.$

Condition $C2$ and $C3$ yield

$\max_{i}\Bigl{\\{}\mu_{p}(T_{i}\bar{A}_{ii}(t)T_{i}^{-1})+\sum_{j\neq
i}\|T_{i}\bar{A}_{ij}(t)T_{j}^{-1}\|_{p}\Bigr{\\}}\leq-\bar{\sigma},$

and

$\max_{i}\Bigl{\\{}\sum_{j}\|T_{i}\bar{B}_{ij}(t)T_{j}^{-1}\|_{p}\Bigr{\\}}\leq\underline{\sigma},$

for some $0\leq\underline{\sigma}<\bar{\sigma}<+\infty$. This implies that

$\mu_{G}(TA(t)T^{-1})+\|TB(t)T^{-1}\|_{G}\leq-\bar{\sigma}+\underline{\sigma}:=-\sigma,$

and Lemma 2 then yields

$\displaystyle|\tilde{e}(t)|_{G}\leq$
$\displaystyle\sup_{t_{0}-\tau_{\max}\leq s\leq
t_{0}}|\tilde{e}(s)|_{G}e^{-\hat{\lambda}(t-t_{0})}+\frac{\|T\|_{G}}{\bar{\sigma}-\underline{\sigma}}\max_{i}\|\tilde{w}_{i}(\cdot)\|_{\mathcal{L}_{\infty}^{p}},$

with $\hat{\lambda}$ defined as in the statement of the proposition. Since
$\tilde{e}(t)=Te(t)$ we get $|e(t)|_{G}\leq\|T^{-1}\|_{G}|\tilde{e}(t)|_{G}$
and $|\tilde{e}(t)|_{G}\leq\|T\|_{G}|e(t)|_{G}$. We also notice that the
definition of $\tilde{w}_{i}(t)$ implies that
$\|\tilde{w}_{i}(\cdot)\|_{\mathcal{L}_{\infty}^{p}}=\|w_{i}(\cdot)\|_{\mathcal{L}_{\infty}^{p}}$.
Hence

$\displaystyle|e(t)|_{G}\leq\|T^{-1}\|_{G}\|T\|_{G}\biggl{(}\sup_{t_{0}-\tau_{\max}\leq
s\leq t_{0}}|e(s)|_{G}e^{-\hat{\lambda}(t-t_{0})}$
$\displaystyle\phantom{=}+\frac{1}{\bar{\sigma}-\underline{\sigma}}\max_{i}\|w_{i}(\cdot)\|_{\mathcal{L}_{\infty}^{p}}\biggr{)}.$

We note that
$|e_{i}(t)|_{p}=|[x_{i}^{\mathsf{T}}(t)-x_{i}^{\ast\mathsf{T}}(t),\zeta_{i,1}^{\mathsf{T}}(t),\ldots,\zeta_{i,m}^{\mathsf{T}}(t)|_{p}\allowbreak\geq|[x_{i}^{\mathsf{T}}(t)-x_{i}^{\ast\mathsf{T}}(t),0_{1\times
n},\ldots,0_{1\times n}]|_{p}=|x_{i}(t)-x_{i}^{\ast}(t)|_{p},$ and
$|e_{i}(t_{0})|_{p}=|[x_{i}^{\mathsf{T}}(t_{0})-x_{i}^{\ast\mathsf{T}}(t_{0}),\zeta_{i,1}^{\mathsf{T}}(t_{0}),\ldots,\zeta_{i,m}^{\mathsf{T}}(t_{0})]|_{p}\leq|x_{i}(t_{0})-x_{i}^{\ast}(t_{0})|_{p}+|\zeta_{i,1}(t_{0})|_{p}+\cdots+|\zeta_{i,m}(t_{0})|_{p}$.
Hence, $|x_{i}(t)-x_{i}^{\ast}(t)|_{p}\leq|e_{i}(t)|_{p}$ and
$|e_{i}(t_{0})|_{p}\leq|x_{i}(t_{0})-x_{i}^{\ast}(t_{0})|_{p}+\sum_{k=1}^{m}|\zeta_{i,k}(t_{0})|_{p}$.
We then finally obtain the upper bound of the state deviation

$\displaystyle\begin{split}&\max_{i}|x_{i}(t)-x_{i}^{\ast}(t)|_{p}\leq\frac{\kappa_{G}(T)}{\bar{\sigma}-\underline{\sigma}}\max_{i}\|w_{i}(\cdot)\|_{\mathcal{L}_{\infty}^{p}}\\\
&+\kappa_{G}(T)e^{-\hat{\lambda}(t-t_{0})}\max_{i}\sup_{t_{0}-\tau_{\max}\leq
s\leq t_{0}}|x_{i}(s)-x_{i}^{\ast}(s)|_{p}\\\
&+\kappa_{G}(T)e^{-\hat{\lambda}(t-t_{0})}\max_{i}\sup_{t_{0}-\tau_{\max}\leq
s\leq t_{0}}\sum_{k=1}^{m}|r_{i,k}(s)\\\
&+\sum_{b=0}^{m-k}\frac{(m-1-b)!}{(m-k-b)!}\cdot\bar{d}_{i,m-1-b}\cdot
s^{m-k-b}|_{p},\forall N.\end{split}$

##  Appendix B Recasting the fulfillment of $C2$ and $C3$ as an optimisation
problem

$\begin{split}\bar{B}_{ii}(t)=|\mathcal{N}_{i}|\left[\begin{matrix}\bar{g}_{0}I_{2}&0_{2}&0_{2}\\\
\bar{g}_{1}I_{2}&0_{2}&0_{2}\\\
\bar{g}_{2}I_{2}&0_{2}&0_{2}\end{matrix}\right]\frac{\partial\tanh(\eta_{j}-\eta_{i}-\delta_{ji}^{\ast})}{\partial\eta_{i}},\bar{B}_{ij}(t)=\left[\begin{matrix}\bar{g}_{0}I_{2}&0_{2}&0_{2}\\\
\bar{g}_{1}I_{2}&0_{2}&0_{2}\\\
\bar{g}_{2}I_{2}&0_{2}&0_{2}\end{matrix}\right]\frac{\partial\tanh(\eta_{j}-\eta_{i}-\delta_{ji}^{\ast})}{\partial\eta_{j}}.\end{split}$
(21)

The optimisation problem in (15) was obtained by noticing that $C2$ and $C3$
can be fulfilled by solving the following optimisation problem:

$\displaystyle\begin{split}&\min_{\bar{\xi}}\ \mathcal{J}\\\ &s.t.\quad
k_{0}\geq 0,k_{1}\geq 0,k_{2}\geq 0,k_{0}^{(\tau)}\geq 0,k_{1}^{(\tau)}\geq
0,k_{2}^{(\tau)}\geq 0,\\\
&\phantom{s.t.\quad}k^{\psi}>0,k_{0}+k_{0}^{(\tau)}>0,k_{1}+k_{1}^{(\tau)}>0,k_{2}+k_{2}^{(\tau)}>0,\\\
&\phantom{s.t.\quad}\bar{\sigma}>0,\underline{\sigma}\geq
0,\bar{\sigma}-\underline{\sigma}>0,\mu_{2}(\bar{T}\bar{A}_{ii}\bar{T}^{-1})\leq-\bar{\sigma},\\\
&\phantom{s.t.\quad}\sum_{j\in\mathcal{N}_{i}}\|\bar{T}\bar{B}_{ij}(t)\bar{T}^{-1}\|_{2}+\|\bar{T}\bar{B}_{ii}(t)\bar{T}^{-1}\|_{2}\leq\underline{\sigma},\end{split}$
(22)

where the decision variables are
$\bar{\xi}:=[k_{0},k_{1},k_{2},k_{0}^{(\tau)},k_{1}^{(\tau)},\allowbreak
k_{2}^{(\tau)},k^{\psi},\underline{\sigma},\bar{\sigma},\alpha_{1},\alpha_{2}]$
and the cost is defined as in Section 5. The matrices
$\bar{B}_{ii},\bar{B}_{ij}$ are given in (21) and $\bar{A}_{ii}$ is given in
(14), in accordance with Proposition 1 while the transformation matrix
$\bar{T}$ is also given in (14). In order to find the control gains, we
propose to solve the optimisation problem for fixed $\alpha_{1},\alpha_{2}$.
Further, in order to obtain a suitable formulation for the optimisation, we
recast the constraints in (22) as LMIs as follows. First, by definition,
$\mu_{2}(\bar{T}\bar{A}_{ii}\bar{T}^{-1})\leq-\bar{\sigma}$ is equivalent to
$[\bar{T}\bar{A}_{ii}\bar{T}^{-1}]_{s}\preceq-\bar{\sigma}I_{6}$. Moreover,
the constraint
$\sum_{j\in\mathcal{N}_{i}}\|\bar{T}\bar{B}_{ij}(t)\bar{T}^{-1}\|_{2}+\|\bar{T}\bar{B}_{ii}(t)\bar{T}^{-1}\|_{2}\leq\underline{\sigma}$
is satisfied if we impose that
$\|\bar{T}\bar{B}_{ii}(t)\bar{T}^{-1}\|_{2}\leq\frac{\underline{\sigma}}{2}$
and, simultaneously,
$\|\bar{T}\bar{B}_{ij}(t)\bar{T}^{-1}\|_{2}\leq\frac{\underline{\sigma}}{2\bar{N}},\forall
j\in\mathcal{N}_{i}$. In turn, since
$-1\leq\frac{\partial\tanh(\eta_{j}-\eta_{i}-\delta_{ji}^{\ast})}{\partial\eta_{i}}\leq
0$ and $|\mathcal{N}_{i}|\leq\bar{N}$ we have (by means of the absolutely
homogeneous property for matrix norms) that
$\|\bar{T}\bar{B}_{ii}(t)\bar{T}^{-1}\|_{2}\leq\|\bar{T}\tilde{B}_{ii}\bar{T}^{-1}\|_{2}$
with $\tilde{B}_{ii}$ defined in (14). Analogously, we have that
$\|\bar{T}\bar{B}_{ij}(t)\bar{T}^{-1}\|_{2}\leq\|\bar{T}\tilde{B}_{ij}\bar{T}^{-1}\|_{2}$,
with $\tilde{B}_{ij}$ defined in (14). Hence, the constraints on the norm in
(22) are satisfied if
$\|\bar{T}\tilde{B}_{ii}\bar{T}^{-1}\|_{2}\leq\frac{\underline{\sigma}}{2}$
and
$\|\bar{T}\tilde{B}_{ij}\bar{T}^{-1}\|_{2}\leq\frac{\underline{\sigma}}{2\bar{N}},\forall
j\in\mathcal{N}_{i}$ which, following (Boyd and Vandenberghe, 2004, Example
4.6.3), can be written as
$\left(\frac{\underline{\sigma}}{2\bar{N}}\right)^{2}I_{6}-(\bar{T}\tilde{B}_{ij}\bar{T}^{-1})^{\mathsf{T}}(\bar{T}\tilde{B}_{ij}\bar{T}^{-1})\succeq
0$ and
$\left(\frac{\underline{\sigma}}{2}\right)^{2}I_{6}-(\bar{T}\tilde{B}_{ii}\bar{T}^{-1})^{\mathsf{T}}(\bar{T}\tilde{B}_{ii}\bar{T}^{-1})\succeq
0$. Now, by means of Schur complement (Horn and Johnson, 1985, Theorem 7.7.7),
this pair of inequalities is equivalent to

$\left[\begin{matrix}\frac{\underline{\sigma}}{2\bar{N}}I_{6}&(\bar{T}\tilde{B}_{ij}\bar{T}^{-1})^{\mathsf{T}}\\\
\bar{T}\tilde{B}_{ij}\bar{T}^{-1}&\frac{\underline{\sigma}}{2\bar{N}}I_{6}\end{matrix}\right]\succeq
0,$
$\left[\begin{matrix}\frac{\underline{\sigma}}{2}I_{6}&(\bar{T}\tilde{B}_{ii}\bar{T}^{-1})^{\mathsf{T}}\\\
\bar{T}\tilde{B}_{ii}\bar{T}^{-1}&\frac{\underline{\sigma}}{2}I_{6}\end{matrix}\right]\succeq
0.$

Since the cost is quadratic and the constraints are LMIs, the optimisation
problem in (15) is convex for fixed $\alpha$’s.

## References

  * Abolfazli et al. (2021) Abolfazli, E., Besselink, B., Charalambous, T., 2021. On time headway selection in platoons under the MPF topology in the presence of communication delays. IEEE Transactions on Intelligent Transportation Systems, 1–14.
  * Aminzare and Sontag (2014) Aminzare, Z., Sontag, E. D., 2014. Contraction methods for nonlinear systems: A brief introduction and some open problems. In: 53rd IEEE Conference on Decision and Control. IEEE, pp. 3835–3847.
  * Andreasson et al. (2014) Andreasson, M., Dimarogonas, D. V., Sandberg, H., Johansson, K. H., 2014. Distributed control of networked dynamical systems: Static feedback, integral action and consensus. IEEE Transactions on Automatic Control 59 (7), 1750–1764.
  * Besselink and Johansson (2017) Besselink, B., Johansson, K. H., 2017. String stability and a delay-based spacing policy for vehicle platoons subject to disturbances. IEEE Transactions on Automatic Control 62 (9), 4376–4391.
  * Besselink and Knorn (2018) Besselink, B., Knorn, S., 2018. Scalable input-to-state stability for performance analysis of large-scale networks. IEEE Control Systems Letters 2 (3), 507–512.
  * Boyd and Vandenberghe (2004) Boyd, S., Vandenberghe, L., 2004. Convex optimization. Cambridge university press.
  * Burbano Lombana and di Bernardo (2016) Burbano Lombana, D. A., di Bernardo, M., 2016. Multiplex PI control for consensus in networks of heterogeneous linear agents. Automatica 67, 310–320.
  * Centorrino et al. (2022) Centorrino, V., Bullo, F., Russo, G., 2022. Contraction analysis of Hopfield neural networks with Hebbian learning. arXiv preprint arXiv:2204.05382.
  * Davydov et al. (2021) Davydov, A., Proskurnikov, A. V., Bullo, F., 2021. Non-Euclidean contractivity of recurrent neural networks. arXiv preprint arXiv:2110.08298.
  * De Domenico et al. (2013) De Domenico, M., Solé-Ribalta, A., Cozzo, E., Kivelä, M., Moreno, Y., Porter, M. A., Gómez, S., Arenas, A., 2013. Mathematical formulation of multilayer networks. Physical Review X 3 (4), 041022.
  * Desoer and Haneda (1972) Desoer, C., Haneda, H., 1972. The measure of a matrix as a tool to analyze computer algorithms for circuit analysis. IEEE Transactions on Circuit Theory 19 (5), 480–486.
  * di Bernardo et al. (2016) di Bernardo, M., Fiore, D., Russo, G., Scafuti, F., 2016. Convergence, consensus and synchronization of complex networks via contraction theory. Complex Systems and Networks, 313–339.
  * di Bernardo et al. (2015) di Bernardo, M., Salvi, A., Santini, S., 2015. Distributed consensus strategy for platooning of vehicles in the presence of time-varying heterogeneous communication delays. IEEE Transactions on Intelligent Transportation Systems 16 (1), 102–112.
  * Dörfler and Bullo (2014) Dörfler, F., Bullo, F., 2014. Synchronization in complex networks of phase oscillators: A survey. Automatica 50 (6), 1539–1564.
  * Fan et al. (2021) Fan, H., Ramirez, A., Sipahi, R., 2021. Consensus stability of a large scale robotic network under input and transmission delays. IEEE Transactions on Control of Network Systems, 1–1.
  * Feng et al. (2019) Feng, S., Zhang, Y., Li, S. E., Cao, Z., Liu, H. X., Li, L., 2019. String stability for vehicular platoon control: Definitions and analysis methods. Annual Reviews in Control 47, 81–97.
  * Freeman et al. (2006) Freeman, R. A., Yang, P., Lynch, K. M., 2006. Stability and convergence properties of dynamic average consensus estimators. In: 45th IEEE Conference on Decision and Control. IEEE, pp. 338–343.
  * Gao et al. (2008) Gao, H., Chen, T., Lam, J., 2008. A new delay system approach to network-based control. Automatica 44 (1), 39–52.
  * Giaccagli et al. (2021) Giaccagli, M., Astolfi, D., Andrieu, V., Marconi, L., 2021. Sufficient conditions for global integral action via incremental forwarding for input-affine nonlinear systems. IEEE Transactions on Automatic Control, 1–1.
  * Gomez et al. (2013) Gomez, S., Diaz-Guilera, A., Gomez-Gardenes, J., Perez-Vicente, C. J., Moreno, Y., Arenas, A., 2013. Diffusion dynamics on multiplex networks. Physical review letters 110 (2), 028701.
  * Horn and Johnson (1985) Horn, R. A., Johnson, C. R., 1985. Matrix Analysis. Cambridge University Press.
  * Jafarpour et al. (2021) Jafarpour, S., Davydov, A., Proskurnikov, A., Bullo, F., 2021. Robust implicit networks via non-Euclidean contractions. Advances in Neural Information Processing Systems 34.
  * Kim et al. (2010) Kim, K.-S., Rew, K.-H., Kim, S., 2010. Disturbance observer for estimating higher order disturbances in time series expansion. IEEE Transactions on Automatic Control 55 (8), 1905–1911.
  * Knorn and Besselink (2020) Knorn, S., Besselink, B., 2020. Scalable robustness of interconnected systems subject to structural changes. IFAC-PapersOnLine 53 (2), 3373–3378.
  * Knorn et al. (2014a) Knorn, S., Donaire, A., Agüero, J. C., Middleton, R. H., 2014a. Passivity-based control for multi-vehicle systems subject to string constraints. Automatica 50 (12), 3224–3230.
  * Knorn et al. (2014b) Knorn, S., Donaire, A., Agüero, J. C., Middleton, R. H., 2014b. Passivity-based control for multi-vehicle systems subject to string constraints. Automatica 50 (12), 3224–3230.
  * Lawton et al. (2003) Lawton, J., Beard, R., Young, B., 2003. A decentralized approach to formation maneuvers. IEEE Transactions on Robotics and Automation 19 (6), 933–941.
  * Li et al. (2019) Li, D., Ge, S. S., He, W., Ma, G., Xie, L., 2019. Multilayer formation control of multi-agent systems. Automatica 109, 108558.
  * Li et al. (2013) Li, W., Chen, Z., Liu, Z., 2013. Leader-following formation control for second-order multiagent systems with time-varying delay and nonlinear dynamics. Nonlinear Dynamics 72 (4), 803–812.
  * Liu et al. (2019) Liu, H., Tian, Y., Lewis, F. L., Wan, Y., Valavanis, K. P., 2019. Robust formation tracking control for multiple quadrotors under aggressive maneuvers. Automatica 105, 179–185.
  * Liu et al. (2018) Liu, X., Ge, S. S., Goh, C.-H., 2018. Vision-based leader–follower formation control of multiagents with visibility constraints. IEEE Transactions on Control Systems Technology 27 (3), 1326–1333.
  * Lohmiller and Slotine (1998) Lohmiller, W., Slotine, J.-J., 1998. On contraction analysis for non-linear systems. Automatica 34 (6), 683–696.
  * Mirabilio et al. (2021) Mirabilio, M., Iovine, A., De Santis, E., Di Benedetto, M. D., Pola, G., 2021. String stability of a vehicular platoon with the use of macroscopic information. IEEE Transactions on Intelligent Transportation Systems 22 (9), 5861–5873.
  * Monteil and Russo (2017) Monteil, J., Russo, G., 2017. On the design of nonlinear distributed control protocols for platooning systems. IEEE Control Systems Letters 1 (1), 140–145.
  * Monteil et al. (2019) Monteil, J., Russo, G., Shorten, R., 2019. On $\mathcal{L}_{\infty}$ string stability of nonlinear bidirectional asymmetric heterogeneous platoon systems. Automatica 105, 198–205.
  * Mucha et al. (2010) Mucha, P. J., Richardson, T., Macon, K., Porter, M. A., Onnela, J.-P., 2010. Community structure in time-dependent, multiscale, and multiplex networks. Science 328 (5980), 876–878.
  * Ofir et al. (2022) Ofir, R., Bullo, F., Margaliot, M., 2022. Minimum effort decentralized control design for contracting network systems. arXiv preprint arXiv:2203.10392.
  * Olfati-Saber and Murray (2004) Olfati-Saber, R., Murray, R., 2004. Consensus problems in networks of agents with switching topology and time-delays. IEEE Transactions on Automatic Control 49 (9), 1520–1533.
  * Panagou and Kyriakopoulos (2011) Panagou, D., Kyriakopoulos, K. J., 2011. Switching control approach for the robust practical stabilization of a unicycle-like marine vehicle under non-vanishing perturbations. In: 2011 IEEE International Conference on Robotics and Automation. pp. 1525–1530.
  * Pant et al. (2002) Pant, A., Seiler, P., Hedrick, K., 2002. Mesh stability of look-ahead interconnected systems. IEEE Transactions on Automatic Control 3, 403 – 407.
  * Park et al. (2012) Park, G., Joo, Y., Shim, H., Back, J., 2012. Rejection of polynomial-in-time disturbances via disturbance observer with guaranteed robust stability. In: 51st IEEE Conference on Decision and Control. pp. 949–954.
  * Ploeg et al. (2014) Ploeg, J., van de Wouw, N., Nijmeijer, H., 2014. $\mathcal{L}_{p}$ string stability of cascaded systems: Application to vehicle platooning. IEEE Transactions on Control Systems Technology 22 (2), 786–793.
  * Qian and Del Vecchio (2018) Qian, Y., Del Vecchio, D., 2018. Realizing integral control in living cells: how to overcome leaky integration due to dilution? Journal of The Royal Society Interface 15 (139), 20170902.
  * Revay and Manchester (2020) Revay, M., Manchester, I., 2020. Contracting implicit recurrent neural networks: Stable models with improved trainability. In: 2nd Conference on Learning for Dynamics and Control. Vol. 120 of Proceedings of Machine Learning Research. PMLR, pp. 393–403.
  * Robotarium (2022) Robotarium, 2022. Robotarium Matlab Simulator.
URL https://tinyurl.com/3rajpnep

  * Russo et al. (2010a) Russo, G., di Bernardo, M., Sontag, E. D., 2010a. Global entrainment of transcriptional systems to periodic inputs. PLoS computational biology 6 (4), e1000739.
  * Russo et al. (2010b) Russo, G., di Bernardo, M., Sontag, E. D., 2010b. Stability of networked systems: A multi-scale approach using contraction. In: 49th IEEE Conference on Decision and Control. pp. 6559–6564.
  * Russo and Wirth (2022) Russo, G., Wirth, F., 2022. Matrix measures, stability and contraction theory for dynamical systems on time scales. Discrete and Continuous Dynamical Systems - B 27 (6), 3345–3374.
  * Seiler et al. (1999) Seiler, P., Pant, A., Hedrick, J., 1999. Preliminary investigation of mesh stability for linear systems. In: ASME International Mechanical Engineering Congress and Exposition. Vol. 16349. American Society of Mechanical Engineers, pp. 359–364.
  * Seyboth and Allgöwer (2015) Seyboth, G. S., Allgöwer, F., 2015. Output synchronization of linear multi-agent systems under constant disturbances via distributed integral action. In: 2015 American Control Conference. IEEE, pp. 62–67.
  * Shiromoto et al. (2019) Shiromoto, H. S., Revay, M., Manchester, I. R., 2019. Distributed nonlinear control design using separable control contraction metrics. IEEE Transactions on Control of Network Systems 6 (4), 1281–1290.
  * Siljak (2011) Siljak, D. D., 2011. Decentralized control of complex systems. Courier Corporation.
  * Silva et al. (2021) Silva, G. F., Donaire, A., McFadyen, A., Ford, J. J., 2021. String stable integral control design for vehicle platoons with disturbances. Automatica 127, 109542.
  * Silva et al. (2022) Silva, G. F., Donaire, A., Seron, M. M., McFadyen, A., Ford, J., 2022. String stability in microgrids using frequency controlled inverter chains. IEEE Control Systems Letters 6, 1484–1489.
  * Simpson-Porco (2021) Simpson-Porco, J. W., 2021. Analysis and synthesis of low-gain integral controllers for nonlinear systems. IEEE Transactions on Automatic Control 66 (9), 4148–4159.
  * Sridhar and Govindarasu (2014) Sridhar, S., Govindarasu, M., 2014. Model-based attack detection and mitigation for automatic generation control. IEEE Transactions on Smart Grid 5 (2), 580–591.
  * Stüdli et al. (2017) Stüdli, S., Seron, M., Middleton, R., 2017. From vehicular platoons to general networked systems: String stability and related concepts. Annual Reviews in Control 44, 157–172.
  * Swaroop and Hedrick (1996) Swaroop, D., Hedrick, J., 1996. String stability of interconnected systems. IEEE Transactions on Automatic Control 41 (3), 349–357.
  * Tanner et al. (2002) Tanner, H., Pappas, G., Kumar, V., 2002. Input-to-state stability on formation graphs. In: 41st IEEE Conference on Decision and Control, 2002. Vol. 3. pp. 2439–2444 vol.3.
  * Tanner et al. (2004) Tanner, H., Pappas, G., Kumar, V., 2004. Leader-to-formation stability. IEEE Transactions on Robotics and Automation 20 (3), 443–455.
  * Tran et al. (2020) Tran, D., Yucelen, T., Pasiliao, E. L., 2020. Formation control with multiplex information networks. IEEE Transactions on Control Systems Technology 28 (2), 462–476.
  * Tsukamoto et al. (2021) Tsukamoto, H., Chung, S.-J., Slotine, J.-J. E., 2021. Contraction theory for nonlinear stability analysis and learning-based control: A tutorial overview. Annual Reviews in Control 52, 135–169.
  * Wang and Slotine (2006) Wang, W., Slotine, J.-J., 2006. Contraction analysis of time-delayed communications and group cooperation. IEEE Transactions on Automatic Control 51 (4), 712–717.
  * Wang et al. (2013) Wang, X., Saberi, A., Stoorvogel, A. A., Grip, H. F., Yang, T., 2013. Consensus in the network with uniform constant communication delay. Automatica 49 (8), 2461–2467.
  * Wen et al. (2008) Wen, L., Yu, Y., Wang, W., 2008. Generalized halanay inequalities for dissipativity of volterra functional differential equations. Journal of Mathematical Analysis and Applications 347 (1), 169–178.
  * Wilson et al. (2020) Wilson, S., Glotfelter, P., Wang, L., Mayya, S., Notomista, G., Mote, M., Egerstedt, M., 2020. The robotarium: Globally impactful opportunities, challenges, and lessons learned in remote-access, distributed control of multirobot systems. IEEE Control Systems Magazine 40 (1), 26–44.
  * Xie and Russo (2022) Xie, S., Russo, G., 2022. On the design of scalable networks rejecting first order disturbances. In: NecSys 2022. See https://arxiv.org/abs/2202.07638 for an extended version of the conference paper with preliminary proofs.
  * Xie et al. (2021) Xie, S., Russo, G., Middleton, R. H., 2021. Scalability in nonlinear network systems affected by delays and disturbances. IEEE Transactions on Control of Network Systems 8 (3), 1128–1138.
  * Xie and Liu (2013) Xie, X.-J., Liu, L., 2013. A homogeneous domination approach to state feedback of stochastic high-order nonlinear systems with time-varying delay. IEEE Transactions on Automatic Control 58 (2), 494–499.