# Robust synchronization via maximal monotone couplings

Felix Miranda-Villatoro

► To cite this version:

Felix Miranda-Villatoro. Robust synchronization via maximal monotone couplings. ECC 2023 - European Control Conference, Jun 2023, Bucarest (Romania), Romania. pp.1-8. hal-03799388v2

HAL Id: hal-03799388

https://hal.science/hal-03799388v2

Submitted on 25 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
# Robust synchronization via set-valued maximal monotone couplings

Félix A. Miranda-Villatoro<fnref n="1" />

*Abstract*— We explore the use of set-valued coupling laws for the design of robust synchronized behaviors in networks of dynamical systems. Under an incremental dissipativity context, it is shown that coupling systems via maximal monotone mappings leads to synchronization that is robust against matched disturbances. Additionally, it is shown that perfect synchronization of heterogeneous networks with persistent matched disturbances is attained with finite coupling strength but infinite incremental gain on the coupling maps. The real-life implementation of the proposed controllers is studied under the context of practical synchronization via Yosida regularizations. Simulations illustrate the effectiveness of the proposed methods.

## I. INTRODUCTION

The study of interacting dynamical systems reaching a common uniform behavior has received numerous attention from the control community in recent years. The problem has been broadly addressed under several scenarios and under different conditions, see e.g., [1], [2], [3], [4], [5], [6] and references therein. Robust synchronization in the presence of heterogeneities, modeled as disturbances, occupies a special place of interest due to its broad application in several fields. Among the studies dealing with robust synchronization one finds [7], [8], [9]. For instance, in [8] the authors studied practical synchronization of perturbed agents coupled via linear diffusion. It is shown there that the ultimate bound of the synchronization error is inversely proportional to the coupling gain, relying on high-gain couplings as a way of enhancing precision. Nonsmooth coupling laws were analyzed in [10], [7], where sufficient conditions for practical synchronization were provided. Though it is left open if there is any advantage on using nonsmooth coupling schemes, since Assumption 3 in [7] also holds in the case of smooth coupling maps.

The main aim of this paper is to investigate the advantages of set-valued coupling, as well as, issues concerning its implementation. The second goal is to show the important role of monotonicity as unifying framework for robust synchronization strategies such as set-valued sliding-mode control. To this end, we analyze the robust synchronization of systems against matched persistent disturbances, where each agent is described by a Lur'e-type system, whereas the coupling strategy is characterized by a static set-valued map. First, theoretical results are presented regarding perfect synchronization in presence of matched disturbances. It is shown that perfect synchronization is achievable with a special family of *dry friction set-valued couplings*.

Afterwards we turn our attention towards implementation issues. Even though some nonsmooth (discontinuous) synchronization strategies have been recently proposed in [7], [10], there are still open issues concerning their implementation. For instance, there are many applications where discontinuous controllers are impermissible because of their intrinsic connection with chattering [11]. For set-valued controllers, the problem translates into developing appropriate selection strategies that will compensate for the disturbances. In this regard, inspired by analog electrical circuits, we change the static set-valued couplings by dynamic ones described by linear complementarity systems. It is shown that the latter achieves synchronization with arbitrary precision (practical synchronization) and allows for selection strategies that are independent of the parameters of the agents.

This paper is organized as follows. Preliminaries are covered in Section II, whereas Section III presents the problem to deal with and related results on well-posedness of dynamical systems with set-valued feedback laws. Section IV deals with the perfect asymptotic synchronization of perturbed systems via the design of set-valued maximal monotone coupling laws, and the effects of regularization are analyzed in Section V where numerical results are presented for an heterogeneous network of FitzHugh-Nagumo oscillators. Finally, conclusions take place at the end of the paper.

### Notation

Let $\langle \cdot, \cdot \rangle$ denote the Euclidean inner product in $\mathbb{R}^l$ and $\| \cdot \|$ the corresponding Euclidean norm. The set $B_r^l := \{q \in \mathbb{R}^l \| \|q\| \le r\}$ denotes the closed ball with radius $r$ and center at zero in $\mathbb{R}^l$. When the dimension of the space is clear from the context we will denote the closed ball simply as $B_r$. The interior of a set $S \subset \mathbb{R}^l$ is denoted as int $S$, whereas the relative interior is denoted as rint $S$. dist($s, S$) denotes the conventional distance function from a point $s$ to a closed convex set $S$. The matrix $I_l$ denotes the identity matrix in $\mathbb{R}^{l \times l}$. For a matrix $A \in \mathbb{R}^{l \times l}$, $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ denote the minimum and maximum eigenvalue of $A$, respectively. The vector $1_N$ denotes the vector of ones in $\mathbb{R}^N$.

## II. PRELIMINARIES

### A. Elements from graph theory

The interconnection structure of a network of $N$ agents is usually modeled by a graph $G(V, E)$, where $V = \{v_1, \dots, v_N\}$ is the set of vertices, $v_i \in V$ represents the $i$-th agent, $E \subset V \times V$ is the set of edges. If $\{v_i, v_j\} \in E$, then the vertices $v_i$ and $v_j$ are *adjacent*. In addition, each edge $\{v_i, v_j\} \in E$ is *incident* with the vertices $v_i$ and $v_j$. When the set of vertices and edges is clear from the context

<footnote n="1">¹Félix A. Miranda-Villatoro is with the Université Grenoble-Alpes, LJK, and INRIA Grenoble-Rhône-Alpes Research Center, 38000 Grenoble, France. felix.miranda-villatoro@inria.fr</footnote>
we will denote the graph simply as $\mathcal{G}$. Note that, in the set notation used, $\{v_i, v_j\} = \{v_j, v_i\}$, that is, the graph under consideration is *undirected*. Let $v_i, v_j \in \mathcal{V}$, a path of length $m$ from $v_i$ to $v_j$ is a sequence of $m+1$ vertices $\{v_k\}_{k=0}^m \subseteq \mathcal{V}$ such that $\nu_0 = v_i$, $\nu_m = v_j$, and $\{\nu_k, \nu_{k+1}\} \in \mathcal{E}$ for $k \in \{0, \dots, m-1\}$. The graph $\mathcal{G}(\mathcal{V}, \mathcal{E})$ is *connected* if for any pair of distinct vertices $(v_i, v_j) \in \mathcal{V} \times \mathcal{V}$ there exists a path from $v_i$ to $v_j$. In what follows we set $N = |\mathcal{V}|$ and $E = |\mathcal{E}|$, the number of vertices and edges in the network, respectively. A subgraph $\mathcal{G}'(\mathcal{V}', \mathcal{E}')$ of $\mathcal{G}(\mathcal{V}, \mathcal{E})$, denoted as $\mathcal{G}' \subseteq \mathcal{G}$, is any graph such that $\mathcal{V}' \subseteq \mathcal{V}$ and $\mathcal{E}' \subseteq \mathcal{E}$. A subgraph of $\mathcal{G}$ is called *spanning subgraph* if $\mathcal{V}' = \mathcal{V}$ and $\mathcal{E}' \subseteq \mathcal{E}$. For each edge $\epsilon_k = \{v_i, v_j\} \in \mathcal{E}$ a sign to each end of $\epsilon_k$ is assigned. Such sign assignation will provide an *orientation* to the graph $\mathcal{G}(\mathcal{V}, \mathcal{E})$. Along all the manuscript, it is assumed that an orientation has been chosen and it is fixed. Thus, the *oriented incidence matrix* $\Theta \in \mathbb{R}^{N \times E}$ is given as, see e.g., [12],

$$ [\Theta]_{i,k} = \begin{cases} +1, & \text{if } v_i \text{ is the positive end of } \epsilon_k; \\ -1, & \text{if } v_i \text{ is the negative end of } \epsilon_k; \\ 0, & \text{otherwise.} \end{cases} $$

If the graph $\mathcal{G}$ has $N$ vertices and $c$ connected components, then any associated incidence matrix has rank $N-c$. Moreover, the graph Laplacian $L = \Theta\Theta^\top$ and $\Theta^\top\mathbf{1}_N = 0$.

### B. Maximal monotone maps

