# Incremental stabilization and multi-agent synchronization of discrete-time nonlinear systems

Samuele Zoboli, Daniele Astolfi, Vincent Andrieu, Giacomo Casadei, Luca Zaccarian

## ► To cite this version:

Samuele Zoboli, Daniele Astolfi, Vincent Andrieu, Giacomo Casadei, Luca Zaccarian. Incremental stabilization and multi-agent synchronization of discrete-time nonlinear systems. 2024. hal-04444190v3

HAL Id: hal-04444190

https://hal.science/hal-04444190v3

Preprint submitted on 14 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Incremental stabilization and multi-agent synchronization of discrete-time nonlinear systems

Samuele Zoboli<sup>a</sup>, Daniele Astolfi<sup>c</sup>, Vincent Andrieu<sup>c</sup>, Giacomo Casadei<sup>d</sup>, Luca Zaccarian<sup>b</sup>

## Abstract

In this paper, we propose a novel distributed state-feedback design for robust synchronization of networks of identical discrete-time nonlinear agents under generic time-invariant communication graphs. We focus on the class of almost differentiable (possibly time-varying) dynamics that are linear in the input. By generalizing results on synchronization of linear agents, we build strong links between the solution to the synchronization problem in the linear and nonlinear framework. This is also enabled by the introduction of new results on design of incrementally stabilizing controllers based on contraction analysis. Finally, we propose numerically tractable sufficient conditions for the synchronization of networks of non-smooth Lur'e systems.

**Keywords:** Discrete-time, nonlinear systems, synchronization, multi-agent systems, contraction, incremental stability.

## 1. Introduction

Multi-agent systems control has attracted a lot of attention from our community. Many modern control problems can be formalized as a network of interacting agents aiming at achieving some sort of agreement [20]. In this paper, we focus on the problem of state synchronization of a network of homogeneous systems (namely, all the agents are identical), described by a time-varying discrete-time nonlinear system that is linear in the control input. The problem of synchronization presents mature results in the linear framework, especially for continuous-time agent dynamics [21, 27]. For discrete-time systems, major contributions can be found in [6, 26, 32], and a small gain theorem for interconnected systems has recently appeared in [28]. However, in this setting, further investigations are still required to understand the relationship between the structure of a network (i.e., the eigenvalues and eigenvectors associated with its representation) and the possibility of finding a suitable synchronizing control. As a matter of fact, the structure of the communication graph has a significant impact on discrete-time networks compared to continuous-time ones. Only recently a novel implicit formulation allowed overcoming some of these structural limitations, proposing a result comparable to the

continuous-time framework [3]. Avoiding implicit computations, graph normalization approaches like [13, 32] have proven to be effective, yet they may not always be feasible. Additionally, existing solutions for nonlinear dynamics are limited to specific agent structures, such as Lur'e system forms [30] or linear systems with saturated inputs [4]. Moreover, they commonly employ observer design [14, 23] or data-based optimization techniques [10].

In this work, we propose solutions to the synchronization problem based on discrete-time contraction analysis [16, 18, 19, 35] and incremental stability [1, 9, 12, 30, 33, 34]. There are two main motivations for this choice. Firstly, contraction analysis allows us to study nonlinear systems via linear systems-like arguments. Hence, we can take inspiration from the well-established linear systems literature [32] and provide a link between the two scenarios. Secondly, incremental stability easily translates to synchronization of homogeneous networks. As a matter of fact, trajectories originated from incrementally stable dynamics “forget” their initial conditions, while offering strong robustness properties. In particular, in a network where the agents are homogeneous, each one can be considered as a singular trajectory of the same system starting from different initial conditions. As distances between trajectories of incrementally stable systems asymptotically decrease to zero, by designing a distributed controller making the network dynamics incrementally stable we indirectly obtain robust state synchronization<fnref n="1" />.

The three main contributions of this paper can be sum-

<footnote><sup>a</sup>LAAS-CNRS, Université de Toulouse, UPS, Toulouse, France</footnote>

<footnote><sup>b</sup>LAAS-CNRS, Université de Toulouse, CNRS, Toulouse, France</footnote>

<footnote><sup>c</sup>Université Claude Bernard Lyon 1, CNRS, LAGEPP UMR 5007, 43 boulevard du 11 novembre 1918, F-69100, Villeurbanne, France</footnote>

<footnote><sup>d</sup>Université Grenoble Alpes, CNRS, Inria, Grenoble INP, GIPSA-lab, F-38000 Grenoble, France.</footnote>

<footnote>*Research supported by ANR via grants DELICIO, number ANR-19-CE23-0006, and OLYMPIA, number ANR-23-CE48-0006.
Email addresses: samuele.zoboli@laas.fr (Samuele Zoboli), daniele.astolfi@univ-lyon1.fr (Daniele Astolfi), vincent.andrieu@univ-lyon1.fr (Vincent Andrieu), giacomo.casadei@univ-grenoble-alpes.fr (Giacomo Casadei), zaccarian@laas.fr (Luca Zaccarian)</footnote>

<footnote n="1">¹It is worth noting that in the context of synchronization it is</footnote>
marized as follows:

(i) We present new results for the robust synchronization of identical discrete-time linear systems based on the Modified Algebraic Riccati Inequality (MARI) [4, 31, 36, 42]. Unlike approaches that rely on normalized Laplacians or agent-specific gains [13, 26, 32], we assume all controllers are identical. Our method, which includes the normalized Laplacian case as a sub-case, yields bounds linking the graph's connectivity with the simultaneous stabilizability of the agents.

(ii) We revisit the problem of incremental stabilization of a nonlinear system and we provide a new closed-form state-feedback controller guaranteeing uniform exponential incremental stability properties for the closed-loop system. The proposed feedback is a direct generalization of a linear LQR design. The conditions are based on contraction analysis and apply to dynamics that are time-varying and not differentiable everywhere, generalizing existing results [7, 9, 19, 35]. Note that while other works address a similar problem using control contraction metrics [37, 38], their solutions are non-constructive, as they rely on solving nested optimization problems at each time step.

(iii) We exploit the new results on contraction analysis and incremental stability to provide new controller designs for robust synchronization of identical discrete-time nonlinear agents. This is done under the assumptions of almost differentiable dynamics, linear input vector field, and generic connected communication graphs (i.e., possibly directed and weighted). As in the linear case, we relate the connectivity properties of the graph to the simultaneous stabilizability properties of the agents.

The paper is organized as follows. Section 2 introduces some preliminary concepts of graph theory and presents the problem of exponential synchronization. Section 3 presents the results for the linear framework. Section 4 extends the solution to nonlinear agents. Section 5 provides an LMI-based method for controller design. Numerical examples are illustrated in Section 6. Finally, Section 7 concludes the paper while the Appendix contains the proofs of our results.

**Notation.** We denote by $\mathbb{R}$ the set of real numbers, by $\mathbb{R}_{\ge 0}$ the set of non-negative real numbers, by $\mathbb{N}$ the set of non-negative integers and by $\mathbb{C}$ the set of complex numbers. Given a complex number $z \in \mathbb{C}$, we use $\Re(z)$ and $\Im(z)$ to identify its real and imaginary parts respectively. $|\cdot|$ denotes the standard Euclidean norm, $|x|_{\mathcal{M}}$ is a distance function between any point $x \in \mathbb{R}^n$ and a closed set $\mathcal{M} \subset \mathbb{R}^n$, namely, $|x|_{\mathcal{M}} := \inf_{z \in \mathcal{M}} |x-z|$. Given two symmetric matrices $A, B \in \mathbb{R}^{n \times n}$, we say $A \preceq B$ if the matrix $A - B$ is negative semidefinite. Similarly, we say $A \succ B$ if $A - B$ is positive definite. Moreover, we denote by $\mathbb{S}_{\ge 0}^n$ (resp. $\mathbb{S}_{\le 0}^n$) the set of symmetric positive definite (resp. semi-definite) real matrices of dimension $n$. Given a set of

vectors $x_1, \dots, x_N$ with $N \in \mathbb{N}$ and arbitrary dimensions, we define $\text{col}(x_1, \dots, x_N) = [x_1^\top \cdots x_N^\top]^\top$. Given a square matrix $A \in \mathbb{R}^{n \times n}$, we use $\text{spec } A$ to represent its spectrum. Given a set of square matrices $A_1, \dots, A_N$ with $N \in \mathbb{N}$ we identify by $\text{diag}(A_1, \dots, A_N)$ the block diagonal matrix whose diagonal blocks are $A_1, \dots, A_N$ respectively. Given a positive integer $n \in \mathbb{N}$, we identify by $I_n$ the identity matrix of dimension $n$. We use **1**, resp. **0**, to identify a column vector of 1s, resp. 0s, of appropriate dimension. The symbol $\otimes$ denotes the Kronecker product, which satisfies $(A \otimes B)^\top = A^\top \otimes B^\top$, $(A \otimes B)(C \otimes D) = (AC \otimes BD)$ and for $A, B, C$ symmetric, $C \succeq 0$, $A \preceq B \implies (A \otimes C) \preceq (B \otimes C)$ and $(C \otimes A) \preceq (C \otimes B)$. For a square matrix $A \in \mathbb{R}^{n \times n}$, we define $\text{Sym}(A) := \frac{1}{2}(A^\top + A)$. Finally, we denote with $\text{co}$ the convex hull.

## 2. Problem statement

In this section, we first present some preliminary concepts related to graph theory and multi-agent networks. Then, we formulate the problem of exponential synchronization.

### 2.1. Highlights on graph theory

In multi-agent systems, a communication graph is typically described by a triplet $\mathcal{G} = \{\mathcal{V}, \mathcal{E}, \mathcal{A}\}$ where $\mathcal{V} = \{v_1, v_2, \dots, v_N\}$ is a set of $N \subset \mathbb{N}$ vertices (or nodes), $\mathcal{E} \subset \mathcal{V} \times \mathcal{V}$ is the set of edges $\epsilon_{jh}$ modeling the interconnection between such nodes, and $\mathcal{A} \in \mathbb{R}^{N \times N}$ is the adjacency matrix, whose entries $a_{jh} \ge 0$ weight the flow of information from vertex $j$ to vertex $h$ and $a_{jh} > 0 \iff (j, h) \in \mathcal{E}$. We denote by $L \in \mathbb{R}^{N \times N}$ the Laplacian matrix of the graph, defined as

$$ \ell_{jh} = -a_{jh} \quad \text{for } j \neq h, \qquad \ell_{jh} = \sum_{i=1}^{N} a_{ji} \quad \text{for } j = h, $$

where $\ell_{jh}$ is the $(j, h)$-th entry of $L$. We denote with $\mathcal{N}_i$ the set of in-neighbors of node $i$, i.e. the set $\mathcal{N}_i := \{j \in \{1, \dots, N\} | \epsilon_{ij} \in \mathcal{E}\}$. In this paper, we consider directed, weighted graphs. As such, the Laplacian matrix is not assumed to be diagonalizable and admits complex eigenvalues. We also suppose the graph contains a directed spanning tree. Without loss of generality, we assume such a directed spanning tree to be rooted at node 1. Consequently, we identify the Laplacian of the network as

$$ L = \begin{pmatrix} L_{11} & L_{12} \\ L_{21} & L_{22} \end{pmatrix} \qquad (1) $$

where $L_{11}$ is a scalar, $L_{12}$ is a $N-1$ row vector, $L_{21}$ is a $N-1$ column vector and $L_{22}$ is a $(N-1) \times (N-1)$ matrix. Based on these assumptions, we recall the following result from the literature, see e.g., [15, 25].

**Lemma 1 ([25])** The Laplacian matrix $L$ has a simple eigenvalue 0 and all the other eigenvalues have positive real parts if and only if the directed graph has a directed spanning tree.

<footnote>more natural to invoke the notion of incremental stability than the one of convergent system, e.g., [16, 35]. Indeed, state synchronization may be achieved on an unbounded trajectory.</footnote>
## 2.2. The problem of multi-agent synchronization

Consider a network of discrete-time agents, where the dynamics of each node is described by a nonlinear, possibly time-varying, difference equation of the form

$$x_i^+ = f(x_i, t) + Bu_i + w_i, \quad i = 1, \dots, N, \quad (2)$$

where $f: \mathbb{R}^n \times \mathbb{N} \to \mathbb{R}^n$, $B \in \mathbb{R}^{n \times m}$ is full column rank, $x_i \in \mathbb{R}^n$, $u_i \in \mathbb{R}^m$ and $w_i \in \mathbb{R}^n$ represent the state, the control input and an unmeasured disturbance of node $i$ at timestep $t \in \mathbb{N}$, respectively, and $x_i^+ \in \mathbb{R}^n$ represents the state of node $i$ at timestep $t+1$. Note that $w_i \in \mathbb{R}^n$ can be any (not necessarily bounded) input sequence and can represent, for instance, an external disturbance affecting the dynamics or some (bounded) model uncertainty. We define the state of the entire network $\mathbf{x} \in \mathbb{R}^{Nn}$ and the entire disturbance $\mathbf{w} \in \mathbb{R}^{Nn}$ as

$$\mathbf{x} := \mathrm{col}(x_1, \dots, x_N), \quad \mathbf{w} := \mathrm{col}(w_1, \dots, w_N). \quad (3)$$

Our objective is to design a distributed state-feedback control law of the form

$$u_i = \sum_{j \in \mathcal{N}_i} a_{ij} (\alpha(x_i, t) - \alpha(x_j, t)) = \sum_{j=1}^{N} \ell_{ij} \alpha(x_j, t) \quad (4)$$

for all $i = 1, \dots, N$, where the function $\alpha: \mathbb{R}^n \times \mathbb{N} \to \mathbb{R}^m$ is our sole design parameter and it stabilizes the dynamics (2) on the synchronization manifold $\mathcal{M}$ defined as

$$\mathcal{M} := \{\mathbf{x} \in \mathbb{R}^{Nn} \mid x_i = x_j, \text{ for all } i, j \in \{1, \dots, N\}\}, \quad (5)$$

where the states of all the agents of the network agree with each other. Note that the $i$-th agent can only use the state $x_j$ of its neighbors $j \in \mathcal{N}_i$ alongside its own local information $x_i$. Moreover, the control action $u_i$ is equal to zero on the synchronization manifold in the absence of disturbances, i.e. when $\mathbf{w} = 0$. In other words, when synchronization is achieved, no correction term is needed for each individual agent. Hence, independent stabilization of all agents on an agreed equilibrium point is not a valid solution in general. We formalize our synchronization problem as follows.

**Problem 1 (Robust network synchronization)** The distributed feedback control law (4) solves the robust synchronization problem for the agents (2) if there exist a function $\alpha: \mathbb{R}^n \times \mathbb{N} \to \mathbb{R}^m$ and real numbers $c \ge 1$, $\rho \in (0, 1)$ and $\gamma \ge 0$ such that, for any initial condition $(\mathbf{x}(t_0), t_0) \in \mathbb{R}^{Nn} \times \mathbb{N}$, and any disturbance sequences $t \mapsto w_i(t)$, $i \in \{1, \dots, N\}$, the solutions to the closed-loop system

$$x_i^+ = f(x_i, t) + B \sum_{j=1}^{N} \ell_{ij} \alpha(x_j, t) + w_i, \quad i = 1, \dots, N.$$

