# Robust synchronization of heterogeneous networks of nonlinear time-varying dynamical systems using set-valued linear complementarity couplings

Quang Hung Pham, Félix Miranda-Villatoro, Bernard Brogliato

► To cite this version:

Quang Hung Pham, Félix Miranda-Villatoro, Bernard Brogliato. Robust synchronization of heterogeneous networks of nonlinear time-varying dynamical systems using set-valued linear complementarity couplings. 2025. hal-05235089

HAL Id: hal-05235089

https://inria.hal.science/hal-05235089v1

Preprint submitted on 2 Sep 2025

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 of heterogeneous networks of
nonlinear time-varying dynamical systems using
set-valued linear complementarity couplings

Quang Hung Pham * Félix A. Miranda-Villatoro *

Bernard Brogliato *

September 2, 2025

Abstract

This work is largely concerned with the synchronization of nonlinear time-
varying agents, using static couplings designed as linear complementarity con-
straints (hence possibly set-valued operators). The closed-loop system takes the
form of linear complementarity systems with a specific structure. Closed-loop
well-posedness, equilibrium points, exponential and finite-time stability, and ro-
bustness with respect to uncertainties are carefully studied. Linear complemen-
tarity problems and variational inequalities are the main analysis tools. Pas-
sivity and maximal monotonicity of operators arguments are used. The time-
discretization and practical implementation are studied in detail. Several ex-
amples with simulations illustrate the theoretical findings. It is shown that the
proposed framework allows to consider new nonsmooth feedback controllers.

# 1 Introduction

Complex networks have emerged as a central topic in systems and control research over the past three decades, due to their widespread occurrence in technological [1, 2], natural [3, 4, 5], and social systems [6, 7, 8], see details in [9, 10, 11] and the references therein. It is pointed out in [9] that a complex network is typically composed of three main components:

* A collection of agents (nodes), each characterized by its individual dynamics.

* A coupling protocol (or coupling law) that define how agents are interconnected and how they influence each other.

<footnote n="1">¹Univ. Grenoble Alpes, INRIA, CNRS, Grenoble INP, LJK, 38000 Grenoble, France. quang-hung.pham@inria.fr; felix.miranda-villatoro@inria.fr; bernard.brogliato@inria.fr</footnote>
* A graph that describes the structure of the network by specifying which nodes are connected to one another.

One of the fundamental problems in the study of complex networks is the problem of synchronization, which, informally speaking, concerns designing the network (including both its structure and coupling protocols) such that all agents connected within the network asymptotically reach the same state or behavior over time, despite starting from different initial states. The study of synchronization traces back to Christiaan Huygens (1665), who observed that two pendulum clocks mounted on the same support would naturally synchronize their oscillations. This initial experimental insight inspired the development of mathematical models of coupled oscillators by A.T. Winfree [12] (1960s) and Y. Kuramoto [13] (1970s). In the 1990s, L.M. Pecora and T.L. Carroll were among the first to investigate synchronization in chaotic systems [14], establishing the foundational concepts for analyzing network synchronization. Later, D.J. Watts and S.H. Strogatz [15] examined how network structures, including small-world and complex networks, affect synchronization phenomena, providing a framework that has shaped later research on emergent behaviors in networks. A deep understanding of underlying synchronization mechanisms is essential, not only for advancing theoretical knowledge in nonlinear dynamics and control, but for designing robust and adaptive strategies that guarantee the desired performance and stability of complex networks in practice. Synchronization has widespread applications across a variety of domains, ranging from biology and neuroscience to engineering, computer science, and the social and economic sciences; a detailed discussion can be found in [16] and its cited references. In the majority of previous studies on the synchronization problem, agents in a network are commonly assumed to have individual dynamics described by the following input-output dynamical systems:

$$ \begin{cases} \dot{x}_i(t) = f_i(t, x_i(t)) + B_i u_i(t) \\ y_i(t) = C_i x_i(t) \end{cases}, \quad \forall i = 1, \dots, N, \quad (1) $$

where $x_i(t) \in \mathbb{R}^n$, $u_i(t) \in \mathbb{R}^q$, $y_i(t) \in \mathbb{R}^p$ are the state, the input and the output of agent $i$, respectively, $f_i : \mathbb{R}_+ \times \mathbb{R}^n \to \mathbb{R}^n$ is a single-valued function, $B_i, C_i$ are matrices of appropriate dimensions, and $N$ is the number of agents in the network.

When all $f_i, B_i$, and $C_i$ in (1) are identical for all $i = 1, \dots, N$, the network of agents (1) is referred to as a homogeneous network. Synchronization of such networks has been extensively investigated under various types of coupling laws and graph topologies, including static linear diffusive coupling law [17], dynamic coupling law [18], as well as static undirected connected graphs [19], weighted directed graphs [20], among others; for further details, see [21, 22, 23, 24, 25, 26, 27, 28, 29], and references therein. For instance, the case where $f_i \equiv f : \mathbb{R}_+ \times \mathbb{R}^n \to \mathbb{R}^n$ and $B_i = C_i = I_n$, for all $i = 1, \dots, N$, is considered in [19]. In that work, under a given connected and undirected network graph, it is proved that if $f$ is globally Lipschitz in the second variable with Lipschitz
constant $\alpha$, and the linear coupling gain is greater than $\max(\alpha, \alpha/\lambda_2)$ where $\lambda_2$ is the second smallest eigenvalue of the Laplacian matrix associated to the graph, then the network achieves global asymptotic synchronization. Another instance where $B_i$ and $C_i$ are identical, $f_i(t, x_i) = Ax_i$, for all $i = \overline{1, \overline{N}}$, for some matrix $A \in \mathbb{R}^{n \times n}$, is treated in [24]. In that context, a directed, balanced graph (i.e., the in-degree and the out-degree of each node are equal), and the linear diffusive couplings with time-varying gain are utilized for achieving network synchronization under the assumption that $(A, B, C)$ is observable and passive (i.e., $\exists P = P^\top > 0 : PA + A^\top P \le 0, PB = C^\top$).

In addition, several studies, such as [30, 31], achieve state synchronization via output synchronization. Studies on homogeneous network synchronization are still ongoing, aiming at a significant reduction of coupling energy cost and output noise sensitivity. For example, [32] propose a partial-state coupling with static and adaptive nonlinear protocols that achieve synchronization of homogeneous networks under mild conditions and further contribute to reducing coupling energy through avoiding detailed network information. In [33], and [34], the authors investigate nonlinear couplings expressed as definite integrals of nonlinear gain functions, showing that this approach achieves more efficient synchronization while requiring less input energy and exhibiting lower sensitivity to measurement noise than conventional linear couplings.

In practice, uncertainties and perturbations always occur and impact the agents in the network at arbitrary time. These may originate from the agents' individual dynamics, inaccuracies of outputs and inputs, or exogenous factors that are generally unknown. In this context, the network synchronization becomes highly challenging. Motivated by this practical challenge, studies in the last decade have almost entirely turned their focus to investigating synchronization in heterogeneous networks. The most of studies on such a network consider agents of the form (1) where $f_i$, $B_i$ and $C_i$ are not necessarily identical for all $i = \overline{1, \overline{N}}$. One of the earliest studies in this direction is [18], which investigated a network of agents (1) characterized by identical nominal linear dynamics, while allowing for uncertainties in the outputs. Specifically, the perturbed dynamics systems of agents in [18] are described by:

$$ \left\{ \begin{array}{l} \dot{x}_i(t) = Ax_i(t) + Bu_i(t), \\ y_i(t) = Cx_i(t) + d_i(t), \quad \forall i = 1, \dots, N, \\ z_i(t) = u_i(t), \\ d_i(t) = \Delta_i z_i(t), \end{array} \right. \qquad (2) $$

where $x_i(t) \in \mathbb{R}^n$, $u_i(t) \in \mathbb{R}^q$, $y_i(t) \in \mathbb{R}^p$, and $d_i(t) \in \mathbb{R}^p$ are the state, the input, the output and the perturbation of agent $i$, respectively; $A, B$, and $C$ are matrices of appropriate dimensions, $\Delta_i$ is a proper and stable rational transfer matrix, and $N$ is the number of agents in the network. Using an observer-based dynamic protocol and a connected undirected graph, under the assumption that $A$ has no eigenvalues on
the imaginary axis, [18] asserts that the heterogeneous network of agents (2) achieves global asymptotic synchronization for all $Δ_i$ such that $||Δ_i||_∞ ≤ η$ for all $i = 1, N$, and for some $η > 0$, where $||·||_∞$ denotes the infinity norm. The cases where the matrix A can have imaginary axis eigenvalues, and where additive perturbations are nonlinear, are also addressed in [18]. Another study is conducted in [35], where more complex conditions are considered to achieve global asymptotic synchronization.

Notably, it has been shown in [36] that the asymptotic synchronization problem can be rigorously reformulated as the problem of asymptotic stability of a synchronization manifold. While this manifold always exists in homogeneous networks, it may fail to exist or, if it does, it may be very difficult to verify and establish the asymptotic stability of such a manifold in heterogeneous networks. This leads to the majority of previous studies concentrating on practical synchronization (that is, all agents asymptotically converge to a same state or behavior within a small, acceptable bound) rather than on asymptotic synchronization. Extensive studies have been carried out on the practical synchronization of networks of agents (1). For instance, [37, 38, 39, 40] investigate (1) with $B_i = C_i = I_n$ for all $i = 1, N$; [36, 41] consider (1) where $B_i, C_i$ are identical for all $i = 1, N$ such that $C_iB_i$ is similar to a positive diagonal matrix (i.e, there exists $U$ such that $U^{-1}CBU$ is diagonal matrix with positive entries); and [42] considers (1) where $B_i$ or $C_i$ are identical for all $i = 1, N$. In all these works, the linear diffusive coupling is used; additionally, a nonlinear coupling is also proposed in [39]. Extensions to agent dynamics beyond the structure of (1) have been reported in [43, 44, 45].

According to [46], synchronization is currently a well-established research area; nevertheless, one of the biggest problems that still needs to be tackled is how to ensure synchronization in heterogeneous networks, where the heterogeneity may arise from external perturbations or uncertainties of the agents. In addition, as highlighted in [38, 41], the use of solely linear diffusive coupling protocols without local controllers makes asymptotic synchronization in heterogeneous networks unattainable. These raise a critical design challenge: to construct a nonlinear coupling protocol under which the heterogeneous network achieves asymptotic synchronization.

To the best of the authors' knowledge, existing research works addressing the design of nonsmooth nonlinear coupling protocols that achieve asymptotic synchronization in heterogeneous networks of agents (1) are limited to [47], [48], and [49]. In particular, in [47], a distributed discontinuous coupling protocol is proposed, constructed by adding a discontinuous signum-function-based coupling layer to a traditional linear diffusive one. Using this coupling protocol, the network of agents (1) where $B_i = C_i = I_n$ for all $i = 1, N$, achieves global asymptotic synchronization under the assumptions that the network is uniformly ultimately bounded to a ball $rB$ for some $r > 0$, $f_i$ satisfies the QUAD condition in the ball $rB$ (a generalization of Lipschitz continuity and one-sided Lipschitz continuity) for all $i = 1, N$, the network graph is connected and undirected,
and the coupling gains exceed certain thresholds determined by the agents dynamics, the network topology, the mismatch between individual and average vector fields, and the so-called minimum density of the network graph.

Inspired by [47], the author in [49] proposed a class of nonsmooth set-valued coupling laws (or coupling protocols) characterized by maximal monotone operators, which directly extends the coupling protocols used in [47] by replacing the discontinuous signum-function-based coupling with one based on a maximal monotone operator. In this context, the agents (1) are considered with $f_i(t, x_i) = A_0x_i - A_1\varphi(A_2x_i) + B\xi_i(t, x_i)$, $B_i = B$, and $C_i = C$ for all $i = \overline{1, N}$, for some matrices $A_0, A_1, A_2, B, C$ of appropriate dimensions, and Lipschitz functions $\varphi: \mathbb{R}^m \to \mathbb{R}^m$, and $\xi_i: \mathbb{R}_+ \times \mathbb{R}^n \to \mathbb{R}^q$. It is shown in [49] that networks of such agents achieve global asymptotic (resp. finite-time) synchronization under an incremental sector condition on $\varphi$, a dry friction condition, technical conditions related to passivity, and threshold conditions on the coupling gains. By following a different approach inspired by the so-called funnel controller, a class of nonlinear coupling laws referred to as the funnel coupling laws, is introduced in [48]. Under this type of coupling law, [48] established sufficient conditions that guarantee asymptotic synchronization for networks of (1) where $B_i = C_i = I_n$ for all $i = \overline{1, N}$. The problem of finite-time synchronization is also addressed in [48].

Although nonsmooth and nonlinear coupling protocols have been proposed to achieve asymptotic synchronization in heterogeneous networks, the practical implementation of such coupling protocols under unknown perturbations remains an open problem. The difficulties of this problem arise from the following two factors:

* The discontinuity of nonsmooth nonlinear coupling protocols.

* The extreme sensitivity of the set-valued feedback controllers, as reported in [50].

For instance, sign-coupling protocols used in [47], while theoretically guaranteeing that networks of (1) achieve asymptotic synchronization, encounter difficulties in practical implementation due to their discontinuity at the origin. In particular, perfect regulation requires knowing the exact values (selections) of the control $u_i(t)$ when $y_i(t) - y_j(t) = 0$ for all $i, j = \overline{1, N}$, in order to compensate for heterogeneities in the network and to maintain synchronization. However, the discontinuity of the signum function at the origin prevents this, leading to incorrect control selections and causing it to exhibit potentially hazardous chattering.

In this work, we present a number of significant advancements over prior work, which can be summarized as follows:

1. We consider an extended class of agents, in the sense that their dynamics directly generalize those given in (1) by allowing the presence of a feedthrough term (or direct feedthrough term) in the output of agents.
2. We propose a novel class of nonsmooth nonlinear set-valued coupling laws, which enables the heterogeneous networks to achieve global exponential (finite-time) synchronization. This class has been shown to encompass both the conventional linear diffusive coupling protocol, and the sign-coupling protocol used in [47]. It is also a particular subclass of the class of coupling laws introduced in [49].

3. We provide algorithms specifically designed to compute discrete-time controllers within a practical implementation framework, allowing the controllers to avoid exhibiting hazardous chattering. This is crucial in practical applications.

**Notation:** Let $\langle \cdot, \cdot \rangle$ denote the Euclidean inner product. Denote by $\| \cdot \|$, and $\| \cdot \|_1$ the spectral norm and the 1-norm (or column norm), respectively. Let $S \subseteq \mathbb{R}^m$ be a nonempty set, $S^0$ is its least-norm element (if it exists), $\text{int}(S)$ is its interior, $\text{rint}(S)$ is its relative interior, $\text{cl}(S)$ is its closure, $\text{aff}(S)$ is its affine hull, and $\text{proj}_S(\cdot)$ is the projection map on $S$. Let $S^* = \{y \in \mathbb{R}^m : \langle y, x \rangle \ge 0, \forall x \in S\}$. Let $f : \mathbb{R}^n \to \mathbb{R}^m$ be a function, $\text{dom } f = \{x \in \mathbb{R}^n : f(x)$ is defined$\}$ is its domain. Let $F : \mathbb{R}^n \rightrightarrows \mathbb{R}^m$ be a set-valued mapping (or operator), $\text{Dom } F = \{x \in \mathbb{R}^n : F(x) \neq \emptyset\}$ is its domain, $\text{rge } F = \{y \in \mathbb{R}^m : \exists x \in \mathbb{R}^n \text{ with } y \in F(x)\}$ is its range, $\text{Graph } F = \{(x, y) \in \mathbb{R}^n \times \mathbb{R}^m : y \in F(x)\}$ is its graph, and $F^{-1} : \mathbb{R}^m \rightrightarrows \mathbb{R}^n$ is its inverse mapping. Let $M \in \mathbb{R}^{m \times n}$ be a matrix, $M^\top$ is its transpose, $M^\dagger$ is its pseudo-inverse, $\ker M$ is its null space, $\text{rge } M$ is its range, $\lambda_{\min}(M)$ and $\lambda_{\max}(M)$ are its smallest and largest eigenvalues (if $m = n$), respectively. Denote by $\text{diag}(D_1, D_2, \dots, D_n)$ the block diagonal matrix in $\mathbb{R}^{nm \times nm}$ with $D_1, D_2, \dots, D_n \in \mathbb{R}^{m \times m}$ the elements on the diagonal. The Kronecker product is denoted as $\otimes$. A matrix $D \in \mathbb{R}^{m \times m}$ is positive definite (resp. positive semi-definite), denoted by $D \succ 0$ (resp. $D \succcurlyeq 0$) if $x^\top Dx > 0$ (resp. $x^\top Dx \ge 0$) for all $x \in \mathbb{R}^m \setminus \{0\}$. A matrix $D$ is copositive-plus if $x^\top Dx \ge 0$ for all $x \in \mathbb{R}_+^m$ and the following implication holds: $[x^\top Dx = 0, x \ge 0] \Rightarrow [(D + D^\top)x = 0]$. LCS, TVLCS, TVCS, LCP are for linear complementarity systems, time-varying linear complementarity systems, time-varying complementarity systems, and linear complementarity problems, respectively. Given $q \in \mathbb{R}^m$ and $D \in \mathbb{R}^{m \times m}$, LCP($q$, $D$) denotes the problem: find $\lambda$ such that: $0 \le \lambda \perp D\lambda + q \ge 0$, and SOL($q$, $D$) is its solution set. We denote $\mathbf{0}_n = (0, 0, \dots, 0)^\top \in \mathbb{R}^n$ and $\mathbf{1}_n = (1, 1, \dots, 1)^\top \in \mathbb{R}^n$. Denote by $I_n$ the identity matrix in $\mathbb{R}^{n \times n}$. The symbol $\subsetneq$ denotes a strict subset, and $\subsetneq, \subset$ denote subsets without distinguishing whether the sets are equal or not. Denote by $\mathcal{B}_\epsilon = \{\gamma \in \mathbb{R}^n : \|\gamma\| \le \epsilon\}$ the closed ball centered at the origin with radius $\epsilon$ in $\mathbb{R}^n$, $\mathcal{B}(x, \epsilon) = \{\gamma \in \mathbb{R}^n : \|\gamma - x\| \le \epsilon\}$ the closed ball centered at $x$ with radius $\epsilon$ in $\mathbb{R}^n$.
# 1.1 Preliminaries

## 1.1.1 Basic concepts of graph theory

The following concepts can be found in [51, 52]. An undirected graph $\mathcal{G}$ is an ordered pair $(\mathcal{V}, \mathcal{E})$, denoted as $\mathcal{G} = (\mathcal{V}, \mathcal{E})$, where $\mathcal{V}$ is a set of vertices and $\mathcal{E}$ is a set of unordered pairs of vertices, called edges (undirected). Denote by $|\mathcal{V}|, |\mathcal{E}|$ the number of elements of $\mathcal{V}, \mathcal{E}$, respectively. For each $k \in \{1, ..., |\mathcal{E}|\}$, the edge $\epsilon_k$ in $\mathcal{E}$ is uniquely defined by an unordered pair of vertices, i.e., there exists a unique pair $(i, j) \equiv (j, i)$ where $i, j \in \mathcal{V}$, such that $\epsilon_k = (i, j)$. If several edges are created by the same pairs of vertices, then they are called multiple edges. An edge with identical endpoints is called a loop. A simple graph is a graph without loops and multiple edges. A path from vertex $i$ to vertex $j$, for some $i, j \in \mathcal{V}$, is an ordered sequence of vertices with starting point $i$, ending point $j$ and any pair of consecutive vertices in the sequence is an edge of the graph. A graph is connected if there exists a path between any two vertices. On each edge $\epsilon_k \in \mathcal{E}$ of the undirected graph $\mathcal{G}$, for some $k \in \{1, ..., |\mathcal{E}|\}$, its end points are assigned by + (positive) or - (negative). This allows us to define an oriented incidence matrix $\Theta = [\Theta_{ik}] \in \mathbb{R}^{|\mathcal{V}| \times |\mathcal{E}|}$ associated with $\mathcal{G}$, which is given by the following, with $i \in \mathcal{V}, \epsilon_k \in \mathcal{E}$:

$$ \Theta_{ik} = \begin{cases} +1 & , \text{ if } i \text{ is assigned by } + \text{ on } \epsilon_k, \\ -1 & , \text{ if } i \text{ is assigned by } - \text{ on } \epsilon_k, \\ 0 & , \text{ otherwise.} \end{cases} $$

All oriented incidence matrices of a connected graph have rank $|\mathcal{V}| - 1$. An example of such matrices is illustrated in Figure 1.

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

$$ \Theta = \begin{bmatrix} \epsilon_1 & \epsilon_2 & \epsilon_3 & \epsilon_4 & \epsilon_5 \\ 1 & 1 & 1 & 0 & 0 \\ -1 & 0 & 0 & 1 & 0 \\ 0 & -1 & 0 & -1 & 1 \\ 0 & 0 & -1 & 0 & -1 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix} \begin{matrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \end{matrix} $$

<caption>Figure 1: (Left) An undirected graph with an orientation selected. (Right) Associated oriented incidence matrix $\Theta$.</caption>

## 1.1.2 Maximal monotone and proximal operators

Let us recall some useful definitions, see [53, 54, 55, 56, 57].

**Definition 1.** A set-valued operator $F: \mathbb{R}^m \to \mathbb{R}^m$ is:

(i) monotone if $\langle u - v, x - y \rangle \ge 0, \forall (x, u), (y, v) \in \text{Graph } F$.
(ii) maximal monotone if its graph is not properly contained in the graph of any other monotone map.

(ii) uniformly monotone with modulus $\Phi: \mathbb{R}_+ \to [0, +\infty]$ if $\Phi$ is increasing, vanishes only at 0, and

$$ \langle u - v, x - y \rangle \geq \Phi(\|x - y\|), \quad \forall (x, u), (y, v) \in \text{Graph } F. $$

(iii) strongly monotone with constant $\beta > 0$ if

$$ \langle u - v, x - y \rangle \geq \beta \|x - y\|^2, \quad \forall (x, u), (y, v) \in \text{Graph } F. $$

The operator $F$ is said to be uniformly (resp. strongly) maximal monotone if $F$ is maximal monotone and uniformly (resp. strongly) monotone.

**Definition 2.** Let $f: \mathbb{R}^m \to (-\infty, +\infty]$ be proper (i.e., $-\infty < f \not\equiv +\infty$). The subdifferential of $f$ is the set-valued operator $\partial f: \mathbb{R}^m \rightrightarrows \mathbb{R}^m$,

$$ x \mapsto \partial f(x) := \{u \in \mathbb{R}^m : f(y) \geq f(x) + \langle u, y - x \rangle, \forall y \in \mathbb{R}^m \}. $$

**Definition 3.** A function $f: \mathbb{R}^m \to [-\infty, +\infty]$ is lower semicontinuous at $\bar{x}$ if $f(\bar{x}) \le \liminf_{x \to \bar{x}} f(x)$, and lower semicontinuous on $\mathbb{R}^m$ if this holds for every $\bar{x} \in \mathbb{R}^m$.

**Definition 4.** Let $S \subseteq \mathbb{R}^m$ be a convex nonempty set. The normal cone to $S$ at $x \in S$ is defined by $N_S(x) := \{p \in \mathbb{R}^m : \langle p, x - y \rangle \ge 0, \forall y \in S\}$.

If $S := \mathbb{R}_+^m$, then the operator $N_{\mathbb{R}_+^m}$ is maximal monotone, and for any $x, y \in \mathbb{R}^m$:

$$ 0 \le x \perp y \ge 0 \quad \Leftrightarrow \quad x \in -N_{\mathbb{R}_+^m}(y) \quad \Leftrightarrow \quad y \in -N_{\mathbb{R}_+^m}(x). \quad (3) $$

**Definition 5.** Let $F: \mathbb{R}^m \rightrightarrows \mathbb{R}^m$ be a set-valued operator. The operator $J_F: \mathbb{R}^m \rightrightarrows \mathbb{R}^m$ defined by $J_F := (Id + F)^{-1}$, is called the resolvent of $F$, where $Id: \mathbb{R}^m \to \mathbb{R}^m$ denotes the identity mapping. Let $f: \mathbb{R}^m \to (-\infty, +\infty]$ be a proper function. The proximal operator $\text{prox}_f: \mathbb{R}^m \rightrightarrows \mathbb{R}^m$ of $f$ is defined by $\text{prox}_f := (Id + \partial f)^{-1}$.

If $F$ is a maximal monotone operator, then $J_F$ is a single-valued map with $\text{Dom}(J_F) = \mathbb{R}^m$, and maximal monotone. If $f: \mathbb{R}^m \to (-\infty, +\infty]$ be a proper lower semicontinuous convex function, then the following hold:

1. $\partial f$ is maximal monotone.

2. $\text{prox}_f$ is single-valued mapping with $\text{Dom}(\text{prox}_f) = \mathbb{R}^m$, and maximal monotone.

3. For all $\gamma > 0$, for all $x \in \mathbb{R}^m$: $x = \text{prox}_{\gamma f}(x) + \gamma \text{prox}_{\frac{1}{\gamma} f^*} \left( \frac{x}{\gamma} \right)$.

4. For all $x \in \mathbb{R}^m$: $\text{prox}_f(x) = \arg\min_{y \in \mathbb{R}^m} \left( f(y) + \frac{1}{2} \|y - x\|^2 \right)$.

**Remark 1.** It is known that passivity and monotonicity are closely related [55, 58, 53, 59, 60, 61, 62]. Conditions which are reminiscent from passivity are used throughout this work to enable us to use maximal monotonicity arguments, see (18), (38), (47), (124), item 3 in Proposition 6, item 2.2 in Theorem 3, and items 2.2 and 2.3 in Theorem 5.
### 1.1.3 QUAD Condition

**Definition 6.** (QUADness) A vector field $f: \mathbb{R}_+ \times \mathbb{R}^n \to \mathbb{R}^n$ is said to be QUAD($P, Q$) in $\Omega \subseteq \mathbb{R}^n$ if there exist symmetric matrices $P, Q \in \mathbb{R}^{n \times n}$ such that:

$$ (x - y)^{\top} P (f(t, x) - f(t, y)) \leq (x - y)^{\top} Q (x - y), \quad \forall x, y \in \Omega, \forall t \in \mathbb{R}_+. $$

The QUAD condition has been widely used in previous studies on network synchronization (e.g., [63], [?], [64], [19], [65], [66], and references therein). Let $\mathcal{OL}$ denote the class of One-sided Lipschitz continuous functions, i.e., the set of all functions $f: \mathbb{R}^n \to \mathbb{R}^n$ that satisfy: $\langle x - y, f(x) - f(y) \rangle \leq \nu \|x - y\|^2$, for all $x, y \in \mathbb{R}^n$, where $\nu \in \mathbb{R}$ is called the one-sided Lipschitz constant of $f$ (see [67]). Clearly the class of QUAD functions, denoted by $\mathcal{Q}$, is a natural extension of $\mathcal{OL}$, which in turn extends the class of Lipschitz continuous functions $\mathcal{L}$. Thus, the inclusion relations can be expressed as:

$$ \mathcal{L} \subseteq \mathcal{OL} \subseteq \mathcal{Q}. $$

The sufficient conditions for a differentiable function $f$ to be QUAD are provided in [47, Proposition 2], and is recalled below.

**Proposition 1.** If a function $f: \mathbb{R}^n \to \mathbb{R}^n$ has an upper bounded Jacobian in $\Omega \subseteq \mathbb{R}^n$, in the sense that for all $x \in \Omega$:

$$ \frac{\partial f_i(x)}{\partial x_i} \leq \eta_{ii}, \quad \text{and} \quad \left| \frac{\partial f_i(x)}{\partial x_j} \right| \leq \eta_{ij}, \quad \forall i \neq j, i, j = 1, \dots, n, $$

then $f$ is QUAD($I_n, \mathcal{Q}$) in $\Omega$ with $\mathcal{Q}$ being diagonal and $Q_{ii} = \eta_{ii} + \sum_{j=1, j \neq i}^{n} \frac{\eta_{ij} + \eta_{ji}}{2}$.

### 1.1.4 Linear Complementarity Problems (LCP), Linear Complementarity Systems (LCS), and Time-Varying Linear Complementarity Systems (TVLCS)

Let us recall the definitions of LCP, LCS, TVLCS and the fundamental result of LCP theory [56].

**Definition 7.** Given a matrix $Q \in \mathbb{R}^{n \times n}$ and a vector $q \in \mathbb{R}^n$. The linear complementarity problem (LCP) is to find a vector $z \in \mathbb{R}^n$ such that $0 \leq z \perp Qz + q \geq 0$, that is,

$$ z \geq 0, \quad Qz + q \geq 0, \quad z^{\top}(Qz + q) = 0, $$

or to show that no such vector $z$ exists. This problem is abbreviated LCP($q, \mathcal{Q}$). Its solution set is denoted as SOL($q, \mathcal{Q}$).

**Theorem 1.** The LCP($q, \mathcal{Q}$) has a unique solution for all $q \in \mathbb{R}^n$ if and only if $\mathcal{Q}$ is P-matrix (i.e., all its principal minors are positive). Moreover, such a solution is a piecewise linear function in $q \in \mathbb{R}^n$.
Note that every positive definite matrix is a **P-matrix**, and a symmetric **P-matrix** is positive definite. Some key results on the existence of solutions for LCP (which may not be unique) can be found in [56, Theorem 3.8.6, Corollary 3.8.10]. Let us now introduce LCS and TVLCS.

**Definition 8.** A Linear Complementarity System (LCS) is defined as:

$$ egin{cases} 
\dot{x}(t) = Ax(t) + B\lambda(t) + Eu(t), \\
w(t) = Cx(t) + D\lambda(t) + Fu(t), \\
0 \le \lambda(t) \perp w(t) \ge 0,
\end{cases} 
\quad (4) $$

where $x(t) \in \mathbb{R}^n$, $w(t), \lambda(t) \in \mathbb{R}^m$, $u(t) \in \mathbb{R}^u$, and the matrices $A, B, C, D, E, F$ are constant with suitable dimensions. A Time-Varying Linear Complementarity System (TVLCS) is an extension of LCS, defined as:

$$ egin{cases} 
\dot{x}(t) = A(t)x(t) + B(t)\lambda(t) + E(t)u(t), \\
w(t) = C(t)x(t) + D(t)\lambda(t) + F(t)u(t), \\
0 \le \lambda(t) \perp w(t) \ge 0,
\end{cases} $$

where $x(t) \in \mathbb{R}^n$, $w(t), \lambda(t) \in \mathbb{R}^m$, $u(t) \in \mathbb{R}^u$, and the matrix functions $A(t), B(t), C(t), D(t), E(t), F(t)$ are constant with suitable dimensions.

An extension of LCS is LCCS (for Linear Cone Complementarity Systems), in which the complementarity constraint is given by $\mathcal{K} \in \lambda(t) \perp w(t) \in \mathcal{K}^*$, where $\mathcal{K} \subseteq \mathbb{R}^m$ is a closed convex cone. Meanwhile, TVLCCS (for Time-Varying Linear Cone Complementarity Systems) extends TVLCS by defining the complementarity constraint as $\mathcal{K}(t) \in \lambda(t) \perp w(t) \in \mathcal{K}^*(t)$, where $\mathcal{K}(t) \subseteq \mathbb{R}^m$ is a time-varying closed convex cone. A summary of well-posedness results for LCS and LCCS is available in [68, 69], see also [70, 71]. The well-posedness of TVLCS and TVLCCS is given in [69]. When the linear term $Ax$ in (4) is replaced by $A(t, x)$, where $A$ may depend on time $t$ and be nonlinear in $x$, we simply speak of a Time-Varying Complementarity System (TVCS).

## 1.2 Network model description

In this work, we concentrate on studying networks whose structure is represented by a connected, simple, undirected graph $\mathcal{G} = (\mathcal{V}, \mathcal{E})$, and their agents are dynamical systems of the form:

$$ egin{cases} 
\dot{x}_i(t) = f_i(t, x_i(t)) + B_i u_i(t), \\
y_i(t) = C_i x_i(t) + D_i u_i(t), 
\end{cases} 
\quad \forall i \in \{1, \dots, |\mathcal{V}|\}, 
\quad (5) $$

in which, $x_i(t) \in \mathbb{R}^n$, $u_i(t) \in \mathbb{R}^q$, $y_i(t) \in \mathbb{R}^p$ are the state, input and output of the agent (vertex) $i$ in the network, respectively, $f_i: [0, +\infty) \times \mathbb{R}^n \to \mathbb{R}^n$ is a single-valued function, and $B_i, C_i, D_i$ are matrices with compatible dimensions.
**Remark 2.** In this work we consider possible nonzero $D_i$ matrices. To the best of authors' knowledge this is an original contribution.

The novel couplings formed by the LCPs structured as the following are used to connect the agents in a specific configuration, where $v_{ij}(t)$ and $\eta_{ij}(t)$ play the role of the input and the output of the coupling that connects agents $i$ and $j$, respectively.

$$ \eta_{ij}(t) = H_{ij}\lambda_{ij}(t), \text{ for } i < j, (i, j) \in \mathcal{E}, \qquad (6a) $$

$$ w_{ij}(t) = M_{ij}\lambda_{ij}(t) + N_{ij}v_{ij}(t) + r_{ij}, \qquad (6b) $$

$$ 0 \le \lambda_{ij}(t) \perp w_{ij}(t) \ge 0, \qquad (6c) $$

where $\lambda_{ij}(t), w_{ij}(t) \in \mathbb{R}^m$ are called complementarity variables, $H_{ij}, M_{ij}, N_{ij}$ are matrices that have compatible dimensions, and $r_{ij} \in \mathbb{R}^m$ is a constant vector. Let us now describe the interconnection law as follows:

$$ \left\{ \begin{array}{l} u_i(t) = \sum_{j \in \mathcal{N}_i} u_{ij}(t) \\ u_{ij}(t) = -u_{ji}(t) = \eta_{ij}(t), \\ v_{ij}(t) = y_i(t) - y_j(t), \end{array} \right. \qquad (7) $$

where $\mathcal{N}_i \subset \{1, ..., |\mathcal{V}|\}$ is the set of neighbours of agent $i$ in $\mathcal{G}$. Finally, the overall dynamics of the network is achieved by combining the three systems (5), (6), and (7). Figure 2 depicts a subsystem of the overall dynamical system of the network that describes the interaction between two agents $i$, and $j$ through their coupling.

<figure><img src="image_3.png" /><figcaption>Figure 2: A subsystem of the dynamical system with nonsmooth nonlinearities in feedback.</figcaption></figure>

## 1.3 Overall dynamics in compact form

Let us now consider the dynamical systems in (5). Let $x(t) := [x_1(t)^\top \quad x_2(t)^\top \quad \dots \quad x_{|\mathcal{V}|}(t)^\top]^\top \in \mathbb{R}^{n|\mathcal{V}|}$, $F$ be a single-valued function given by:

$$ F(\cdot, \cdot) : [0, +\infty) \times \mathbb{R}^{n|\mathcal{V}|} \rightarrow \mathbb{R}^{n|\mathcal{V}|} \\ (t, x) \quad \mapsto F(t, x) := [f_1(t, x_1)^\top \quad f_2(t, x_2)^\top \quad \dots \quad f_{|\mathcal{V}|}(t, x_{|\mathcal{V}|})^\top]^\top, \qquad (8) $$

$u(t) := [u_1(t)^\top \quad u_2(t)^\top \quad \dots \quad u_{|\mathcal{V}|}(t)^\top]^\top \in \mathbb{R}^{q|\mathcal{V}|}$, $y(t) := [y_1(t)^\top \quad y_2(t)^\top \quad \dots \quad y_{|\mathcal{V}|}(t)^\top]^\top \in \mathbb{R}^{p|\mathcal{V}|}$, and $B := \text{diag}(B_1, \dots, B_{|\mathcal{V}|})$ be a block diagonal matrix with all matrices $B_i$ on
the main diagonal. A similar notation holds for the block diagonal matrices $C := \text{diag}(C_1, ..., C_{|v|})$, and $D := \text{diag}(D_1, ..., D_{|v|})$. Then the overall dynamics of agents is rewritten as:

$$ \begin{cases} \dot{x}(t) = F(t, x(t)) + Bu(t), \\ y(t) = Cx(t) + Du(t). \end{cases} \qquad (9) $$

Let us consider the couplings in (6). The number of edges of $\mathcal{G}$ is the number of couplings in the network. Let $\eta(t) := [\eta_{i_1 j_1}(t)^\top \cdots \eta_{i_2 j_2}(t)^\top]^\top \in \mathbb{R}^{q|\mathcal{E}|}$ for $i_1 < j_1$ and $i_2 < j_2$ such that $(i_1, j_1), (i_2, j_2) \in \mathcal{E}$, be the stack vector of all vectors $\eta_{ij}(t)$. A similar notation following the order of indices in $\eta(t)$, applies to the stack vectors $v(t), \lambda(t), w(t), r$. Then the structure of all couplings in (6) is rewritten in compact form as:

$$ \left\{ \begin{array}{l} \eta(t) = H\lambda(t), \\ w(t) = M\lambda(t) + Nv(t) + r, \\ 0 \le \lambda(t) \perp w(t) \ge 0, \end{array} \right. \qquad (10) $$

in which, $H \in \mathbb{R}^{q|\mathcal{E}| \times m|\mathcal{E}|}$ is a block diagonal matrix where all matrices $H_{ij}$ on the main diagonal in an order consistent with the order of $\eta_{ij}(t)$ in the form of $\eta(t)$. A similar notation holds for the block diagonal matrices $M \in \mathbb{R}^{m|\mathcal{E}| \times m|\mathcal{E}|}$, $N \in \mathbb{R}^{m|\mathcal{E}| \times p|\mathcal{E}|}$. Under the above settings, the interconnection laws in (7) are rewritten in compact form as:

$$ \left\{ \begin{array}{l} u(t) = (\Theta \otimes I_q)\eta(t), \\ v(t) = (\Theta^\top \otimes I_p)y(t), \end{array} \right. \qquad (11) $$

where $I_q \in \mathbb{R}^{q \times q}$, $I_p \in \mathbb{R}^{p \times p}$ are identity matrices, and $\Theta$ is the oriented incidence matrix of $\mathcal{G}$, characterized by the fact that in each column of $\Theta$ the position of 1 is always higher than the position of -1, that is, the following holds for all $k \in \{1, ..., |\mathcal{E}|\}$:

$$ [\Theta_{ik} = 1, \Theta_{jk} = -1] \quad \Rightarrow \quad [i < j]. \qquad (12) $$

It follows from (9), (10), (11) that the overall dynamics of the coupled network is the TVCS:

$$ \dot{x}(t) = F(t, x(t)) + B(\Theta \otimes I_q)H\lambda(t), \qquad (13a) $$

$$ w(t) = N(\Theta^\top \otimes I_p)Cx(t) + [N(\Theta^\top \otimes I_p)D(\Theta \otimes I_q)H + M]\lambda(t) + r, \qquad (13b) $$

$$ 0 \le \lambda(t) \perp w(t) \ge 0. \qquad (13c) $$

Using the equivalent form of the complementarity constraint given in (3), (13b) and (13c) are rewritten as follows:

$$ \lambda(t) \in - \left( \mathcal{N}_{\mathbb{R}_+^{m|\mathcal{E}|-r}}^{-1} + \Gamma \right)^{-1} \left( N(\Theta^\top \otimes I_p)Cx(t) \right), \qquad (14) $$

where $\Gamma := N(\Theta^\top \otimes I_p)D(\Theta \otimes I_q)H + M$. Let us introduce the operator $\mathbf{\Phi}$, known as a Lur'e operator [53, 55, 61, 69]:

$$ \mathbf{\Phi}: \mathbb{R}^{m|\mathcal{E}|} \to \mathbb{R}^{m|\mathcal{E}|}, \quad v \mapsto \left( \mathcal{N}_{\mathbb{R}_+^{m|\mathcal{E}|-r}}^{-1} + \Gamma \right)^{-1}(v). \qquad (15) $$
Substitute $\lambda(t)$ into (13a), then (13) is equivalent to the following differential inclusion (DI):

$$ \dot{x}(t) \in F(t, x(t)) - B(\Theta \otimes I_q)H\Phi(N(\Theta^\top \otimes I_p)Cx(t)). \quad (16) $$

**Remark 3.** When $\Gamma \succ 0$, $\Phi$ is a single-valued function with domain equal to $\mathbb{R}^{m|\mathcal{E}|}$, and Lipschitz continuous [61, Proposition 1]. When $\Gamma \succcurlyeq 0$, $\Phi$ is a maximal monotone operator [72, Lemma 1] with domain equal to $(\text{SOL}(0, \Gamma))^* - r$, using [69, Lemma 2.3]. In addition, using [55, p.1054], [73, (2.4.14)-p.95], and [54, Example 6.23, Corollary 11.25] then

$$ (\text{SOL}(0, \Gamma))^* = \mathbb{R}_+^{m|\mathcal{E}|} + \text{rge}(\Gamma + \Gamma^\top) + \{\Gamma^\top \gamma \mid \gamma \in \mathbb{R}_+^{m|\mathcal{E}|}\}. \quad (17) $$

If $\Gamma = \Gamma^\top$, then this reduces to $\text{Dom } \Phi = \mathbb{R}_+^{m|\mathcal{E}|} + \text{rge}(\Gamma) - r$ and we recover [74].

## 2 Well-posedness of the TVCS (13) and Problem Formulation

The following theorem guarantees the well-posedness of (13), which is obtained by mimicking a transformation method for LCS proposed in [72, Section 3.2] to rewrite (16) as a DI of a special form, and then applies [75, Theorem 2.1] to this DI. Let us recall that $F(t, x) := [f_1(t, x_1)^\top \ f_2(t, x_2)^\top \ \dots \ f_{|V|}(t, x_{|V|})^\top]^\top$ where $x = [x_1^\top \ x_2^\top \ \dots \ x_{|V|}^\top]^\top \in \mathbb{R}^{n|V|}$.

**Theorem 2.** Suppose that for each $i = 1, |V|$, $f_i : \mathbb{R}_+ \times \mathbb{R}^n \to \mathbb{R}^n$ is locally integrable on $\mathbb{R}_+$ in the first variable and Lipschitz continuous in the second variable, $\Gamma \succcurlyeq 0$, and the following hold:

1. There exists $P \in \mathbb{R}^{n|V| \times n|V|}$ such that $P = P^\top \succ 0$, and

$$ \ker (\Gamma + \Gamma^\top) \subseteq \ker (PB(\Theta \otimes I_q)H - C^\top(\Theta \otimes I_p)N^\top), \quad (18) $$

2. $\text{rge}(N(\Theta^\top \otimes I_p)C) \cap \text{rint}((\text{SOL}(0, \Gamma))^* - r) \neq \emptyset$,

where $\Gamma = N(\Theta^\top \otimes I_p)D(\Theta \otimes I_q)H + M$, and $(\text{SOL}(0, \Gamma))^*$ defined by (17). Then, (13) has a unique absolutely continuous solution for any initial condition $x(0) = x_0 \in (N(\Theta^\top \otimes I_p)C)^{-1}((\text{SOL}(0, \Gamma))^* - r)$.

*Proof.* Since $\Gamma \succ 0$, and items 1, 2 hold, mimicking [72, Section 3.2] (or [69, Section 2.2.1]), the DI (16) is rewritten as follows:

$$ \dot{z}(t) \in h(t, z(t)) - \Psi(z(t)), \quad (19) $$

where $z(t) := Rx(t)$ with $R = R^\top \succ 0$ is a square root of $P$ in item 1, and $\Psi(\cdot)$ is maximal monotone defined as:

$$ \Psi(z(t)) := R^{-1}C^\top(\Theta \otimes I_p)N^\top\Phi(N(\Theta^\top \otimes I_p)CR^{-1}z(t)), $$
$h(t, z(t)) := RF(t, R^{-1}z(t)) + (R^{-1}C^{\top}(\Theta \otimes I_p)N^{\top} - RB(\Theta \otimes I_q)H) \lambda^{\text{im}}(R^{-1}z(t))$, with $\lambda^{\text{im}}(\cdot) = \text{proj}_{\text{rge}(\Gamma+\Gamma^{\top})}(\mathbf{\Phi}(N(\Theta^{\top} \otimes I_p)C\cdot))$ is a Lipschitz continuous function, and $\text{dom } h(t, \cdot) = \text{Dom } \mathbf{\Psi} = R(\text{dom } \lambda^{\text{im}}) = R(N(\Theta^{\top} \otimes I_p)C)^{-1}(\text{Dom } \mathbf{\Phi})$ is a closed set, see [69]. Let us check conditions of Theorem B1 (a special case of [75, Theorem 2.1]) on (19). Firstly, since $f_i(\cdot, x)$ is measurable for each fixed $x \in \mathbb{R}^n$, for all $i = \overline{1, |\mathcal{V}|}$, $F(\cdot, R^{-1}z)$ is measurable, therefore $h(\cdot, z)$ is measurable for each fixed $z \in \mathbb{R}^n$. In addition, $h(t, \cdot)$ is Lipschitz continuous for each fixed $t \in \mathbb{R}_+$ and from there, measurability follows. Indeed, for each fixed $t \in \mathbb{R}_+$, for any $\alpha = [\alpha_1^{\top} \dots \alpha_{|\mathcal{V}|}^{\top}]^{\top}$, $\beta = [\beta_1^{\top} \dots \beta_{|\mathcal{V}|}^{\top}]^{\top} \in \text{dom } h(t, \cdot)$, we have:

$\|h(t, \alpha) - h(t, \beta)\|_1 \le \|R\|_1 \|F(t, R^{-1}\alpha) - F(t, R^{-1}\beta)\|_1 + \|\hat{R}\|_1 \|\lambda^{\text{im}}(R^{-1}\alpha) - \lambda^{\text{im}}(R^{-1}\beta)\|_1$, where $\hat{R} = R^{-1}C^{\top}(\Theta \otimes I_p)N^{\top} - RB(\Theta \otimes I_q)H$. Let $c_{ij} \in \mathbb{R}^n$ be a vector in the matrix

of the form $R^{-1} = \begin{bmatrix} c_{11}^{\top} & \dots & c_{1|\mathcal{V}|}^{\top} \\ \vdots & \ddots & \vdots \\ c_{n|\mathcal{V}|1}^{\top} & \dots & c_{n|\mathcal{V}||\mathcal{V}|}^{\top} \end{bmatrix} \in \mathbb{R}^{n|\mathcal{V}| \times n|\mathcal{V}|}$. Then

$$R^{-1}\alpha = \left[ \sum_{i=1}^{|\mathcal{V}|} c_{1i}^{\top} \alpha_i \dots \sum_{i=1}^{|\mathcal{V}|} c_{n|\mathcal{V}|i}^{\top} \alpha_i \right]^{\top} \in \mathbb{R}^{n|\mathcal{V}|},$$

and similarly for $R^{-1}\beta$. Therefore, the following holds:

$$\|F(t, R^{-1}\alpha) - F(t, R^{-1}\beta)\|_1 = \sum_{j=1}^{|\mathcal{V}|} \left\| f_j \left( t, \sum_{i=1}^{|\mathcal{V}|} c_{ji}^{\top} \alpha_i \right) - f_j \left( t, \sum_{i=1}^{|\mathcal{V}|} c_{ji}^{\top} \beta_i \right) \right\|_1 \\ \le c \sum_{j=1}^{|\mathcal{V}|} L_j(t) \|\alpha - \beta\|_1$$

where $L_j(t)$ is a Lipschitz constant of $f_i(t, \cdot)$, and $c = \max_{j=1,...,|\mathcal{V}|} \max_{i=1,...,|\mathcal{V}|} \|c_{ji}\|_1$. Consequently,

$$\|h(t, \alpha) - h(t, \beta)\|_1 \le L_h(t) \|\alpha - \beta\|_1,$$

where $L_0, L_h(t) := c \|R\|_1 \sum_{j=1}^{|\mathcal{V}|} L_j(t) + L_0 \|\hat{R}\|_1 \|R^{-1}\|_1$ are Lipschitz constants of $\lambda^{\text{im}}(\cdot)$ and $h(t, \cdot)$, respectively. Moreover, for all $t \in \mathbb{R}_+$ and $z \in \text{dom } h(t, \cdot)$:

$$\|h(t, z)\|_1 \le \|h(t, z) - h(t, z^*)\|_1 + \|h(t, z^*)\|_1 \le L_h(t)\|z\|_1 + (L_h(t)\|z^*\|_1 + \|h(t, z^*)\|_1) \\ = \gamma_1(t)\|z\|_1 + \gamma_2(t),$$

where $z^* \in \text{dom } h(t, \cdot)$ such that $\|z^*\|_1 < +\infty$ and $\|\lambda^{\text{im}}(R^{-1}z^*)\|_1 < +\infty$, $\gamma_1(\cdot) := L_h(\cdot)$, $\gamma_2(\cdot) := L_h(\cdot)\|z^*\|_1 + \|h(\cdot, z^*)\|_1$ are locally integrable on $\mathbb{R}_+$ since $L_h(\cdot)$ is uniformly bounded on $\mathbb{R}_+$ and $f_i$ is locally integrable in the first variable on $\mathbb{R}_+$ for all $i = \overline{1, |\mathcal{V}|}$. Therefore, (19) has a unique AC solution for any initial condition $z(0) = z_0 \in \text{Dom } \mathbf{\Psi} = R(N(\Theta^{\top} \otimes I_p)C)^{-1}(\text{Dom } \mathbf{\Phi})$, where $\text{Dom } \mathbf{\Phi} = (\text{SOL}(0, \Gamma))^* - r$, see Remark 3. This is equivalent to: (13) has a unique absolutely continuous solution for any initial condition $x(0) = x_0 \in (N(\Theta^{\top} \otimes I_p)C)^{-1}((\text{SOL}(0, \Gamma))^* - r)$. $\square$
**Remark 4.** It follows from (11), (10) and (9) that the following holds:

$$ 
\begin{aligned}
u(t) = (\Theta \otimes I_q)\eta(t) &= (\Theta \otimes I_q)H\lambda(t) \\
&\in (\Theta \otimes I_q)HSOL(Nv(t) + r, M) \\
&\equiv (\Theta \otimes I_q)HSOL(N(\Theta^\top \otimes I_p)(Cx(t) + Du(t)) + r, M),
\end{aligned}
$$

where the right-hand side may be a function of $u(t)$ in general. Therefore, the input $u(t)$ in (11) is implicitly represented. When inserted into the overall dynamics system in (13), its existence is guaranteed by Theorem 2. Indeed, the well-posedness of (13) is guaranteed by Theorem 2, which means that there exists $\lambda(t)$ satisfying (13b), and (13c). Hence, $u(t) = (\Theta \otimes I_q)H\lambda(t)$ is defined. In addition, such $\lambda(t)$ may not be unique, so neither can $u(t)$.

The main objective of this work is to find sufficient conditions to guarantee that all agents in the network, as defined in Section 1.2, reach a common state. This behavior of the network is referred to as synchronization, defined as follows.

**Definition 9.** (Global exponential synchronization). The network of agents (5) with the coupling (6) and the interconnection law (7) is globally exponentially synchronized if and only if for all admissible initial states $x(0) = [x_1(0)^\top \dots x_{|V|}(0)^\top]^\top$, there exist constants $\gamma > 0$, and $K(x(0)) > 0$ (depending on $x(0)$) such that:

$$ \|x_i(t) - x_j(t)\| \le K(x(0))e^{-\gamma t}, \quad \forall i, j = 1, \dots, |V|. $$

**Definition 10.** (Global finite-time synchronization). The network of agents (5) with the coupling (6) and the interconnection law (7) is globally finite-time synchronized if and only if for all admissible initial states $x(0) = [x_1(0)^\top \dots x_{|V|}(0)^\top]^\top$, there exists a finite $T(x(0)) > 0$ (depending on $x(0)$) such that $\|x_i(t) - x_j(t)\| = 0$, for all $t \ge T$, for all $i, j = 1, \dots, |V|$.

Let $\Theta_c$ be an oriented incidence matrix associated with the complete graph $K_{|V|} = (V, \mathcal{E}_c)$, where $|\mathcal{E}_c| = \frac{|V|(|V|-1)}{2}$, as in (20) below. Let us consider the error state $e(t)$, which is defined by $e(t) := (\Theta_c^\top \otimes I_n)x(t)$ with $x(\cdot)$ a solution of (13). Denote $e(t) = [e_{12}(t)^\top \dots e_{ij}(t)^\top \dots e_{|V|-1,|V|}(t)^\top]^\top$ with $1 \le i < j \le |V|$, then $e_{ij}(t) = x_i(t) - x_j(t)$.

$$ 
\Theta_c = \underbrace{\begin{bmatrix} 1 & 1 & \cdots & 0 & 0 & \cdots & 0 & 0 \\ -1 & 0 & \cdots & 1 & 1 & \cdots & 0 & 0 \\ 0 & -1 & \cdots & -1 & 0 & \cdots & 0 & 0 \\ 0 & 0 & \cdots & 0 & -1 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 0 & 0 & \cdots & 0 & 1 \\ 0 & 0 & \cdots & 0 & 0 & \cdots & -1 & -1 \end{bmatrix}}_{|\mathcal{E}_c| = \frac{|V|(|V|-1)}{2} \text{ columns}} \quad |V| \text{ rows, } e(t) = \underbrace{\begin{bmatrix} x_1(t) - x_2(t) \\ x_1(t) - x_3(t) \\ \vdots \\ x_2(t) - x_3(t) \\ x_2(t) - x_4(t) \\ \vdots \\ x_{|V|-1}(t) - x_{|V|}(t) \end{bmatrix}}_{|V| \text{ rows, } e(t) = \begin{bmatrix} x_1(t) - x_2(t) \\ x_1(t) - x_3(t) \\ \vdots \\ x_2(t) - x_3(t) \\ x_2(t) - x_4(t) \\ \vdots \\ x_{|V|-1}(t) - x_{|V|}(t) \end{bmatrix}}. \quad (20) $$
It follows from (13) that the dynamics of the error state is given by:

$$ \dot{e}(t) = (\Theta_c^\top \otimes I_n)F(t, x(t)) + (\Theta_c^\top \otimes I_n)B(\Theta \otimes I_q)H\lambda(t), \quad (21a) $$

$$ w(t) = N(\Theta_c^\top \otimes I_p)Cx(t) + \Gamma\lambda(t) + r, \quad (21b) $$

$$ 0 \le \lambda(t) \perp w(t) \ge 0, \quad (21c) $$

where $x(\cdot)$ is a solution of the TVCS (13). Intuitively, the network is globally exponentially (resp. finite-time) synchronized if the system (21) has an equilibrium point at the origin which is globally exponentially (resp. finite-time) stable. Indeed, assuming that the origin is the globally exponentially (resp. finite-time) stable equilibrium point of (21), there exist positive constants $k, \gamma$ (resp. $T(e(0)) < +\infty$) such that $\|e(t)\| \le k\|e(0)\|e^{-\gamma t}$ (resp. $e(t) = 0$, for all $t \ge T(e(0))$), for all initial error state $e(0)$, using [76, Definition 4.5] (resp. [77, Definition 8]). In addition, since each element of $e(t)$ has the form $x_i(t) - x_j(t)$, for some $i, j \in \{1, ..., |V|\}$, it follows that $\|x_i(t) - x_j(t)\| \le K(x(0))e^{-\gamma t}$ (resp. $\|x_i(t) - x_j(t)\| = 0$, for all $t \ge T(x(0))$), for all $i, j = 1, ..., |V|$, with $K(x(0)) = k\|(\Theta_c^\top \otimes I_n)x(0)\|$.

## 3 Analysis of the Error Dynamics (21) Equilibria

As mentioned in the previous section, the existence of an equilibrium point at the origin in the error dynamics (21) is a prerequisite for conducting further investigations into the synchronization of the network. In fact, proving the existence of an equilibrium point for a general system is a challenging problem. In most of the previous studies on network synchronization (e.g., [49], [?], [?]), the error dynamics are typically modeled by an ordinary differential equation (ODE) whose solution existence is ensured either by classical ODE theory (in the smooth case) or by Filippov's theory (in the nonsmooth case). In that context, the (unique) existence of an equilibrium point at the origin in the error dynamics has been implicitly guaranteed through assumptions that ensure:

1. Solutions $e(\cdot)$ to the error dynamics exist for all initial states $e(0)$ within the entire state space.

2. The existence of a Lyapunov function $V(e(\cdot))$ whose time derivative is negative for all $e(t) \ne 0$. (Lyapunov method)

In this work, the error dynamics is defined by the TVCS (21), which is a combination of the ODE (21a) and the complementarity constraint (21c). If the well-posedness of (21) is guaranteed, the Lyapunov method can be applied to prove that the origin is the equilibrium point of (21). This is established in Theorem 3. However, in a general analysis of the equilibria of (21), the Lyapunov method requires strong assumptions that may not be essential, or even impractical in cases where (21) does not have an equilibrium point at the origin, admits multiple equilibrium points, or when the origin is not globally exponentially stable. Therefore, in this section, rather than using the Lyapunov method, we directly investigate the existence of equilibrium points of (21)
(not limited to the origin), using three approaches: one is via the LCPs, another is with generalized equations, and a third is by using a dry friction condition. These approaches allow us to broaden the investigation of TVCS equilibria.

## 3.1 Equilibrium Points of (21) with an LCP Approach

Following this approach, the original definition of an equilibrium point of (21) is presented below, where SOL($q, Q$) denotes the set of all solutions to the LCP($q, Q$) as defined in Definition 7.

**Definition 11.** A point $e^* \in \mathbb{R}^{n|\mathcal{E}_c|}$ is called an equilibrium point of the system (21) if and only if for all $t \ge 0$, there exists $x^*(t) \in \mathbb{R}^{n|\mathcal{V}|}$ and a multiplier $\lambda^*(t) \in \text{SOL}(N(\Theta^\top \otimes I_p)Cx^*(t) + r, \Gamma)$ such that:

$$ (\Theta_c^\top \otimes I_n)F(t, x^*(t)) + (\Theta_c^\top \otimes I_n)B(\Theta \otimes I_q)H\lambda^*(t) = 0, \quad (22) $$

where $(x^*(t), \lambda^*(t))$ solves the TVCS in (13). Moreover, $e^* = (\Theta_c^\top \otimes I_n)x^*(t)$.

In Definition 11, it is noteworthy that $x^*$ and $\lambda^*$ may be time-varying, while $e^*$ is a constant. Denote by $\Lambda_i \in \mathbb{R}^{n \times m|\mathcal{E}|}$ the matrices such that $B(\Theta \otimes I_q)H = [\Lambda_1^\top \dots \Lambda_{|\mathcal{V}|}^\top]^\top$, for $i = 1, \dots, |\mathcal{V}|$. Definition 11 can be rewritten in the following equivalent forms:

$$ \begin{array}{l} (21) \text{ has an equilibrium point} \\ e^* = (\Theta_c^\top \otimes I_n)x^*(t). \end{array} \quad \Leftrightarrow \quad \left\{ \begin{array}{l} \exists (x^*(t), \lambda^*(t)) \text{ which solves (13):} \\ F(t, x^*(t)) + B(\Theta \otimes I_q)H\lambda^*(t) \in \ker(\Theta_c^\top \otimes I_n), \forall t \ge 0. \\ \\ \exists (x^*(t), \lambda^*(t)) \text{ which solves (13): for all } i, j = 1, \dots, |\mathcal{V}|, \\ f_i(t, x_i^*(t)) - f_j(t, x_j^*(t)) = (\Lambda_j - \Lambda_i)\lambda^*(t), \forall t \ge 0. \end{array} \right. \quad (23) $$

**Remark 5.** For an equilibrium point $e^* = (\Theta_c^\top \otimes I_n)x^*(t)$, using [78, Theorem 6], then

$$ x^*(t) = (\Theta_c^\top \otimes I_n)^\dagger e^* + v(t), \quad \text{for some } v(t) \in \ker(\Theta_c^\top \otimes I_n), $$

where $(\Theta_c^\top \otimes I_n)^\dagger = \frac{1}{|\mathcal{V}|}(\Theta_c \otimes I_n)$ by Lemma B1. Since $\ker(\Theta_c^\top \otimes I_n) = \{[\alpha(t)^\top \dots \alpha(t)^\top]^\top : \alpha(t) \in \mathbb{R}^n\}$, therefore

$$ x_i^*(t) = \alpha(t) + p_i, \quad \text{for some } \alpha: \mathbb{R}_+ \to \mathbb{R}^n, \text{ and } p_i \in \mathbb{R}^n, \text{ for all } i = 1, \dots, |\mathcal{V}|. $$

Without loss of generality, we can consider $p_i = 0$ for some $i = 1, \dots, |\mathcal{V}|$.

Let us recall that the problem discussed in this section is whether the system (21) admits any equilibrium point, and if so, what are they? To address this problem we need to find a solution $x^*(\cdot)$ of (13) such that (23) holds, or to show that no such $x^*(\cdot)$ exists. Note that the existence of a solution for (13) under Theorem 2 is not sufficient to indicate the existence of such an $x^*(\cdot)$. Let $\widehat{NC} = [(N_{ij}(C_i - C_j))^\top \dots (N_{kl}(C_k - C_l))^\top]^\top$ for $i < j$ and $k < l$ such that $(i, j), (k, l) \in \mathcal{E}$. The following proposition provides some necessary conditions for (21) to have an equilibrium point.
**Proposition 2.** Suppose that (21) has an equilibrium point $e^*$. The following statements hold:

1. Let $e^* = [e_{12}^{*\top} \dots e_{ij}^{*\top} \dots e_{|V|-1,|V|}^{*\top}]^\top$ with $1 \le i < j \le |V|$. Then
$$e_{ij}^* = e_{i|V|}^* - e_{j|V|}^*, \quad \forall 1 \le i < j < |V|$$

2. There exist a function $\alpha: \mathbb{R}_+ \to \mathbb{R}^n$ and a tuple $(p_1, p_2, \dots, p_{|V|})$ where $p_i = 0$ for some $i$:

2a. $f_i(t, \alpha(t) + p_i) - f_j(t, \alpha(t) + p_j) \in \text{rge}(\Lambda_i - \Lambda_j)$, $\forall t \ge 0$, $\forall i, j = 1, \dots, |V|$.

2b. Let $p = [p_1^\top \dots p_{|V|}^\top]^\top$. If $\Gamma \succcurlyeq 0$ (not necessarily symmetric), then for all $t \ge 0$,

$$\widehat{NC}\alpha(t) + N(\Theta^\top \otimes I_p)Cp + r \in \mathbb{R}_+^{m|\mathcal{E}|} + \text{rge}(\Gamma + \Gamma^\top) + \{\Gamma^\top \gamma \mid \gamma \in \mathbb{R}_+^{m|\mathcal{E}|}\}.$$

*Proof.* [1.] Since $e^*$ is the equilibrium point of (21), it can be deduced that there exists $x^*(t) = [x_1^*(t)^\top \dots x_{|V|}^*(t)^\top]^\top$ such that $e^* = (\Theta_c^\top \otimes I_n)x^*(t)$. For all $1 \le i < j < |V|$, using the definition of $\Theta_c$ in (20), then $e_{ij}^* = x_i^*(t) - x_j^*(t)$. Thus,

$$e_{ij}^* = x_i^*(t) - x_{|V|}^*(t) - x_j^*(t) + x_{|V|}^*(t) = e_{i|V|}^* - e_{j|V|}^*.$$

[2.] Using the equivalence (23), there exists $x^*(t)$ which solves (13) such that:

$$e^* = (\Theta_c^\top \otimes I_n)x^*(t) \qquad (24a)$$

$$\exists \lambda^*(t) \in \text{SOL}(N(\Theta^\top \otimes I_p)Cx^*(t) + r, \Gamma): \qquad (24b)$$

$$f_i(t, x_i^*(t)) - f_j(t, x_j^*(t)) = (\Lambda_j - \Lambda_i)\lambda^*(t), \quad \forall t \ge 0, \forall i, j = 1, \dots, |V|. \qquad (24c)$$

Using (24a) and Remark 5, there exist $\alpha(t)$ and a tuple $(p_1, \dots, p_{|V|})$ where $p_i = 0$ for some $i$ such that $x_i^*(t) = \alpha(t) + p_i$ for all $i = 1, \dots, |V|$. Therefore, (24c) implies that item 2a holds. In addition, it follows from (24b) that

$$\text{SOL}(\widehat{NC}\alpha(t) + N(\Theta^\top \otimes I_p)Cp + r, \Gamma) \neq \emptyset$$

since $N(\Theta^\top \otimes I_p)Cx^*(t) = \widehat{NC}\alpha(t) + N(\Theta^\top \otimes I_p)Cp$. This is equivalent to $\widehat{NC}\alpha(t) + N(\Theta^\top \otimes I_p)Cp + r \in (\text{SOL}(0, \Gamma))^*$ if $\Gamma \succcurlyeq 0$ [56, Corollary 3.8.10, Exercise 3.12.1]. Using Remark 3, the proof is complete. $\square$

Sufficient conditions to ensure that (21) has an equilibrium point are given below. Moreover these results provide a method for calculating some equilibrium points of (21).

**Proposition 3.** Suppose that there exists a tuple $(p_1, \dots, p_{|V|})$ such that $f_i(t, p_i) = 0$ for all $i = 1, \dots, |V|$, and $0 \in \text{SOL}(N(\Theta \otimes I_p)Cp + r, \Gamma)$ where $p = [p_1^\top \dots p_{|V|}^\top]^\top$. Then $(\Theta_c^\top \otimes I_n)p$ is an equilibrium point of (21).
Proof. Let $x^*(t) = p$ and $\lambda^*(t) = 0$, then $(x^*(t), \lambda^*(t))$ solves (13). In addition, $f_i(t, p_i) - f_j(t, p_j) = 0 = (\Lambda_j - \Lambda_i)\lambda^*(t)$ for all $i, j = 1, ..., |V|$. Using (23), then $(\Theta_c^\top \otimes I_n)x^*(t) = (\Theta_c^\top \otimes I_n)p$ is an equilibrium point of (21). $\square$

The following corollary is a direct consequence of Proposition 3.

**Corollary 1.** Suppose that there exists $p \in \mathbb{R}^n$ such that $f_i(t, p) = 0$ for all $i = 1, ..., |V|$, and $0 \in \text{SOL}(\widehat{NC}p + r, \Gamma)$. Then the origin is an equilibrium point of (21).

An alternative sufficient condition is provided in the following proposition, which applies to cases where the conditions in Corollary 1 may not be satisfied, but still ensures that the origin is an equilibrium point of (21).

**Proposition 4.** Suppose that there exist $\alpha(t) \in \mathbb{R}^n$ and $\lambda^*(t) \in \text{SOL}(\widehat{NC}\alpha(t) + r, \Gamma)$ such that:

$$f_i(t, \alpha(t)) + \Lambda_i \lambda^*(t) = \dot{\alpha}(t), \quad \forall t \ge 0, \quad \forall i = 1, ..., |V|. \quad (25)$$

Then the origin is an equilibrium point of (21).

Proof. Let $x^*(t) = [x_1^*(t)^\top \ x_2^*(t)^\top \ ... \ x_{|V|}^*(t)^\top]^\top$ where $x_i^*(t) = \alpha(t)$ for all $i = 1, ..., |V|$. It follows from (25) that:

$$\begin{cases} (x^*(t), \lambda^*(t)) \text{ solves } (13), \text{ since } \widehat{NC}\alpha(t) = N(\Theta^\top \otimes I_p)Cx^*(t). \\ f_i(t, x_i^*(t)) - f_j(t, x_j^*(t)) = (\Lambda_j - \Lambda_i)\lambda^*(t), \quad \forall t \ge 0, \quad \forall i, j = 1, ..., |V|. \end{cases}$$

Using (23), the origin is an equilibrium point of (21). $\square$

The following corollary is a consequence of Proposition 4.

**Corollary 2.** Suppose that $f_i \equiv f$ for all $i = 1, ..., |V|$, and the ODE: $\dot{x}(t) = f(t, x(t))$ has a solution. Assume further that $N_{ij}C_i = N_{ij}C_j$ (see (5) and (6b)) and $r_{ij} \ge 0$, for all $i < j$ such that $(i, j) \in \mathcal{E}$. Then the origin is an equilibrium point of (21).

Proof. Let $\alpha(t)$ be a solution of the ODE: $\dot{x}(t) = f(t, x(t))$, and $\lambda^*(t) \equiv 0$. It can be checked that Proposition 4 is satisfied since $\widehat{NC} = 0$ and $0 \in \text{SOL}(r, \Gamma)$ since $r \ge 0$. $\square$

## 3.2 Equilibrium Points as Generalized Equations Solutions

Noting that (21b) and (21c) are equivalent to (14), the error dynamics (21) is rewritten as the DI:

$$\dot{e}(t) \in (\Theta_c^\top \otimes I_n)F(t, x(t)) - (\Theta_c^\top \otimes I_n)B(\Theta \otimes I_q)H\Phi(N(\Theta^\top \otimes I_p)Cx(t)), \quad (26)$$

where $\Phi(\cdot)$ is defined in (15). Following this approach, a point $e^* = (\Theta_c^\top \otimes I_n)x^*(t) \in \mathbb{R}^{|\mathcal{E}|}$ is an equilibrium point of (21) if and only if there exists $x^*(t)$ which solves (13) and such that for all $t \ge 0$,

$$0 \in (\Theta_c^\top \otimes I_n)F(t, x^*(t)) - (\Theta_c^\top \otimes I_n)B(\Theta \otimes I_q)H\Phi(N(\Theta^\top \otimes I_p)Cx^*(t)). \quad (27)$$
The prerequisite for the existence of an equilibrium point of (21) is that (27) has a solution, which is analyzed in this section. The following lemma provides a special structure of the coupling that ensures (27) has a solution. Moreover, this structure will be used consistently throughout the remainder of this work, as it plays a critical role in proving the existence of an equilibrium point at the origin for error dynamics (21) under the dry friction condition (Proposition 13), as well as in establishing the global exponential (finite-time) synchronization of the network (Theorem 3, Corollary 6).

**Lemma 1.** Consider the agents in (5) and the couplings in (6). Suppose that $q = p$ and for each pair $(i, j) \in \mathcal{V} \times \mathcal{V}$ where $i < j$ and $(i, j) \in \mathcal{E}$, there exist matrices $W_{ij} \in \mathbb{R}^{p \times p}$, $H'_{ij} \in \mathbb{R}^{p \times m'}$, $M'_{ij} \in \mathbb{R}^{m' \times m'}$, $N'_{ij} \in \mathbb{R}^{m' \times p}$, $r'_{ij} \in \mathbb{R}^{m'}$ where $m' = m - 2p$ such that:

$$ H_{ij} = [H'_{ij} | W_{ij}(I_p \otimes [1 - 1])] , \quad M_{ij} = \begin{bmatrix} M'_{ij} & 0 \\ 0 & I_{2p} \end{bmatrix} , \quad N_{ij} = \begin{bmatrix} N'_{ij} \\ I_p \otimes \begin{bmatrix} 1 \\ -1 \end{bmatrix} \end{bmatrix} , \quad r_{ij} = \begin{bmatrix} r'_{ij} \\ \mathbf{0}_{2p} \end{bmatrix} . \qquad (28) $$

Then, the interconnection law in (7) is rewritten as:

$$ \left\{ \begin{aligned} u_i(t) &= \sum_{j \in \mathcal{N}_i} u_{ij}(t), \\ u_{ij}(t) &= -u_{ji}(t) = \eta'_{ij}(t) - W_{ij}v_{ij}(t), \\ v_{ij}(t) &= y_i(t) - y_j(t), \end{aligned} \right. \qquad (29) $$

where $\eta'_{ij}(t)$ is defined similarly to (6), using the coefficient matrices $H'_{ij}$, $M'_{ij}$, $N'_{ij}$, $r'_{ij}$.

*Proof.* Recall that for each pair $(i, j) \in \mathcal{V} \times \mathcal{V}$ where $i < j$ and $(i, j) \in \mathcal{E}$, the coupling defined in (6) has the form:

$$ \left\{ \begin{aligned} \eta_{ij}(t) &= H_{ij}\lambda_{ij}(t), \\ w_{ij}(t) &= M_{ij}\lambda_{ij}(t) + N_{ij}v_{ij}(t) + r_{ij}, \\ 0 &\le \lambda_{ij}(t) \perp w_{ij}(t) \ge 0. \end{aligned} \right. \qquad (30) $$

Under condition (28), the system (30) is rewritten as:

$$ \left\{ \begin{aligned} \eta_{ij}(t) &= H'_{ij}\lambda'_{ij}(t) + W_{ij}(I_p \otimes [1-1])\lambda''_{ij}(t), \\ w'_{ij}(t) &= M'_{ij}\lambda'_{ij}(t) + N'_{ij}v_{ij}(t) + r'_{ij}, \\ w''_{ij}(t) &= \lambda''_{ij}(t) + (I_p \otimes [1-1]^\top)v_{ij}(t), \\ 0 &\le \lambda'_{ij}(t) \perp w'_{ij}(t) \ge 0, \\ 0 &\le \lambda''_{ij}(t) \perp w''_{ij}(t) \ge 0, \end{aligned} \right. \qquad (31) $$

where $\lambda'_{ij}(t), w'_{ij}(t) \in \mathbb{R}^{m'}$ and $\lambda''_{ij}(t), w''_{ij}(t) \in \mathbb{R}^{2p}$ such that $\lambda_{ij}(t) = [\lambda'_{ij}(t)^\top \ \lambda''_{ij}(t)^\top]^\top$, and $w_{ij}(t) = [w'_{ij}(t)^\top \ w''_{ij}(t)^\top]^\top$. Set $\eta'_{ij}(t) = H'_{ij}\lambda'_{ij}(t)$ and $\eta''_{ij}(t) = (I_p \otimes [-1 \ 1])\lambda''_{ij}(t)$,
then $\eta_{ij}(t) = \eta'_{ij}(t) - W_{ij}\eta''_{ij}(t)$. In addition, the set of equations (31) can be split into two subsets as follows:

$$ egin{cases} 
\eta'_{ij}(t) = H'_{ij}\lambda'_{ij}(t), \\
w'_{ij}(t) = M'_{ij}\lambda'_{ij}(t) + N'_{ij}v_{ij}(t) + r'_{ij}, \\
0 \le \lambda'_{ij}(t) \perp w'_{ij}(t) \ge 0,
\end{cases} 
\quad ; 
\quad 
\begin{cases} 
\eta''_{ij}(t) = (I_p \otimes [-1 \ 1])\lambda''_{ij}(t), \\
w''_{ij}(t) = \lambda''_{ij}(t) + (I_p \otimes [1 \ -1]^\top)v_{ij}(t), \\
0 \le \lambda''_{ij}(t) \perp w''_{ij}(t) \ge 0.
\end{cases} 
\tag{32} $$

It can be checked that $\eta''_{ij}(t) = v_{ij}(t)$. Indeed, let $\eta''_{ij}(t) = [\eta''_{ij,1}(t) \dots \eta''_{ij,p}(t)]^\top \in \mathbb{R}^p$, $v_{ij}(t) = [v_{ij,1}(t) \dots v_{ij,p}(t)]^\top \in \mathbb{R}^p$, $\lambda''_{ij}(t) = [\lambda''_{ij,1}(t)^\top \dots \lambda''_{ij,p}(t)^\top]^\top \in \mathbb{R}^{2p}$, $w''_{ij}(t) = [w''_{ij,1}(t)^\top \dots w''_{ij,p}(t)^\top]^\top \in \mathbb{R}^{2p}$. It follows from (32) that for all $k = 1, \dots, p$:

$$ egin{cases} 
\eta''_{ij,k}(t) = [-1 \ 1]\lambda''_{ij,k}(t), \\
w''_{ij,k}(t) = \lambda''_{ij,k}(t) + [1 \ -1]^\top v_{ij,k}(t), \\
0 \le \lambda''_{ij,k}(t) \perp w''_{ij,k}(t) \ge 0.
\end{cases} $$

Solving these LCP, we obtain that $\eta''_{ij,k}(t) = v_{ij,k}(t)$ for all $k = 1, \dots, p$. Therefore $\eta''_{ij}(t) = v_{ij}(t)$, and thus $\eta_{ij}(t) = \eta'_{ij}(t) - W_{ij}v_{ij}(t)$. It follows that the input $u_i$ in (7) can be reformulated as (29). The proof is complete. $\square$

In the remainder of this section, we assume that (28) holds. Using this specially structured coupling, it can be proved that the inclusion (27) has a solution $x^*(t)$ under certain sufficient conditions. To clarify this, let us follow the steps outlined below. Firstly, using Lemma 1, proceeding similarly to (11), it follows from (29) that:

$$ egin{cases} 
u(t) = (\Theta \otimes I_p)\eta'(t) - (\Theta \otimes I_p)Wv(t) = (\Theta \otimes I_p)\eta'(t) - (\Theta \otimes I_p)W(\Theta^\top \otimes I_p)y(t), \\ 
v(t) = (\Theta^\top \otimes I_p)y(t), 
\end{cases} 
\tag{33} $$

in which, $\eta'(t)$ is defined as follows:

$$ egin{cases} 
\eta'(t) = H'\lambda'(t), \\
w'(t) = M'\lambda'(t) + N'v(t) + r', \\
0 \le \lambda'(t) \perp w'(t) \ge 0,
\end{cases} 
\tag{34} $$

where $\eta'(t), \lambda'(t), w'(t)$ and $W \in \mathbb{R}^{p|\mathcal{E}| \times p|\mathcal{E}|}$, $H' \in \mathbb{R}^{p|\mathcal{E}| \times m'|\mathcal{E}|}$, $M' \in \mathbb{R}^{m'|\mathcal{E}| \times m'|\mathcal{E}|}$, $N' \in \mathbb{R}^{m'|\mathcal{E}| \times p|\mathcal{E}|}$, $r' \in \mathbb{R}^{m'|\mathcal{E}|}$ are defined similarly to $\eta(t), \lambda(t), w(t)$ and $H, M, N, r$ in Section 1.3, respectively. Using (9), the following is obtained:

$$ u(t) = (\Theta \otimes I_p)\eta'(t) - (\Theta \otimes I_p)W(\Theta^\top \otimes I_p)(Cx(t) + Du(t)). \tag{35} $$

Let $\Xi := I_{p|\mathcal{E}|} + (\Theta \otimes I_p)W(\Theta^\top \otimes I_p)D$, and $\Delta := I_{p|\mathcal{E}|} + (\Theta^\top \otimes I_p)D(\Theta \otimes I_p)W$. The following assumption is considered:

**Assumption 1.** The matrix $\Delta$ is invertible.
Under this assumption, $\Xi^{-1}$ exists [73, Corollary 2.8.8]. Therefore, the input $u(t)$ in (35) is rewritten as follows:

$$u(t) = \Xi^{-1}(\Theta \otimes I_p)\eta'(t) - \Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)Cx(t). \quad (36)$$

By combining (9), (36), (34), (33) the overall dynamics of the entire network is given as the TVCS:

$$\left\{
\begin{array}{l}
\dot{x}(t) = F(t, x(t)) - B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)Cx(t) + B\Xi^{-1}(\Theta \otimes I_p)H'\lambda'(t), \\
w'(t) = N'\Delta^{-1}(\Theta^{\top} \otimes I_p)Cx(t) + \Gamma'\lambda'(t) + r', \\
0 \le \lambda'(t) \perp w'(t) \ge 0,
\end{array}
\right.
\quad (37)$$

where $\Gamma' := N'(\Theta^{\top} \otimes I_p)D\Xi^{-1}(\Theta \otimes I_p)H' + M'$, and $\Delta^{-1} = I_{p|\mathcal{E}|} - (\Theta^{\top} \otimes I_p)D\Xi^{-1}(\Theta \otimes I_p)W$, using [73, Corollary 2.8.8]. In addition, (37) is equivalent to (13).

**Remark 6.** When $W$ and $D$ are symmetric matrices, then $\Xi^{-1}(\Theta \otimes I_p) = (\Theta \otimes I_p)(\Delta^{-1})^{\top}$, and the matrices $D\Xi^{-1}$, $(\Delta^{-1})^{\top}W$ are symmetric. Indeed, since $\Xi(\Theta \otimes I_p) = (\Theta \otimes I_p)(I_{p|\mathcal{E}|} + W(\Theta^{\top} \otimes I_p)D(\Theta \otimes I_p)) = (\Theta \otimes I_p)\Delta^{\top}$, the following holds:

$$\Xi^{-1}(\Theta \otimes I_p) = \Xi^{-1}(\Theta \otimes I_p)\Delta^{\top}(\Delta^{-1})^{\top} = \Xi^{-1}\Xi(\Theta \otimes I_p)(\Delta^{-1})^{\top} = (\Theta \otimes I_p)(\Delta^{-1})^{\top}.$$

In addition, using [73, Corollary 2.8.8], the following hold:

$$\begin{aligned}
(D\Xi^{-1})^{\top} &= (\Xi^{\top})^{-1}D = (I_{p|\mathcal{V}|} + D(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p))^{-1}D \\
&= (I_{p|\mathcal{V}|} - D(I_{p|\mathcal{V}|} + (\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)D)^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p))D \\
&= D - D\Xi^{-1}(\Xi - I_{p|\mathcal{V}|}) = D\Xi^{-1}.
\end{aligned}$$

$$\begin{aligned}
((\Delta^{-1})^{\top} W)^{\top} &= W \Delta^{-1} = W(I_{p|\mathcal{E}|} + (\Theta^{\top} \otimes I_p)D(\Theta \otimes I_p)W)^{-1} \\
&= W(I_{p|\mathcal{E}|} - (\Theta^{\top} \otimes I_p)D(\Theta \otimes I_p)(I_{p|\mathcal{E}|} + W(\Theta^{\top} \otimes I_p)D(\Theta \otimes I_p))^{-1}W) \\
&= W - W(\Theta^{\top} \otimes I_p)D(\Theta \otimes I_p)(\Delta^{\top})^{-1}W = W - (\Delta^{\top} - I_{p|\mathcal{E}|})(\Delta^{\top})^{-1}W \\
&= (\Delta^{-1})^{\top}W.
\end{aligned}$$

Let us introduce Assumptions 2-4, which play an important role in the main results of this section.

**Assumption 2.** $\Gamma' \ge 0$.

**Assumption 3.** There exists $P' \in \mathbb{R}^{n|\mathcal{V}| \times n|\mathcal{V}|}$ such that $P' = P'^{\top} > 0$, and

$$\ker (\Gamma' + (\Gamma')^{\top}) \subseteq \ker (P' B \Xi^{-1} (\Theta \otimes I_p) H' - C^{\top} (\Theta \otimes I_p) (\Delta^{-1})^{\top} N'^{\top}). \quad (38)$$

**Assumption 4.** rge $(N'\Delta^{-1}(\Theta^{\top} \otimes I_p)C) \cap \text{rint}((\text{SOL}(0, \Gamma'))^{\star} - r') \neq \emptyset$.

The well-posedness of (37) is established in the following corollary as a trivial adaptation of Theorem 2, which was originally proved for (13), to the closely related (37).
**Corollary 3.** Suppose that for each $i = 1, \dots, |V|$, $f_i: \mathbb{R}_+ \times \mathbb{R}^n \to \mathbb{R}^n$ is locally integrable on $\mathbb{R}_+$ in the first variable and Lipschitz continuous in the second variable, and Assumptions 1 through 4 hold. Then, (37) has a unique absolutely continuous solution $x(\cdot)$ for any initial condition $x(0) = x_0 \in (N'\Delta^{-1}(\Theta^\top \otimes I_p)C)^{-1}((\text{SOL}(0, \Gamma'))^* - r')$.

Let us now introduce the operator $\Phi'$, similar to $\Phi$ in (15), defined as follows, where $\Gamma'$ is given after (37):

$$ \Phi': \mathbb{R}^{m'|\mathcal{E}|} \to \mathbb{R}^{m'|\mathcal{E}|}, \quad v \mapsto \left( \mathcal{N}_{\mathbb{R}_+^{m'|\mathcal{E}|-r'}}^{-1} + \Gamma' \right)^{-1}(v). \quad (39) $$

**Remark 7.** In a similar manner to Remark 3, when Assumption 2 holds (i.e., $\Gamma' \succeq 0$), $\Phi'$ is a maximal monotone operator with domain $\text{Dom } \Phi' = (\text{SOL}(0, \Gamma'))^* - r'$, where

$$ (\text{SOL}(0, \Gamma'))^* = \mathbb{R}_+^{m'|\mathcal{E}|} + \text{rge}(\Gamma' + \Gamma'^{\top}) + \{\Gamma'^{\top}\gamma \mid \gamma \in \mathbb{R}_+^{m'|\mathcal{E}|}\}. \quad (40) $$

If $\Gamma' = \Gamma'^{\top}$, then this reduces to $\text{Dom } \Phi' = \mathbb{R}_+^{m'|\mathcal{E}|} + \text{rge}(\Gamma') - r'$.

Proceeding similarly to Section 1.3, (37) is rewritten as:

$$ \begin{aligned} \dot{x}(t) \in F(t, x(t)) - B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)Cx(t) \\ - B\Xi^{-1}(\Theta \otimes I_p)H'\Phi'(N'\Delta^{-1}(\Theta^{\top} \otimes I_p)Cx(t)), \end{aligned} \quad (41) $$

which is equivalent to (16). Therefore, the error dynamics (26) (or equivalently, (21)) is rewritten as:

$$ \dot{e}(t) \in (\Theta_c^\top \otimes I_n) [F(t, x(t)) - B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)Cx(t) \\ - B\Xi^{-1}(\Theta \otimes I_p)H'\Phi'(N'\Delta^{-1}(\Theta^\top \otimes I_p)Cx(t))] \quad (42) $$

Under Assumption 1 and special structure (28), the inclusion (27) can be rewritten in the equivalent form:

$$ 0 \in (\Theta_c^\top \otimes I_n) [F(t, x^*(t)) - B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)Cx^*(t) \\ - B\Xi^{-1}(\Theta \otimes I_p)H'\Phi'(N'\Delta^{-1}(\Theta^\top \otimes I_p)Cx^*(t))] \quad (43) $$

Let us introduce $\lambda^{\text{im}}(\cdot) := \text{proj}_{\text{rge}(\Gamma'+\Gamma'^{\top})}(\Phi'(N'\Delta^{-1}(\Theta^\top \otimes I_p)C))$, and the following function:

$$ \tilde{F}(t, \cdot) := -P'F(t, \cdot) + P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)C(\cdot) + (P'B\Xi^{-1}(\Theta \otimes I_p)H' \\ - C^\top(\Theta \otimes I_p)(\Delta^{-1})^\top(N')^\top)\lambda^{\text{im}}(\cdot), $$

where $P'$ is given in Assumption 3. The following proposition holds.

**Proposition 5.** Consider (43) and its equivalent form (27). Suppose that Assumptions 1-4 hold, and one of the following holds (A or B or C) for all $t \ge 0$:

A. $\tilde{F}(t, \cdot)$ is strongly maximal monotone.

B. $\tilde{F}(t, \cdot)$ is uniformly maximal monotone with a supercoercive modulus.
C. There exists a convex lower semicontinuous function $g(t, \cdot) : \mathbb{R}^{n|\nu|} \to \mathbb{R}$ such that $\tilde{F}(t, \cdot) = \partial g(t, \cdot)$, and $\text{rge}(\tilde{F}(t, \cdot)) = \mathbb{R}^{n|\nu|}$.

Then (43) has a solution, thus so does (27).

Proof. Since Assumption 1 holds, $\Xi^{-1}$ exists. For each $t \ge 0$, let us consider the following inclusion:

$$0 \in F(t, x) - B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)Cx \\ - B\Xi^{-1}(\Theta \otimes I_p)H'\Phi'(N'\Delta^{-1}(\Theta^{\top} \otimes I_p)Cx). \quad (44)$$

Proceeding similarly to the proof of Theorem 2, under Assumptions 2, 3, 4, this inclusion is rewritten as follows:

$$0 \in \tilde{F}(t, x) + \Psi(x), \quad (45)$$

where $\Psi(\cdot) = C^{\top}(\Theta \otimes I_p)(\Delta^{-1})^{\top}(N')^{\top}\Phi'(N'\Delta^{-1}(\Theta^{\top} \otimes I_p)C\cdot)$ is a maximal monotone operator, and

$$\tilde{F}(t, \cdot) = -P'F(t, \cdot) + P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)C(\cdot) + (P'B\Xi^{-1}(\Theta \otimes I_p)H' \\ - C^{\top}(\Theta \otimes I_p)(\Delta^{-1})^{\top}(N')^{\top})\lambda^{\text{im}'}(\cdot),$$

with $\lambda^{\text{im}'}(\cdot) := \text{proj}_{\text{rge}(\Gamma'+\Gamma'')}(\Phi'(N'\Delta^{-1}(\Theta^{\top} \otimes I_p)C\cdot))$ is Lipschitz continuous [72, Lemma 3]. Note that $\text{rint}(\text{Dom } \Psi) \neq \emptyset$. Indeed, since Assumption 4 holds, there exists $z \in \mathbb{R}^{n|\nu|}$ such that:

$$N'\Delta^{-1}(\Theta^{\top} \otimes I_p)Cz \in \text{rint}((\text{SOL}(0, \Gamma'))^{\star} - r').$$

Since Assumption 2 holds, $\text{Dom } \Phi' = (\text{SOL}(0, \Gamma'))^{\star} - r'$, using Remark 7. Therefore,