A set-valued map<fnref n="1" /> $\mathbf{M}: \mathbb{R}^l \rightrightarrows \mathbb{R}^l$ maps points from $\mathbb{R}^l$ to subsets of $\mathbb{R}^l$. The set $\text{dom}\mathbf{M} := \{\eta \in \mathbb{R}^l | \mathbf{M}(\eta) \neq \emptyset\} \subseteq \mathbb{R}^l$ denotes the domain of $\mathbf{M}$; the set $\text{gph}(\mathbf{M}) := \{(\eta, \vartheta) \in \mathbb{R}^l \times \mathbb{R}^l | \vartheta \in \mathbf{M}(\eta)\}$ denotes the graph of $\mathbf{M}$; and the set $\text{rge}\mathbf{M} := \{\vartheta \in \mathbb{R}^l | \exists \eta \in \mathbb{R}^l \text{ such that } \vartheta \in \mathbf{M}(\eta)\}$ denotes the range of $\mathbf{M}$.

A set-valued map $\mathbf{M}$ is monotone if for any two pairs $(\eta_1, \vartheta_1), (\eta_2, \vartheta_2) \in \text{gph}(\mathbf{M})$,

$$ \langle \eta_1 - \eta_2, \vartheta_1 - \vartheta_2 \rangle \ge 0. $$

In addition, $\mathbf{M}$ is called maximal monotone if it is monotone and its graph is not strictly contained in the graph of any other monotone map. For instance, maximal monotone maps comprise the subdifferentials, $\partial f$, of convex, proper, lower semicontinuous functions [13],

$$ \partial f(\eta) = \{\vartheta \in X | \langle \vartheta, \zeta - \eta \rangle \le f(\zeta) - f(\eta), \text{ for all } \zeta \in \mathbb{R}^l\}. $$

Special cases include the set-valued sign map **Sgn**: $\mathbb{R} \rightrightarrows [-1, 1]$, such that **Sgn**(0) = [-1, 1] and **Sgn**(η) = {$η/|η|$} for $η \neq 0$ and the normal cone to a convex set $S \in \mathbb{R}^l$ given as, $\mathbb{N}_S(η) = \{\vartheta \in \mathbb{R}^l | \langle \vartheta, \sigma - \eta \rangle \le 0, \text{ for all } \sigma \in S\}$.

For a maximal monotone $\mathbf{M}$: $\text{dom}\mathbf{M} \rightrightarrows \mathbb{R}^l$, the so-called *Yosida approximation* of index $\varepsilon > 0$, $\mathcal{Y}_M^\varepsilon: \mathbb{R}^l \to \mathbb{R}^l$, is the single-valued, Lipschitz continuous function,

$$ \vartheta \mapsto \frac{1}{\varepsilon}(\vartheta - \mathcal{J}_{\varepsilon\mathbf{M}}(\vartheta)), \quad (1) $$

where $\mathcal{J}_{\varepsilon\mathbf{M}}: \mathbb{R}^l \to \mathbb{R}^l$ is the so-called resolvant of $\varepsilon\mathbf{M}$, that is,

$$ \mathcal{J}_{\varepsilon\mathbf{M}} := (I + \varepsilon\mathbf{M})^{-1}, \quad (2) $$

which is single-valued and firmly non-expansive [13]. It follows from (1) and (2) that for any $\varepsilon > 0$, and any $\vartheta \in \mathbb{R}^l$,

$$ \vartheta = \mathcal{J}_{\varepsilon\mathbf{M}}(\vartheta) + \varepsilon \mathcal{Y}_{\mathbf{M}}^{\varepsilon}(\vartheta), \quad (3) $$

and,

$$ \mathcal{Y}_{\mathbf{M}}^{\varepsilon}(\vartheta) \in \mathbf{M}(\mathcal{J}_{\varepsilon\mathbf{M}}(\vartheta)). \quad (4) $$

Moreover, it follows from [13, Lemma 12.14] that

$$ \mathcal{Y}_{\mathbf{M}}^{\varepsilon}(\vartheta) = \mathcal{J}_{\frac{1}{\varepsilon}\mathbf{M}-1}\left(\frac{\vartheta}{\varepsilon}\right). \quad (5) $$

### III. PROBLEM FORMULATION AND WELL-POSEDNESS OF LUR'E NETWORKS WITH MONOTONE COUPLING

We consider an ensemble of $N$ dynamical systems, where the $k$-th agent (denoted as $v_k \in \mathcal{V}$), is a perturbed Lur'e-type system given by

$$ v_k: \begin{cases} \dot{x}_k(t) = Ax_k(t) - B_1\varphi(y_k(t)) \\ +B_2(u_k(t) + \xi(t, x(t))), \\ y_k(t) = C_1x_k(t), \\ w_k(t) = C_2x_k(t). \end{cases}, \quad (6) $$

where $x_k(t) \in \mathbb{R}^n$ is the state of the $k$-th agent; $y_k(t) \in \mathbb{R}^m$ is an (internal) variable used only by the $k$-th agent; $u_k(t), w_k(t) \in \mathbb{R}^p$ are, respectively, the input and output used to exchange information with other agents; and $\xi_k: \mathbb{R} \times \mathbb{R}^n \to \mathbb{R}^m$ is an unknown function accounting for unmodeled dynamic effects, as well as, external disturbances affecting the system. The map $\varphi: \mathbb{R}^m \to \mathbb{R}^m$ is assumed Lipschitz continuous, satisfying extra conditions stated below. Finally, the matrices $A, B_i$, and $C_i$, $i \in \{1, 2\}$, are constant and of the appropriate dimensions. In this work the following standing assumptions are considered.

**Assumption 1:** Each function $\xi_k: \mathbb{R} \times \mathbb{R}^n \to \mathbb{R}^m$ is Lipschitz continuous in the second argument and the map $t \mapsto \xi_k(t, x(t))$ is measurable and uniformly bounded in $\mathcal{L}^\infty[\mathbb{R}_+; \mathbb{R}^m]$, for all $k \in \{1, \dots, N\}$.

**Assumption 2:** The map $\varphi: \mathbb{R}^m \to \mathbb{R}^m$ is such that for any two $\eta, \tilde{\eta} \in \mathbb{R}^m$, the following incremental sector condition holds:

$$ \begin{bmatrix} \Delta\eta \\ \Delta\varphi \end{bmatrix}^\top \begin{bmatrix} K_1^\top K_2 + K_2^\top K_1 & (K_1 + K_2)^\top \\ K_1 + K_2 & 2I_m \end{bmatrix} \begin{bmatrix} \Delta\eta \\ \Delta\varphi \end{bmatrix} \le 0 \quad (7) $$

where $\Delta\eta = \eta - \tilde{\eta}$, $\Delta\varphi = \varphi(\eta) - \varphi(\tilde{\eta})$, $K_1, K_2 \in \mathbb{R}^{m \times m}$ are diagonal matrices such that $0 < K_2 - K_1$.

Functions $\varphi$ satisfying the incremental sector condition (7) include “N”-shape and incrementally passive maps. Thus, the family of Lur'e systems into consideration is general enough to show bistable, oscillatory or chaotic behavior for instance.

<footnote n="1">¹i.e., a non-necessarily single-valued map.</footnote>
The overall dynamics is written in compact form as