satisfy for all $t \ge t_0$

$$|\mathbf{x}(t)|_{\mathcal{M}} \le c \rho^{t-t_0} |\mathbf{x}(t_0)|_{\mathcal{M}} + \sup_{\substack{t \in [t_0, t] \\ i, j \in [1, N]}} \gamma |w_i(t) - w_j(t)|, \quad (6)$$

where $\mathcal{M}$ is defined in (5).

**Remark 1** Although (4) assumes that each agent has access to its own and its neighbors complete state information, an output-feedback scenario with the inclusion of an observer and reduction-type arguments [8] could be considered due to the robust convergence required by (6). For the sake of clarity and ease in notation, we focus on the full-state information problem.

## 3. Synchronization of Linear Systems

It is well known that the problem of linear synchronization boils down to a problem of simultaneous stabilization, see, e.g., [32, 41]. In this section, we reinterpret these results into our framework and present the discrete synchronization problem to non-diagonalizable Laplacian matrices. We start by recalling necessary and sufficient conditions for state synchronization. Then, we derive novel conditions for multi-agent synchronization under any Laplacian whose eigenvalues belong to a given compact set by exploiting the gain margin properties of Riccati-based designs. While closely related results on linear discrete-time synchronization exist, our formulation provides a novel specific viewpoint that can be easily generalized to the nonlinear setting (see Section 4), and that is of independent interest also in the linear setting.

We start by presenting a general result for network synchronization for linear systems. It is shown that the existence of a common control law for systems associated with each non-zero eigenvalue of the Laplacian is equivalent to solving the synchronization problem. Similar results have appeared multiple times in the literature, e.g., [21, Theorem 13],[13, Lemma 1], [32, Lemma 1] for the normalized Laplacian scenario and [40, Theorem 3.2]. However, we present a reformulated proof that will be instrumental for the analysis of the nonlinear framework.

**Theorem 1** Consider a network of $N \in \mathbb{N}$ agents described by dynamics

$$x_i^+ = Ax_i + Bu_i + w_i \quad (7)$$

The diffusive control law $u_i = K \sum_{j=1}^{N} \ell_{ij} x_j$ solves Problem 1 if and only if the following two conditions hold:

(i) matrix $A$ is Schur-Cohn stable$^2$, or, if $A$ is not Schur-Cohn stable, the (possibly directed, weighted) interconnection graph $G$ has a directed spanning tree;

(ii) the gain $K$ is such that matrix $(A+\lambda_i BK)$ is Schur-Cohn stable for any $\lambda_i \in \mathrm{spec} L \setminus \{0\}$.

The proof of Theorem 1 can be found in Appendix A.1. Note that the results of Theorem 1 are not constructive. To obtain a set of constructive conditions, we first recall the design procedure for pencil matrices stabilization. Then, in Section 3.2, we apply this result to the case of networks.

<footnote>$^2$A matrix is said to be Schur-Cohn stable, or simply Schur stable, if all of its eigenvalues are inside the unitary disk.</footnote>
### 3.1. The pencil matrix for discrete-time linear systems

Consider a discrete-time linear system described by

$$x^+ = Ax + Bu, \quad u = Kx, \quad (8)$$

with $x \in \mathbb{R}^n$, $u \in \mathbb{R}^m$ and, without loss of generality, $B$ is assumed to be full-column rank. The goal is to find a gain matrix $K \in \mathbb{R}^{m \times n}$ such that the complex closed-loop matrix $(A + \lambda BK)$ is Schur-Cohn for some complex numbers $\lambda$. Inspired by [29, Definition 3.13], we formally define this notion as follows.

#### Definition 1 (Complex gain margin for LTI systems)

The matrix $K$ is said to have a complex gain margin with radius $r > 0$ and center $c$ if $A + \lambda BK$ is Schur-Cohn for any $\lambda$ in $\{\lambda \in \mathbb{C} : |\lambda - c| \le r\}$.

To find the complex gain margin of a matrix $K$, we propose a solution based on the discrete-time Modified Algebraic Riccati Inequality (MARI) [4, 31, 36] defined as

$$A^{\top}PA - \sigma A^{\top}PB(R + B^{\top}PB)^{-1}B^{\top}PA \le \rho P, \quad (9)$$

where $R \in \mathbb{S}_{\ge 0}^m$, $P \in \mathbb{S}_{\ge 0}^n$ and generally $\sigma \in (0, 1]$, $\rho \in (0, 1)$. A differential version of the MARI (9) will be the fundamental tool allowing the derivation of our nonlinear results. Therefore, we briefly discuss some of the properties of such an inequality. Note that, since $B$ is assumed to be full column rank, the matrix $R + B^\top PB$ is positive definite and, consequently, invertible. The main difference between the MARI (9) and the more common discrete-time Algebraic Riccati Inequality (DARI)

$$A^{\top}PA - A^{\top}PB(R + B^{\top}PB)^{-1}B^{\top}PA + Q \prec P \quad (10)$$

lies in the presence of the scalar $\sigma$. First, note that if $R \succ 0$, then the positive semi-definite matrix $Q$ can be embedded in the right-hand side of (9) by exploiting $\rho P \prec P - Q$, which holds for a suitable $\rho \in (0, 1)$ as long as $P - Q \succ 0$. Inequality $P - Q \succ 0$ holds when $R \in \mathbb{S}_{\ge 0}^n$ because one can rearrange (10) as $A^\top[P - PB(R + B^\top PB)^{-1}B^\top P]A \prec P - Q$ and applying the Woodbury matrix identity

$$(A + BD^{-1}C)^{-1} = A^{-1} - A^{-1}B(D + CA^{-1}B)^{-1}CA^{-1}, \quad (11)$$

with $A^{-1} = P$, $D = R$, $C = B^\top$, $B = B$, yields to

$$0 \prec A^{\top}(P^{-1} + BR^{-1}B^{\top})^{-1}A \prec P - Q. \quad (12)$$

Therefore, the DARI (10) is a special case of the MARI (9) when $\sigma = 1$ and $R$ is positive definite. As such, the MARI allows for an extra degree of freedom. Its role is to weigh the impact of the input on the solution to the inequality. In other words, the smaller the $\sigma$, the less we can rely on the input to stabilize the system.

We next reformulate [42, Theorem 2], showing that the degree of freedom offered by the MARI (9) allows obtaining sufficient conditions for the existence of a state feedback gain $K$ solving the simultaneous stabilization problem. This result expands on the findings of [4, Theorem 2], by taking into account the interplay between $R$ and $P$ in the MARI (10), thus allowing for larger certified gain margins.

**Proposition 1** Let the pair $(A, B)$ be stabilizable and $R \in \mathbb{S}_{\ge 0}^m$. Let $P \in \mathbb{S}_{\ge 0}^n$ be a solution to the MARI (9) for some $\sigma \in (0, 1]$ and for some $\rho \in (0, 1)$. Then the matrix

$$K = -(R + B^{\top}PB)^{-1}B^{\top}PA, \quad (13)$$

has a complex gain margin center and radius

$$c = 1 + \frac{\lambda_{\min}(R)}{\lambda_{\max}(B^{\top}PB)}, \quad r = \sqrt{c(c - \sigma)}. \quad (14)$$

Moreover, the following bound holds

$$c - \sqrt{c(c - \sigma)} \le 1 - \sqrt{1 - \sigma} \le 1 + \sqrt{1 - \sigma} \le c + \sqrt{c(c - \sigma)} \quad (15)$$

The proof of Proposition 1 is omitted as it can be derived from [42, Theorem 2].

**Remark 2** Continuous-time systems can achieve infinite-gain margins<fnref n="3" /> because their stable region (the left-half plane) is unbounded [11, 29]. In contrast, since the stable region for discrete-time systems (the unit disc) is bounded, an overly large gain can cause instability, resulting in only finite-gain margins [24].

### 3.2. Main result on robust linear synchronization

By relying on the result of Proposition 1, we now provide a set of constructive conditions for Theorem 1. These conditions will be generalized to the nonlinear systems framework in Section 4 and thus provide an essential intuition on the proposed synchronization strategy.

We consider a network of systems (7) and combine the results of Theorem 1 and Proposition 1 to design the state feedback gain $K$ inducing synchronization over general time-invariant graphs. To this end, define the following quantities:

$$\eta_i := \left( \frac{|\lambda_i|}{\Re(\lambda_i)} \right)^2 = 1 + \left( \frac{\Im(\lambda_i)}{\Re(\lambda_i)} \right)^2, \quad i = 2, \dots, N, \quad (16a)$$

$$\bar{\eta} := \max\{\eta_2, \dots, \eta_N\}, \quad \underline{\eta} := \min\{\eta_2, \dots, \eta_N\}, \quad (16b)$$

$$\bar{\lambda} := \max_{i \in \{2, \dots, N\}} \Re(\lambda_i), \quad \lambda := \min_{i \in \{2, \dots, N\}} \Re(\lambda_i), \quad (16c)$$

for the non-zero eigenvalues $\lambda_i$, $i = 2, \dots, N$, of a Laplacian matrix $L$. Our MARI-based design is effective whenever the following inclusion holds for the graph-induced quantities (16) and the MARI parameters $\sigma \in (0, 1]$ and $c$ as in (14):

$$\frac{\bar{\eta}\sigma}{c} \in \left(0, 1 - \frac{(\bar{\eta}\bar{\lambda} - \underline{\eta}\underline{\lambda})^2}{(\bar{\eta}\bar{\lambda} + \underline{\eta}\underline{\lambda})^2}\right]. \quad (17)$$

The next lemma establishes a useful implication of (17).

<footnote n="3">³Namely, if $K$ stabilizes the pair $(A, B)$, then $\kappa K$ is stabilizing feedback for any $\kappa \ge 1$ [11, 29].</footnote>
**Lemma 2** Consider the quantities $\eta, \bar{\eta}, \lambda, \bar{\lambda}$ in (16), $c$ in (14) and let $\sigma \in (0, 1]$. The following interval of the real axis

$$ \mathcal{K} := \left[ \frac{c - \sqrt{c(c - \bar{\eta}\sigma)}}{\bar{\eta}\lambda}, \frac{c + \sqrt{c(c - \bar{\eta}\sigma)}}{\bar{\eta}\bar{\lambda}} \right] \quad (18) $$

is nonempty if and only if (17) holds.

The proof of Lemma 2 is given in Appendix A.2.

We are ready to present the main result on robust synchronization of linear systems.

**Theorem 2** Consider the network (7) and suppose that $L$ is a Laplacian matrix describing a directed and weighted communication graph with a directed spanning tree. Let $R \in \mathbb{S}_{\ge 0}^n$ and suppose that there exists $P \in \mathbb{S}_{\ge 0}^n$ such that (9) holds for a selection of $\sigma$ satisfying (17) with $\bar{\eta}, \eta, \bar{\lambda}, \lambda$ defined in (16). Then, the distributed control law $u_i$ in (4), with $\alpha(x) = \kappa K x$ and $K$ selected as in (13), solves Problem 1 for any scalar gain $\kappa \in \mathcal{K}$ as defined in (18).

The proof of Theorem 2 is presented in Appendix A.3.

We now highlight the importance of bounds (17) on $\sigma$ in Proposition 2. First, differently from the continuous-time scenario [15], the bounds on the scalar gain $\kappa$ depend on the imaginary part of the Laplacian eigenvalues via $\eta$ and $\bar{\eta}$. This is expected, as discrete-time stability requires the eigenvalues to lay inside the unit disc, thus imposing bounds on both their real and imaginary parts. Differently, continuous-time stability requires the eigenvalue to lay in the negative half-plane, which constrains their real parts only. Hence, in the case where there is at least one complex eigenvalue, definitions (16) imply $\bar{\eta} > 1$.

As a consequence, it is necessary that $\sigma < c\bar{\eta}^{-1}$ for a real solution to the square roots in (18) to exist, thus possibly excluding the DARI scenario $\sigma = 1$. In the case of real eigenvalues, $\eta_i = 1$ for all $i = 1, \dots, N$. Then, (17) and (18) simplify to

$$ \begin{aligned} \frac{\sigma}{c} &\in \left(0, 1 - \frac{(\bar{\lambda} - \lambda)^2}{(\bar{\lambda} + \lambda)^2}\right] \\ \kappa &\in \left[\frac{c - \sqrt{c(c - \sigma)}}{\bar{\lambda}}, \frac{c + \sqrt{c(c - \sigma)}}{\bar{\lambda}}\right], \end{aligned} $$

By (15), this last bound includes the results in [6], where all eigenvalues of $L$ are supposed to be real. Finally, we highlight that smaller values of $\sigma$ lead to robust synchronization over a broader range of graphs.

**Remark 3** In the continuous-time case, thanks to the infinite-gain margin property (see Remark 2), the communication graph imposes only a lower bound to the gain $\kappa$, see e.g. [15, 22].

# 4. Synchronization of Nonlinear Systems

As shown in Section 3, the problem of synchronization can be interpreted as a simultaneous stabilization problem by means of an opportune change of coordinates. In order to extend the results of Theorem 1 to the nonlinear case, in this section, we first focus on the general properties a single nonlinear system has to fulfill. Then, we exploit these properties to recast the results of Section 3.2 to networks of nonlinear agents.

## 4.1. The case of single nonlinear discrete agents

We consider time-varying discrete-time nonlinear system of the form

$$ x^+ = \varphi(x, t) + w \quad (19) $$

where the function $\varphi: \mathbb{R}^n \times \mathbb{N} \to \mathbb{R}^n$ is such that the following mild property holds.

**Property 1** Function $\varphi: \mathbb{R}^n \times \mathbb{N} \to \mathbb{R}^n$ is continuous in its first argument and there exists a (possibly unbounded) set of matrices $\mathcal{D}\varphi \subset \mathbb{R}^{n \times n}$ such that, for each $x_a, x_b \in \mathbb{R}^n$ and each $t \in \mathbb{N}$, there exists an integrable function $\psi: [0, 1] \to \mathbb{R}^{n \times n}$ satisfying $\psi(s) \in \mathcal{D}\varphi, \forall s \in [0, 1]$ and

$$ \varphi(x_a, t) - \varphi(x_b, t) = \int_0^1 \psi(s) ds (x_a - x_b). \quad (20) $$

The above definition allows considering a wide class of dynamical systems. First, note that when $n = 1$, Property 1 boils down to the requirement of $\varphi$ being absolutely continuous. Trivially, such a class of systems includes continuously differentiable ones with $\mathcal{D}\varphi$ containing all of their Jacobians. Moreover, Property 1 includes functions that are differentiable almost everywhere (i.e. everywhere but on a set of measure zero), such as piecewise smooth and Lipschitz functions. In this case, $\mathcal{D}\varphi$ contains all the possible Clarke generalized gradients [5]. As a particular case, for linear systems of the form (8), $\mathcal{D}\varphi = \{A\}$. When moving to nonlinear systems, this allows the inclusion of some useful nonlinearities, such as saturations and arctangents, by selecting $\mathcal{D}\varphi$ as the vertices of the convex hull of all possible Jacobians. Finally, in the time-varying scenario, the set $\mathcal{D}\varphi$ contains all possible generalized Jacobians of $\psi$ for all times. Hence, this set is possibly unbounded and does not depend on $t$. Instead, function $\psi$ depends on $t$, $x_a$, and $x_b$, but we omit this dependence for compact notation, and only indicate its dependence on $s$.

