item 3.2 holds. Moreover, it can be checked that item 4 holds for $\mu = 0.01$ since:

$$2Q - \gamma\lambda_2(L)\bar{C}^\top\bar{C} + \mu\bar{P}' = \begin{bmatrix} 1 - \frac{1}{4}\gamma\lambda_2(L) + \frac{1}{2}\mu & 0 \\ 0 & \frac{1}{2}\mu - 0.01 \end{bmatrix} = \begin{bmatrix} -21.008 & 0 \\ 0 & -0.005 \end{bmatrix} \lesssim 0,$$

where the gain $\gamma = 30$, $\lambda_2(L) \approx 2.9344$, using Remark 11. The eigenvalues of the matrix $(\frac{1}{4}(A(t) + A(t)^\top) - Q)$ are shown in Figure 17 to illustrate its negative definiteness.

8. Item 3.3 holds with $\xi = 2$ and item 5 holds for all $\rho \in [4.45, 5)$, using Remark 11.

Since all conditions of Theorem 3 hold, $e^* = 0_{|\varepsilon_c|}$ is the unique equilibrium point of the error dynamics and it is globally exponentially stable. This means that the network achieves global exponential synchronization. In this example, the set of initial states satisfying Theorem 3 is the entire space $\mathbb{R}^{2|\nu|}$.

**Numerical simulation of the network:** Let us consider the Semi-Backward Euler scheme (81) with time step $h = 10^{-2}s$, $\delta_1 = 4$, $\delta_2 = 2$, $\delta_3 = 1$, and $\Theta \in \mathbb{R}^{30 \times 129}$ a oriented incidence matrix associated with $\mathcal{G}$ that satisfies (12). The initial state $x(0)$, and the parameters $\gamma_{i,1}$, $\gamma_{i,2}$, $\delta_{4,i}$, $\delta_{5,i}$ are given in Appendix C.2. The numerical simulations of the network are presented in Figures 15, and 16. They have been implemented in Python 3 using CVXPY Tool (as given in Appendix C.1) to solve LCPs. The simulation scripts can be found in Appendix C.2. In Figure 15, the trajectory with the larger amplitude corresponds to $x_{i,1}(t)$ after synchronization, while the other corresponds to $x_{i,2}(t)$. Figure 16 illustrates the convergence behavior of the error state $e(t) = (\Theta^\top \otimes I_2)x(t)$ on a logarithmic scale. 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. This behavior is further confirmed in Figure 15, where the trajectories $x_{i,1}(t)$ achieve synchronization within a very short time (approximately 0.02s), while $x_{i,2}(t)$ reach synchronization within around 5s. 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 16, $e(t)$ can be regarded as zero for all $t \ge 40s$, within numerical precision.
<figure><img src="image_16.png" /><figcaption>Figure 15: Numerical simulations of the states $x_{i,1}(t)$, $x_{i,2}(t)$ in Example 4, using Semi-Backward Euler scheme (81) and CVXPY Tool given in Appendix C.1 to solve LCPs, with time step $h = 10^{-2}s$, $\delta_1 = 4$, $\delta_2 = 2$, $\delta_3 = 1$, the initial state $x(0)$ and the parameters $\gamma_{i,1}, \gamma_{i,2}, \delta_{4,i}, \delta_{5,i}$ given in Appendix C.2.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16: Numerical simulations of the log-scaled norm of the error state $e(t)$ in Example 4 where $e(t) = (\Theta_c^\top \otimes I_2)x(t)$, using Semi-Backward Euler scheme (81) and CVXPY Tool given in Appendix C.1 to solve LCPs, with time step $h = 10^{-2}s$, $\delta_1 = 4$, $\delta_2 = 2$, $\delta_3 = 1$, the initial state $x(0)$ and the parameters $\gamma_{i,1}, \gamma_{i,2}, \delta_{4,i}, \delta_{5,i}$ given in Appendix C.2.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 17: Circuit diagram of agent i in Example 4 (left figure), and the trajectories of the eigenvalues of the matrix $\frac{1}{4}(A(t) + A(t)^T) - Q$ (right figure).</figcaption></figure>

**The structure of the couplings:** Recall from (7) that $v_{ij} = y_i - y_j = \bar{C}(x_i - x_j)$. Using Lemma 1, the couplings (6) in this example are given by:

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

Let $\eta'_{ij}(t) = [3 - 2]\lambda'_{ij}(t)$, then $\eta'_{ij} : \mathbb{R} \Rightarrow \mathbb{R}$ is a set-valued mapping in the variable $v_{ij}$, defined by:

$$ v_{ij} \mapsto \eta'_{ij}(v_{ij}) = \begin{cases} [-5, 15], & \text{if } v_{ij} = 0 \\ -v_{ij} - 5, & \text{if } v_{ij} > 0 \\ -9v_{ij} + 15, & \text{if } v_{ij} < 0 \end{cases}. $$

The graphs of $\eta'_{ij}$ and $\eta_{ij}$ are given in Figure 18, which depicts the relationship between the input $v_{ij}(t)$ and the output $\eta_{ij}(t)$ of the coupling. Note that the axes in Figure 18 use different unit scales to better illustrate the graph.

<figure><img src="image_19.png" /><figcaption>Figure 18: Input-output relation at the coupling between agents i and j in Example 4.</figcaption></figure>
## 5.5 Fourth Example: Chua's Circuits

The following example illustrates Corollary 6, in which a network of Chua's circuits is considered. Chua's circuits are known to exhibit complex chaotic behavior, such as the Double Scroll Attractor - two symmetric spirals representing aperiodic, sensitive, yet bounded dynamics typical of nonlinear chaotic systems [89]. The main difference between this example and Examples 2, 3 and 4 lies in the fact that the agents under consideration are nonlinear, exhibit complex chaotic behaviors. In addition, similarly to Examples 2, 3, the matrices $B_i$ (defined in (5)) corresponding to the inputs $u_i(t)$ in this example are identical and invertible, which enables us to take into account full-state perturbations.

**Example 5.** (Chua's Circuits Synchronization). Let us consider ten Chua's circuits as in Figure 19. For each $i = 1, \overline{1,10}$, let $x_{i,1}, x_{i,2}$ be the voltages across capacitors $C_1, C_2$, respectively, and $x_{i,3}$ the current through the inductor $L$; $y_{i,1}, y_{i,2}$ be the output voltages, and $y_{i,3}$ be the output current; $(u_{i,1} + \xi_{i,1,2}), (u_{i,2} + \xi_{i,2,2})$ be the input currents, and $(u_{i,3} + \xi_{i,3,2})$ be the input voltage; $\xi_{i,1,1}, \xi_{i,2,1}$ be the currents of the AC sources, and $\xi_{i,3,1}$ be the voltage of the DC source. In practice, the quantities $\xi_{i,1,2}, \xi_{i,2,2}, \xi_{i,3,2}$ represent perturbations originating from the inputs, while $\xi_{i,1,1}, \xi_{i,2,1}, \xi_{i,3,1}$ correspond to perturbations caused by uncertainties related to the agents. Therefore, they could be time-varying and dependent on $x_{i,1}, x_{i,2}, x_{i,3}$. This practical framework will be examined in Example 8 when studying practical implementations. Meanwhile, within the current framework, these represent measurable components that contribute to increasing the network's heterogeneity. In this example, these components are informally referred to as known perturbations of the network.

<figure><img src="image_20.png" /><figcaption>Figure 19: Circuit diagram of agent $i$ in Example 5.</figcaption></figure>

Following Kirchhoff's laws (see details in Appendix A.2), the dynamics of each Chua's
circuit is given by:

$$
\left\{
\begin{aligned}
\dot{x}_i(t) &= \underbrace{\begin{bmatrix} -\frac{1}{RC_1} & \frac{1}{RC_1} & 0 \\ \frac{1}{RC_2} & -\frac{1}{RC_2} & -\frac{1}{C_2} \\ 0 & \frac{1}{L} & 0 \end{bmatrix}}_{A_1} x_i(t) + \underbrace{\begin{bmatrix} -\frac{1}{C_1} \\ 0 \\ 0 \end{bmatrix}}_{A_2} g(x_{i,1}(t)) \\
&\quad + \underbrace{\begin{bmatrix} \frac{1}{C_1} & 0 & 0 \\ 0 & \frac{1}{C_2} & 0 \\ 0 & 0 & \frac{1}{L} \end{bmatrix}}_{\bar{B}} (u_i(t) + \xi_i(t, x_i(t))) \\
y_i(t) &= \underbrace{I_3}_{\bar{C}} x_i(t),
\end{aligned}
\right.
\qquad (95)
$$

for $i = 1, ..., 10$, where $x_i = [x_{i,1} \ x_{i,2} \ x_{i,3}]^\top$; $u_i = [u_{i,1} \ u_{i,2} \ u_{i,3}]^\top$; $\xi_i = [\xi_{i,1} \ \xi_{i,2} \ \xi_{i,3}]^\top$ with $\xi_{i,j} = \xi_{i,j,1} + \xi_{i,j,2}$, for $j = 1, 2, 3$; $y_i = [y_{i,1} \ y_{i,2} \ y_{i,3}]^\top$; $C_1, C_2$ are capacitances, $R$ is the resistance, $L$ is the inductance; and $g: \mathbb{R} \to \mathbb{R}$ is a nonlinear function that models the voltage-current characteristic of the nonlinear resistor, with voltage as the input and current as the output.

Clearly the systems (95) have the form (5) where $f_i(t, x_i) = f(x_i) + \bar{B}\xi_i(t, x_i)$ with $f(x_i) = A_1x_i + A_2g(x_{i,1})$, $B_i = \bar{B}$, $C_i = \bar{C}$, $D_i = \mathbf{0}_{3\times3}$, for all $i = 1, ..., 10$, and $A_1, A_2, \bar{B}, \bar{C}$ are given in (95). In this example, let $\mathbf{R} = \mathbf{C}_2 = 1$, $\mathbf{C}_1 = \frac{1}{16}$, $L = \frac{1}{25}$, the function $g$ be defined by $g: \mathbb{R} \to \mathbb{R}$, $v \mapsto g(v) := -\frac{1}{2}v - \frac{7}{20}(|v+1| - |v-1|)$, and the known perturbations $\xi_{i,1}(t, x_i)$, $\xi_{i,2}(t, x_i)$, and $\xi_{i,3}(t, x_i)$ be defined by:

$$
\left\{
\begin{aligned}
\xi_{i,1}(t, x_i) &= \gamma_{i,1,1} \sin(2t + 0.1x_{i,1}) \cos(-2t) + \gamma_{i,1,2} \sin(x_{i,2}t - 1) \cos(t - x_{i,1}), \\
\xi_{i,2}(t, x_i) &= \gamma_{i,2,1} \cos(-t + x_{i,2} - 0.3x_{i,1}) \sin(-x_{i,1}) + \gamma_{i,2,2} \sin(-0.2t - x_{i,3}), \\
\xi_{i,3}(t, x_i) &= \gamma_{i,3,1} \sin(0.5t - 0.3x_{i,2}x_{i,3}) \cos(0.5t) + \gamma_{i,3,2} \cos(t + x_{i,3} + 2x_{i,2}),
\end{aligned}
\right.
\qquad (96)
$$

for $i = 1, ..., 10$, where $\gamma_{i,j,k} \in [-0.1, 0.1]$, $\forall i = 1, ..., 10$, $\forall j = 1, 2, 3$, $\forall k = 1, 2$. A network with the structure of a Newman-Watts-Strogatz small-world type [88], where each vertex is connected to its $k = 2$ nearest neighbors and with a probability of extra connections $p = 0.7$, see Figure 20, is considered.
<figure><img src="image_21.png" /><figcaption>Figure 20: Undirected graph $\mathcal{G} = (\mathcal{V}, \mathcal{E})$ illustrating the structure of a Newman-Watts-Strogatz small-world type used in Example 4, where each vertex is connected to its $k=2$ nearest neighbors and with a probability of extra connections $p = 0.7$. This configuration corresponds $|\mathcal{V}| = 10$, $|\mathcal{E}| = 16$. 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 \mathcal{V}$, $i < j$, and $(i, j) \in \mathcal{E}$, we have:

$$W_{ij} = 5I_3, \quad H'_{ij} = I_3 \otimes [1 \ -3], \quad N'_{ij} = (H'_{ij})^\top, \quad M'_{ij} = I_3 \otimes \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}, \quad \text{and } r'_{ij} = \mathbf{1}_3 \otimes \begin{bmatrix} -4 \\ -4 \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_{|\mathcal{E}|}$ and $\Gamma' = M' = \text{diag}(M'_{ij}) \succ 0$, Assumptions 1 and 2 hold.

3. It can be checked that Assumption 3 holds with $P' = (I_{|\mathcal{V}|} \otimes \text{diag}(\frac{1}{16}, \frac{1}{25}))$. Therefore, item 2 holds with $\bar{P}' = \text{diag}(\frac{1}{16}, \frac{1}{25}) \succ 0$.

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

5. It can be checked that $f_k$ is locally integrable on $\mathbb{R}_+$ in the first variable and Lipschitz continuous in the second variable for all $k = \frac{1}{10}$. From all the above, the conditions in Corollary 3 hold.

6. Using Remark 10 and Proposition 9, it can be checked that Assumption 6 holds for all $\rho < 4$.

7. The function $f$, defined after (95), is $\text{QUAD}(\bar{P}', Q)$ with $Q = \begin{bmatrix} 0.2 & 1 & 0 \\ 1 & -1 & 0 \\ 0 & 0 & 0 \end{bmatrix}$. Indeed, since $-g$ is a Lipschitz continuous function with Lipschitz constant 1.2,
it is QUAD(1, 1.2). For all $\alpha_1 = [\alpha_{1,1} \ \alpha_{1,2} \ \alpha_{1,3}]^\top$, $\alpha_2 = [\alpha_{2,1} \ \alpha_{2,2} \ \alpha_{2,3}]^\top \in \mathbb{R}^3$, the following holds:

$$ (\alpha_1 - \alpha_2)^\top \bar{P}'(f(\alpha_1) - f(\alpha_2)) = (\alpha_1 - \alpha_2)^\top \bar{P}' A_1 (\alpha_1 - \alpha_2) \\ + (\alpha_{1,1} - \alpha_{2,1})^\top (-g(\alpha_{1,1}) + g(\alpha_{2,1})) \\ \leq (\alpha_1 - \alpha_2)^\top (\bar{P}' A_1 + \text{diag}(1.2, 0, 0))(\alpha_1 - \alpha_2) \\ = (\alpha_1 - \alpha_2)^\top Q (\alpha_1 - \alpha_2). $$

Therefore, item 3.2 holds. Moreover, item 4 holds for $\mu = 0.01$ since:

$$ 2Q - \gamma \lambda_2(L) \bar{C}^\top \bar{C} + \mu \bar{P}' = \begin{bmatrix} 0.4 - \gamma \lambda_2(L) + \frac{\mu}{16} & 2 & 0 \\ 2 & -2 - \gamma \lambda_2(L) + \mu & 0 \\ 0 & 0 & -\gamma \lambda_2(L) + \frac{\mu}{25} \end{bmatrix} \le 0, $$

where the gain $\gamma = 10$, $\lambda_2(L) \approx 0.9963$, using Remark 11.

8. It can be checked that item 3.3 holds with $\xi = 0.4\sqrt{3}$ and item 5 holds for all $\rho \in [1.48, 4)$, using Remark 11.

Since all conditions of Theorem 3, and Corollary 6 hold, $e^* = \mathbf{0}_{|\varepsilon_c|}$ is the unique equilibrium point of the error dynamics corresponding to the network's dynamics, and it is globally stable in finite time with the upperbound $t^* = \frac{2}{4.464\rho-6.5726}V(e(0))$, where $e(0) = (\Theta_c^\top \otimes I_3)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 satisfying Theorem 3 is the entire space $\mathbb{R}^{3|\nu|}$.

**Numerical simulation of the network:** Let us consider Semi-Backward Euler scheme (81) with time step $h = 5 \times 10^{-4}s$, and $\Theta \in \mathbb{R}^{10 \times 16}$ a oriented incidence matrix associated with $G$ that satisfies (12). The initial state $x(0)$ and the parameters $\gamma_{i,j,k}$ are defined in Appendix C.2. The numerical simulations of the network are presented in Figures 21, 22, and 23. They have been implemented in Python 3 using CVXPY Tool (as given in Appendix C.1) to solve LCPs. The simulation scripts can be found in Appendix C.2. Figure 21 illustrates the three-dimensional states $x_i$, showing the interrelationship between its components and providing a global view of the behavior of the Chua's circuits in 3D space. In the top two figures of Figure 21, it is clear that, prior to the Chua's circuits being interconnected into a network, each circuit exhibits significantly different behaviors and displays Double Scroll Attractors, which are characteristic manifestations of complex chaotic dynamics. The three-dimensional states of each circuit, when connected into a network, become synchronized as shown in the bottom two figures of Figure 21. Here, we can observe trajectories starting from different initial points converging into a new trajectory (the blue curve), which is referred to as the synchronized trajectory. Figure 22 illustrates the time evolution of the states $x_{i,1}(t)$, $x_{i,2}(t)$, and $x_{i,3}(t)$, indicating that all these states reach synchronization
in roughly 0.03 seconds. In addition, in Figure 23, the linear decay of the log-scaled norm of the error state $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 23, $e(t)$ can be regarded as zero for all $t \ge 0.03$, within numerical precision. This illustrates Corollary 6, which states that the network achieves global synchronization within a finite time.

<figure><img src="image_22.png" /><figcaption>Figure 21: Numerical simulations of the three-dimensional states $x_i = [x_{i,1} \ x_{i,2} \ x_{i,3}]^\top$ of Chua's circuits in Example 5 prior to interconnection (top two figures) and post interconnection (bottom two figures), using Semi-Backward Euler scheme (81) with time step $h = 5 \times 10^{-4}$s, and CVXPY Tool given in Appendix C.1 to solve LCPs. The initial state $x(0)$ and the parameters $\gamma_{i,j,k}$ are defined in Appendix C.2.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 22: Numerical simulations of the states $x_{i,1}(t)$, $x_{i,2}(t)$, $x_{i,3}(t)$ in Example 5 in different coordinate regions, using Semi-Backward Euler scheme (81) with time step $h = 5 \times 10^{-4}$s, and CVXPY Tool given in Appendix C.1 to solve LCPs. The initial state $x(0)$ and the parameters $\gamma_{i,j,k}$ are defined in Appendix C.2.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 23: Numerical simulations of the log-scaled norm of the error state $e(t)$ in Example 5 where $e(t) = (\Theta_c^\top \otimes I_3)x(t)$, using Semi-Backward Euler scheme (81) with time step $h = 5 \times 10^{-4}$s, and CVXPY Tool given in Appendix C.1 to solve LCPs. The initial state $x(0)$ and the parameters $\gamma_{i,j,k}$ are defined in Appendix C.2.</figcaption></figure>

# 6 Practical Implementation

In this section, we focus on designing a discrete-time controller $u(t)$ to ensure the synchronization of the given systems of the form of (5). This is crucial for practical applications, and it is known that set-valued feedback controllers (like sliding-mode inputs) are extremely sensitive to the chosen discretization method [50]. Recall that $\Gamma' = N'(\Theta^\top \otimes I_p)D\Xi^{-1}(\Theta \otimes I_p)H' + M'$, $\Xi = I_{p|V|} + (\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$. All conditions in Lemma 1, item 3.1 in Theorem 3, and Assumptions 1, 2, 3', 4 are considered to hold throughout this section. Let us consider the following lemma.

**Lemma 4.** Suppose that all conditions in Lemma 1, and Assumption 1 hold. Then the input $u(t)$ in (11) is rewritten as:

$$u(t) \in -\tilde{\mathcal{M}}(Cx(t)) - \Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)Cx(t), \quad (97)$$

where $\tilde{\mathcal{M}}(Cx(t)) := \Xi^{-1}(\Theta \otimes I_p)H'\Phi'(N'\Delta^{-1}(\Theta^{\top} \otimes I_p)Cx(t))$, the state $x(t)$ and the operator $\Phi'$ are defined in (9) and (39), respectively.

*Proof.* Since all conditions in Lemma 1 and Assumption 1 hold, it follows from (36) that:

$$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 (98)$$

where $x(t)$ is the state defined in (9), $\eta'(t) = H'\lambda'(t)$, and $\lambda'(t)$ is defined in (34). In addition, it follows from (37) that $\lambda'(t)$ is a solution of the following LCP:

$$\left\{ \begin{array}{l} 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 (99)$$

Let $\Phi'$ be defined as in (39). Proceeding similarly to (14), the LCP (99) is equivalent to:

$$\lambda'(t) \in -\Phi'(N'\Delta^{-1}(\Theta^{\top} \otimes I_p)Cx(t)).$$

From the above, the input $u(t)$ is rewritten equivalently as:

$$u(t) \in -\Xi^{-1}(\Theta \otimes I_p)H'\Phi'(N'\Delta^{-1}(\Theta^{\top} \otimes I_p)Cx(t)) - \Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)Cx(t).$$

This is equivalent to: $u(t) \in -\tilde{\mathcal{M}}(Cx(t)) - \Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)Cx(t)$. $\square$

Recall that in this section, item 3.1 in Theorem 3 is assumed to hold. Therefore $F(t, x)$ (defined in (8)) can be rewritten in the form:

$$F(t, x) = F_{\text{nom}}(t, x) + B\xi_{\text{per}}(t, x), \quad \forall x = [x_1^{\top} \dots x_{|V|}^{\top}]^{\top} \in \mathbb{R}^{n|V|}, \forall t \ge 0,$$

where $F_{\text{nom}}(t, x) := [f(t, x_1)^{\top} \dots f(t, x_{|V|})^{\top}]^{\top} + B\xi_{\text{nom}}(t, x)$ with the function $f$ defined in Theorem 3, and $\xi_{\text{nom}}, \xi_{\text{per}}: \mathbb{R}_+ \times \mathbb{R}^{n|V|} \to \mathbb{R}^{p|V|}$ are functions such that:

$$\xi_{\text{nom}}(t, x) + \xi_{\text{per}}(t, x) = [\xi_1(t, x_1)^{\top} \dots \xi_{|V|}(t, x_{|V|})^{\top}]^{\top}, \quad \forall x = [x_1^{\top} \dots x_{|V|}^{\top}]^{\top}, \forall t \ge 0,$$

with $\xi_i$ defined in Theorem 3 for all $i = 1, \dots, |V|$. In the framework of this section, $F_{\text{nom}}$ denotes the nominal vector field, and $\xi_{\text{per}}$ denotes the perturbations. From Lemma 4 and the compact form of agents in (9), the dynamics of the network is rewritten as:

$$\left\{ \begin{array}{l} \dot{z}(t) = RF_{\text{nom}}(t, R^{-1}z(t)) + RB\xi_{\text{per}}(t, R^{-1}z(t)) + RBu(t), \\ u(t) \in -\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)CR^{-1}z(t) - \tilde{\mathcal{M}}(CR^{-1}z(t)), \end{array} \right. \quad (100)$$
where $z(t) := Rx(t)$ with $R = R^\top > 0$ a square root of $P'$ given in Assumption 3'. Using Forward (full-explicit) Euler scheme for the plant (the agents), and Semi-Backward (semi-implicit) Euler scheme for the input, the time-discretization of (100), (assuming that perturbations are known, i.e., $\xi_{per}(t, R^{-1}z)$ is known) is given by:

$$ \left\{ \begin{array}{l} z_{k+1} = z_k + h(RF_{nom}(t_k, R^{-1}z_k) + RB\xi_{per}(t_k, R^{-1}z_k) + RBu_k), \\ u_k \in -\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)CR^{-1}z_k - \tilde{M}(CR^{-1}z_{k+1}), \end{array} \right. \qquad (101) $$

where $h > 0$ is the time step, $(t_0, t_1, ..., t_n)$ is a partition of $[0, T]$ for some $T > 0$, and $t_{k+1} - t_k = h$, $n = \frac{T}{h}$, $z_k := z(t_k)$, and $u_k := u(t_k)$, for all $k = \overline{0, n}$. In practice perturbations $\xi_{per}$ are unknown, hence (101) cannot be implemented. Following ideas introduced in [90] (see [50, 91] for details and more references on discrete-time realization and implementation of set-valued feedback controllers), let us consider the virtual state $\tilde{z}_{k+1}$ corresponding to the nominal system, defined by:

$$ \left\{ \begin{array}{l} \tilde{z}_{k+1} = z_k + h(RF_{nom}(t_k, R^{-1}z_k) + RBu_k), \\ u_k \in -\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)CR^{-1}z_k - \tilde{M}(CR^{-1}\tilde{z}_{k+1}). \end{array} \right. \qquad (102) $$

Let $\mathcal{M}(\cdot) := RB\tilde{M}(CR^{-1}(\cdot))$, then the system (102) is rewritten as:

$$ \tilde{z}_{k+1} \in z_k + hRF_{nom}(t_k, R^{-1}z_k) - hRB\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)CR^{-1}z_k - h\mathcal{M}(\tilde{z}_{k+1}). \qquad (103) $$

The inclusion (103) implies that:

$$ \tilde{z}_{k+1} \in \mathcal{J}_{h\mathcal{M}}(z_k + hRF_{nom}(t_k, R^{-1}z_k) - hRB\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)CR^{-1}z_k), \qquad (104) $$

where resolvents are defined as in Definition 5. In addition, it follows from (102) that:

$$ Bu_k = \frac{1}{h}R^{-1}(\tilde{z}_{k+1} - z_k - hRF_{nom}(t_k, R^{-1}z_k)). \qquad (105) $$

Recall that the matrix $B \in \mathbb{R}^{n|\mathcal{V}| \times p|\mathcal{V}|}$ is defined in (9).

**Lemma 5.** Suppose that Assumptions 2, 3', 4 hold, and $\text{rank}(B) = p|\mathcal{V}|$. Then $\mathcal{M}$ defined after (102) is a maximal monotone operator, and the following hold:

$$ \tilde{z}_{k+1} = \mathcal{J}_{h\mathcal{M}}(z_k + hRF_{nom}(t_k, R^{-1}z_k) - hRB\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)CR^{-1}z_k), \qquad (106) $$

$$ u_k = \frac{1}{h}(B^\top B)^{-1}B^\top R^{-1}(\tilde{z}_{k+1} - z_k - hRF_{nom}(t_k, R^{-1}z_k)), \qquad (107) $$

where $\tilde{z}_{k+1}$, and $u_k$ are defined in (102).

*Proof.* Let $G := N'\Delta^{-1}(\Theta^\top \otimes I_p)CR^{-1}$. Since Assumption 3' holds, the operator $\mathcal{M}$ defined after (102) is redefined by $\mathcal{M}(x) = G^\top \Phi'(Gx)$, for all $x \in \mathbb{R}^{n|\mathcal{V}|}$, where the operator $\Phi'$ is defined in (39). Since Assumption 2 holds, $\Phi'$ is maximal monotone,
and Dom $\Phi' = (\text{SOL}(0, \Gamma'))^* - r'$, see Remark 7. Since Assumption 4 holds, $\mathcal{M}$ is maximal monotone [54, Theorem 12.43]. Therefore, $\mathcal{J}_{h\mathcal{M}}$ is a single-valued mapping with $\text{Dom}(\mathcal{J}_{h\mathcal{M}}) = \mathbb{R}^{n|\mathcal{V}|}$ (see Section 1.1.2), and (106) is inferred from (104). Since $\text{rank}(B) = p|\mathcal{V}|$, using [73, Fact 2.15.2], equation (107) is derived from (105). $\square$

Note that if items 1, and 2 in Theorem 3 hold (i.e., $B_i = \bar{B}$, $C_i = \bar{C}$ for all $i \in \mathcal{V}$, and $P'B = C^\top$, for some matrices $\bar{B}$, $\bar{C}$ and $P' = P'^\top > 0$), then the condition $\text{rank}(B) = p|\mathcal{V}|$ in Lemma 5 is equivalent to the crucial condition that $C_iB_i$ is similar to a diagonal matrix with positive entries for all $i \in \mathcal{V}$, as stated in [36, 92, 93]. Let us recall that:

$$ \mathcal{M} := G^\top \circ \Phi' \circ G, \qquad (108) $$

where $G = N'\Delta^{-1}(\Theta^\top \otimes I_p)CR^{-1} \in \mathbb{R}^{m'|\mathcal{E}| \times n'|\mathcal{V}|}$, and $\Phi'$ is defined in (39). Recall that $z_k = Rx_k$ where $x_k := x(t_k)$ is the state of all agents at time $t = t_k$ defined in (9). In what follows, we compute the controller $u_k$ based on the given state $x_k$ (resp. on the state $z_k = Rx_k$) so that the network achieves global exponential synchronization, using two approaches:

1. Computing $\mathcal{J}_{h\mathcal{M}}$ by solving an optimization problem (OP) when $\Gamma' = (\Gamma')^\top$.

2. Solving an LCP (applied in the case where $\Gamma'$ is not necessarily symmetric)

## 6.1 Controller Design Based on Resolvent Computation

In this subsection, the assumption $\Gamma' = (\Gamma')^\top$ is assumed to hold. The following lemma provides a method for computing the resolvent $\mathcal{J}_{h\mathcal{M}}$ where $\mathcal{M}$ is defined in (108). Consequently, the controller $u_k$ will be derived from Lemma 5, and an algorithm is provided to compute $u_k$.

**Lemma 6.** Assume that Assumptions 2, 4, 6 hold, and $\Gamma' = (\Gamma')^\top$. Then, for each $z \in \mathbb{R}^{n|\mathcal{V}|}$:

$$ \mathcal{J}_{h\mathcal{M}}(z) = z + hG^\top \arg\min_{x \in \mathbb{R}_+^{n|\mathcal{V}|}} \left( \frac{1}{2h} x^\top \Gamma' x + \frac{1}{h} x^\top r' + \frac{1}{2} \|G^\top x + \frac{z}{h}\|^2 \right), $$

where $\mathcal{M}$, $G$ and $r'$ are defined in (108) and in (34), respectively.

*Proof.* Since $\Gamma' = (\Gamma')^\top$ and Assumption 2 holds, using [94, Section 3.2.1, p.6], $\Phi' = \partial g^*$, where $g^*$ is the conjugate function of the proper lower semicontinuous convex function $g: \mathbb{R}^{m'|\mathcal{E}|} \to \mathbb{R} \cup \{+\infty\}$ defined by:

$$ g(x) = \sup_{v \in \mathbb{R}_+^{m'|\mathcal{E}|}} \langle v, x \rangle + \frac{1}{2} x^\top \Gamma' x - x^\top r', \quad \forall x \in \mathbb{R}^{m'|\mathcal{E}|}, $$

and $\text{Dom}\,\Phi' \subseteq \text{dom}(g^*) \subseteq \text{cl}(\text{Dom}\,\Phi')$. In addition, it follows from Remark 7 that $\text{Dom}\,\Phi' = (\text{SOL}(0, \Gamma'))^* - r'$ is a closed convex set. This implies that $\text{dom}(g^*) = (\text{SOL}(0, \Gamma'))^* - r'$. Using [57, Theorem 13.37, Proposition 13.13], it can be checked that
$g^*$ is also a proper lower semicontinuous convex function. Since Assumption 4 holds, $\mathcal{M} = G^\top \circ \partial g^* \circ G = \partial(g^* \circ G)$, where $g^* \circ G$ is a proper lower semicontinuous convex function, using [57, Corollary 16.53, Proposition 6.19 - (x)]. Consequently, $\mathcal{J}_{h\mathcal{M}} = \text{prox}_{h(g^* \circ G)}$. Using item 3 presented immediately prior to Remark 1, the following hold for all $z \in \mathbb{R}^{n|\mathcal{V}|}$:

$$ \mathcal{J}_{h\mathcal{M}}(z) = \text{prox}_{h(g^* \circ G)}(z) = z - h\text{prox}_{\frac{1}{h}(G^* \circ G)^*}\left(\frac{z}{h}\right) = z - h\text{prox}_{\frac{1}{h}G^\top \triangleright g}\left(\frac{z}{h}\right), \quad (109) $$

where $(G^\top \triangleright g)(y) := \inf\{g(x) \mid G^\top x = y\}$ with $\text{dom}(G^\top \triangleright g) = \text{rg}(G^\top)$, and the last equality follows from [57, Proposition 13.24, Theorem 13.37]. In addition, it can be checked that for all $y \in \mathbb{R}^{n|\mathcal{V}|}$:

$$ \begin{aligned} (G^\top \triangleright g)(y) &= \inf\{g(x) \mid G^\top x = y\} = \inf_{G^\top x = y} \left(\frac{1}{2}x^\top \Gamma'x - x^\top r'\right) \\ &= \inf_{G^\top x = -y} \left(\frac{1}{2}x^\top \Gamma'x + x^\top r'\right). \end{aligned} $$

Since Assumptions 2, and 6 hold, it follows from item 1.1 in Proposition 7 that $r' \in (\text{SOL}(0, \Gamma'))^*$. Therefore, LCP($r', \Gamma'$) has a solution using [56, Theorem 3.8.6, Exercise 3.12.1]. In addition, since $\Gamma' = (\Gamma')^\top \succcurlyeq 0$, the following quadratic program:

$$ \min_{x \ge 0} \left( \frac{1}{2} x^\top \Gamma' x + x^\top r' \right) \qquad (110) $$

has a finite optimal value, see [56, p.5] or Appendix C.1. This implies that for a fixed $y \in \text{rg}(G^\top)$, the function $x \mapsto \frac{1}{2}x^\top \Gamma'x + x^\top r'$ is bounded below on the set $\{x \in \mathbb{R}^{m'|\mathcal{E}|} : x \ge 0, G^\top x = -y\}$. Using [56, Theorem 2.8.1], it can be verified that:

$$ \inf_{G^\top x = -y} \left( \frac{1}{2} x^\top \Gamma' x + x^\top r' \right) = \min_{G^\top x = -y} \left( \frac{1}{2} x^\top \Gamma' x + x^\top r' \right), \quad \forall y \in \text{rg}(G^\top). \quad (111) $$

Using item 4 presented immediately prior to Remark 1 and (111), the following hold:

$$ \begin{aligned} \text{prox}_{\frac{1}{h}G^\top \triangleright g}\left(\frac{z}{h}\right) &= \arg\min_{y \in \mathbb{R}^{n|\mathcal{V}|}} \left(\frac{1}{h}(G^\top \triangleright g)(y) + \frac{1}{2}\|y - \frac{z}{h}\|^2\right) \\ &= \arg\min_{y \in \text{rg}(G^\top)} \left( \min_{G^\top x = -y} \left( \frac{1}{2h}x^\top \Gamma'x + \frac{1}{h}x^\top r' \right) + \frac{1}{2}\|y - \frac{z}{h}\|^2 \right) \\ &= \arg\min_{y \in \text{rg}(G^\top)} \min_{G^\top x = -y} \Upsilon(x) = \arg\min_{y \in \text{rg}(G^\top)} H(y), \end{aligned} \quad (112) $$

where $\Upsilon(x) := \frac{1}{2h}x^\top \Gamma'x + \frac{1}{h}x^\top r' + \frac{1}{2}\|G^\top x + \frac{z}{h}\|^2$, and $H(y) := \min_{G^\top x = -y} \Upsilon(x)$. Analogously to the explanation given just before (111), it can be checked that the optimization problem $\min_{x \ge 0} \Upsilon(x)$ admits a minimizer. Let $x^* = \arg\min_{x \ge 0} \Upsilon(x)$, then it always holds that $H(y) \ge \Upsilon(x^*)$, for all $y \in \text{rg}(G^\top)$. On the other hand, we have $H(-G^\top x^*) \le \Upsilon(x^*)$ since $x^* \in \{x : x \ge 0, G^\top x = -(-G^\top x^*)\}$. Therefore $H(-G^\top x^*) = \Upsilon(x^*) \le H(y)$ for all $y \in \text{rg}(G^\top)$. This means that $\arg\min_{y \in \text{rg}(G^\top)} H(y) = -G^\top x^*$. The conclusion is derived from (109), (112). $\square$
The following proposition is an immediate consequence of Lemmas 5, and 6.

**Proposition 14.** Suppose that all conditions in Lemmas 5, and 6 hold, then $u_k$ in (107) is rewritten as:

$$u_k = -\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)CR^{-1}z_k + \Xi^{-1}(\Theta \otimes I_p)H'\kappa(\bar{z}_k),$$

where $\kappa(\bar{z}_k) := \arg \min_{x \in \mathbb{R}_+^{n|V|}} \left( \frac{1}{2h}x^{\top}\Gamma'x + \frac{1}{h}x^{\top}r' + \frac{1}{2}\|G^{\top}x + \frac{\bar{z}_k}{h}\|^2 \right)$, and $\bar{z}_k$ is defined by:

$$\bar{z}_k := z_k + hRF_{nom}(t_k, R^{-1}z_k) - hRB\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)CR^{-1}z_k.$$