$$ \left\{ \begin{array}{l} \dot{x}(t) = (I_N \otimes A)x(t) - (I_N \otimes B_1)\Phi(y(t)) \\ \qquad + (I_N \otimes B_2)(u(t) + \xi(t, x(t))), \\ y(t) = (I_N \otimes C_1)x(t), \\ w(t) = (I_N \otimes C_2)x(t), \end{array} \right. \qquad (8) $$

where $x(t) \in \mathbb{R}^{Nn}$, $y(t) \in \mathbb{R}^{Nm}$, $u(t)$, $w(t)$, $\xi(t, x(t)) \in \mathbb{R}^{Np}$ are the vector variables with components $x(t) = [x_1(t)^\top, \dots, x_N(t)^\top]^\top$. A similar notation holds for the signals $u(t)$, $w(t)$, $y(t)$ and $\xi(t)$. The aggregated map $\Phi: \mathbb{R}^{Nm} \to \mathbb{R}^{Nm}$ is such that $y \mapsto [\varphi(y_1)^\top, \dots, \varphi(y_N)^\top]^\top$.

Let $\mathcal{G}(\mathcal{V}, \mathcal{E})$ be a connected graph encoding the connection structure of the network, and let $\Theta \in \mathbb{R}^{N \times E}$ be an oriented incidence matrix of $\mathcal{G}$. Then $\Theta$ has rank $N-1$ and null $(\Theta^\top) = \text{rge}(\mathbf{1}_N)$. Intuitively, asymptotic synchronization is achieved when there is a coupling law $u(t)$ such that

$$ \text{dist}(x(t); \text{rge}(\mathbf{1}_N \otimes I_n)) \to 0 \text{ as } t \to \infty, $$

that is, the mismatch $x_i(t) - x_j(t) \to 0$ as $t \to \infty$, for all $i, j \in \{1, \dots, N\}$. Equivalently, the network (8) achieves global (semi-global, practical) asymptotic synchronization if and only if the zero solution associated to the dynamics of the spatial increments $\Delta x := (\Theta^\top \otimes I_n)x \in \mathbb{R}^{En}$ is globally (semi-globally, practically) asymptotically stable.

**Problem formulation:** Let $\mathcal{G}(\mathcal{V}, \mathcal{E})$ be a given undirected and connected graph fixing the network structure. Our target consists in designing a coupling law $u(t)$ that only uses the information of neighbors indicated by the graph $\mathcal{G}(\mathcal{V}, \mathcal{E})$, such that the network (8) achieves robust synchronization in the presence of matched disturbances coming from external signals and/or uncertain parameters in the individual models.

It is well known that under the presence of persistence disturbances, such as $\xi$ in (8), only practical synchronization is possible. Moreover, the coupling gain is the parameter controlling the size of the ultimate bound, so that $\text{dist}(x; S) \to 0$ as the coupling gain increases up to infinity, see e.g. [8]. In order to achieve robust synchronization with bounded coupling gains, we focus on coupling laws characterized by maximal monotone maps. Specifically, we consider the coupling law,

$$ u(t) \in -(\Theta_a W \Theta_a^\top \otimes I_p)w(t) - \gamma(\Theta_b \otimes I_p)\mathbf{M}((\Theta_b^\top \otimes I_p)w(t)), \quad (9) $$

where $\gamma > 0$; $W \in \mathbb{R}^{E_a \times E_a}$ is a diagonal matrix with $\min([W]_{i,i}) = \beta > 0$; $\mathbf{M}: \text{dom}\mathbf{M} \subseteq \mathbb{R}^{E_b p} \Rightarrow \mathbb{R}^{E_b p}$ is such that $\mathbf{M} = \mathbf{M}_1 \times \dots \times \mathbf{M}_{E_b}$, and $\mathbf{M}_j: \text{dom}\mathbf{M}_j \subseteq \mathbb{R}^p \Rightarrow \mathbb{R}^p$ is maximal monotone, for $j \in \{1, \dots, E_b\}$. The matrices $\Theta_a \in \mathbb{R}^{N \times E_a}$ and $\Theta_b \in \mathbb{R}^{N \times E_b}$ are, non-necessarily equal, oriented incidence matrices, both with rank $N-1$ associated, respectively, to the spanning subgraphs $\mathcal{G}_a(\mathcal{V}, \mathcal{E}_a) \subseteq \mathcal{G}$ and $\mathcal{G}_b(\mathcal{V}, \mathcal{E}_b) \subseteq \mathcal{G}$. That is, the coupling law (9) induces a multi-edge structure between vertices, (sometimes known as multi-layer or multiplex network, see e.g., [14]), where the sub-network $\mathcal{G}_a$ is driven by linear diffusion, whereas $\mathcal{G}_b$ is driven by the set-valued coupling law. Even though we are

considering only two layers, the developments that follow easily extend to the case of several layers. An additional constraint concerning $\mathbf{M}$ is stated below in Section IV. The well-posedness of the interconnected system (8)-(9) is a consequence of the following theorem.

**Theorem 1:** Consider the following differential inclusion,

$$ \dot{x}(t) \in h(t, x(t)) - \mathbf{F}(x(t)) \qquad (10) $$

where $\mathbf{F}: \text{dom}\mathbf{F} \subseteq \mathbb{R}^l \Rightarrow \mathbb{R}^l$ is a maximal monotone map, $h: I \times \text{dom}\mathbf{F} \to \mathbb{R}^l$ is measurable and such that for all $x \in \text{dom}\mathbf{F}$,

$$ \|h(t, x)\| \le \gamma_1(t)\|x\| + \gamma_2(t), $$

for some functions $\gamma_i \in L_{\text{loc}}^1(I)$. Then, for any initial condition $x(t_0) = x_0 \in \text{dom}\mathbf{F}$, there exists at least one absolutely continuous function $\chi(t, x_0)$ satisfying (10) for almost all times $t \in I$. Moreover, if $h(t, \cdot)$ is Lipschitz continuous then the solution is unique.

For a proof of Theorem 1 the reader is addressed to [15, Theorem 2.1]. The following corollary is an adaptation of Theorem 2 in [16] and state sufficient conditions guaranteeing the existence and uniqueness of solutions for (8).

**Corollary 1:** Assume that there exists a matrix $P = P^\top > 0$ such that $PB_2 = C_2^\top$ and $rge(\Theta_b^\top \otimes C_2 P^{-1/2}) \cap \text{rint}(\text{dom}\mathbf{M}) \neq \emptyset$. Then, for any initial condition $x(t_0) = x_0$ satisfying $(\Theta_b^\top \otimes C_2)x_0 \in \text{dom}\mathbf{M}$ there exists a unique absolutely continuous function $\chi(t, x_0)$ satisfying (8)-(9) for almost all times $t \ge t_0$.

**Proof:** Consider the change of coordinates $z = (I_N \otimes P^{1/2})x$, then (8)-(9) becomes:

$$ \dot{z} \in h(t, z) - \mathbf{F}(z), $$

where

$$ \begin{aligned} h(t, z) = & (I_N \otimes P^{1/2} A P^{-1/2})z \\ & - (I_N \otimes P^{1/2} B_1)\Phi\left((I_N \otimes C_1 P^{-1/2})z\right) \\ & + (I_N \otimes P^{1/2} B_2)\xi\left(t, (I_N \otimes P^{-1/2})z\right) \\ & - \left(\Theta_a W \Theta_a^\top \otimes P^{1/2} B_2 C_2 P^{-1/2}\right)z; \\ \mathbf{F} = & (\Theta_b \otimes P^{1/2} B_2) \circ \mathbf{M} \circ (\Theta_b^\top \otimes C_2 P^{-1/2}). \end{aligned} $$

Since $h(t, \cdot)$ is Lipschitz continuous, then for any $x \in \text{dom}\mathbf{F}$,

$$ \|h(t, x)\| \le L_h \|x\| + \|h(t, 0)\|, $$

and it satisfies the conditions stated in Theorem 1. Now, it follows from the hypothesis in $P$ and [13, Theorem 12.43] that $\mathbf{F}$ is also maximal monotone, and the conclusion follows as a direct consequence of Theorem 1. ■

## IV. PERFECT SYNCHRONIZATION VIA SET-VALUED MAXIMAL MONOTONE COUPLINGS

In this section it is shown that the coupling strategy (9) achieves asymptotic synchronization with a finite coupling $\gamma$ in the presence of matched disturbances. To that end, we consider the following assumption,
**Assumption 3:** There exists $\rho_M > 0$ such that M in (9) satisfies

$$ \mathcal{B}_{\rho_M} \subset \text{int} \mathbf{M}(0). \quad (11) $$

Note that (11) implies that $\mathbf{M}(0)$ is neither empty nor a singleton. In particular, linear coupling maps are excluded in the layer $\mathcal{G}_b$. Assumption 3 is also known as a *dry friction* condition in the literature of contact dynamics [17]. The following lemma will be useful in the developments that follow.

**Lemma 1:** Let $\mathbf{M}$ satisfies Assumption 3. Then for any $(\eta, \vartheta) \in \text{gph}(\mathbf{M})$, $\eta^\top \vartheta \ge \rho_M \|\eta\|$.

**Proof:** Let $(\eta, \vartheta) \in \text{gph}(\mathbf{M})$. It follows from the assumption on $\mathbf{M}(0)$ that $(0, \vartheta) \in \text{gph}(\mathbf{M})$ for all $\vartheta \in \mathcal{B}_{\rho_M}$. Hence, monotonicity of $\mathbf{M}$ implies that,

$$ \langle \eta, \vartheta - \hat{\vartheta} \rangle \ge 0 \text{ for all } \hat{\vartheta} \in \mathcal{B}_{\rho_M}. $$

Therefore, $\langle \eta, \vartheta \rangle \ge \sup_{\hat{\vartheta} \in \mathcal{B}_{\rho_M}} \langle \eta, \hat{\vartheta} \rangle = \rho_M \|\eta\|$, and the implication follows. ■

The following theorem shows how maximal monotone set-valued maps satisfying Assumption 3 achieve perfect asymptotic synchronization with finite coupling gain $\gamma$.

**Theorem 2:** Let Assumptions 1-2 and 3 hold. Let the initial condition $x(t_0) = x_0$ be such that $(\Theta_b^\top \otimes C_2)x_0 \in \text{dom}\mathbf{M}$, and let $0 < \bar{\xi}$ be such that $\|\xi_k(t, x(t))\| \le \bar{\xi}_k$ for all $k \in \{1, \dots, N\}$ and for almost all $t \ge t_0$. If

1) there exist $\mu > 0$ and a matrix $P = P^\top > 0$ satisfying the conditions of Corollary 1, together with