We now explore the design of stabilizers showing gain margin properties in the nonlinear framework by means of incremental input-to-state stability ($\delta$ISS) arguments. Typically, $\delta$ISS is obtained via incremental Lyapunov functions [1, 7, 9, 33, 34]. However, in [35, Theorem 15], the equivalence between uniform global exponentially $\delta$ISS and global contractivity is shown for continuously differentiable discrete-time dynamics. Hence, we aim at exploiting contraction to obtain $\delta$ISS.
## 4.2. Sufficient conditions for exponential δISS

While the first in-depth analysis of the relation between contraction, incremental stability and convergence in discrete-time appeared in [35], to the authors' knowledge the first results date back to [9]. We now generalize these existing results to the framework of non-smooth dynamics whose vector fields satisfy Property 1. Moreover, we extend recent advances of [7] to time-varying dynamics. We start by recalling the definition of δISS for discrete-time systems, see e.g. [1, 7, 34].

**Definition 2 (Exponential δISS) System (19)** is globally uniformly Incrementally Input-to-State Stable with exponential convergence rate (exponentially δISS) if there exist $c, \gamma > 0$ and $\rho \in [0, 1)$ such that, for all $t \ge t_0$ with $t_0 \in \mathbb{N}$ and for any initial states $x_1(t_0), x_2(t_0)$ and any pair of disturbance sequences $t \mapsto w_1(t), t \mapsto w_2(t)$, the resulting solutions $x_1(t), x_2(t)$ of (19) satisfy

$$|x_1(t) - x_2(t)| \le cp^{t-t_0}|x_1(t_0) - x_2(t_0)| + \sup_{t \in [t_0, t]} \gamma|w_1(t) - w_2(t)|. \quad (21)$$

We now introduce a sufficient δISS condition that extends the results in [7, Theorem 2], [35, Theorem 14], [9, Theorem 6.1] to the case of time-varying non-smooth vector fields satisfying Property 1.

**Lemma 3** Consider system (19) and suppose that $\varphi$ satisfies Property 1 with a specific set-valued map $\mathcal{D}\varphi$, and there exists $P \in \mathbb{S}_{>0}^n$ and $\rho \in (0, 1)$ satisfying

$$J^\top PJ \le \rho^2 P, \quad \forall J \in \mathcal{D}\varphi. \quad (22)$$

Then, system (19) is exponentially δISS according to Definition 2.

The proof of Lemma 3 is given in Appendix B.1.

Equipped with sufficient conditions for contraction of non-smooth dynamics, we conclude this subsection with some pedagogical examples, providing a useful insight into the applicability of the result.

**Example 1** Consider a system of the form

$$x^+ = \varphi(x) = \text{sat}_r(Ax), \quad (23)$$

where $r \in \mathbb{R}_{>0}^n$ and the vector saturation function $\text{sat}(\cdot)$ has components $\text{sat}_i(\cdot) := \max(\min(\cdot, r_i), -r_i)$. It can be easily verified that, for all $x \in \mathbb{R}^n$, the generalized Jacobian of $\varphi$ [5] satisfies $\partial\varphi(x) \subset \text{co}\{\Delta A, \Delta \in \Delta\}$, where $\Delta := \{\Delta = \text{diag}(\delta_1, \dots, \delta_n) : \delta_i \in \{0, 1\}, \forall i = 1, \dots, n\}$ is a finite set of matrices representing the vertices of a polytope. More generally, let $\mathcal{V}$ be a set of matrices $\mathcal{V} := \{A_1, \dots, A_v\}$ with $v \in \mathbb{N}$, such that $\partial\varphi(x) \in \text{co}\{\mathcal{V}\}$ for all $x \in \mathbb{R}^n$. Then, it suffices to verify (22) on $\mathcal{V}$ and convexity of the equivalent formulation

$$\begin{pmatrix} \rho^2 P & J^\top \\ J & P^{-1} \end{pmatrix} \succeq 0$$

(obtained via a Schur complement) ensures that (22) holds for $\mathcal{D}\varphi = \text{co}\{\mathcal{V}\}$. Similar reasonings can be followed for smooth monotone saturation-like functions, such as arctangents or hyperbolic tangents.

**Example 2** The incremental stability of neural networks is a growing research area [2, 7, 39]. Our contraction analysis tools can be applied to derive such properties. For example, a multilayer perceptron (MLP) with $L \in \mathbb{N}$ layers and ReLU activations has the following dynamics

$$\left\{ \begin{array}{l} x^+ = y_L \\ y_\ell = W_\ell \nu(y_{\ell-1}) + b_\ell, \quad \ell = 1, \dots, L \\ y_0 = x \end{array} \right. \quad (24)$$

with the ReLU function $\nu(\cdot)$ applied component-wise, i.e., $\nu(x)$ has components $\nu_i(x_i) := \max(0, x_i)$, $y_\ell, b_\ell \in \mathbb{R}^{n_\ell}$ and $W_\ell \in \mathbb{R}^{n_\ell \times n_{\ell-1}}$. Denoting by $x \mapsto \varphi(x)$ the function satisfying $y_L = \varphi(x)$ recursively defined in (24), by the chain rule [5, Theorem 2.6.6], for all $x \in \mathbb{R}^n$, we have

$$\partial\varphi(x) \subset \text{co}\{W_L\Delta_{L-1}W_{L-1}\dots\Delta_1 W_1, \\ \Delta_i \in \Delta, \forall i = 1, \dots, L-1\},$$

with $\Delta$ defined as in Example 1. Proceeding as in Example 1, we can conclude exponential δISS properties of (24) by checking (26) on a set $\mathcal{V}$ satisfying $\partial\varphi(x) \in \text{co}\{\mathcal{V}\}$ for all $x \in \mathbb{R}^n$. Similar results extend to more complex recurrent neural networks, as shown in [7].

## 4.3. Nonlinear robust feedback with gain margin design

Paralleling the linear derivation in Theorem 1, we now exploit the results of Lemma 3 to design a feedback stabilizers $u = \alpha(x, t)$ inducing exponential δISS (as defined in Definition 2) with respect to $w$ and for a nonlinear system of the form

$$x^+ = f(x, t) + Bu + w, \quad (25)$$

where $f: \mathbb{R}^n \times \mathbb{N} \to \mathbb{R}^n$ satisfies Property 1 and $B$ is full column rank. Again, the result we obtain establishes a gain margin property as defined below. We emphasize the similarity with (2), anticipating that this construction will be used in Section 4.4 to solve the robust nonlinear synchronization problem.

**Definition 3 (Gain margin for nonlinear systems)**

A function $\alpha: \mathbb{R}^n \times \mathbb{N} \to \mathbb{R}^m$ is a δISS feedback with gain margin of center $c \in \mathbb{R}$ and radius $r > 0$ if, for any real number $\kappa \in [c - r, c + r]$, system (25) with $u = \kappa\alpha(x, t)$ is exponentially δISS with respect to $w$.

The following proposition can be seen as a non-smooth nonlinear counterpart of Proposition 1. We now state our first result for the nonlinear framework.
**Proposition 2** Let $R \in \mathbb{S}_{\ge 0}^m$ and assume that $f$ in (25) satisfies Property 1 for some $\mathcal{D}f \subset \mathbb{R}^{n \times n}$. Moreover, suppose that there exists $P \in \mathbb{S}_{\ge 0}^n$ satisfying

$$J^\top Q J \le \rho^2 P, \quad \forall J \in \mathcal{D}f, \qquad (26a)$$

$$Q := P - \sigma P B (R + B^\top P B)^{-1} B^\top P, \qquad (26b)$$

for some $\rho \in (0, 1)$ and $\sigma \in (0, 1]$. Then for the system (25) the function

$$u = \alpha(x, t) = -\kappa (R + B^\top P B)^{-1} B^\top P f(x, t) \quad (27)$$

is a $\delta$ISS feedback with gain margin of center and radius

$$c = 1 + \frac{\lambda_{\min}(R)}{\lambda_{\max}(B^\top P B)}, \quad r = \sqrt{c(c - \sigma)}. \qquad (28)$$

The proof of Proposition 2 can be found in Appendix B.2. Note that for linear dynamics, the feedback (27) coincides with (13). Similar to the linear case, taking into account the relation between the choice of $R$ and the solution $P$ of (26) via the value $c$ in (28) allows for improved gain margins. We also remark that condition (26) is difficult to handle for two main reasons: i) it represents an infinite set of inequalities, ii) it is nonlinear in the variable $P$. We postpone to Section 5 a discussion on possible workarounds to these issues.

**Remark 4** The parameter $\sigma$ in (26b) measures the system's open-loop $\delta$ISS properties, analogous to its role in the linear MARI (9). Its sign determines if the system requires a stabilizing input ($\sigma > 0$), is already contracting ($\sigma = 0$), or is robust to destabilizing inputs ($\sigma < 0$). Inequality (26) can therefore be interpreted as a nonlinear version of the MARI.

**Remark 5** Mirroring linear MARI-based designs, our nonlinear $\delta$ISS approach provides a guaranteed gain margin. Any gain $\kappa \in [c - \sqrt{c(c - \sigma)}, c + \sqrt{c(c - \sigma)}]$ is stabilizing, a property analogous to findings in [11]. As with MARI, this range is centered at $c \ge 1$. Notably, the margin vanishes ($\kappa = 1$) when $\sigma = 1$ and $R$ is singular.

**Remark 6** Note that (26b) implies $Q \in \mathbb{S}_{\ge 0}^n$ if $R \in \mathbb{S}_{\ge 0}^m$. Indeed, by (26b) and (11) with $A = P$, $B = B$, $C = B^\top$ and $\mathcal{D} = R$, the invertibility of $R$ yields

$$Q = (1 - \sigma)P + \sigma(P^{-1} + B R^{-1} B^\top)^{-1}.$$

Since $\sigma \in (0, 1]$ and $P > 0$, matrix $Q$ is a $\sigma$-governed linear interpolation between positive definite matrices and, thus, it is positive.

### 4.4. Main result on nonlinear synchronization

We now exploit the $\delta$ISS results of Proposition 2 to derive a solution to the nonlinear multi-agent robust synchronization problem. The idea is to exploit contraction properties to show convergence to the synchronization manifold. The link between Proposition 1 and Proposition 2 is

evident. However, the Jordan transformation which is the basis of the proof of Theorem 1 and, consequently, of Theorem 2 cannot be easily applied in the nonlinear scenario. Thus, in order to exploit Proposition 2 for synchronization, we introduce the following technical lemma, which will be used in the proof to design an appropriate transformation.

**Lemma 4** Let the weighted graph $\mathcal{G} = \{\mathcal{V}, \mathcal{E}, \mathcal{A}\}$ be directed, with Laplacian $L$ and $L_{11}, L_{12}$ defined as in (1). Moreover, suppose $\mathcal{G}$ has a directed spanning tree. Then, there exist $M \in \mathbb{S}_{\ge 0}^{N-1}$ and constants $\underline{m}, \bar{m}, \underline{\mu}, \bar{\mu} > 0$, $\rho_M \in (0, 1]$ such that

$$\underline{m} I_{N-1} \le M \le \bar{m} I_{N-1}, \quad \rho_M \le \frac{\underline{m}}{\bar{m}} \qquad (29a)$$

$$M(L_{22} - \mathbf{1} L_{12}) + (L_{22} - \mathbf{1} L_{12})^\top M \ge 2 \underline{\mu} M \qquad (29b)$$

$$(L_{22} - \mathbf{1} L_{12})^\top M (L_{22} - \mathbf{1} L_{12}) \le \bar{\mu}^2 M. \qquad (29c)$$

The proof of Lemma 4 is given in Appendix B.3.

We now present the following main result on network synchronization of nonlinear systems.

**Theorem 3** Consider the network (2) and suppose that $f$ satisfies Property 1 for some $\mathcal{D}f \subset \mathbb{R}^{n \times n}$ and $L$ is a Laplacian matrix describing a directed and weighted communication graph with a directed spanning tree. Let $\rho \in (0, 1)$ and $\sigma \in (0, 1]$ satisfy

$$\rho \le \rho_M, \quad \sigma \le \frac{c}{\zeta}, \quad \zeta := \left(\frac{\bar{\mu}}{\underline{\mu}}\right)^2 \qquad (30)$$

with $\rho_M, \underline{\mu}, \bar{\mu}$ as in Lemma 4 and $c$ as in (28). If, for some $R \in \mathbb{S}_{\ge 0}^m$, there exists $P \in \mathbb{S}_{\ge 0}^n$ satisfying (26a), (26b), then, the distributed control law $u_i = \sum_{j=1}^N \ell_{ij}\alpha(x_j, t)$ in (4) with $\alpha$ defined as in (27) and $\kappa$ satisfying

$$\kappa \in \left[ \frac{c - \sqrt{c(c - \zeta\sigma)}}{\zeta\underline{\mu}}, \frac{c + \sqrt{c(c - \zeta\sigma)}}{\zeta\underline{\mu}} \right], \qquad (31)$$

solves Problem 1 for the network (2), namely, (6) holds.

The proof of Theorem 3 can be found in Appendix B.4.

The difference between the allowed set of gains (31) and (18) arises from the Lyapunov-like conditions (29). Indeed, while Theorem 3 builds on Lemma 4, (18) follows from spectral arguments. Nevertheless, while more conservative, the matrix inequalities (29) can be related to spectral properties of $L$, as presented in the next remark.

**Remark 7** The contraction inequality (26a) in Theorem 3 is tightly related to the structure of the Laplacian matrix and its eigenvalues. For an undirected, leader-connected network graph, the Laplacian $L$ in (1) satisfies $L_{11} = 0, L_{12} = 0$, and $L_{22}$ is symmetric positive definite. Choosing $M = \text{Id}_{N-1}$, $\underline{\mu}$ (resp. $\bar{\mu}$) in Lemma 4 can be selected as the smallest (resp. largest) eigenvalue of $L_{22}$, recovering a result similar to the linear scenario.
Consequently, admissible $\sigma$ values relate to the condition number of $L_{22}$. Furthermore, with $M = \text{Id}_{N-1}$, the contraction rate $\rho$ disentangles from the network structure ($\underline{m} = \overline{m} = 1$), meaning condition (30) imposes no constraints on $\rho \in (0,1)$, as in Lemma 3.

## 5. GEVPs for exponential $\delta$ISS

We now discuss numerically efficient formulations of the results of Section 4, and apply them to specific classes of systems. LMI-based conditions for robust stabilization are a valuable tool for control design for discrete-time nonlinear systems, see e.g. [12, 19, 30, 43]. Hence, inspired by these works and recent LMI approaches for solving MARI inequalities [31, 36], we propose LMI-based conditions to obtain the solution of the MARI-like inequality (26). Our result provides a viable solution to the design problem of robustly synchronizing controllers.

First, we introduce an equivalent formulation for (26). The parameters of the proposed reformulation can be obtained by solving a generalized eigenvalue problem (GEVP). Then, we focus our attention on the case where the set of possible open-loop Jacobians of the system dynamics (25) is polytopic. Finally, we target the specific case of Lur'e systems.

### 5.1. Formulation as a GEVP

We start by reformulating Proposition 2 as an LMI problem. This provides convex analysis conditions for constructing matrix $P$.