The relationship of $\kappa(\bar{z}_k)$ with the complementarity variable $\lambda'_{k+1}$ of the network is shown in Remark 13. Recall that $z_k = Rx_k$. Algorithm 2 allows us to compute the control input $u_k$ for a given $x_k$ so that the network achieves global exponential synchronization, via the solution of an optimization problem $\kappa(\bar{z}_k)$ defined in Proposition 14, which can be found using CVXPY Tool in Python, see Appendix C.1. This algorithm is employed in Examples 7, 8.

**Algorithm 2** Compute $u_k$ using the Resolvent Approach

**Require:** State $x_k$, time step $h$

**Ensure:** $u_k$

$$1: \quad \bar{z}_k \leftarrow (Rx_k + hRF_{nom}(t_k, x_k) - hRB\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)Cx_k)$$

$$2: \quad \kappa(\bar{z}_k) \leftarrow \arg \min_{x \ge 0} \left( \frac{1}{2h}x^{\top}\Gamma'x + \frac{1}{h}x^{\top}r' + \frac{1}{2}\|G^{\top}x + \frac{\bar{z}_k}{h}\|^2 \right)$$

$$3: \quad u_k \leftarrow -\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)Cx_k + \Xi^{-1}(\Theta \otimes I_p)H'\kappa(\bar{z}_k)$$

$$4: \quad \textbf{return} \ u_k$$

## 6.2 Solving an LCP for Controller Calculation

Recall from Lemma 5 that for a given state $z_k$, the controller $u_k$ can be computed via the result of $\tilde{z}_{k+1}$, which is obtained by computing the resolvent $J_{h\mathcal{M}}$. This approach was thoroughly analyzed in Section 5.1 under the assumptions that $\Gamma' = (\Gamma')^{\top}$, and Assumption 6 holds. In this subsection, we use the same computation of $u_k$ as given in Lemma 5, but with a different method for computing $\tilde{z}_{k+1}$. By this approach, $\Gamma'$ is not necessarily assumed symmetric. Throughout the rest of this subsection, all conditions in Lemma 5 are assumed to hold. Since $\mathcal{M} = G^{\top} \circ \mathbf{\Phi}' \circ G$ (see (108)), then the inclusion (103) is rewritten as:

$$\tilde{z}_{k+1} \in z_k + hRF_{nom}(t_k, R^{-1}z_k) - hRB\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)CR^{-1}z_k - hG^{\top}\mathbf{\Phi}'(G\tilde{z}_{k+1}). \quad (113)$$

Proceeding similarly to the proof of Lemma 4, the inclusion (113) is rewritten as:

$$\left\{ \begin{array}{l} \tilde{z}_{k+1} = z_k + hR(F_{nom}(t_k, R^{-1}z_k) - hRB\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)CR^{-1}z_k + hG^{\top}\lambda'_{k+1}, \\ w'_{k+1} = G\tilde{z}_{k+1} + \Gamma'\lambda'_{k+1} + r', \\ 0 \le \lambda'_{k+1} \perp w'_{k+1} \ge 0, \end{array} \right. \quad (114)$$
where $\lambda'_{k+1}, w'_{k+1} \in \mathbb{R}^{m'|\epsilon|}$ are complementarity variables. After some manipulations, (114) can be rewritten as:

$$ \left\{ \begin{array}{l} \tilde{z}_{k+1} = \bar{z}_k + hG^\top \lambda'_{k+1}, \\ w'_{k+1} = (\Gamma' + hGG^\top)\lambda'_{k+1} + G\bar{z}_k + r', \\ 0 \le \lambda'_{k+1} \perp w'_{k+1} \ge 0, \end{array} \right. \qquad (115) $$

where $\bar{z}_k := z_k + hR(F_{nom}(t_k, R^{-1}z_k) - hRB\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)CR^{-1}z_k$. The following proposition is a direct consequence of Lemma 5 and (115).

**Proposition 15.** Suppose that all conditions in Lemma 5 hold, then $u_k$ in (107) is rewritten as:

$$ u_k = -\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)CR^{-1}z_k + \Xi^{-1}(\Theta \otimes I_p)H'\lambda'_{k+1}, $$

where $\lambda'_{k+1} = \text{SOL}(G\bar{z}_k + r', \Gamma' + hGG^\top)$, and $\bar{z}_k$ is defined by:

$$ \bar{z}_k := z_k + hRF_{nom}(t_k, R^{-1}z_k) - hRB\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)CR^{-1}z_k. $$

**Remark 13.** Proposition 15 is an extension of Proposition 14. Indeed, when $\Gamma' = \Gamma^\top$, the value of $\lambda'_{k+1}$ in Proposition 15 is exactly equal to the value of $\kappa(\bar{z}_k)$ in Proposition 14, using the equivalence of an LCP and an optimization problem in the case the LCP matrix is symmetric and positive semidefinite (see Appendix C.1).

Recall that $z_k = Rx_k$. Algorithm 3 allows us to compute the control input $u_k$ for a given $x_k$ so that the network achieves global exponential synchronization, via the solution $\lambda'_{k+1}$ of an LCP defined in Proposition 15, which can be found using LCP Solvers in Python, see Appendix C.1. This algorithm is employed in Example 6.

**Algorithm 3** Compute $u_k$ using the LCP Approach

**Require:** State $x_k$, time step h

**Ensure:** $u_k$

1: $\bar{z}_k \leftarrow (Rx_k + hRF_{nom}(t_k, x_k) - hRB\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)Cx_k)$

2: $\lambda'_{k+1} \leftarrow \text{SOL}(G\bar{z}_k + r', \Gamma' + hGG^\top)$

3: $u_k \leftarrow -\Xi^{-1}(\Theta \otimes I_p)W(\Theta^\top \otimes I_p)Cx_k + \Xi^{-1}(\Theta \otimes I_p)H'\lambda'_{k+1}$

4: **return** $u_k$

## 6.3 Examples and Numerical Simulations

**Example 6.** Let us consider Example 2. In this framework, unknown perturbations are taken into account, and practical implementation is carried out. Numerical simulations of the control inputs, the states, the unknown perturbations, the outputs in the network with time step $h = 10^{-2}s$, and the log-scaled norms of error states with time steps $h = 10^{-2}s, 10^{-3}s, 10^{-4}s$, are depicted in Figures 24, 25, 27, 26. They utilize Algorithm 3 with LCP Solver given in Appendix C.1 to solve LCPs, and Semi-Backward
Euler scheme (101) with initial state $x(0) = [4 \ 1 \ -2]^\top$ and the parameters of $f, \xi_i$ as given in (89). The numerical simulations are implemented on Google Colaboratory using Python 3 (see details in Appendices C.1, C.2). The simulation scripts can be found in Appendix C.2.

Figure 24 shows that the control inputs $u_i(t)$ designed in this example do not exhibit the hazardous chattering phenomenon as mentioned in [50]. As illustrated in Figure 26, the error state $e(t)$ does not converge to zero but remains bounded by a threshold that depends on the time step $h$, and is observed to satisfy $\|e(t)\| \le h$, for all $t \ge 2$. As $h$ decreases, the norm of the error state $e(t)$ also decreases proportionally. This implies that improving the synchronization performance (i.e., reducing the magnitude of $e(t)$) requires the control signals must be transmitted at increasingly shorter time intervals. Furthermore, the difference in the norm of the error state $e(t)$ between the ideal case with known perturbations shown in Figure 7 (Example 2) and the realistic case with unknown perturbations depicted in Figure 26 can be explained by the interpretation of $\xi_i(t, x_i)$ given in (82), as discussed in the introduction of Section 5. Specifically, in Example 2, $\xi_i(t, x_i)$ was treated as a part of the vector field $f_i$, and Figure 7 merely illustrates that a network composed of agents with distinct vector fields (i.e., heterogeneous) can achieve exponential synchronization. In contrast, in this example, $\xi_i(t, x_i)$ acts as unknown perturbations that directly affects the state $x(t)$ at each time, and the discrete-time controller $u(t)$ is designed to compensate for these perturbations and to ensure that the network achieves practical synchronization in discrete-time.

<figure><img src="image_25.png" /><figcaption>Figure 24: Numerical simulations of the control inputs $u_i(t)$ in Example 6, using Algorithm 3 and Semi-Backward Euler scheme (101) with time step $h = 10^{-2}s$, initial state $x(0) = [4 \ 1 \ -2]^\top$, and the parameters of $f, \xi_i$ given in (89). LCP Solver given in Appendix C.1 is used to solve LCPs in Algorithm 3.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure 25: Numerical simulations of the states $x_i(t)$ (left figure) and the perturbations $\xi_i(t, x_i)$ (right figure) in Example 6, using Algorithm 3 and Semi-Backward Euler scheme (101) with time step $h = 10^{-2}s$, initial state $x(0) = [4 \ 1 \ -2]^T$, and the parameters of $f, \xi_i$ given in (89). LCP Solver given in Appendix C.1 is used to solve LCPs in Algorithm 3.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Figure 26: Numerical simulations of the log-scaled norms of the error states $e_{12}(t)$, $e_{13}(t)$, $e_{23}(t)$ and $e(t)$ in Example 6 where $e(t) = (\Theta_c^\top \otimes I_1)x(t)$, using Algorithm 3 and Semi-Backward Euler scheme (101) with time steps $h = 10^{-2}s$ (top figure), $h = 10^{-3}s$ (middle figure), $h = 10^{-4}s$ (bottom figure), initial state $x(0) = [4 \ 1 \ -2]^T$, and the parameters of $f, \xi_i$ given in (89). LCP Solver given in Appendix C.1 is used to solve LCPs in Algorithm 3.</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Figure 27: Numerical simulations of the outputs $y_i(t)$ in Example 6, using Algorithm 3 and Semi-Backward Euler scheme (101) with time step $h = 10^{-2}s$, initial state $x(0) = [41-2]^T$, and the parameters of $f, \xi_i$ given in (89). LCP Solver given in Appendix C.1 is used to solve LCPs in Algorithm 3.</figcaption></figure>

**Example 7.** Let us consider Example 4. In this framework, unknown perturbations are taken into account, and practical implementation is carried out. Numerical simulations of the control inputs, the states, the unknown perturbations in the network with time step $h = 10^{-2}s$, and the log-scaled norms of error states with time steps $h = 10^{-2}s, 10^{-3}s$ are depicted in Figures 28, 29, 30. Algorithm 2 with CVXPY Tool given in Appendix C.1 to solve OPs, and Semi-Backward Euler scheme (101) are used. The initial state $x(0)$, and the parameters are given in Appendix C.2. The numerical simulations are implemented on Google Colaboratory using Python 3 (see details in Appendices C.1, C.2). The simulation scripts can be found in Appendix C.2.

Let $e_{ij,1}(t) := x_{i,1}(t) - x_{j,1}(t)$, and $e_{ij,2}(t) := x_{i,2}(t) - x_{j,2}(t)$, for all $t \ge 0$, for all $i \ne j$, $i, j = 1, \dots, 30$. These are the entries of the error state $e(t)$ defined in (20). As illustrated in Figure 30, the error states $e(t)$, $e_{ij,1}(t)$, and $e_{ij,2}(t)$ in this example do not converge to zero but remain bounded by thresholds that depend on the time step $h$, decrease as $h$ decreases, and are observed to approximately satisfy that for all $i \ne j$, $i, j = 1, \dots, 30$:

$$ \begin{cases} |e_{ij,1}(t)| \le 3h, \quad \forall t \ge 1. \\ |e_{ij,2}(t)| \le 0.5h, \quad \|e(t)\| \le 20h, \quad \forall t \ge 10. \end{cases} $$

The difference in the norm of the error state $e(t)$ between the ideal case with known perturbations (shown in Figure 16 in Example 4) and the realistic case with unknown perturbations (depicted in the top two figure of Figure 30), can be explained similarly to Example 6 via the role of $\bar{\xi}_{i,1}(t, x_i)$, $\bar{\xi}_{i,2}(t, x_i)$ given in (93).

Figure 28 shows that the discrete-time controllers $u_i(t)$ designed in this example do not exhibit the hazardous chattering phenomenon as reported in [50]. This is crucial for practical applications. At the beginning of Section 5.1 and Section 6, it is recalled that discretization methods significantly influence the chattering phenomenon in the design of set-valued feedback controllers, where Semi-Backward Euler scheme effectively suppresses such phenomenon, while Forward Euler scheme fails to do so. This fact is
once again confirmed in this example. Indeed, the time-discretization of the network dynamics (100) using Forward Euler scheme, is given by:

$$ \begin{cases} z_{k+1} = z_k + h(RF_{nom}(t_k, R^{-1}z_k) + RB\xi_{per}(t_k, R^{-1}z_k) + RBu_k), \\ u_k \in -\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)CR^{-1}z_k - \tilde{\mathcal{M}}(CR^{-1}z_k), \end{cases} $$

where $\tilde{\mathcal{M}}$ is defined in Lemma 4, $z_k = Rx_k$ with $R = R^{\top} > 0$ a square root of $P'$ given in Assumption 3', and $x_k$ the state of the network at time $t = t_k$. Proceeding similarly to the proof of Lemma 4, it can be rewritten as:

$$ \left\{ \begin{aligned} z_{k+1} &= z_k + h(RF_{nom}(t_k, R^{-1}z_k) + RB\xi_{per}(t_k, R^{-1}z_k) + RBu_k), \\ u_k &= -\Xi^{-1}(\Theta \otimes I_p)W(\Theta^{\top} \otimes I_p)CR^{-1}z_k + \Xi^{-1}(\Theta \otimes I_p)H'\lambda'_k, \\ 0 &\le \lambda'_k \perp w'_k = Gz_k + \Gamma'\lambda'_k + r' \ge 0, \end{aligned} \right. \quad (116) $$

where $\lambda'_k, w'_k \in \mathbb{R}^{m'|\mathcal{E}|}$ are complementarity variables. Under the same initial state $x(0)$ and parameters (given in Appendix C.2), the numerical simulations of the control inputs $u_i(t)$, the states $x_i(t)$, and the log-scaled norm of the error state $e(t)$ in the network when using Forward Euler scheme (116) with time step $h = 10^{-2}s$, are depicted in Figure 31. This figure shows that the discrete-time controllers $u_i(t)$ exhibit chattering, in contrast to the controllers $u_i(t)$ in Figure 28 where Semi-Backward Euler scheme is used for time-discretization of (100). Consequently, the network cannot achieve synchronization because $u_i(t)$ is inaccurately computed. These analyses clarify why Semi-Backward Euler scheme (101) is selected for the study of practical implementation in Section 6.

<figure><img src="image_29.png" /><figcaption>Figure 28: Numerical simulations of the control inputs $u_i(t)$ in Example 7 in different coordinate regions, using Algorithm 2 and Semi-Backward Euler scheme (101) with time step $h = 10^{-2}s$. The initial state $x(0)$, and the parameters are given in Appendix C.2. CVXPY Tool given in Appendix C.1 is used to solve OPs in Algorithm 2.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure 29: Numerical simulations of the states $x_{i,1}(t)$, $x_{i,2}(t)$ (left figure) and the perturbations $\xi_i(t, x_i)$ given in (93) (right figure) in Example 7. Algorithm 2 and Semi-Backward Euler scheme (101) are used with time step $h = 10^{-2}s$. The initial state $x(0)$, and the parameters are given in Appendix C.2. CVXPY Tool given in Appendix C.1 is used to solve OPs in Algorithm 2.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure 30: Numerical simulations of the log-scaled norms of the error states $e(t)$, $e_{ij,1}(t)$ and $e_{ij,2}(t)$ in Example 7 where $e(t) = (\Theta_c^\top \otimes I_2)x(t)$, $e_{ij,1}(t) = x_{i,1}(t) - x_{j,1}(t)$, and $e_{ij,2}(t) = x_{i,2}(t) - x_{j,2}(t)$. Three figures on the left using the time step $h = 10^{-2}s$, and the rest using the time step $h = 10^{-3}s$. Algorithm 2 with CVXPY Tool given in Appendix C.1 to solve OPs, and Semi-Backward Euler scheme (101) are used. The initial state $x(0)$, and the parameters are given in Appendix C.2.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Figure 31: Numerical simulations of the control inputs $u_i(t)$, the states $x_{i,1}(t), x_{i,2}(t)$ and the log-scaled norm of the error state $e(t)$ in Example 7 where $e(t) = (\Theta_c^\top \otimes I_1)x(t)$. Forward (fully-explicit) Euler scheme (116) with LCP Solver given in Appendix C.1 to solve LCPs, is used with time step $h = 10^{-2}s$. The initial state $x(0)$, and the parameters are given in Appendix C.2.</figcaption></figure>

**Example 8.** Let us consider Example 5. In this framework, unknown perturbations are taken into account, and practical implementation is carried out. Numerical simulations of the control inputs, the states in the network, and the log-scaled norms of error states are depicted in Figures 32, 33, using Algorithm 2 with CVXPY Tool given in Appendix C.1 to solve OPs, and Semi-Backward Euler scheme (101) with time step $h = 4 \times 10^{-4}s$. They are implemented on Google Colaboratory using Python 3 (see details in Appendices C.1, C.2). The simulation scripts can be found in Appendix C.2.

Figure 32 shows that the discrete-time controllers $u_i(t)$ designed in this example do not exhibit the hazardous chattering phenomenon as reported in [50]. This is crucial for practical applications. Let $e_{ij,1}(t) := x_{i,1}(t) - x_{j,1}(t)$, $e_{ij,2}(t) := x_{i,2}(t) - x_{j,2}(t)$, and $e_{ij,3}(t) := x_{i,3}(t) - x_{j,3}(t)$, for all $t \ge 0$, for all $i \ne j$, $i, j = 1, \dots, 10$, where $x_i(t) = [x_{i,1}(t) \ x_{i,2}(t) \ x_{i,3}(t)]^\top$ is the state of agent i in the network. These are the entries of the error state $e(t) = (\Theta_c^\top \otimes I_3)x(t)$ where $x(t) = [x_1(t) \ \dots \ x_{10}(t)]^\top$ is the state of the network. As illustrated in Figure 33 (four figures on the right), the error states $e(t)$, $e_{ij,1}(t)$, $e_{ij,2}(t)$, $e_{ij,3}(t)$ in this example do not converge to zero but remain bounded by thresholds that depend on the time step $h$, decrease as $h$ decreases, and are observed to approximately satisfy that for all $i \ne j$, $i, j = 1, \dots, 10$:

$$|e_{ij,1}(t)| \le 5h, \quad |e_{ij,2}(t)| \le 0.5h, \quad |e_{ij,3}(t)| \le 10h, \quad \|e(t)\| \le 25h, \quad \forall t \ge 1.$$

The difference in the norm of the error state $e(t)$ between the ideal case with known perturbations (shown in Figure 23 in Example 5) and the realistic case with unknown perturbations (depicted in the top right figure in Figure 33), can be explained similarly
to Example 6 via the role of $\xi_i$ defined after (95). Figure 33 show that the network achieves practical synchronization in discrete-time in a very short time, and the state after synchronization exhibits Double Scroll Attractor behavior, which is a type of complex chaotic behavior.

Since the state of the network exhibits complex chaotic behavior, they change rapidly, complexly, and are highly sensitive to small differences in the initial conditions. Therefore, if time step $h$ is not sufficiently small, the simulation results will no longer be accurate since the network's variations can be missed, and the state of the network may diverge substantially from the true trajectory. This is illustrated in Figure 34. This figure shown that, with the same data, Algorithm 2, and the Semi-Backward Euler scheme (101), increasing the time step $h$ from $4 \times 10^{-4}$ to $5 \times 10^{-4}$, although the network remains practically synchronized in discrete-time, the state $x(t)$ has diverged and no longer exhibits Double Scroll Attractor behavior. This consequently causes the discrete-time controllers $u_i(t)$ to exhibit hazardous chattering phenomenon.

<figure><img src="image_33.png" /><figcaption>Figure 32: Numerical simulations of the control inputs $u_i(t) = [u_{i,1}(t) \ u_{i,2}(t) \ u_{i,3}(t)]^\top$ in Example 8 in different coordinate regions, using Algorithm 2 with CVXPY Tool given in Appendix C.1 to solve OPs, and Semi-Backward Euler scheme (101) with time step $h = 4 \times 10^{-4}$s. The initial state $x(0)$, and the parameters are given in Appendix C.2.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>Figure 33: Numerical simulations of the states $x_i(t) = [x_{i,1}(t) \ x_{i,2}(t) \ x_{i,3}(t)]^\top$ (four figures on the left) and the log-scaled norms of the error states $e(t), e_{ij,1}(t), e_{ij,2}(t)$ and $e_{ij,3}(t)$ in Example 5 (four figures on the right), where $e(t) = (\Theta_c^\top \otimes I_3)x(t), e_{ij,1}(t) = x_{i,1}(t) - x_{j,1}(t), e_{ij,2}(t) = x_{i,2}(t) - x_{j,2}(t), \text{ and } e_{ij,3}(t) = x_{i,3}(t) - x_{j,3}(t)$. Algorithm 2 with CVXPY Tool given in Appendix C.1 to solve OPs, and Semi-Backward Euler scheme (101) are used with time step $h = 4 \times 10^{-4}$s. The initial state $x(0)$, and the parameters are given in Appendix C.2.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>Figure 34: Numerical simulations of the control inputs $u_i(t) = [u_{i,1}(t) \ u_{i,2}(t) \ u_{i,3}(t)]^\top$ and the states $x_i(t)$ in Example 8, using Algorithm 2 with CVXPY Tool given in Appendix C.1 to solve OPs, and Semi-Backward Euler scheme (101) with time step $h = 5 \times 10^{-4}$s. The initial state $x(0)$, and the parameters are given in Appendix C.2.</figcaption></figure>

**Remark 14.** In Examples 4, 5, the quantities $\bar{\xi}_{i,2}$, $\xi_{i,j,2}$ are introduced as input perturbations affecting agent i in practice, and are subsequently used in Examples 7, 8 for practical implementation, respectively. They are modeled in (93) and (96) as functions solely of time and the local state $x_i$ of agent i. However, these input perturbations can be generalized to depend on time and the entire network state $x = [x_1^\top \dots x_{|V|}^\top]^\top$, thus providing a more realistic representation. For instance, the following generalized input perturbation can be considered in Examples 4, 7:

$$ \bar{\xi}_{i,2}(t, x) = 0.3 \cos(-2t + 0.2x_{i,2}) \sum_{j \in \mathcal{N}_i} \sin(\gamma_{i,j}x_{j,1} + x_{j,2} - 3t^2), \text{ for some } \gamma_{i,j} \in \mathbb{R}, $$

where $\mathcal{N}_i \subseteq \{1, \dots, 30\}$ is the set of neighbours of agent i in the network.

## 6.4 Evaluation of Computational Efficiency: Algorithms 2, 3, and LCP Solvers when $\Gamma' = (\Gamma')^\top$

In this section, we evaluate the execution time required to perform the simulations of Examples 6, 7, and 8 in practical implementation. Since the matrix $\Gamma'$ in these examples is symmetric, they are well-suited for comparing the computational performance of Algorithms 2, and 3 in terms of execution time. Note that Algorithm 2 employs the CVXPY Tool given in Appendix C.1 to solve OPs, while Algorithm 3 utilizes either LCP solver or the CVXPY Tool, both given in Appendix C.1, to solve LCPs. In what follows, the comparisons will be based on the following two criteria:

1. Resolvent approach (Algorithm 2) vs. LCP approach (Algorithm 3)

2. LCP solver used in Algorithm 3 (CVXPY Tool vs. LCP SOLVER)
All numerical simulations of Examples 6, 7, 8, are implemented on Google Collaboratory using Python version 3.11.12, where computations take place on Google's servers (specifically, on a virtual machine within the Google Cloud infrastructure). Table 5 provides the computation time required to obtain the numerical simulation of the network over the time interval [0, 20] seconds with different time steps $h$. Therefore, the number of steps (iterations) in the simulation process varies for each case of the time step $h$. For example, when $h = 10^{-2}$ is considered, it means that the number of iterations performed is 2000. Some observations from Table 5 are as follows:

* Resolvent Approach has a shorter computation time than LCP Approach in Example 7, but a longer one in Example 6, and an approximately equal one in Example 8.

* In LCP Approach, using CVXPY Tool to solve LCPs yields a shorter execution time than using LCP SOLVER in Examples 7, and 8, but a longer one in Example 6.

Therefore, we can conclude that the choice between Resolvent Approach and LCP Approach should depend on each specific case in order to make a decision. A similar remark can be made regarding the selection of the LCP solver in LCP Approach.

<table><caption>Table 5: Computational times for both approaches on the time interval [0, 20] second.</caption><thead><tr><th rowspan="2">Cases</th><th rowspan="2"></th><th rowspan="2">Results<br>(CVXPY Tool)</th><th colspan="2">LCP Approach</th></tr><tr><th>CVXPY Tool</th><th>LCP SOLVER</th></tr></thead><tbody><tr><td rowspan="2">Example 6</td><td>$h = 10^{-2}$ (s)</td><td>5.51 (s)</td><td>4.92 (s)</td><td>0.91 (s)</td></tr><tr><td>$h = 10^{-3}$ (s)</td><td>58.7 (s)</td><td>54.8 (s)</td><td>10.1 (s)</td></tr><tr><td rowspan="2">Example 7</td><td>$h = 10^{-2}$ (s)</td><td>6.83 (s)</td><td>14.8 (s)</td><td>395.9 (s)</td></tr><tr><td>$h = 10^{-3}$ (s)</td><td>74.2 (s)</td><td>153.8 (s)</td><td>3555 (s)</td></tr><tr><td rowspan="2">Example 8</td><td>$h = 5 \times 10^{-4}$ (s)</td><td>128.6 (s)</td><td>128.4 (s)</td><td>388.4 (s)</td></tr><tr><td>$h = 4 \times 10^{-4}$ (s)</td><td>166.3 (s)</td><td>165.8 (s)</td><td>482.6 (s)</td></tr></tbody></table>

Since the time step $h = 10^{-2}$s is used in Examples 6, 7, and $h = 4 \times 10^{-4}$s is used in Example 8, the number of steps (iterations) performed to simulate the control input $u(t)$ on the entire time interval [0, 20] seconds in Examples 6, 7, 8 are 2000, 2000, and 50000, respectively. The statistical summary of computation times for each discrete step over all executed steps is presented in Table 6, including metrics such as the maximum (except initial time), minimum, average, and 99% of data (i.e., the value below which 99% of the data fall). The values in this table suggest the approximate time required in practice for the computing unit to process state information and transmit the control signal $u_k$ back to the network, which guarantees the practical synchronization of the network in discrete-time. For instance, in Example 7, the computer requires on average approximately 0.0032 seconds from the moment it receives the state $x_k$ until it transmits the control signal $u_k$ back to the network. This also means that
numerical simulations of the network in Example 7 are meaningful only for time steps $h \ge 0.0032s$. The comparisons between the two approaches and the two LCP solvers in LCP Approach, based on the information given in Table 6, coincide with the observations from Table 5 presented above. In addition, it is important to note that execution times presented in Tables 5, and 6, may vary depending on various other factors, such as computer's capabilities, the optimization level of the code in each algorithm, or the internet connection quality if the computations are performed on Google Colaboratory.

**Conclusion:** It is inferred from this preliminary study, that both approaches are worth investigating and that future studies (in particular experimental ones) have to be led to determine which algorithm should be implemented.

<table><caption>Table 6: Statistical summary of computation time per discrete step of control input $u(t)$ on the time interval $[0, 20]$ seconds. The number of steps performed in Example 6 is 2000, in Example 7 is 2000, and in Example 8 is 50000. EIT: except initial time.</caption><thead><tr><th rowspan="2">Cases</th><th rowspan="2"></th><th rowspan="2">Results<br>(CVXPY Tool)</th><th colspan="2">LCP Approach</th></tr><tr><th>CVXPY Tool</th><th>LCP SOLVER</th></tr></thead><tbody><tr><td rowspan="4">Example 6</td><td>Max (EIT)</td><td>0.0097 (s)</td><td>0.0067 (s)</td><td>0.0016 (s)</td></tr><tr><td>Min</td><td>0.002</td><td>0.0019</td><td>0.0002</td></tr><tr><td>Average</td><td>0.0026</td><td>0.0023</td><td>0.0004</td></tr><tr><td>99% of Data</td><td>0.0045</td><td>0.0039</td><td>0.001</td></tr><tr><td rowspan="4">Example 7</td><td>Max (EIT)</td><td>0.0099</td><td>0.0224</td><td>0.434</td></tr><tr><td>Min</td><td>0.0026</td><td>0.0059</td><td>0.079</td></tr><tr><td>Average</td><td>0.0032</td><td>0.0071</td><td>0.198</td></tr><tr><td>99% of Data</td><td>0.0054</td><td>0.0122</td><td>0.353</td></tr><tr><td rowspan="4">Example 8</td><td>Max (EIT)</td><td>0.0339</td><td>0.0389</td><td>0.0694</td></tr><tr><td>Min</td><td>0.0021</td><td>0.0022</td><td>0.0054</td></tr><tr><td>Average</td><td>0.0028</td><td>0.0029</td><td>0.0092</td></tr><tr><td>99% of Data</td><td>0.0052</td><td>0.0062</td><td>0.0174</td></tr></tbody></table>

# 7 Various Extensions and Perspectives

## 7.1 All data of agents are time-varying

In this section, let us consider networks whose structure is represented by a connected, simple, undirected graph $G = (V, E)$ and their agents are time-varying dynamical systems of the form:

$$ \begin{cases} \dot{x}_i(t) = f_i(t, x_i(t)) + B_i(t)u_i(t), \\ y_i(t) = C_i(t)x_i(t) + D_i(t)u_i(t), \end{cases} \qquad (117) $$

where $x_i(t) \in \mathbb{R}^n$, $u_i(t), y_i(t) \in \mathbb{R}^p$ are the state, input and output of the agent $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(\cdot), C_i(\cdot), D_i(\cdot)$ are matrix functions with compatible dimensions. Let us propose the couplings formed by the time-varying linear complementarity problems (TVLCP) structured as the following:

$$
\left\{
\begin{aligned}
\eta_{ij}(t) &= H_{ij}(t)\lambda_{ij}(t), \quad \text{for } i < j, (i, j) \in \mathcal{E}, && (118a) \\
w_{ij}(t) &= M_{ij}(t)\lambda_{ij}(t) + N_{ij}(t)v_{ij}(t) + r_{ij}(t), && (118b) \\
0 &\le \lambda_{ij}(t) \perp w_{ij}(t) \ge 0, && (118c)
\end{aligned}
\right.
$$

where $\lambda_{ij}(t), w_{ij}(t) \in \mathbb{R}^m$ are complementarity variables, $H_{ij}(\cdot), M_{ij}(\cdot), N_{ij}(\cdot)$ and $r_{ij}(\cdot)$ are matrix functions that have compatible dimensions. The agents in (117) and couplings in (118) are extensions of the agents, and couplings given in (5), (6), respectively, where the coefficient matrices are considered to be time-varying. Using the interconnection law given in (7), and mimicking Section 1.3, the overall dynamical system of the entire network is defined by the following TVCS:

$$
\left\{
\begin{aligned}
\dot{x}(t) &= F(t, x(t)) + B(t)(\Theta \otimes I_p)H(t)\lambda(t), \\
w(t) &= N(t)(\Theta^{\top} \otimes I_p)C(t)x(t) + \Gamma(t)\lambda(t) + r(t), \\
0 &\le \lambda(t) \perp w(t) \ge 0,
\end{aligned}
\right.
\qquad (119)
$$

where $\Gamma(t) := N(t)(\Theta^{\top} \otimes I_p)D(t)(\Theta \otimes I_p)H(t) + M(t)$, $x(t) \in \mathbb{R}^{n|\mathcal{V}|}$, $\lambda(t)$, $w(t) \in \mathbb{R}^{m|\mathcal{E}|}$ are defined as in Section 1.3, $F(t, x(t))$ is defined in (8), $\Theta$ is an oriented incidence matrix associated with $\mathcal{G}$ satisfying (12), and $B(t)$, $C(t)$, $D(t)$, $H(t)$, $N(t)$, $M(t)$, $r(t)$ are defined in the same way as $B$, $C$, $D$, $H$, $N$, $M$, $r$ in Section 1.3, but note that the matrices are now replaced by matrix functions. Mimicking Lemma 1, the lemma below is achieved immediately.

**Lemma 7.** Consider the agents in (117) and the couplings in (118). Suppose that for each pair $(i, j) \in \mathcal{V} \times \mathcal{V}$ where $i < j$ and $(i, j) \in \mathcal{E}$, there exist matrix functions $W_{ij}(\cdot)$, $H'_{ij}(\cdot)$, $M'_{ij}(\cdot)$, $N'_{ij}(\cdot)$, and $r'(t)$ such that for all $t \ge 0$:

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

where $W_{ij}(t) \in \mathbb{R}^{p \times p}$, $H'_{ij}(t) \in \mathbb{R}^{p \times m'}$, $M'_{ij}(t) \in \mathbb{R}^{m' \times m'}$, $N'_{ij}(t) \in \mathbb{R}^{m' \times p}$, $r'_{ij}(t) \in \mathbb{R}^{m'}$ with $m' = m - 2p$. Then the interconnection law in (7) is rewritten as:

$$
\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) - W_{ij}(t)v_{ij}(t), \\
v_{ij}(t) = y_i(t) - y_j(t),
\end{array}
\right.
$$
where $\eta'_{ij}(t)$ is defined similarly to (118), but using the matrix functions $H'_{ij}(\cdot)$, $M'_{ij}(\cdot)$, $N'_{ij}(\cdot)$, $r'_{ij}(\cdot)$.

In this subsection, we assume that (120) holds. Let us define the following matrices:

$$
\begin{align*}
\Xi(t) &:= I_{p|\nu|} + (\Theta \otimes I_p)W(t)(\Theta^\top \otimes I_p)D(t), \\
\hat{D}(t) &:= (\Theta^\top \otimes I_p)D(t)(\Theta \otimes I_p), \\
\Delta(t) &:= I_{p|\mathcal{E}|} + \hat{D}(t)W(t), \\
\bar{\Delta}(t) &:= I_{p|\mathcal{E}|} + W(t)\hat{D}(t).
\end{align*}
$$

Note that if $D(t)$ and $W(t)$ are symmetric matrices then $\bar{\Delta}(t) = \Delta^\top(t)$. In addition, the following properties always hold:

$$
\Xi(t)^{-1}(\Theta \otimes I_p) = (\Theta \otimes I_p)\bar{\Delta}(t)^{-1}, \quad (121)
$$

$$
(\Theta^{\top} \otimes I_p) D(t) \Xi(t)^{-1} = \Delta(t)^{-1} (\Theta^{\top} \otimes I_p) D, \quad (122)
$$

using [84, Fact 3.20.6]. An extension of Assumption 1 is considered below:

**Assumption 7.** For all $t \ge 0$, the matrix $\Delta(t)$ is invertible.

Under the structure of couplings as in (120) and Assumption 7, mimicking Section 3.2 and using the properties (121) and (122), the TVCS (119) is rewritten as:

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

where $\tilde{F}(t, x(t)) := F(t, x(t)) - B(t)(\Theta \otimes I_p)\bar{\Delta}(t)^{-1}W(t)(\Theta^\top \otimes I_p)C(t)x(t)$, and

$$
\Gamma'(t) := N'(t)\Delta(t)^{-1}\hat{D}(t)H'(t) + M'(t),
$$

the complementarity variables $\lambda'(t), w'(t) \in \mathbb{R}^{m'|\mathcal{E}|}$ and the matrix functions $W(\cdot)$, $N'(\cdot)$, $H'(\cdot)$, $M'(\cdot)$, $r'(\cdot)$ are defined in the same way as $\lambda'(t)$, $w'(t)$, $W$, $N'$, $H'$, $M'$, $r'$ in Section 3.2, but note that the matrices are now replaced by matrix functions. Let us consider the following assumptions, which are direct extensions of Assumptions 2, 3, 4, respectively.

**Assumption 8.** For all $t \ge 0$, $\Gamma'(t) \ge 0$.

**Assumption 9.** There is a continuously differentiable $P'(t) = P'(t)^\top > 0$ such that for all $t \ge 0$:

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

**Assumption 10.** For all $t \ge 0$,

$$
\text{rge} \left( N'(t) \Delta(t)^{-1} (\Theta^\top \otimes I_p) C(t) \right) \cap \text{rint} \left( (\text{SOL}(0, \Gamma'(t)))^* - r'(t) \right) \neq \emptyset.
$$
Let $T$ be a positive number, and $\mathcal{H}(t) := (\text{SOL}(0, \Gamma'(t)))^*$ for all $t \in [0, T]$. Let us introduce the operator $\mathbf{\Phi}'_{tv}$, known as a time-varying Lur'e operator [69, 74, 94]:

$$ \mathbf{\Phi}'_{tv} : [0, T] \times \mathbb{R}^{m'|\mathcal{E}|} \to \mathbb{R}^{m'|\mathcal{E}|}, \quad (t, v) \mapsto \left( \mathcal{N}_{\mathbb{R}_{+}^{m'|\mathcal{E}|-r'(t)}}^{-1} + \Gamma'(t) \right)^{-1} (v). \quad (125) $$

The following theorem is derived from [69, Theorem 2], with a slight modification, guaranteeing the well-posedness of the TVCS (123) on the time interval $[0, T]$.

**Theorem 4.** Suppose that for each $i = 1, |V|$, $f_i : \mathbb{R}_+ \times \mathbb{R}^n \to \mathbb{R}^n$ is Lipschitz continuous in the second variable, and Assumptions 8, 9, 10 hold. For all $t \in [0, T]$, let $R(t)$ be a square root of $P'(t)$ defined in Assumption 9 such that $R(t) = R(t)^\top > 0$, and $G(t) := N'(t)\Delta(t)^{-1}(\Theta^\top \otimes I_p)C(t)R(t)^{-1}$. Assume further that the following hold:

1. There exists $\hat{Q} : [0, T] \to \mathbb{R}^{k \times m'|\mathcal{E}|}$ such that $\hat{Q}(\cdot)$ is bounded on $[0, T]$, and

(i) for each $t \in [0, T]$, SOL(0, $\Gamma'(t)$) = $\{\hat{Q}(t)^\top\beta$, with $\beta \in \mathbb{R}_+^k\}$, for some $k$,

(ii) $\mathcal{H}(s) \subseteq \mathcal{H}(t)$, for all $0 \le s \le t \le T$.

2. rge($G(t)$) is a constant linear space for all $t \in [0, T]$.

3. rge($G(s)$) $\cap (\mathcal{H}(s) - r'(t)) \neq \emptyset$, for all $0 \le s \le t \le T$.

4. $G(s)^{-1}(\mathcal{H}(t) - r'(t)) \subseteq G(t)^{-1}(\mathcal{H}(t) - r'(t))$, for all $0 \le s \le t \le T$.

5. Mappings $G(\cdot), \Gamma'(\cdot)$ are right-continuous, $\mathcal{H}(\cdot)$, $(\text{rge}(G(\cdot)) + \text{aff}(\mathcal{H}(\cdot)))^\perp$ are inner semi-continuous from right, and $r'(\cdot)$ is absolutely continuous on $[0, T]$.

6. All the matrix functions $B(\cdot), C(\cdot), D(\cdot), H'(\cdot), W(\cdot), N'(\cdot), M'(\cdot)$, and $\Delta(\cdot)^{-1}$ are bounded on $[0, T]$, and for each $i = 1, V$, for each fixed $\alpha \in \text{dom } f_i(t, \cdot)$, there exists a number $c_i(T) > 0$ such that:

$$ \|f_i(t, \alpha)\| \le c_i(T), \quad \forall t \in [0, T]. $$

Then (123) has a unique absolutely continuous solution $x(\cdot)$ on $[0, T]$ for any initial condition given by:

$$ x(0) = x_0 \in (N'(0)\Delta(0)^{-1}(\Theta^\top \otimes I_p)C(0))^{-1}(\mathcal{H}(0) - r'(0)). \quad (126) $$

*Proof.* Since Assumption 9 holds, using [69, Section 2.3.1-(b)], then the TVCS (123) is rewritten as:

$$ \left\{ \begin{array}{l} \dot{z}(t) = \dot{R}(t)R(t)^{-1}z(t) + R(t)\tilde{F}(t, R(t)^{-1}z(t)) + R(t)B(t)(\Theta \otimes I_p)\bar{\Delta}(t)^{-1}H'(t)\lambda'(t), \\ w'(t) = N'(t)\Delta(t)^{-1}(\Theta^\top \otimes I_p)C(t)R(t)^{-1}z(t) + \Gamma'(t)\lambda'(t) + r'(t), \\ 0 \le \lambda'(t) \perp w'(t) \ge 0, \end{array} \right. \quad (127) $$
where $z(t) = R(t)x(t)$ with $R(t) = R(t)^{\top} \succ 0$ be a square root of $P'(t)$ defined in Assumption 9. The TVCS (127) has the form of [69, System (47)], but the linear term $\bar{A}(t)z$ in [69, System (47)] is replaced by $\dot{R}(t)R(t)^{-1}z + R(t)\tilde{F}(t, R(t)^{-1}z)$. Let

$$ \lambda^{\text{im}}(t, R(t)^{-1}z) := \text{proj}_{\text{rge}}(\Gamma'(t)+\Gamma'(t)^{\top}) (\mathbf{\Phi}'_{tv}(t, N'(t)\Delta(t)^{-1}(\Theta^{\top} \otimes I_p)C(t)R(t)^{-1}z)) $$

where the operator $\mathbf{\Phi}'_{tv}$ is given in (125), and

$$ h(t, z) := \dot{R}(t)R(t)^{-1}z + R(t)\tilde{F}(t, R(t)^{-1}z) + \left( R(t)^{-1}C(t)^{\top}(\Theta \otimes I_p)(\Delta(t)^{-1})^{\top}N'(t)^{\top} - R(t)B(t)(\Theta \otimes I_p)\bar{\Delta}(t)^{-1}H'(t) \right) \lambda^{\text{im}}(t, R(t)^{-1}z). $$

In what follows, we apply [69, Theorem 2] to (127). Mimicking the proof of [69, Theorem 2], it can be checked that if all conditions presented in [69, Theorem 2] hold, and $h(t, \cdot)$ is a Lipschitz continuous function and satisfies that there exist $c(T), d(T) > 0$ such that:

$$ \|h(t, z)\| \le c(T) + d(T)\|z\|, \quad \forall z \in \text{dom } h(t, \cdot), \forall 0 \le t \le T, \qquad (128) $$

then the well-posedness of (127) is guaranteed. Let us now verify these conditions for (127). Since items 1 through 6 and Assumptions 8, 9, 10 hold where the identity matrix is a solution of (124), all conditions in [69, Theorem 2] hold. Since $f_i$ is Lipschitz continuous in the second variable for all $i = \overline{1, \mathcal{V}}$ (by assumption), and $\lambda^{\text{im}}(t, \cdot)$ is Lipschitz continuous [69, Lemma 6], mimicking the proof of Theorem 2, it can be checked that $h(t, \cdot)$ is Lipschitz continuous. In addition, mimicking the initial steps of Second part in the proof of [69, Theorem 2] with note that $h(t, z)$ takes the role of $f(t, z)$ in that context, then (128) holds by the Lipschitz continuity of $h(t, \cdot)$ and the validity of item 6. Therefore, the TVCS (127) has a unique absolutely continuous solution $z(\cdot)$ on $[0, T]$ for any initial condition $z(0) \in R(0)(N'(0)\Delta(0)^{-1}(\Theta^{\top} \otimes I_p)C(0))^{-1}(\mathcal{H}(0) - r'(0))$. This is equivalent to (123) has a unique absolutely continuous solution $x(\cdot)$ on $[0, T]$ for any initial condition given in (126). $\square$

A generalization of Assumption 6 is studied, in which it plays a key role in the proof of the main result.

**Assumption 11.** For all $t \ge 0$, there exists $\rho > 0$ such that $\mathcal{B}_{\rho} \subset \text{int}(\mathcal{P}(t, 0))$, where $\mathcal{B}_{\rho} = \{\gamma \in \mathbb{R}^{p|\mathcal{E}|} : \|\gamma\| \le \rho\}$, $\mathcal{P}(t, \cdot) := (N'(t)\Delta(t)^{-1})^{\top}\mathbf{\Phi}'_{tv}(t, N'(t)\Delta(t)^{-1}\cdot)$, and $\mathbf{\Phi}'_{tv}$ is given in (125).

**Lemma 8.** Suppose that Assumptions 8 and 10 hold. Then for all $t \ge 0$, the operator $\mathcal{P}(t, \cdot)$ defined in Assumption 11 is maximal monotone.

*Proof.* Take any $t \ge 0$. Since Assumption 8 holds, $\mathbf{\Phi}'_{tv}(t, \cdot)$ is maximal monotone [72, Lemma 1], and $\text{Dom } \mathbf{\Phi}'_{tv}(t, \cdot) = (\text{SOL}(0, \Gamma'(t)))^* - r'(t)$ [69, Lemma 2]. In addition, it can be verified that rge$(N'(t)\Delta(t)^{-1}) \cap \text{rint}(\text{Dom } \mathbf{\Phi}'_{tv}(t, \cdot)) \ne \emptyset$ due to Assumption 10. Hence the operator $\mathcal{P}(t, \cdot) = (N'(t)\Delta(t)^{-1})^{\top}\mathbf{\Phi}'_{tv}(t, N'(t)\Delta(t)^{-1}\cdot)$ is maximal monotone, using [54, Theorem 12.43]. $\square$
The next result follows directly from Lemma 8 and [49, Lemma 1]. It is needed for Theorem 5 below and is thus recalled.

**Lemma 9.** Let $\mathcal{P}(t, \cdot)$ be given in Assumption 11 and Assumptions 8, 10 and 11 hold. Then for all $t \ge 0$, for any $(\eta, \nu) \in \text{graph}(\mathcal{P}(t, \cdot))$, we have $\eta^\top \nu \ge \rho \|\eta\|$.

Let us consider the error state $e(t) := (\Theta_c^\top \otimes I_n)x(t)$ with $x(\cdot)$ a solution of (123), and $\Theta_c$ defined in (20). It follows from (123) that the dynamics of the error state is given by:

$$ \dot{e}(t) = (\Theta_c^\top \otimes I_n) \tilde{F}(t, x(t)) + (\Theta_c^\top \otimes I_n) B(t) (\Theta \otimes I_p) \bar{\Delta}(t)^{-1} H'(t) \lambda'(t), \quad (129a) $$

$$ w'(t) = N'(t)\Delta(t)^{-1}(\Theta^\top \otimes I_p)C(t)x(t) + \Gamma'(t)\lambda'(t) + r'(t), \quad (129b) $$

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

where $\tilde{F}(t, x(t))$ and $\Gamma'(t)$ are defined after (123). The following theorem provides some sufficient conditions that guarantee that (129) has a unique equilibrium point at the origin, which is globally exponentially stable.

**Theorem 5.** Suppose that all conditions in Lemma 7 and Theorem 4 hold, and Assumption 11 holds. Let the initial condition for the network be

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

Assume further that the following assertions hold:

1. There exist matrix functions $\bar{B}(\cdot) : \mathbb{R}_+ \to \mathbb{R}^{n \times p}$, $\bar{C}(\cdot) : \mathbb{R}_+ \to \mathbb{R}^{p \times n}$ such that $B_i(t) = \bar{B}(t)$, $C_i(t) = \bar{C}(t)$, for all $t \ge 0$, for all $i \in \mathcal{V}$.

2. There exists a continuously differentiable function $\bar{P}' : \mathbb{R}_+ \to \mathbb{R}^{n \times n}$ such that for all $t \ge 0$:

2.1. $c_1 I_n \ll \bar{P}'(t) = \bar{P}'(t)^\top \ll c_2 I_n$, for some $c_1, c_2 > 0$.

2.2. $\bar{P}'(t) \bar{B}(t) = \bar{C}(t)^\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}(t)\xi_i(t, \cdot)$, $\forall t \ge 0$, $\forall i \in \mathcal{V}$.

3.2. There exists a symmetric matrix $Q$ such that $f$ is QUAD($\bar{P}'(t), Q$), that is, $(\gamma_1 - \gamma_2)^\top \bar{P}'(t) (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$, $\forall i, j \in \mathcal{V}$.

4. $\bar{\Delta}(t)^{-1} H'(t) = (\Delta(t)^{-1})^\top N'(t)^\top$, for all $t \ge 0$.

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

$$ \dot{\bar{P}}(t) + 2Q - \gamma \lambda_2(L) \bar{C}(t)^\top \bar{C}(t) + \mu \bar{P}'(t) \ll 0, \forall t \ge 0. $$
6. The threshold value of $\rho$ given in Assumption 11 satisfies: $\rho > \frac{\xi\sqrt{|\mathcal{V}|-1}}{\sqrt{2\lambda_2(L)}}$,

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

Then the origin of the error dynamics (129) is globally exponentially stable.

*Proof.* Recall that the error state $e(t) = (\Theta_c^\top \otimes I_n)x(t)$ with $x(t)$ a solution of the TVCS (123), is governed by the dynamics given in (129). Proceeding similarly to Section 1.3, (129) is rewritten as:

$$ \dot{e}(t) \in (\Theta_c^\top \otimes I_n) (\tilde{F}(t, x(t)) - B(t)(\Theta \otimes I_p)\bar{\Delta}(t)^{-1}H'(t)\mathbf{\Phi}'_{tv}(t, N'(t)\Delta(t)^{-1}(\Theta^\top \otimes I_p)C(t)x(t))), \quad (130) $$

where $\tilde{F}(t, x(t))$ and $\mathbf{\Phi}'_{tv}$ are defined after (123), and in (125), respectively. For simplicity in 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(t, e) = \frac{1}{2} e^{\top} (I_{|\mathcal{E}_c|} \otimes \bar{P}'(t)) e, $$

where $\bar{P}'(t) = \bar{P}'(t)^\top > 0$ defined in item 2. Since item 2.1 holds, we have

$$ c_1 I_{n|\mathcal{E}_c|} \preceq (I_{|\mathcal{E}_c|} \otimes \bar{P}'(t)) \preceq c_2 I_{n|\mathcal{E}_c|}, $$

using [73, Fact 7.4.21, Fact 8.10.3]. Hence $\frac{1}{2}c_1\|e\|^2 \le V(t, e) \le \frac{1}{2}c_2\|e\|^2$, for all $t \ge 0$, for all $e \in \mathbb{R}^{n|\mathcal{E}_c|}$. The time derivative of $V$ along the trajectory of $e(t)$ is given by:

$$ \dot{V}(t, e) = \frac{1}{2} e^{\top} (I_{|\mathcal{E}_c|} \otimes \dot{\bar{P}}'(t)) e + e^{\top} (I_{|\mathcal{E}_c|} \otimes \bar{P}'(t)) \dot{e}. \quad (131) $$

It follows from (130), the definition of $\tilde{F}(t, x(t))$ (see below (123)), items 1, and 4 that:

$$ e^{\top} (I_{|\mathcal{E}_c|} \otimes \bar{P}'(t)) \dot{e} = W_1 - W_2 - W_3, $$

where the terms $W_1, W_2, W_3$ are defined as:

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

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

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

for some $\nu \in \mathbf{\Phi}'_{tv}(t, N'(t)\Delta(t)^{-1}(\Theta^\top \otimes I_p)C(t)x)$. Let us now evaluate $W_1, W_2, W_3$. First, mimicking the evaluates immediately preceding (70), the following is obtained using items 2.2 and 3:

$$ W_1 \le \frac{1}{2} e^{\top} (I_{|\mathcal{E}_c|} \otimes (2Q)) e + \xi \sqrt{|\mathcal{E}_c|} \|(I_{|\mathcal{E}_c|} \otimes \bar{C}(t))e\|. \quad (132) $$

Mimicking the evaluates preceding (71) and using items 1, and 2.2, the following holds:

$$ W_2 \ge \frac{1}{2}\gamma\lambda_2(L) \|(I_{|\mathcal{E}_c|} \otimes \bar{C}(t)^{\top}\bar{C}(t))e\|^2, \quad (133) $$
where $\gamma = \lambda_{min}(\bar{\Delta}(t)^{-1}W(t) + W(t)^{\top}\bar{\Delta}(t)^{-\top})$, and $\lambda_2(L)$ is the second smallest eigenvalue of the Laplacian matrix $L = \Theta\Theta^{\top}$. Since $\nu \in \Phi'_{tv}(t, N'(t)\Delta(t)^{-1}(\Theta^{\top} \otimes I_p)C(t)x)$ and assumptions 8, 10 (part of the conditions in Theorem 4), 11 are assumed to hold, mimicking (76) and using Lemma 9, the following is obtained:

$$W_3 \geq \rho \sqrt{|\mathcal{V}|\lambda_2(L)} \|(I_{|\mathcal{E}_c|} \otimes \bar{C}(t))e\|. \quad (134)$$

Since $|\mathcal{E}_c| = \frac{|\mathcal{V}|(|\mathcal{V}|-1)}{2}$, it follows from all the above and items 5, 6 that:

$$\dot{V}(t, e) \leq \frac{1}{2} e^{\top} \left( I_{|\mathcal{E}_c|} \otimes \dot{\bar{P}}'(t) + I_{|\mathcal{E}_c|} \otimes (2Q) - \gamma \lambda_2(L) (I_{|\mathcal{E}_c|} \otimes \bar{C}(t)^{\top} \bar{C}(t)) \right) e + \left( \xi \sqrt{|\mathcal{E}_c|} - \rho \sqrt{|\mathcal{V}|\lambda_2(L)} \right) \|(I_{|\mathcal{E}_c|} \otimes \bar{C}(t))e\| \leq -\mu V(t, e). \quad (135)$$

Thus the origin of (129) is globally exponentially stable [76, Theorem 4.10]. □

The results in Sections 3.1, 3.2, 3.3.1, 3.3.2, and 4.2 together with a global finite-time synchronization result (in the spirit of Corollary 6), could be extended to networks considered in this subsection. It also can be verified that the network in Example 4, with the time-varying capacitance $C(t)$ replacing the time-invariant capacitance $C$, achieves global exponential synchronization by applying Theorem 5.

## 7.2 Agents are Time-Varying Complementarity Systems

In this subsection, let us consider networks whose structure is represented by a connected, simple, undirected graph $G = (\mathcal{V}, \mathcal{E})$, and their agents are time-varying complementarity systems of the form:

$$\left\{
\begin{array}{l}
\dot{x}_i(t) = f_i(t, x_i(t)) + B_{i,1}\bar{\lambda}_i(t) + B_{i,2}u_i(t) \\
\bar{w}_i(t) = C_{i,1}x_i(t) + D_{i,1}\bar{\lambda}_i(t) + D_{i,2}u_i(t) \\
y_i(t) = C_{i,2}x_i(t) + D_{i,3}u_i(t) \\
0 \le \bar{\lambda}_i(t) \perp \bar{w}_i(t) \ge 0
\end{array}
\right.
, \quad \forall i \in \{1, \dots, |\mathcal{V}|\}, \quad (136)$$

where $x_i(t) \in \mathbb{R}^n$, $\bar{\lambda}_i(t) \in \mathbb{R}^l$, $u_i(t) \in \mathbb{R}^q$, $y_i(t) \in \mathbb{R}^p$ are the state, complementarity variable, 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,1}, B_{i,2}, C_{i,1}, C_{i,2}, D_{i,1}, D_{i,2}, D_{i,3}$ are matrices with compatible dimensions. Proceeding similarly to Section 1.3, the overall dynamics of agents (136) is rewritten as follows:

$$\left\{
\begin{array}{l}
\dot{x}(t) = F(t, x(t)) + B_1\bar{\lambda}(t) + B_2u(t) \\
\bar{w}(t) = C_1x(t) + D_1\bar{\lambda}(t) + D_2u(t) \\
y(t) = C_2x(t) + D_3u(t) \\
0 \le \bar{\lambda}(t) \perp \bar{w}(t) \ge 0,
\end{array}
\right.
\quad (137)$$

in which, $F(\cdot, \cdot)$ is defined in (8), $x(t) = [x_1(t)^{\top} \ x_2(t)^{\top} \ ... \ x_{|\mathcal{V}|}(t)^{\top}]^{\top} \in \mathbb{R}^{n|\mathcal{V}|}$ and a similar notation is used for $u(t) \in \mathbb{R}^{q|\mathcal{V}|}$, $y(t) \in \mathbb{R}^{p|\mathcal{V}|}$, $\bar{w}(t), \bar{\lambda}(t) \in \mathbb{R}^{l|\mathcal{V}|}$; and
$B_1 := \text{diag}(B_{1,1}, ..., B_{|V|,1})$ be a block diagonal matrix; likewise, the same notation applies to the block diagonal matrices $B_2, C_1, C_2, D_1, D_2, D_3$. Using the couplings of the form (6) and the interconnection law (7), then the systems (10) and (11) are still obtained. Proceeding similarly to Section 1.3, it follows from (137), (10) and (11) that the overall dynamics of the entire network is also a TVCS defined by:

$$ egin{cases} 
\dot{x}(t) = F(t, x(t)) + B_1 \bar{\lambda}(t) + B_2(\Theta \otimes I_q)H\lambda(t) \\
\bar{w}(t) = C_1 x(t) + D_1 \bar{\lambda}(t) + D_2(\Theta \otimes I_q)H\lambda(t) \\
w(t) = N(\Theta^{\top} \otimes I_p)C_2 x(t) + [N(\Theta^{\top} \otimes I_p)D_3(\Theta \otimes I_q)H + M]\lambda(t) + r \\
0 \le \bar{w}(t) \perp \bar{\lambda}(t) \ge 0 \\
0 \le w(t) \perp \lambda(t) \ge 0.
\end{cases} 
\quad (138) $$

Let $B := [B_1 B_2(\Theta \otimes I_q)H]$, $C := \begin{bmatrix} C_1 \\ N(\Theta^{\top} \otimes I_p)C_2 \end{bmatrix}$, $\hat{\lambda}(t) := \begin{bmatrix} \bar{\lambda}(t) \\ \lambda(t) \end{bmatrix}$, $\hat{w}(t) := \begin{bmatrix} \bar{w}(t) \\ w(t) \end{bmatrix}$, $\hat{r} := \begin{bmatrix} 0_{l|V|} \\ r \end{bmatrix}$, and $D := \begin{bmatrix} D_1 & D_2(\Theta \otimes I_q)H \\ 0 & \hat{\Gamma} \end{bmatrix}$ where $\hat{\Gamma} := N(\Theta^{\top} \otimes I_p)D_3(\Theta \otimes I_q)H + M$.

Then (138) is rewritten as:

$$ egin{cases} 
\dot{x}(t) = F(t, x(t)) + B\hat{\lambda}(t) \\
\hat{w}(t) = Cx(t) + D\hat{\lambda}(t) + \hat{r} \\
0 \le \hat{w}(t) \perp \hat{\lambda}(t) \ge 0.
\end{cases} 
\quad (139) $$

The well-posedness of (139) is guaranteed by a trivial adaptation of Theorem 2, applied to (139) instead of (13). An example of such an adaptation can be found in Corollary 3. A further extension is when exogenous terms are present in the complementarity constraints in (136). This implies that some time-varying terms are also present in the complementarity constraints in (138), so state jumps could be considered. Considering state jumps in the closed-loop is an important extension that could be studied in the future. Further information on state jumps in LCS and TVLCCS are in [95, 53, 96, 97, 98, 68, 69].

## 7.3 Other Possible Future Works and Extensions

In addition to the two above extensions of this thesis' work, we may mention the following paths that could be tackled in future works:

* Discrete-time schemes: In this thesis we have mainly tackled the design and the numerical simulation of discrete-time schemes (forward Euler, backward Euler, semi-implicit methods). Analysis of stability and convergence (when timestep vanishes) of the discretizations presented in Sections C and 6 should be tackled. A drawback of the backward Euler method is that the whole state vector has to be available to solve the generalized equation to compute the controller (see (102)). It would be quite useful to relax this.
* State-jumps: In a general setting, complementarity systems may be subjected to state discontinuities, see, e.g., [95, 96, 53, 97, 72, 69]. Such jumps are present to prevent the complementarity constraints to be violated at some time instants (either the initial time, of greater times if discontinuous exogenous terms act inside the complementarity constraints, see [68] for a summary). Clearly the influence of state jumps both in the agents (when agents are nonautonomous complementarity systems) or in the couplings (if couplings are extended to nonautonomous dynamical complementarity systems) is worth investigating.

* First-order sweeping process (FOSwP) synchronization: Indeed the framework presented above paves the way towards considering FOSwP (which are a particular type of differential inclusions introduced by J.J. Moreau in 1971/1972, and widely studied in the Applied Mathematics field since then, see [53] for references), following ideas in [60, 96, 53] to transform complementarity systems into FOSwP, using passivity-like conditions (similar to those used throughout this thesis). This could be done when time-varying exogenous terms act in the complementarity constraints (due either to synchronization of complementarity systems that involve such terms, or to some time-varying controller).

* Krasovskii-LaSalle invariance could be applied in the autonomous case (though, such an assumption may seriously restrict the class of agents and of trajectories). Could this help simplifying the feedback controller to guarantee asymptotic synchronization?

* Output-feedback synchronization: Design state observers and close the loop with the estimated state variable. It seems in view of the results about state observer design for linear complementarity systems and extensions [99, 100, 101] that this is an attainable objective.

* Exponential (resp. asymptotic) synchronization of the agents (5) with $q \neq p$ (i.e., the input and output dimensions are different): Recall that the main role of the special structure presented in Lemma 1 is to generate a linear component (with respect to the state $x(t)$) in the control input $u(\cdot)$ given in (9). In this case, $u(\cdot)$ is composed of two parts: a linear component and a nonlinear component, as defined in (36). As indicated in the proof of Theorem 3, the linear component acts on the nominal vector field $f$ of the agents to guarantee trajectories convergence, while the nonlinear component serves to compensate for perturbations $\xi_i$, thereby maintaining the network's synchronization. Therefore, for agents where the nominal dynamics are exponentially (resp. asymptotically) stable, the control input design focuses solely on compensating for perturbations acting on agents. This means that the linear component of the control input $u(\cdot)$ can be relaxed, that is, Lemma 1 is no longer necessary. In this case, the inputs $u_i(\cdot)$ and outputs $y_i(\cdot)$ in (5) may have different dimensions, and the network's stability could be analyzed directly from the TVCS (13) instead of (37).
* Piecewise-smooth systems (PWS) synchronization: A coupling (referred to as *sign coupling*) that combines the linear diffusive coupling with a discontinuous coupling defined via the signum function, is employed in [64] to synchronize networks of PWSs. Since LCP couplings proposed in this work are an extension of both linear diffusive coupling and sign coupling (see (59) in Example 1), they can be considered in the synchronization of PWSs.

# 8 Conclusions

This work proposes a new class of nonsmooth and nonlinear set-valued coupling laws constructed by linear complementarity problems (referred as LCP couplings), which is an extension of the well-known linear diffusive couplings and sign couplings introduced in [47]. These LCP couplings are proven to enable the heterogeneous networks to achieve global exponential (resp. finite-time) synchronization. An algorithm enabling the identification of the structure of such couplings is proposed. Passivity, maximal monotonicity, linear complementarity problems and variational inequalities are the main tools for analysis. The practical implementation is studied in detail with algorithms are designed for different cases. The well-posedness of the closed-loop, and the existence of equilibrium points of the error dynamics are carefully investigated. Several examples with numerical simulations are presented to illustrate the theoretical results and the proposed algorithms. An extension to networks of agents in which all data are time-varying is investigated. This work concludes with a discussion of potential directions for future research.
# Appendices

## A Dynamics of Circuits Based on Kirchhoff's Laws

In this section, we provide the dynamical systems of the electrical circuits considered in Examples 4 and 5, derived by applying Kirchhoff's Current Law (KCL) and Kirchhoff's Voltage Law (KVL).

### A.1 Dynamics of Circuits in Example 4

Let us consider the circuits given in Example 4. Recall that these circuits are described in Figure 35.

<figure><img src="image_36.png" /><figcaption>Figure 35: Circuit diagram of agent *i* in Example 4.</figcaption></figure>

Let $\mathbf{R}(t)$ be the time-varying resistance, $\mathbf{L}(t)$ be the time-varying inductance, and $\mathbf{C}$ be the capacitance; $x_{i,1}$ be the charge of the capacitor, and $x_{i,2}$ be the current through the inductor; $I_{\mathbf{C}}$ be the current through the capacitor; $U_{\mathbf{L}}$, $U_{\mathbf{R}}$, $U_{\mathbf{C}}$ be the voltages across the inductor, the resistor and the capacitor. Let $i$ denote the current in branch AB; $(u_i + \bar{\xi}_{i,2})$ be the input current, and $y_i$ be the output voltage; $I_i^{\mathbf{S}}(t)$, $\bar{\xi}_{i,1}$ be the currents of the AC current sources. In practice, the quantity $\bar{\xi}_{i,2}$ represents a perturbation originating from the input, while $\bar{\xi}_{i,1}$ corresponds to a perturbation caused by uncertainties related to the agents. Therefore, they could be time-varying and dependent on $x_{i,1}$, $x_{i,2}$. This practical framework is addressed in Example 7. Let us now use Kirchhoff's laws to derive the dynamics of each circuit. The following statements hold:

1. Applying KCL at:

1.1. Node A: $-x_{i,2} - i + u_i + \bar{\xi}_{i,2} = 0$.

1.2. Node B: $i + I_i^{\mathbf{S}} + \bar{\xi}_{i,1} - I_{\mathbf{C}} = 0$.

2. Applying KVL around the loop $A \to B \to D \to E$: $U_{\mathbf{C}} - U_{\mathbf{R}} - U_{\mathbf{L}} = 0$.

Since $I_{\mathbf{C}} = \dot{x}_{i,1}$, it follows from items 1.1 and 1.2 that: $\dot{x}_{i,1} = -x_{i,2} + I_i^{\mathbf{S}} + u_i + \bar{\xi}_{i,1} + \bar{\xi}_{i,2}$. Since $U_{\mathbf{L}} = \mathbf{L}(t)\dot{x}_{i,2} + \dot{\mathbf{L}}(t)x_{i,2}$, $U_{\mathbf{R}} = \mathbf{R}(t)x_{i,2}$, and $U_{\mathbf{C}} = \frac{1}{C}x_{i,1}$, it follows from item 2
that: $\mathbf{L}(t)\dot{x}_{i,2} = \frac{1}{C}x_{i,1} - \mathbf{R}(t)x_{i,2} - \dot{\mathbf{L}}(t)x_{i,2}$. Let $x_i := [x_{i,1} \ x_{i,2}]^\top$, and $\bar{\xi}_i := \bar{\xi}_{i,1} + \bar{\xi}_{i,2}$. From the above, the dynamics of each circuit is defined by:

$$ \left\{ \begin{aligned} \dot{x}_i(t) &= \begin{bmatrix} 0 & -1 \\ \frac{1}{C\mathbf{L}(t)} & -\frac{\mathbf{R}(t)}{\mathbf{L}(t)} - \frac{\dot{\mathbf{L}}(t)}{\mathbf{L}(t)} \end{bmatrix} x_i(t) + \begin{bmatrix} \mathbf{I}_i^S(t) \\ 0 \end{bmatrix} + \begin{bmatrix} 1 \\ 0 \end{bmatrix} (u_i(t) + \bar{\xi}_i(t, x_i)) \\ y_i(t) &= \begin{bmatrix} \frac{1}{C} & 0 \end{bmatrix} x_i(t). \end{aligned} \right. $$

## A.2 Dynamics of Circuits in Example 5

Let us consider Chua's circuits given in Example 5. Recall that these circuits are described in Figure 36.

<figure><img src="image_37.png" /><figcaption>Figure 36: Circuit diagram of agent i in Example 5.</figcaption></figure>

Let $C_1, C_2$ be capacitances, $R$ be the resistance, $L$ be the inductance, $g: \mathbb{R} \to \mathbb{R}$ be a nonlinear function that models the voltage-current characteristic of the nonlinear resistor, with voltage as the input and current as the output. Let $I_{C_1}, I_{C_2}, x_{i,3}, I_R, I_{NLR}$ (resp. $x_{i,1}, x_{i,2}, U_L, U_R, U_{NLR}$) be the currents (resp. the voltages) through (resp. across) the capacitors $C_1, C_2$, the inductor $L$, the resistor $R$, and the nonlinear resistor, respectively. Let $I_{b_1}, I_{b_2}$ denote the currents in branches EG and BA, respectively; $y_{i,1}, y_{i,2}$ be the output voltages, and $y_{i,3}$ be the output current; $(u_{i,1} + \xi_{i,1,2}), (u_{i,2} + \xi_{i,2,2})$ be the input currents, and $(u_{i,3} + \xi_{i,3,2})$ be the input voltage; $\xi_{i,1,1}, \xi_{i,2,1}$ be the currents of the AC sources, and $\xi_{i,3,1}$ be the voltage of the DC source. In practice, the quantities $\xi_{i,1,2}, \xi_{i,2,2}, \xi_{i,3,2}$ represent perturbations originating from the inputs, while $\xi_{i,1,1}, \xi_{i,2,1}, \xi_{i,3,1}$ correspond to perturbations caused by uncertainties related to the agents. Therefore, they could be time-varying and dependent on $x_{i,1}, x_{i,2}, x_{i,3}$. This practical framework is addressed in Example 8. Let us now use Kirchhoff's laws to derive the dynamics of each Chua's circuit. The following statements hold:

1. Applying KCL at:

1.1. Node A: $I_{b_2} + \xi_{i,2,1} - x_{i,3} = 0$.

1.2. Node B: $-I_{C_2} - I_{b_2} - I_R + u_{i,2} + \xi_{i,2,2} = 0$.

1.3. Node E: $-I_{C_1} - I_{b_1} + I_R + u_{i,1} + \xi_{i,1,2} = 0$.
1.4. Node G: $I_{b_1} - I_{NLR} + \xi_{i,1,1} = 0$.

2. Applying KVL around:

2.1. The loop $A \to B \to C \to D$: $\xi_{i,3,1} + x_{i,2} + u_{i,3} + \xi_{i,3,2} - U_L = 0$.

2.2. The loop $B \to E \to F \to C$: $U_R + x_{i,1} - x_{i,2} = 0$.

2.3. The loop $E \to G \to H \to F$: $U_{NLR} - x_{i,1} = 0$.

Since $I_{C_1} = C_1\dot{x}_{i,1}$, $I_{NLR} = g(U_{NLR})$, and $I_R = \frac{U_R}{R}$, it follows from items 1.3, 1.4, 2.2, and 2.3 that:

$$ 
\begin{aligned}
C_1 \dot{x}_{i,1} &= -I_{b_1} + I_R + u_{i,1} + \xi_{i,1,2} = -I_{NLR} + \xi_{i,1,1} + I_R + u_{i,1} + \xi_{i,1,2} \\
&= -g(x_{i,1}) + \xi_{i,1,1} + \frac{x_{i,2} - x_{i,1}}{R} + u_{i,1} + \xi_{i,1,2} \\
&= -\frac{1}{R}x_{i,1} + \frac{1}{R}x_{i,2} - g(x_{i,1}) + u_{i,1} + \xi_{i,1,1} + \xi_{i,1,2}.
\end{aligned} 
\quad (i1) $$

Since $I_{C_2} = C_2\dot{x}_{i,2}$ and $I_R = \frac{U_R}{R}$, it follows from items 1.2, 1.1, and 2.2 that:

$$ 
\begin{aligned}
C_2 \dot{x}_{i,2} &= -I_{b_2} - I_R + u_{i,2} + \xi_{i,2,2} = -x_{i,3} + \xi_{1,2,1} + \frac{x_{i,1} - x_{i,2}}{R} + u_{i,2} + \xi_{i,2,2} \\
&= \frac{1}{R}x_{i,1} - \frac{1}{R}x_{i,2} - x_{i,3} + u_{i,2} + \xi_{1,2,1} + \xi_{1,2,2}.
\end{aligned} 
\quad (i2) $$

Since $U_L = L\dot{x}_{i,3}$, it follows from item 2.1 that:

$$ L\dot{x}_{i,3} = x_{i,2} + u_{i,3} + \xi_{i,3,1} + \xi_{i,3,2}. \quad (i3) $$

Let $x_i := [x_{i,1} \ x_{i,2} \ x_{i,3}]^\top$, $u_i := [u_{i,1} \ u_{i,2} \ u_{i,3}]^\top$, $\xi_i := [\xi_{i,1} \ \xi_{i,2} \ \xi_{i,3}]^\top$ where for each $j = 1, 2, 3$, $\xi_{i,j} := \xi_{i,j,1} + \xi_{i,j,2}$, and $y_i := [y_{i,1} \ y_{i,2} \ y_{i,3}]^\top$. From (i1), (i2), (i3), the dynamics of each Chua's circuit is defined by:

$$ 
\left\{
\begin{aligned}
\dot{x}_i(t) &= \begin{bmatrix} -\frac{1}{RC_1} & \frac{1}{RC_1} & 0 \\ \frac{1}{RC_2} & -\frac{1}{RC_2} & -\frac{1}{C_2} \\ 0 & \frac{1}{L} & 0 \end{bmatrix} x_i(t) + \begin{bmatrix} -\frac{1}{C_1} \\ 0 \\ 0 \end{bmatrix} g(x_{i,1}(t)) \\
&\quad + \begin{bmatrix} \frac{1}{C_1} & 0 & 0 \\ 0 & \frac{1}{C_2} & 0 \\ 0 & 0 & \frac{1}{L} \end{bmatrix} (u_i(t) + \xi_i(t, x_i(t))) \\
y_i(t) &= x_i(t).
\end{aligned}
\right.
$$

## B Some useful results

**Theorem B1.** [75, Theorem 2.1] *Considering the following differential inclusion:*

$$ \dot{x}(t) \in h(t, x(t)) - \mathbf{\Psi}(x(t)) \quad (ii1) $$

where $\mathbf{\Psi}: \text{Dom } \mathbf{\Psi} \subseteq \mathbb{R}^l \Rightarrow \mathbb{R}^l$ is a maximal monotone map, $h: I \times \overline{\text{Dom } \mathbf{\Psi}} \to \mathbb{R}^l$ is measurable and for almost all $t$ in $I$ is continuous on $\text{Dom } \mathbf{\Psi}$, and such that for all $x \in \overline{\text{Dom } \mathbf{\Psi}}$:
$$ \|h(t, x)\| \le \gamma_1(t)\|x\| + \gamma_2(t), \text{ for almost all } t \text{ in } I, $$

for some $\gamma_1, \gamma_2$ are locally integrable function on $I$. Then, for any $x(0) = x_0 \in \text{Dom } \Psi$, there exists at least one AC solution to (ii1).

Moreover, if $F(t, \cdot)$ is Lipschitz continuous, the solution in Theorem B1 is unique. This is obtained using the maximal monotonicity of $\Psi(\cdot)$ and Gronwall inequality (see for instance [96, p.10]).

**Proposition B1.** [57, Example 25.13] Let $f: \mathbb{R}^n \to (-\infty, +\infty]$ be proper. Then $\partial f$ is 3\* monotone.

**Proposition B2.** [57, Example 25.15] Let $A: \mathbb{R}^n \to \mathbb{R}^n$ be monotone. Then $A$ is 3\* monotone in each of the following cases:

(i) $A$ is strongly monotone.

(ii) $A$ is uniformly monotone with a supercoercive modulus.

**Proposition B3.** [57, Proposition 22.11] Let $A: \mathbb{R}^n \to \mathbb{R}^n$ be maximally monotone and suppose that one of the following holds:

(i) $A$ is strongly monotone.

(ii) $A$ is uniformly monotone with a supercoercive modulus.

Then $A$ is surjective.

**Proposition B4.** [57, Corollary 25.27] Let $A, B: \mathbb{R}^n \to \mathbb{R}^n$ be monotone operators such that $A+B$ is maximally monotone, $A$ or $B$ is surjective, and one of the following conditions is satisfied:

(i) $A$ and $B$ are 3\* monotone.

(ii) $\text{Dom } A \subset \text{Dom } B$ and $B$ is 3\* monotone.

Then $A+B$ is surjective.

The following lemma will play a role in the analysis of equilibrium points of (21) and their stability.

**Lemma B1.** Let $\Theta \in \mathbb{R}^{|\mathcal{V}|\times|\mathcal{E}|}$ be an oriented incidence matrix of a graph $\mathcal{G} = (\mathcal{V}, \mathcal{E})$, $\Theta_c \in \mathbb{R}^{|\mathcal{V}|\times|\mathcal{E}_c|}$ be an oriented incidence matrix of the complete graph $K_{|\mathcal{V}|} = (\mathcal{V}, \mathcal{E}_c)$, where $\mathcal{E}_c = \frac{|\mathcal{V}|(|\mathcal{V}|-1)}{2}$. Then, the following assertion hold:

1. $\Theta_c \Theta_c^\top \Theta = |\mathcal{V}| \Theta$.

2. $\Theta_c^\dagger = \frac{1}{|\mathcal{V}|} \Theta_c^\top$, where $\Theta_c^\dagger \in \mathbb{R}^{|\mathcal{E}_c|\times|\mathcal{V}|}$ denotes the pseudoinverse (Moore-Penrose inverse) of $\Theta_c$.
Proof. The proof of item 1 can be found in [49, Lemma 17]. Let us prove that item 2 holds. Using item 1, all of the following four criteria (known as the Moore–Penrose conditions) hold for $\Theta_c^\dagger = \frac{1}{|\mathcal{V}|}\Theta_c^\top$:

(a) $\Theta_c \Theta_c^\dagger \Theta_c = \Theta_c$

(b) $\Theta_c^\dagger \Theta_c \Theta_c^\dagger = \Theta_c^\dagger$

(c) $(\Theta_c^\dagger \Theta_c)^\top = \Theta_c^\dagger \Theta_c$

(d) $(\Theta_c \Theta_c^\dagger)^\top = \Theta_c \Theta_c^\dagger$

Since the pseudoinverse of a matrix is unique, the proof is complete. $\square$

**Remark 15.** If $L$ is a Laplacian matrix of a connected unweighted undirected graph $G = (\mathcal{V}, \mathcal{E})$, there exists an oriented incidence matrix $\Theta \in \mathbb{R}^{|\mathcal{V}| \times |\mathcal{E}|}$ such that $L = \Theta \Theta^\top$ [51, Lemma 8.1], and the spectrum of $L$ satisfies:

$$0 = \lambda_1 < \lambda_2 \le \dots \le \lambda_{|\mathcal{V}|}, \qquad (ii2)$$

with exactly one zero eigenvalue, and all other eigenvalues positive. Let $v_1, \dots, v_{|\mathcal{V}|}$ be orthonormal eigenvectors corresponding to the eigenvalues $\lambda_1, \dots, \lambda_{|\mathcal{V}|}$, respectively. Let $R(v) := \frac{v^\top L v}{v^\top v}$, defined for all $v \ne \mathbf{0}_{|\mathcal{V}|}$, which is known as the Rayleigh quotient. For all $v \in \text{span}\{v_2, \dots, v_{|\mathcal{V}|}\}$, there exist $\alpha_2, \dots, \alpha_{|\mathcal{V}|} \in \mathbb{R}$ such that $v = \sum_{i=2}^{|\mathcal{V}|} \alpha_i v_i$. Therefore, the Rayleigh quotient satisfies:

$$R(v) = \frac{(\alpha_2 v_2 + \dots + \alpha_{|\mathcal{V}|} v_{|\mathcal{V}|})^\top (\alpha_2 \lambda_2 v_2 + \dots + \alpha_{|\mathcal{V}|} \lambda_{|\mathcal{V}|} v_{|\mathcal{V}|})}{(\alpha_2 v_2 + \dots + \alpha_{|\mathcal{V}|} v_{|\mathcal{V}|})^\top (\alpha_2 v_2 + \dots + \alpha_{|\mathcal{V}|} v_{|\mathcal{V}|})} = \frac{\sum_{i=2}^{|\mathcal{V}|} \lambda_i \alpha_i^2}{\sum_{i=2}^{|\mathcal{V}|} \alpha_i^2} \ge \lambda_2,$$

where the second equality is derived from the fact that $\{v_2, \dots, v_{|\mathcal{V}|}\}$ form an orthonormal set of vectors, and the inequality follows from (ii2). The inequality becomes an equality when $v = v_2$. Therefore, $\min_{\substack{v \ne \mathbf{0}_{|\mathcal{V}|}, \\ v \in \text{span}\{v_2, \dots, v_{|\mathcal{V}|}\}} R(v) = \lambda_2$. A general result for computing the eigenvalues $\lambda_i$, for all $i = \overline{1, |\mathcal{V}|}$, is given in [85, Theorem 4.2.6], where it is known as the Courant–Fischer minimax theorem. Since $\text{span}\{v_2, \dots, v_{|\mathcal{V}|}\} = (\text{span}\{v_1\})^\top = (\text{span}\{\mathbf{1}_{|\mathcal{V}|}\})^\top$, the following holds:

$$\min_{\substack{v \ne \mathbf{0}_{|\mathcal{V}|}, \\ v^\top \mathbf{1}_{|\mathcal{V}|} = 0}} \frac{v^\top L v}{v^\top v} = \min_{v \in \text{span}\{v_2, \dots, v_{|\mathcal{V}|}\}} \frac{v^\top L v}{v^\top v} = \lambda_2. \qquad (ii3)$$

This can also be found in [19]. Using this fact and $\Theta_c \Theta_c^\top \Theta = |\mathcal{V}| \Theta$ (see Lemma B1), the following chain of inequalities holds for all $x \in \mathbb{R}^{|\mathcal{V}|}$:

$$\begin{aligned} x^\top (L \otimes I_n) x &= \frac{1}{|\mathcal{V}|^2} ((\Theta_c \Theta_c^\top \otimes I_n) x)^\top (L \otimes I_n) ((\Theta_c \Theta_c^\top \otimes I_n) x) \\ &= \frac{1}{|\mathcal{V}|^2} y^\top (L \otimes I_n) y \ge \frac{\lambda_2}{|\mathcal{V}|^2} y^\top y = \frac{\lambda_2}{|\mathcal{V}|} x^\top (\Theta_c \Theta_c^\top \otimes I_n) x, \end{aligned} \qquad (ii4)$$

where $y = (\Theta_c \Theta_c^\top \otimes I_n)x \in \mathbb{R}^{|\mathcal{V}|}$. The inequality (ii4) plays a crucial role in the proofs of Proposition 6, Theorems 3, and 5.
# C Numerical Simulation Setup

This section provides the setup for the numerical simulations conducted in this work, including the initial conditions, parameters, code, and computational tools used.

## C.1 LCP Solver and CVXPY Tool in Python

Solving an LCP plays a crucial role in performing the numerical simulations of the examples in this work using the Semi-Backward Euler scheme. Numerous LCP solvers have been developed across various programming languages. However, in this work, all numerical simulations are conducted on Google Colaboratory - a free platform offered by Google that allows writing and running Python code in an interactive environment without local setup. Therefore, we have employed an LCP solver named **LCP-Solve.py**, which implements Lemke's Complementary Pivot algorithm. The original code was written in the GAUSS programming language by Rob Dittmar, and was later adapted to Python by Enzo Michelangeli. The solver is available as part of the OpenOpt project on GitHub, and is distributed under the MIT license, which permits free use, modification, and distribution. Detailed information is available at here.

In addition, it follows from [56, p.5] that if the matrix $Q$ in the problem LCP($q, Q$) is symmetric and positive semi-definite, then LCP($q, Q$) is completely equivalent to the following quadratic program (QP):

$$ \begin{array}{ll} \text{minimize} & f(\lambda) = \frac{1}{2}\lambda^\top Q \lambda + q^\top \lambda \\ \text{subject to} & \lambda \ge 0. \end{array} $$

This QP problem is efficiently solved using **CVXPY** - a powerful and user-friendly Python library designed for convex optimization problems. It simplifies the formulation of complex QP problems by allowing users to express objectives and constraints in intuitive mathematical syntax, and deliver fast and reliable solutions even for large-scale problems. See more details at: https://www.cvxpy.org/. Consequently, **CVXPY** Tool is useful for solving LCPs, and optimization problems such as the minimization problem in Algorithm 2 in Section 6.1.

## C.2 Code Repository and Supplementary Analysis

All numerical simulations in this work were implemented on Google Colaboratory using Python version 3.11.12, and are stored in a GitLab repository at Inria, accessible via the following link: https://gitlab.inria.fr/qpham/numerical-simulations. Therein, each example is simulated using Semi-Backward Euler scheme, as well as both LCP solvers described in Appendix C.1. The choice of LCP solver depends on the user's requirements. In certain cases, the selection of an appropriate LCP solver can lead to a substantial reduction in execution time, as illustrated in Table 5 in Section 6.4. The initial states and parameters used for simulating Examples 4, 5, 7, 8 are provided below.
**In Example 3:** The parameters $\gamma_{i,1}, \gamma_{i,2}$ are stored in a matrix named `paragamma`, which was generated using independent uniform distributions over $[-9, 9]$ for each $i = 1, 6$, with different random seeds (based on the index $i$) to ensure diversity and reproducibility. Python code is given by:

```python
paragamma = np.zeros((12, 1))
for i in range(0, 12, 2):
    rng = np.random.default rng(seed = 100 + i)
    paragamma[i, 0] = rng .uniform(-9, 9)
    paragamma[i+1, 0] = rng .uniform(-9, 9)
```

**In Examples 4, 7:** The initial condition $x(0)$ was randomly generated in Python, using a uniform distribution in the interval $[-5, 5]$, with a fixed seed for reproducibility. Python code is given by:

```python
rng_x0 = np.random.default rng(seed=1)
x0 = rng_x0 .uniform(-1, 1, (2 * 30, 1)) * 5
```

The parameters $\gamma_{i,1}, \gamma_{i,2}$ (resp. $\delta_{4,i}, \delta_{5,i}$) are stored in a matrix named `paragamma` (resp. `paradelta`), which was generated using independent uniform distributions over $[-0.4, 0.4]$ (resp. $[0, 10]$) for each $i = 1, 30$, with different random seeds (based on the index $i$) to ensure diversity and reproducibility. Python code is given by:

```python
paragamma = np .zeros((30, 2))
for i in range(30):
    rng = np.random.default rng(seed = 100 + i)
    paragamma[i, 0] = rng .uniform(-0.4, 0.4)
    paragamma[i, 1] = rng .uniform(-0.4, 0.4)
```

```python
paradelta = np .zeros((30, 2))
for i in range(30):
    rng = np.random.default rng(seed = 100 + i)
    paradelta[i, 0] = rng .uniform(0, 10)
    paradelta[i, 1] = rng .uniform(0, 10)
```

**In Examples 5, 8:** The initial condition $x(0)$ is given by the following Python code:

```python
x01 = np .array([[0.7], [0], [0]])
x02 = np .array([[0.2], [-0.4], [0.3]])
x03 = np .array([[1.5], [0.1], [-0.2]])
x04 = np .array([[0.3], [-0.5], [0.1]])
x05 = np .array([[0.1], [0.2], [-0.1]])
x06 = np .array([[1], [-0.2], [0]])
x07 = np .array([[0], [0.3], [-0.5]])
```
x08 = np.array([[-0.7], [-0.1], [0.6]])
x09 = np.array([0.2, 0.4, [-0.3]])
x10 = np.array([-1.2, -0.3, [-0.7]])
x0 = np.vstack((x01, x02, x03, x04, x05, x06, x07, x08, x09, x10))

The parameters $\gamma_{i,j,1}$, $\gamma_{i,j,2}$ are stored in a matrix named `paragamma`, which was generated using independent uniform distributions over [-0.1, 0.1] for each $i = 1, \dots, 10$, for each $j = 1, 2, 3$, with different random seeds (based on the index $i$) to ensure diversity and reproducibility. Python code is given by:

paragamma = np(zeros((3 * 10, 2)))
for i in range(0, 3 * 10, 3):
    rng = np.random.defaultenseleseed(seed = 100 + i)
    paragamma[i, 0] = rng. uniform(-0.1, 0.1)
    paragamma[i + 1, 0] = rng. uniform(-0.1, 0.1)
    paragamma[i + 2, 0] = rng. uniform(-0.1, 0.1)
    paragamma[i, 1] = rng. uniform(-0.1, 0.1)
    paragamma[i + 1, 1] = rng. uniform(-0.1, 0.1)
    paragamma[i + 2, 1] = rng. uniform(-0.1, 0.1)

## References

[1] S.N. Dorogovtsev and J.F.F. Mendes. *Evolution of Networks: From Biological Nets to the Internet and WWW*. OUP Oxford, 2013.

[2] A. Barrat, M. Barthélemy, R. Pastor-Satorras, and A. Vespignani. The architecture of complex weighted networks. *Proceedings of the National Academy of Sciences*, 101(11):3747–3752, 2004.

[3] Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem, and Alessandro Vespignani. Epidemic processes in complex networks. *Rev. Mod. Phys.*, 87:925–979, Aug 2015.

[4] Petter Holme and Juan C. Rocha. Networks of climate change: connecting causes and consequences. *Applied Network Science*, 8(1):10, 2023.

[5] Ugo Bastolla, Miguel A. Fortuna, Alberto Pascual-García, Antonio Ferrera, Bartolo Luque, and Jordi Bascompte. The architecture of mutualistic networks minimizes competition and increases biodiversity. *Nature*, 458(7241):1018–1020, 2009.

[6] S. Wasserman and K. Faust. *Social Network Analysis: Methods and Applications*. Social Network Analysis: Methods and Applications. Cambridge University Press, 1994.
[7] J. Scott. *Social Network Analysis: A Handbook*. SAGE Publications, 2000.

[8] Dana L. Haynie. Friendship networks and delinquency: The relative nature of peer delinquency. *Journal of Quantitative Criminology*, 18(2):99-134, 2002.

[9] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang. Complex networks: Structure and dynamics. *Physics Reports*, 424(4):175-308, 2006.

[10] Niloy Ganguly, Andreas Deutsch, and Animesh Mukherjee, editors. *Dynamics On and Of Complex Networks: Applications to Biology, Computer Science, and the Social Sciences*. Modeling and Simulation in Science, Engineering and Technology. Birkhäuser, Boston, MA, 2009.

[11] Luciano da Fontoura Costa, Osvaldo N. Oliveira Jr., Gonzalo Travieso, Francisco Aparecido Rodrigues, Paulino Ribeiro Villas Boas, Lucas Antiqueira, Matheus Palhares Viana, and Luis Enrique Correa Rocha. Analyzing and modeling real-world phenomena with complex networks: a survey of applications. *Advances in Physics*, 60(3):329-412, 2011.

[12] Arthur T. Winfree. Biological rhythms and the behavior of populations of coupled oscillators. *Journal of Theoretical Biology*, 16(1):15-42, 1967.

[13] Yoshiki Kuramoto. Self-entrainment of a population of coupled non-linear oscillators. In Huzihiro Araki, editor, *International Symposium on Mathematical Problems in Theoretical Physics*, pages 420-422, Berlin, Heidelberg, 1975. Springer Berlin Heidelberg.

[14] Louis M. Pecora and Thomas L. Carroll. Synchronization in chaotic systems. *Phys. Rev. Lett.*, 64:821-824, Feb 1990.

[15] Duncan J. Watts and Steven H. Strogatz. Collective dynamics of 'small-world' networks. *Nature*, 393(6684):440-442, 1998.

[16] Alex Arenas, Albert Díaz-Guilera, Jurgen Kurths, Yamir Moreno, and Chang-song Zhou. Synchronization in complex networks. *Physics Reports*, 469(3):93-153, 2008.

[17] Z. Li and G. Chen. Global synchronization and asymptotic stability of complex dynamical networks. *IEEE Transactions on Circuits and Systems II: Express Briefs*, 53(1):28-33, 2006.

[18] H. Trentelman, K. Takaba, and N. Monshizadeh. Robust synchronization of uncertain linear multi-agent systems. *Automatic Control, IEEE Transactions on*, 58:1511-1523, 06 2013.

[19] P. DeLellis, M. Di Bernardo, and G. Russo. On QUAD, Lipschitz, and contracting vector fields for consensus and synchronization of networks. *IEEE Transactions on Circuits and Systems I: Regular Papers*, 58(3):576-583, 2011.
[20] J.H. Seo, H. Shim, and J. Back. Consensus of high-order linear systems using dynamic output feedback compensator: Low gain approach. *Automatica*, 45(11):2659-2664, 2009.

[21] Mario di Bernardo, Davide Fiore, Giovanni Russo, and Francesco Scafuti. *Convergence, Consensus and Synchronization of Complex Networks via Contraction Theory*, pages 313-339. Springer Berlin Heidelberg, Berlin, Heidelberg, 2016.

[22] R. Olfati-Saber and R.M. Murray. Consensus problems in networks of agents with switching topology and time-delays. *IEEE Transactions on Automatic Control*, 49(9):1520-1533, 2004.

[23] Wei Ren and R.W. Beard. Consensus seeking in multiagent systems under dynamically changing interaction topologies. *IEEE Transactions on Automatic Control*, 50(5):655-661, 2005.

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

[25] P. DeLellis, M. di Bernardo, T.E. Gorochowski, and G. Russo. Synchronization and control of complex networks via contraction, adaptation and evolution. *IEEE Circuits and Systems Magazine*, 10(3):64-82, 2010.

[26] Pietro DeLellis, Mario diBernardo, Franco Garofalo, and Maurizio Porfiri. Evolution of complex networks via edge snapping. *IEEE Transactions on Circuits and Systems I: Regular Papers*, 57(8):2132-2143, 2010.

[27] Irinel-Constantin Morărescu, Wim Michiels, and Marc Jungers. Synchronization of coupled nonlinear oscillators with shifted gamma-distributed delays. In *2013 American Control Conference*, pages 6865-6870, 2013.

[28] Jihene Ben Rejeb, Irinel-Constantin Morărescu, and Jamal Daafouz. Control design with guaranteed cost for synchronization in networks of linear singularly perturbed systems. *Automatica*, 91:89-97, 2018.

[29] Anthony Couthures, Vineeth S. Varma, Samson Lasaulce, and Irinel-Constantin Morărescu. Global synchronization of multi-agent systems with nonlinear interactions. *IEEE Control Systems Letters*, 9:354-359, 2025.

[30] A. V. Proskurnikov. The circle criterion for synchronization in nonlinearly coupled networks. *IFAC Proceedings Volumes*, 46(23):737-742, 2013. 9th IFAC Symposium on Nonlinear Control Systems.

[31] A. V. Proskurnikov, F. Zhang, M. Cao, and J. M. A. Scherpen. A general criterion for synchronization of incrementally dissipative nonlinearly coupled agents. In *2015 European Control Conference (ECC)*, pages 581-586, 2015.
[32] Z. Liu, M. Zhang, A. Saberi, and A. A. Stoorvogel. State synchronization of multi-agent systems via static or adaptive nonlinear dynamic protocols. Automatica, 95:316-327, 2018.

[33] Erik Steur, Alexey Pavlov, and Nathan van de Wouw. Design of nonlinear coupling for efficient synchronization in networks of nonlinear systems. In 2023 62nd IEEE Conference on Decision and Control (CDC), pages 2883-2890, 2023.

[34] A. Pavlov, E. Steur, and N. van de Wouw. Nonlinear integral coupling for synchronization in networks of nonlinear systems. Automatica, 140:110202, 2022.

[35] J. Zhao, D.J. Hill, and T. Liu. Synchronization of dynamical networks with nonidentical nodes: Criteria and control. IEEE Transactions on Circuits and Systems I: Regular Papers, 58(3):584-594, 2011.

[36] E. Panteley, A. Loría, and L. Conteville. On practical synchronization of heterogeneous networks of nonlinear systems: application to chaotic systems. In 2015 American Control Conference (ACC), pages 5359-5364, 2015.

[37] Jan Maximilian Montenbruck, Georg S. Seyboth, and Frank Allgöwer. Practical and robust synchronization of systems with additive linear uncertainties. IFAC Proceedings Volumes, 46(23):743-748, 2013. 9th IFAC Symposium on Nonlinear Control Systems.

[38] J.M. Montenbruck, M. Bürger, and F. Allgöwer. Practical synchronization with diffusive couplings. Automatica, 53:235-243, 2015.

[39] P. DeLellis, M. di Bernardo, and D. Liuzza. Convergence and synchronization in heterogeneous networks of smooth and piecewise smooth systems. Automatica, 56:1-11, 2015.

[40] J. Kim, J. Yang, H. Shim, J. S. Kim, and J. H. Seo. Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents. IEEE Transactions on Automatic Control, 61(10):3096-3102, 2016.

[41] E. Panteley and A. Loría. Synchronization and dynamic consensus of heterogeneous networked systems. IEEE Transactions on Automatic Control, 62(8):3758-3773, 2017.

[42] J.G. Lee and H. Shim. A tool for analysis and synthesis of heterogeneous multi-agent systems under rank-deficient coupling. Automatica, 117:108952, 2020.

[43] H. Ahmed, R. Ushirobira, D. Efimov, and W. Perruquetti. Robust synchronization for multistable systems. IEEE Transactions on Automatic Control, 61(6):1625-1630, 2016.
[44] A. V. Proskurnikov and M. Mazo. Simple synchronization protocols for heterogeneous networks: beyond passivity. *IFAC-PapersOnLine*, 50(1):9426–9431, 2017. 20th IFAC World Congress.

[45] Wu Yang, Yan-Wu Wang, Irinel-Constantin Morărescu, Xiao-Kang Liu, and Yue-hua Huang. Fixed-time synchronization of competitive neural networks with multiple time scales. *IEEE Transactions on Neural Networks and Learning Systems*, 33(8):4133–4138, 2022.

[46] Simone Baldi and Paolo Frasca. Adaptive synchronization of unknown heterogeneous agents: An adaptive virtual model reference approach. *Journal of the Franklin Institute*, 356(2):935–955, 2019. Special Issue on Modeling, Analysis and Control of Networked Autonomous Agents.

[47] M. Coraggio, P. DeLellis, and M. di Bernardo. Distributed discontinuous coupling for convergence in heterogeneous networks. *IEEE Control Systems Letters*, 5(3):1037–1042, 2021.

[48] J.G. Lee, S. Trenn, and H. Shim. Synchronization with prescribed transient behavior: Heterogeneous multi-agent systems under funnel coupling. *Automatica*, 141:110276, 2022.

[49] F.A. Miranda-Villatoro. Robust synchronization via set-valued maximal monotone couplings. In *2023 European Control Conference (ECC)*, pages 1–6, 2023.

[50] F. Miranda-Villatoro, F. Castanos, and B. Brogliato. When proximal-point algorithms meet set-valued sliding-mode systems. An optimization point of view of discrete-time sliding modes. *IEEE Control Systems Magazine*, 45(5), 2025.

[51] F. Bullo. *Lectures on Network Systems*. Kindle Direct Publishing, 1.7 edition, 2024.

[52] I.N. Bronshtein, K.A. Semendyayev, G. Musiol, and H. Mühlig. *Handbook of Mathematics*. Springer Berlin Heidelberg, 2007.

[53] B. Brogliato and A. Tanwani. Dynamical systems coupled with monotone set-valued operators: Formalisms, applications, well-posedness, and stability. *SIAM Review*, 62(1):3–129, 2020.

[54] R.T. Rockafellar and R.J.B. Wets. *Variational Analysis*, volume 317 of *Grundlehren der mathematischen Wissenschaften*. Springer, 1998.

[55] M.K. Camlibel, L. Iannelli, and A. Tanwani. Convergence of proximal solutions for evolution inclusions with time-dependent maximal monotone operators. *Mathematical Programming A*, 194:1017–1059, 2021.
[56] R.W. Cottle, J.S. Pang, and R.E. Stone. *The Linear Complementarity Problem*. Classics in Applied Mathematics. Society for Industrial and Applied Mathematics, 1992.

[57] H.H. Bauschke and P.L. Combettes. *Convex Analysis and Monotone Operator Theory in Hilbert Spaces*. CMS Books in Mathematics. Springer, 2011.

[58] M.K. Camlibel and J.M. Schumacher. Linear passive systems and maximal monotone mappings. *Math. Program.*, 157:397–420, 2016.

[59] B. Brogliato, R. Lozano, B. Maschke, and O. Egeland. *Dissipative Systems Analysis and Control. Theory and Applications*. Communications and Control Engineering. Springer Nature Switzerland AG, third edition, 2020. Erratum at https://inria.hal.science/hal-03364924/document.

[60] B. Brogliato. Absolute stability and the Lagrange–Dirichlet theorem with monotone multivalued mappings. *Systems and Control Letters*, 51(5):343–353, 2004.

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

[62] B. Brogliato. Dissipative dynamical systems with set-valued feedback loops: Well-posed set-valued Lur'e dynamical systems. *IEEE Control Systems Magazine*, 42(3):93–114, 2022.

[63] M. Coraggio, P. De Lellis, S. Hogan, and M. Di Bernardo. Synchronization of networks of piecewise-smooth systems. *IEEE Control Systems Letters*, 2(4):653–658, 2018.

[64] M. Coraggio, P. De Lellis, and M. Di Bernardo. Convergence and synchronization in networks of piecewise-smooth systems via distributed discontinuous coupling. *Automatica*, 129:109596, 07 2021.

[65] Wenlian Lu and Tianping Chen. Quad-condition, synchronization, consensus of multiagents, and anti-synchronization of complex networks. *IEEE Transactions on Cybernetics*, 51(6):3384–3388, 2021.

[66] Pietro DeLellis, Mario diBernardo, and Francesco Garofalo. Novel decentralized adaptive strategies for the synchronization of complex networks. *Automatica*, 45(5):1312–1318, 2009.

[67] E. Hairer and G. Wanner. *Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems*. Springer Series in Computational Mathematics. Springer Berlin Heidelberg, 2010.
[68] A. Younes, F. Miranda-Villatoro, and B. Brogliato. Trajectory tracking in linear complementarity systems with and without state jumps: A passivity approach. *Nonlinear Analysis: Hybrid Systems*, 54:101520, 2024.

[69] Q.H. Pham, B. Brogliato, and F. Miranda-Villatoro. Well-posedness of passive time-varying linear cone complementarity systems. *Journal of Convex Analysis*, 32(3):681–722, 2025.

[70] B.K. Le. Well-posedness and nonsmooth Lyapunov pairs for state-dependent maximal monotone differential inclusions. *Optimization*, 69(6):1187–1217, 2020.

[71] B.K. Le. On a class of Lur'e dynamical systems with state-dependent set-valued feedback. *Set-Valued Variational Anal.*, 28(3):537–557, 2020.

[72] A. Tanwani, B. Brogliato, and C. Prieur. Well-posedness and output regulation for implicit time-varying evolution variational inequalities. *SIAM Journal on Control and Optimization*, 56(2):751–781, 2018.

[73] D.S. Bernstein. *Matrix Mathematics: Theory, Facts, and Formulas*. Princeton University Press, second edition, 2011.

[74] Q.-H. Pham and B. Brogliato. Correction to “Analysis of the Implicit Euler Time-Discretization of a Class of Descriptor-Variable Linear Cone Complementarity Systems”, J. Convex Analysis 29/2 (2022) 481–517. *Journal of Convex Analysis*, 31(3):1035–1037, July 2024.

[75] H. Attouch and A. Damlamian. On multivalued evolution equations in Hilbert spaces. *Israel Journal of Mathematics*, 12:373–390, 1972.

[76] H.K. Khalil. *Nonlinear Systems*. Pearson Education. Prentice Hall, 2002.

[77] S. Adly, A. Hantoute, and B.K. Le. Nonsmooth Lur'e dynamical systems in Hilbert spaces. *Set-Valued and Variational Analysis*, 24(1):13–35, 2016.

[78] D.C. Lay, S.R. Lay, and J.J. McDonald. *Linear Algebra and Its Applications*. Pearson Education Ltd, fifth edition, 2020.

[79] N.V. Tran, H.T.T. Le, A.V. Truong, and T.V. Phan. A fixed-time stable forward-backward dynamical system for solving generalized monotone inclusions. *Journal of Applied Mathematics and Computing*, 70(6):5857–5885, 2024.

[80] P.T. Vuong. A second order dynamical system and its discretization for strongly pseudo-monotone variational inequalities. *SIAM Journal on Control and Optimization*, 59(4):2875–2897, 2021.

[81] B. Baji and A. Cabot. An inertial proximal algorithm with dry friction: Finite convergence results. *Set-Valued Analysis*, 14:1–23, 03 2006.
[82] S. Adly, H. Attouch, and A. Cabot. *Finite Time Stabilization of Nonlinear Oscillators Subject to dry Friction*, pages 289–304. 01 2006.

[83] S. Adly and H. Attouch. First-order inertial algorithms involving dry friction damping. *Mathematical Programming*, 193, 01 2021.

[84] D.S. Bernstein. *Scalar, Vector, and Matrix Mathematics: Theory, Facts, and Formulas*. Princeton University Press, NJ, USA, 2018.

[85] R.A. Horn and C.R. Johnson. *Matrix Analysis*. Cambridge University Press, 2 edition, 2012.

[86] V. Acary and B. Brogliato. *Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics*. Lecture Notes in Applied and Computational Mechanics. Springer Berlin Heidelberg, 2008.

[87] V. Acary, O. Bonnefon, and B. Brogliato. *Nonsmooth Modeling and Simulation for Switched Circuits*, volume 69 of *Lecture Notes in Electrical Engineering*. Springer Verlag, October 2010. Erratum at https://inria.hal.science/hal-01311078v2/document.

[88] M.E.J. Newman and D.J. Watts. Renormalization group analysis of the small-world network model. *Physics Letters A*, 263(4):341–346, 1999.

[89] L.O. Chua, C.W. Wu, A. Huang, and Guo-Qun Zhong. A universal circuit for studying and generating chaos. i. routes to chaos. *IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications*, 40(10):732–744, 1993.

[90] V. Acary, B. Brogliato, and Y.V. Orlov. Chattering-free digital sliding-mode control with state observer and disturbance rejection. *IEEE Transactions on Automatic Control*, 57(5):1087–1101, 2012.

[91] B. Brogliato and A. Polyakov. Digital implementation of sliding-mode control via the implicit method: A tutorial. *International Journal of Robust and Nonlinear Control*, 31(9):3528–3586, 2021.

[92] A. Pogromsky and H. Nijmeijer. Cooperative oscillatory behavior of mutually coupled dynamical systems. *IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications*, 48(2):152–162, 2001.

[93] A. Pogromsky, G. Santoboni, and Nijmeijer. H. Partial synchronization: from symmetry towards stability. *Physica D: Nonlinear Phenomena*, 172(1):65–87, 2002.

[94] B. Brogliato. Analysis of the implicit Euler time-discretization of a class of descriptor-variable linear cone complementarity systems. *Journal of Convex Analysis*, 29(2):481–517, 2022.
[95] M.K. Çamlıbel, W.P.M.H. Heemels, and J.M. Schumacher. On linear passive complementarity systems. *European Journal of Control*, 8(3):220–237, 2002.

[96] B. Brogliato and L. Thibault. Existence and uniqueness of solutions for non-autonomous complementarity dynamical systems. *Journal of Convex Analysis*, 17(3-4):961–990, 2010.

[97] M.K. Çamlıbel, L. Iannelli, and F. Vasca. Passivity and complementarity. *Math. Program.*, 145:531–563, 2014.

[98] S. Greenhalgh, V. Acary, and B. Brogliato. On preserving dissipativity properties of linear complementarity dynamical systems with the $\theta$-method. *Numer. Math.*, 125:601–637, 2013.

[99] B. Brogliato and W. P. M. H. Heemels. Observer design for Lur'e systems with multivalued mappings: A passivity approach. *IEEE Transactions on Automatic Control*, 54(8):1996–2001, 2009.

[100] W. P. M. H. Heemels, M. K. Camlibel, J. M. Schumacher, and B. Brogliato. Observer-based control of linear complementarity systems. *International Journal of Robust and Nonlinear Control*, 21(10):1193–1218, 2011.

[101] A. Tanwani, B. Brogliato, and C. Prieur. Stability and observer design for Lur'e systems with multivalued, nonmonotone, time-varying nonlinearities and state jumps. *SIAM Journal on Control and Optimization*, 52(6):3639–3672, 2014.