$$ \begin{bmatrix} \tilde{Q}_1 + \mu P & PB_1 - C_1^\top(K_1 + K_2)^\top \\ * & -2I_m \end{bmatrix} < 0, \quad (12) $$

where $\tilde{Q}_1 = A^\top P + PA - 2\beta\lambda_{2,a}C_2^\top C_2 - C_1^\top(K_1^\top K_2 + K_2^\top K_1)C_1$ and $\lambda_{2,a}$ is the second smallest eigenvalue of the Laplacian $L_a = \Theta_a \Theta_a^\top$; and

2) the gain $\gamma > 0$ is such that

$$ \gamma \rho_M \sqrt{\lambda_{2,b}} > \|\bar{\xi}\|, \quad (13) $$

where $\bar{\xi} = [\bar{\xi}_1, \dots, \bar{\xi}_N]^\top \in \mathbb{R}^N$.

Then, the perturbed network (8) with coupling law (9) achieves global asymptotic synchronization.

**Proof:** Let us consider the Lyapunov function candidate,

$$ V(x) = \frac{1}{2N} x^\top (\Theta_c \Theta_c^\top \otimes P)x, \quad (14) $$

where $P = P^\top > 0$ satisfies all assumptions of the theorem and $\Theta_c \in \mathbb{R}^{N \times E_c}$ is an oriented incidence matrix associated to the complete graph $K_N$, so that $E_c = \frac{N(N-1)}{2}$. Note that $V(x) = 0$ if and only if $x$ lies inside the synchronization manifold $S = \text{rge}(\mathbf{1}_N \otimes I_n)$. Computing the time derivative of $V$ along the trajectories of (8) with coupling (9) yields

$$ \begin{aligned} \frac{d}{dt} V(x) \le \frac{1}{2N} x^\top (\Theta_c \Theta_c^\top \otimes A^\top P + PA)x \\ - \frac{1}{N} x^\top (\Theta_c \Theta_c^\top \otimes PB_1) \Phi(y) + \frac{1}{N} \|(\Theta_c^\top \otimes I_p)w\| \|(\Theta_c^\top \otimes I_p)\xi\| \\ - w^\top (\Theta_a W \Theta_a^\top \otimes I_p)w - \gamma w^\top (\Theta_b \otimes I_p)\vartheta_b \end{aligned} \quad (15) $$

where, $\vartheta_b \in \mathbf{M}((\Theta_b^\top \otimes I_p)w)$ and we have used the fact that $\Theta_c \Theta_c^\top \Theta = (N I_N - \mathbf{1}_N \mathbf{1}_N^\top)\Theta = N\Theta$ for any $\Theta$ satisfying $\Theta^\top \mathbf{1}_N = 0$, to obtain the last two terms. Now, recalling that $\min_i([W]_{i,i}) = \beta$, it follows from the Courant-Fisher minimax theorem [18, Theorem 4.2.6] that,

$$ w^\top (\Theta_a W \Theta_a^\top \otimes I_p)w \ge \frac{\beta \lambda_{2,a}}{N} x^\top (\Theta_c \Theta_c^\top \otimes C_2^\top C_2) x, \quad (16) $$

where $\lambda_{2,a}$ denotes the (unweighted) algebraic connectivity of the subnetwork $\mathcal{G}_a$. Making use of (16) together with Lemma 1, back into (15) leads us to,

$$ \begin{aligned} \frac{d}{dt} V(x) \le \frac{1}{2N} x^\top (\Theta_c \Theta_c^\top \otimes Q_1) x \\ - \frac{1}{N} x^\top (\Theta_c \Theta_c^\top \otimes PB_1) \Phi(y) - \gamma \rho_M \|(\Theta_b^\top \otimes I_p)w\| \\ + \frac{1}{N} \|(\Theta_c^\top \otimes I_p)w\| \|(\Theta_c^\top \otimes I_p)\xi\|, \end{aligned} \quad (17) $$

where $Q_1 := (A^\top P + PA - 2\beta\lambda_{2,a}C_2^\top C_2)$. It follows from simple calculations that

$$ \|(\Theta_c^\top \otimes I_p)\xi\| \le \sqrt{N} \|\bar{\xi}\|. \quad (18) $$

Hence, the substitution of (18) into (17), and using once again the Courant-Fisher minimax theorem, leads us to

$$ \begin{aligned} \frac{d}{dt} V(x) \le \frac{1}{2N} x^\top (\Theta_c \Theta_c^\top \otimes Q_1) x \\ - \frac{1}{N} x^\top (\Theta_c \Theta_c^\top \otimes PB_1) \Phi(y) \\ - \frac{1}{\sqrt{N}} \left( \gamma \rho_M \sqrt{\lambda_{2,b}} - \|\bar{\xi}\| \right) \|(\Theta_c^\top \otimes I_p)w\|. \end{aligned} \quad (19) $$

Now we turn to the incremental sector condition (7). Noticing that each row of $\Theta_c^\top$ has exactly two non-zero entries given by +1 and -1, it follows from (7) that

$$ 0 \le -\frac{1}{2N} \begin{bmatrix} \Delta_c^n x \\ \Delta_c^m \Phi \end{bmatrix}^\top \begin{bmatrix} I_{E_c} \otimes D_1 & I_{E_c} \otimes D_2 \\ I_{E_c} \otimes D_2 & 2I_{E_c m} \end{bmatrix} \begin{bmatrix} \Delta_c^n x \\ \Delta_c^m \Phi \end{bmatrix}, \quad (20) $$

where $\Delta_c^n x = (\Theta_c^\top \otimes I_n)x$, $\Delta_c^m \Phi = -(\Theta_c^\top \otimes I_m)\Phi(y)$, $D_1 = C_1^\top(K_1^\top K_2 + K_2^\top K_1)C_1$, and $D_2 = (K_1 + K_2)C_1$. Hence, the addition of (19) and (20) leads us to,

$$ \begin{aligned} \frac{d}{dt} V(x) \le \frac{1}{2N} \begin{bmatrix} \Delta_c^n x \\ \Delta_c^m \Phi \end{bmatrix}^\top \begin{bmatrix} I_{E_c} \otimes \tilde{Q}_1 & I_{E_c} \otimes \tilde{Q}_2 \\ I_{E_c} \otimes \tilde{Q}_2^\top & -2I_{E_c m} \end{bmatrix} \begin{bmatrix} \Delta_c^n x \\ \Delta_c^m \Phi \end{bmatrix} \\ - \frac{1}{\sqrt{N}} \left( \gamma \rho_M \sqrt{\lambda_{2,b}} - \|\bar{\xi}\| \right) \|(\Theta_c^\top \otimes I_p)w\|, \end{aligned} \quad (21) $$

where $\tilde{Q}_1$ is as in (12) and $\tilde{Q}_2 = PB_1 - D_2^\top$. It thus follows from (12) and (13) that $\frac{d}{dt}V(x) \le -\mu V(x)$. Finally, recalling that $V(x) = 0$ if and only if $x \in \text{rge}(\mathbf{1}_N \otimes I_n)$, global asymptotic synchronization (with a rate of convergence of at least $\mu$) follows. This concludes the proof. ■