**Proposition 3** Let $\sigma \in (0,1]$, $\rho \in (0,1)$. The following LMIs in the decision variables $W, \Sigma$ and $U$

$$W \succ 0, \quad \Sigma \succ 0, \quad U \succ 0, \quad \begin{pmatrix} \rho W & W J^\top \\ J W & \rho \Sigma \end{pmatrix} \succeq 0, \quad (32a)$$

$$\begin{pmatrix} W + \sigma B U B^\top - \Sigma & B U B^\top \\ B U B^\top & \frac{1}{\sigma} \left( \frac{W}{1-\sigma} + B U B^\top \right) \end{pmatrix} \succeq 0, \quad (32b)$$

hold if and only if conditions (26) hold with $P = W^{-1} \succ 0$ and $R = U^{-1} \in S_{\succ 0}^m$. Moreover, (32) is a generalized eigenvalue problem in $(\sigma, \rho)$, namely, if it is feasible for $(\bar{\sigma}, \bar{\rho})$, then it is feasible for any $(\sigma, \rho)$ such that $\sigma \ge \bar{\sigma}$, $\rho \ge \bar{\rho}$ and, conversely, if it is infeasible for $(\underline{\sigma}, \underline{\rho})$, then it is infeasible for any $(\sigma, \rho)$ such that $\sigma \le \underline{\sigma}$, $\rho \le \underline{\rho}$.

The proof of Proposition 3 can be found in Appendix B.5. By continuity arguments, if $W, \Sigma, U$ solve (32) for some $J, \rho, \sigma$, then there exist scalars $\hat{\rho} \in (\rho, 1)$, $\hat{\sigma} \in (\sigma, 1]$ and (typically small) $\delta > 0$, such that (32) holds for all $\hat{J}$ satisfying $|J - \hat{J}| \le \delta$. In other words, controller (27) provides closed-loop robustness to (small) uncertainties in $J$, and consequently in $Df$. In fact, (32) allows easily addressing uncertain nonlinearities whose (generalized) Jacobian is known to be bounded in an interval (see Section 5.2 and the example in Section 6.1). Moreover, the fact that

(26) is a GEVP implies that optimized values of $\sigma$ and $\rho$ can be estimated using iterative methods, such as bi-section or Newton-Rhapson iterations. This allows establishing performance/robustness tradeoffs associated with the fact that $\rho$ governs the guaranteed convergence rate, through (6), while $\sigma$ governs the guaranteed robustness, through (31). In fact, note that lower values for $\sigma$ in (31) imply that less control effort is required to incrementally stabilize the system (because the lower bound on $\kappa$ approaches zero), in addition to enlarging the gain margin radius $(\zeta\mu)^{-1}\sqrt{c(c - \zeta\sigma)}$.

Combined with Proposition 2, Proposition 3 requires the satisfaction of (32) for all $J \in Df$. This may turn out to be impractical, as $Df$ could be infinite-dimensional. However, under some additional assumptions on system (25), we can follow a polytopic approach similar to the one in Examples 1 and 2. Hence, we propose the following result addressing the case where the open-loop system Jacobian belongs to a polytopic set defined by a finite number of vertices.

**Corollary 1** Let $R \in S_{\succ 0}^m$ and assume that $f$ in (25) satisfies Property 1 for some $Df \subset \mathbb{R}^{n \times n}$. Moreover, suppose there exists a finite set of matrices $V := \{A_1, \dots, A_v\} \subset \mathbb{R}^{n \times n}$ such that $Df \subseteq co\{V\}$. If there exist matrices $W, \Sigma \in S_{\succ 0}^n$ and scalars $\rho \in (0,1)$, $\sigma \in (0,1]$ satisfying (32) for all $J \in V$, the control law $u = \alpha(x, t)$ with $\alpha$ defined in (27) and $P = W^{-1}$ makes the closed-loop exponentially $\delta$ISS with respect to $w$ with gain margin of center and radius (28).

### 5.2. Lur'e systems

We further specialize our result to the case of Lur'e systems. Namely, we now consider nonlinear discrete-time systems of the form

$$x^+ = f(x) + Bu = Ax + F\phi(Cx) + Bu, \quad (33)$$

where $C \in \mathbb{R}^{p \times n}$, $F \in \mathbb{R}^{n \times p}$ and the square nonlinearity $\phi: \mathbb{R}^p \to \mathbb{R}^p$ is a pool of $p$, possibly different, feedback nonlinear elements $\phi(y) := \text{diag}(\phi_1(y_1), \dots, \phi_p(y_p))$ whose components $\phi_i, i = 1, \dots, p$, satisfy Property 1 for some intervals $D\phi_i \subset \mathbb{R}$, $i = 1, \dots, p$. We assume that each function $\psi_i$ belongs to an incremental sector $[0, \omega_i]$, with $\omega_i \ge 0$, in the following classical sense:

$$(\phi_i(s_1) - \phi_i(s_2))(\phi_i(s_1) - \phi_i(s_2) - \omega_i(s_1 - s_2)) \le 0, \quad (34)$$

for all $s_1, s_2 \in \mathbb{R}$. Conditions of the form (34) are useful to tackle uncertainties in the slope of the nonlinearity $\phi$. By the non-smooth mean value theorem [5, Theorem 2.3.7], we may combine bounds (34) into

$$2\text{Sym}\{J_{\phi}S(J_{\phi} - \Omega)\} \le 0 \quad (35)$$

which holds for all diagonal $J_{\phi} \in D\phi = \text{diag}(D\phi_1, \dots, D\phi_p)$, for any diagonal $S \in S_{\succ 0}^p$ and for some diagonal $\Omega = \text{diag}(\omega_1, \dots, \omega_p) \in S_{\succ 0}^p$. We then have the following result.
**Proposition 4** Let $R \in \mathbb{S}_{\ge 0}^m$ and suppose that $\phi$ in (33) satisfies (34) for some $\Omega = \text{diag}(\omega_1, \dots, \omega_p)$. If there exist symmetric matrices $W, \Sigma$, a diagonal matrix $S \in \mathbb{S}_{\ge 0}^p$ and scalars $\rho \in (0, 1), \sigma \in (0, 1]$ satisfying (32b) and

$$ W \succ 0, \quad \Sigma \succ 0, \quad \begin{pmatrix} \rho W & W A^\top & W C^\top \\ A W & \rho \Sigma & -F \Omega^\top S \\ C W & -S \Omega F^\top & 2S \end{pmatrix} \succeq 0, \quad (36) $$

then the control law $u = \alpha(x, t)$ with $\alpha$ defined in (27) and $P = W^{-1}$ makes the closed-loop (33)-(27) exponentially $\delta$ISS with respect to $w$ with gain margin of center and radius (28).

The proof of Proposition 4 can be found in Appendix B.6.

**Remark 8** Our result in (36) differs from related work in several important ways. Unlike [17, Eq. (13)], which handles classical stabilization for time-invariant systems, our Proposition 4 ensures incremental stability for time-varying systems with a guaranteed convergence rate via $\rho$. Furthermore, our Riccati-based design induces robustness margins of incremental stability that are key for achieving synchronization and decouples the plant and controller analysis, in contrast to the unified LMI structure in [17]. Finally, compared to the analysis in [16, Section IV.B], our work focuses on control design rather than analysis. Nonetheless, the more general incremental Lyapunov functions used therein represent a promising direction for future work.

# 6. Numerical example

## 6.1. Polytopic approach

Consider the planar dynamics for agent $i$ in (2)

$$ \underbrace{\begin{pmatrix} x_a \\ x_b \end{pmatrix}}_{x_i^+} = \underbrace{\begin{pmatrix} f_1(x_a, x_b, t) \\ f_2(x_a, x_b, t) \end{pmatrix}}_{f(x_i, t)} + \underbrace{\begin{pmatrix} 1 \\ 1 \end{pmatrix} \sin(t) \sin(x_a)}_{\tilde{f}(x_i, t)} + \underbrace{\begin{pmatrix} 1 \\ 0 \end{pmatrix}}_{B} u_i, $$

with $x_i = \text{col}(x_a, x_b)$, where we omit the index “i” from $x_a, x_b$ for simplicity, and

$$ f_1(x_a, x_b, t) = k_1 x_a + k_2 x_b + k_5 e^{-x_a^2} + k_6 \sin(x_b) + k_8 (1 + \sin(\omega_1 t)) \ln(1 + \max(0, x_a)) + k_9 \tanh(x_b), $$

$$ f_2(x_a, x_b, t) = k_3 x_a + k_4 e^{-x_b^2} + k_7 \cos(x_a) + k_{10} \arctan(x_a) + 2 \sin(\sin(t)\omega_2 t), $$

and positive parameters $k_1 = 1.1, k_2 = 0.02, k_3 = 0.085, k_4 = 0.4\sqrt{2e}, k_5 = \frac{0.3\sqrt{2e}}{2}, k_6 = 0.1, k_7 = 0.005, k_8 = 0.05, k_9 = 0.03, k_{10} = 0.01, \omega_1 = \pi/8, \omega_2 = \sqrt{e}$. To highlight the robustness against unmodeled dynamics inherent in the results of Section 4, we design $u$ using only the function $f$, thereby omitting $\tilde{f}$ from the computation of $\alpha$ in (27).

<figure><img src="image_2.png" /><figcaption>Figure 1: Communication graphs considered in Section 6.</figcaption></figure>

Since $\ln(1 + \max(0, x_a))$ is Lipschitz continuous on any compact set, $f$ satisfies Property 1. Moreover, $f_a$ is not differentiable on the line $x_a = 0$, due to the max function and because $k_8 \neq 0$. For $x_a < 0$, we have

$$ J_1 := \frac{\partial f}{\partial x}(x, t) = \begin{pmatrix} j_1 & j_2 \\ j_3 & j_4 \end{pmatrix}. $$

where $j_1 = k_1 - 2k_5 x_a e^{-x_a^2}$, $j_2 = k_2 + k_6 \cos(x_b) + k_9 \text{sech}(x_b)^2$, $j_3 = k_3 - k_7 \sin(x_a) + \frac{k_{10}}{1+x_a^2}$, $j_4 = -2k_4 x_b e^{-x_b^2}$. For $x_a > 0$, we obtain

$$ J_2 := \frac{\partial f}{\partial x}(x, t) = \begin{pmatrix} j_1 + \frac{k_8(1+\sin(\omega_1 t))}{1+x_a} & j_2 \\ j_3 & j_4 \end{pmatrix}. $$

As commented after (20), since $f$ is Lipschitz continuous, we can use Clarke's generalized gradient and select

$$ \mathcal{D}f = \bigcup_{\alpha \in [0, 2k_8]} \left\{ \begin{pmatrix} j_1 + \frac{\alpha \operatorname{sgn}(\max\{0, x_a\})}{1+\max\{0, x_a\}} & j_2 \\ j_3 & j_4 \end{pmatrix} \right\} $$

with $\operatorname{sgn}(0) = 0$. Note that $\mathcal{D}f$ has an infinite number of elements, but it is a bounded set. In fact, each $J = (\begin{pmatrix} j_1 & j_2 \\ j_3 & j_4 \end{pmatrix}) \in \mathcal{D}f$ belongs to the polytope

$$ \underline{j}_1 \le \underline{j}_1 \le \bar{j}_1, \quad \underline{j}_2 \le \underline{j}_2 \le \bar{j}_2, \quad \underline{j}_3 \le \underline{j}_3 \le \bar{j}_3, \quad \underline{j}_4 \le \underline{j}_4 \le \bar{j}_4 $$

with the following selections, exploiting the bound $-\frac{1}{\sqrt{2e}} \le ye^{-y^2} \le \frac{1}{\sqrt{2e}}$ for all $y \in \mathbb{R}$,

$$ \underline{j}_1 = k_1 - \frac{2k_5}{\sqrt{2e}} = 0.8, \quad \bar{j}_1 = k_1 + \frac{2k_5}{\sqrt{2e}} + 2k_8 = 1.5, $$

$$ \underline{j}_2 = k_2 - k_6 = -0.08, \quad \bar{j}_2 = k_2 + k_6 + k_9 = 0.15, $$

$$ \underline{j}_3 = k_3 - k_7 = 0.08, \quad \bar{j}_3 = k_3 + k_7 + k_{10} = 0.1, $$

$$ \underline{j}_4 = -\frac{2k_4}{\sqrt{2e}} = -0.8, \quad \bar{j}_4 = \frac{2k_4}{\sqrt{2e}} = 0.8. $$

Therefore, we leverage Proposition 3 and solve (32) for all $J$ at the vertices of the polytope. Specifically, we fix $\rho = 0.9$ and $\sigma = 0.7$ and solve 20 LMIs: the three strict positivity constraints in (32a), $2^4 = 16$ non-strict LMIs in
<figure><img src="image_3.png" /><figcaption>Figure 2: Time evolution of the multi-agent system trajectories and average error norms from the randomly sampled initial condition $x(0) = \text{col}(9.9343, -2.7653, 12.9538, 30.4606, -4.6831, -4.6827, 31.5843, 15.3487, -9.3895, 10.8512)$. (a) State trajectories of the open-loop multi-agent system without noise; (b) State trajectories of the closed-loop multi-agent system with noise; (c) Average agents' error norms wrt Agent 1, with and without noise. Shadowed regions represent variance.</figcaption></figure>

(32a) corresponding to the polytope vertices, and one LMI (32b). We obtain the solution $U = 15.4705$ and

$$W = \begin{pmatrix} 0.0619 & -0.0018 \\ -0.0018 & 0.0857 \end{pmatrix}, \quad \Sigma = \begin{pmatrix} 0.2025 & -0.0006 \\ -0.0006 & 0.0847 \end{pmatrix}.$$

We consider $N = 5$ agents interconnected by the graph in Figure 1a represented by the Laplacian matrix

$$L = \begin{pmatrix} 0.1 & 0 & 0 & -0.1 & 0 \\ -3 & 3.1 & 0 & 0 & -0.1 \\ -2.9 & -0.1 & 3 & 0 & 0 \\ -2.7 & 0 & -0.1 & 2.8 & 0 \\ -3 & 0 & 0 & 0 & 3 \end{pmatrix}.$$

This graph has a directed spanning tree rooted in node 1. Then, by Lemma 4, inequalities (29) hold with $\bar{m} = 2.4414$, $\bar{m} = 2.7092$, $\underline{\mu} = 2.7174$, $\bar{\mu} = 3.1378$, and

$$M = \begin{pmatrix} 2.5530 & 0.0561 & 0.0034 & -0.1052 \\ 0.0561 & 2.5132 & 0.0155 & 0.0075 \\ 0.0034 & 0.0155 & 2.5389 & -0.0099 \\ -0.1052 & 0.0075 & -0.0099 & 2.6316 \end{pmatrix}.$$

With the above solution, conditions (30) and (31) hold with $\rho_M = 0.9011$, $c = 1.0040$, $c/\zeta = 0.7530$, and $\kappa \in [0.2036, 0.3506]$. In the simulation, the gain is selected as $\kappa = 0.2771$. We simulate the multi-agent system with perturbed dynamics

$$x_i^+ = f(x_i, t) + \tilde{f}(x_i, t) + Bu_i + w_i, \quad i = 1, \dots, 5, \quad (37)$$