$$z \in (N'\Delta^{-1}(\Theta^{\top} \otimes I_p)C)^{-1} (\text{rint}(\text{Dom } \Phi')).$$

It can be checked that $\text{Dom } \Phi'$ is a convex set with note that $(\text{SOL}(0, \Gamma'))^{\star}$ is given in (40). Using [54, Proposition 2.44, (b)], the following holds:

$$z \in \text{rint} \left( (N'\Delta^{-1}(\Theta^{\top} \otimes I_p)C)^{-1} (\text{Dom } \Phi') \right) = \text{rint}(\text{Dom } \Psi).$$

Thus $\text{rint}(\text{Dom } \Psi) \neq \emptyset$. In addition, since $\text{dom } \tilde{F}(t, \cdot) = \text{Dom } \Psi$, the following holds:

$$\text{rint}(\text{dom } \tilde{F}(t, \cdot)) \cap \text{rint}(\text{Dom } \Psi) \neq \emptyset. \quad (46)$$

[If A holds (resp. B holds)]: Since (46) holds, $\tilde{F}(t, \cdot) + \Psi(\cdot)$ is a maximal monotone operator [54, Corollary 12.44]. Using Propositions B2-(i), B3-(i) (resp. B2-(ii), B3-(ii)), $\tilde{F}(t, \cdot)$ is 3$^{\star}$ monotone and surjective. Therefore, $\tilde{F}(t, \cdot) + \Psi(\cdot)$ is surjective, using Proposition B4-(ii). This implies that (45) has a unique solution $x^{\star}(t)$. Since (45) is equivalent to (44), $x^{\star}(t)$ is the unique solution of (44), which yields $x^{\star}(t)$ is a solution of (43). Since (43) is equivalent to (27), $x^{\star}(t)$ is a solution of (27).

[If C holds]: Using [57, Theorem 20.25], $\tilde{F}(t, \cdot)$ is a maximal monotone operator. It follows from (46) and [54, Corollary 12.44] that $\tilde{F}(t, \cdot) + \Psi(\cdot)$ is maximal monotone. In addition, $\tilde{F}(t, \cdot)$ is surjective since $\text{rge}(\tilde{F}(t, \cdot)) = \mathbb{R}^{n|\nu|}$, and is 3$^{\star}$ monotone, using Proposition B1. Therefore, $\tilde{F}(t, \cdot) + \Psi(\cdot)$ is surjective, using Proposition B4-(ii). Mimicking the above cases, (27) has a solution. $\square$
Let us introduce the following assumption, which is a special case of Assumption 3.

**Assumption 3'.** There exists $P' \in \mathbb{R}^{n|\nu| \times n|\nu|}$ such that $P' = P'^{\top} > 0$, and

$$P'B\Xi^{-1}(\Theta \otimes I_p)H' = C^{\top}(\Theta \otimes I_p)(\Delta^{-1})^{\top}N'^{\top}. \quad (47)$$

**Remark 8.** When Assumption 3' holds, the function $\tilde{F}(t, \cdot)$ in Proposition 5 is rewritten as:

$$\tilde{F}(t, \cdot) = -P'F(t, \cdot) + P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)C(\cdot).$$

In this case, item **A** (resp. item **B**) in Proposition 5 holds if for all $t \ge 0$, $\text{dom } F(t, \cdot) = \mathbb{R}^{n|\nu|}$, $F(t, \cdot)$ is continuous, and $\tilde{F}(t, \cdot)$ is strongly monotone (resp. uniformly monotone with a supercoercive modulus), using [54, Example 12.7]. In addition, it also follows from [79, Corollary 3.4] that if $F(t, \cdot)$ is Lipschitz continuous, and $\tilde{F}(t, \cdot)$ is strongly monotone, then item **A** in Proposition 5 holds.

The following corollary provides sufficient conditions to ensure that the strong monotonicity of the function $\tilde{F}(t, \cdot)$ in Proposition 5 is verified when Assumption 3' holds.

**Corollary 4.** Let $P'$ be given in Assumption 3'. Suppose that one of the following items holds:

1. For each $t \ge 0$, there exists a convex function $g(t, \cdot) : \mathbb{R}^{n|\nu|} \to \mathbb{R}$ such that $-P'F(t, \cdot) = \partial g(t, \cdot)$, and

$$(P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)C) + (P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)C)^{\top} > 0.$$

2. There exist a symmetric matrix $Q \in \mathbb{R}^{n|\nu| \times n|\nu|}$ such that $F$ is QUAD($P'$, $Q$) (see Definition 6), and the following holds for some $\mu > 0$:

$$(P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)C - Q) + (P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)C - Q)^{\top} \\ \geq \mu I_{n|\nu|}.$$

Then the function $\tilde{F}(t, \cdot)$ in Proposition 5 is strongly monotone for each $t \ge 0$.

**Proof.** **[If item 1 holds]**: For each $t \ge 0$, since $\partial g(t, \cdot)$ is a monotone operator [57, Example 20.3], $-P'F(t, \cdot)$ is monotone. The matrix inequality in item 1 implies that $P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)C(\cdot)$ is strongly monotone. Therefore, $-P'F(t, \cdot) + P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)C$ is strongly monotone. The conclusion follows from Remark 8.