Note that both, the matrix $W$ and the gain $\gamma$ are bounded, as opposed to linear designs, where perfect synchronization can only be attained with infinite gains [8], [9]. By contrast, set-valued maps can achieve perfect asymptotic synchronization with finite coupling strength. However, Assumption 3
implies that the coupling map has infinite incremental gain. In the following subsection we study how the implementable controller has an intrinsic limitation in precision that comes from the regularization and not from the magnitude of the coupling gains.

**Corollary 2:** Let $p = n$ and $C_2 = I_n$. If all the assumptions of Theorem 2 hold, then the network achieves perfect synchronization in finite time with the upper bound

$$t^* \le \frac{\sqrt{2\lambda_{\max}(P)V(x_0)}}{\gamma\rho_M\sqrt{\lambda_{2,b}} - \|\bar{\xi}\|}. \quad (22)$$

*Proof:* It follows directly from the assumptions and the integration of (21) in the interval $[0, t]$. ■ The conditions in Theorem 2 can be relaxed to get a semi-global version as stated in the following corollary. In what follows we impose the following detectability assumption as appears in [19].

**Assumption 4:** The incremental dynamics is asymptotically zero-state detectable. That is, for any two systems $v_a, v_b$ such that the output mismatch $w_a(t) - w_b(t) \to 0$, whereas the inputs $\tilde{u}_r(t) = u_r(t) + \xi_r(t, x_r(t))$, $r \in \{a, b\}$ satisfy $\tilde{u}_a(t) \to 0$, $\tilde{u}_b(t) \to 0$ implies $x_a(t) - x_b(t) \to 0$.

**Corollary 3:** Let Assumptions 1-4 and 3 hold. Moreover, let the initial condition $x(t_0) = x_0 \in \mathbb{R}^{Nn}$ be such that $(\Theta_b^\top \otimes C_2)x_0 \in \text{dom}M$ and for some $c > 0$ $x_0 \in \text{int}\Omega_c$ where,

$$\Omega_c := \{x \in \mathbb{R}^{Nn} | V(x) \le c\}. \quad (23)$$

If

1) there exist matrices $R = R^\top$, such that $\lambda_{\max}(R) > 0$, and $P = P^\top > 0$ satisfying the conditions of Corollary 1 together with

$$\begin{bmatrix} \tilde{Q}_1 - C_2^\top R C_2 & PB_1 - C_1^\top (K_1 + K_2)^\top \\ * & -2I_m \end{bmatrix} \preceq 0, \quad (24)$$

where $\tilde{Q}_1$ is the same as in Theorem 2; and

2) the gain $\gamma > 0$ is such that

$$\gamma \rho_M \sqrt{\lambda_{2,b}} > \|\bar{\xi}\| + \frac{\lambda_{\max}(R) \|C_2\| \sqrt{c}}{\sqrt{2\lambda_{\min}(P)}} \quad (25)$$

Then the perturbed network (8) with monotone coupling (9) achieves semi-global asymptotic synchronization.

*Proof:* Taking the same Lyapunov function candidate as in the proof of Theorem 2 we arrive at (21). Hence, it follows from (24) that

$$\frac{d}{dt} V(x) \le \frac{\lambda_{\max}(R)}{2N} x^\top (\Theta_c \Theta_c^\top \otimes C_2^\top C_2) x - \frac{1}{\sqrt{N}} (\gamma \rho_M \sqrt{\lambda_{2,b}} - \|\bar{\xi}\|) \|(\Theta_c^\top \otimes I_p)w\|, \quad (26)$$

By assumption $x_0 \in \text{int}\Omega_c$. Hence, there exists $t_1 > t_0$ such that for all $t \in [t_0, t_1]$, $x(t) \in \Omega_c$. Thus, for $t \in [t_0, t_1]$ one has that

$$\frac{1}{2N} x^\top (\Theta_c \Theta_c^\top \otimes C_2 C_2^\top) x \le \frac{c \|C_2\|^2}{\lambda_{\min}(P)}. \quad (27)$$

Consequently, for $t \in [t_0, t_1]$, the derivative of V satisfies

$$\frac{d}{dt} V(x) \le -\frac{1}{\sqrt{N}} \left( \gamma \rho_M \sqrt{\lambda_{2,b}} - \|\bar{\xi}\| - \frac{\lambda_{\max}(R) \|C_2\| \sqrt{c}}{\sqrt{2\lambda_{\min}(P)}} \right) \|(\Theta_c^\top \otimes I_p)w\| \quad (28)$$

It follows from (25) that for $t \in [t_0, t_1]$, V is a non-increasing function and therefore $x(t_1) \in \text{int}\Omega_c$. Following an induction argument, we conclude that for all $t \ge t_0$, $x(t) \in \Omega_c$, that is, $\Omega_c$ is a positively invariant set of the closed-loop (8)-(9). Finally, it follows from (28) that $(\Theta_c^\top \otimes I_p)w \to 0$ and the zero-state detectability assumption guarantees that $(\Theta_c^\top \otimes I_p)x \to 0$. That is, $\text{dist}(x(t); S) \to 0$ whenever $x(t) \in B_{s(c)} \subset \Omega_c$, where $s(c) = \frac{c}{N\lambda_{\min}(P)}$, and semi-global asymptotic synchronization follows. ■

Corollary 3 allows for designs when the subgraph $G_a$ does not provide sufficient energy dissipation to achieve global synchronization. In such a case, the set-valued subnetwork induced by $G_b$ is used to compensate the lack of dissipation by achieving a trade-off between the size of the region of attraction, estimated by $\Omega_c$, and the size of the gain $\gamma$.

It is also worth to remark that if $\lambda_{\max}(R) \le 0$ in (24) then (26) together with the zero-state detectability assumption imply that the network (8) with coupling law (9) achieves global asymptotic synchronization.

## V. IMPLEMENTATION OF SET-VALUED SIGN-LIKE COUPLINGS

For achieving perfect regulation, as stated in Theorem 2 or Corollary 3, it is necessary to know the *exact* values (selection) of the control $u$ in the set $M(0)$ that will counteract the effects caused by the matched disturbances. Even though such scenario is unreal, the usefulness of Theorem 2 amounts to guide us in the design of regularized controllers achieving *approximate* selections. Our choice over approximate selections is motivated by the mechanisms envisaged for implementing the coupling signals using analog circuitry. For instance, let us consider the analog circuit in Figure 1 where each diode satisfies an *ideal* complementarity relation as,

$$0 \le I_{D_k} \perp V^* - V_{D_k} \ge 0, \quad (29)$$

where the notation $0 \le a \perp b \ge 0$ stands for the following three conditions: i) $a \ge 0$, ii) $b \ge 0$ and iii) $ab = 0$, $I_{D_k}$, $V_{D_k}$ denote the current through and the voltage across the k-th diode, respectively, and $V^*$ denotes the activation voltage of the diode. In real-life circuits, each diode will have

<figure><img src="image_2.png" /><figcaption>Fig. 1. Analog circuit implementing an approximation of the set-valued coupling law $u_{(j,i)} = -u_{(i,j)} \in \gamma \text{Sgn}(w_i - w_j)$.</figcaption></figure>
parasitic dynamic and resistance effects which are taken into account via the resistance R and the inductance L. It follows from Kirchhoff's laws that the coupling circuit in Figure 1 satisfies,

$$-u_{(i,j)} = u_{(j,i)} = z_1 - \gamma = \gamma - z_2 \quad (30a)$$

$$w_i - w_j = \vartheta + R(z_1 - z_2) + L(\dot{z}_1 - \dot{z}_2) \quad (30b)$$

$$\vartheta = V_{D_1} - V_{D_2} \quad (30c)$$

$$0 \le z_1 \perp V^* - V_{D_1} \ge 0 \quad (30d)$$

$$0 \le z_2 \perp V^* - V_{D_2} \ge 0 \quad (30e)$$

It follows from (30d)-(30e) and (30a) that, (regardless of the value of $z_1$ and $z_2$), $-u_{(i,j)} = u_{(j,i)} \in [-\gamma, \gamma]$. Moreover, we have the following three cases,