where the random noise $w_i$ is sampled from a uniform distribution $w_i \sim \mathcal{U}(-\bar{w}, \bar{w}) \times \mathcal{U}(-\bar{w}, \bar{w})$ with bounds $\bar{w} = 0.5$ for all $i = 1, \dots, 5$. Simulation results are shown in Figure 2 for a randomly sampled initial condition.

Figure 2a shows that the single agent dynamics are open-loop unstable even when $w_i = 0$ for all $i$, and the network does to synchronize. Figure 2b depicts the closed-loop trajectories of (37). Despite the unmodeled dynamics and the additive noise, the network achieves synchronization to a nontrivial (unbounded) trajectory. Figure 2c details the average of $|x_i - x_1|$ over $i = 2, \dots, 5$ alongside

the corresponding variance. The blue line illustrates the exponential error decay when $w_i = 0$ for all $i = 1, \dots, 5$, while the orange one shows the $\delta$ISS robustness properties yielding ultimately bounded inter-agent errors.

## 6.2. Lur'e systems

We consider now a network of $N = 6$ agents connected according to the graph in Figure 1b represented by

$$L = 0.2 * \begin{pmatrix} 2 & -1 & 0 & 0 & 0 & -1 \\ -1 & 3 & 0 & -1 & 0 & -1 \\ 0 & 0 & 3 & -1 & -1 & -1 \\ -2 & 0 & 0 & 4 & -1 & -1 \\ 0 & -1 & -1 & 0 & 3 & -1 \\ 0 & -1 & -1 & 0 & -1 & 3 \end{pmatrix},$$

and evolving according to the Lur'e dynamics (33) where

$$A = \begin{pmatrix} 1.1 & 0.1 \\ -0.3 & 0.5 \end{pmatrix}, \quad B = \begin{pmatrix} 0.2 \\ 0.3 \end{pmatrix}, \\ C = (1 \ -1), \quad F^{\top} = (-0.1 \ 0.6),$$

and $\phi(\cdot) = \text{sat}_3(\cdot) = \max(\min(\cdot, 3), -3)$. It is simple to verify that $J_{\phi} \in \{0, 1\}$. Unlike the example in Section 6.1, the infinite set of LMIs in (32a) can be replaced by the four conditions in (36). Consequently, we need to solve only 5 LMIs, corresponding to (36) and (32b). Conditions (35), (36) and (32b) are solved with $\rho = \rho_M - 0.001 = 0.999$, $\sigma = \zeta^{-1} + 0.001 = 0.101$, $S = \Omega = 1$ and provide $U = 2.8981$ and

$$W = \begin{pmatrix} 0.0455 & -0.2400 \\ -0.2400 & 1.3805 \end{pmatrix}, \quad \Sigma = \begin{pmatrix} 0.0458 & -0.2396 \\ -0.2396 & 1.3807 \end{pmatrix}.$$

As $c = 1.0205$, the recovered values satisfy (30) and (31). Similar to the example in Section 6.1, we simulate the multi-agent system with perturbed dynamics

$$x_i^+ = Ax_i + F\phi(Cx_i) + Bu_i + 0.2\sin(t)\sin(x_i) + w_i, \quad (38)$$

where $\sin(x_i) = \text{col}(\sin(x_a), \sin(x_b))$ and the random noise $w_i$ is sampled from a uniform distribution $w_i \sim \mathcal{U}(-\bar{w}, \bar{w}) \times \mathcal{U}(-\bar{w}, \bar{w})$ with bounds $\bar{w} = 0.5$ for all $i = 1, \dots, 6$. Simulation results are shown in Figure 3 for a randomly sampled initial condition.
<figure><img src="image_4.png" /><figcaption>Figure 3: Time evolution of the Lur'e multi-agent system trajectories and average error norms from the randomly sampled initial condition $\mathbf{x}(0) = \text{col}(88.2026, 20.0079, 48.9369, 112.0447, 93.3779, -48.8639, 47.5044, -7.568, -5.1609, 20.5299, 7.2022, 72.7137). (a) State trajectories of the open-loop multi-agent system without noise; (b) State trajectories of the closed-loop multi-agent system with noise; (c) Average agents' error norms wrt Agent 1, with and without noise. Shadowed regions represent variance.</figcaption></figure>

## 7. Conclusions

In this paper, we explored the design of incrementally stabilizing feedback laws for discrete-time nonlinear systems and their application to multi-agent synchronization under generic connected communication graphs. Starting from the linear scenario, we provided constructive designs for robust stabilization and sufficient convex conditions for network synchronization. We exploited new contraction analysis results and focused on Euclidean metrics. The analysis is focused on input-linear systems. Future work will focus on generalizing these results to systems with different input-vector fields, potentially using non-quadratic Lyapunov functions or Riemannian metrics. We also plan to explore more precise robustness bounds, incorporate learning-based techniques to reduce reliance on model knowledge, and analyze the impact of discrete-time observers to relax full-state information assumptions.

**Acknowledgement.** We thank Mattia Giaccagli for providing the main steps of the proof of Lemma 4.

## Appendix A. Proofs of the Linear results

### Appendix A.1. Proof of Theorem 1

*Sufficiency.* Using the Kronecker notation, the closed-loop network dynamics can be written as

$$ \mathbf{x}^{+} = ((\mathbf{I}_N \otimes \mathbf{A}) + (\mathbf{L} \otimes \mathbf{BK}))\mathbf{x} + \mathbf{w}, \quad (A.1) $$

with $\mathbf{x}$ and $\mathbf{w}$ defined in (3). To show convergence to the synchronization manifold $\mathcal{M}$, define a virtual leader $z$ as the root node of the directed spanning tree. Without loss of generality, assume $z := x_1$. Consider now the associated transformation

$$ T := \begin{pmatrix} 1 & 0 \\ -\mathbf{1} & \mathbf{I}_{N-1} \end{pmatrix}, \quad T^{-1} = \begin{pmatrix} 1 & 0 \\ \mathbf{1} & \mathbf{I}_{N-1} \end{pmatrix}. \quad (A.2) $$

According to the partitioning in (1), we have

$$ \mathcal{L} := \mathbf{T}\mathbf{L}\mathbf{T}^{-1} = \begin{pmatrix} 0 & L_{12} \\ \mathbf{0} & \tilde{L} \end{pmatrix}, \quad \tilde{L} := L_{22} - \mathbf{1}L_{12}, \quad (A.3) $$

where we exploited $L\mathbf{1} = \mathbf{0}$ (see Section 2.1). As a consequence, defining $N-1$ error coordinates as $\mathbf{e} := \text{col}(e_2, \dots, e_N) \in \mathbb{R}^{(N-1)n}$ with $e_i := x_i - z$, we obtain $\mathbf{e} = (0 \ \mathbf{I}_{N-1})\mathbf{T}\mathbf{x}$ and

$$ \begin{aligned} \mathbf{e}^{+} &= (\mathbf{0} \quad \mathbf{I}_{N-1}) \\ &\quad \times \left( ((TT^{-1} \otimes \mathbf{A}) + (\mathcal{L} \otimes \mathbf{BK})) \begin{pmatrix} z \\ \mathbf{e} \end{pmatrix} + (\mathbf{T} \otimes \mathbf{I}_n)\mathbf{w} \right) \\ &= A_{\text{cl}}\mathbf{e} + \tilde{\mathbf{w}}, \end{aligned} \quad (A.4) $$

$$ A_{\text{cl}} := (\mathbf{I}_{N-1} \otimes \mathbf{A}) + (\tilde{\mathbf{L}} \otimes \mathbf{BK}), \quad (A.5) $$

where $\tilde{\mathbf{w}} := \text{col}(\tilde{w}_2, \dots, \tilde{w}_N) \in \mathbb{R}^{Nn}$ with $\tilde{w}_i := w_i - w_1$. If $A_{\text{cl}}$ is Schur-Cohn stable, the use of standard arguments for linear systems yields

$$ |\mathbf{e}(t)| \le c \rho^{t-t_0} |\mathbf{e}(t_0)| + \sup_{t \in [t_0, t]} \gamma |\tilde{\mathbf{w}}(t)|, \quad (A.6) $$

for some $c, \gamma > 0$ and $\rho \in (0, 1)$. Since $T$ in (A.2) has a bounded norm, we have, for some $c_1 > 0$,

$$ \begin{aligned} |\mathbf{e}|^2 &= \left| \begin{pmatrix} z \\ \mathbf{e} \end{pmatrix} - \begin{pmatrix} z \\ \mathbf{0} \end{pmatrix} \right|^2 = \inf_{z^* \in \mathbb{R}^n} \left| \begin{pmatrix} z \\ \mathbf{e} \end{pmatrix} - \begin{pmatrix} z^* \\ \mathbf{0} \end{pmatrix} \right|^2 \\ &= \inf_{\mathbf{x}^* \in \mathcal{M}} |(\mathbf{T} \otimes \mathbf{I}_n)\mathbf{x} - (\mathbf{T} \otimes \mathbf{I}_n)\mathbf{x}^*|^2 \\ &\le \inf_{\mathbf{x}^* \in \mathcal{M}} c_1^{-1} |\mathbf{x} - \mathbf{x}^*|^2 = c_1^{-1} |\mathbf{x}|_{\mathcal{M}}^2 \end{aligned} $$

Also, since $T$ is invertible, for some $c_2 > 0$,

$$ \begin{aligned} |\mathbf{x}|_{\mathcal{M}}^2 &= \inf_{\mathbf{x}^* \in \mathcal{M}} |\mathbf{x} - \mathbf{x}^*|^2 \\ &= \inf_{\mathbf{x}^* \in \mathcal{M}} |(\mathbf{T}^{-1} \otimes \mathbf{I}_n)(\mathbf{T} \otimes \mathbf{I}_n)(\mathbf{x} - \mathbf{x}^*)|^2 \\ &\le \inf_{z^* \in \mathbb{R}^n} c_2 \left| \begin{pmatrix} z \\ \mathbf{e} \end{pmatrix} - \begin{pmatrix} z^* \\ \mathbf{0} \end{pmatrix} \right|^2 = c_2 |\mathbf{e}|^2. \end{aligned} $$

Then, we obtain the following relations

$$ \sqrt{c_1}|\mathbf{e}| \le |\mathbf{x}|_{\mathcal{M}} \le \sqrt{c_2}|\mathbf{e}|, \quad |\tilde{\mathbf{w}}| \le \sup_{i,j \in [1,N]} c_3|w_i - w_j|, \quad (A.7) $$
for some $c_1, c_2, c_3 > 0$. As a consequence, if $A_{cl}$ is Schur-Cohn stable, one obtains robust synchronization as in Problem 1. Therefore, in the rest of the proof we set $\tilde{w} = 0$ and we show that $A_{cl}$ in (A.5) is Schur-Cohn stable. Let $T_J \in \mathbb{C}^{(N-1)\times(N-1)}$ be a transformation such that $\tilde{L}_J := T_J(L_{22} - \mathbf{1} L_{12})T_J^{-1}$ is in Jordan canonical form. By defining the resulting closed-loop matrix after change of coordinates

$$ \hat{A}_{cl} = (T_J \otimes I_{N-1}) A_{cl} (T_J^{-1} \otimes I_{N-1}), $$

Schur-Cohn stability of $\hat{A}_{cl}$ implies Schur-Cohn stability of $A_{cl}$. By the properties of the Kronecker product and (A.5), we have

$$ \hat{A}_{cl} = (I_{N-1} \otimes A) + (\tilde{L}_J \otimes BK). \quad (A.8) $$

Since $\tilde{L}_J$ is in its Jordan form, the former matrix is block triangular with diagonal block equal to $(A + \lambda BK)$ with $\lambda$ in spec $\tilde{L}_J$.

Hence, Schur-Cohn stability of $\hat{A}_{cl}$ holds if and only if the complex matrix $(A + \lambda BK)$ is Schur-Cohn stable for all $\lambda \in \text{spec } \tilde{L}_J$. Due to the similarity transformations, spec $\tilde{L}_J = \text{spec } (L_{22} - \mathbf{1} L_{12}) = \text{spec } L \setminus \{\lambda_1\}$, where $\lambda_1 = 0$ is the eigenvalue associated to the eigenvector $\mathbf{1}$. Consider now the case in item (i) where the graph has a directed spanning tree. Then, by Lemma 1, $L$ has only one zero eigenvalue and the gain $K$ is such that $(A + \lambda BK)$ is Schur-Cohn stable for all $\lambda \in \text{spec } \tilde{L}_J$, which implies Schur-Cohn stability of $A_{cl}$ in (A.5). If instead $A$ is Schur-Cohn in item (i), then $(A + \lambda BK)$ is Schur-Cohn stable for all eigenvalues $\lambda$ of $L$ (including the zero ones) and $A_{cl}$ in (A.5) is exponentially stable.

*Necessity.* Consider a Laplacian matrix of the form (1). Following the lines of the sufficiency proof, the error dynamics between agents and a virtual leader are described by (A.4). We first study the connectivity requirement. Suppose that synchronization is achieved, $A$ is unstable and at least one agent is not connected. Without loss of generality, assume $x_1$ to be such a node. Since it is not connected, the Laplacian takes the form

$$ \mathcal{L}' = \begin{pmatrix} 0 & \mathbf{0}^{\top} \\ \mathbf{0} & \tilde{L}' \end{pmatrix}, $$

where $\tilde{L}'$ is the Laplacian matrix of the connected portion of the graph. Then, by (A.4) with $\tilde{w} = 0$, we have