**[If item 2 holds]**: Let any $x, y \in \mathbb{R}^{n|\nu|}$, since $F$ is QUAD($P'$, $Q$), we have:

$$\langle -P'(F(t, x) - F(t, y)), x - y \rangle \ge (x - y)^{\top}(-Q)(x - y).$$

It follows from the matrix inequality in item 2 that for all $t \ge 0$, the following holds:

$$\langle -P'(F(t, x) - F(t, y)) + P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)C(x - y), x - y \rangle \ge \frac{1}{2}\mu\|x - y\|^2.$$

Therefore, $-P'F(t, \cdot) + P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)C(\cdot)$ is strongly monotone. The conclusion follows from Remark 8. □
Recall that $\lambda^{\text{im}'}(\cdot) := \text{proj}_{\text{rge}(\Gamma'+\Gamma'^{\top})}(\mathbf{\Phi}'(N'\Delta^{-1}(\Theta^{\top}\otimes I_p)C\cdot))$. The following proposition will be useful when conditions stated in Proposition 5 do not apply.

**Proposition 6.** Suppose that for each $t \ge 0$, for all $i = \overline{1,|V|}$, $f_i(t, \cdot)$ is Lipschitz continuous, dom $f_i(t, \cdot) = \mathbb{R}^n$, Assumptions 1 through 4 holds, and for each pair $(i, j)$ satisfying $i < j$ and $(i, j) \in \mathcal{E}$, the following conditions hold:

1. $B_i W_{ij} = B_j W_{ij}$,

2. $B_i W_{ij} B_i^\top = B_i W_{ij} B_j^\top$,

3. $P'B(\Theta \otimes I_p) = C^\top(\Theta \otimes I_p)$, where $P'$ is given in Assumption 3,

4. $\lambda_{\min}(\widetilde{W} + \widetilde{W}^\top)\lambda_2(L)\lambda_{\min}^2(P') > 2(L_F\|P'\| + \|K_0\| + L_\lambda\|K_1\|)$,

where $\widetilde{W} := \text{diag}(B_i W_{ij} B_i^\top)$; $\lambda_2(L)$ is the second smallest eigenvalue of the Laplacian matrix $L = \Theta\Theta^\top$; $L_F := \max_{t \in \mathbb{R}_+} L_F(t)$ with $L_F(t)$ a Lipschitz constant of $F(t, \cdot)$; $L_\lambda$ is a Lipschitz constant of $\lambda^{\text{im}'}(\cdot)$; $K_0 := P'B(\Theta \otimes I_p)W\Delta^{-1}(\Theta^\top \otimes I_p)D(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)C$; $K_1 := P'B\Xi^{-1}(\Theta \otimes I_p)H' - C^\top(\Theta \otimes I_p)(\Delta^{-1})^\top N'^\top$. Then (43) has a solution $x^*(\cdot)$, thus so does (27).

*Proof.* Since item 4 holds, there exists $h > 0$ such that:

$$ \frac{L_F \|P'\| + \|K_0\| + L_\lambda \|K_1\|}{|V| \lambda_{\min}^2(P')} < h \le \frac{\lambda_{\min}(\widetilde{W} + \widetilde{W}^\top)\lambda_2(L)}{2|V|} $$

For each $t \ge 0$, let $v(t, \cdot)$ be the function $v(t, x) := [\alpha(t, \bar{x})^\top \dots \alpha(t, \bar{x})^\top]^\top \in \ker(\Theta_c^\top \otimes I_n)$, where $\alpha(t, \bar{x}) = h\bar{x} \in \mathbb{R}^n$, with $\bar{x} = (\mathbf{1}_{|V|}^\top \otimes I_n)P'x \in \mathbb{R}^n$. Let us consider the following inclusion:

$$ \begin{aligned} 0 \in F(t, x) - v(t, x) - B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)Cx \\ - B\Xi^{-1}(\Theta \otimes I_p)H'\mathbf{\Phi}'(N'\Delta^{-1}(\Theta^\top \otimes I_p)Cx). \end{aligned} \quad (48) $$

Proceeding similarly to the proof of Theorem 2, under Assumptions 2-4, this inclusion is rewritten as:

$$ 0 \in \tilde{F}(t, x) + \mathbf{\Psi}(x), \quad (49) $$

where $\mathbf{\Psi}(\cdot) = C^\top(\Theta \otimes I_p)(\Delta^{-1})^\top(N')^\top\mathbf{\Phi}'(N'\Delta^{-1}(\Theta^\top \otimes I_p)C\cdot)$ is a maximal monotone operator, and

$$ \begin{aligned} \tilde{F}(t, \cdot) = -P'F(t, \cdot) + P'v(t, \cdot) + P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)C(\cdot) \\ + (P'B\Xi^{-1}(\Theta \otimes I_p)H' - C^\top(\Theta \otimes I_p)(\Delta^{-1})^\top(N')^\top)\lambda^{\text{im}'}(\cdot), \end{aligned} $$

with $\lambda^{\text{im}'}(\cdot) := \text{proj}_{\text{rge}(\Gamma'+\Gamma'^{\top})}(\mathbf{\Phi}'(N'\Delta^{-1}(\Theta^\top \otimes I_p)C\cdot))$ is Lipschitz continuous [72, Lemma 3]. Since $f_i(t, \cdot)$ is Lipschitz continuous for all $i = \overline{1,|V|}$, proceeding similarly to the proof of Theorem 2, $F(t, \cdot)$ is Lipschitz continuous. Let any $x_1, x_2 \in \mathbb{R}^{|V|}$, then the following inequalities hold:

$$ -\langle P'(F(t, x_1) - F(t, x_2), x_1 - x_2 \rangle \ge -L_F(t) \|P'\| \|x_1 - x_2\|^2 \ge -L_F \|P'\| \|x_1 - x_2\|^2, \quad (50) $$
where $L_F(t)$ is a Lipschitz constant of $F(t, \cdot)$ and $L_F = \max_{t \in \mathbb{R}_+} L_F(t)$. Also:

$$ -\langle K_0(x_1 - x_2), x_1 - x_2 \rangle \geq -\|K_0\|\|x_1 - x_2\|^2, \quad (51) $$

$$ \langle K_1(\lambda^{im}(x_1) - \lambda^{im}(x_2)), x_1 - x_2 \rangle \geq -L_{\lambda} \|K_1\| \|x_1 - x_2\|^2. \quad (52) $$

Using item 3, the following holds:

$$ \begin{aligned} \langle P'B(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)C(x_1 - x_2), x_1 - x_2 \rangle & = \langle P'B(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)B^\top P'^\top(x_1 - x_2), x_1 - x_2 \rangle \\ & = \langle P'(\Theta \otimes I_n)\widetilde{W}(\Theta^\top \otimes I_n)P'(x_1 - x_2), x_1 - x_2 \rangle \\ & \geq \frac{\lambda_{min}(\widetilde{W} + \widetilde{W}^\top)}{2}(x_1 - x_2)^\top P'(\Theta\Theta^\top \otimes I_n)P'(x_1 - x_2) \\ & \geq \frac{\lambda_2(L)\lambda_{min}(\widetilde{W} + \widetilde{W}^\top)}{2|V|}(x_1 - x_2)^\top P'(\Theta_c\Theta_c^\top \otimes I_n)P'(x_1 - x_2) \\ & = \frac{\lambda_2(L)\lambda_{min}(\widetilde{W} + \widetilde{W}^\top)}{2|V|}\|(\Theta_c^\top \otimes I_n)P'(x_1 - x_2)\|^2, \end{aligned} \quad (53) $$

in which, the second equality follows from items 1 and 2, and the second inequality follows from (ii4) in Remark 15 (see Appendix B). In addition, we have:

$$ \begin{aligned} \langle P'(v(t, x_1) - v(t, x_2)), x_1 - x_2 \rangle &= \langle v(t, x_1) - v(t, x_2), P'(x_1 - x_2) \rangle \\ &= (\alpha(t, \bar{x}_1) - \alpha(t, \bar{x}_2))^\top (\mathbf{1}_{|V|}^\top \otimes I_n) P'(x_1 - x_2) \\ &= (\alpha(t, \bar{x}_1) - \alpha(t, \bar{x}_2))^\top (\bar{x}_1 - \bar{x}_2) \\ &= h \|\bar{x}_1 - \bar{x}_2\|^2 \\ &= h(x_1 - x_2)^\top P' (\mathbf{1}_{|V|} \mathbf{1}_{|V|}^\top \otimes I_n) P'(x_1 - x_2) \\ &= h(x_1 - x_2)^\top P' ( (|V|I_{|V|} - \Theta_c \Theta_c^\top) \otimes I_n) P'(x_1 - x_2) \\ &\geq h|V| \lambda_{min}^2(P') \|x_1 - x_2\|^2 - h \|(\Theta_c^\top \otimes I_n)P'(x_1 - x_2)\|^2, \end{aligned} \quad (54) $$

where the last equality follows from the fact that $I_{|V|} - \frac{1}{|V|} \mathbf{1}_{|V|} \mathbf{1}_{|V|}^\top = \frac{1}{|V|} \Theta_c \Theta_c^\top$. Using [73, Corollary 2.8.8] and the definition of $\Xi$ defined after (36), $\Xi^{-1}$ is given by $\Xi^{-1} = I_{p|V|} - (\Theta \otimes I_p)W\Delta^{-1}(\Theta^\top \otimes I_p)D$. Therefore, it follows from (50)-(54) that:

$$ \langle \tilde{F}(t, x_1) - \tilde{F}(t, x_2), x_1 - x_2 \rangle \geq \epsilon \|x_1 - x_2\|^2, \quad (55) $$

where $\epsilon := h|V|\lambda_{min}^2(P') - L_F\|P'\| - \|K_0\| - L_\lambda\|K_1\| > 0$. In addition, it can be checked that $v(t, \cdot)$ is Lipschitz continuous. Thus, $\tilde{F}(t, \cdot)$ is strongly monotone and Lipschitz continuous. The inclusion (49) has a unique solution $x^*(\cdot)$, using [79, Corollary 3.4]. Since (49) is equivalent to (48), $x^*(\cdot)$ is a solution of (48), so does (43). The conclusion follows from the equivalence of (43) and (27). $\square$

The following corollary provides necessary conditions for items in Proposition 6 to hold.

**Corollary 5.** Suppose that items 1, 2, 4 in Proposition 6 hold. Then the following hold:
1. For all $i \in \mathcal{V}$, $\ker(B_i^\top) = \{0\}$.

2. $p \ge n$, where $p, n$ are dimensions of the input and the state in (5), respectively.

**Proof.** 1. Assume contrarily that there exists $i \in \mathcal{V}$ such that $\ker(B_i^\top) \ne \{0\}$. Then there exists a vector $v = [\mathbf{0}_n^\top \dots \hat{v}^\top \dots \mathbf{0}_n^\top]^\top \in \mathbb{R}^{n|\mathcal{V}|}$ with $\hat{v} \in \ker(B_i^\top) \setminus \{0\}$ such that $(\tilde{W} + \tilde{W}^\top)v = 0$, where $\tilde{W}$ is defined in Proposition 6. Therefore, 0 is an eigenvalue of the matrix $\tilde{W} + \tilde{W}^\top$. This implies that the smallest eigenvalue of $\tilde{W} + \tilde{W}^\top$ satisfies $\lambda_{\min}(\tilde{W} + \tilde{W}^\top) \le 0$, which contradicts item 4 in Proposition 6.

2. Let any $i \in \mathcal{V}$. As indicated above, if items 1, 2, 4 in Proposition 6 hold, then $\ker(B_i^\top) = \{0\}$. This means that there does not exist a nonzero vector $x = [x_1 \dots x_n]^\top \in \mathbb{R}^n$ such that:

$$x_1 b_i^1 + x_2 b_i^2 + \dots + x_n b_i^n = 0,$$

where $B_i^\top = [b_i^1 \dots b_i^n] \in \mathbb{R}^{p \times n}$. It follows that $n$ vectors $b_i^1, \dots, b_i^n \in \mathbb{R}^p$ are linearly independent. In fact, in the space $\mathbb{R}^p$, there can be at most $p$ linearly independent vectors. Therefore, $p \ge n$. $\square$

**Remark 9.** Proposition 5 mainly relies on the structure of the functions $f_i(t, \cdot)$ while relaxing some conditions on the matrices $B_i$ and $C_i$, such as items 1, and 2 in Corollary 5. In contrast, Proposition 6 mainly relies on the structure of the matrices $B_i$ and $C_i$, while relaxing the conditions on the $f_i(t, \cdot)$, such as the QUAD property. There are many algorithms that allow us to find the solution $x^*(\cdot)$ of (49) where $\tilde{F}(t, \cdot)$ is strongly monotone and $\Psi$ is maximal monotone, see, e.g., [79, 80].

**Example 1.** Let us consider the following three systems for some $d_1, d_2, d_3 \in \mathbb{R}$:

$$\left\{ \begin{array}{l} \dot{x}_1 = -x_1 - 2 + u_1 \\ y_1 = x_1 + d_1 u_1 \end{array} \right. ; \left\{ \begin{array}{l} \dot{x}_2 = -(1+t)x_2 + u_2 \\ y_2 = x_2 + d_2 u_2 \end{array} \right. ; \left\{ \begin{array}{l} \dot{x}_3 = 1 - \sin(\frac{1}{t+1}x_3) + u_3 \\ y_3 = x_3 + d_3 u_3 \end{array} \right. . \tag{56}$$

Each agent in (56) has the form (5). They are an extension of the systems studied for synchronization in [?], where the considered coupling is a special case of the LCP coupling (a combination of linear and sign couplings). In this example, we will consider a network of the three agents in (56), in the case $d_1 = d_3 = 1, d_2 = 0$, with the graph shown in Figure 3:

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

$$\Theta = \begin{bmatrix} 1 & 1 \\ -1 & 0 \\ 0 & -1 \end{bmatrix}$$

<caption>Figure 3: Graph of the network in Example 1, and its oriented incidence matrix $\Theta$.</caption>
Suppose that two couplings in the network have the form (6) with the following data:

$$
\left\{
\begin{array}{l}
H_{12} = [1 - 1 2 - 2], \quad M_{12} = \operatorname{diag}(\mathbf{1}_2 \mathbf{1}_2^\top, I_2), \quad N_{12} = [1 - 1 1 - 1]^\top, \\
r_{12} = \left[\frac{-1}{2} \frac{-1}{2} 0 0\right]^\top, \quad H_{13} = [0 0 - 1 1], \quad M_{13} = I_4, \quad N_{13} = [0 0 1 - 1]^\top, \quad r_{13} = \mathbf{0}_4.
\end{array}
\right.
$$

It can be checked that these couplings have the special form (28) where $H'_{12} = [1 - 1]$,
$H'_{13} = \mathbf{0}_2^\top$, $N'_{12} = [1 - 1]^\top$, $N'_{13} = \mathbf{0}_2$, $W_{12} = 2, W_{13} = -1$, $M'_{12} = \mathbf{1}_2 \mathbf{1}_2^\top$, $M'_{13} = I_2$,
$r'_{12} = [\frac{-1}{2} \quad \frac{-1}{2}]^\top$, $r'_{13} = \mathbf{0}_2$. It follows from (37) that the overall dynamics of the network
is defined as the TVCS:

$$
\left\{
\begin{array}{l}
\dot{x}(t) = \begin{bmatrix} -x_1 - 2 \\ -(1+t)x_2 \\ 1 - \sin(\frac{1}{t+1}x_3) \end{bmatrix} - \begin{bmatrix} 1 & 0 & -1 \\ 0 & 2 & -2 \\ -1 & -2 & 3 \end{bmatrix} x(t) + \begin{bmatrix} 0 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \end{bmatrix} \lambda'(t), \\
\\
w'(t) = \begin{bmatrix} 0 & 0 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \end{bmatrix}^\top x(t) + \begin{bmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \lambda'(t) + \begin{bmatrix} \frac{-1}{2} \\ \frac{-1}{2} \\ 0 \\ 0 \end{bmatrix}, \\
\\
0 \le \lambda'(t) \perp w'(t) \ge 0.
\end{array}
\right.
\qquad (57)
$$

Let us check the validity of assumptions in Proposition 5.

1. Since $\Delta = \begin{bmatrix} 3 & -1 \\ 2 & -1 \end{bmatrix}$ and $\Gamma' = \operatorname{diag}(\begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}, I_2) \succcurlyeq 0$, Assumptions 1, 2 hold.

2. It can be checked that Assumption 3' holds with $P' = I_3$. Therefore, Assumption 3 holds.

3. Using (40) in Remark 7, $(\text{SOL}(0, \Gamma'))^* = \mathbb{R}^4$. Therefore, Assumption 4 holds.

4. Using Proposition 1, it holds that $F$ is QUAD($P'$, $Q$) where $Q = \operatorname{diag}(-1, -1, 1)$ and $P' = I_3$. Moreover, $P'B\Xi^{-1}\Theta W\Theta^\top C - Q = \begin{bmatrix} 2 & 0 & -1 \\ 0 & 3 & -2 \\ -1 & -2 & 2 \end{bmatrix} \succcurlyeq 0.09I_3$.

Using item 2 in Corollary 4 and Remark 8, all conditions in Proposition 5 are satisfied.
Then (27) has a solution $x^*(\cdot)$. Indeed, let $x^* = [-1 \ 0 \ 0]^\top$, and $\lambda'^* = [\frac{1}{4} \ \frac{1}{4} \ 0 \ 0]^\top$. It can
be checked that $(x^*, \lambda'^*)$ solves the TVCS (57), and therefore $x^*$ solves (43) and (13).
Following Proposition 5, since (43) is equivalent to (27), $x^*$ is a solution of (27). In
addition, since $x^*$ solves the TVCS (13), the point $e^* = [-1 \ -1 \ 0]^\top$ is an equilibrium
point of the error dynamics (21), using the definition given prior to (27). Note that
the well-posedness of (57) is guaranteed by Corollary 3.
**The structure of the couplings:** The couplings in this example are rewritten as:

$$ \left\{ \begin{array}{l} \eta_{12}(t) = [1 - 1]\lambda'_{12}(t) - 2v_{12}(t) \\ 0 \le \lambda'_{12}(t) \perp w'_{12}(t) = \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix} \lambda'_{12}(t) + \begin{bmatrix} 1 \\ -1 \end{bmatrix} v_{12}(t) + \begin{bmatrix} -\frac{1}{2} \\ -\frac{1}{2} \end{bmatrix} \ge 0. \end{array} \right. , \quad (58) $$

and $\eta_{13}(t) = v_{13}(t)$, using Lemma 1. In addition, it can be checked that:

$$ \eta'_{12}(t) = [1 - 1]\lambda'_{12}(t) \in -v_{12}(t) - \frac{1}{2}\operatorname{sgn}(v_{12}(t)). $$

Therefore, the couplings (6) used in this example are:

$$ \left\{ \begin{array}{l} \eta_{12}(t) \in -3v_{12}(t) - \frac{1}{2}\operatorname{sgn}(v_{12}(t)), \\ \eta_{13}(t) = v_{13}(t). \end{array} \right. \quad (59) $$

Table 1 presents the appropriate linear coupling matrices for this example for various triplets $(d_1, d_2, d_3)$. In these cases, all data of two couplings except for $W_{12}, W_{13}$, are retained as in the above case, that is, the case $(d_1, d_2, d_3) = (1, 0, 1)$. Denote $\Gamma' = \operatorname{diag}(\tilde{\Gamma}, I_2)$, and let $\mu$ be a positive number such that $P'B\Xi^{-1}\Theta W\Theta^\top C - Q \ge \mu I_3$. The values of $\tilde{\Gamma}, \mu$ can be found in Table 1.

<table><caption>Table 1: Linear coupling matrices, along with the corresponding results (Example 1).</caption><thead><tr><th></th><th>(0,0,0)</th><th>(0,1,0)</th><th>(0,0,1)</th><th>(0,1,1)</th><th>(0,0,-1)</th></tr></thead><tbody><tr><th>Results</th><th colspan="5" style="text-align:center;">(d<sub>1</sub>, d<sub>2</sub>, d<sub>3</sub>)</th></tr><tr><th>W<sub>12</sub></th><td>2</td><td>4</td><td>1</td><td>2</td><td>2</td></tr><tr><th>W<sub>13</sub></th><td>3</td><td>4</td><td>-1.2</td><td>-1.3</td><td>0.9</td></tr><tr><th>Γ̃</th><td><b>[1 1]</b><br><b>[1 1]</b></td><td><b>[6/5 4/5]</b><br><b>[4/5 6/5]</b></td><td><b>[1 1]</b><br><b>[1 1]</b></td><td><b>[4/3 2/3]</b><br><b>[2/3 4/3]</b></td><td><b>[1 1]</b><br><b>[1 1]</b></td></tr><tr><th>(SOL(0, Γ'))*</th><td><b>I</b>R<sup>4</sup></td><td><b>I</b>R<sup>4</sup></td><td><b>I</b>R<sup>4</sup></td><td><b>I</b>R<sup>4</sup></td><td><b>I</b>R<sup>4</sup></td></tr><tr><th>μ</th><td>0.04</td><td>0.035</td><td>0.105</td><td>0.03</td><td>0.19</td></tr><tr><th>Results</th><th colspan="5" style="text-align:center;">(d<sub>1</sub>, d<sub>2</sub>, d<sub>3</sub>)</th></tr><tr><th>W<sub>12</sub></th><td>5</td><td>1</td><td>2</td><td>8</td><td>10</td></tr><tr><th>W<sub>13</sub></th><td>5</td><td>-2</td><td>-1</td><td>-0.7</td><td>-0.115</td></tr><tr><th>Γ̃</th><td><b>[12/11 10/11]</b><br><b>[10/11 12/11]</b></td><td><b>[1 1]</b><br><b>[1 1]</b></td><td><b>[1 1]</b><br><b>[1 1]</b></td><td><b>[13/12 11/12]</b><br><b>[11/12 13/12]</b></td><td><b>[131/120 109/120]</b><br><b>[109/120 131/120]</b></td></tr><tr><th>(SOL(0, Γ'))*</th><td><b>I</b>R<sup>4</sup></td><td><b>I</b>R<sup>4</sup></td><td><b>I</b>R<sup>4</sup></td><td><b>I</b>R<sup>4</sup></td><td><b>I</b>R<sup>4</sup></td></tr><tr><th>μ</th><td>0.03</td><td>0.2</td><td>0.095</td><td>0.27</td><td>0.185</td></tr></tbody></table>

### 3.3 Equilibrium Point at the Origin of (21) under Dry Friction Condition

In this subsection, the couplings (6) with a special structure given in Lemma 1 are considered. Let us introduce the following assumption, which plays a crucial role in
proving the main results of this work.

**Assumption 6.** There exists $\rho > 0$ such that $\mathcal{B}_\rho \subset \text{int}(\mathcal{P}(0))$, where $\mathcal{B}_\rho = \{\gamma \in \mathbb{R}^{p|\mathcal{E}|} : \|\gamma\| \le \rho\}$, and $\mathcal{P}(\cdot) := (N'\Delta^{-1})^\top \mathbf{\Phi}'(N'\Delta^{-1})$.

This assumption is known as a condition of dry friction [81, 49]. For instance, the case $\mathcal{P} = \partial\varphi$ with $\varphi: \mathbb{R}^n \to \mathbb{R}$ being a convex function, is considered in [82, 83]. When $D=0$, it can be seen that $\Delta = I_{p|\mathcal{E}|}$ and the validity of Assumption 6 depends solely on a part of the structure of the couplings. The following proposition provides some necessary conditions for the validity of Assumption 6.

**Proposition 7.** Let Assumption 6 hold. Then $\Gamma'$ is not a $\mathbf{P}$-matrix and the following assertions are true:

1. If Assumption 2 also holds, then:

1.1. $r' \in (\text{SOL}(0, \Gamma'))^*$ in (34),

1.2. $\text{int}(\mathcal{P}(0)) = (N'\Delta^{-1})^\top \text{rint}(\mathbf{\Phi}'(0))$, where $\mathcal{P}$ is given in Assumption 6,

1.3. if either $r' \ne 0$ or $(\Gamma' + \Gamma'^\top) \ne 0$, then $\text{int}(\mathbf{\Phi}'(0)) = \emptyset$.

2. If $D, W$ are symmetric and $H' = N'^\top$, then

$$ \ker((N'\Delta^{-1})^\top) \cap \mathbf{\Phi}'(0) \cap (-\text{SOL}(r', M')) \ne \emptyset. $$

*Proof.* Proceeding similarly to (14), it can be checked that $\mathbf{\Phi}'(v) = -\text{SOL}(v+r', \Gamma')$, for all $v \in \mathbb{R}^{m'|\mathcal{E}|}$. Since Assumption 6 holds, $\mathcal{P}(0)$ is neither an empty set nor a singleton, and $\mathcal{B}_\rho \in \mathcal{P}(0)$, where $\mathcal{P}(\cdot) := (N'\Delta^{-1})^\top \mathbf{\Phi}'(N'\Delta^{-1})$, for some $\rho > 0$. This implies that $\Gamma'$ is not a $\mathbf{P}$-matrix, using Theorem 1, $\mathbf{\Phi}'(0)$ is not empty and $\text{int}(\mathcal{P}(0)) = \text{rint}(\mathcal{P}(0))$. Let us consider the following cases:

[Item 1.] Suppose that Assumption 2 holds.

Item 1.1: Using [56, Exercise 3.12.1, Corollary 3.8.10, Theorem 3.8.6], the nonemptiness of the set $\mathbf{\Phi}'(0)$ implies that $r' \in (\text{SOL}(0, \Gamma'))^*$. Since $\Gamma' \ge 0$, and $\mathbf{\Phi}'(0) \ne \emptyset$, let $\bar{z}$ be an arbitrary element of $\mathbf{\Phi}'(0)$. Using [56, Theorem 3.1.7], $\mathbf{\Phi}'(0)$ is a closed convex polyhedral set defined by:

$$ \mathbf{\Phi}'(0) = \{z \in \mathbb{R}^{m'|\mathcal{E}|} : r' - \Gamma'z \ge 0, r'^\top(z + \bar{z}) = 0, (\Gamma' + \Gamma'^\top)(z + \bar{z}) = 0\}. $$

Let $S_1 := \{z \in \mathbb{R}^{m'|\mathcal{E}|} : r'^\top(z + \bar{z}) = 0\}$, and $S_2 := \{z \in \mathbb{R}^{m'|\mathcal{E}|} : (\Gamma' + \Gamma'^\top)(z + \bar{z}) = 0\}$. Then,

$$ \mathbf{\Phi}'(0) \subseteq S_1 \cap S_2. \qquad (60) $$

Item 1.2: Using [54, Proposition 2.44], we obtain $\text{rint}(\mathcal{P}(0)) = (N'\Delta^{-1})^\top(\text{rint} \mathbf{\Phi}'(0))$. The conclusion follows.

Item 1.3: Let us consider the affine hull of $\mathbf{\Phi}'(0)$ in the following cases:

* If $r' \ne 0$, the set $S_1 = \{z \in \mathbb{R}^{m'|\mathcal{E}|} : r'^\top z = -r'^\top \bar{z}\} = -\bar{z} + \{z \in \mathbb{R}^{m'|\mathcal{E}|} : r'^\top z = 0\}$ is an affine space and $S_1 \subsetneq \mathbb{R}^{m'|\mathcal{E}|}$. It follows from (60) that

$$ \text{aff}(\mathbf{\Phi}'(0)) \subseteq S_1 \subsetneq \mathbb{R}^{m'|\mathcal{E}|}. $$
• If $(\Gamma' + \Gamma'^{\top}) \neq 0$, $\ker(\Gamma' + \Gamma'^{\top}) \subsetneq \mathbb{R}^{m'|\epsilon|}$. In addition, since $S_2 = -\bar{z} + \ker(\Gamma' + \Gamma'^{\top})$ is an affine space, $S_2 \subsetneq \mathbb{R}^{m'|\epsilon|}$. It follows from (60) that $\mathrm{aff}(\Phi'(0)) \subseteq S_2 \subsetneq \mathbb{R}^{m'|\epsilon|}$.

Hence $\mathrm{aff}(\Phi'(0)) \subsetneq \mathbb{R}^{m'|\epsilon|}$. Let us prove that $\mathrm{int}(\Phi'(0)) = \emptyset$. Suppose, on the contrary, that there exists $x \in \mathrm{int}(\Phi'(0))$. Since $\mathrm{int}(\Phi'(0))$ is an open set, there exists a closed ball $\mathcal{B}(x, \epsilon)$ satisfying $\mathcal{B}(x, \epsilon) \subseteq \mathrm{int}(\Phi'(0))$ for some $\epsilon > 0$. It follows that:

$$ \mathbb{R}^{m'|\epsilon|} = \mathrm{aff}(\mathcal{B}(x, \epsilon)) \subseteq \mathrm{aff}(\mathrm{int}(\Phi'(0))) \subseteq \mathrm{aff}(\Phi'(0)) \subsetneq \mathbb{R}^{m'|\epsilon|}, $$

which is absurd. The conclusion follows.

[Item 2.] Suppose that $D, W$ are symmetric and $H' = N'^{\top}$. Since $\mathcal{B}_{\rho} \subset \mathcal{P}(0)$, $0 \in \mathcal{P}(0)$. There exists $x \in \Phi'(0)$ such that $(N'\Delta^{-1})^{\top}x = 0$. This implies that $x \in \ker((N'\Delta^{-1})^{\top}) \cap \Phi'(0)$. In addition, using Remark 6, the following holds:

$$ \begin{aligned} \Gamma'x &= N'(\Theta^{\top} \otimes I_p)D\Xi^{-1}(\Theta \otimes I_p)H'x + M'x \\ &= N'(\Theta^{\top} \otimes I_p)D(\Theta \otimes I_p)(N'\Delta^{-1})^{\top}x + M'x = M'x. \end{aligned} $$

This implies that $x \in -\mathrm{SOL}(r', M')$, since $x \in \Phi'(0) = -\mathrm{SOL}(r', \Gamma')$. From all of the above, the conclusion follows. $\square$

As mentioned in Proposition 7, a necessary condition for Assumption 6 to hold is that $\Gamma'$ is not a $\mathbf{P}$-matrix. In addition, Assumption 2 is necessary to ensure the well-posedness of the overall dynamics (see Corollary 3), so $\Gamma' \ge 0$ is required. The following proposition provide sufficient conditions to ensure such properties of $\Gamma'$.

**Proposition 8.** Suppose that $D, W, M'$ are symmetric matrices, $H' = N'^{\top}$, and one of the following holds:

1. $D\Xi^{-1} \ge 0$, $M' \ge 0$, and $\ker(D_s(\Theta \otimes I_p)H') \cap \ker M' \ne \{0\}$, where $D_s = D_s^{\top} = (D\Xi^{-1})^{\frac{1}{2}} \ge 0$.

2. $\Gamma' \ge 0$ and $\ker H' \cap \ker M' \ne \{0\}$.

Then $\Gamma' \ge 0$ and it is not a $\mathbf{P}$-matrix.

*Proof.* Since $D, W$ are symmetric matrices, so $D\Xi^{-1}$ is symmetric, using Remark 6. In addition, since $N' = H'^{\top}$ and $M'$ is symmetric, $\Gamma'$ is a symmetric matrix.

[If item 1 holds] Since $D\Xi^{-1} \ge 0$, there exists a matrix $\hat{D} = \hat{D}^{\top} \ge 0$ such that $\hat{D}^2 = D\Xi^{-1}$. The following holds:

$$ \Gamma' = N'(\Theta^{\top} \otimes I_p)D\Xi^{-1}(\Theta \otimes I_p)H' + M' = (\hat{D}(\Theta \otimes I_p)H')^{\top}(\hat{D}(\Theta \otimes I_p)H') + M'. $$

Using [73, Theorem 2.4.3], $\ker((\hat{D}(\Theta \otimes I_p)H')^{\top}(\hat{D}(\Theta \otimes I_p)H')) = \ker(\hat{D}(\Theta \otimes I_p)H')$. In addition, by the assumption, $\ker(\hat{D}(\Theta \otimes I_p)H') \cap \ker M' \ne \{0\}$. Therefore, there exists $x \in \mathbb{R}^{m'|\epsilon|} \setminus \{0\}$ such that $\Gamma'x = 0$. This means that zero is an eigenvalue of $\Gamma'$. It follows from the symmetry of $\Gamma'$ that $\Gamma'$ is not a $\mathbf{P}$-matrix, using [56, Theorem 3.3.4 (c)]. Moreover, since $M' \ge 0$, so $\Gamma' \ge 0$.

[If item 2 holds] Since $\ker H' \cap \ker M' \ne \{0\}$, there exists $x \in \mathbb{R}^{m'|\epsilon|} \setminus \{0\}$ such that $\Gamma'x = 0$. Proceeding similarly to the first case, $\Gamma'$ is not a $\mathbf{P}$-matrix. By assumption, $\Gamma' \ge 0$. $\square$
### 3.3.1 Methodology for Coupling Design under Assumption 6 when $D = 0$

Since all conditions in Lemma 1 are assumed to hold, $q = p$ where $q$ and $p$ are the dimensions of the input $u_i(t)$ and the output $y_i(t)$ defined in (5), respectively. The following proposition provides an explicit design for $N'$, $M'$, $r'$ in the coupling that ensures Assumption 6 holds when $D = 0$ and the dimension of the input $u_i(t)$ defined in (5) is one (i.e., $q = p = 1$). The case $q = p > 1$ (resp. $D \neq 0$) is addressed in Remark 10 (resp. Proposition 10). Recall that $N'$, $M'$, $r'$ are given in (34).

**Proposition 9.** Suppose that $D = 0$, all conditions in Lemma 1 hold, and $q = p = 1$.

Assume further that for all pairs $(i, j)$ where $i < j$, $(i, j) \in \mathcal{E}$, then $N'_{ij} = \begin{bmatrix} n_{ij,1} \\ n_{ij,2} \end{bmatrix}$,

$M'_{ij} = \begin{bmatrix} m_{ij,1} & m_{ij,2} \\ m_{ij,2} & m_{ij,3} \end{bmatrix}$, and $r'_{ij} = \begin{bmatrix} r_{ij,1} \\ r_{ij,2} \end{bmatrix}$ are such that:

1. Structure of $N'_{ij}$: $n_{ij,1} > 0$ and $n_{ij,2} < 0$.

2. Structure of $M'_{ij}$: $m_{ij,1}, m_{ij,3} \ge 0, (m_{ij,1})^2 + (m_{ij,3})^2 > 0$, and $m_{ij,1}m_{ij,3} = (m_{ij,2})^2$.

3. Structure of $r'_{ij}$: defined as in Table 2.

Then $-N'_{ij}^\top \text{SOL}(r'_{ij}, M'_{ij}) \ni \mathcal{B}_{\rho_{ij}}$, Assumptions 2 and 6 hold for all $\rho < \min_{i<j, (i,j) \in \mathcal{E}} \rho_{ij}$, where $\rho_{ij}$ is given in Table 3.

<table><caption>Table 2: Structure of $r'_{ij}$ based on the given structures of $N'_{ij}$ and $M'_{ij}$.</caption><thead><tr><th rowspan="3">Cases</th><th rowspan="3">$m_{ij,1} = 0$</th><th rowspan="3">$m_{ij,3} = 0$</th><th colspan="3">$m_{ij,1} > 0, m_{ij,3} > 0$</th></tr><tr><th rowspan="2">$m_{ij,2} > 0$</th><th colspan="2">$m_{ij,2} < 0$</th></tr><tr><th>$n_{ij,1} > n_{ij,2} \frac{m_{ij,1}}{m_{ij,2}}$</th><th>$n_{ij,1} < n_{ij,2} \frac{m_{ij,1}}{m_{ij,2}}$</th></tr></thead><tbody><tr><th rowspan="2">$r'_{ij}$</th><th>$r_{ij,1} = 0$</th><th>$r_{ij,1} < 0$</th><th>$r_{ij,1} < 0$</th><th>$r_{ij,1} > 0$</th><th>$r_{ij,1} < 0$</th></tr><tr><th>$r_{ij,2} < 0$</th><th>$r_{ij,2} = 0$</th><th>$r_{ij,2} = \frac{m_{ij,2}}{m_{ij,1}} r_{ij,1}$</th><th>$r_{ij,2} = \frac{m_{ij,2}}{m_{ij,1}} r_{ij,1}$</th><th>$r_{ij,2} = \frac{m_{ij,2}}{m_{ij,1}} r_{ij,1}$</th></tr></tbody></table>

<table><caption>Table 3: Value $\rho_{ij}$ associated with the given structures of $N'_{ij}$, $M'_{ij}$, $r'_{ij}$.</caption><thead><tr><th rowspan="3">Cases</th><th rowspan="3">$m_{ij,1} = 0$</th><th rowspan="3">$m_{ij,3} = 0$</th><th colspan="3">$m_{ij,1} > 0, m_{ij,3} > 0$</th></tr><tr><th rowspan="2">$m_{ij,2} > 0$</th><th colspan="2">$m_{ij,2} < 0$</th></tr><tr><th>$n_{ij,1} > n_{ij,2} \frac{m_{ij,1}}{m_{ij,2}}$</th><th>$n_{ij,1} < n_{ij,2} \frac{m_{ij,1}}{m_{ij,2}}$</th></tr></thead><tbody><tr><th>$\rho_{ij}$</th><td>$\frac{n_{ij,2}r_{ij,2}}{m_{ij,3}}$</td><td>$\frac{|n_{ij,1}r_{ij,1}|}{m_{ij,1}}$</td><td>$\min\left(\frac{n_{ij,2}r_{ij,1}}{m_{ij,2}}, \frac{|n_{ij,1}r_{ij,1}|}{m_{ij,1}}\right)$</td><td>$\frac{n_{ij,2}r_{ij,1}}{m_{ij,2}}$</td><td>$\frac{|n_{ij,1}r_{ij,1}|}{m_{ij,1}}$</td></tr></tbody></table>

*Proof.* Since $D = 0$, $\Delta = I_{p|\mathcal{E}|}$ and the operator $\mathcal{P}$ in Assumption 6 is defined by $\mathcal{P}(\cdot) = N'^\top \mathbf{\Phi}'(N')$. Using (3), it can be checked that $\mathbf{\Phi}'(v) = -\text{SOL}(v + r', M')$, for
all $v \in \mathbb{R}^{m'|\mathcal{E}|}$. Therefore, $\mathcal{P}(0) = -N^T \text{SOL}(r', M')$. For each pair $(i, j) \in \mathcal{V} \times \mathcal{V}$ where $i < j$, $(i, j) \in \mathcal{E}$, let us consider the following LCP:

$$ \text{LCP}(r'_{ij}, M'_{ij}) : \quad 0 \le \lambda_{ij} \perp M'_{ij} \lambda_{ij} + r'_{ij} \ge 0, $$

where $\lambda_{ij} = [\lambda_{ij,1} \ \lambda_{ij,2}]^T \in \mathbb{R}^2$ and $M'_{ij}, r'_{ij}$ are given in items 2 and 3. It can be checked that there exists $\rho_{ij} > 0$ such that $\mathcal{B}_{\rho_{ij}} \in -N'_{ij}^T \text{SOL}(r'_{ij}, M'_{ij})$. Indeed, let us consider the following cases:

1. $[m_{ij,1} = 0, m_{ij,3} > 0$, and $r_{ij,1} = 0, r_{ij,2} < 0]$ The problem LCP($r'_{ij}, M'_{ij}$) is rewritten as:

$$ 0 \le \begin{bmatrix} \lambda_{ij,1} \\ \lambda_{ij,2} \end{bmatrix} \perp \begin{bmatrix} 0 \\ m_{ij,3}\lambda_{ij,2} + r_{ij,2} \end{bmatrix} \ge 0. $$

It follows that $\text{SOL}(r'_{ij}, M'_{ij}) = [\mathbb{R}_+ \left(-\frac{r_{ij,2}}{m_{ij,3}}\right)]^T$. Let $\rho_{ij} = \frac{n_{ij,2}r_{ij,2}}{m_{ij,3}} > 0$, it follows from item 1 that:

$$ -N'_{ij}^T \text{SOL}(r'_{ij}, M'_{ij}) = -n_{ij,1} \mathbb{R}_+ + \frac{n_{ij,2}r_{ij,2}}{m_{ij,3}} = \mathbb{R}_- + \rho_{ij} = (-\infty, \rho_{ij}] \ni \mathcal{B}_{\rho_{ij}}. $$

2. $[m_{ij,3} = 0, m_{ij,1} > 0$, and $r_{ij,2} = 0, r_{ij,1} < 0]$ The problem LCP($r'_{ij}, M'_{ij}$) is rewritten as:

$$ 0 \le \begin{bmatrix} \lambda_{ij,1} \\ \lambda_{ij,2} \end{bmatrix} \perp \begin{bmatrix} m_{ij,1}\lambda_{ij,1} + r_{ij,1} \\ 0 \end{bmatrix} \ge 0. $$

Thus $\text{SOL}(r'_{ij}, M'_{ij}) = \left[ \left(-\frac{r_{ij,1}}{m_{ij,1}}\right) \mathbb{R}_+ \right]^T$. Let $\rho_{ij} = \left|\frac{n_{ij,1}r_{ij,1}}{m_{ij,1}}\right| > 0$, it follows from assumption item 1 that:

$$ -N'_{ij}^T \text{SOL}(r'_{ij}, M'_{ij}) = -n_{ij,2} \mathbb{R}_+ + \frac{n_{ij,1}r_{ij,1}}{m_{ij,1}} = \mathbb{R}_+ - \rho_{ij} = [-\rho_{ij}, +\infty) \ni \mathcal{B}_{\rho_{ij}}. $$

3. $[m_{ij,1}, m_{ij,3}, m_{ij,2} > 0$, and $r_{ij,1} < 0, r_{ij,2} = \frac{m_{ij,2}}{m_{ij,1}}r_{ij,1}]$ The LCP($r'_{ij}, M'_{ij}$) is rewritten as:

$$ 0 \le \begin{bmatrix} \lambda_{ij,1} \\ \lambda_{ij,2} \end{bmatrix} \perp \begin{bmatrix} m_{ij,1}\lambda_{ij,1} + m_{ij,2}\lambda_{ij,2} + r_{ij,1} \\ m_{ij,2}\lambda_{ij,1} + m_{ij,3}\lambda_{ij,2} + \frac{m_{ij,2}}{m_{ij,1}}r_{ij,1} \end{bmatrix} \ge 0. \quad (61) $$

Since $m_{ij,2} > 0, m_{ij,3} > 0$ by assumption and $m_{ij,3} = \frac{(m_{ij,2})^2}{m_{ij,1}}$ as assumption item 2 holds, the LCP (61) is equivalent to:

$$ 0 \le \begin{bmatrix} \lambda_{ij,1} \\ \lambda_{ij,2} \end{bmatrix} \perp \begin{bmatrix} \frac{m_{ij,1}}{m_{ij,2}}\lambda_{ij,1} + \lambda_{ij,2} + \frac{r_{ij,1}}{m_{ij,2}} \\ \frac{m_{ij,1}}{m_{ij,2}}\lambda_{ij,1} + \lambda_{ij,2} + \frac{r_{ij,1}}{m_{ij,2}} \end{bmatrix} \ge 0. \quad (62) $$
Since $\frac{r_{ij,1}}{m_{ij,2}} < 0$, $(\lambda_{ij,1})^2 + (\lambda_{ij,2})^2 > 0$. It follows from the LCP (62) that: $\lambda_{ij,2} = -\frac{m_{ij,1}}{m_{ij,2}}\lambda_{ij,1} - \frac{r_{ij,1}}{m_{ij,2}}$. The set of solutions of the LCP (62) is given by:

$$ \mathrm{SOL}(r'_{ij}, M'_{ij}) = \left\{ \begin{bmatrix} \lambda_{ij,1} \\ \lambda_{ij,2} \end{bmatrix} \in \mathbb{R}_+^2 : \lambda_{ij,2} = -\frac{m_{ij,1}}{m_{ij,2}}\lambda_{ij,1} - \frac{r_{ij,1}}{m_{ij,2}} \right\}. \quad (63) $$

Let $\rho_{ij} = \min\left(\frac{n_{ij,2}r_{ij,1}}{m_{ij,2}}, \left|\frac{n_{ij,1}r_{ij,1}}{m_{ij,1}}\right|\right) > 0$. Since item 1 holds, $n_{ij,2}\frac{m_{ij,1}}{m_{ij,2}} - n_{ij,1} < 0$. Thus the following holds:

$$ -N_{ij}^{\top} \mathrm{SOL}(r'_{ij}, M'_{ij}) = -n_{ij,1}\lambda_{ij,1} - n_{ij,2}\lambda_{ij,2} = \left(n_{ij,2}\frac{m_{ij,1}}{m_{ij,2}} - n_{ij,1}\right)\lambda_{ij,1} + \frac{n_{ij,2}r_{ij,1}}{m_{ij,2}}. $$

It can be checked that $\lambda_{ij,1} \in [0, -\frac{r_{ij,1}}{m_{ij,1}}]$. This implies that:

$$ -N_{ij}^{\top} \mathrm{SOL}(r'_{ij}, M'_{ij}) \ge \left(n_{ij,2}\frac{m_{ij,1}}{m_{ij,2}} - n_{ij,1}\right) \left(-\frac{r_{ij,1}}{m_{ij,1}}\right) + n_{ij,2}\frac{r_{ij,1}}{m_{ij,2}} = \frac{n_{ij,1}r_{ij,1}}{m_{ij,1}}. $$

$$ -N_{ij}^{\prime\top} \mathrm{SOL}(r'_{ij}, M'_{ij}) \le \frac{n_{ij,2}r_{ij,1}}{m_{ij,2}}. $$

$$ \text{Hence, } -N_{ij}^{\prime\top} \mathrm{SOL}(r'_{ij}, M'_{ij}) = \left[ \frac{n_{ij,1}r_{ij,1}}{m_{ij,1}}, \frac{n_{ij,2}r_{ij,1}}{m_{ij,2}} \right] \ni \mathcal{B}_{\rho_{ij}}. $$

4. $[m_{ij,1}, m_{ij,3}, -m_{ij,2} > 0, n_{ij,1} > n_{ij,2}\frac{m_{ij,1}}{m_{ij,2}}, \text{ and } r_{ij,1} > 0, r_{ij,2} = \frac{m_{ij,2}}{m_{ij,1}}r_{ij,1}]$ The problem LCP($r'_{ij}, M'_{ij}$) is rewritten as (61). Since $-m_{ij,2}, m_{ij,3} > 0$ by assumption and $m_{ij,3} = \frac{(m_{ij,2})^2}{m_{ij,1}}$ from item 2, the LCP (61) is equivalent to:

$$ 0 \le \begin{bmatrix} \lambda_{ij,1} \\ \lambda_{ij,2} \end{bmatrix} \perp \begin{bmatrix} -\frac{m_{ij,1}}{m_{ij,2}}\lambda_{ij,1} - \lambda_{ij,2} - \frac{r_{ij,1}}{m_{ij,2}} \\ \frac{m_{ij,1}}{m_{ij,2}}\lambda_{ij,1} + \lambda_{ij,2} + \frac{r_{ij,1}}{m_{ij,2}} \end{bmatrix} \ge 0. \quad (64) $$

It follows that $\lambda_{ij,2} = -\frac{m_{ij,1}}{m_{ij,2}}\lambda_{ij,1} - \frac{r_{ij,1}}{m_{ij,2}} > 0$. The set of solutions of the LCP (64) is given by (63). Let $\rho_{ij} = \frac{n_{ij,2}r_{ij,1}}{m_{ij,2}} > 0$. Since $n_{ij,2}\frac{m_{ij,1}}{m_{ij,2}} - n_{ij,1} < 0$, the following holds:

$$ \begin{aligned} -N_{ij}^{\top} \mathrm{SOL}(r'_{ij}, M'_{ij}) &= \left(n_{ij,2}\frac{m_{ij,1}}{m_{ij,2}} - n_{ij,1}\right)\lambda_{ij,1} + \frac{n_{ij,2}r_{ij,1}}{m_{ij,2}} \\ &= \mathbb{R}_- + \frac{n_{ij,2}r_{ij,1}}{m_{ij,2}} \\ &= (-\infty, \rho_{ij}] \ni \mathcal{B}_{\rho_{ij}}. \end{aligned} $$

5. $[m_{ij,1}, m_{ij,3}, -m_{ij,2} > 0, n_{ij,1} < n_{ij,2}\frac{m_{ij,1}}{m_{ij,2}}, \text{ and } r_{ij,1} < 0, r_{ij,2} = \frac{m_{ij,2}}{m_{ij,1}}r_{ij,1}]$ Following a similar procedure to item 4 in the proof, the problem LCP($r'_{ij}, M'_{ij}$) is rewritten as (64). It also implies that $\lambda_{ij,2} = -\frac{m_{ij,1}}{m_{ij,2}}\lambda_{ij,1} - \frac{r_{ij,1}}{m_{ij,2}}$. Therefore $\lambda_1 \ge -\frac{r_{ij,1}}{m_{ij,1}}$ and the set of solution of this LCP is given by (63). Let $\rho_{ij} = \left|\frac{n_{ij,1}r_{ij,1}}{m_{ij,1}}\right| > 0$. Since $n_{ij,2}\frac{m_{ij,1}}{m_{ij,2}} - n_{ij,1} > 0$, the following hold:

$$ \begin{aligned} -N_{ij}^{\top} \mathrm{SOL}(r'_{ij}, M'_{ij}) &= \left(n_{ij,2}\frac{m_{ij,1}}{m_{ij,2}} - n_{ij,1}\right)\lambda_{ij,1} + \frac{n_{ij,2}r_{ij,1}}{m_{ij,2}} \\ &\ge \left(n_{ij,2}\frac{m_{ij,1}}{m_{ij,2}} - n_{ij,1}\right)\left(-\frac{r_{ij,1}}{m_{ij,1}}\right) + n_{ij,2}\frac{r_{ij,1}}{m_{ij,2}} = \frac{n_{ij,1}r_{ij,1}}{m_{ij,1}}. \end{aligned} $$

$$ \text{Hence, } -N_{ij}^{\prime\top} \mathrm{SOL}(r'_{ij}, M'_{ij}) = [-\rho_{ij}, +\infty) \ni \mathcal{B}_{\rho_{ij}}. $$
It follows from the form of $N', M', r'$ defined in (34) that:

$$ \begin{aligned} \mathcal{P}(0) = -N'^{\top}\mathrm{SOL}(r', M') &= \begin{bmatrix} -N'_{ij}^{\top} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & -N'_{kl}^{\top} \end{bmatrix} \begin{bmatrix} \mathrm{SOL}(r'_{ij}, M'_{ij}) \\ \vdots \\ \mathrm{SOL}(r'_{kl}, M'_{kl}) \end{bmatrix} \\ &= \begin{bmatrix} -N'_{ij}^{\top}\mathrm{SOL}(r'_{ij}, M'_{ij}) \\ \vdots \\ -N'_{kl}^{\top}\mathrm{SOL}(r'_{kl}, M'_{kl}) \end{bmatrix}, \end{aligned} $$

for some pairs $(i, j), (k, l) \in \mathcal{V} \times \mathcal{V}$ where $i < j, k < l$, and $(i, j), (k, l) \in \mathcal{E}$. Let $\rho_{min} = \min_{i<j,(i,j) \in \mathcal{E}} \rho_{ij}$, then $\mathcal{P}(0) \ni \mathcal{B}_{\rho_{min}}$. Hence $\mathrm{int}(\mathcal{P}(0)) \ni \mathcal{B}_{\rho}$ for all $\rho < \rho_{min}$. This means that Assumption 6 holds for all $\rho < \rho_{min}$. In addition, from item 2, it can be checked that the symmetric matrix $M'_{ij}$ is positive semidefinite for all pair $(i, j)$ where $i < j$, and $(i, j) \in \mathcal{E}$. Therefore $\Gamma' = M' \succcurlyeq 0$. This means that Assumption 2 holds. $\square$

**Remark 10.** When $D = 0$ and the dimension of the input $u_i(t)$ defined in (5) is greater than one (i.e., $q = p > 1$), for each pair $(i, j) \in \mathcal{V} \times \mathcal{V}$ where $i < j$, and $(i, j) \in \mathcal{E}$, let $N'_{ij} = \mathrm{diag}(N'_{ij,1}, \dots, N'_{ij,p})$, $M'_{ij} = \mathrm{diag}(M'_{ij,1}, \dots, M'_{ij,p})$, and $r'_{ij} = [r'_{ij,1}^{\top}, \dots, r'_{ij,p}^{\top}]^{\top}$ such that for all $k = 1, \dots, p$, the matrices $N'_{ij,k}$, $M'_{ij,k}$ and $r'_{ij,k}$ satisfy items 1, 2 and 3 in Proposition 9, respectively. Then an identical conclusion to that of Proposition 9 is obtained, where Assumption 6 holds for all $\rho < \min_{i<j,(i,j) \in \mathcal{E}} \rho_{ij,k}$, where $\rho_{ij,k}$ is given as in Table 3, but here it corresponds to the structures of $N'_{ij,k}$, $M'_{ij,k}$ and $r'_{ij,k}$. Note that in this case, $m' = 2p = 2q$, where $m'$ defined in Lemma 1.

Note that Proposition 9 only provides a partial structure of the coupling given in (34), in which the data $H'$ is missing. However, the condition on $H'$ is implicitly encoded in Assumption 3, where it guarantees the well-posedness of the network dynamics (37). The method for computing $H'$ to ensure that Assumption 3 holds, will be presented in Algorithm 1 in Section 4.2. The symmetry of the matrix $M'_{ij}$ in the above proposition can be relaxed. In that case, by following a similar approach to Proposition 9, it is possible to construct a class of matrices $N'_{ij}$, $M'_{ij}$, and $r'_{ij}$ which satisfy Assumption 6, but the details are omitted to avoid excessive calculations.

### 3.3.2 Methodology for Coupling Design under Assumption 6 when $D \neq 0$

Let $\hat{D} := (\Theta^{\top} \otimes I_p)D(\Theta \otimes I_p)$. Recall that $\Xi = I_{p|\mathcal{V}|} + (\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)D$ and $\Delta = I_{p|\mathcal{E}|} + \hat{D}W$. The proposition below is a direct extension of Proposition 9, and Remark 10 to the case $D \neq 0$, where the network structure is constrained by a condition imposed on the oriented incidence matrix $\Theta$.

**Proposition 10.** Suppose that all conditions in Lemma 1 hold with $m' = 2p$, and the following conditions hold:
1. $N'\Delta^{-1} = \text{diag}(N'_1, ..., N'_{|\mathcal{E}|})$, where $N'_k = \text{diag}(N'_{k,1}, ..., N'_{k,p})$ with $N'_{k,l}$ satisfying item 1 in Proposition 9, for all $k = \overline{1, |\mathcal{E}|}, l = \overline{1, p}$.

2. $N'\Delta^{-1}\hat{D}H' + M' = \text{diag}(\Gamma'_1, ..., \Gamma'_{|\mathcal{E}|})$, where $\Gamma'_k = \text{diag}(\Gamma'_{k,1}, ..., \Gamma'_{k,p})$ with $\Gamma'_{k,l}$ satisfying item 2 in Proposition 9, for all $k = \overline{1, |\mathcal{E}|}, l = \overline{1, p}$.

3. $r' = [r'_1^\top ... r'_{|\mathcal{E}|}^\top]^\top$, where $r'_k = [r'_{k,1}^\top ... r'_{k,p}^\top]^\top$ with $r'_{k,l}$ satisfying item 3 in Proposition 9, associated with the structures of $\Gamma'_{k,l}$, and $N'_{k,l}$, for all $k = \overline{1, |\mathcal{E}|}, l = \overline{1, p}$.

Then Assumption 2, and Assumption 6 hold for all $\rho < \min_{\substack{l=1,...,p \\ k=1,...,|\mathcal{E}|}} \rho_{k,l}$, where $\rho_{k,l}$ is given in Table 3, and is associated with the structures of $\Gamma'_{k,l}$, $N'_{k,l}$ and $r'_{k,l}$.

*Proof.* Using [84, Fact 3.20.6], it can be verified that:

$$N'(\Theta^\top \otimes I_p)D\Xi^{-1}(\Theta \otimes I_p)H' = N'\Delta^{-1}\hat{D}H'$$

Using the definition of $\Gamma'$ (see after (37)) and item 2, $\Gamma' = \text{diag}(\Gamma'_1, ..., \Gamma'_{|\mathcal{E}|}) \ge 0$. Hence Assumption 2 holds. Using (3), it can be checked that $\Phi'(0) = -\text{SOL}(r', \Gamma')$. Therefore it follows from items 1, 2, 3 that:

$$ \begin{aligned} \mathcal{P}(0) = -(N'\Delta^{-1})^\top \text{SOL}(r', \Gamma') &= \begin{bmatrix} -N'_1^\top \text{SOL}(r'_1, \Gamma'_1) \\ \vdots \\ -N'_{|\mathcal{E}|}^\top \text{SOL}(r'_{|\mathcal{E}|}, \Gamma'_{|\mathcal{E}|}) \end{bmatrix}, \\ N'_k^\top \text{SOL}(r'_k, \Gamma'_k) &= \begin{bmatrix} N'_{k,1}^\top \text{SOL}(r'_{k,1}, \Gamma'_{k,1}) \\ \vdots \\ N'_{k,p}^\top \text{SOL}(r'_{k,p}, \Gamma'_{k,p}) \end{bmatrix}, \forall k = 1, ..., |\mathcal{E}|, \end{aligned} $$

where $\mathcal{P}$ is defined in Assumption 6. For all $k = \overline{1, |\mathcal{E}|}$ and $l = \overline{1, p}$, since $N'_{k,l}$, $\Gamma'_{k,l}$, $r'_{k,l}$ satisfy items 1, 2, and 3 in Proposition 9 (by assumption), respectively, proceeding similarly to the proof of Proposition 9, $\mathcal{B}_{\rho_{k,l}} \in -N'_{k,l}^\top \text{SOL}(r'_{k,l}, \Gamma'_{k,l})$, where $\rho_{k,l}$ is given in Table 3, and is associated with the structures of $\Gamma'_{k,l}$, $N'_{k,l}$ and $r'_{k,l}$. The conclusion follows. $\square$

Recall that the matrix $H'$ (defined in (34)) is block diagonal with each block on the diagonal being a matrix in $\mathbb{R}^{p \times m'}$. From this, we observe that if item 1 in Proposition 10 holds and the matrix $\hat{D}$ is block diagonal with each block on the diagonal being a matrix in $\mathbb{R}^{p \times p}$, then there always exists $M'$ such that item 2 in Proposition 10 holds. The following proposition provides some sufficient conditions to ensure that item 1 in Proposition 10 holds, and it plays a role in the coupling design algorithm, presented in Section 4.2, to ensure the synchronization of the network.

**Proposition 11.** Suppose that $\hat{D} = (\Theta^\top \otimes I_p)D(\Theta \otimes I_p)$ is block diagonal, that is $\hat{D} = \text{diag}(\hat{D}_1, ..., \hat{D}_{|\mathcal{E}|})$ such that $\hat{D}_i \in \mathbb{R}^{p \times p}$ for all $i = 1, ..., |\mathcal{E}|$. The following assertions hold:
1. If $W = \text{diag}(W_1, ..., W_{|\mathcal{E}|})$ such that for all $k = 1, ..., |\mathcal{E}|$:

(a) $(I_p + W_k \hat{D}_k)^{-1}$ exists,

(b) $\lambda_{min}^k := \lambda_{min}((I_p + W_k \hat{D}_k)^{-1}W_k + W_k^\top (I_p + W_k \hat{D}_k)^{-\top}) > 0$.

Then Assumption 1 holds, and $\lambda_{min}(\bar{\Delta}^{-1}W + W^\top \bar{\Delta}^{-\top}) = \min_{k=1,...,|\mathcal{E}|} (\lambda_{min}^k) > 0$, where $\bar{\Delta} := I_{p|\mathcal{E}|} + W \hat{D}$.

2. If Assumption 1 holds with $W = \text{diag}(W_1, ..., W_{|\mathcal{E}|})$, and $N' = \text{diag}(N'_1, ..., N'_{|\mathcal{E}|})$ such that:

$$ \begin{array}{ll} N'_{k}(I_p + \hat{D}_k W_k)^{-1} = \text{diag}(N'_{k,1}, ..., N'_{k,p}), & \forall k = 1, ..., |\mathcal{E}|, \\ N'_{k,l} \text{ satisfies item 1 in Proposition 9}, & \forall k = 1, ..., |\mathcal{E}|, \forall l = 1, ..., p. \end{array} $$

Then item 1 in Proposition 10 holds.

Proof. [Item 1.] Since $W = \text{diag}(W_1, ..., W_{|\mathcal{E}|})$ (see definition of $W$ after (34)), and $\hat{D} = \text{diag}(\hat{D}_1, ..., \hat{D}_{|\mathcal{E}|})$ where $W_k, \hat{D}_k \in \mathbb{R}^{p \times p}$ for all $k = 1, ..., |\mathcal{E}|$, the following holds:

$$ \bar{\Delta} = I_p + W \hat{D} = \text{diag}((I_p + W_1 \hat{D}_1), ..., (I_p + W_{|\mathcal{E}|} \hat{D}_{|\mathcal{E}|})). $$

Using [84, Lemma 3.9.2], it follows from the assumption (a) that $\bar{\Delta}$ is invertible and

$$ \bar{\Delta}^{-1} = \text{diag}((I_p + W_1 \hat{D}_1)^{-1}, ..., (I_p + W_{|\mathcal{E}|} \hat{D}_{|\mathcal{E}|})^{-1}). $$

This implies that Assumption 1 holds [84, Fact 3.20.6], and the following holds:

$$ \bar{\Delta}^{-1} W = \text{diag}((I_p + W_1 \hat{D}_1)^{-1} W_1, ..., (I_p + W_{|\mathcal{E}|} \hat{D}_{|\mathcal{E}|})^{-1} W_{|\mathcal{E}|}). $$

Using [85, Problem 1.1.P4], the conclusion follows from the assumption (b).

[Item 2.] Proceeding similarly to the proof of item 1, we obtain:

$$ N' \Delta^{-1} = N'(I_{p|\mathcal{E}|} + \hat{D}W)^{-1} = \text{diag}(N'_1(I_p + \hat{D}_1 W_1)^{-1}, ..., N'_{|\mathcal{E}|}(I_p + \hat{D}_{|\mathcal{E}|}W_{|\mathcal{E}|})^{-1}), $$

The conclusion follows from the assumption. $\square$

Note that the matrix $\bar{\Delta}$ in item 1 in Proposition 11 will play an important role in Theorem 3 (main result), and $\bar{\Delta} = \Delta^\top$ if $D$ and $W$ are symmetric matrices. Similar to the remark at the end of Section 3.3.1, Propositions 10 and 11 only provide a partial structure of the coupling (comprising $W, N', M', r'$) given in (34), in which the data $H'$ is missing. It can be seen from the proof of Proposition 10 that $\Gamma' = N'\Delta^{-1}\hat{D}H' + M'$. Therefore, the conditions on $H'$ are implicitly encoded in item 2 in Proposition 10 and Assumptions 2, 3, 4 (through the matrix $\Gamma'$), which guarantee the well-posedness of the network dynamics (37). The method for computing $H'$ to ensure that such conditions hold, will be presented in Algorithm 1 in Section 4.2.
If the research framework is extended to the case where the network structure is subject to design (i.e., $\Theta$ can be freely chosen), the following proposition asserts that there always exists an oriented incidence matrix $\Theta$ associated with the graph $G = (V, E)$ (characterizing the network structure) satisfying (12) such that the matrix $\hat{D}$ is block diagonal with blocks on the diagonal in $\mathbb{R}^{p \times p}$. It also provides a rule for constructing such $\Theta$, presented immediately after the proof of the proposition.

**Proposition 12.** Suppose that there exists $\bar{n} \in \{1, ..., |V|\}$ such that $D_i = 0$ for all $i = 1, ..., \bar{n}$. Then there always exists an oriented incidence matrix $\Theta$ satisfying (12) such that the matrix $\hat{D}$ is block diagonal with blocks on the diagonal in $\mathbb{R}^{p \times p}$.

*Proof.* Let $G_a = (V_a, E_a)$ be a connected, simple, undirected graph with $V_a = \{1, ..., \bar{n}\}$. Let $\Theta^\top := \begin{bmatrix} \Theta_a^\top & 0 \\ Z & -I_{|V|-\bar{n}} \end{bmatrix}$, where $\Theta_a$ is an oriented incidence matrix associated with the graph $G_a$ satisfying (12) and $Z$ is a matrix in which each row contains exactly one entry equal to 1 and all other entries are zero. Clearly $\Theta \in \mathbb{R}^{|V| \times |E|}$ with $|E| = |E_a| + |V| - \bar{n}$, is an oriented incidence matrix associated with the graph $G$ (defined before (5)) and satisfies (12). Let us now check that $\hat{D}$ is a block diagonal matrix. Indeed, we have:

$$ \begin{aligned} \hat{D} = (\Theta^\top \otimes I_p) D(\Theta \otimes I_p) &= \begin{bmatrix} (\Theta_a^\top \otimes I_p) & 0 \\ (Z \otimes I_p) & -I_{p(|V|-\bar{n})} \end{bmatrix} \begin{bmatrix} 0 & 0 \\ 0 & \bar{D} \end{bmatrix} \begin{bmatrix} (\Theta_a \otimes I_p) & (Z^\top \otimes I_p) \\ 0 & -I_{p(|V|-\bar{n})} \end{bmatrix} \\ &= \begin{bmatrix} 0 & 0 \\ 0 & -\bar{D} \end{bmatrix} \begin{bmatrix} (\Theta_a \otimes I_p) & (Z^\top \otimes I_p) \\ 0 & -I_{p(|V|-\bar{n})} \end{bmatrix} = \begin{bmatrix} 0 & 0 \\ 0 & \bar{D} \end{bmatrix}, \end{aligned} $$

where $\bar{D} = \text{diag}(D_{\bar{n}+1}, ..., D_{|V|})$ and $0$ denotes zero matrices of appropriate dimensions. The proof is complete. $\square$

Note that the assumption “$D_i = 0, \forall i = \overline{1, \bar{n}}$” in Proposition 11 is made without loss of generality, in the sense that the agents (5) can always be labeled from the outset according to this rule. The proof of Proposition 11 provides a rule for building the network structure to ensure that $\hat{D}$ is a block diagonal matrix, stated as follows:

The agents (5) labeled from 1 to $\bar{n}$ can be connected in any manner, as long as they constitute a connected, simple, undirected graph. Meanwhile, each of the remaining agents is allowed to connect uniquely to one agent among these $\bar{n}$ agents.

### 3.3.3 Existence of the Origin Equilibrium of (21) under Assumption 6

The following proposition considers a particular class of heterogeneous networks (where all matrices $B_i, C_i$ in agents (5) are identical, and the vector fields differ only by a term matched with the input), and provides several sufficient conditions to guarantee the existence of an equilibrium point at the origin in the error dynamics corresponding to these network. Assumption 6 (dry friction condition) plays a central role in this proposition. It is worth noting that all the assumptions in this proposition are part
of the assumptions used in the main result (Theorem 3) presented in Section 4.1. A relationship between this proposition and Theorem 3 is shown in Remark 12.

**Proposition 13.** Suppose that all conditions in Lemma 1 hold, Assumptions 1, and 3' hold, and the following assertions hold:

1. There exist $\bar{B} \in \mathbb{R}^{n \times p}$, $\bar{C} \in \mathbb{R}^{p \times n}$ such that $B_i = \bar{B}, C_i = \bar{C}$, for all $i \in \mathcal{V}$.

2. There exists a positive definite matrix $\bar{P}' = \bar{P}'^\top \in \mathbb{R}^{n \times n}$ such that $\bar{P}'\bar{B} = \bar{C}^\top$ and $P' = (I_{|\mathcal{V}|} \otimes \bar{P}')$, where $P'$ is given in Assumption 3'.

3. There exist functions $f: \mathbb{R}_+ \times \mathbb{R}^n \to \mathbb{R}^n$ and $\xi_i: \mathbb{R}_+ \times \mathbb{R}^n \to \mathbb{R}^p$ such that:

3.1. $f_i(t, \cdot) = f(t, \cdot) + \bar{B}\xi_i(t, \cdot)$, $\forall t \ge 0$, $\forall i \in \mathcal{V}$.

3.2. There exists $\xi > 0$ such that: $\|\xi_i(t, \alpha) - \xi_j(t, \beta)\| \le \xi$, $\forall t \ge 0$, $\forall \alpha, \beta \in \mathbb{R}^n$, and $\forall i, j \in \mathcal{V}$.

Assume further that $f$ and $\xi_i$ are locally integrable in the first variable and Lipschitz continuous in the second variable for all $i = \overline{1, |\mathcal{V}|}$, and Assumption 6 holds with

$$ \rho > \frac{\xi \sqrt{|\mathcal{V}| - 1}}{\lambda_{\min}(\Theta_t^\top \Theta_t)}, \qquad (65) $$

for some oriented incidence matrix $\Theta_t$ corresponding to a spanning tree of the graph $\mathcal{G} = (\mathcal{V}, \mathcal{E})$. Then the error dynamics (21) has an equilibrium point at the origin.

*Proof.* Since items 1 and 3.1 hold, the dynamical system (37) is rewritten as:

$$ \begin{cases} \dot{x}(t) = \bar{F}(t, x(t)) + B\xi(t, x(t)) - B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)Cx(t) \\ \qquad + B\Xi^{-1}(\Theta \otimes I_p)H'\lambda'(t) \end{cases} \qquad (66a) $$

$$ w'(t) = N'\Delta^{-1}(\Theta^\top \otimes I_p)Cx(t) + \Gamma'\lambda'(t) + r' \qquad (66b) $$

$$ 0 \le \lambda'(t) \perp w'(t) \ge 0, \qquad (66c) $$

where $\bar{F}(t, x) := [f(t, x_1)^\top \dots f(t, x_{|\mathcal{V}|})^\top]^\top$, and $\xi(t, x) := [\xi_1(t, x_1)^\top \dots \xi_{|\mathcal{V}|}(t, x_{|\mathcal{V}|})^\top]^\top$. Since Assumption 3', items 1 and 2 hold, (66a) is equivalent to:

$$ P' \dot{x}(t) = P' \bar{F}(t, x(t)) - P' B \Xi^{-1}(\Theta \otimes I_p) W(\Theta^\top \otimes I_p) C x(t) + C^\top (\xi(t, x(t)) + (\Theta \otimes I_p)(N' \Delta^{-1})^\top \lambda'(t)) $$

where $P' \succ 0$ is given in Assumption 3'. Using item 2, it can be rewritten as:

$$ \begin{aligned} \dot{x}(t) = \bar{F}(t, x(t)) - B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)Cx(t) + B(\xi(t, x(t)) \\ + (\Theta \otimes I_p)(N'\Delta^{-1})^\top \lambda'(t)). \end{aligned} \qquad (67) $$

Note that $\lambda'(t) \in -\Phi'(N'\Delta^{-1}(\Theta^\top \otimes I_p)Cx(t))$ (see (41)). Since $f$ and $\xi_i$ are locally integrable in the first variable and Lipschitz continuous in the second variable for all
$i = \overline{1, |V|}$ (by assumption), so is $(f + \frac{1}{p|V|} \bar{B} \mathbf{1}_p \mathbf{1}_p^\top \sum_{i=1}^{|V|} \xi_i)$. Using [69, Theorem A1], the ODE:

$$ \dot{\bar{x}}(t) = f(t, \bar{x}) + \frac{1}{p|V|} \bar{B} \mathbf{1}_p \mathbf{1}_p^\top \sum_{i=1}^{|V|} \xi_i(t, \bar{x}) $$

has a unique locally absolutely continuous solution on $\mathbb{R}_+$ for any $\bar{x}(0) \in \mathbb{R}^n$. Let $\bar{x}(\cdot)$ be a solution of this ODE, and $x^*(t) := [\bar{x}(t)^\top \dots \bar{x}(t)^\top]^\top$ for all $t \ge 0$. Since item 1 holds, $Cx^* = [(\bar{C}\bar{x})^\top \dots (\bar{C}\bar{x})^\top]^\top$. Therefore, $(\Theta^\top \otimes I_p)Cx^* = \mathbf{0}_{p|\mathcal{E}|}$. This implies that $B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)Cx^* = \mathbf{0}_{n|V|}$, and $N'\Delta^{-1}(\Theta^\top \otimes I_p)Cx^* = \mathbf{0}_{m'|\mathcal{E}|}$. In what follows, we prove that $x^*(\cdot)$ is a solution of (66). To this end, it suffices to show the existence of $\lambda^*(t) \in -\mathbf{\Phi}'(N'\Delta^{-1}(\Theta^\top \otimes I_p)Cx^*(t)) \equiv -\mathbf{\Phi}'(0)$ such that $(x^*(t), \lambda^*(t))$ solves (67). By assumption, there is a spanning tree $\mathcal{G}_t = (\mathcal{V}, \mathcal{E}_t)$ of $\mathcal{G} = (\mathcal{V}, \mathcal{E})$ and an oriented incidence matrix $\Theta_t$ associated to $\mathcal{G}_t$ such that (65) holds. Note that $|\mathcal{E}_t| = |V| - 1$. In fact, since $\mathcal{G}_t$ is connected, $\mathrm{rank}(\Theta_t^\top \Theta_t) = \mathrm{rank}(\Theta_t) = |V| - 1$, using [73, Corollary 2.5.3]. Therefore, $\Theta_t^\top \Theta_t$ is a symmetric and positive definite matrix. Let

$$ \bar{h} := -(\Theta_t^\top \Theta_t \otimes I_p)^{-1}(\Theta_t^\top \otimes I_p)\xi(t, x^*) \in \mathbb{R}^{p|\mathcal{E}_t|}. $$

It follows from item 3.2 and (65) that $\|\bar{h}\| \le \rho$ with note that $\|(\Theta_t^\top \Theta_t \otimes I_p)^{-1}\| = (\lambda_{min}(\Theta_t^\top \Theta_t))^{-1}$. Let $\bar{h} := [h_1^\top \dots h_{|\mathcal{E}_t|}^\top]^\top$ with $h_i \in \mathbb{R}^p$ for all $i = 1, |\mathcal{E}_t|$. Let $h$ be the vector formed by the vectors $h_i$ for $i = 1, |\mathcal{E}_t|$ and $(|\mathcal{E}| - |\mathcal{E}_t|)$ vectors $\mathbf{0}_p$, such that $(\Theta \otimes I_p)h = (\Theta_t \otimes I_p)\bar{h}$. Clearly $h \in \mathcal{B}_\rho = \{\gamma \in \mathbb{R}^{p|\mathcal{E}|} : \|\gamma\| \le \rho\}$ since $\|\bar{h}\| \le \rho$. Since Assumption 6 holds, there exists $\lambda^* \in -\mathbf{\Phi}'(0)$ such that $(N'\Delta^{-1})^\top \lambda^* = h$. Therefore, the following hold:

$$ \begin{aligned} \xi(t, x^*) + (\Theta \otimes I_p)(N'\Delta^{-1})^\top \lambda^* &= \xi(t, x^*) + (\Theta \otimes I_p)h = \xi(t, x^*) + (\Theta_t \otimes I_p)\bar{h} \\ &= (I_{p|V|} - (\Theta_t \otimes I_p)(\Theta_t^\top \Theta_t \otimes I_p)^{-1}(\Theta_t^\top \otimes I_p)) \xi(t, x^*). \end{aligned} $$

It follows from [73, Fact 3.13.6, Eq (3.5.1)] that $I_{p|V|} - (\Theta_t \otimes I_p)(\Theta_t^\top \Theta_t \otimes I_p)^{-1}(\Theta_t^\top \otimes I_p)$ is the projection onto $\mathrm{ker}(\Theta_t^\top \otimes I_p)$. In addition, using [78, Theorem 8], the projection of $\xi(t, x^*)$ onto $\mathrm{ker}(\Theta_t^\top \otimes I_p) \equiv \mathrm{span}\{\mathbf{1}_{p|V|}\}$, is computed as:

$$ \begin{aligned} (I_{p|V|} - (\Theta_t \otimes I_p)(\Theta_t^\top \Theta_t \otimes I_p)^{-1}(\Theta_t^\top \otimes I_p)) \xi(t, x^*) &= \frac{\langle \xi(t, x^*), \mathbf{1}_{p|V|} \rangle}{\langle \mathbf{1}_{p|V|}, \mathbf{1}_{p|V|} \rangle} \mathbf{1}_{p|V|} \\ &= \frac{1}{p|V|} \begin{bmatrix} \mathbf{1}_p \mathbf{1}_p^\top \sum_{i=1}^{|V|} \xi_i(t, x^*) \\ \vdots \\ \mathbf{1}_p \mathbf{1}_p^\top \sum_{i=1}^{|V|} \xi_i(t, x^*) \end{bmatrix}. \end{aligned} $$

From all the above, $(x^*(t), \lambda^*(t))$ solves (67) (or equivalently, TVCS (37)). Since all condition in Lemma 1 and Assumption 1 hold, (37) is equivalent to (13) (see Section 3.2). Therefore, $x^*(\cdot)$ is a solution of (13). In addition, since $x^*(t) \in \mathrm{ker}(\Theta_c^\top \otimes I_n)$, the error dynamics (21) has an equilibrium point at the origin. The proof is complete. $\square$
# 4 Global synchronization of heterogeneous networks

In this section, we provide sufficient conditions to guarantee that the error dynamics (21) has a unique equilibrium point at the origin, which is globally exponentially (resp. finite-time) stable. As mentioned in Section 2, this guarantees the global exponential (resp. finite-time) synchronization of the heterogeneous network formed by agents (5), couplings (6) and the interconnection law (7). Throughout this section, the couplings (6) with a special structure given in Lemma 1 are considered. Let us consider the following lemma.

**Lemma 2.** Suppose that Assumptions 2 and 4 hold. Then the operator $\mathcal{P}$ defined in Assumption 6 is maximal monotone.

*Proof.* Since Assumption 2 holds, it follows from Remark 7 that $\Phi' = (\text{SOL}(0, \Gamma'))^* - r'$. In addition, it can be verified that $\text{rge}(N'\Delta^{-1}) \cap \text{rint}(\text{Dom} \Phi') \neq \emptyset$ due to Assumption 4. Therefore, the operator $\mathcal{P}(\cdot) = (\Delta^{-1})^\top N'\Phi'(N'\Delta^{-1}\cdot)$ is maximal monotone, using [54, Theorem 12.43]. $\square$

The next result follows directly from Lemma 2 and [49, Lemma 1]. It is needed for Theorem 3 below and is thus recalled.

**Lemma 3.** Let $\mathcal{P}$ be given in Assumption 6 and Assumptions 2, 4 and 6 hold. Then $\eta^\top \nu \ge \rho \|\eta\|$, for all $(\eta, \nu) \in \text{graph}(\mathcal{P})$.

## 4.1 Main results

The main result of this work is presented in the following theorem, where the global exponential stability of the equilibrium point at the origin of the error dynamics (21) is established. Later, finite-time stability will be proved. Recall that $\Phi'(0) = -\text{SOL}(r', \Gamma')$ (see proof of Proposition 7), $\Xi = I_{p|\nu|} + (\Theta \otimes I_p)W(\Theta^\top \otimes I_p)D$, $\Delta = I_{p|\mathcal{E}|} + \hat{D}W$, and $\bar{\Delta} = I_{p|\mathcal{E}|} + W\hat{D}$, where $\hat{D} = (\Theta^\top \otimes I_p)D(\Theta \otimes I_p)$.

**Theorem 3.** Suppose that all conditions in Lemma 1 and Corollary 3 hold, and Assumption 6 holds. Let the initial condition for the network be

$$x(0) \in (N'\Delta^{-1}(\Theta^\top \otimes I_p)C)^{-1}((\text{SOL}(0, \Gamma'))^* - r').$$

Assume further that the following assertions hold:

1. There exist $\bar{B} \in \mathbb{R}^{n \times p}$, $\bar{C} \in \mathbb{R}^{p \times n}$ such that $B_i = \bar{B}$, $C_i = \bar{C}$, for all $i \in \mathcal{V}$.

2. There exists a positive definite matrix $\bar{P}' = \bar{P}'^\top \in \mathbb{R}^{n \times n}$ such that:

2.1. $P' = (I_{|\nu|} \otimes \bar{P}')$, where $P'$ is given in Assumption 3.

2.2. $\bar{P}'\bar{B} = \bar{C}^\top$.

3. There exist functions $f: \mathbb{R}_+ \times \mathbb{R}^n \to \mathbb{R}^n$ and $\xi_i: \mathbb{R}_+ \times \mathbb{R}^n \to \mathbb{R}^p$ such that:
3.1. $f_i(t, \cdot) = f(t, \cdot) + \bar{B}\xi_i(t, \cdot), \forall t \ge 0, \forall i \in \mathcal{V}$.

3.2. There is a symmetric matrix $Q \in \mathbb{R}^{n \times n}$ such that $f$ is QUAD($\bar{P}', Q$), i.e.,
$$ (\gamma_1 - \gamma_2)^\top \bar{P}' (f(t, \gamma_1) - f(t, \gamma_2)) \le (\gamma_1 - \gamma_2)^\top Q (\gamma_1 - \gamma_2), \forall \gamma_1, \gamma_2 \in \mathbb{R}^n, \forall t \ge 0. $$

3.3. There exists $\xi > 0$ such that: $\|\xi_i(t, \alpha) - \xi_j(t, \beta)\| \le \xi, \forall t \ge 0, \forall \alpha, \beta \in \mathbb{R}^n$, and $\forall i, j \in \mathcal{V}$.

4. There exists $\mu > 0$ such that the gain $\gamma = \lambda_{min}(\bar{\Delta}^{-1}W + W^\top\bar{\Delta}^{-\top})$ satisfies:

$$ 2Q - (\gamma - c_0)\lambda_2(L)\bar{C}^\top\bar{C} + \mu\bar{P}' \ll 0. $$

5. The threshold value of $\rho$ given in Assumption 6 satisfies:

$$ \rho > \frac{\xi\sqrt{|\mathcal{V}| - 1}}{\sqrt{2\lambda_2(L)}} + c_1, $$

where $\lambda_2(L)$ is the second smallest eigenvalue of the Laplacian matrix $L = \Theta\Theta^\top$,

$$ c_0 = \frac{4\sqrt{h}}{\lambda_{min}^{\Gamma'}} \|N'\Delta^{-1}\| \|(\Delta^{-1})^\top(N')^\top - \bar{\Delta}^{-1}H'\|, $$

$$ c_1 = c \|(\Delta^{-1})^\top(N')^\top - \bar{\Delta}^{-1}H'\|, $$

with $c$ the norm of the least norm element of the set $\Phi'(0)$, $h$ and $\lambda_{min}^{\Gamma'}$ the rank and the smallest nonzero eigenvalue of $\Gamma' + (\Gamma')^\top$, respectively.

Then the origin of the error dynamics (42), which is equivalent to (21), is globally exponentially stable.

In fact, the state space of the network's dynamics (37) (or equivalently (13) when all conditions in Lemma 1 and Assumption 1 hold) is restricted to the set $S := (N'\Delta^{-1}(\Theta^\top \otimes I_p)C)^{-1}((\text{SOL}(0, \Gamma'))^* - r')$, which is positive invariant. Therefore, the term “global” in the global exponential stability of the origin is understood as global within the set $S$.

*Proof.* Since all conditions in Corollary 3 hold and the initial state $x(0)$ of the network satisfies $x(0) \in (N'\Delta^{-1}(\Theta^\top \otimes I_p)C)^{-1}((\text{SOL}(0, \Gamma'))^* - r')$, the TVCS (37) has a unique absolutely continuous solution $x(\cdot)$ on $\mathbb{R}_+$, using Corollary 3. In addition, it follows that Assumptions 1 through 4 hold. Since the couplings of the network have the special structure (28) (all conditions in Lemma 1 hold) and Assumption 1 holds, the TVCS (13) is equivalent to (37), see Section 3.2. Therefore, the error state (defined prior to (20)) can be rewritten as $e(t) := (\Theta_c^\top \otimes I_n)x(t)$ where $x(\cdot)$ is a solution of the system (37) (or equivalently, the differential inclusion (41)) instead of (13). This means that the dynamics of the error state $e(t)$ is defined by (42), which is equivalent to (21). For simplicity of notation, we will omit the dependence on $t$ for the variables $e(t)$ and $x(t)$. Let us consider the following candidate Lyapunov function:

$$ V(e) = \frac{1}{2} e^\top (I_{|\epsilon_c|} \otimes \bar{P}') e, $$
with $\bar{P}'$ defined in item 2. Note that since $\bar{P}' = \bar{P}'^{\top} \succ 0$, $(I_{\varepsilon_c} \otimes \bar{P}') \succ 0$, using [73, Fact 7.4.21]. The time derivative of $V$ along trajectories of (42) is given by:

$$ 
\begin{aligned}
\dot{V}(e) &= x^{\top}(\Theta_c \Theta_c^{\top} \otimes \bar{P}') (F(t, x) - B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)Cx - B\Xi^{-1}(\Theta \otimes I_p)H'\nu) \\
&= x^{\top}(\Theta_c \Theta_c^{\top} \otimes I_n)(P'F(t, x) - P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)Cx \\
&\qquad - P'B\Xi^{-1}(\Theta \otimes I_p)H'\nu),
\end{aligned}
\quad (68) $$

for some $\nu \in \Phi'(N'\Delta^{-1}(\Theta^{\top} \otimes I_p)Cx)$, and the second equality follows from item 2.1. Proceeding similarly to [72, Section 3.2] (or [69, Section 2.2.1]), (68) is rewritten as follows:

$$ 
\begin{aligned}
\dot{V}(e) &= x^{\top}(\Theta_c \Theta_c^{\top} \otimes I_n) \left[ P'F(t, x) - P'B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)Cx \right. \\
&\qquad \left. + (C^{\top}(\Theta \otimes I_p)(\Delta^{-1})^{\top}(N')^{\top} - P'B\Xi^{-1}(\Theta \otimes I_p)H') \lambda^{\text{im}}(x) \right. \\
&\qquad \left. - C^{\top}(\Theta \otimes I_p)(\Delta^{-1})^{\top}(N')^{\top} \nu \right] \\
&= W_1 - W_2 + W_3 - W_4,
\end{aligned}
\quad (69) $$

where $\lambda^{\text{im}}(\cdot) := \text{proj}_{\text{rge}(\Gamma'+\Gamma'\top)}(\Phi'(N'\Delta^{-1}(\Theta^{\top} \otimes I_p)C\cdot))$ is a Lipschitz continuous single-valued function [72, Lemma 3], and the terms $W_1, W_2, W_3, W_4$ are defined as follows:

$$ W_1 := x^{\top}(\Theta_c \Theta_c^{\top} \otimes \bar{P}')F(t, x). $$

$$ W_2 := x^{\top}(\Theta_c \Theta_c^{\top} \otimes \bar{P}')B\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)Cx. $$

$$ W_3 := x^{\top}(\Theta_c \Theta_c^{\top} \otimes I_n) (C^{\top}(\Theta \otimes I_p)(\Delta^{-1})^{\top}(N')^{\top} - P'B\Xi^{-1}(\Theta \otimes I_p)H') \lambda^{\text{im}}(x). $$

$$ W_4 := x^{\top}(\Theta_c \Theta_c^{\top} \otimes I_n)C^{\top}(\Theta \otimes I_p)(\Delta^{-1})^{\top}N'\nu. $$

Let us now evaluate $W_1, W_2, W_3, W_4$. First, using items 2.2 and 3, the following holds:

$$ 
\begin{aligned}
W_1 &= \sum_{1 \le i < j \le |V|} (x_i - x_j)^{\top} \bar{P}' (f_i(t, x_i) - f_j(t, x_j)) \\
&= \sum_{1 \le i < j \le |V|} (x_i - x_j)^{\top} \bar{P}' [ (f(t, x_i) - f(t, x_j)) + \bar{B}(\xi_i(t, x_i) - \xi_j(t, x_j))] \\
&\le \sum_{1 \le i < j \le |V|} (x_i - x_j)^{\top} Q (x_i - x_j) + \langle (\Theta_c^{\top} \otimes \bar{C})x, (\Theta_c^{\top} \otimes I_p)\xi(t, x) \rangle \\
&\le \sum_{1 \le i < j \le |V|} (x_i - x_j)^{\top} Q (x_i - x_j) + \| (I_{|\varepsilon_c|} \otimes \bar{C})e \| \| (\Theta_c^{\top} \otimes I_p)\xi(t, x) \| \\
&\le e^{\top} (I_{|\varepsilon_c|} \otimes Q) e + \xi \sqrt{|\varepsilon_c|} \| (I_{|\varepsilon_c|} \otimes \bar{C})e \|,
\end{aligned}
\quad (70) $$

where the second equality follows from item 3.1, and the first inequality follows from items 2.2, 3.2, and 3.3. Using items 1, 2.2, [84, Fact 3.20.6], the following is true:

$$ 
\begin{aligned}
W_2 &= x^{\top}(\Theta_c \Theta_c^{\top} \otimes \bar{P}'\bar{B})(\Theta \otimes I_p)\bar{\Delta}^{-1}W(\Theta^{\top} \otimes \bar{C})x \\
&= x^{\top}(\Theta_c \Theta_c^{\top} \Theta \otimes \bar{C}^{\top})\bar{\Delta}^{-1}W(\Theta^{\top} \otimes \bar{C})x = |V|x^{\top}(\Theta \otimes \bar{C}^{\top})\bar{\Delta}^{-1}W(\Theta^{\top} \otimes \bar{C})x. \\
&\ge \frac{1}{2}\gamma\|V\|\|(\Theta^{\top} \otimes \bar{C})x\|^2 \ge \frac{1}{2}\gamma\lambda_2(L)\|(I_{|\varepsilon_c|} \otimes \bar{C})e\|^2,
\end{aligned}
\quad (71) $$
where $\gamma = \lambda_{min}(\bar{\Delta}^{-1}W + W^{\top}\bar{\Delta}^{-\top})$, the third equality follows from Lemma B1 (see Appendix B), and the last inequality stems from:

$$ \begin{aligned} \|(\Theta^{\top} \otimes \bar{C})x\|^2 &= x^{\top}(\Theta\Theta^{\top} \otimes \bar{C}^{\top}\bar{C})x \\ &\geq \frac{\lambda_2(L)}{|\mathcal{V}|}x^{\top}(\Theta_c\Theta_c^{\top} \otimes \bar{C}^{\top}\bar{C})x \\ &= \frac{\lambda_2(L)}{|\mathcal{V}|}\|(I_{|\mathcal{E}_c|} \otimes \bar{C})e\|^2, \end{aligned} \quad (72) $$

where the first inequality follows from (ii4) in Remark 15 (see Appendix B) with $\lambda_2(L)$ the second smallest eigenvalue of the Laplacian matrix $L = \Theta\Theta^{\top}$. By assumption in item 4 (with $c_1$ existing), the matrix $\Gamma' + (\Gamma')^{\top}$ has a nonzero eigenvalue $\lambda_{min}^{\Gamma'}$. Therefore, $\Gamma' + (\Gamma')^{\top} \neq 0$. In addition, since $\Gamma' + (\Gamma')^{\top} \in \mathbb{R}^{m'|\mathcal{E}| \times m'|\mathcal{E}|}$ is a symmetric positive semidefinite matrix, there exists a matrix $V_{\Gamma'}$ such that $V_{\Gamma'}V_{\Gamma'}^{\top} = \Gamma' + (\Gamma')^{\top}$, $V_{\Gamma'}^{\top}V_{\Gamma'} \succ 0$, and $\|(V_{\Gamma'}^{\top}V_{\Gamma'})^{-1}\|_1 = (\lambda_{min}^{\Gamma'})^{-1} > 0$, where $\lambda_{min}^{\Gamma'}$ is the smallest eigenvalue of $\Gamma' + (\Gamma')^{\top}$. Indeed, using [73, Theorem 5.6.4], there exists a unitary matrix $S$ such that:

$$ \Gamma' + (\Gamma')^{\top} = S \begin{bmatrix} \Lambda & \mathbf{0}_{h \times (m'|\mathcal{E}|-h)} \\ \mathbf{0}_{(m'|\mathcal{E}|-h) \times h} & \mathbf{0}_{(m'|\mathcal{E}|-h) \times (m'|\mathcal{E}|-h)} \end{bmatrix} S^{\top} = S_1 \Lambda S_1^{\top} = V_{\Gamma'} V_{\Gamma'}^{\top}, $$

where $V_{\Gamma'} = S_1\Lambda^{\frac{1}{2}}$, $h = \text{rank}(\Gamma' + \Gamma^{\top})$, $\Lambda = \text{diag}(\lambda_1, ..., \lambda_h)$, with $\lambda_i$ the positive eigenvalues of $(\Gamma' + \Gamma^{\top})$ and $S = [S_1 S_2]$. Since $S$ is a unitary matrix, so $S_1^{\top}S_1 = I_h$ and then $V_{\Gamma'}^{\top}V_{\Gamma'} = \Lambda$. Therefore,

$$ \|(V_{\Gamma'}^{\top} V_{\Gamma'})^{-1}\| \leq \sqrt{h} \|(V_{\Gamma'}^{\top} V_{\Gamma'})^{-1}\|_1 = \sqrt{h} (\lambda_{min}^{\Gamma'})^{-1}, \quad (73) $$

using [73, Fact 9.8.12-(vii)]. In addition, since Assumption 6 holds, $\Phi'(0)$ is defined. Mimicking the proof of [72, Lemma 3], and using item 1, the following is obtained:

$$ \begin{aligned} \|\lambda^{\text{im'}}(x)\| &\leq \|\lambda^{\text{im'}}(x) - \lambda^{\text{im}}(0)\| + \|\lambda^{\text{im}}(0)\| \\ &\leq 2\|(V_{\Gamma'}^{\top}V_{\Gamma'})^{-1}\|\|N'\Delta^{-1}\|\|(\Theta^{\top} \otimes \bar{C})x\| + \|\lambda^{\text{im}}(0)\|, \\ &\leq \bar{c}_0\|(\Theta^{\top} \otimes \bar{C})x\| + \bar{c}_1, \end{aligned} \quad (74) $$

where the last inequality follows from (73) with $\bar{c}_0 := 2\sqrt{h}(\lambda_{min}^{\Gamma'})^{-1}\|N'\Delta^{-1}\|$, and $\bar{c}_1 := \|\lambda^{\text{im}}(0)\|$ the norm of the least norm element of the set $\Phi'(0)$, see [72, Lemma 3]. Using items 1, 2.2, and [84, Fact 3.20.6], the following holds:

$$ \begin{aligned} W_3 &= x^{\top}(\Theta_c\Theta_c^{\top}\Theta \otimes \bar{C}^{\top})\left((\Delta^{-1})^{\top}(N')^{\top} - \bar{\Delta}^{-1}H'\right)\lambda^{\text{im}}(x) \\ &= |\mathcal{V}|x^{\top}(\Theta \otimes \bar{C}^{\top})\left((\Delta^{-1})^{\top}(N')^{\top} - \bar{\Delta}^{-1}H'\right)\lambda^{\text{im}}(x) \\ &\leq |\mathcal{V}|\|(\Delta^{-1})^{\top}(N')^{\top} - \bar{\Delta}^{-1}H'\|\|(\Theta^{\top} \otimes \bar{C})x\|\|\lambda^{\text{im}}(x)\| \\ &\leq |\mathcal{V}|\|(\Delta^{-1})^{\top}(N')^{\top} - \bar{\Delta}^{-1}H'\|\|(\Theta^{\top} \otimes \bar{C})x\| (\bar{c}_0\|(\Theta^{\top} \otimes \bar{C})x\| + \bar{c}_1) \\ &= \bar{c}_0|\mathcal{V}|\|(\Delta^{-1})^{\top}(N')^{\top} - \bar{\Delta}^{-1}H'\|\|(\Theta^{\top} \otimes \bar{C})x\|^2 \\ &\quad + \bar{c}_1|\mathcal{V}|\|(\Delta^{-1})^{\top}(N')^{\top} - \bar{\Delta}^{-1}H'\|\|(\Theta^{\top} \otimes \bar{C})x\| \\ &\leq \bar{c}_0\lambda_2(L)\|(\Delta^{-1})^{\top}(N')^{\top} - \bar{\Delta}^{-1}H'\|\|(I_{|\mathcal{E}_c|} \otimes \bar{C})e\|^2 \\ &\quad + \bar{c}_1\sqrt{\lambda_2(L)|\mathcal{V}|}\|(\Delta^{-1})^{\top}(N')^{\top} - \bar{\Delta}^{-1}H'\|\|(I_{|\mathcal{E}_c|} \otimes \bar{C})e\|, \\ &= \frac{1}{2}\lambda_2(L)c_0\|(I_{|\mathcal{E}_c|} \otimes \bar{C})e\|^2 + c_1\sqrt{\lambda_2(L)|\mathcal{V}|}\|(I_{|\mathcal{E}_c|} \otimes \bar{C})e\|, \end{aligned} \quad (75) $$
where the second equality, the second inequality and the last inequality follow from Lemma B1, (74) and (72), respectively. Since item 1 holds, $\nu \in \Phi'(N'\Delta^{-1}(\Theta^\top \otimes \bar{C})x)$ and the following holds:

$$ 
\begin{aligned}
W_4 &= x^\top (\Theta_c \Theta_c^\top \Theta \otimes \bar{C}^\top) (\Delta^{-1})^\top N'^\top \nu = |\mathcal{V}| ((\Theta^\top \otimes \bar{C})x)^\top (\Delta^{-1})^\top N'^\top \nu \\
    &\geq \rho |\mathcal{V}| \|(\Theta^\top \otimes \bar{C})x\| \geq \rho \sqrt{|\mathcal{V}|\lambda_2(L)} \|(I_{|\mathcal{E}_c|} \otimes \bar{C})e\|, \quad (76)
\end{aligned}
$$

where the first equality follows from Lemma B1, the first inequality is inferred from Lemma 3 since Assumptions 2, 4, and 6 are satisfied, and the last inequality follows from (72). Since $|\mathcal{E}_c| = \frac{|\mathcal{V}|(|\mathcal{V}|-1)}{2}$, it follows from all the above and items 4, 5 that:

$$ 
\begin{aligned}
\dot{V}(e) &\leq \frac{1}{2} e^\top \left( (I_{|\mathcal{E}_c|} \otimes 2Q) - (\gamma - c_0)\lambda_2(L) (I_{|\mathcal{E}_c|} \otimes \bar{C}^\top \bar{C}) \right) e \\
&\quad + \left( \xi \sqrt{|\mathcal{E}_c|} + c_1 \sqrt{\lambda_2(L)|\mathcal{V}|} - \rho \sqrt{|\mathcal{V}|\lambda_2(L)} \right) \|(I_{|\mathcal{E}_c|} \otimes \bar{C})e\| \leq -\mu V(e). \quad (77)
\end{aligned}
$$

Thus the origin of the error dynamics system (42) (or equivalent to (21)) is globally exponentially stable. $\square$

The following corollary shows that synchronization occurs in finite time.

**Corollary 6 (Full-state coupling).** Suppose that all conditions in Theorem 3 hold, and $\bar{C} = I_n$. Then the network achieves global synchronization in finite time with the upperbound:

$$ t^* := \frac{2\sqrt{\lambda_{max}(\bar{P}')V(e(0))}}{(\rho - c_1)\sqrt{2|\mathcal{V}|\lambda_2(L)} - \xi\sqrt{|\mathcal{V}|(|\mathcal{V}| - 1)}. $$

*Proof.* Proceeding similarly to the proof of Theorem 3, the inequality (77) can be rewritten as:

$$ 
\begin{aligned}
\dot{V}(e) &\leq \frac{1}{2} e^{\top} \left( (I_{|\mathcal{E}_c|} \otimes 2Q) - (\gamma - c_0)\lambda_2(L) I_{n|\mathcal{E}_c|} \right) e + \left( \xi \sqrt{|\mathcal{E}_c|} + (c_1 - \rho) \sqrt{|\mathcal{V}|\lambda_2(L)} \right) \|e\| \\
&\leq - \left( (\rho - c_1) \sqrt{|\mathcal{V}|\lambda_2(L)} - \xi \sqrt{|\mathcal{E}_c|} \right) \sqrt{\frac{2}{\lambda_{max}(\bar{P}')}} \sqrt{V(e)}, \quad (78)
\end{aligned}
$$

where the last inequality follows from $V(e) = \frac{1}{2}e^\top(I_{|\mathcal{E}_c|} \otimes \bar{P}')e \leq \frac{1}{2}\lambda_{max}(\bar{P}')\|e\|^2$. Since $|\mathcal{E}_c| = \frac{|\mathcal{V}|(|\mathcal{V}|-1)}{2}$, the result follows from direct integration of (78). $\square$

When $f_i(t, x_i)$ in (5) reduces to $f_i(x_i)$ for all $i \in \mathcal{V}$, then the error dynamics (42), (21) are autonomous systems. In that case, some conditions in Theorem 3 could be relaxed by using the Krasovskii-LaSalle invariance principle.

**Remark 11.** If $(\Gamma' + \Gamma'^\top) = 0$, or $\bar{\Delta}^{-1}H' = (N'\Delta^{-1})^\top$, or more generally Assumption 3' holds, the term $W_3$ in the proof of Theorem 3 vanishes. Then items 4 and 5 in Theorem 3 and the upperbound in Corollary 6 are rewritten as follows:

[Item 4.] There exists $\mu > 0$ such that the gain $\gamma = \lambda_{min}(\bar{\Delta}^{-1}W + W^\top\bar{\Delta}^{-\top})$ satisfies:

$$ 2Q - \gamma\lambda_2(L)\bar{C}^\top\bar{C} + \mu\bar{P}' \leq 0. $$
[Item 5.] The threshold value of $\rho$ given in Assumption 6 satisfies:

$$ \rho > \frac{\xi \sqrt{|\mathcal{V}|-1}}{\sqrt{2\lambda_2(L)}}. $$

$$ \text{[Upperbound]. } t^* := \frac{2\sqrt{\lambda_{\max}(\bar{P}')V(e(0))}}{\rho\sqrt{2|\mathcal{V}|\lambda_2(L)} - \xi\sqrt{|\mathcal{V}|(|\mathcal{V}|-1)}}. $$

**Remark 12.** It is worth noting that Proposition 13 in Section 3.3 guarantees the existence of an equilibrium point at the origin of the error dynamics (21) under a certain set of assumptions. Theorem 3 builds upon this result by employing a largely similar set of assumptions with only minor differences in a few items (e.g., threshold for the value $\rho$ defined in Assumption 6) and introducing some additional conditions (e.g., items 3.2 and 4). These enhancements allow Theorem 3 to not only ensure existence, but also establish the global exponential stability of the equilibrium point at the origin of (21).

## 4.2 Coupling Design Algorithm for Network Synchronization

In this subsection, items 1, 2, 3 in Theorem 3 are assumed to hold, and $Q, \bar{C}, D, \bar{P}', \xi$ (defined in (9) and Theorem 3, characterizing the agents (5)) and $\Theta, |\mathcal{V}|$ (characterizing the network structure) are assumed to be known. The design of the LCP couplings (6) to achieve global exponential (resp. finite-time) synchronization of the agents (5) using the interconnection law (7), is generally highly challenging. This subsection summarizes the main results presented in the previous sections into an algorithm that facilitates a more straightforward and efficient design of such couplings.

Recall that the network structure is represented by a connected, simple, undirected graph $\mathcal{G} = (\mathcal{V}, \mathcal{E})$, and is characterized by an oriented incidence matrix $\Theta$ associated with $\mathcal{G}$. The matrices $H, M, N, r$ (defined after (10)) characterize the LCP couplings (6). Due to the special structure of $H, M, N, r$ proposed in Lemma 1 and consistently used throughout this work, the LCP couplings (6) are now characterized by the matrices $W, H', M', N', r'$ defined after (34). It can be seen that Assumptions 1, 2, 3, 4, 6, and items 4, 5 in Theorem 3 are conditions related to the LCP couplings (6), play roles in ensuring the global exponential (resp. finite-time) synchronization of the network, and they are difficult to verify. The following algorithm provides the steps to compute the matrices $W, H', M', N', r'$ to ensure that such conditions are satisfied. A brief proof is given later. Recall that $\hat{D} = (\Theta^{\top} \otimes I_p)D(\Theta \otimes I_p) \in \mathbb{R}^{p|\mathcal{E}| \times p|\mathcal{E}|}$, $\bar{\Delta} = I_{p|\mathcal{E}|} + W\hat{D}$, and $\Delta = I_{p|\mathcal{E}|} + \hat{D}W$.
**Algorithm 1** Calculation of W, N', H', M', r' to satisfy Assumptions 1, 2, 3, 4, 6, and items 4, 5 in Theorem 3.

**Require:** The matrices $\Theta, Q, \bar{C}, D, \bar{P}',$ and the constants $|\mathcal{V}|, \xi$.

1: Check whether the matrix $\hat{D}$ is block diagonal with each block on the diagonal being a matrix in $\mathbb{R}^{p \times p}$. If yes, continue; otherwise, stop.

2: Compute $W$ satisfying the conditions in item 1 in Proposition 11 such that item 4 in Theorem 3 holds with $c_0 = 0$. If such $W$ exists, continue; otherwise, stop.

3: Compute $N'$ such that the conditions in item 2 in Proposition 11 hold. If such $N'$ exists, continue; otherwise, stop.

4: Compute $H' = \bar{\Delta}(\Delta^{-1})^{\top}(N')^{\top}$.

5: Compute $M'$ such that the condition in item 2 in Proposition 10 holds and $\text{rgc}(N'\Delta^{-1}(\Theta^{\top} \otimes I_p)C) \cap \text{rint}((\text{SOL}(0, \Gamma'))^*) \neq \emptyset$. If such $M'$ exists, continue; otherwise, stop.

6: Compute $r'$ such that the condition in item 3 in Proposition 10, and item 5 in Theorem 3 hold with $c_1 = 0$. If such $r'$ exists, the algorithm terminates successfully.

The design of $H'$ in step 4 in Algorithm 1 is explained as follows. In Theorem 3, the conditions on $H'$ are implicitly encoded in Assumptions 2, 3, 4, 6 through the matrix $\Gamma' = N'\Delta^{-1}\hat{D}H' + M'$, and in items 4, 5 through the constants $c_0, c_1$. However, if steps 1, 2, 3 in the Algorithm 1 are successfully executed, then $H'$ no longer plays a role in $\Gamma'$ since there always exists a matrix $M'$ that makes $\Gamma'$ satisfy the desired properties, such as item 2 in Proposition 11. Therefore, the conditions on $H'$ are solely reflected through items 4, 5 and Assumption 3, which is recalled as:

$$ \ker (\Gamma' + (\Gamma')^{\top}) \subseteq \ker (P'B\Xi^{-1}(\Theta \otimes I_p)H' - C^{\top}(\Theta \otimes I_p)(\Delta^{-1})^{\top}N'^{\top}), $$

for some $P' = P'^{\top} \succ 0$. Since items 1, 2 in Theorem 3 are assumed to hold, and $\Xi^{-1}(\Theta \otimes I_p) = (\Theta \otimes I_p)\bar{\Delta}^{-1}$ using [84, Fact 3.20.6], Assumption 3 is rewritten as:

$$ \ker (\Gamma' + (\Gamma')^{\top}) \subseteq \ker (C^{\top}(\Theta \otimes I_p)(\bar{\Delta}^{-1}H' - (\Delta^{-1})^{\top}N'^{\top})). \quad (79) $$

It can be seen that if $\bar{\Delta}^{-1}H' = (\Delta^{-1})^{\top}N'^{\top}$, then (79) and item 4 hold (since step 2 succeeded), and simultaneously ensures that step 6 always succeeds. This provides us with the option of designing $H' = \bar{\Delta}(\Delta^{-1})^{\top}N'^{\top}$. Note that when the matrix $DW$ is symmetric, $H' = N'^{\top}$ since $\bar{\Delta} = \Delta = I_{p|\mathcal{E}|}$.

The brief proof showing that the successful execution of Algorithm 1 implies that Assumptions 1, 2, 3, 4, 6, and items 4, 5 in Theorem 3 are satisfied, is presented as follows:

1. Since $W$ in step 2 exists, Assumption 1 holds (see item 1 in Proposition 11).

2. Since $N', H'$ in steps 3, 4 exist, Assumption 3 holds, and $c_0 = c_1 = 0$ where $c_0, c_1$ are defined in Theorem 3.
3. Since $N', M'$, and $r'$ in steps 3, 5, 6 exist, respectively, then Assumptions 2, 6, 4 hold (see Propositions 10, 11, and item 1.1 in Proposition 7), and items 4, 5 in Theorem 3 hold since $c_0 = c_1 = 0$.

Algorithm 1 is used to design the LCP couplings (6) for all examples in the next section.

# 5 Examples and Numerical Simulations

As mentioned in Section 1, most previous studies on the synchronization of heterogeneous networks, that is, networks in which the agents (5) possess pairwise distinct vector fields $f_i$ (e.g., [38], [36], [41]), have achieved practical synchronization. In such cases, the error state $e(t)$ (defined in (20)) remains bounded but does not converge to zero, in contrast to exponential synchronization. The main objective of this section is to provide examples that illustrate the theoretical results presented in Section 4.1, where a heterogeneous network is shown to reach global exponential (or finite-time) synchronization. In this section, the terms $\bar{B}\xi_i$ appearing in item 3.1 in Theorem 3 are interpreted as quantifiable representations of the mismatches between the vector fields of agents (5). With precise knowledge of $\bar{B}\xi_i$, one can design control inputs $u_i(t)$ (defined in (5)) that theoretically guarantee global exponential (or finite-time) synchronization of the network. When considering practical scenarios where $\bar{B}\xi_i$ are unknown perturbations, but the bounds on the differences between $\xi_i$ and $\xi_j$ for all $i, j = 1, |V|$ can be measured (i.e., item 3.3 in Theorem 3 is satisfied), the design of a discrete-time controller $u_i(t)$ that still ensures network synchronization is referred to as a practical implementation and it is investigated in detail in Section 6.

## 5.1 Introduction to Numerical Schemes

In each example of this section, numerical simulations of the closed-loop system described by the TVCS (37) are carried out to show the exponential (or finite-time) synchronization of the states of agents (5) once interconnected via LCP couplings (6) and the interconnection law (7). It is useful to illustrate the behavior of the network's continuous-time dynamics. The discretization and simulation of the TVCSs (or of sliding-mode systems) require specific numerical schemes, which are well established in Contact Mechanics (see [86], [87] for details). Common schemes employed for discretizing complementarity systems include the Forward (fully-explicit) Euler scheme, the Backward (fully-implicit) Euler scheme, and the Semi-Backward (semi-implicit) Euler scheme. The selection of an appropriate scheme for simulation must be based on its main features, which are briefly outlined as follows:

1. The Forward Euler scheme is simple and computationally inexpensive, making it suitable for non-stiff systems. However, it has low accuracy and poor stability when applied to stiff systems, often requiring very small time steps to ensure reliable results. It is also well-known that such schemes yield numerical chattering effects that should be avoided.
2. The Backward Euler scheme provides strong numerical stability even with large time steps, making it well-suited for stiff systems. Concerning nonsmooth systems with sliding modes, such methods are known to correctly handle sliding motions without any numerical chattering. This is important both for numerical simulations, but also for practical implementation of set-valued controllers. However, at each step, it requires solving a nonlinear equation, increasing computational cost and complexity.

3. The Semi-Backward Euler scheme (which in our case consists of discretizing implicitly the set-valued terms only) combines the advantages of both Forward and Backward Euler schemes. It evaluates some terms using current values and others using future values, improving stability over Forward Euler while being computationally less demanding than Backward Euler.

Based on the above analysis, in this section, we discretize the TVCS (37) using the Semi-Backward Euler scheme. Specifically, let $[0, T]$, $T > 0$, be the interval that we need to obtain the solution, $(t_0, t_1, ..., t_l)$ be a partition of $[0, T]$ with $t_{k+1} - t_k = h > 0$ the time step, $l = \frac{T}{h}$, and $x_k := x(t_k)$, $\lambda'_k := \lambda'(t_k)$, and $w'_k := w'(t_k)$, for all $k = 0, \overline{l}$. Mimicking ideas for discretizing LCS [87], a time-discretization of TVCS (37) using Semi-Backward Euler scheme, is given by:

$$ \left\{ \begin{array}{l} x_{k+1} = x_k + hF(t_k, x_k) - hB\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)Cx_k + hB\Xi^{-1}(\Theta \otimes I_p)H'\lambda'_{k+1}, \\ 0 \le \lambda'_{k+1} \perp w'_{k+1} = N'\Delta^{-1}(\Theta^\top \otimes I_p)Cx_{k+1} + \Gamma'\lambda'_{k+1} + r' \ge 0. \end{array} \right. \qquad (80) $$

After some algebraic manipulations, (80) can be rewritten as:

$$ \left\{ \begin{array}{l} x_{k+1} = x_k + hF(t_k, x_k) - hB\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)Cx_k + hB\Xi^{-1}(\Theta \otimes I_p)H'\lambda'_{k+1}, \\ 0 \le \lambda'_{k+1} \perp \bar{\Gamma}'\lambda'_{k+1} + q_k \ge 0, \end{array} \right. \qquad (81) $$

where $\bar{\Gamma}' := \Gamma' + hN'\Delta^{-1}(\Theta^\top \otimes I_p)CB\Xi^{-1}(\Theta \otimes I_p)H'$, and

$$ q_k := N'\Delta^{-1}(\Theta^\top \otimes I_p)C(x_k + hF(t_k, x_k) - hB\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)Cx_k) + r', \forall k = 0, \overline{l}. $$

Therefore, given the initial state $x(0)$, at each step $k + 1$ of the simulation process, the problem LCP($q_k$, $\bar{\Gamma}'$) is solved to compute $\lambda'_{k+1}$, which is then used to update the state $x_{k+1}$.

## 5.2 First Example

The following example illustrates Corollary 6 in Section 4.1. In addition, it shows that synchronization of agents (5) can be achieved even when $D_i \neq 0$ or $D_i \ll 0$, for $i \in \mathcal{V}$. This aspect, emphasized as an original contribution, see Remark 2. The important role of Assumption 6 in Theorem 3 and Corollary 6 is also discussed in this example. The LCP couplings in this example are designed using Algorithm 1.
**Example 2.** Let us consider the following three systems:

$$ \left\{ \begin{array}{l} \dot{x}_i = f(t, x_i) + \xi_i(t, x_i) + u_i, \quad \forall i = 1, 2, 3, \\ y_i = x_i + d_i u_i \end{array} \right. \qquad (82) $$

where $f(t, x_i) = f(t)x_i$ with $f(t) = a_1 \cos(a_2t + a_3)$, and the mismatches between agents are given as follows:

$$ \left\{ \begin{array}{l} \xi_1(t, x_1) = \gamma_{1,1} \sin(2tx_1 + 1) \cos(0.5t - x_1) + \gamma_{1,2} \cos(14tx_1 + 1) \sin(-4t), \\ \xi_2(t, x_2) = \gamma_{2,1} \sin(20t + x_2) \cos(3t - 1) + \gamma_{2,2} \cos(0.2t + 1) \sin(-tx_2), \\ \xi_3(t, x_3) = \gamma_{3,1} \sin(12tx_3 - 2) \cos(2t + 4) + \gamma_{3,2} \cos(4tx_3) \sin(-4t - 1), \end{array} \right. $$

where $a_1 \in [-0.4, 0.4]$, $a_2, a_3 \in \mathbb{R}$, $d_i \in \mathbb{R}$, and $\gamma_{i,1}, \gamma_{i,2} \in [-0.2, 0.2]$, for all $i = 1, 2, 3$. In this example, we consider $d_1 = 0, d_2 = d_3 = 1$, and a network with the graph in Figure 4.

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

$$ \Theta = \begin{bmatrix} 1 & 1 \\ -1 & 0 \\ 0 & -1 \end{bmatrix} $$

<caption>Figure 4: Graph of the network in Example 2, and its oriented incidence matrix $\Theta$.</caption>

Let us consider two couplings as in (28), where for all $i, j \in \mathcal{V}$, $i < j$, and $(i, j) \in \mathcal{E}$, we have:

$$ W_{ij} = 1, \quad H'_{ij} = \begin{bmatrix} 2 & -1 \end{bmatrix}, \quad N'_{ij} = \begin{bmatrix} 2 \\ -1 \end{bmatrix}, \quad M'_{ij} = \begin{bmatrix} -1 & 0 \\ 0 & 0.5 \end{bmatrix}, \quad \text{and} \quad r'_{ij} = \begin{bmatrix} 2 \\ -2 \end{bmatrix}. \quad (83) $$

It follows from (37) that the dynamics of the network is defined as the TVCS:

$$ \left\{ \begin{array}{l} \dot{x}(t) = \begin{bmatrix} f_1(t, x_1) \\ f_2(t, x_2) \\ f_3(t, x_3) \end{bmatrix} - \begin{bmatrix} 1 & -0.5 & -0.5 \\ -0.5 & 0.5 & 0 \\ -0.5 & 0 & 0.5 \end{bmatrix} x(t) + \begin{bmatrix} 1 & -0.5 & 1 & -0.5 \\ -1 & 0.5 & 0 & 0 \\ 0 & 0 & -1 & 0.5 \end{bmatrix} \lambda'(t), \\ \\ w'(t) = \begin{bmatrix} 1 & -1 & 0 \\ -0.5 & 0.5 & 0 \\ 1 & 0 & -1 \\ -0.5 & 0 & 0.5 \end{bmatrix} x(t) + \begin{bmatrix} 1 & -1 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 0 & 0 & 1 & -1 \\ 0 & 0 & -1 & 1 \end{bmatrix} \lambda'(t) + \begin{bmatrix} 2 \\ -2 \\ 2 \\ -2 \end{bmatrix}, \\ \\ 0 \le \lambda'(t) \perp w'(t) \ge 0, \end{array} \right. \qquad (84) $$
where $f_i(t, x_i) = f(t)x_i + \xi_i(t, x_i)$, for all $i = 1, 2, 3$, and $\lambda'(t), w'(t) \in \mathbb{R}^4$ are the complementarity variables defined by:

$\lambda'(t) = [\lambda'_{12,1}(t) \ \lambda'_{12,2}(t) \ \lambda'_{13,1}(t) \ \lambda'_{13,2}(t)]^\top$; $w'(t) = [w'_{12,1}(t) \ w'_{12,2}(t) \ w'_{13,1}(t) \ w'_{13,2}(t)]^\top$.

Set $\lambda'_{ij}(t) = [\lambda'_{ij,1}(t) \ \lambda'_{ij,2}(t)]^\top$, and $w'_{ij}(t) = [w'_{ij,1}(t) \ w'_{ij,2}(t)]^\top$, for $(i, j) \in \{(1, 2), (1, 3)\}$.

**The structure of the couplings:** Recall from (7) that $v_{12} = y_1 - y_2 = (x_1 - x_2) + (d_1u_1 - d_2u_2)$. Using Lemma 1, the couplings (6) in this example are given by:

$$ \left\{ \begin{array}{l} \eta_{12}(t) = [2 - 1]\lambda'_{12}(t) - v_{12}(t), \\ 0 \le \lambda'_{12}(t) \perp w'_{12}(t) = \begin{bmatrix} -1 & 0 \\ 0 & 0.5 \end{bmatrix} \lambda'_{12}(t) + \begin{bmatrix} 2 \\ -1 \end{bmatrix} v_{12}(t) + \begin{bmatrix} 2 \\ -2 \end{bmatrix} \ge 0, \end{array} \right. \qquad (85) $$

and $\eta_{13}(t)$ has the same structure as $\eta_{12}(t)$. Let $\eta'_{12}(t) = [2 - 1]\lambda'_{12}(t)$, then $\eta'_{12}$ is a set-valued mapping in the variable $v_{12}$, defined by:

$$ \eta'_{12} : [-1, +\infty) \Rightarrow \mathbb{R}, \quad v_{12} \mapsto \{2v_{12}, -2(v_{12} + 2)\}. $$

Its graph is given in Figure 5a, which depicts the relationship between the input $v_{12}(t)$ and the output $\eta_{12}(t)$ of the coupling. On the other hand, let $\bar{v}_{12} = \frac{1}{2}(x_1 - x_2)$ where $x_1, x_2$ are the states of the agents 1, 2 in (82). It follows from (85) and (84) that the following holds:

$$ \left\{ \begin{array}{l} \eta'_{12}(t) = [2 - 1]\lambda'_{12}(t), \\ 0 \le \lambda'_{12}(t) \perp w'_{12}(t) = \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \lambda'_{12}(t) + \begin{bmatrix} 2 \\ -1 \end{bmatrix} \bar{v}_{12}(t) + \begin{bmatrix} 2 \\ -2 \end{bmatrix} \ge 0. \end{array} \right. \qquad (86) $$

Therefore, it can be checked that $\eta'_{12} : [0, +\infty) \Rightarrow \mathbb{R}$, $\bar{v}_{12} \mapsto \begin{cases} [-2, +\infty), & \text{if } \bar{v}_{12} = 0 \\ -\bar{v}_{12} - 2, & \text{if } \bar{v}_{12} > 0 \end{cases}$.

Since the input $u$ defined in (11) and (36) is identical, it follows that:

$$ \eta_{12} = \frac{1}{2}\eta'_{12} - \bar{v}_{12} = \begin{cases} [-1, +\infty), & \text{if } \bar{v}_{12} = 0 \\ -\frac{3}{2}\bar{v}_{12} - 1, & \text{if } \bar{v}_{12} > 0 \end{cases}. $$

This implies that for every $\bar{v}_{12}(t) \in \mathbb{R}_+$ (which relates only to the states $x_1, x_2$), the set of possible output values of the coupling between agents 1, 2 is always closed and convex. The same holds for the coupling between agents 1, 3. In practical implementation, the closedness and convexity of the set of output values of the couplings ensures the existence and uniqueness of the input of the agents $u(t)$ corresponding to a given appropriate state $x(t)$ at any time $t = t_k$, where $u(t)$ and $x(t)$ are defined in (9). The graphs of mappings $\eta'_{12}(\bar{v}_{12})$, $\eta_{12}(\bar{v}_{12})$ are depicted in Figure 5b.
<figure><img src="image_6.png" /><figcaption>Figure 5: Input-output relation at the coupling between agents 1 and 2 in Example 2.</figcaption></figure>

**Existence and uniqueness of the equilibrium point:** In this part, let us consider the nominal system of (84) where $\xi_i$ is not taken into account. Let $x^*(t) = \mathbf{0}_3$, and $\lambda^*(t) = \mathbf{0}_4$. It can be checked that $(x^*(t), \lambda^*(t))$ solves this nominal system, hence the point $e^* = \mathbf{0}_3$ is an equilibrium point of the error dynamics corresponding to the nominal system, using Definition 11. Moreover, it is the unique equilibrium point. Indeed, suppose on the contrary that there exists a nonzero equilibrium point $e^*$. It follows from Remark 5 that there exist $x^*(t) = [x_1^*(t) \ x_2^*(t) \ x_3^*(t)]^\top \in \mathbb{R}^3$, $\alpha: \mathbb{R}_+ \to \mathbb{R}$, and $p_2, p_3 \in \mathbb{R}$ such that:

$$ \left\{ \begin{array}{l} x^*(t) \text{ is a solution of the nominal system of (84),} \\ x_1^*(t) = \alpha(t), \ x_2^*(t) = \alpha(t) + p_2, \ x_3^*(t) = \alpha(t) + p_3, \\ \mathbf{0}_3 \neq e^* = (\Theta_c^\top \otimes I_n)x^*(t). \end{array} \right. \qquad (87) $$

Since $\dot{x}_1^*(t) = \dot{x}_2^*(t) = \dot{x}_3^*(t)$, it follows from (87) and the nominal system of (84) that for all $t \ge 0$:

$$ f(t)p_2 - \frac{1}{2}p_2 - \frac{1}{2}(p_2 + p_3) - 2\lambda'_{12,1}(t) + \lambda'_{12,2}(t) - \lambda'_{13,1}(t) + \frac{1}{2}\lambda'_{13,2}(t) = 0, \quad (88a) $$

$$ f(t)p_3 - \frac{1}{2}p_3 - \frac{1}{2}(p_2 + p_3) - \lambda'_{12,1}(t) + \frac{1}{2}\lambda'_{12,2}(t) - 2\lambda'_{13,1}(t) + \lambda'_{13,2}(t) = 0, \quad (88b) $$

$$ f(t)(p_2 - p_3) - \frac{1}{2}(p_2 - p_3) - \lambda'_{12,1}(t) + \frac{1}{2}\lambda'_{12,2}(t) + \lambda'_{13,1}(t) - \frac{1}{2}\lambda'_{13,2}(t) = 0, \quad (88c) $$
where $\lambda'_{12,1}(t), \lambda'_{12,2}(t), \lambda'_{13,1}(t)$ and $\lambda'_{13,2}(t)$ satisfies the following LCPs:

$$ \left\{ \begin{array}{l} 0 \le \begin{bmatrix} \lambda'_{12,1}(t) \\ \lambda'_{12,2}(t) \end{bmatrix} \perp \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \begin{bmatrix} \lambda'_{12,1}(t) \\ \lambda'_{12,2}(t) \end{bmatrix} + \begin{bmatrix} -1 \\ \frac{1}{2} \end{bmatrix} p_2 + \begin{bmatrix} 2 \\ -2 \end{bmatrix} \ge 0, \\ 0 \le \begin{bmatrix} \lambda'_{13,1}(t) \\ \lambda'_{13,2}(t) \end{bmatrix} \perp \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \begin{bmatrix} \lambda'_{13,1}(t) \\ \lambda'_{13,2}(t) \end{bmatrix} + \begin{bmatrix} -1 \\ \frac{1}{2} \end{bmatrix} p_3 + \begin{bmatrix} 2 \\ -2 \end{bmatrix} \ge 0. \end{array} \right. $$

Proceeding similarly to the computation of $\eta'_{12}$ after (86), then the following hold:

$$ 2\lambda'_{12,1}(t) - \lambda'_{12,2}(t) = \begin{cases} (-\infty, 2], & \text{if } p_2 = 0 \\ \frac{1}{2}p_2 - 2, & \text{if } p_2 < 0 \end{cases}, $$

$$ 2\lambda'_{13,1}(t) - \lambda'_{13,2}(t) = \begin{cases} (-\infty, 2], & \text{if } p_3 = 0 \\ \frac{1}{2}p_3 - 2, & \text{if } p_3 < 0 \end{cases}. $$

The following reasoning shows that $p_2 = p_3 = 0$, which contradicts (87).

1. If $p_2 < 0, p_3 < 0$, and $p_2 \neq p_3$, it follows from (88c) that for all $t \ge 0$, $(f(t) - \frac{3}{4})(p_2 - p_3) = 0$. This leads to a contradiction.

2. If $p_2 < 0, p_3 < 0$, and $p_2 = p_3$, it follows from (88a) that for all $t \ge 0$, $f(t)p_2 - \frac{9}{4}p_2 + 3 = 0$. This leads to a contradiction.

3. If $p_2 = 0, p_3 < 0$ (resp. $p_2 < 0, p_3 = 0$), it follows from (88a) (resp. (88b)), (88c) that the following holds for all $t \ge 0$:

$$ \left( \text{resp.} \begin{cases} -\frac{3}{4}p_3 + 1 = 2\lambda'_{12,1} - \lambda'_{12,2}, \\ -p_3f(t) + \frac{3}{4}p_3 - 1 = \lambda'_{12,1} - \frac{1}{2}\lambda'_{12,2}. \end{cases} \right) \begin{pmatrix} -\frac{3}{4}p_2 + 1 = 2\lambda'_{13,1} - \lambda'_{13,2} \\ p_2f(t) - \frac{3}{4}p_2 + 1 = -\lambda'_{13,1} + \frac{1}{2}\lambda'_{13,2}. \end{pmatrix} $$

This leads to a contradiction.

The existence and uniqueness of the equilibrium point at the origin of the error dynamics corresponding to the heterogeneous network dynamics (84) (i.e., $\xi_i$ is taken into account) are implicitly established by Theorem 3, whose assumptions will be verified in the next part.

**Global synchronization:** Let us check the validity of all assumptions in Theorem 3.

1. From the network's structure, all conditions in Lemma 1, and items 1, 3.1 hold.

2. Since $\Delta = \begin{bmatrix} 2 & 0 \\ 0 & 2 \end{bmatrix}$ and $\Gamma' = \text{diag}(\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}, \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}) \ge 0$, Assumptions 1 and 2 hold.
3. It can be checked that Assumption 3' holds with $P' = I_3$. Therefore, Assumption 3 holds, and item 2 holds with $\bar{P}' = 1$.

4. Using (40) in Remark 7, we have

$$(\mathrm{SOL}(0, \Gamma'))^* = \{[q_1 \ q_2 \ q_3 \ q_4]^\top \in \mathbb{R}^4 : q_1 + q_2 \ge 0, q_3 + q_4 \ge 0\}.$$

Moreover, $(\mathrm{SOL}(0, \Gamma'))^* - r' = (\mathrm{SOL}(0, \Gamma'))^*$, and the following holds:

$$\mathrm{rge}(N'\Delta^{-1}(\Theta^\top \otimes I_p)C) = \{[(x-y)(0.5y-0.5x)(x-z)(0.5z-0.5x)]^\top, \forall x, y, z \in \mathbb{R}\}.$$

Since $[2 - 1 \ 2 - 1]^\top \in \mathrm{rge}(N'\Delta^{-1}(\Theta^\top \otimes I_p)C) \cap \mathrm{rint}((\mathrm{SOL}(0, \Gamma'))^* - r')$, Assumption 4 holds.

5. Obviously $f_i$ is locally integrable on $\mathbb{R}_+$ in the first variable and Lipschitz continuous in the second variable, for all $i = 1, 2, 3$. From all the above, the conditions in Corollary 3 are satisfied.

6. Since $N'\Delta^{-1} = \begin{bmatrix} 1 & -0.5 & 0 & 0 \\ 0 & 0 & 1 & -0.5 \end{bmatrix}^\top$, $\Gamma' = \mathrm{diag}(\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}, \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix})$, and $r' = [2 \ -2 \ 2 \ -2]^\top$, so Assumption 6 holds for all $\rho < 1$, using Proposition 10.

7. It can be checked that $f(t, \cdot) = a_1 \cos(a_2t + a_3)(\cdot)$ is QUAD($\bar{P}', Q$) where $\bar{P}' = 1$, $Q = |a_1|$. Therefore, item 3.2 holds. Moreover, item 4 holds for all $\mu \in (0, 0.2)$ since:

$$2Q - 2\gamma\lambda_2(L) + \mu = 2|a_1| - 2 \cdot 0.5 \cdot 1 + \mu \le -0.2 + \mu < 0,$$

where the gain $\gamma = 0.5$, $\lambda_2(L) = 1$, and $|a_1| \le 0.4$, using Remark 11.

8. It can be checked that item 3.3 holds with $\xi = 0.8$ and item 5 holds for all $\rho \in (0.8, 1)$, using Remark 11.

Since all conditions of Theorem 3, Corollary 6 hold, the equilibrium point $e^* = \mathbf{0}_3$ is globally stable in finite time with the upperbound $t^* = \frac{2\sqrt{V(e(0))}}{(\rho-0.8)\sqrt{6}}$, where $e(0) = \Theta_c^\top x(0)$ with $x(0)$ the initial state. This means that the network achieves global synchronization in finite time. In this example, the set of initial states that satisfy Theorem 3 is given by:

$$S = \{[x_1 \ x_2 \ x_3]^\top \in \mathbb{R}^3 : x_1 \ge x_2, x_1 \ge x_3\}.$$

**Numerical simulation of the TVCS (84):** Let us consider the Semi-Backward Euler scheme as given in (81) with time step $h = 10^{-3}s$, $x(0) = [4 \ 1 \ -2]^\top$, and the parameters of $f(t)$, $\xi_1(t, x)$, $\xi_2(t, x)$, $\xi_3(t, x)$ are chosen as below:

$$a_1 = 0.4, a_2 = 10, a_3 = \frac{\pi}{2}; \gamma_{1,1} = \gamma_{3,1} = -0.1, \gamma_{1,2} = \gamma_{2,1} = \gamma_{3,2} = 0.2, \gamma_{2,2} = -0.2. \tag{89}$$

The numerical simulations of the network (84) are presented in Figures 6, and 7. They have been implemented in Python 3 using an LCP Solver given in Appendix C.1. The
simulation scripts can be found in Appendix C.2. Note that, in numerical computations in Python 3, due to floating-point rounding errors, it is common to treat any number with a magnitude smaller than $10^{-13}$ as effectively zero. Therefore, based on the right figure in Figure 6 and the bottom figure in Figure 7, we can see that the network exhibits exponential synchronization, that is, the error state $e(\cdot)$ converges to zero at an exponential rate. The value $t^*$ in Corollary 6 is defined by $t^* = 21.32s$ (i.e., $\rho = 0.999$), where it means that the network achieved global synchronization within the finite time interval $[0, 21.32]$ seconds. This threshold is theoretical, but in practice, the network may synchronize faster than this predicted time. The behavior of the complementarity variables is also depicted in the top figure of Figure 7, where at certain time instances, jump behaviors occur to prevent violation of the linear complementarity constraint in (81). Similar to what was shown in [41, p.3763], the synchronized state $x_{syn}(t)$ in this example is governed by the following dynamics:

$$ \dot{x}_{syn}(t) = \frac{1}{3} \sum_{i=1}^{3} (f(t, x_{syn}(t)) + \xi_{i}(t, x_{syn}(t))), \text{ with an initial state } x_{syn}(0) = \frac{1}{3} \sum_{i=1}^{3} x_{i}(0). $$

**Crucial role of Assumption 6 (Dry Friction Condition):** In this example, we also show that Assumption 6 plays a crucial role in Theorem 3 to ensure the global exponential synchronization of the network. To this end, let us reconsider the couplings given in (83), where $r'_{ij}$ is replaced by $r'_{ij} = [2 \ 2]^T$. Proceeding similarly to the part named **Global synchronization** above, it can be checked that all assumptions in Theorem 3 hold, except Assumption 6 and its related condition in item 5 (which is meaningless if Assumption 6 is violated). The numerical simulations of the states and the norms of the error states in this case are depicted in Figure 8, using Semi-Backward Euler scheme (81) and LCP Solver given in Appendix C.1 to solve LCPs, with time step $h = 10^{-3}s$, initial state $x(0) = [4 \ 1 \ -2]^T$, and the parameters of $f, \xi_i$ defined in (89). Figure 8 shows that the states $x_i(t)$ in the network are not synchronized (left figure), and equivalently that the error states fail to converge to zero over time (right figure). This leads to the conclusion that in this example, in the absence of Assumption 6, the conclusion of Theorem 3 do not hold.

<figure><img src="image_7.png" /><figcaption>Figure 6: Numerical simulations of the states $x_i(t)$ in Example 2 for the disconnected case (left figure) and for the connected case (right figure), using Semi-Backward Euler scheme (81) and LCP Solver given in Appendix C.1 to solve LCPs, with time step $h = 10^{-3}s$, initial state $x(0) = [4 \ 1 \ -2]^T$, and the parameters of $f, \xi_i$ defined in (89).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Numerical simulations of the complementarity variables $\lambda'_{ij,k}(t)$ in the TVCS (84) (top figure), and the log-scaled norms of the error states in Example 2 (bottom figure) where $e(t) = (\Theta_c^\top \otimes I_1)x(t)$. Time step $h = 10^{-3}s$, initial state $x(0) = [4 \ 1 \ -2]^\top$, the parameters of $f, \xi_i$ are defined in (89), using Semi-Backward Euler scheme (81) and LCP Solver given in Appendix C.1 to solve LCPs.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Numerical simulations of the states $x_i(t)$ (left figure) and the log-scaled norms of the error states (right figure) in Example 2 with $e(t) = (\Theta_c^\top \otimes I_1)x(t)$, when agents are interconnected via the couplings defined in (83) with $r'_{ij}$ replaced by $r'_{ij} = [2 \ 2]^\top$. Semi-Backward Euler scheme (81) and LCP Solver given in Appendix C.1 to solve LCPs are used, with time step $h = 10^{-3}s$, initial state $x(0) = [4 \ 1 \ -2]^\top$, and the parameters of $f, \xi_i$ defined in (89).</figcaption></figure>

**The various cases of ($d_1, d_2, d_3$):** In the various cases of the triplet ($d_1, d_2, d_3$), it can be designed couplings that ensure the synchronization network. Indeed, let $H'_{ij}, N'_{ij}, r'_{ij}$ be exactly as given in (83), and the matrices $W_{12}, W_{13}, M'_{12}, M'_{13}$ be as given in Table 4. The couplings of the form (28), generated by these matrices, guarantee global finite-time synchronization of the network. Recall that $\gamma = \lambda_{min}(\bar{\Delta}^{-1}W + W^\top\bar{\Delta}^{-\top})$ represents the coupling gain, the matrices $\Gamma'$ and $\Delta$ are defined after (37), and before Assumption 1, respectively. In the cases considered in Table 4, $\Gamma' = \text{diag}(\Gamma'_{12}, \Gamma'_{13})$ where $\Gamma'_{12}, \Gamma'_{13}$ can be found in Table 4. The values of $\Gamma', \Delta, \gamma$, and $\Theta$ given in Table 4, can be used to verify the global finite-time synchronization of the network of agents (82) via Theorem 3.
<table><caption>Table 4: Matrices $W_{12}$, $W_{13}$, and $M'_{12}$, $M'_{13}$ are used for LCP couplings to achieve global finite-time synchronization among systems (82) in the various cases of ($d_1, d_2, d_3$), where $\Gamma' = \text{diag}(\Gamma'_{12}, \Gamma'_{13})$, $\Theta$, $\Delta = \bar{\Delta}^\top$, and $\gamma$ are used to verify the conditions in Theorem 3.</caption><thead><tr><th>Results<br>(d<sub>1</sub>, d<sub>2</sub>, d<sub>3</sub>)</th><th>(0,0,0)</th><th>(0,-1,0)</th><th>(1,0,0)</th><th>(0,-1,1)</th><th>(0,-1,-1)</th></tr></thead><tbody><tr><td>W<sub>12</sub></td><td>1</td><td><sup>1</sup>&frasl;<sub>2</sub></td><td>1</td><td><sup>1</sup>&frasl;<sub>2</sub></td><td><sup>1</sup>&frasl;<sub>2</sub></td></tr><tr><td>W<sub>13</sub></td><td>1</td><td>3</td><td>1</td><td>3</td><td><sup>1</sup>&frasl;<sub>2</sub></td></tr><tr><td>Θ</td><td><table><tbody><tr><td>1</td><td>1</td></tr><tr><td>-1</td><td>0</td></tr><tr><td>0</td><td>-1</td></tr></tbody></table></td><td><table><tbody><tr><td>1</td><td>1</td></tr><tr><td>-1</td><td>0</td></tr><tr><td>0</td><td>-1</td></tr></tbody></table></td><td><table><tbody><tr><td>1</td><td>0</td></tr><tr><td>-1</td><td>1</td></tr><tr><td>0</td><td>-1</td></tr></tbody></table></td><td><table><tbody><tr><td>1</td><td>1</td></tr><tr><td>-1</td><td>0</td></tr><tr><td>0</td><td>-1</td></tr></tbody></table></td><td><table><tbody><tr><td>1</td><td>1</td></tr><tr><td>-1</td><td>0</td></tr><tr><td>0</td><td>-1</td></tr></tbody></table></td></tr><tr><td>Δ = Δ<sup>T</sup></td><td><table><tbody><tr><td>1</td><td>0</td></tr><tr><td>0</td><td>1</td></tr></tbody></table></td><td><table><tbody><tr><td><sup>1</sup>&frasl;<sub>2</sub></td><td>0</td></tr><tr><td>0</td><td>1</td></tr></tbody></table></td><td><table><tbody><tr><td>2</td><td>0</td></tr><tr><td>0</td><td>1</td></tr></tbody></table></td><td><table><tbody><tr><td><sup>1</sup>&frasl;<sub>2</sub></td><td>0</td></tr><tr><td>0</td><td>4</td></tr></tbody></table></td><td><table><tbody><tr><td><sup>1</sup>&frasl;<sub>2</sub></td><td>0</td></tr><tr><td>0</td><td><sup>1</sup>&frasl;<sub>2</sub></td></tr></tbody></table></td></tr><tr><td>γ</td><td>2</td><td>2</td><td>1</td><td><sup>3</sup>&frasl;<sub>2</sub></td><td>2</td></tr><tr><td>Γ'<sub>12</sub></td><td><table><tbody><tr><td>0</td><td>0</td></tr><tr><td>0</td><td>0</td></tr></tbody></table></td><td><table><tbody><tr><td>-8</td><td>4</td></tr><tr><td>4</td><td>-2</td></tr></tbody></table></td><td><table><tbody><tr><td>2</td><td>-1</td></tr><tr><td>-1</td><td><sup>1</sup>&frasl;<sub>2</sub></td></tr></tbody></table></td><td><table><tbody><tr><td>-8</td><td>4</td></tr><tr><td>4</td><td>-2</td></tr></tbody></table></td><td><table><tbody><tr><td>-8</td><td>4</td></tr><tr><td>4</td><td>-2</td></tr></tbody></table></td></tr><tr><td>Γ'<sub>13</sub></td><td><table><tbody><tr><td>0</td><td>0</td></tr><tr><td>0</td><td>0</td></tr></tbody></table></td><td><table><tbody><tr><td>0</td><td>0</td></tr><tr><td>0</td><td>0</td></tr></tbody></table></td><td><table><tbody><tr><td>0</td><td>0</td></tr><tr><td>0</td><td>0</td></tr></tbody></table></td><td><table><tbody><tr><td>1</td><td>-<sup>1</sup>&frasl;<sub>2</sub></td></tr><tr><td>-<sup>1</sup>&frasl;<sub>2</sub></td><td><sup>1</sup>&frasl;<sub>4</sub></td></tr></tbody></table></td><td><table><tbody><tr><td>-8</td><td>4</td></tr><tr><td>4</td><td>-2</td></tr></tbody></table></td></tr><tr><td>M'<sub>12</sub></td><td><table><tbody><tr><td>1</td><td>-1</td></tr><tr><td>-1</td><td>1</td></tr></tbody></table></td><td><table><tbody><tr><td>9</td><td>-5</td></tr><tr><td>-5</td><td>3</td></tr></tbody></table></td><td><table><tbody><tr><td>-1</td><td>0</td></tr><tr><td>0</td><td><sup>1</sup>&frasl;<sub>2</sub></td></tr></tbody></table></td><td><table><tbody><tr><td>9</td><td>-5</td></tr><tr><td>-5</td><td>3</td></tr></tbody></table></td><td><table><tbody><tr><td>9</td><td>-5</td></tr><tr><td>-5</td><td>3</td></tr></tbody></table></td></tr><tr><td>M'<sub>13</sub></td><td><table><tbody><tr><td>1</td><td>-1</td></tr><tr><td>-1</td><td>1</td></tr></tbody></table></td><td><table><tbody><tr><td>1</td><td>-1</td></tr><tr><td>-1</td><td>1</td></tr></tbody></table></td><td><table><tbody><tr><td>1</td><td>-1</td></tr><tr><td>-1</td><td>1</td></tr></tbody></table></td><td><table><tbody><tr><td>0</td><td>-<sup>1</sup>&frasl;<sub>2</sub></td></tr><tr><td>-<sup>1</sup>&frasl;<sub>2</sub></td><td><sup>3</sup>&frasl;<sub>4</sub></td></tr></tbody></table></td><td><table><tbody><tr><td>9</td><td>-5</td></tr><tr><td>-5</td><td>3</td></tr></tbody></table></td></tr></tbody></table>

## 5.3 Second Example

The following example illustrates Corollary 6 in Section 4.1. The difference between this example and Example 2 is that the inputs and outputs of each agents are considered in $\mathbb{R}^2$, with matrices $D_i \neq 0$ and $D_i$ not necessarily symmetric, positive semi-definite or negative semi-definite, for $i \in \mathcal{V}$. The LCP couplings in this example are designed using Algorithm 1, where $H' \neq (N')^\top$. In addition, a structure of the network in this example illustrates and confirms that the construction rules of networks when $D_i \neq 0$ given at the end of Section 3.3.2 can ensure the global finite-time synchronization of the network with the couplings defined by Algorithm 1.

**Example 3.** Let us consider the following six systems:

$$ \left\{ \begin{array}{ll} \dot{x}_i = f(t, x_i) + \xi_i(t, x_i) + 30u_i & , \quad \forall i = 1, \dots, 6, \\ y_i = x_i + D_i u_i & \end{array} \right. \qquad (90) $$

where $f(t, x_i) = A(t)x_i$ with $A(t) = \begin{bmatrix} 0 & 10 \\ -0.1 - 0.6 \sin(10\pi t) & -0.05 \end{bmatrix}$, $D_2 = \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$,
$D_3 = \begin{bmatrix} -2 & 1 \\ 0 & 2 \end{bmatrix}, D_4 = \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}, D_1 = D_5 = D_6 = 0$, and the mismatches between agents are given by:

$\xi_i(t, x_i) = [\gamma_{i,1} \sin(tx_{1,1} + 1 - x_{1,2}) \cos(0.5t - x_{1,2}) \quad \gamma_{i,2} \cos(4t + x_{1,2} - x_{i,1}) \sin(t + 1)]^\top$
where $\gamma_{i,1}, \gamma_{i,2} \in [-9, 9]$, for all $i = 1, \dots, 6$. In this example, we consider a network with the graph in Figure 9.

<figure><img src="image_10.png" /><figcaption>Figure 9: Graph of the network in Example 3, and its oriented incidence matrix $\Theta$.</figcaption></figure>

Let us consider six couplings as in (28), which are defined as follows:

$$ \bullet W_{25} = W_{15} = W_{45} = W_{56} = W_{16} = \mathbf{I}_2, \text{ and } W_{13} = \text{diag}(0.25, 1). $$

$$ \bullet N'_{25} = \begin{bmatrix} 4 & -2 \\ -2 & 1 \\ -2 & 4 \\ 1 & -2 \end{bmatrix}, N'_{45} = \begin{bmatrix} 4 & 2 \\ -2 & -1 \\ 2 & 4 \\ -1 & -2 \end{bmatrix}, N'_{13} = \begin{bmatrix} 1 & 2 \\ -0.5 & -1 \\ 0 & 6 \\ 0 & -3 \end{bmatrix}, N'_{15} = \begin{bmatrix} 2 & 0 \\ -1 & 0 \\ 0 & 2 \\ 0 & -1 \end{bmatrix}, $$
and $N'_{16} = N'_{56} = N'_{15}$.

$$ \bullet H'_{13} = \begin{bmatrix} 1 & -0.5 & 0.5 & -0.25 \\ 0 & 0 & 6 & -3 \end{bmatrix}, H'_{25} = N'_{25}^\top, H'_{15} = N'_{15}^\top, H'_{45} = N'_{45}^\top, H'_{56} = N'_{56}^\top, $$
and $H'_{16} = N'_{25}^\top$.

$$ \bullet M'_{25} = \begin{bmatrix} -11 & 5 & 12 & -6 \\ 5 & -2 & -6 & 3 \\ 12 & -6 & -11 & 5 \\ -6 & 3 & 5 & -2 \end{bmatrix}, M'_{15} = M'_{56} = M'_{16} = \begin{bmatrix} 1 & -1 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 0 & 0 & 1 & -1 \\ 0 & 0 & -1 & 1 \end{bmatrix}, $$
$$ M'_{45} = \begin{bmatrix} -11 & 5 & -12 & 6 \\ 5 & -2 & 6 & -3 \\ -12 & 6 & -11 & 5 \\ 6 & -3 & 5 & -2 \end{bmatrix}, \text{ and } M'_{13} = \begin{bmatrix} 5 & -3 & -10 & 5 \\ -3 & 2 & 5 & -2.5 \\ 0 & 0 & -23 & 11 \\ 0 & 0 & 11 & -5 \end{bmatrix}. $$

$$ \bullet r_{25} = r_{15} = r_{45} = r_{56} = r_{16} = r_{13} = \begin{bmatrix} 2 & -2 & 2 & -2 \end{bmatrix}^\top. $$

Let $W = \text{diag}(W_{25}, W_{15}, W_{45}, W_{56}, W_{16}, W_{13})$. A similar notation holds for $H', N', M'$.

Let $r' = [r_{25}^\top \ r_{15}^\top \ r_{45}^\top \ r_{56}^\top \ r_{16}^\top \ r_{13}^\top]^\top$. These couplings were determined using Algorithm 1. Let us check the validity of all assumptions in Theorem 3.
1. From the network's structure, all conditions in Lemma 1, and items 1, 3.1 hold.

2. Since $\Delta = \text{diag} \left( \begin{bmatrix} 2 & -1 \\ -1 & 2 \end{bmatrix}, \mathbf{I}_2, \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}, \mathbf{I}_2, \begin{bmatrix} 0.5 & 1 \\ 0 & 3 \end{bmatrix}, \mathbf{I}_2 \right)$ and the matrix $\Gamma'$ is block diagonal with blocks given by $\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$. Hence Assumptions 1 and 2 hold.

3. It can be checked that Assumption 3' holds with $P' = \frac{1}{30}I_{12}$. Hence Assumption 3 holds, and item 2 holds with $\bar{P}' = \frac{1}{30}I_2$. This implies that $c_2 = 0$ where $c_2$ defined in Theorem 3.

4. Since $\Gamma'$ in this example is similar to $\Gamma'$ in Example 1, proceeding similarly to Example 1, it can be checked that Assumption 4 holds.

5. Obviously $f_i := f+\xi_i$ is locally integrable on $\mathbb{R}_+$ in the first variable and Lipschitz continuous in the second variable, for all $i = 1, ..., 6$. From all the above, the conditions in Corollary 3 are satisfied.

6. Since $N'\Delta^{-1} = \text{diag} \left( \begin{bmatrix} 2 \\ -1 \end{bmatrix}, \dots, \begin{bmatrix} 2 \\ -1 \end{bmatrix} \right)$, $\Gamma' = \text{diag} \left( \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}, \dots, \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \right)$, and $r' = [2 \ 2 \ \dots \ 2 \ 2]^\top$, Assumption 6 holds for all $\rho < 2$ using Proposition 10.

7. Using Proposition 1, it can be checked that $f(t, \cdot) = A(t)(\cdot)$ is QUAD($\bar{P}', Q$) where $\bar{P}' = \frac{1}{30}$, $Q = \frac{1}{30}\text{diag}(5.35, 5.3)$. Therefore, item 3.2 holds. Moreover, item 4 also holds for $c_2 = 0$ and $\mu = 0.01$ since:

$$2Q - \gamma\lambda_2(L)\mathbf{I}_2 + \frac{1}{30}\mu\mathbf{I}_2 = \begin{bmatrix} -0.0203 & 0 \\ 0 & -0.0237 \end{bmatrix} \preceq 0,$$

where the gain $\gamma = 0.5976$, $\lambda_2(L) = 0.6314$.

8. It can be checked that item 3.3 holds with $\xi = 0.6\sqrt{2}$ and item 5 holds for all $\rho \in (1.69, 2)$ where $c_2 = 0$.

Since all conditions of Theorem 3, Corollary 6 hold, the network achieves global synchronization in finite time. The set of initial states that satisfy Theorem 3 is given by:

$$S = \{[x_1^\top \dots x_6^\top]^\top \in \mathbb{R}^{12} : x_1 \ge x_3, x_1 \ge x_5 \ge x_6, x_2 \ge x_5, x_4 \ge x_5\},$$

where $x_i \ge x_j$ means that $x_{i,1} - x_{j,1} \ge 0$ and $x_{i,2} - x_{j,2} \ge 0$, where $x_i = [x_{i,1} \ x_{i,2}]^\top$, $x_j = [x_{j,1} \ x_{j,2}]^\top$, for $i, j \in \{1, ..., 6\}$ and $i \neq j$.

**Numerical simulation of the network:** Let us consider the Semi-Backward Euler scheme (81) with time step $h = 10^{-3}s$, the initial state $x(0)$ is given as below:

$$x(0) = \begin{bmatrix} 1 & 3 & 2 & 2 & -\frac{1}{5} & -\frac{1}{2} & \frac{3}{2} & \frac{3}{2} & -\frac{1}{2} & 1 & -1 & -\frac{3}{2} \end{bmatrix}^\top, \quad (91)$$
and the parameters $\gamma_{i,1}$, $\gamma_{i,2}$ are given in Appendix C.2. The numerical simulations of the network are presented in Figures 10, 11, 12, and 13. They have been implemented in Python 3 using an LCP Solver (as given in Appendix C.1) to solve LCPs. The simulation scripts can be found in Appendix C.2. Comparing the states $x_i(t)$ before and after the connection, as depicted in Figures 10, 11, respectively, it can be seen that they have been synchronized within a very short time (approximately 0.02s). The behavior of the mismatches (known perturbations) $\xi_i$ is also depicted in Figure 12 with an oscillation amplitude in [-9,9], corresponding to the amplitude of the parameters $\gamma_{i,1}$ and $\gamma_{i,2}$ defined before Figure 9. Figure 13 illustrates the convergence behavior of the error state $e(t) = (\Theta^\top \otimes I_2)x(t)$ on a logarithmic scale, where the linear decay of the log-scaled norm of $e(t)$ in this semi-log plot indicates that it converges to zero at an exponential rate. It is also worth noting that, in numerical computations using Python 3, due to floating-point rounding limitations, any number with magnitude smaller than $10^{-13}$ is typically considered to be numerically zero. Therefore, based on the numerical simulation results in Figure 13, $e(t)$ can be regarded as zero for all $t \ge 0.02s$, within numerical precision.

<figure><img src="image_11.png" /><figcaption>Figure 10: Numerical simulations of the states $x_i(t) = [x_{i,1}(t) \ x_{i,2}(t)]^\top$ in Example 3 for the disconnected case, using Forward Euler scheme with time step $h = 10^{-3}s$, initial state $x(0)$ as given in (91), and the parameters $\gamma_{i,1}$, $\gamma_{i,2}$ as given in Appendix C.2.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Numerical simulations of the states $x_i(t) = [x_{i,1}(t) \ x_{i,2}(t)]^\top$ in Example 3 for the connected case in different coordinate regions, using Semi-Backward Euler scheme (81) and LCP Solver given in Appendix C.1 to solve LCPs, with time step $h = 10^{-3}s$, initial state $x(0)$ as given in (91), and the parameters $\gamma_{i,1}$, $\gamma_{i,2}$ as given in Appendix C.2.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 12: Numerical simulations of the mismatches $\xi_i(t) = [\xi_{i,1}(t) \xi_{i,2}(t)]^\top$ in Example 3, using Semi-Backward Euler scheme (81) with time step $h = 10^{-3}s$, and LCP Solver given in Appendix C.1 to solve LCPs. The initial state $x(0)$ and the parameters of $\xi_i$ are given in (91), and in Appendix C.2, respectively.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: Numerical simulations of the log-scaled norm of the error state $e(t)$ in Example 3 where $e(t) = (\Theta_c^\top \otimes I_2)x(t)$, using Semi-Backward Euler scheme (81) with time step $h = 10^{-2}s$, and LCP Solver given in Appendix C.1 to solve LCPs. The initial state $x(0)$ and the parameters of $\xi_i$ are given in (91), and in Appendix C.2, respectively.</figcaption></figure>

## 5.4 Third Example: Circuits

The following example illustrates Theorem 3, whose core idea is to show that the couplings (6) can be designed to synchronize a large-scale network of agents (5). Although Example 4 below considers only 30 agents, it can be extended to larger cases. For a heterogeneous network of agents given in Example 4 with 300 agents and graphs containing over 2,000 edges, it can be checked that the global exponential synchronization of the network is achieved, using Semi-Backward Euler scheme (81) and the CVXPY Tool (for solving the LCPs) given in Appendix C.1.

**Example 4 (Circuits synchronization).** Let us consider 30 time-varying electrical circuits as in Figure 17. For each $i = 1, 30$, let $x_{i,1}$ be the charge of the capacitor, and $x_{i,2}$ be the current through the inductor; $(u_i + \bar{\xi}_{i,2})$ be the input current, and $y_i$ be the output voltage; $I_i^S(t)$, $\bar{\xi}_{i,1}$ be the currents of the AC current sources. In practice, $\bar{\xi}_{i,1}$ is considered as a type of accumulated perturbation resulting from the uncertainty of individual agents, while $\bar{\xi}_{i,2}$ represents a perturbation arising from the input. Therefore, they could be time-varying and dependent on $x_{i,1}$, $x_{i,2}$. This practical framework will be examined in Example 7 when studying practical implementations. Meanwhile, within the current framework, these represent measurable components that contribute
to increasing the network's heterogeneity, alongside the heterogeneity of the AC current sources. In this example, these components are informally referred to as known perturbations of the network. Following Kirchhoff's laws (see details in Appendix A.1), the dynamics of each circuit is given by:

$$ egin{cases} 
\dot{x}_i(t) = \underbrace{\begin{bmatrix} 0 & -1 \\ \frac{1}{CL(t)} & -\frac{R(t)}{L(t)} - \frac{\dot{L}(t)}{L(t)} \end{bmatrix}}_{A(t)} x_i(t) + \underbrace{\begin{bmatrix} I_i^S(t) \\ 0 \end{bmatrix}}_{\bar{B}} + \underbrace{\begin{bmatrix} 1 \\ 0 \end{bmatrix}}_{\bar{B}} (u_i(t) + \bar{\xi}_i(t, x_i)) \\ 
y_i(t) = \underbrace{\begin{bmatrix} \frac{1}{C} & 0 \end{bmatrix}}_{\bar{C}} x_i(t), 
\end{cases} \quad (92) $$

for $i = 1, ..., 30$, where $x_i = [x_{i,1} \ x_{i,2}]^\top$, $\bar{\xi}_i = \bar{\xi}_{i,1} + \bar{\xi}_{i,2}$, and the other variables are defined as:

$$ egin{cases} 
\text{Time-varying inductance: } \mathbf{L}(t) = \mathbf{L} + \alpha \sin(0.2t + 2) \cos(3t - 1). \\ 
\text{Time-varying resistance: } \mathbf{R}(t) = \mathbf{R} + \beta \sin(2t - 3) \cos(-4t + 0.1). \\ 
\text{AC current sources: } \mathbf{I}_i^S(t) = \delta_1 \sin(\delta_2 t + \delta_3) + 0.2 \cos(\delta_{4,i} t - 1) \sin(\delta_{5,i} t + 2). \\ 
\text{Nominal inductance, resistance, and capacitance: } \mathbf{L}, \mathbf{R}, \mathbf{C} > 0. \\ 
\text{Perturbations: } \bar{\xi}_{i,1}(t, x_i) = \gamma_{i,1} \sin(0.4t + x_{i,1} - 2x_{i,2}); \\ 
\bar{\xi}_{i,2}(t, x_i) = \gamma_{i,2} \cos(-2t + 0.2x_{i,2}). \\ 
\text{Parameters: } \alpha, \beta \in [-\frac{2}{5}, \frac{2}{5}] \text{ and } |\alpha| < \mathbf{L}, |\beta| < \mathbf{R}; \delta_j \in \mathbb{R}, \text{ for all } j = 1, 2, 3; \\ 
\delta_{4,i}, \delta_{5,i} \in [0, 10], \text{ and } \gamma_{i,1}, \gamma_{i,2} \in [-\frac{2}{5}, \frac{2}{5}], \text{ for all } i = \overline{1, 30}. 
\end{cases} \quad (93) $$

Clearly the systems in (92) have the form of (5) where $f_i(t, x_i) = A(t)x_i + \bar{B}\xi_i(t, x_i)$, $\xi_i(t, x_i) = \bar{\xi}_i(t, x_i) + \mathbf{I}_i^S(t)$, $B_i = \bar{B}$, $C_i = \bar{C}$, $D_i = 0$, for all $i = \overline{1, 30}$, and $A(t)$, $\bar{B}$, $\bar{C}$ are given in (92). In this example, let $\mathbf{C} = 2$, $\mathbf{L} = \mathbf{R} = \frac{3}{2}$, $\alpha = \beta = \frac{2}{5}$. A network with the structure of a Newman-Watts-Strogatz small-world type [88], where each vertex is connected to its $k = 6$ nearest neighbors and with a probability of extra connections $p = 0.5$, see Figure 14, is considered.
<figure><img src="image_15.png" /><figcaption>Figure 14: Undirected graph $G = (V, E)$ illustrating the structure of a Newman-Watts-Strogatz small-world type used in Example 4, where each vertex is connected to its $k = 6$ nearest neighbors and with a probability of extra connections $p = 0.5$. This configuration corresponds $|V| = 30, |E| = 129$. This figure is generated by `networkx.newman_watts_strogatz_graph` in Python 3 with `seed = 1`.</figcaption></figure>

Let us consider couplings as in (28), where for all $i, j \in V, i < j$, and $(i, j) \in E$, we have:

$$W_{ij} = 15, \quad H'_{ij} = [3 \ -2], \quad N'_{ij} = \begin{bmatrix} 3 & 2 \\ -2 & 2 \end{bmatrix}, \quad M'_{ij} = \begin{bmatrix} 1 & 2 \\ 2 & 4 \end{bmatrix}, \quad \text{and} \quad r'_{ij} = \begin{bmatrix} -5 \\ -10 \end{bmatrix}.$$

The dynamics of the network is obtained from (37). Let us now check the validity of all assumptions in Theorem 3.

1. From the network's structure, all conditions in Lemma 1, and items 1, 3.1 hold.

2. Since $\Delta = I_{|E|}$ and $\Gamma' = M' = \text{diag}(M'_{ij}) \ge 0$, Assumptions 1 and 2 hold.

3. It can be checked that Assumption 3' holds with $P' = \frac{1}{2}I_{2|V|}$. Therefore, Assumption 3 holds, and item 2 holds with $\bar{P}' = \frac{1}{2}I_2$.

4. Using (40) in Remark 7, $(\text{SOL}(0, \Gamma'))^* = (\text{SOL}(0, \Gamma'))^* - r' = \mathbb{R}^{2|E|}$. Therefore, Assumption 4 holds.

5. Obviously $f_i$ is locally integrable on $\mathbb{R}_+$ in the first variable and Lipschitz continuous in the second variable, for all $i = 1, 30$. From all the above, the conditions in Corollary 3 are satisfied.

6. Using Proposition 9, it can be checked that Assumption 6 holds for all $\rho < 5$.

7. Let $Q = \text{diag}(0.5, -0.05)$. It can be checked that $\frac{1}{4}(A(t)+A(t)^\top) \ll Q$ for all $t \ge 0$. Therefore, $f(t, \cdot) = A(t)(\cdot)$ is QUAD($\bar{P}', Q$) where $\bar{P}' = \frac{1}{2}I_2$. This implies that