1) $u_{(j,i)} = -\gamma$. It follows from (30a) that $z_2 > 0$ and (30e) implies that $V^* - V_{D_2} = 0$. Hence, (30c) and (30d) imply that $\vartheta \in (-\infty, 0]$.

2) $u_{(j,i)} = \gamma$. It follows from (30a) that $z_1 > 0$ and (30d) implies that $V^* - V_{D_1} = 0$. Hence, (30c) and (30e) imply that $\vartheta \in [0, +\infty)$.

3) $-\gamma < u_{(j,i)} < \gamma$. It follows from (30a) that $z_1 > 0$ and $z_2 > 0$. Hence the complementarity conditions (30d)-(30e) imply that $V^* - V_{D_1} = V^* - V_{D_2} = 0$ and it follows from (30c) that $\vartheta = 0$.

It is not difficult to see that the three cases above characterize the relation $\vartheta \in \mathbf{N}_{[-\gamma, \gamma]}(u_{(j,i)})$. Equivalently, $u_{(j,i)} \in \gamma \mathbf{Sgn}(\vartheta)$. Hence, it follows from (30a)-(30b), that

$$2L\dot{u}_{(j,i)} \in -2Ru_{(j,i)} - \mathbf{N}_{[-\gamma, \gamma]}(u_{(j,i)}) + (w_i - w_j) \quad (31)$$

Since $\mathbf{N}_{[-\gamma, \gamma]}$ is a maximal monotone operator, then the dynamics (31) is well-possed. Notice that, in the ideal case where $R = 0$ and $L = 0$, (31) yields

$$u_{(j,i)}^{\text{ideal}} \in \gamma \mathbf{Sgn}(w_i - w_j),$$

and we retrieve the ideal set-valued sign coupling. The dynamics (31) has a unique equilibrium $u_{(j,i)}^*$, driven by the difference $w_i - w_j$ and satisfying,

$$\frac{w_i - w_j}{2R} \in \left( I + \frac{1}{2R} \mathbf{N}_{[-\gamma, \gamma]} \right) (u_{(j,i)}^*).$$

Equivalently,

$$u_{(j,i)}^* = \mathcal{J}_{\frac{1}{2R}\mathbf{N}_{[-\gamma, \gamma]}} \left( \frac{w_i - w_j}{2R} \right).$$

Thus, it follows from (5) that

$$u_{(j,i)}^* = \mathcal{Y}_{\gamma \mathbf{Sgn}}^{2R}(w_i - w_j). \quad (32)$$

Moreover, the dynamics (31) is incrementally stable for the any two trajectories with the same external input $w_i - w_j$. Indeed, setting $\Delta u_{(j,i)} = u_{(j,i)}^a - u_{(j,i)}^b$ for any two trajectories $u_{(j,i)}^k$ of (31), $k \in \{a, b\}$, it follows that $V(\Delta u_{(j,i)}) = L(\Delta u_{(j,i)})^2$ satisfies

$$\begin{aligned} \frac{d}{dt} V(\Delta u_{(j,i)}) &= \Delta u_{(j,i)}(-2R\Delta u_{(j,i)} - (\vartheta_a - \vartheta_b)) \\ &\le -2R(\Delta u_{(j,i)})^2 = -2\frac{R}{L}V(\Delta u_{(j,i)}). \quad (33) \end{aligned}$$

Hence, (31) is (exponentially) contracting with rate $\lambda = \frac{R}{L}$. Integration on both sides of (33) leads us to,

$$\|u_{(j,i)}^a(t) - u_{(j,i)}^b(t)\| \le e^{-\frac{R}{L}t} \|u_{(j,i)}^a(0) - u_{(j,i)}^b(0)\|.$$

Recalling that for all time $t > 0$ $u_{(j,i)}(t) \in [-\gamma, \gamma]$, it follows that

$$u_{(j,i)}(t) \in u_{(j,i)}^*(t) + 2\gamma e^{-\frac{R}{L}t} \mathcal{B}_1. \quad (34)$$

Therefore, the coupling circuit in Figure 1 converges exponentially fast towards the Yosida approximation in (32). In what follows we consider the case where the couplings on $\mathcal{G}_b$ are scalar systems satisfying (31). It follows from (34) that it is equivalent to consider coupling laws given by

$$u = -(\Theta_a W \Theta_a^\top \otimes I_p)w - \gamma(\Theta_b \otimes I_p)\mathcal{Y}_M^\varepsilon ((\Theta_b^\top \otimes I_p)w), \quad (35)$$

where $\varepsilon = 2R > 0$, $\mathbf{M}: \mathbb{R}^{E_b} \to \mathbb{R}^{E_b}$ is the componentwise signum multifunction such that $\eta \mapsto \mathbf{Sgn}(\eta_1) \times \dots \times \mathbf{Sgn}(\eta_{E_b})$ and the vector of matched perturbations is as

$$\tilde{\xi}(t, x) = \xi(t, x) + 2\gamma\delta_b e^{-\frac{R}{L}t}\zeta(t),$$

where $\zeta(t) \in \mathcal{B}_1$ and $\delta_b$ is the maximum degree of nodes in $\mathcal{G}_b$.

**Lemma 2:** Let $\mathbf{M}$ be a maximal monotone map satisfying Assumption 3. Then, for any $\varepsilon > 0$, $\|\mathcal{Y}_M^\varepsilon(\vartheta)\| < \rho_M$, if and only if, $\vartheta \in \text{int}(\varepsilon\mathcal{B}_{\rho_M})$.

**Proof:** Let $\vartheta \in \text{int}(\varepsilon\mathcal{B}_{\rho_M})$. It follows from (2) that $\mathcal{J}_{\varepsilon\mathbf{M}}(\vartheta) = 0$. Hence, (3) implies that

$$\|\mathcal{Y}_M^\varepsilon(\vartheta)\| = \frac{\|\vartheta - \mathcal{J}_{\varepsilon\mathbf{M}}(\vartheta)\|}{\varepsilon} = \frac{\|\vartheta\|}{\varepsilon} < \rho_M.$$

For the converse, let $\vartheta \notin \text{int}(\varepsilon\mathcal{B}_{\rho_M})$. In particular, $\vartheta \ne 0$, $\varepsilon\rho_M \frac{\vartheta}{\|\vartheta\|} \in \varepsilon\mathcal{B}_{\rho_M} \subset \varepsilon\mathbf{M}(0)$, and $\mathcal{J}_{\varepsilon\mathbf{M}}(\varepsilon\rho_M \frac{\vartheta}{\|\vartheta\|}) = 0$. Thus, the nonexpansiveness property of the resolvant implies that

$$\begin{aligned} \|\mathcal{J}_{\varepsilon\mathbf{M}}(\vartheta)\| &= \left\| \mathcal{J}_{\varepsilon\mathbf{M}}(\vartheta) - \mathcal{J}_{\varepsilon\mathbf{M}}\left(\varepsilon\rho_M \frac{\vartheta}{\|\vartheta\|}\right) \right\| \\ &\le \left\| \vartheta - \varepsilon\rho_M \frac{\vartheta}{\|\vartheta\|} \right\| = \|\vartheta\| - \varepsilon\rho_M. \quad (36) \end{aligned}$$

Therefore,

$$\begin{aligned} \|\mathcal{Y}_M^\varepsilon(\vartheta)\| &= \frac{1}{\varepsilon}\|\vartheta - \mathcal{J}_{\varepsilon\mathbf{M}}(\vartheta)\| \ge \frac{1}{\varepsilon}(\|\vartheta\| - \|\mathcal{J}_{\varepsilon\mathbf{M}}(\vartheta)\|) \\ &\ge \frac{1}{\varepsilon}(\|\vartheta\| - \|\vartheta\| + \varepsilon\rho_M) = \rho_M, \end{aligned}$$

and the conclusion follows.

**Corollary 4:** Let all the assumptions of Theorem 2 hold with $p = 1$ and consider the coupling (31). Then the asymptotic behavior of the network (8) with coupling (31) is practically synchronized. Moreover, an estimation of the ultimate bound is given by the set