$$ \mathbf{e}^{+} = ((I_{N-1} \otimes A) + (\tilde{L}' \otimes BK))\mathbf{e}. $$

Notice that $L'$ describes a graph with a directed spanning tree. Then, by Lemma 1, it has one zero eigenvalue. By performing similar steps to the ones in the sufficiency proof, we define the transformed closed-loop matrix

$$ \hat{A}_{cl} = (I_{N-1} \otimes A) + (\tilde{L}'_J \otimes BK), $$

where $\tilde{L}'_J$ is in Jordan form. Note that $\hat{A}_{cl}$ is Schur-Cohn stable if and only if the complex matrix $(A + \lambda BK)$ is

Schur-Cohn stable for all $\lambda \in \text{spec } \tilde{L}'_J$. However, spec $\tilde{L}'_J$ includes a zero eigenvalue. Hence, $\hat{A}_{cl}$ is stable if and only if $A$ is Schur-Cohn stable, showing the first item by establishing contradiction.

We now prove the necessity of item (ii). If the agents are synchronized, the $e$ subsystem in (A.4) is asymptotically stable and the matrix $\hat{A}_{cl}$ in (A.8) is Schur-Cohn. Since $\tilde{L}_J$ contains all the nonzero eigenvalues of $L$ and $\hat{A}_{cl}$ is block-upper triangular, item (ii) must hold, and this concludes the proof.

## Appendix A.2. Proof of Lemma 2

It is trivial that (17) implies $\bar{\eta}\sigma \le c$ and then the square root in (2) is well-defined. To complete the proof, we show that the bound

$$ \frac{c - \sqrt{c(c - \bar{\eta}\sigma)}}{\bar{\eta}\lambda} \le \frac{c + \sqrt{c(c - \bar{\eta}\sigma)}}{\bar{\eta}\bar{\lambda}} \quad (A.9) $$

holds if and only if (17) is satisfied, namely if and only if

$$ \frac{\bar{\eta}\sigma}{c} \le 1 - \frac{(\bar{\eta}\bar{\lambda} - \underline{\eta}\underline{\lambda})^2}{(\bar{\eta}\bar{\lambda} + \underline{\eta}\underline{\lambda})^2} = \frac{4\bar{\eta}\underline{\eta}\bar{\lambda}\underline{\lambda}}{(\bar{\eta}\bar{\lambda} + \underline{\eta}\underline{\lambda})^2} = \left( \frac{2\bar{\eta}\bar{\lambda}}{\bar{\eta}\bar{\lambda} + \underline{\eta}\underline{\lambda}} \right)^2 \frac{\underline{\eta}\underline{\lambda}}{\bar{\eta}\bar{\lambda}}, $$

which, due to the positivity of the squared term, is equivalent to

$$ \left( \frac{\bar{\eta}\bar{\lambda} + \underline{\eta}\underline{\lambda}}{2\bar{\eta}\bar{\lambda}} \right)^2 \bar{\eta}\sigma - c \frac{\underline{\eta}\underline{\lambda}}{\bar{\eta}\bar{\lambda}} \le 0. \quad (A.10) $$

Thus, we must show that (A.9) $\iff$ (A.10). By the lower bound of (17), $\bar{\eta}\sigma > 0$. Then, multiplying (A.10) by $\bar{\eta}\sigma$ paired with addition and subtraction of $c(c - \bar{\eta}\sigma)$ at the right-hand side yields the equivalent inequality

$$ \begin{aligned} c(c - \bar{\eta}\sigma) &\ge \left( \frac{\bar{\eta}\bar{\lambda} + \underline{\eta}\underline{\lambda}}{2\bar{\eta}\bar{\lambda}} \right)^2 (\bar{\eta}\sigma)^2 + c^2 - \left( 1 + \frac{\underline{\eta}\underline{\lambda}}{\bar{\eta}\bar{\lambda}} \right) \bar{\eta}\sigma c \\ &\ge \left( c - \frac{\bar{\eta}\bar{\lambda} + \underline{\eta}\underline{\lambda}}{2\bar{\eta}\bar{\lambda}} \right)^2 \bar{\eta}\sigma^2. \end{aligned} \quad (A.11) $$

By taking the square root, (A.11) is equivalent to

$$ \sqrt{c(c - \bar{\eta}\sigma)} \ge c - \frac{\bar{\eta}\bar{\lambda} + \underline{\eta}\underline{\lambda}}{2\bar{\eta}\bar{\lambda}} \bar{\eta}\sigma, \quad (A.12) $$

where the right-hand side is non-negative since $\bar{\eta}\underline{\lambda} \le \bar{\eta}\bar{\lambda}$ and $\bar{\eta}\sigma \le c$. Exploiting the expansion $c\bar{\eta}\sigma = (c - \sqrt{c(c - \bar{\eta}\sigma)})(c + \sqrt{c(c - \bar{\eta}\sigma)}) > 0$, inequality (A.12) is equivalent to

$$ \begin{aligned} \frac{\bar{\eta}\underline{\lambda}}{\bar{\eta}\bar{\lambda}} &\ge \left( c - \sqrt{c(c - \bar{\eta}\sigma)} \right) \frac{2}{\bar{\eta}\sigma} - 1 \\ &\ge \frac{c^2 + c(c - \bar{\eta}\sigma) - 2c\sqrt{c(c - \bar{\eta}\sigma)}}{c\bar{\eta}\sigma} \\ &\ge \frac{\left( c - \sqrt{c(c - \bar{\eta}\sigma)} \right)^2}{(c - \sqrt{c(c - \bar{\eta}\sigma)})(c + \sqrt{c(c - \bar{\eta}\sigma)})} \end{aligned} $$

which coincides with (A.9), thus completing the proof.
## Appendix A.3. Proof of Theorem 2

By Theorem 1, Problem 1 is solved (equivalently, (6)) holds) if the matrices $(A + \lambda_i \kappa BK)$ are Schur-Cohn for all $\lambda_i \in \text{spec } L \setminus \emptyset$. By Proposition 1, each one of these matrices is Schur-Cohn if $|\kappa \lambda_i - c|^2 \le c(c - \sigma)$. By expanding the norm, we conclude that the closed-loop matrix associated to $\lambda_i$ is Schur-Cohn if $c\sigma - 2c\kappa \Re(\lambda_i) + \kappa^2|\lambda_i|^2 \le 0$. Solving for $\kappa$ and recalling the definition of $\eta_i$ in (16), we obtain robust synchronization if

$$ \kappa \in \left[ \frac{c - \sqrt{c(c - \eta_i \sigma)}}{\eta_i \Re(\lambda_i)}, \frac{c + \sqrt{c(c - \eta_i \sigma)}}{\eta_i \Re(\lambda_i)} \right] \quad \forall i = 2, \dots, N, \qquad (A.13) $$

because we simultaneously stabilize all the closed-loop matrices. First, note that from (16c) we have $\eta_i \ge 1$. Moreover, $\eta_i < \infty$, because all eigenvalues $\lambda_i$ have positive real part. Then, since $\sigma > 0$ and $c \ge 1$, for any $i = 2, \dots, N$ it holds that

$$ \frac{c - \sqrt{c(c - \bar{\eta} \sigma)}}{\bar{\eta} \bar{\lambda}} \le \frac{c - \sqrt{c(c - \eta_i \sigma)}}{\eta_i \Re(\lambda_i)} \le \frac{c - \sqrt{c(c - \bar{\eta} \sigma)}}{\bar{\eta} \bar{\lambda}}. $$

Consequently, for any $\kappa \in \mathcal{K}$ as per (18), condition (A.13) holds and (6) holds, as to be proven.

## Appendix B. Proofs of the nonlinear results

### Appendix B.1. Proof of Lemma 3

Consider the candidate Lyapunov function $V: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}_{\ge 0}$ defined as

$$ V(x_1, x_2) := (x_1 - x_2)^\top P (x_1 - x_2), $$

for any two states $x_1, x_2 \in \mathbb{R}^n$. Given any selection of $x_1, x_2$ and $w_1, w_2 \in \mathbb{R}^n$, define function $\tilde{\Phi}: \mathbb{R} \to \mathbb{R}^n$ as

$$ \tilde{\Phi}(s) = \varphi(sx_1 + (1-s)x_2, t) + sw_1 + (1-s)w_2. $$

We have

$$ x_1^+ - x_2^+ = \tilde{\Phi}(1) - \tilde{\Phi}(0) = \varphi(x_1, t) - \varphi(x_2, t) + w_1 - w_2. $$

In view of (20), we obtain

$$ \begin{aligned} V^+ &= V(x_1^+, x_2^+) \\ &= (x_1^+ - x_2^+)^\top P \left[ \int_0^1 \psi(s) ds (x_1 - x_2) + w_1 - w_2 \right] \end{aligned} $$

for some $\psi(s) \in \mathcal{D}\varphi$, for all $s \in [0, 1]$. Then, adding and subtracting $\rho^2 V(x_1, x_2)$ and $V(x_1^+, x_2^+)$ to the right-hand side of the previous equation yields

$$ \begin{aligned} V(x_1^+, x_2^+) - \rho^2 V(x_1, x_2) &= 2(x_1^+ - x_2^+)^\top P \left[ \int_0^1 \psi(s) ds (x_1 - x_2) + w_1 - w_2 \right] \\ &\quad - (x_1^+ - x_2^+)^\top P(x_1^+ - x_2^+) \int_0^1 ds \\ &\quad - \rho^2 (x_1 - x_2)^\top P(x_1 - x_2) \int_0^1 ds \\ &= \int_0^1 \xi^\top \Upsilon(s) \xi ds + 2(x_1^+ - x_2^+)^\top P(w_1 - w_2) \end{aligned} $$

where we defined $\xi := \text{col}(x_1 - x_2, x_1^+ - x_2^+)$ and

$$ \Upsilon(s) := \begin{pmatrix} -\rho^2 P & \psi^\top(s) P \\ P\psi(s) & -P \end{pmatrix}. $$

Applying Schur complement and inequality (22) we obtain $\Upsilon(s) \preceq 0$ for all $s \in [0, 1]$. As a consequence, since $\rho \in (0, 1)$, we obtain

$$ V(x_1^+, x_2^+) - \rho^2 V(x_1, x_2) \le 2(x_1^+ - x_2^+)^\top P(w_1 - w_2). $$

By the generalized Young's inequality and by considering the decomposition $P = \sqrt{P}^\top \sqrt{P}$, we have

$$ \begin{aligned} 2(x_1^+ - x_2^+)^\top P(w_1 - w_2) &\le (1 - \rho) V(x_1^+, x_2^+) + \frac{1}{1 - \rho} (w_1 - w_2)^\top P(w_1 - w_2). \end{aligned} $$

Then, by combining the previous inequalities we obtain

$$ V(x_1^+, x_2^+) - \rho V(x_1, x_2) \le \frac{|P|}{\rho(1-\rho)} |w_1 - w_2|^2. $$

As $\rho \in (0, 1)$ and $P \succ 0$, the function $V$ is a dissipation-form incremental Lyapunov function [34, Definition 7]. Then, the result follows by [34, Theorem 8]. Finally, by using standard arguments (i.e. [35, Theorem 14]) one can conclude the exponential behavior of solutions.

### Appendix B.2. Proof of Proposition 2

For the sake of compactness, let us start by defining

$$ Y = Y^\top := (R + B^\top P B)^{-1}, \quad \Omega = I_n - \kappa B Y B^\top P. $$

Since $B$ is assumed to be full column rank, the matrix $R + B^\top P B$ is invertible, and $Y$ exists. Then, Lemma 3 states that the closed-loop (25), (27), which can be written as (19) with

$$ \varphi(x, t) = f(x, t) + B\alpha(x, t) = \Omega f(x, t), \quad (B.1) $$

is exponentially $\delta$ISS if

$$ J^\top \Omega^\top P \Omega J \preceq \rho^2 P, \quad \forall J \in \mathcal{D}f. \quad (B.2) $$

By expanding the left-hand side in (B.2) and by adding and subtracting $\sigma J^\top P BY B^\top PJ$, due to (B.2) we obtain the equality

$$ \begin{aligned} J^\top \Omega^\top P \Omega J &= J^\top P J - \sigma J^\top P BY B^\top P J \\ &\quad + (\sigma - 2\kappa) J^\top P BY B^\top P J + \kappa^2 J^\top P BY B^\top P BY B^\top P J, \end{aligned} $$

where we note that the first two terms at the right-hand side coincide with $J^\top Q J$.

Note that, since $R \succeq 0$ and $B$ is full rank, by selecting $c$ as in (28), we have $Y^{-1} \succeq c B^\top P B$. Then inequality (26) implies, for all $J \in \mathcal{D}f$

$$ \begin{aligned} J^\top \Omega^\top P \Omega J &\preceq \rho^2 P + \\ &J^\top P BY ((\sigma - 2\kappa)Y^{-1} + \kappa^2 B^\top P B) Y B^\top P J \\ &\preceq \rho^2 P + (c^{-1}\kappa^2 - 2\kappa + \sigma) J^\top P BY B^\top P J. \end{aligned} $$
Since Y is positive definite, (B.2) holds and Lemma 3 applies if $\kappa^2 - 2c\kappa + c\sigma \le 0$, which holds if and only if

$$c - \sqrt{c(c - \sigma)} \le \kappa \le c + \sqrt{c(c - \sigma)},$$

concluding the proof.

### Appendix B.3. Proof of Lemma 4

Since the graph has a directed spanning tree, Lemma 1 ensures that the Laplacian $L$, as in (1), has one zero eigenvalue and $N-1$ eigenvalues with positive real part. Consider again the transformation (A.3) and the Laplacian $\mathcal{L} = TLT^{-1}$ defined as in (A.2). Since $T$ is full rank, by similarity transformation spec $\mathcal{L}$ = spec $L$, namely, it has one zero eigenvalue and $N-1$ eigenvalues with positive real part. Then, due to the block-triangular structure of $\mathcal{L}$, all the eigenvalues of $\tilde{L} = L_{22} - \mathbf{1}L_{12}$ have positive real part. By the Lyapunov equation we conclude that there exists $M \in \mathbb{S}_{>0}^{N-1}$ satisfying $-M\tilde{L} - \tilde{L}^\top M = -I_{N-1} \prec 0$ and a sufficiently small scalar $\mu > 0$ such that (29b) holds.

We now move to the other inequalities in Lemma 4. Since $M \in \mathbb{S}_{>0}^{N-1}$, (29a) and (29c) trivially hold with of $\bar{m}$ and $\bar{m}$ being the smallest and largest eigenvalues of $M$ respectively, $\rho_M = \bar{m}\bar{m}^{-1}$ and a sufficiently large $\bar{\mu}$.

### Appendix B.4. Proof of Theorem 3

Mimicking the linear framework, we show convergence to the synchronization manifold $\mathcal{M}$ by focusing our analysis on the error between agents. If these error dynamics are robustly stable (ISS) with respect to the perturbation signal generated by the errors between the different signals $w_i$, then Problem 1 is solved. Bearing in mind the steps of the proof of Theorem 1, without loss of generality, we define a virtual leader $z = x_1$ and define $N-1$ error coordinates with respect to such a leader node, $\mathbf{e} := \text{col}(e_2, \dots, e_N) \in \mathbb{R}^{Nn}$ with $e_i := x_i - z$ for all $i = 1, \dots, N$. Similarly, we define the incremental disturbance $\tilde{\mathbf{w}} := \text{col}(\tilde{w}_2, \dots, \tilde{w}_N)$ with $\tilde{w}_i := w_i - w_1$. The error dynamics are described, for all $i = 2, \dots, N$, by

$$e_i^+ = f(z+e_i, \mathbf{t}) - f(z, \mathbf{t}) + B \sum_{j=1}^{N} (\ell_{ij} - \ell_{1j}) \alpha(z+e_j, \mathbf{t}) + \tilde{w}_i. \quad (B.3)$$

Since by definition of the Laplacian entries $\sum_{j=1}^N \ell_{ij} = 0$ for any $i$, we can subtract $B \sum_{j=1}^N (\ell_{ij} - \ell_{1j}) \alpha(z, \mathbf{t}) = 0$ from the right-hand side so that (B.3) becomes

$$\begin{align*}
e_i^+ &= \tilde{f}(e_i, \mathbf{t}) + B \sum_{j=2}^{N} \tilde{\ell}_{ij} \tilde{\alpha}(e_j, \mathbf{t}) + \tilde{w}_i \\
\tilde{f}(e_i, \mathbf{t}) &:= f(z+e_i, \mathbf{t}) - f(z, \mathbf{t}) \tag{B.4} \\
\tilde{\alpha}(e_j, \mathbf{t}) &:= \alpha(z+e_j, \mathbf{t}) - \alpha(z, \mathbf{t}) \\
\tilde{\ell}_{ij} &:= \ell_{ij} - \ell_{1j}.
\end{align*}$$

Overall, the closed-loop system can be written in compact form as

$$\mathbf{e}^+ = \varphi(\mathbf{e}, \mathbf{t}) + \tilde{\mathbf{w}}, \quad (B.5)$$

where we defined

$$\varphi(\mathbf{e}, \mathbf{t}) := \begin{pmatrix} \tilde{f}(e_2, \mathbf{t}) + B \sum_{j=2}^{N} \tilde{\ell}_{2j} \tilde{\alpha}(e_j, \mathbf{t}) \\ \vdots \\ \tilde{f}(e_N, \mathbf{t}) + B \sum_{j=2}^{N} \tilde{\ell}_{Nj} \tilde{\alpha}(e_j, \mathbf{t}) \end{pmatrix}. \quad (B.6)$$

Now, select the following candidate Lyapunov function

$$V(\mathbf{e}) = \mathbf{e}^\top (M \otimes P) \mathbf{e}, \quad (B.7)$$

with $M$ defined in Lemma 4. Note that, due to the properties of the Kronecker product, since $M \succ 0$ and $P \succ 0$, $M \otimes P$ is symmetric and positive-definite. Now, for each value of $z = x_1$ and $\mathbf{e} = \text{col}(e_2, \dots, e_N)$, define the function $F_t: \mathbb{R} \to \mathbb{R}^{Nn}$ as

$$F_t(s) := \begin{pmatrix} \tilde{f}_s(s, e_2, \mathbf{t}) + B \sum_{j=2}^{N} \tilde{\ell}_{2j} \tilde{\alpha}_s(s, e_j, \mathbf{t}) \\ \vdots \\ \tilde{f}_s(s, e_N, \mathbf{t}) + B \sum_{j=2}^{N} \tilde{\ell}_{Nj} \tilde{\alpha}_s(s, e_j, \mathbf{t}) \end{pmatrix}, \quad (B.8)$$

parametrized by $t \in \mathbb{N}$, with the definitions

$$\tilde{f}_s(s, e_i, \mathbf{t}) := f(z + s e_i, \mathbf{t}) - f(z, \mathbf{t}), \quad (B.9a)$$

$$\begin{align*}
\tilde{\alpha}_s(s, e_j, \mathbf{t}) &:= \alpha(z + s e_j, \mathbf{t}) - \alpha(z, \mathbf{t}), \quad (B.9b) \\
&= -\kappa Y B^\top P \tilde{f}_s(s, e_i, \mathbf{t}),
\end{align*}$$

where we used (27) and $Y := (R + B^\top PB)^{-1}$. From (B.8)-(B.9) we have $F_t(0) = 0$ and from (20) we get

$$\varphi(\mathbf{e}, \mathbf{t}) = F_t(1) = F_t(1) - F_t(0) = \int_0^1 \partial F(s) ds \mathbf{e}, \quad (B.10)$$

where $\tilde{\psi}_i(s) \in \mathcal{D}f$ for all $i = 2, \dots, N$ and

$$\partial F(s) := [I_{(N-1)n} - \kappa (\tilde{L} \otimes BYB^\top P)] \Psi(s), \quad (B.11)$$

$$\Psi(s) := \operatorname{diag}(\tilde{\psi}_2(s) \dots, \tilde{\psi}_N(s)), \quad (B.12)$$

are obtained from (B.8) by proceeding as in (A.4) and with $\tilde{L}$ defined as in (A.3). Since $V(\mathbf{e}^+) = 2(\mathbf{e}^+)^\top (M \otimes P) \mathbf{e}^+ - V(\mathbf{e}^+)$, subtracting $\rho V(\mathbf{e})$ on both sides and combining (B.5) with (B.10), we obtain

$$\begin{align*}
V(\mathbf{e}^+) - \rho V(\mathbf{e}) &= 2 (\mathbf{e}^+)^\top (M \otimes P) \int_0^1 \partial F(s) ds \mathbf{e} \\
&\quad - [(\mathbf{e}^+)^\top (M \otimes P) \mathbf{e}^+ + \rho \mathbf{e}^\top (M \otimes P) \mathbf{e}] \int_0^1 ds \\
&\quad + 2 (\mathbf{e}^+)^\top (M \otimes P) \tilde{\mathbf{w}}.
\end{align*}$$

Then, by collecting everything under the integral and defining the extended error vector $\xi = \text{col}(\mathbf{e}, \mathbf{e}^+)$ we obtain

$$V(\mathbf{e}^+) - \rho V(\mathbf{e}) = - \int_0^1 \xi^\top \Upsilon(s) \xi ds + 2(\mathbf{e}^+)^\top (M \otimes PB) \tilde{\mathbf{w}}, \quad (B.13a)$$
$$
\Upsilon(s) := \begin{pmatrix} \rho(M \otimes P) & -\partial F^\top(s)(M \otimes P) \\ -(M \otimes P)\partial F(s) & (M \otimes P) \end{pmatrix}. \tag{B.13b}
$$

Since $P \succ 0$, $M \succ 0$, and $M \otimes P$ is invertible, we can study
the positive definiteness of $\Upsilon(s)$ via its Schur complement

$$
\hat{\Upsilon}(s) = \rho(M \otimes P) - \partial F^{\top}(s)(M \otimes P)\partial F(s).
$$

By using the definition of $\partial F(s)$ in (B.11) and the prop-
erties of Kronecker products, we obtain

$$
\hat{\Upsilon}(s) = \rho T_a - \Psi(s)^{\top} (T_a + T_b + T_b^{\top} + T_c) \Psi(s), \quad (B.14)
$$

where we defined $T_a := M \otimes P$ and

$$
\begin{align*}
T_b &:= -\kappa(M\tilde{L}) \otimes (PBYB^\top P), \\
T_c &:= \kappa^2(\tilde{L}^\top M\tilde{L}) \otimes (PBYB^\top PBYB^\top P).
\end{align*}
$$

Noting that $PBYB^\top P$ is symmetric and applying Lemma 4 and using the properties of the Kronecker product, we obtain the following inequality

$$
2\operatorname{Sym}\{T_b\} \preceq -2\kappa\underline{\mu}(M \otimes PBYB^\top P). \quad (B.15)
$$

Similarly, by exploiting the Kronecker product and by us-
ing again Lemma 4, we get

$$
T_c \preceq \kappa^2 \bar{\mu}^2 (M \otimes PBYB^\top PBYB^\top P). \quad (B.16)
$$

Note that, since $R \succeq 0$ and $B$ is full rank, there always
exists a scalar $c \ge 1$ such that $Y^{-1} \succeq cB^\top PB$, specifically
$c$ as in (28). Hence, bound (B.16) leads to

$$
T_c \preceq \kappa^2 \bar{\mu}^2 c^{-1} (M \otimes PBYB^\top P). \quad (B.17)
$$

Using (B.15) and (B.17), matrix $\hat{\Upsilon}$ in (B.14) can be
bounded as

$$
\begin{align}
\hat{\Upsilon}(s) &\succeq \rho(M \otimes P) - \Psi(s)^{\top} (M \otimes \bar{P}) \Psi(s) \tag{B.18} \\
\bar{P} &= P + (c^{-1}\kappa^2\bar{\mu}^2 - 2\kappa\underline{\mu}) PBYB^{\top} P \nonumber
\end{align}
$$

Now, consider $\bar{P}$. By addition and subtraction, it can be
rewritten as

$$
\bar{P} = P - \sigma PBYB^\top P + (\kappa^2 \bar{\mu}^2 c^{-1} - 2\kappa \underline{\mu} + \sigma) PBYB^\top P.
$$

Then, proceeding as in the proof of Theorem 2 (see equa-
tion (A.13)), if $\kappa^2 \bar{\mu}^2 - 2c\kappa \underline{\mu} + c\sigma \le 0$, namely if

$$
\frac{\mu}{\bar{\mu}^2} \left( c - \sqrt{c^2 - \frac{\bar{\mu}^2}{\mu^2} \sigma c} \right) \le \kappa \le \frac{\mu}{\bar{\mu}^2} \left( c + \sqrt{c^2 - \frac{\bar{\mu}^2}{\mu^2} \sigma c} \right),
$$

which holds due to the selection in (31), we obtain

$$
\bar{P} \preceq P - \sigma PBYB^\top P = Q, \quad (B.19)
$$

with $Q$ defined in (26b). Using (29a) from Lemma 4 and
(B.19), $\hat{\Upsilon}(s)$ in (B.18) satisfies

$$
\hat{\Upsilon}(s) \succeq \underline{m}\rho(\mathrm{I}_{N-1} \otimes P) - \bar{m}\Psi(s)^{\top}(\mathrm{I}_{N-1} \otimes Q)\Psi(s).
$$

Recalling from (B.12) the block-diagonal structure of $\Psi(s)$
and exploiting (26), we obtain

$$
\begin{align*}
\hat{\Upsilon}(s) &\succeq \underline{m}\rho(\mathrm{I}_{N-1} \otimes P) - \bar{m} \operatorname{diag}(\{\tilde{\psi}_i(s)^\top Q \tilde{\psi}_i(s)\}_{i=2}^N) \\
&\succeq \underline{m}\rho(\mathrm{I}_{N-1} \otimes P) - \bar{m} \operatorname{diag}(\{\rho^2 P\}_{i=2}^N) \\
&= \underline{m}\rho(\mathrm{I}_{N-1} \otimes P) - \bar{m} \rho^2(\mathrm{I}_{N-1} \otimes P) \\
&\succeq \rho(\underline{m} - \bar{m} \rho)(\mathrm{I}_{N-1} \otimes P) \\
&\succeq \rho(\underline{m} - \bar{m} \rho_M)(\mathrm{I}_{N-1} \otimes P) \succeq 0,
\end{align*}
$$

where we used $0 < \rho \le \rho_M = \underline{m}\bar{m}^{-1}$. Since $\hat{\Upsilon}(s) \succeq 0$ for each $s \in [0, 1]$, we conclude that also $\Upsilon(s)$ defined in (B.13b) satisfies $\Upsilon(s) \succeq 0$ for all $s \in [0, 1]$, and (B.13a) implies

$$
V(\mathbf{e}^{+}) - \rho V(\mathbf{e}) \leq 2(\mathbf{e}^{+})^{\top} (M \otimes P) \tilde{\mathbf{w}} \quad (B.20)
$$

By the generalized Young's inequality and by consider-
ing the factorization $M \otimes P = \sqrt{M \otimes P}^{\top} \sqrt{M \otimes P} =$
$(\sqrt{M \otimes P})^2$ (with $\sqrt{M \otimes P}$ denoting the unique positive
square root of $M \otimes P \succ 0$), we have

$$
\begin{align*}
2(\mathbf{e}^{+})^{\top} (M \otimes P) \tilde{\mathbf{w}} &= 2(\mathbf{e}^{+})^{\top} (\sqrt{M \otimes P})^{2} \tilde{\mathbf{w}} \\
&\leq (1 - \sqrt{\rho}) V(\mathbf{e}^{+}) + \frac{1}{1 - \sqrt{\rho}} \tilde{\mathbf{w}}^{\top} (M \otimes P) \tilde{\mathbf{w}}.
\end{align*}
$$

Then, since $\rho \in (0, 1)$, inequality (B.20) implies

$$
V(\mathbf{e}^{+}) - \sqrt{\rho}V(\mathbf{e}) \leq \frac{1}{\sqrt{\rho}(1 - \sqrt{\rho})} \tilde{\mathbf{w}}^{\top}(M \otimes P)\tilde{\mathbf{w}},
$$

thus proving exponential ISS properties of the e dynamics
due to the quadratic form of (B.7). Finally, similarly to
the linear scenario of Theorem 1, relations (A.7) hold and
robust synchronization as in Problem 1 is obtained, thus
concluding the proof.

Appendix B.5. *Proof of Proposition 3*

Consider the last LMI of (32). Since $W \succ 0$ and $U \succ 0$,
then its (2, 2) entry is positive definite. Then, by the Schur
complement, (32b) holds if and only if

$$
W + \sigma B U B^{\top} - \Sigma \succeq \sigma(1 - \sigma) B U B^{\top} (W + (1 - \sigma) B U B^{\top})^{-1} B U B^{\top},
$$

which can be rearranged as

$$
\Sigma \preceq W + \sigma B \left( U - (1 - \sigma)UB^{\top} \times (W + (1 - \sigma)BUB^{\top})^{-1}BU \right) B^{\top}. \quad (B.21)
$$

By applying the identity (11) with $\mathcal{A} = U^{-1}$, $\mathcal{B} = \sqrt{1 - \sigma B^{\top}}$, $\mathcal{C} = \sqrt{1 - \sigma B}$, $\mathcal{D} = W$, inequality (B.21) is equivalent to

$$
\Sigma \preceq W + \sigma B (U^{-1} + (1 - \sigma) B^{\top} W^{-1} B)^{-1} B^{\top}.
$$
Left and right multiplication of both sides by $W^{-1} \succ 0$ yields the equivalent condition

$$W^{-1}\Sigma W^{-1} \preceq W^{-1} - \sigma W^{-1}B(-U^{-1} + B^{\top}W^{-1}B) + \sigma B^{\top}W^{-1}B)^{-1}B^{\top}W^{-1}. \quad (B.22)$$