$$\Omega_{\varepsilon}(\|\bar{\xi}\|) = \left\{ x \in \mathbb{R}^{Nn} \mid V(x) \le \varepsilon \sqrt{\frac{N}{\lambda_{2,b}} \frac{\|\bar{\xi}\|}{\mu}} \right\}, \quad (37)$$

where $\mu$ and $\bar{\xi}$ are the same as in Theorem 2.
*Proof:* Taking the Lyapunov function candidate (14), simple computations, similar to those made in the proof of Theorem 2, lead us to

$$ \frac{d}{dt} V(x) \leq -\frac{\mu}{2N} x^{\top} (\Theta_c \Theta_c^{\top} \otimes P) x \\ - \gamma w^{\top} (\Theta_b \otimes I_p) \mathcal{Y}_{\mathbf{M}}^{\varepsilon} ((\Theta_b^{\top} \otimes I_p)w) \\ + \left( \frac{\|\bar{\xi}\| + \kappa(t)}{\sqrt{N}} \right) \|(\Theta_c^{\top} \otimes I_p)w\|, \quad (38) $$

where $\kappa(t) = 2\gamma\sqrt{N}\delta_b e^{-\frac{R}{L}t}$. Now, the substitution of (3) into (38) yields,

$$ \frac{d}{dt} V(x) \leq -\mu V(x) + \left( \frac{\|\bar{\xi}\| + \kappa(t)}{\sqrt{N}} \right) \|(\Theta_c^\top \otimes I_p)w\| \\ - \varepsilon\gamma \|\mathcal{Y}_M^\varepsilon((\Theta_b^\top \otimes I_p)w)\|^2 \\ - \gamma \rho_M \|\mathcal{J}_{\varepsilon M}((\Theta_b^\top \otimes I_p)w)\|, \quad (39) $$

where we have used Lemma 1 and (4) to get the last term. We divide the analysis into two cases.

a) : Assume $\|(\Theta_b^\top \otimes I_p)w\| \notin \varepsilon\mathcal{B}_{\rho_M}$. Hence, using once again (3) together with the triangle inequality in the second term of (39), leads us to

$$ \frac{d}{dt} V(x) \leq -\mu V(x) \\ - \varepsilon\gamma \left( \|\mathcal{Y}_{\mathbf{M}}^{\varepsilon}((\Theta_b^\top \otimes I_p)w)\| - \rho_{\mathbf{M}} \right) \|\mathcal{Y}_{\mathbf{M}}^{\varepsilon}((\Theta_b^\top \otimes I_p)w)\| \\ - \frac{\gamma \rho_{\mathbf{M}} \sqrt{\lambda_{2,b}} - \|\bar{\xi}\| - \kappa(t)}{\sqrt{N}} \|(\Theta_c^\top \otimes I_p)w\|. \quad (40) $$

It follows from Lemma 2 and (13) that the right-hand side of (40) is strictly negative for $t > t^*$ where

$$ t^* > -\frac{L}{R} \log \left( \frac{\gamma \sqrt{\lambda_{2,b}} - \|\bar{\xi}\|}{2\gamma \delta_b \sqrt{N}} \right), \quad (41) $$

with $\frac{L}{R} \ll 1$ so that the time $t^*$ remains “small” in comparison to the temporal scale of the agents.

b) : Now assume that $(\Theta_b^\top \otimes I_p)w \in \varepsilon\mathcal{B}_{\rho_M}$ but $x \notin \Omega_\varepsilon(\|\bar{\xi}\| + \kappa(t))$. In such case it follows from (2) and Assumption 3 that $\mathcal{J}_{\varepsilon M}((\Theta_b^\top \otimes I_p)w) = 0$, and after applying the Courant-Fisher minimax theorem (39) becomes,

$$ \frac{d}{dt} V(x) \leq -\mu \left( V(x) - \varepsilon \sqrt{\frac{N}{\lambda_{2,b}}} \frac{\|\bar{\xi}\| + \kappa(t)}{\mu} \right) \\ - \frac{\gamma}{\varepsilon} \|(\Theta_b^\top \otimes I_p)w\|^2. \quad (42) $$

Hence if $x \notin \Omega_\varepsilon(\|\bar{\xi}\| + \kappa(t))$ then (42) is strictly negative. It follows that for all $t > t^*$ both, (40) and (42), are strictly negative whenever $x \notin \Omega_\varepsilon(\|\bar{\xi}\| + \kappa(t))$. Therefore, $\mathrm{dist}(x; \Omega_\varepsilon(\|\bar{\xi}\|)) \to 0$ as $t \to \infty$ and global practical synchronization follows. ■

From the proof of Corollary 4 it becomes clear that in order to have high precision in the presence of matched disturbances, it is better to look for a Yosida approximation with small index $\varepsilon$, rather than increasing the gain $\gamma$. Moreover, note that the dynamics in (31) is independent of the agent

parameters, as opposed to the network (8)-(9), allowing thus for distributed selection schemes via implicit methods [20].

*Example 1:* As an illustration, let us consider an heterogeneous group of $N = 32$ FitzHugh-Nagumo oscillators with agent dynamics (6) and parameters given by,

$$ A = \begin{bmatrix} 0 & -2 \\ \frac{1}{\delta} & \frac{0.1}{\delta} \end{bmatrix}, B_1 = B_2 = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, C_1 = \begin{bmatrix} 1 & 0 \end{bmatrix}, \quad (43) $$

where $\delta = 0.05$, the matrix $C_2$ is defined below, and each agent $v_k$ is closed with a distinct nonlinear feedback $\varphi_k: \mathbb{R} \to \mathbb{R}$, $\varphi_k(\eta) = \eta^3 - \alpha_k \eta$, where each constant $\alpha_k$ is selected in a random way from the interval $[2.5, 10]$. In addition, each agent is affected by external disturbances of the form $\xi_k = 2 \cos(\beta_k t) \sin(x_{k,1}) + 2 \sin(kx_{k,1}) \cos(\sqrt{k}x_{k,2})$, where each constant $\beta_k$ is selected randomly from the interval $[15, 20]$. Note that extra disturbances will appear due to the heterogeneity between agents. Indeed, it follows from simple computations that for any $\{i, j\} \in \mathcal{E}_c$

$$ \varphi_i(y_i) - \varphi_j(y_j) = \bar{\varphi}_{i,j}(y_i) - \bar{\varphi}_{i,j}(y_j) + \frac{(\alpha_j - \alpha_i)(y_i + y_j)}{2}, \quad (44) $$

where $\bar{\varphi}_{i,j} = \frac{\varphi_i + \varphi_j}{2}$. Moreover, for any $\{i, j\} \in \mathcal{E}_c$, the map $\bar{\varphi}_{i,j}$ satisfies (7) in $\mathcal{B}_s$ for some $s > 0$, whenever, $K_1 \le \min_{\{i,j\} \in \mathcal{E}_c\} \{-\bar{\alpha}_{i,j}\}$ and $\max_{\{i,j\} \in \mathcal{E}_c\} \{3s^2 - \bar{\alpha}_{i,j}\} \le K_2$, where $\bar{\alpha}_{i,j} = \frac{\alpha_i + \alpha_j}{2}$. Taking the same Lyapunov function candidate and following the same steps as in the proof of Theorem 2, *mutatis mutandis*, it follows that if $P$ satisfies (12) and $\gamma$ satisfies

$$ \gamma \rho_{\mathbf{M}} \sqrt{\lambda_{2,b}} \geq \|\bar{\xi}\| + r(s), \quad (45) $$

where $r(s) = \max_{\{i,j\} \in \mathcal{E}_c} |\alpha_j - \alpha_i|$ $s$, then semi-global practical synchronization follows. The same conclusion (with different ultimate bound and region of attraction) is obtained if we consider instead the assumptions of Corollary 3. In that case, $\gamma$ must satisfy,

$$ \gamma \rho_{\mathbf{M}} \sqrt{\lambda_{2,b}} \geq \|\bar{\xi}\| + \frac{\lambda_{\max}(R) \|C_2\| \sqrt{c}}{\sqrt{2\lambda_{\min}(P)}} + r(s) \quad (46) $$

where $c > 0$ is such that (23) holds. The information available to each agent is described by a so-called small-world network of type Newman-Watts-Strogatz, where each vertex is connected to its $k = 8$ nearest neighbors and with probability of extra connections of $p = 0.45$, see [21] for details. The coupling configuration is as follows: $\mathcal{G}_a$ is the same as $\mathcal{G}$ with unitary weights, that is $W = I_{E_a}$, and the subnetwork $\mathcal{G}_b$ constitutes the ring of vertices from agent $v_1$ up to agent $v_N$ with coupling circuit satisfying (31). Thus, the nonlinear coupling acts only on a subnetwork of $\mathcal{G}$. Also, it follows from the network configuration that, $\lambda_{2,a} = 4.36$ and $\sqrt{\lambda_{2,b}} = 0.196$.

We now look for the output matrix $C_2$ such that (12) holds. For the case $C_2 = C_1$, the relaxed LMI (24) has a solution with a minimal value of $R = 736.14$ and $\lambda_{\min}(P) = 0.1$, leading to a large value of the gain $\gamma$ according to (46). Thus, in order to reduce the aforementioned conservativeness, we
consider $C_2$ in (12) as another variable. After standard operations (a loop transformation changing the incremental sector of $\varphi$ to the incremental sector $[0, K_2 - K_1]$, together with a congruence transformation in (12)), we arrive at the following LMI,

$$ \begin{bmatrix} Q_{1,1} & B_1 - P^{-1}C_1^\top \tilde{K}^\top \\ B_1^\top - \tilde{K}C_1 P^{-1} & -2I_m \end{bmatrix} \preceq 0, $$

where $Q_{1,1} = P^{-1}(A+B_1K_1C_1)^\top + (A+B_1K_1C_1)P^{-1} - \beta\lambda_{2,a}B_2B_2^\top + \mu P^{-1}$ and $\tilde{K} = K_2 - K_1$. By setting $C_2 = B_2^\top P$ we recover the output matrix. In this case, with the parameters aforementioned, we obtain that $C_2 = [8.53, -0.21]$ with $\mu = 0.1$. Therefore, setting $\gamma > 5.2(\|\tilde{\xi}\| + r(s))$ guarantees the semiglobal practical synchronization of the perturbed network.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Time evolution of states $x_{k,i}(t)$, $i \in \{1, 2\}$, of Example 1 with nonsmooth coupling (31) with $R = 10^{-2}\Omega$, $L = 0.1\mu H$, $W = I_{E_a}$, and $\gamma = 25$. The black line represents the average behavior $\bar{x}(t)$.</figcaption></figure>

Figure 2 shows the time trajectories of the state variables with the dynamic coupling (31) with $R = 10^{-2}\Omega$, $L = 0.1\mu H$, and $\gamma = 25$. The sum of squares error signal, $e_{sos}(t) := \frac{1}{N}x^\top(\Theta_c\Theta_c^\top \otimes I_n)x$ is displayed in Figure 3 for four different values of resistance $R$, verifying the claimed practical synchronization of the network. All simulations were performed using a backward (implicit) Euler method for (31) and a forward (explicit) Euler method for the agent dynamics (6) with a sampling time of $h = 100\mu s$.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Sum of squares error for the network of Example 1 at four different regularizations: a) $R = 10^{-2}\Omega$ – dashed gray line; b) $R = 10^{-3}\Omega$ – dashed black line; c) $R = 10^{-4}\Omega$ – continuous black line; and d) $R = 10^{-5}\Omega$ – continuous gray line. In all cases $L = 0.1\mu H$, $W = I_{E_a}$, and $\gamma = 25$.</figcaption></figure>

## VI. CONCLUSIONS AND FURTHER RESEARCH

The problem of robust synchronization against matched perturbations was dealt in this paper. It was shown that max-

imal monotone maps with dry friction achieve perfect regulation in the presence of the aforementioned disturbances. Afterwards a dynamic coupling is proposed for implementing the coupling laws. It is shown that practical synchronization is attained and precision improves as the index of the Yosida approximation approaches zero. Further research considers more general dynamic nonsmooth couplings generalizing (31).

## REFERENCES

[1] P. DeLellis, M. di Bernardo, T. E. Gorochowski, and G. Russo, “Synchronization and control of complex networks via contraction, adaptation and evolution,” *IEEE Circuits and Systems Magazine*, pp. 64–82, 2010.

[2] R. Olfati-Saber, J. A. Fax, and R. M. Murray, “Consensus and cooperation in networked multi-agent systems,” *Proceeding of the IEEE*, vol. 95, no. 1, pp. 215–233, 2007.

[3] E. Panteley and A. Loría, “Synchronization and dynamic consensus of heterogeneous networked systems,” *IEEE Transactions on Automatic Control*, vol. 62, no. 8, pp. 3758–3773, 2017.

[4] L. Scardovi and R. Sepulchre, “Synchronization in networks of identical linear systems,” *Automatica*, vol. 45, pp. 2557–2562, 2009.

[5] H. L. Trentelman, K. Takaba, and N. Monshizadeh, “Robust synchronization of uncertain linear multi-agent systems,” *International Journal of Robust and Nonlinear Control*, vol. 58, no. 6, pp. 1511–1523, 2013.

[6] W. Wang and J.-J. E. Slotine, “On partial contraction analysis for coupled nonlinear oscillators,” *Biological Cybernetics*, vol. 92, p. 2004, 2004.

[7] P. DeLellis, M. di Bernardo, and D. Liuzza, “Convergence and synchronization in heterogeneous networks of smooth and piecewise smooth systems,” *Automatica*, vol. 56, pp. 1–11, 2015.

[8] J. M. Montenbruck, M. Bürger, and F. Allgöwer, “Practical synchronization with diffusive couplings,” *Automatica*, vol. 53, pp. 235–243, 2015.

[9] J. Zhao, D. J. Hill, and T. Liu, “Global bounded synchronization of general dynamical networks with nonidentical nodes,” *IEEE Transactions on Automatic Control*, vol. 57, no. 10, pp. 2656–2662, 2012.

[10] M. Coraggio, P. DeLellis, and M. di Bernardo, “Convergence and synchronization in networks of piecewise-smooth systems via distributed discontinuous coupling,” *Automatica*, vol. 129, pp. 1–9, 2021.

[11] V. Utkin, J. Guldner, and J. Shi, *Sliding Mode Control in Electro-Mechanical Systems*. CRC press, 2009, vol. 34.

[12] C. Godsil and G. Royle, *Algebraic Graph Theory*. New York, USA: Springer, 2001.

[13] R. T. Rockafellar and R. J.-B. Wets, *Variational Analysis*. Springer-Verlag, 2009.

[14] S. Boccaletti, G. Bianconi, R. Criado, C. I. del Genio, J. Gómez-Gardeñez, M. Romance, I. Sendiña Nadal, Z. Wang, and M. Zanin, “The structure and dynamics of multilayer networks,” *Physics Reports*, vol. 544, pp. 1–122, 2014.

[15] H. Attouch and A. Damlamian, “On multivalued evolution equations in Hilbert spaces,” *Israel Journal of Mathematics*, vol. 12, pp. 373–390, 1972.

[16] B. Brogliato and D. Goeleven, “Well-posedness, stability and invariance results for a class of multivalued Lur'e dynamical systems,” *Nonlinear Analysis: Theory, Methods and Applications*, vol. 74, no. 1, pp. 195–212, 2011.

[17] B. Baji and A. Cabot, “An inertial proximal algorithm with dry friction: finite convergence results,” *Set-Valued Analysis*, vol. 14, no. 1, pp. 1–23, 2006.

[18] R. A. Horn and C. R. Johnson, *Matrix Analysis*, 2nd ed. New York, USA: Cambridge University Press, 2013.

[19] A. V. Proskurnikov, F. Zhang, M. Cao, and J. M. A. Scherpen, “A general criterion for synchronization of incrementally dissipative nonlinear coupling agents,” in *European Control Conference (ECC)*, Linz, Austria, 2015, pp. 581–586.

[20] J. Bastien and M. Schatzman, “Numerical precision for differential inclusions with uniqueness,” *ESAIM: Mathematical Modelling and Numerical Analysis*, vol. 36, no. 3, pp. 427–460, 2002.

[21] M. E. J. Newman and D. J. Watts, “Renormalization group analysis of the small-world network model,” *Physics Letters A*, vol. 263, pp. 341–346, 1999.