Once again, by applying identity (11) with $A = W$, $B = \sqrt{\sigma}B$, $C = \sqrt{\sigma}B^{\top}$, $D = -(U^{-1} + B^{\top}W^{-1}B)$, inequality (B.22) is equivalent to

$$W^{-1}\Sigma W^{-1} \preceq (W - \sigma B(U^{-1} + B^{\top}W^{-1}B)^{-1}B^{\top})^{-1}. \quad (B.23)$$

By left and right multiplying both sides by $W$, we obtain the equivalent inequality

$$\Sigma \preceq W(W - \sigma B(U^{-1} + B^{\top}W^{-1}B)^{-1}B^{\top})^{-1}W. \quad (B.23)$$

Since $(ABC)^{-1} = C^{-1}B^{-1}A^{-1}$ for any invertible matrices $A, B, C$, inequality (B.23) is equivalent to

$$\Sigma \preceq (W^{-1}(W - \sigma B(U^{-1} + B^{\top}W^{-1}B)^{-1}B^{\top})W^{-1})^{-1} \preceq (P - \sigma PB(R + B^{\top}PB)^{-1}B^{\top}P)^{-1} = Q^{-1}, \quad (B.24)$$

where we used $P = W^{-1}$, $R = U^{-1}$ and the definition of $Q$ in (26b). Consider now the right LMI in (32a). By left and right multiplication by the matrix

$$\tilde{T} = \begin{pmatrix} W^{-1} & 0 \\ 0 & I_n \end{pmatrix} = \begin{pmatrix} P & 0 \\ 0 & I_n \end{pmatrix},$$

we have

$$\begin{pmatrix} \rho P & J^{\top} \\ J & \rho \Sigma \end{pmatrix} = \begin{pmatrix} \rho W^{-1} & J^{\top} \\ J & \rho \Sigma \end{pmatrix} \succeq 0.$$

Then, by the Schur complement, we obtain the equivalent condition $\rho W^{-1} - \rho^{-1}J^{\top}\Sigma^{-1}J \succeq 0$, which, with the selection $P = W^{-1}$, can be written as

$$J^{\top}\Sigma^{-1}J \preceq \rho^2 P. \quad (B.25)$$

Summarizing, we proved the equivalence of (32) with the four inequalities $P \succ 0$, $\Sigma \succ 0$, (B.24) and (B.25), where we emphasize that, under (32), $Q^{-1}$ exists due to the positive definiteness of $Q$ implied by $R \succ 0$, as established in Remark 6. More specifically, (32) is equivalent to

$$P \succ 0, \quad \Sigma \succ 0, \quad \Sigma^{-1} \succeq Q, \quad J^{\top}\Sigma^{-1}J \preceq \rho^2 P. \quad (B.26)$$

To complete the first part of the proof, we show that (B.26) is equivalent to (26). If (26) holds, then $Q \succ 0$ and (B.26) holds with $\Sigma = Q^{-1}$. If (B.26) holds, then

$$J^{\top}QJ \preceq J^{\top}\Sigma^{-1}J \preceq \rho^2 P,$$

thus completing the first part of the proof. To prove that (32) is a GEVP in $(\sigma, \rho, R)$, let us denote by $\bar{W}, \bar{\Sigma}$ the solution of (32) with $(\bar{\sigma}, \bar{\rho}, \bar{R})$. If $\sigma \ge \bar{\sigma} > 0$, $0 \prec R \prec \bar{R}$, we obtain

$$\bar{Q} = \bar{P} - \bar{\sigma}\bar{P}B(\bar{R} + B^{\top}\bar{P}B)^{-1}B^{\top}\bar{P} \preceq \bar{P} - \sigma\bar{P}B(R + B^{\top}\bar{P}B)^{-1}B^{\top}\bar{P} = Q,$$

where $\bar{P} = \bar{W}^{-1} \succ 0$. Then, by (B.26) and since $\bar{\Sigma}^{-1} \succeq \bar{Q}$ and $\rho \ge \bar{\rho}$, the following inequalities hold

$$\bar{P} \succ 0, \quad \bar{\Sigma} \succ 0, \quad \bar{\Sigma}^{-1} \succeq Q, \quad J^{\top}\bar{\Sigma}^{-1}J \preceq \bar{\rho}^2\bar{P}.$$

Due to the equivalence between inequalities (B.26) and (32), we conclude that $P = \bar{P}$ is solution to (32) with $(\sigma, \rho, R)$. Similar reasoning proves infeasibility of (B.26) for any $(\sigma, \rho, R)$ such that $\sigma \le \bar{\sigma}$, $\rho \le \bar{\rho}$, $R \succeq \bar{R}$ if (B.26) is infeasible for $(\sigma, \rho, R)$.

## Appendix B.6. Proof of Proposition 4

Due to [5, Theorem 2.3.7], (35) holds for any diagonal $S \in \mathbb{S}_{\ge 0}^p$ and any $J_{\phi} \in \mathcal{D}\phi$. Define $\Lambda^{\top} := \begin{pmatrix} I_n & 0 \\ 0 & I_n F J_{\phi} \end{pmatrix}$ with any diagonal $J_{\phi} \in \mathcal{D}\phi$. Equation (35) implies

$$\Lambda^{\top} \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -F\Omega^{\top}S \\ 0 & -S\Omega F^{\top} & 2S \end{pmatrix} \Lambda = \Lambda^{\top}\Pi\Lambda \preceq 0. \quad (B.27)$$

Consider now (36), which implies

$$\Lambda^{\top} \begin{pmatrix} \rho W & WA^{\top} & WC^{\top} \\ AW & \rho \Sigma & -F \Omega^{\top} S \\ CW & -S \Omega F^{\top} & 2S \end{pmatrix} \Lambda = \Lambda^{\top}(\Xi + \Pi)\Lambda \succeq 0,$$

where we defined

$$\Xi := \begin{pmatrix} \rho W & WA^{\top} & WC^{\top} \\ AW & \rho \Sigma & 0 \\ CW & 0 & 0 \end{pmatrix}. \quad (B.28)$$

By (B.27), we have $\Lambda^{\top}\Xi\Lambda \succeq \Lambda^{\top}(\Xi + \Pi)\Lambda \succeq 0$, thus showing $\Lambda^{\top}\Xi\Lambda \succeq 0$. Then, the expansion of the product leads to

$$\Lambda^{\top}\Xi\Lambda = \begin{pmatrix} \rho W & W(A + FJ_{\phi}C)^{\top} \\ (A + FJ_{\phi}C)W & \rho \Sigma \end{pmatrix} \succeq 0,$$

for all $J_{\phi} \in \mathcal{D}\phi$. By Proposition 3, the assumption (32b) implies that conditions (26) hold with $P = W^{-1} \succ 0$ for any $J \in \mathcal{D}f = A + F\mathcal{D}\phi C$. The proof is concluded by Proposition 2.

## References

[1] F. Bayer, M. Bürger, and F. Allgöwer. Discrete-time incremental ISS: A framework for robust NMPC. In *European Control Conference*, pages 2068–2073, 2013.

[2] F. Bonassi, M. Farina, and R. Scattolini. Stability of discrete-time feed-forward neural networks in NARX configuration. *IFAC-PapersOnLine*, 54(7):547–552, 2021.

[3] F. Cacace, M. Mattioni, S. Monaco, and D. Normand-Cyrot. Consensus and multi-consensus for discrete-time LTI systems. *Automatica*, 166:111718, 2024.

[4] M. Z. Q. Chen, L. Zhang, H. Su, and G. Chen. Stabilizing solution and parameter dependence of modified algebraic Riccati equation with application to discrete-time network synchronization. *IEEE Transactions on Automatic Control*, 61(1):228–233, 2016.

[5] F. H. Clarke. *Optimization and nonsmooth analysis*. SIAM, 1990.
[6] A. Cristofaro and M. Mattioni. Hybrid consensus for multi-agent systems with time-driven jumps. *Nonlinear Analysis: Hybrid Systems*, 43:101113, 2021.

[7] W. D'Amico, A. La Bella, and M. Farina. An incremental input-to-state stability condition for a generic class of recurrent neural networks. *IEEE Transactions on Automatic Control*, 2024.

[8] M. I. El-Hawwary and M. Maggiore. Reduction theorems for stability of closed sets with application to backstepping control design. *Automatica*, 49(1):214–222, 2013.

[9] V. Fromion and G. Scorletti. The behaviour of incrementally stable discrete time systems. *American Control Conference*, 6:4563–4567, 1999.

[10] H. Fu, X. Chen, W. Wang, and M. Wu. Data-based optimal synchronization control for discrete-time nonlinear heterogeneous multiagent systems. *IEEE Transactions on Cybernetics*, 2020.

[11] M. Giaccagli, V. Andrieu, S. Tarbouriech, and D. Astolfi. Infinite gain margin, contraction and optimality: an LMI-based design. *European Journal of Control*, 68:100685, 2022.

[12] M. E. Gilmore, C. Guiver, and H. Logemann. Semi-global incremental input-to-state stability of discrete-time Lur'e systems. *Systems & Control Letters*, 136:104593, 2020.

[13] K. Hengster-Movric, K. You, F. Lewis, and L. Xie. Synchronization of discrete-time multi-agent systems on graphs using riccati design. *Automatica*, 49(2):414–423, 2013.

[14] H. J. C. Huijberts, T. Lilge, and H. Nijmeijer. Nonlinear discrete-time synchronization via extended observers. *International Journal of Bifurcation and Chaos*, 11(07):1997–2006, 2001.

[15] A. Isidori. Coordination and consensus of linear systems. In *Lectures in Feedback Design for Multivariable Systems*, pages 135–163. Springer, 2017.

[16] M. Jungers, M. F. Shakib, and N. van De Wouw. Discrete-time convergent nonlinear systems. *IEEE Transactions on Automatic Control*, 2024.

[17] M. Jungers and N. van de Wouw. Discrete-Time Nonlinear Control Systems. In *Reference Module in Materials Science and Materials Engineering*. Elsevier, 2025.

[18] Y. Kawano and Y. Hosoe. Contraction analysis of discrete-time stochastic systems. *IEEE Transactions on Automatic Control*, pages 1–16, 2023.

[19] Y. Kawano and K. Kashima. An LMI framework for contraction-based nonlinear control design by derivatives of Gaussian process regression. *Automatica*, 151:110928, 2023.

[20] S. Knorn, Z. Chen, and R. H. Middleton. Overview: Collective control of multiagent systems. *IEEE Transactions on Control of Network Systems*, 3(4):334–347, 2015.

[21] Z. Li and Z. Duan. *Cooperative Control of Multi-Agent Systems: A Consensus Region Approach*. CRC Press, 2015.

[22] Z. Li, G. Wen, Z. Duan, and W. Ren. Designing fully distributed consensus protocols for linear multi-agent systems with directed graphs. *IEEE Transactions on Automatic Control*, 60(4):1152–1157, 2014.

[23] T. Lilge. Nonlinear discrete-time observers for synchronization problems. In *New Directions in nonlinear observer design*, pages 491–510. Springer, 1999.

[24] A. Moreschini, M. Bin, A. Astolfi, and T. Parisini. A generalized passivity theory over abstract time domains. *IEEE Transactions on Automatic Control*, 2024.

[25] W. Ren and R. W. Beard. Consensus seeking in multiagent systems under dynamically changing interaction topologies. *IEEE Transactions on automatic control*, 50(5):655–661, 2005.

[26] A. Saberi, A. A. Stoorvogel, M. Zhang, and P. Sannuti. Synchronization of Multi-Agent Systems in the Presence of Disturbances and Delays. Springer Nature, 2022.

[27] L. Scardovi and R. Sepulchre. Synchronization in networks of identical linear systems. *Automatica*, 45(11):2557–2562, 2009.

[28] I. Schimperna and L. Magni. On incremental input-to-state stability of interconnected discrete-time systems. *Automatica*, 179:112406, 2025.

[29] R. Sepulchre, M. Jankovic, and P. V. Kokotovic. *Constructive nonlinear control*. Springer Science & Business Media, 2012.

[30] D. D. Šiljak and D. M. Stipanovic. Robust stabilization of nonlinear systems: The LMI approach. *Mathematical problems in Engineering*, 6(5):461–493, 2000.

[31] B. Sinopoli, L. Schenato, M. Franceschetti, K. Poolla, M. I. Jordan, and S. S. Sastry. Kalman filtering with intermittent observations. *IEEE Transactions on Automatic Control*, 49(9):1453–1464, 2004.

[32] A. A. Stoorvogel, A. Saberi, M. Zhang, and Z. Liu. Solvability conditions and design for synchronization of discrete-time multiagent systems. *International Journal of Robust and Nonlinear Control*, 28(4):1381–1401, 2017.

[33] H. Y. Sutrisno, S. Trenn, and B. Jayawardhana. Nonlinear singular switched systems in discrete-time: Solution theory and incremental stability under restricted switching signals. In *62nd IEEE Conference on Decision and Control*, pages 914–919, 2023.

[34] D. N. Tran, B. S. Rüffer, and C. M. Kellett. Incremental stability properties for discrete-time systems. In *IEEE 55th Conference on Decision and Control*, pages 477–482, 2016.

[35] D. N. Tran, B. S. Rüffer, and C. M. Kellett. Convergence properties for discrete-time nonlinear systems. *IEEE Transactions on Automatic Control*, 64(8):3415–3422, 2018.

[36] F. J. Vargas and R. A. González. On the existence of a stabilizing solution of modified algebraic Riccati equations in terms of standard algebraic Riccati equations and linear matrix inequalities. *IEEE Control Systems Letters*, 4(1):91–96, 2019.

[37] L. Wei, R. McCloy, and J. Bao. Control contraction metric synthesis for discrete-time nonlinear systems. *IFAC-PapersOnLine*, 54(3):661–666, 2021.

[38] L. Wei, R. McCloy, and J. Bao. Discrete-time contraction-based control of nonlinear systems with parametric uncertainties using neural networks. *Computers & Chemical Engineering*, 166:107962, 2022.

[39] H. Yin, P. Seiler, and M. Arcak. Stability analysis using quadratic constraints for systems with neural network controllers. *IEEE Transactions on Automatic Control*, 67(4):1980–1987, 2022.

[40] K. You and L. Xie. Network topology and communication data rate for consensusability of discrete-time multi-agent systems. *IEEE Transactions on Automatic Control*, 56(10):2262–2275, 2011.

[41] N. Zaupa, G. Giordano, I. Queinnec, S. Tarbouriech, and L. Zaccarian. Equivalent conditions for the synchronization of identical linear systems over arbitrary interconnections. *European Journal of Control*, page 101099, 2024.

[42] L. Zhang, M. Z. Q. Chen, Y. Zou, and G. Chen. Gain margin and Lyapunov analysis of discrete-time network synchronization via riccati design. *International Journal of Robust and Nonlinear Control*, 2023.

[43] Z. Zuo, J. Wang, and L. Huang. Robust stabilization for nonlinear discrete-time systems. *International Journal of Control*, 77(4):384–388, 2004.

<figure><img src="image_5.png" /></figure>

**Samuele Zoboli** obtained his Ph.D. in Control Theory in 2023 from University of Lyon 1, France, and since 2025 he is a CNRS researcher at LAAS-CNRS, Toulouse, France. His research interests include stabilization of discrete-time nonlinear systems, multi-agent systems, reinforcement learning and control-applied artificial intelligence.

<figure><img src="image_6.png" /></figure>

**Daniele Astolfi** obtained a joint Ph.D. degree in Control Theory from the University of Bologna, Italy, and from Mines ParisTech, France, in 2016. Since 2018, he is a CNRS Researcher at LAGEPP, Lyon, France. His research interests include observer design, feedback stabilization and output regulation for nonlinear systems.
systems, networked control systems, hy-
brid systems, and multi-agent systems. He serves as an as-
sociate editor of the IFAC journal Automatica since 2018 and
European Journal of Control since 2023.

<figure><img src="image_7.png" /></figure>

**Vincent Andrieu** is a Senior Researcher (Directeur de Recherche) at the French National Center for Scientific Research (CNRS). He is currently affiliated with the LAGEPP laboratory at the University of Lyon 1, which he joined in 2010. He received his Ph.D. in control theory from Ecole des Mines de Paris (Mines Paris - PSL) in 2005 and a graduate degree in applied mathematics from INSA de Rouen in 2001. His main research interests are in the feedback stabilization of controlled dynamical nonlinear systems and state estimation problems. Since 2018 he is an associate editor of the IEEE Transactions on Automatic Control, and senior editor for System & Control Letters.

<figure><img src="image_8.png" /></figure>

**Giacomo Casadei** is an associate professor in Université Grenoble Alpes (IUT) and a member of team DANCE in GIPSA-Lab since 2024. He received the Master degree in Automation in 2012 and the PhD degree in Automation and Operative Research in 2016 at the University of Bologna, Italy. From 2016 to 2018, he has been postdoc in the NeCS team, a joint team of GIPSA-Lab (CNRS) and INRIA. From 2018 to 2024, he has held an associate professor position in Ecole Centrale de Lyon. His research interests include modeling and control of networks, with a particular focus on synchronization of nonlinear systems and the use of nonlinear control techniques in networks.

<figure><img src="image_9.png" /></figure>

**Luca Zaccarian** has been an Assistant Professor and then an Associate Professor with the University of Roma Tor Vergata, Rome, Italy, since 2000. In 2011, he became Directeur de Recherche at the Laboratoire d'Analyse et d'Architecture of Systems-CNRS, Toulouse, France. Since 2013, he is also a part-time Professor with the University of Trento, Italy. He is senior editor of Automatica.