Félix A Miranda-Villatoro
email: felix.miranda-villatoro@inria.fr

Robust output convergence of heterogeneous networks via nonsmooth hard-threshold couplings *

We study a set-valued maximal monotone coupling law achieving robust output convergence in heterogeneous networks of dynamical systems with uncertainties and persistent disturbances. The coupling consists of an adaptable strategy built from normal cones to convex time-dependent sets (hard-threshold maps). To guarantee the convergence of the output mismatches to a neighborhood of the origin, only connectivity of the intrinsic graph is required (knowledge of the graph algebraic connectivity is not required), whereas only the output of the associated systems is used. Numerical simulations illustrate the effectiveness of the proposed coupling scheme.

Introduction

The study of interacting dynamical systems showing synchronized trajectories has got much attention in recent years, see e.g., [START_REF] Arenas | Synchronization in complex networks[END_REF][START_REF] Delellis | Synchronization and control of complex networks via contraction, adaptation and evolution[END_REF][START_REF] Dörfler | Synchronization in complex networks of phase oscillators: a survey[END_REF][START_REF] Gambuzza | Distributed control of synchronization of a group of network nodes[END_REF][START_REF] Ghosh | The synchronized dynamics of time-varying networks[END_REF][START_REF] Russo | Symmetries, stability, and control in nonlinear systems and networks[END_REF] and references therein. Such intense focus has been motivated by the wide range of applications coming from almost all fields of science: power grids in engineering [START_REF] Motter | Spontaneous synchrony in power-grid networks[END_REF]; biological and artificial neural networks in neuroscience and computer science [START_REF] Lodi | Analyzing synchronized clusters in neuron networks[END_REF]; gene regulatory networks in biology [START_REF] Jiang | Finite-time stochastic synchronization of genetic regulatory networks[END_REF]; distributed resource allocation in operations research [START_REF] Zhu | Distributed resource allocation over directed graphs via continuous-time algorithms[END_REF]; and opinion formation of individuals in sociology [START_REF] Bizyaeva | Nonlinear opinion dynamics with tunable sensitivity[END_REF]; are just some examples from the increasing list of applications reported in the literature. In the control and systems community, the property of synchronization has become a pivotal property, studied from different perspectives and under different contexts. At first, ideal situations were considered, where disturbances and uncertainties were absent [START_REF] Pogromsky | Passivity based design of synchronizing systems[END_REF][START_REF] Pecora | Master stability functions for synchronized coupled systems[END_REF][START_REF] Proskurnikov | A general criterion for synchronization of incrementally dissipative nonlinear coupling agents[END_REF][START_REF] Wang | On partial contraction analysis for coupled nonlinear oscillators[END_REF]. In such contexts, linear, and more generally, smooth approaches, (such as the master-stability function and contraction theory, see e.g., [START_REF] Aminzare | Synchronization of diffusively-connected nonlinear systems: results based on contractions with respect to general norms[END_REF][START_REF] Trentelman | Robust synchronization of uncertain linear mutiagent systems[END_REF][START_REF] Pecora | Master stability functions for synchronized coupled systems[END_REF]), are the dominant tools used.

Nowadays, the complexity of the studied networks has increased, as more complicated dynamics are considered together with a lack of global uniformity on the models describing the individual systems, [START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF][START_REF] Kim | Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents[END_REF][START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]. Such heterogeneity in the network may be present because of different reasons. For instance, each system in the network may be described by a mathematical model that is different from that of its neighbors, with the possibility of having a state space of different dimension. In such a case, full synchronization is not possible. Nevertheless, synchronization of some of the state variables (output synchronization) may still be attainable, depending on the structure of the systems and their interactions. In other cases, even if each system belongs to the same model class, variations or uncertainties in the values of the parameters give rise to heterogeneity. Also, heterogeneity may arise due to the presence of external persistent disturbances affecting some of the agents. It is clear that when all the aforementioned sources are combined into the network, the analysis of such systems becomes a challenging task.

Some recent works dealing with the synchronization of heterogeneous networks include [START_REF] Casadei | Synchronization in networks of identical nonlinear systems via dynamic dead zones[END_REF][START_REF] Coraggio | Distributed discontinuous coupling for convergence in heterogeneous networks[END_REF][START_REF] Kim | Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents[END_REF][START_REF] Lee | Synchronization with prescribed transient behavior: heterogeneous multi-agent systems under funnel coupling[END_REF][START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF][START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF][START_REF] Zhao | Global bounded synchronization of general dynamical networks with nonidentical nodes[END_REF]. These works focus on either of the following types of interconnections: i) diffuse coupling; ii) nonlinear (high-gain) time-varying coupling; or iii) discontinuous coupling. For the case of diffusive coupling, it is shown in [START_REF] Kim | Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents[END_REF][START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF][START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF][START_REF] Zhao | Global bounded synchronization of general dynamical networks with nonidentical nodes[END_REF] that under suitable assumptions (such as connectivity of the network, and semipassivity or QUAD-ness of the vector fields), the synchronization error diminishes as the amplitude of the coupling signal increases. It is shown that, in general, an infinite interconnection signal is needed in order to force the exact invariance of the synchronization manifold, so that asymptotic synchronization cannot be attained in practice. For the case of nonlinear coupling, [START_REF] Casadei | Synchronization in networks of identical nonlinear systems via dynamic dead zones[END_REF][START_REF] Lee | Synchronization with prescribed transient behavior: heterogeneous multi-agent systems under funnel coupling[END_REF] showed, under similar assumptions as in the diffusive case, that a better performance is obtained with high-gain time-varying interconnections. In [START_REF] Lee | Synchronization with prescribed transient behavior: heterogeneous multi-agent systems under funnel coupling[END_REF] nonlinear coupling strategies, inspired from the literature of funnel control, are studied. It is shown there that such strategies achieve high-precision in the presence of heterogeneity. However, the coupling applied requires the knowledge of the full state. Similar interaction laws where proposed in [START_REF] Casadei | Synchronization in networks of identical nonlinear systems via dynamic dead zones[END_REF] regarding adaptable dead-zone maps for the case of identical systems under external noise perturbations. In such a case, the dead-zone induces an implicit funnel for the synchronization error, and the adaptation mechanism drives the funnel towards a neighborhood of the origin. Showing, once again high precision for the synchronized trajectories. Finally, for the case of discontinuous coupling, [START_REF] Coraggio | Distributed discontinuous coupling for convergence in heterogeneous networks[END_REF] established interesting properties of the network that are not shared by their smooth counterparts. Indeed, [START_REF] Coraggio | Distributed discontinuous coupling for convergence in heterogeneous networks[END_REF] showed that discontinuous couplings can achieve the asymptotic convergence of the full state with finite coupling strength. Their approach is reminiscent of sliding-mode control techniques for lumped systems, where the synchronization manifold is seen, to some extent, as a sliding surface for the synchronization error. However, it is noteworthy that the asymptotic convergence in [START_REF] Coraggio | Distributed discontinuous coupling for convergence in heterogeneous networks[END_REF] holds only at theoretical level, as the coupling is implemented in a regularized fashion to avoid the so called chattering effect. Thus, adding a boundary layer around the synchronization manifold leading to practical synchronization.

Notably, nonsmooth control techniques have shown a remarkable performance when dealing with uncertainties and disturbances, endowing the closed-loop with interesting properties such as, finite-time convergence and model reduction, see e.g., [START_REF] Levant | Higher-order sliding modes, differentiation and output-feedback control[END_REF][START_REF] Miranda-Villatoro | Robust output regulation of strongly linear passive systems with multivalued maximally monotone controls[END_REF][START_REF] Utkin | Sliding Mode in Electro-Mechanical Systems[END_REF]. In contrast, to achieve such distinctive performance, special attention must be put at implementation level, as without an appropriate implementation, those techniques may lead to the appearance of chattering, which leads to degradation in performance and the useful life of components. Regarding networks of systems, the prevalent source of nonsmoothness concerns discontinuous (switched) and impulsive couplings, these have been addressed in recent studies such as [START_REF] Chen | Synchronization with on-off coupling: role of time scales in network dynamics[END_REF][START_REF] Coraggio | Convergence and synchronization in networks of piecewise-smooth systems via distributed discontinuous coupling[END_REF][START_REF] Guan | Synchronization of complex dynamical networks with time-varying delays via impulsive distributed control[END_REF][START_REF] Jin | Quasi-synchronization of heterogeneous Lur'e networks with uncertain parameters and impulsive effect[END_REF][START_REF] Phillips | Robust distributed synchronization of networked linear systems with intermittent information[END_REF][START_REF] Stilwell | Sufficient conditions for fast switching synchronization in time-varying network topologies[END_REF][START_REF] Zhao | Synchronization of complex dynamical networks with switching topology: A switched system point of view[END_REF][START_REF] Zhou | Synchronization in slowly switching networks of coupled oscillators[END_REF]. In there, the main concern regards the synchronization of systems under interactions that change abruptly, or in cases where the communication channels put constraints on the coupling between systems. However, besides such cases, very little is known about the performance, and the implementation, of more general nonsmooth strategies. Thus, it becomes natural to consider the robust output synchronization problem under nonsmooth couplings, as a way to counteract all the effects caused by disturbances and uncertainties affecting the network.

In this paper we study the robust output convergence of heterogeneous networks using tools from the theory of differential inclusions and convex analysis. The term convergence is used in this paper, instead of the term "synchronization", as the amplitude of the coupling signals is not constrained to be small [START_REF] Pikovsky | Synchronization: A universal concept in nonlinear sciences[END_REF]Section 1.2.1]. First, some theoretical results are presented regarding perfect output convergence in presence of uncertainties and persistent disturbances. Later, practical convergence is studied by considering a discrete-time implementation of the ideal coupling. The proposed coupling scheme has a strong connection with funnel coupling [START_REF] Lee | Synchronization with prescribed transient behavior: heterogeneous multi-agent systems under funnel coupling[END_REF], adaptable dead-zone coupling [START_REF] Casadei | Synchronization in networks of identical nonlinear systems via dynamic dead zones[END_REF], constrained differential inclusions [START_REF] Papageorgiou | Differential inclusions with state constraints[END_REF], and perturbed Moreau's sweeping processes [START_REF] Moreau | Evolution problem associated with a moving convex set in a Hilbert space[END_REF]. Intuitively, it can be seen as a nonsmooth generalization of the funnel coupling in [START_REF] Lee | Synchronization with prescribed transient behavior: heterogeneous multi-agent systems under funnel coupling[END_REF] and the dead-zone coupling in [START_REF] Casadei | Synchronization in networks of identical nonlinear systems via dynamic dead zones[END_REF]. It consists in applying a correction at each time that the mismatches between the outputs of neighbors reach the boundary of a control set S(t), so that S(t) is rendered positively invariant. Such correction terms are implemented via hard-thresholds to the set S(t), so that the coupling action is inactive whenever the mismatches reside in the interior of S(t). Additionally, S(t) is designed such that it contracts asymptotically towards zero at a rate depending on the aforementioned mismatches. In this way, output convergence is easily achieved by means of high-gain couplings. The main problem that arises with the proposed strategy concerns the well-posedness (existence and uniqueness of absolutely continuous solutions) of the interconnected systems, as the operators characterizing the hardthreshold mapping are set-valued and locally unbounded on the boundary of the set S(t). That is, the associated differential inclusion is not of Filippov type [START_REF] Shtessel | Sliding Mode Control and Observation[END_REF]Definition 2.2]. Thus, special attention is paid to the well-posedness of the problem and sufficient conditions for the existence of solutions are presented in Theorem 6. Finally, it is shown that the proposed strategy is robust against parametric uncertainties, as well as, external disturbances. The main contributions of the paper are summarized as follows:

• The proof of existence of absolutely continuous solutions for set-valued heterogeneous networks with nonsmooth hard-threshold couplings, where the dimension of the state at each vertex is not necessarily the same.

• The robust output convergence in the presence of unmatched, persistent disturbances.

• The presentation of two numerical approaches, (one centralized and one fully distributed), for the implementation of the proposed nonsmooth schemes in digital computers.

• The extension of the funnel-coupling proposed in [START_REF] Lee | Synchronization with prescribed transient behavior: heterogeneous multi-agent systems under funnel coupling[END_REF] to the multivariable, nonsmooth case.

The paper is organized as follows. The next section recalls some results from convex analysis and settles the notation used all over the paper. Section 3 formulates the problem and introduces the proposed coupling strategy. The formal proof concerning the well-posedness of the coupled system is deferred to Section 6. Section 4 studies the asymptotic properties of the interconnected system in a robust framework, whereas, Section 5 presents two specific numerical approaches for the discrete-time implementation of the proposed schemes in digital computers. The paper ends with Section 7 where conclusions are presented.

Notation and preliminaries

The Euclidean inner product is represented as ⟨ • , • ⟩ and the associated norm as ∥ • ∥. The p-norm of a vector is denoted as ∥ • ∥ p , for p ∈ [1, +∞]. For the sake of simplicity we drop the subindex for the 2-norm. In the cases when the argument of the norm is a matrix, A ∈ R r×l , the induced norm is considered, that is ∥A∥ = sup ∥x∥=1 ∥Ax∥. The null and range subspaces of A are represented as, null A and rge A, respectively. The matrix A † ∈ R l×r denotes the Moore-Penrose pseudoinverse of A. The identity matrix is denoted as I l ∈ R l×l , whereas Id : R l → R l denotes the identity map. The set B l denotes the closed unit ball in R l with center at zero (in the cases where the dimension is clear from the context we will drop the subindex l). For any two nonempty sets U, V ⊆ R l , U + V = {u + v|u ∈ U, v ∈ V}, and AU = u∈U {Au}. The interior of a set U is denoted as int U. The spaces L 1 ([0, T ]; R l ) and L ∞ ([0, T ]; R l ) correspond, respectively, to the Lebesgue spaces of all absolutely integrable and essentially bounded functions, from [0, T ] into R l . For a sequence {f k } k∈N ⊂ L 1 ([0, T ]; R l ), the notation f k ⇀ f denotes convergence in the weak topology of L 1 ([0, T ]; R l ), whereas if {f k } k∈N ⊂ L ∞ ([0, T ]; R l ), then the notation f k * ⇀ f denotes convergence in the weak* topology of L ∞ ([0, T ]; R l ). The reader is addressed to [START_REF] Brézis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]Chapter 3] for further details on these notions of convergence.

Elements from convex analysis

Let U, V ⊂ R l be two nonempty, closed, convex sets. The distance from a point x ∈ R l to U is given as dist(x, U) = min y∈U ∥x -y∥ = ∥xproj(U; x)∥ , where proj(U; x) denotes the projection of x onto the set U, i.e., proj(U; x) = arg min y∈U ∥x -y∥ .

(

) 1
Let φ : R l → R ∪ {+∞} be a proper, convex, lower semicontinuous function. The convex subdifferential of φ at x is defined as the set-valued map

∂φ(x) := {w ∈ R l |⟨w, y -x⟩ ≤ φ(y) -φ(x), for all y ∈ R l } .
One of the most important properties of convex subdifferentials, that is central in the developments that follow, concerns its maximal monotonicity. A set-valued map1 M : R l ⇒ R l is monotone (in the sense of Minty-Browder) if for any two pairs (x 1 , w 1 ), (x 2 , w 2 )

∈ gph M := {(x, w) ∈ R l × R l |w ∈ M(x)}, the inequality ⟨x 1 -x 2 , w 1 -w 2 ⟩ ≥ 0 , holds.
In addition, M is maximal monotone, if it is monotone and its graph is not strictly contained in the graph of any other monotone map. Maximal monotonicity guarantees that, for any r > 0, the map (Id +rM) (

) 2
The proximal map to φ is related to the convex subdifferential of φ via the following expression, see e.g., [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Theorem 16.34], prox(φ; x) = (Id +∂φ)

-1 (x) . (3)
The right-hand side of (3) is also called the resolvent of ∂φ at x, see e.g., [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Chapter 23]. It is clear from (1) and ( 2) that the proximal map is a generalization of the projection onto a closed, convex set, as

proj(U; x) = prox(ψ U ; x) ,
where ψ U : R l → R ∪ {+∞} is the indicator function of the set U, such that ψ U (x) = 0 for all x ∈ U and ψ U (x) = +∞ otherwise. The set-valued map N(U; • ) : U ⇒ R l denotes the normal cone to U, given by

N(U; x) = ∂ψ U (x) = w ∈ R l |⟨w, y -x⟩ ≤ 0, for all y ∈ U .
Thus, it follows from (3) that proj(U;

x) = (Id +N(U; • )) -1 (x) (4)
The function σ(U; u) = sup p∈U ⟨p, u⟩ denotes the support function of the set U at u ∈ U.

Elements from graph theory

A graph G(V, E) is a mathematical structure consisting of a set of vertices V = {ν 1 , . . . , ν N } where ν i ∈ V represents the i-th vertex, together with its connections, represented via a set of edges E ⊂ V × V. So that, if {ν i , ν j } ∈ E, then the vertices ν i and ν j are connected (adjacent). Throughout the paper adjacency between vertices ν i and ν j is also denoted as ν i ∼ ν j . In addition, each edge {ν i , ν j } ∈ E is incident with the vertices ν i and ν j . When the set of vertices and edges is clear from the context we will denote the graph simply as G.

Note that, in the set notation used, {ν i , ν j } = {ν j , ν i }, that is, the graph under consideration is undirected. Let ν i , ν j ∈ V, a path of length m from ν i to ν j is a sequence of m + 1 vertices {ν k } m k=0 ⊆ V such that ν 0 = ν i , ν m = ν j , and {ν k , ν k+1 } ∈ E for k ∈ {0, . . . , m -1}. The graph G(V, E) is connected if for any pair of distinct vertices (ν i , ν j ) ∈ V × V there exists a path from ν i to ν j . For each edge ϵ k = {ν i , ν j } ∈ E a sign to each end of ϵ k is assigned. Such sign assignation will provide an orientation to the graph G(V, E). Along all the manuscript, it is assumed that an orientation has been chosen and it is fixed. Thus, the oriented incidence matrix Θ ∈ R |V|×|E| is given as, see e.g., [START_REF] Godsil | Algebraic Graph Theory[END_REF],

[Θ] i,k =      +1, if ν i is the positive end of ϵ k ; -1, if ν i is the negative end of ϵ k ; 0, otherwise.
If the graph G has |V| vertices and c connected components, then any associated incidence matrix has rank |V|c. Moreover, the graph Laplacian L = ΘΘ ⊺ and Θ ⊺ 1 N = 0.

Problem formulation and proposed coupling strategy

Let G(V, E) be an undirected and connected graph, with vertices ν i ∈ V and edges {ν i , ν j } ∈ E, characterizing the connection structure between vertices. Each vertex is associated with a nonlinear system of the form,

ν i : ẋi (t) = f i (x i (t)) + B i u i (t) + G i ζ i (t, x i (t)) ; y i (t) = C i x i (t) , (5)
where, for each i ∈ {1, . . . , |V|}, x i (t) ∈ R ni denotes the state of the i-th system at time t; u i (t), y i (t) ∈ R m denote the variables available for interconnection with other systems, as indicated by the graph G. The term ζ i (t, x i (t)) ∈ R pi is in general unknown, as it takes into account parametric uncertainties of the model, as well as, external disturbances affecting the i-th system. Finally, all the matrices are constant and of the appropriate dimensions, whereas each f i : R ni → R ni is assumed Lipschitz continuous. Throughout the paper it is assumed that each function

ζ i ( • , x i ) is measurable, whereas ζ i (t, •
) is a function of bounded variation, such that there exists a set-valued map Z i : R × R ni ⇒ R pi , that is measurable in its first argument and upper-semicontinuous in its second argument, with compact and convex images and such that ζ i (t, x i (t)) ∈ Z i (t, x i (t)) for almost all times. In addition, the following assumption on the norm of ζ i is considered.

Assumption 1. For each function ζ i there exist non-negative constants M i,1 and M i,2 such that the following inequality holds for all t ≥ 0,

∥ζ i (t, x(t))∥ ≤ M i,1 ∥x(t)∥ + M i,2 .
As stated in the introduction, we will focus on the study of a particular nonsmooth coupling law, described below, enforcing the robust agreement (against the heterogeneity of the network and external disturbances ζ i ) of all vertex outputs y i (t) towards a single common trajectory or a small neighborhood of it. These two properties are formalized in the following definition.

Definition 1. The family of systems ( 5) achieves robust asymptotic output convergence if there exist coupling laws u i (t) such that for any two indices i, j ∈ {1, . . . , |V|}

lim t→∞ ∥y i (t) -y j (t)∥ = 0 .
Likewise, it achieves practical output convergence if for any ε > 0 there exist coupling laws u i (t, ε), i ∈ {1, . . . , |V|}, such that for any two indices i, j ∈ {1, . . . , |V|} lim t→∞ dist (y i (t)y j (t); εB) = 0 .

Problem formulation Given a family of systems (5) and a graph interconnection G(V, E), our target consists on designing coupling laws u i (t) = u i (t, y j1 (t), . . . , y jp i (t)), (where each j k is such that ν j k ∼ ν i and i ∈ {1, . . . |V|}), such that robust asymptotic output convergence holds in the presence of non-vanishing disturbances.

Notice that for each vertex ν i , the external disturbances ζ i are not necessarily matched with the coupling input, that is rge G i ̸ ⊆ rge B i for some i ∈ {1, . . . |V|}. In addition, the dimension of the state space may also be different from vertex to vertex leading to a heterogeneous dynamical network. In such context, it is well-known that smooth coupling, such as diffusive coupling, cannot solve the robust convergence problem exactly, since it can provide at most practical convergence with coupling gains growing towards infinity as the mismatches between trajectories approach a neighborhood of the origin, see e.g., [START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF][START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]. On the other hand, nonsmooth control laws are able of removing, theoretically at least, all the affections caused by uncertainties and external disturbances, whilst maintaining all control signals in a bounded region, see e.g. [START_REF] Coraggio | Distributed discontinuous coupling for convergence in heterogeneous networks[END_REF][START_REF] Miranda-Villatoro | Robust output regulation of strongly linear passive systems with multivalued maximally monotone controls[END_REF][START_REF] Utkin | Sliding Mode in Electro-Mechanical Systems[END_REF].

Motivated by the robustness properties of the set-valued hard-thresholds controllers studied in [START_REF] Miranda-Villatoro | Robust output regulation of strongly linear passive systems with multivalued maximally monotone controls[END_REF], it is proposed to use the following set-valued and time-dependent coupling law

-u i (t) ∈ j∼i N(S (i,j) (t); y i (t) -y j (t)) , (6)
where j∼i indicates that the sum is taken over all indices j such that the vertex ν j is adjacent to vertex ν i and N(S (i,j) (t); • ) : S (i,j) (t) ⇒ R m denotes the normal cone to the time-dependent set S (i,j) (t). At each time instant t ≥ 0, each set S (i,j) (t) = S (j,i) (t) in ( 6) is given as

S (i,j) (t) = r (i,j) (t)C (i,j) ⊂ R m , (7)
where the sets C (i,j) are constant and satisfy the following standing assumption.

Assumption 2. Each set C (i,j) ⊂ R m is nonempty, compact, convex, symmetric (i.e., C (i,j) = -C (i,j)
), and such that 0 ∈ int C (i,j) .

Note that the symmetry condition on C (i,j) guarantees that N(S (i,j) (t); y i (t)-y j (t)) = -N(S (j,i) (t); y j (t)y i (t)), so that the coupling function in ( 6) is undirected. Figure 2 below depicts the evolution of the graph of the normal cone over time for the case when 7) is used to control the size of each set S (i,j) (t) and it obeys

C i,j = [-1, 1]. The variable r (i,j) (t) ∈ R in (
ṙ(i,j) (t) = -γ(∥y i (t) -y j (t)∥)r (i,j) (t) , (8)
with r (i,j) (0) > 0, and γ : R + → R + is a class-K function, that is, γ is strictly increasing and γ(0) = 0. In addition, γ is chosen such that there are constants Mγ , M γ ≥ 0 such that for any two outputs y i (t), y j (t),

γ(∥y i (t) -y j (t)∥) ≤ Mγ ∥y i (t) -y j (t)∥ + M γ . (9)
Note that with such initial condition, r (i,j) (t) > 0 for all t ∈ [0, +∞) and

r (i,j) (t) = e -t 0 γ(∥yi(τ )-yj (τ )∥)dτ r (i,j) (0) .
Thus, for each t ≥ 0, each set S (i,j) (t) is compact and convex, and it contracts asymptotically to the singleton {0} as time grows towards infinity.

The following example shows that the proposed coupling strategy ( 6) is not only a mathematical abstraction, but it can emerge, for instance, in the context of electrical circuits.

Example 3. Let us consider the nonsmooth diode network in Figure 1 carrying out the interconnection between vertices ν i and ν j . Let us assume that each diode satisfies an ideal complementarity condition as (see e.g., [1, Section 1.1])

ν i ν j u (i,j) ID 1 - + s(t) -+ s(t) ID 2 u (j,i) + - yi + - yj
0 ≤ I D k ⊥ V * -V D k ≥ 0 ,
where the notation 0 ≤ P ⊥ Q ≥ 0 is the short form of the following three conditions: i) P ≥ 0, ii) Q ≥ 0, iii) P Q = 0; I D k denotes the current flowing through the k-th diode D k ; V D k denotes the voltage across the terminals of D k ; and V * denotes the threshold voltage of D k for k ∈ {1, 2}. The ideal complementarity condition captures the behavior of the ideal diode in an simple way, since either: the current through the diode is zero and the voltage is lesser or equal to V * , or the voltage is equal to V * and the current is nonnegative. Setting each voltage controlled source in Figure 1 as s(t) = r (i,j) (t) -V * , it follows from Kirchhoff 's laws that (in what follows the explicit dependence on time is omitted in all variables for the sake of simplicity)

u (i,j) = -u (j,i) = I D2 -I D1 ;
(10a)

y j -y i = s + V D2 = -s -V D1 ; (10b) 0 ≤ I D1 ⊥ V * -V D1 ≥ 0 ; (10c) 0 ≤ I D2 ⊥ V * -V D2 ≥ 0 . ( 10d
)
The substitution of (10b) into the complementarity conditions (10c)-(10d) yields

0 ≤ I D1 ⊥ r (i,j) + y j -y i ≥ 0 ; (11a) 0 ≤ I D2 ⊥ r (i,j) + y i -y j ≥ 0 . ( 11b
)
It is clear from [START_REF] Byrnes | Passivity, feedback equivalence, and the global stabilization of minimum phase nonlinear systems[END_REF] that for any value of the currents I D1 and I D2 , the mismatch between the outputs of the i-th and j-th vertices satisfies -r (i,j) ≤ y jy i ≤ r (i,j) .

Hence,

1. If y jy i = -r (i,j) , then it follows from (11b) and (10a) that I D2 = 0 and u (i,j) = -I D1 ≤ 0.

2. If y jy i = r i,j , then it follows from (11a) and (10a) that I D1 = 0 and u (i,j) = I D2 ≥ 0.

3. If -r (i,j) < y jy i < r (i,j) , then I D1 = I D2 = 0 and u (i,j) = 0.

Putting all cases together we arrive at the expression

-u (i,j) ∈ N([-r (i,j) , r (i,j) ]; y i -y j ) .
Figure 2 displays the time evolution of the graph of the normal cone to the set r (i,j) (t)[-1, 1] for some function r (i,j) so that r (i,j) (t) → 0 as t ↑ ∞. Finally, joining all vertices of G(V, E) with coupling circuits as the one in Figure 1 produces a coupling law of the form (6) as

u i = j∼i u (i,j) , (12)
so that (12) indeed coincides with (6) for the scalar case where each set It is worth to emphasize that the general coupling strategy ( 6) is not exclusively implementable via analog electrical circuits. Indeed, in Section 5 we present two algorithms for the implementation of such coupling laws in a digital computer. Before that, the robustness of the coupled network ( 5)-( 6) is studied in the following section.

S (i,j) (t) = r (i,j) (t)[-1, 1]. w t N(S(t); w(t)) -r(t0) r(t0) -r(t1) r(t1) -r(t2) r(t2) -r(t3) r(t3)

Robust output convergence of heterogeneous networks

Let Θ ∈ R |V|×|E| be the associated oriented incidence matrix of the graph G(V, E). The complete network with nonsmooth coupling [START_REF] Bastien | Convergence order of implicit Euler numerical scheme for maximal monotone differential inclusions[END_REF] is written in compact form as

ẋ(t) ∈ F (x(t)) + Bu(t) + ḠZ(t, x(t)) (13a) ṙ(t) = -Γ((Θ ⊺ ⊗ I m )y(t))r(t) (13b) y(t) = Cx(t) (13c) -u(t) ∈ (Θ ⊗ I m )N (S(r(t)); (Θ ⊺ ⊗ I m )y(t)) (13d) where x(t) ∈ R n , n = |V| i=1 n i is the aggregated vector of all vertex states; F : R n → R n is such that x(t) → [f 1 (x 1 (t)) ⊺ , . . . , f |V| (x |V| (t)) ⊺ ] ⊺ ; r(t) ∈ R |E| contains all of the edge variables r(t) = [r (i1,j1) (t), . . . , r (i |E| ,j |E| ) (t)] ⊺ for {ν i k , ν j k } ∈ E; and Z : R × R p ⇒ R p , where p = |V| i p i , is the set-valued map Z(t, x(t)) = Z 1 (t, x 1 (t)) × • • • × Z |V| (t, x |V| (t)). The matrices B ∈ R n×|V|m , C ∈ R |V|m×n and Ḡ ∈ R n×p , are block diagonal given by B = Diag(B 1 , . . . , B |V| ), C = Diag(C 1 , . . . , C |V| ), and Ḡ = Diag(G 1 , . . . , G |V| ). The function Γ : R |E|m → R |E|×|E| is such that Γ (Θ ⊺ ⊗ I m ) Cx(t) =    γ(∥[(Θ ⊺ ⊗ I m ) Cx(t)] 1 ∥) . . . γ(∥[(Θ ⊺ ⊗ I m ) Cx(t)] |E| ∥)    .
Finally, the set S(r(t)) is such that

S(r(t)) = × (i,j) i<j {νi,νj }∈E r (i,j) (t)C (i,j) ⊂ R |E|m (14)
that is, S(r(t)) is the Cartesian product between all the sets indexed by the edges of the graph G.

We first formulate sufficient conditions for the existence of solutions of ( 13) and we defer the formal proof up to Section 6. By a solution of ( 13) we mean a pair of absolutely continuous functions x : R + → R n and r : R + → R |E| , satisfying the inclusion [START_REF] Castaing | Convex Analysis and Measurable Mutifunctions[END_REF] for almost all times t. To that end, we impose the following assumption.

Assumption 4. There exists a non-singular matrix R ∈ R n×n such that

R B = R -⊺ C⊺ . (15)
Assumption 4 constraints the relative degree (with respect to the input-output pair (u(t), y(t))) of the network (13) to {1, . . . , 1}, which is a necessary condition for passivity of the connected system [START_REF] Castaing | Convex Analysis and Measurable Mutifunctions[END_REF], see e.g., [START_REF] Byrnes | Passivity, feedback equivalence, and the global stabilization of minimum phase nonlinear systems[END_REF]. However, since at this point the only assumption we make on the vector fields f i concerns their Lipschitz continuity, the composed network is not necessarily passive. This lack of passivity allows us to consider the robust agreement of systems with interesting behaviors, such as self-sustained oscillations or chaotic solutions. It is easy to see that if for each i ∈ {1, . . . , |V|} there exists a symmetric, positive definite matrix,

P i ∈ R ni×ni , such that P i B i = C ⊺ i , then R = Diag(P 1/2 1 , . . . , P 1/2
|V| ) satisfies Assumption 4. Notice that in a more general setting we may have some systems for which there is no symmetric, positive definite matrix, P i , satisfying [START_REF] Coraggio | Convergence and synchronization in networks of piecewise-smooth systems via distributed discontinuous coupling[END_REF]. In addition to Assumption 4 we also consider the following standing assumption. Assumption 5. The set of admissible initial conditions is compact so that there is a constant R > 0, such that ∥r(0)∥ ≤ R, and

P i B i = C ⊺ i . Nevertheless, it might still exists R ∈ R n satisfying
(Θ ⊺ ⊗ I m )y(0) ∈ S(r(0)) .
Assumption 5 is necessary for the existence of absolutely continuous solutions. Otherwise, a state jump at t = 0 will take place so that (Θ ⊺ ⊗ I m )y(0

+ ) = lim t→0 + (Θ ⊺ ⊗ I m )y(t) ∈ S(r(0 + )).
We are now ready to state the main result concerning the existence of solutions of [START_REF] Moreau | Evolution problem associated with a moving convex set in a Hilbert space[END_REF]. The proof of Theorem 6 is postponed until Section 6. We now change focus to the consequences of Theorem 6 concerning the robust output convergence problem.

Remark 1. Note that in Theorem 6 the final time T is arbitrary but finite, and some upper-bounds in the proof depend on such T . In terms of the behavior, this translates into that fact that, even though there are not finite escape times, it is possible for the trajectories to grow unbounded as time evolves. An alternative approach consists in considering extra conditions, regarding the vertex dynamics, in order to guarantee the existence of a compact, positively invariant region, having in this way an uniform upper-bound on the state variables. For instance, if each vertex is assumed semipassive [START_REF] Pogromsky | Passivity based design of synchronizing systems[END_REF], then the trajectories of ( 13) are guaranteed to be uniformly bounded in the whole domain [0, +∞). Namely, the i-th vertex is semipassive if there exist ρ > 0 and a continuously differentiable, radially unbounded function, V i : R ni → R + such that for any admissible input u i (t),

d dt V (x i (t)) ≤ y i (t) ⊺ u i (t) -H i (x i (t))
where

H i : R ni → R is such that H i (x i (t)) ≥ M > 0 whenever ∥x i (t)∥ > ρ.
Thus, the derivative of the positive definite function

V (x) = |V| i=1 V i (x i ) is d dt V (x) ≤ |V| i=1 y ⊺ i u i - |V| i=1 H i (x i ) ≤ - |V| i=1 j∼i y ⊺ i N(S(r (i,j) (t)); y i -y j ) - |V| i=1 H i (x i ) = -y ⊺ (Θ ⊗ I m )N(S(r(t)); (Θ ⊺ ⊗ I m )y) - |V| i=1 H i (x i ) ≤ - |V| i=1 H i (x i ) , (16)
where the last inequality above follows from the maximal monotonicity of the normal cone map. Since the right-hand side of ( 16) is strictly negative whenever ∥x i ∥ > ρ for some i ∈ {1, . . . , |V|}, and the function V is radially unbounded, then all trajectories converge to the largest invariant region contained in the compact set

Ω L = {x ∈ R n |V (x) ≤ L} ,
where L is such that |V|ρB ⊂ Ω L , and boundedness of trajectories of (13) follows.

The well-posedness result in Theorem 6 conveys important consequences regarding the output convergence of general heterogeneous networks in which each agent may be subject to external persistent disturbances. Namely, well-posedness of the network [START_REF] Castaing | Convex Analysis and Measurable Mutifunctions[END_REF] implies that (Θ ⊺ ⊗I m )y(t) ∈ S(r(t)) for all times t ≥ 0. Moreover, it follows from the dynamics of r(t) in ( 13) together with the strict positivity of Γ that r(t) → 0 as t → ∞. Therefore,

(Θ ⊺ ⊗ I m )y(t) → 0 as t → ∞ .
The following corollary is thus an immediate consequence of Theorem 6.

Corollary 7. Under the assumptions of Theorem 6, the coupling law (6) achieves the robust asymptotic output convergence of the family of systems [START_REF] Arenas | Synchronization in complex networks[END_REF].

Note that, in the general heterogeneous case, the full convergence of the state is not guaranteed since, in principle, the dimension of the state of each individual system might differ from that of its neighbors. Nevertheless, in the cases when all individual systems have the same state dimension, then extra conditions can guarantee full convergence.

Definition 2. The network ( 13) is asymptotically zero-state detectable if for any two systems ν a , ν b such that the output mismatch y ay b → 0, whereas the terms ũs (t) = B s u s (t) + G s ζ(t, x s (t)), s ∈ {a, b} satisfy, ũa (t) → 0, ũb (t) → 0 implies x a (t)x b (t) → 0.

Corollary 8. Let all assumptions of Theorem 6 hold. If in addition the dimension of the state is the same for each individual system and the incremental dynamics of the network is asymptotically zero-state detectable, then the coupling (6) achieves the full-state convergence of the family of systems [START_REF] Arenas | Synchronization in complex networks[END_REF].

The results concerning asymptotic output convergence of the network are valid only when the ideal setvalued map in ( 6) is used. However, in practice, it is not possible to have a real-world implementation of the ideal coupling as depicted in Section 3, since it requires the use of ideal components. In real-life experiments there are losses due to parasitic resistance effects and unmodeled dynamics. In the next section an implementable coupling law is proposed via the implicit discretization of the continuous-time coupling [START_REF] Bastien | Convergence order of implicit Euler numerical scheme for maximal monotone differential inclusions[END_REF], so that the convergence of the output towards a unique trajectory is maintained with a precision depending on the sampling step time h (practical output convergence).

Numerical implementation of the coupling strategy

Nowadays, due to the great performance and accessibility of digital electronic devices, it may be useful to implement the coupling law discussed above using a discrete-time scheme in digital microcomputers. To make such digital implementation, it becomes necessary to pay extra attention to the discretization used, as it is well known, explicit discretization schemes are prone to the appearance of chattering, which leads to degradation of closed-loop performance or even loss of stability, see e.g., [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems[END_REF][START_REF] Bastien | Convergence order of implicit Euler numerical scheme for maximal monotone differential inclusions[END_REF]. In this section we study two approaches for computing a discrete-time coupling law based on (6).

Implementation via implicit discretization

The discretization considered in this subsection is based on the discretization scheme presented in [START_REF] Acary | Implicit Euler numerical scheme and chattering-free implementation of sliding mode systems[END_REF], where the set-valued component is discretized implicitly and the disturbances are ignored for the selection process. It is similar to the discretization used in the proof of Theorem 6 (see Section 6.2), but in the present context, we are interested only in the output dynamics, as the full state is assumed unknown. Concretely, in order to compute a selection of the generalized equation ( 6) we consider the output dynamics

ẏ(t) = C F (x(t)) + Ḡζ(t, x(t)) + C Bu(t) (17a) ṙ(t) = -Γ ((Θ ⊺ ⊗ I m )y(t)) r(t) (17b) -u(t) ∈ (Θ ⊗ I m )N (S(r(t)); (Θ ⊺ ⊗ I m )y(t)) , (17c)
where x(t) is the state of the network (13) at time t ≥ 0. The discretized output dynamics is

y k+1 = ỹk+1 + h C(F (x k ) -F ( C † y k ) + Ḡζ k ) (18a) ỹk+1 = y k + h CF ( C † y k ) + h C Bu k+1 (18b)
r k+1 = I + hΓ (Θ ⊺ ⊗ I m )y k -1 r k (18c) -u k+1 ∈ (Θ ⊗ I m )N(S(r k+1 ); (Θ ⊺ ⊗ I m )ỹ k+1 ) . (18d)
The new variable ỹk+1 in [START_REF] Delellis | On QUAD, Lipschitz, and contracting vector fields for consensus and synchronization of networks[END_REF] plays the role of a nominal output. Its role consists in making the selection strategy independent of unknown data, (i.e., independent of C(F (x k ) + Ḡζ k )). Note that, in the cases when C is invertible, as for instance C = I n , the difference F (x k ) -F ( C † y k ) = 0 for all x k and all y k . It is also noteworthy to recall that, in general, the complete state is not available, so that the term F (x k ) is unknown.

A possible way to reduce the level of uncertainty, consists in the design of state-observers as is done in [START_REF] Ahmed | Output global oscillatory synchronisation of heterogeneous systems[END_REF][START_REF] Murguia | Network synchronization by dynamic diffusive coupling[END_REF].

In this work the coupling strategy is static, so that no observer design is discussed. This, in order to keep the numerical implementation of the coupling law as simple as possible.

The following corollary establishes the well-posedness of the implicitly defined coupling (18b)-(18d), as well as, its effectiveness for solving the practical output convergence problem in a discrete-time context. Corollary 9. Let all assumptions of Theorem 6 hold. For the closed-loop heterogeneous network (18) the coupling action is given explicitly as

u k+1 = - 1 h Λ -1 Id -proj(L † S(r k+1 ); • ) (L † LΛ -1 (y k + h CF ( C † y k ))) , (19)
where L = (Θ ⊺ ⊗ I m )Λ and Λ = ( C B)

1 2
. Moreover, the "nominal ouput" ỹk+1 ∈ S(r k+1 ) for all k ≥ 0 and if the function F and ζ k are uniformly bounded then the network (18)-( 19) achieves practical output convergence.

Proof. It follows from Assumption 4 that the product C B is symmetric and positive definite. Let Λ = Λ ⊺ ≻ 0 be such that Λ 2 = C B. Then, the change of variables w k = Λ -1 y k , yields,

w k+1 = wk+1 + hΛ -1 φ k (20a) wk+1 = w k + hp k + hv k+1 (20b) r k+1 = Γ(w k )r k (20c) -v k+1 ∈ L ⊺ N(S(r k+1 ); L wk+1 ) (20d)
where

v k+1 = Λu k+1 , φ k = C(F (x k )-F ( C † Λw k )+ Ḡζ k ), p k = Λ -1 CF ( C † Λw k ), Γ(w k ) = I + hΓ Lw k -1
and L = (Θ ⊺ ⊗ I m )Λ. Note that p k is assumed to be known, as it depends only on nominal parameters and the measurable output y k = Λw k . Setting w k ⊥ = L † Lw k and w k = (I -L † L)w k as the projections onto (null L) ⊥ and null L respectively, (similarly for p k ⊥ and p k ), it follows that

w k+1 = wk+1 + hΛ -1 φ k (21a) wk+1 ⊥ = w k ⊥ + hp k ⊥ + hv k+1 (21b) wk+1 = w k + hp k (21c) r k+1 = Γ(w k ⊥ )r k (21d) -v k+1 ∈ N( Ŝ(r k+1 ); wk+1 ⊥ ) (21e)
where Ŝ(r k+1 ) = {s ∈ R |V|m |Ls ∈ S(r k+1 )} = {s ∈ R |V|m |s ∈ L † S(r k+1 )}. It follows from (21b) and (21e) that

w k ⊥ + hp k ⊥ -wk+1 ⊥ ∈ hN( Ŝ(r k+1 ); wk+1 ⊥ ) Thus, it follows from [7, Theorem 3.14] that wk+1 ⊥ = proj Ŝ(r k+1 ); w k ⊥ + hp k ⊥ , (22)
and the substitution of ( 22) back into (21b) gives us an explicit expression for v k+1 as

v k+1 = - 1 h Id -proj( Ŝ(r k+1 ); • ) (w k ⊥ + hp k ⊥ ) ,
from where [START_REF] Dörfler | Synchronization in complex networks of phase oscillators: a survey[END_REF] follows. Also note that [START_REF] Ghosh | The synchronized dynamics of time-varying networks[END_REF] implies that wk+1

⊥ ∈ Ŝ(r k+1 ), that is L wk+1 ⊥ = L wk+1 = (Θ ⊺ ⊗ I m )ỹ k+1 ∈ S(r k+1 ) ,
and (18a) leads us to

(Θ ⊺ ⊗ I m )(y k+1 -h C(F (x k ) -F ( C † y k ) + Ḡζ k )) ∈ S(r k+1 ) . (23)
If the map F and ζ k , are uniformly bounded then there is a finite M > 0 such that

(Θ ⊺ ⊗ I m )y k+1 ∈ S(r k+1 ) + hM B .
So that the error is dependent of the sampling time h. Finally, as (21d) is Schur stable, then practical output convergence follows. This concludes the proof.

Remark 2. The assumption regarding the uniform boundedness of F and ζ, is sufficient but not necessary. For instance, if F is not uniformly bounded but Lipschitz continuous instead, then it follows from ( 23) that

(Θ ⊺ ⊗ I m )y k+1 ∈ S(r k+1 ) + h(M 1 ∥x k ∥ + M 2 )B ,
for some finite M 1 , M 2 ≥ 0, so that the outputs remain close to each other whenever the trajectories rest in an invariant compact region and the sampling time h is small enough. Further assumptions on the individual vector fields, such as, semipassivity [START_REF] Pogromsky | Passivity based design of synchronizing systems[END_REF], incremental dissipativeness [START_REF] Proskurnikov | A general criterion for synchronization of incrementally dissipative nonlinear coupling agents[END_REF], or QUAD-ness of the vector field [START_REF] Delellis | On QUAD, Lipschitz, and contracting vector fields for consensus and synchronization of networks[END_REF], are common to ensure the existence of a compact invariant region of the state, allowing to relax the uniform boundedness of both F and ζ, see Remark 1 in Section 4.

It is worth to point out that the implicit discretization scheme in [START_REF] Delellis | On QUAD, Lipschitz, and contracting vector fields for consensus and synchronization of networks[END_REF] is centralized, as it uses the pseudoinverse of L, so that in order to compute the individual input u k+1 to each agent, it is necessary to have access to the output of the entire network at each time step. The following subsection presents a simple way of generating a distributed approach in order to achieve practical output convergence of the network, via a simple regularization of the ideal coupling. Example 10. Let us consider an heterogenous small-world network [START_REF] Newman | Renormalization group analysis of the small-world network model[END_REF] consisting of 16 vertices, where each vertex is of one of the following four different type of systems: A) Chua chaotic system, B) FitzHugh-Nagumo oscillator, C) Lorentz chaotic system, and D) Rössler chaotic system. Figure 3 displays the network configuration.

Each vertex is of the form

ẋi (t) = (A i + ∆A i )x i (t) + (E i + ∆E i )h i (t) + B i u i (t) g i (t) = (F i + ∆F i )x i (t) y i (t) = C i x i (t)
Let A, B, C, D, be the sets of indices for systems of type A), B), C), and D) respectively. The set of nominal parameters, A i , B i , C i , E i and F i , is shown in Table 1 for each class of vertex system i ∈ A ∪ B ∪ C ∪ D. Note that the dimension of the internal states may vary from vertex to vertex. So that in this case full convergence is not possible to achieve. The terms ∆A i , ∆E i , ∆F i denote parametric uncertainties for each model class. Each uncertainty term is taken as a sample from a standard distribution, of appropriate dimensions, with zero mean and standard deviation σ = 0.1. The function γ : R + → R + used for setting the rate of descend of r k is set as,

γ(η) = η, if η > δ ; 0, otherwise .
where δ = 0.5. In this setting, the input to each vertex is computed in a centralized fashion via [START_REF] Dörfler | Synchronization in complex networks of phase oscillators: a survey[END_REF]. Figure 4 shows the state x i (t), the output y i (t), and the input u i (t) of each vertex in the network with a sampling time of h = 10ms.

Implementation via regularization

The coupling law ( 6) is largely inspired by the electrical network in Figure 1. Thus, with the end of obtaining a distributed numerical scheme we consider a regularized version of the ideal electrical coupling discussed in Example 3. Specifically, let us consider the circuit shown in Figure 5. After simple computations, similar to those shown in Example 3, we obtain that the regularized coupling circuit obeys -u (i,j) = u (j,i) ∈ N(S (i,j) , y iy j + Ru (i,j) ) . [START_REF] Dörfler | Synchronization in complex networks of phase oscillators: a survey[END_REF].

Ai Bi Ci Ei Fi hi i ∈ A   -8.8 8.8 0 1 -1 1 0 -15 0     1 0 0   B ⊺ i   8.8 0 0   B ⊺ i tanh(2gi) + 0.7gi i ∈ B 0 -2 20 2 0 1 B ⊺ i 1 0 E ⊺ i -g 3 i + 5.5gi i ∈ C   -10 10 0 28 -1 0 0 0 -8 3     1 1 1   B ⊺ i   0 0 1 0 0 1   - -x1x3 x1x2 i ∈ D   -0.05 -0.5 -1 1 0.133 0 0 0 -1     1 1 0     1 5 0   ⊺   0 0 1     1 0 0   ⊺ max{0, 15(gi -3)}
Table 1: Nominal parameters for the three classes of systems for the network in Figure 3.

ν i ν j u (i,j) R 2 ID 1 -+ s(t) -+ s(t) ID 2 R 2 u (j,i) + - yi + - yj

Figure Regularized coupling circuit

Recalling that N(S, y) = ∂ψ S (y) is a convex cone it follows that for any R > 0 -Ru (i,j) ∈ N(S (i,j) ; y iy j + Ru (i,j) ) .

(25)

Hence, y iy j ∈ (Id +N(S (i,j) ; • ))(y iy j + Ru (i,j) ), and from (4) we retrieve the explicit expression for u (i,j) as

u (j,i) = -u (i,j) = 1 R Id -proj(S (i,j) ; • ) (y i -y j ) . ( 26
)
The coupling ( 26) is single-valued and it is shown below that it achieves practical output convergence for values of R sufficiently small. Figure 6 depicts the coupling function [START_REF] Jiang | Finite-time stochastic synchronization of genetic regulatory networks[END_REF] for the scalar case.

w t 1 Rt (Id -proj(S(t); • )) (w(t)) -r(t0) r(t0) -r(t1) r(t1) -r(t2) r(t2) -r(t3) r(t3) R t0 1 1 R t0
Figure 6: Time evolution of the graph of the coupling law in [START_REF] Jiang | Finite-time stochastic synchronization of genetic regulatory networks[END_REF]. In this case as t ↑ ∞, R t ↓ 0 and the graph of the coupling map approaches the graph of the normal to the set S(t), see Figure 2.

Note that, ( 26) is independent of the parameters of the vertices dynamics, as only the associated outputs and the set S (i,j) are needed. Thus, ( 26) can be implemented in a distributed fashion. It is also important to remark that ( 26) is a Lipschitz continuous function of the output mismatch y iy j and therefore the network ( 13)-( 26) does not have finite-escape times, that is, it is well-posed in the entire domain [0, +∞). Moreover, in the limit as R ↓ 0, the network is also well-posed, as in that case, the coupling (24) coincides with the original coupling ( 6) via [START_REF] Casadei | Synchronization in networks of identical nonlinear systems via dynamic dead zones[END_REF], see Theorem 6 and the remark after it.

Corollary 11. Under the assumptions of Theorem 6, the regularized coupling

u k+1 i = j∼i u k+1 (i,j) , (27)
where u k+1 (i,j) is given by

-u k+1 (i,j) = u k+1 (j,i) = 1 R k (i,j) Id -proj(S(r k+1 (i,j) ); • ) (y k i -y k j ) , (28)
and the regularization parameters

R k (i,j) R + satisfy R k (i,j) = R k (j,i) := min hδ 1 ∥u k (i,j) ∥ , δ 2 , (29)
for some δ i > 0 i ∈ {1, 2}, achieves the practical output convergence of (13) whenever the time step h is sufficiently small.

Proof. It follows from ( 28) that

y k i -y k j + R k (i,j) u k+1 (i,j) = proj(S(r k+1 (i,j) ); y k i -y k j ) .
Hence, ∥R k (i,j) u k+1 (i,j) ∥ = dist(y k iy k j ; S(r k+1 (i,j) )). Consequently,

y k i -y k j ∈ S(r k+1 (i,j) ) + R k (i,j) ∥u k+1 (i,j) ∥B ⊂ S(r k+1 (i,j) ) + hδ 1 ∥u k+1 (i,j) ∥ ∥u k (i,j) ∥ B .
In the limit as h ↓ 0, ∥u k+1 (i,j) ∥ ∥u k (i,j) ∥ → 1, (notice that such limit is well-defined, since in the limit we retrieve the strategy (6) whose well-posedness is guaranteed by Theorem 6), so that for h sufficiently small the output mismatch y k i -y k j lies in a neighborhood of S(r k+1 (i,j) ) and practical output convergence follows. This concludes the proof.

Example 12. Let us consider again the heterogeneous network of systems described in Example 10 above. This time each coupling is implemented in a distributed fashion using [START_REF] Kim | Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents[END_REF]. The nominal parameters of each vertex are the same as shown in Table 1 and the function γ describing the rate of descend of the variables r (i,j) also remains unchanged. The disturbances are taken in the same manner as in Example 10. In this example, the sampling time is decreased to h = 2.5ms, (as for a sampling time of h = 10ms, the coupling (28) fails to achieve the ultimate boundedness of the mismatch y k i -y k j ). Note that, in this case, no knowledge regarding the vector field F is used for the computation of the coupling signals. The regularization parameters R k (i,j) are set as in [START_REF] Kurdila | Convex Functional Analysis. Systems & Control: Foundations & Applications[END_REF] with hδ 1 = 25 and δ 2 = 0.1. Figure 7 shows the time trajectories of the state, control input, and output of each vertex system, illustrating the practical output convergence of the heterogeneous network of Figure 3. Finally, Figure 8 displays the sum of square error signals for the coupling schemes [START_REF] Dörfler | Synchronization in complex networks of phase oscillators: a survey[END_REF] and (28) under the same sampling rate. As expected, there is a trade-off between centrality and precision, since the centralized scheme [START_REF] Dörfler | Synchronization in complex networks of phase oscillators: a survey[END_REF] leads to a lower error compared to the fully distributed scheme [START_REF] Kim | Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents[END_REF]. j∼i ∥y i (t k )y j (t k )∥ 2 for coupling schemes [START_REF] Dörfler | Synchronization in complex networks of phase oscillators: a survey[END_REF](black) and (28) (gray) with a common sampling period of h = 2.5ms.

Existence of solutions of the coupled network

This section presents the proof of Theorem 6 regarding existence of absolutely continuous solutions of the differential inclusion [START_REF] Castaing | Convex Analysis and Measurable Mutifunctions[END_REF]. The proof is an adaptation and an extension of the proof in [START_REF] Haddad | Differential inclusions governed by a state dependent sweeping process[END_REF] and is presented here for completeness. The proof is divided in several subsections. First, a change of coordinates and an useful decomposition of the dynamics is performed, putting the model into a more suitable structure for its analysis. Second, a sequence of approximate solutions is constructed based on its associated discrete-time model. Third, it is shown that these sequences of approximate solutions converge in appropriated spaces. Finally, it is shown that the limit functions are indeed trajectories of the continuous-time model.

State transformation and model decomposition

Let us consider the following change of variables

z(t) = Rx(t) , (30)
where R satisfies Assumption 4. Thus, ( 13) is transformed into

   ż(t) ∈ F (z(t)) -Θ⊺ N S(r(t)); Θz(t) + (I N ⊗ RG)Z(t, (I |V| ⊗ R -1 )z(t)) , ṙ(t) = -Γ Θz(t) r(t) , (31)
where

, F = R • F • R -1 and Θ = (Θ ⊺ ⊗ I m ) CR -1 = (Θ ⊺ ⊗ I m ) B⊺ R ⊺ ∈ R |E|m×n .
It follows from Assumption 2 that 0 ∈ int S(r(t)) for all t ≥ 0. Consequently, for any ϑ ∈ null Θ, Θϑ ∈ S(r(t)) for all t ≥ 0. Hence, rge( Θ) ∩ int S(r(t)) ̸ = ∅ for all t, and it follows from [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Proposition 16.32] that Θ⊺ N(S(r(t));

Θz(t)) = ∂(ψ S(r(t)) • Θ)(z(t)) .
It follows easily from the definition of normal cone and the chain rule for convex subdifferentials that Θ⊺ N(S(r(t));

Θz(t)) = N( S(r(t)); z(t)) , (32)
where

S(r(t)) = ϑ ∈ R n | Θϑ ∈ S(r(t)) ⊂ R n , (33)
is closed and convex. Hence, (31) transforms into

ż(t) ∈ F (z(t)) -N( S(r(t)); z(t)) + Z(t, z(t)) ṙ(t) = -Γ( Θz(t))r(t) (34)
where

Z(t, z(t)) = R ḠZ(t, R -1 z(t)).
In what follows we denote the orthogonal projection of any η ∈ R n onto the subspace null Θ as η , and the complementary projection onto (null Θ) ⊥ as η ⊥ . That is, let Π be the matrix representation of the projection onto null Θ. Then, for any η ∈ R n , η = η ⊥ + η , where η ⊥ = (I -Π)η and η = Πη.

Hence, letting z(t) = z ⊥ (t) + z (t), ( 34) is rewritten as, ż⊥ (t) ∈ F⊥ (z(t)) -N( S(r(t)); z ⊥ (t)) + Z⊥ (t, z(t)) (35a)

(t) ∈ F (z(t)) + Z (t, z(t)) (35b) ṙ(t) = -Γ( Θz ⊥ (t))r(t) (35c)
where we have used the fact that Π Θ⊺ = 0, for Π a projector onto null Θ, together with (32).

Construction of approximate solutions

Clearly, that well-posedness of ( 13) is equivalent to that of [START_REF] Moreau | Evolution problem associated with a moving convex set in a Hilbert space[END_REF]. Thus, we continue by discretizing in time the dynamics (35) using a semi-implicit Euler method.

z k+1 ⊥ ∈ z k ⊥ + hq ⊥ (t k , z k ) -hN( S(r k+1 ); z k+1 ⊥ ) (36a)
z k+1 = z k + hq (t k , z k ) (36b)
r k+1 = I + hΓ( Θz k ⊥ ) -1 r k (36c)
where h := t k+1t k > 0 is the time-step, and

q(t k , z k ) = F (z k ) + ζ(t k , z k ) , so that ζ(t k , z k ) ∈ Z(t k , z k ) and q(t k , z k ) = q⊥ (t k , z k ) + q (t k , z k ).
It thus follows from the maximal monotonicity of the subdifferential map and Minty's theorem [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Theorem 21.1], that there exists a unique selection z k+1 ⊥ satisfying (36a). Indeed, the rearrangement of terms in (36a) and the use of (4) (see also [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Theorem 3.14]), gives an explicit expression for z k+1 ⊥ satisfying (36a) as

z k+1 ⊥ = proj S(r k+1 ); z k ⊥ + hq ⊥ (t k , z k ) = proj Θ † S(r k+1 ); z k ⊥ + hq ⊥ (t k , z k ) , (37)
where the second equality is a consequence of Proposition 13-ii) in A. It is worth to remark that the selection (37) is unique and it was possible to compute because of the implicit discretization of the set-value term in (36a). From [START_REF] Murguia | Network synchronization by dynamic diffusive coupling[END_REF] and in view of Assumption 5 and Proposition 13 it follows that z k ⊥ ∈ Θ † S(r k ), and z k ∈ S(r k ) for all k ≥ 0 .

It is also worth mentioning that, as Γ( Θz k ⊥ ) is a diagonal matrix with non-negative entries, the inverse in (36c) is well-defined. Moreover, since Γ( Θ(z ⊥ k )) is a non-negative diagonal matrix, it follows that all eigenvalues of I + hΓ( Θ(z k ⊥ )) lie outside the unitary ball. Recalling that for matrices we consider the induced norm, it is clear that

∥(I + hΓ( Θ(z k ⊥ )) -1 ∥ ≤ 1 for all z k ∈ R n . (38)
Hence, the recursion on r k , given by (36c), is Schur stable regardless of the value of z k . Let us consider the following family of piecewise linear functions parametrized by h.

z h (t) = z k + t -t k h (z k+1 -z k ) for t ∈ [t k , t k+1 ) , (39)
r h (t) = r k + t -t k h (r k+1 -r k ) for t ∈ [t k , t k+1 ) . (40)
It is clear from ( 39)-( 40) that z h and r h are differentiable almost everywhere in [0, T ]. In what follows it is shown that for h sufficiently small the sequences {z h } h>0 , {r h } h>0 , { żh } h>0 , { ṙh } h>0 , converge in suitable spaces.

Convergence of approximate solutions

Let 0 = t 0 < t 1 < • • • < t N = T . First, show that the derivatives żh and ṙh are bounded for t ∈ [0, T ].

To that end, let t ∈ [0, T ] \ {t 0 , t 1 , . . . , t N } and consider the following inequality

∥ żh (t)∥ = 1 h ∥z k+1 -z k ∥ ≤ 1 h ∥z k+1 -z k -hq(t k , z k )∥ + ∥q(t k , z k )∥ . (41)
Let us focus on the first term on the right-hand side of the inequality. It follows from ( 36)- [START_REF] Murguia | Network synchronization by dynamic diffusive coupling[END_REF] and the definition of projection that

∥z k+1 -z k -hq(t k , z k )∥ = ∥z k+1 ⊥ + z k+1 -z k -hq(t k , z k )∥ = ∥ proj S(r k+1 ); z k ⊥ + hq ⊥ (t k , z k ) + z k + hq (t k , z k ) -z k -hq(t k , z k )∥ = ∥ proj S(r k+1 ); z k ⊥ + hq ⊥ (t k , z k ) -z k ⊥ -hq ⊥ (t k , z k )∥ = dist(z k ⊥ + hq ⊥ (t k , z k ), S(r k+1 )) ≤ dist(z k ⊥ , S(r k+1 )) + h∥q ⊥ (t k , z k ∥ = dist(z k ⊥ , S(r k+1 )) -dist(z k ⊥ , S(r k )) + h∥q ⊥ (t k , z k )∥ ,
where we used the fact that z k ⊥ ∈ S(r k ) on the last line. Making use of Lemma 14 in A yields the upper bound

∥z k+1 -z k -hq(t k , z k )∥ ≤ L S ∥r k+1 -r k ∥ + h∥q ⊥ (t k , z k )∥ , (42)
and the substitution of ( 42) into (41) yields the estimate

∥ żh (t)∥ ≤ L S h ∥r k+1 -r k ∥ + 2∥q(t k , z k )∥ ≤ L S ∥ ṙh (t)∥ + 2(L F + Mζ )∥z k ∥ + 2(∥ F (0)∥ + M ζ ) , (43)
where L F denotes the Lipschitz constant of F and ∥ ζ(t k , z k )∥ ≤ Mζ ∥z k ∥ + M ζ follows from Assumption 1.

On the other hand, it follows from (36c) and ( 40) that

∥ ṙh (t)∥ = Γ( Θz k ⊥ )r k+1 = Γ( Θz k ⊥ )(I + hΓ( Θz k ⊥ )) -1 r k ≤ ∥Γ( Θz k ⊥ )∥ k r=1 (I + hΓ( Θz k ⊥ )) -1 ∥r(0)∥ .
Note that [START_REF] Bizyaeva | Nonlinear opinion dynamics with tunable sensitivity[END_REF] implies that there exist positive constants MΓ , M Γ such that

∥Γ( Θz k ⊥ )∥ ≤ MΓ ∥z k ∥ + M Γ .
Therefore,

∥ ṙh (t)∥ ≤ R( MΓ ∥z k ∥ + M Γ ) , ( 44
)
where R is given by Assumption 5. Thus, the substitution of ( 44) into [START_REF] Pikovsky | Synchronization: A universal concept in nonlinear sciences[END_REF], and by setting

α 1 := L S RM Γ + 2 ∥ F (0)∥ + M ζ , α 2 := 2L F + 2 Mζ + R M , yields ∥ żh (t)∥ ≤ α 1 + α 2 ∥z k ∥ = α 1 + α 2 ∥z h (t) -(t -t k ) żh (t)∥ ≤ α 1 + α 2 (t -t k )∥ żh (t)∥ + α 2 ∥z h (t)∥ ≤ α 1 + α 2 h∥ żh (t)∥ + α 2 ∥z h (t)∥ . (45)
Let 0 < δ < 1, then for all 0 < h < h * := 1-δ α2 it follows that

∥ żh (t)∥ ≤ α 1 1 -α h + α 2 1 -α 2 h ∥z h (t)∥ ≤ α 1 δ + α 2 δ ∥z h (t)∥ . (46)
The integration on both sides of (46) from t = 0 up to t = T , yields the estimate

∥z h (t) -z h (0)∥ = t 0 żh (τ )dτ ≤ t 0 ∥ żh (τ )dτ ∥ ≤ α 1 δ t + α 2 δ t 0 ∥z h (τ )∥dτ . (47)
Hence,

∥z h (t)∥ ≤ α 1 δ t + ∥z(0)∥ + α 2 δ t 0 ∥z h (τ )∥dτ .
By Gronwall's inequality, (see e.g., [53, Proposition 3.10]), it follows that

∥z h (t)∥ ≤ α 1 δ t + ∥z(0)∥ e α 2 δ t ≤ α 1 δ T + ∥z(0)∥ e α 2 δ T =: M T , (48)
and therefore,

∥ żh (t)∥ ≤ α 1 + α 2 M T δ , (49)
∥ ṙh (t)∥ ≤ R( MΓ M T + M Γ ) , (50)
∥ ζh (t, z h (t))∥ ≤ Mζ M T + M ζ . (51)
Consequently, {z h } 0<h<h , that there exists a subsequence, still denoted as { żh }, such that żh ⇀ ż weakly in L 1 ([0, T ]; R n ) as h ↓ 0. The same arguments apply for the sequences {r h } 0<h<h * and { ṙh } 0<h<h * . That is, r h → r uniformly and ṙh ⇀ ṙ weakly in L 1 ([0, T ]; R |E|m ) as h ↓ 0. Also, as ζh ∈ L ∞ ([0, T ]; R p ), it follows from the Banach-Alaoglu-Bourbaki theorem [START_REF] Brézis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]Theorem 3.16], that ζh * ⇀ ζ in the weak*-topology of L ∞ ([0, T ]; R p ). Recalling that the uniform limit of sequences of continuous functions is again continuous, and since the sequence of derivatives { żh } and { ṙh } are equi-integrable, it follows that both limit functions, z and r, are absolutely continuous in [0, T ], [START_REF] Fitzpatrick | Absolute continuity of a functions and uniform integrability of its divided differences[END_REF].

Compliance of the limit functions with the continuous time model

It rests to show that the limit functions of the previous subsection satisfy the differential inclusion [START_REF] Moreau | Evolution problem associated with a moving convex set in a Hilbert space[END_REF] for almost all t ∈ [0, T ]. To that end, let τ h (t) = t k and τ

+ h (t) = t k+1 for t ∈ [t k , t k + h). Notice that lim h→0 τ h (t) = lim h→0 τ + h (t) = t. Hence, for all t ∈ [0, T ], z k ⊥ = [z h (τ h (t))] ⊥ , r k = r h (τ h (t)), and it follows from Lemma 14 in A that dist(z ⊥ (t), Θ † S(r(t))) ≤ dist(z ⊥ (t), Θ † S(r h (τ h (t)))) + L S ∥r(t) -r h (τ h (t))∥ ≤ ∥z ⊥ (t) -[z h (τ h (t))] ⊥ ∥ + L S ∥r(t) -r h (τ h (t))∥ ,
Taking the limit as h ↓ 0 on both sides of the inequality leads us to z ⊥ (t) ∈ Θ † S(r(t)), for all t ∈ [0, T ] ,

and Proposition 13 implies that

z(t) = z ⊥ (t) + z (t) ∈ S(r(t)) for all t ∈ [0, T ] . (52)
Now it follows from ( 48)-( 49), Assumption 1, and the continuity of , that there exists a constant 0 < MT < ∞ such that

- z k+1 -z k h + q(t k , z k ) ≤ MT , for all t k ∈ [0, T ] . (53)
On the other hand, [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]Example 16.49] together with Proposition 13 yield, MT N( S(r k+1 );

z k+1 ⊥ ) ∩ B = MT ∂ dist(z k+1 ⊥ , S(r k+1 )) = MT ∂ dist(z k+1 ⊥ , Θ † S(r k+1 )) (54)
Hence, [START_REF] Zeidler | Nonlinear Functional Analysis and its Applications I: Fixed-Point Theorems[END_REF] together with [START_REF] Motter | Spontaneous synchrony in power-grid networks[END_REF] imply that

- z k+1 ⊥ -z k ⊥ h + q⊥ (t k , z k ) ∈ MT ∂ dist(z k+1 ⊥ , Θ † S(r k+1 )) , Therefore, -[ żh (τ h (t))] ⊥ + hq ⊥ (τ h (t), z h (τ h (t))) ∈ MT ∂ dist([z h (τ + h (t))] ⊥ , Θ † S(r h (τ + h (t)))) , As z h → z and r h → r uniformly as h ↓ 0, then also [z h (τ + h (t))] ⊥ → z ⊥ (t) and r h (τ + h (t)) → r(t) pointwise for all t ∈ [0, T ] and letting l ∈ L 1 ([0, T ]; R n ) be arbitrary, it follows from Lemma 15 that Ω ⟨-ż⊥ (t) + q⊥ (t, z(t)), l(t)⟩dt = lim h↓0 Ω ⟨-[ żh (τ h (t))] ⊥ + q⊥ (τ h (t), z h (τ h (t))), l(t)⟩dt ≤ lim sup h↓0 Ω ⟨-[ żh (τ h (t))] ⊥ + q⊥ (τ h (t), z h (τ h (t))), l(t)⟩dt ≤ Ω lim sup h↓0 σ MT ∂ dist([z h (τ + h (t))] ⊥ , Θ † S(r h (τ + h (t)))); l(t) dt ≤ Ω σ MT ∂ dist(z ⊥ (t); Θ † S(r(t))); l(t) dt . (55)
Since z ⊥ (t) ∈ Θ † S(r(t)) for all t ∈ [0, T ], it follows from [START_REF] Castaing | Convex Analysis and Measurable Mutifunctions[END_REF]Proposition III.35] that for almost all t ∈ [0, T ],

-[ ż(t)] ⊥ + F⊥ (z(t)) + ζ⊥ (t, z(t)) ∈ MT ∂ dist(z ⊥ (t); Θ † S(r(t))) ⊂ N( S(r(t)); z ⊥ (t)) . Analogously, for any l ∈ L ∞ ([0, T ]; R n ) Ω ⟨ ζ(t, z(t)), l(t)⟩dt = lim h↓0 Ω ⟨ ζh (τ h (t), z h (τ h (t))), l(t)⟩dt ≤ lim sup h↓0 Ω ⟨ ζh (τ h (t), z h (τ h (t))), l(t)⟩dt ≤ Ω lim sup h↓0 σ Z(τ h (t), z h (τ h (t)))); l(t) dt ≤ Ω σ Z(t,
(t) ∈ L 1 ([0, T ]; R n ) Ω ⟨ ż (t), l(t)⟩dt = lim h→0 Ω ⟨ F (z h (τ h (t))) + ζ (τ h (t), z (τ h (t))), l(t)⟩dt = Ω ⟨ F (z(t)) + ζ (t, z(t)), l(t)⟩dt , (56)
Ω ⟨ ṙ(t), l(t)⟩dt = lim h→0 Ω ⟨-Γ( Θz h (τ h (t)))r h (τ + h (t)), l(t)⟩dt = Ω ⟨-Γ( Θz ⊥ (t))r(t), l(t)⟩dt . (57)
Therefore, for almost all t ∈ [0, T ],

ż (t) = F (z(t)) + ζ (t, z(t)) ∈ F (z(t)) + Z (t, z(t)) , ṙ(t) = -Γ( Θz ⊥ (t))r(t) .
This concludes the proof.

Conclusions and further research

A set-valued coupling protocol achieving robust output convergence in networks of heterogeneous systems is presented. The general case in which each vertex system may have a different dynamical model, with possible different state space, is considered. The proposed scheme uses only the information from the outputs of vertices. It provides a generalization of the funnel coupling strategy in [START_REF] Lee | Synchronization with prescribed transient behavior: heterogeneous multi-agent systems under funnel coupling[END_REF] to the multivariable, set-valued setting, and it is also an extension of the adaptable dead-zone coupling in [START_REF] Casadei | Synchronization in networks of identical nonlinear systems via dynamic dead zones[END_REF]. Two numerical algorithms for the implementation of the coupling law in digital computers are also discussed. The former is based on the implicit discretization of the proposed set-valued coupling leading to a centralized coupling scheme in discrete time. The latter is based on a regularization from the analog circuit and it follows a black-box approach, as no parameter from the nominal dynamics of the vertices is necessary to compute the coupling signal, but it requires a faster sampling frequency in order to maintain the error bounded. Further extensions include the analysis of splitting algorithms to get a distributed way of approximating the coupling [START_REF] Dörfler | Synchronization in complex networks of phase oscillators: a survey[END_REF].

A Some technical results

In this appendix section some technical results, used mainly in the proof of Theorem 6, are presented. Proof. First notice that S(ρ) is nonempty as, by construction, for any ρ ∈ R |E| + , 0 ∈ S(ρ), so that 0 ∈ S(ρ). Moreover, since Θ as a continuous (linear) mapping from R n to R |E|m , then for any ρ ∈ int R |E| + , there exists ε > 0 such that εB ⊂ S(ρ), and int S(ρ) ̸ = ∅. Now, to show that (58) holds, notice that ϑ ∈ S(ρ) if, and only if, there is v ∈ S(ρ) such that Θϑ = v. Moreover, the last equation holds if, and only if, (see e.g., [START_REF] Bernstein | Matrix Mathematics: Theory, Facts, and Formulas[END_REF]Proposition 6.1.7]), ϑ ∈ Θ † v + null Θ , from where (58) easily follows. Now, item i) follows from (58), the fact that rge Θ † = (null Θ) ⊥ [8, Proposition 6.1.6-vii)], and the uniqueness of the splitting β = β ⊥ + β For the second item, let ϑ ∈ N( S(ρ); β ⊥ ). Hence, from (58) we have that 0 ≥ ⟨ϑ, ηβ ⊥ ⟩ for all η ∈ S(ρ) ⇔0 ≥ ⟨ϑ, η + η ⊥β ⊥ ⟩ for all η ∈ null Θ, η ⊥ ∈ Θ † S(ρ) .

Since [START_REF] Lodi | Analyzing synchronized clusters in neuron networks[END_REF] implies that ϑ ∈ rge( Θ⊺ ) = (null Θ) ⊥ , it follows that ⟨ϑ, η ⟩ = 0 for all η ∈ null Θ and then (59) implies that N( S(ρ); β ⊥ ) ⊆ N( Θ † S(ρ); β ⊥ )∩(null Θ) ⊥ . For the converse, take ϑ ∈ N( Θ † S(ρ); β ⊥ )∩(null Θ) ⊥ . Then, ⟨ϑ, η +η ⊥ -β ⊥ ⟩ = ⟨ϑ, η ⊥ -β ⊥ ⟩ for any η ∈ null Θ, and the equivalence (59) implies that N( S(ρ); β ⊥ ) ⊇ N( Θ † S(ρ); β ⊥ ) ∩ (null Θ) ⊥ . This concludes the proof.

Let U, V ⊂ R l be two nonempty, compact sets. The Pompeiu-Hausdorff distance between U and V, denoted as d H (U, V), is given as

d H (U, V) = sup x∈U ∪V | dist(x, U) -dist(x, V)| , (60)
or equivalently d H (U, V) = inf {ε ≥ 0| U ⊆ V + εB l , and V ⊆ U + εB l } ,

Thus, it follows from (60), that for any η ∈ R l ,

| dist(η, U) -dist(η, V)| ≤ d H (U, V) . (62)
The Pompeiu-Hausdorff distance measures how far two sets are from each other and it provides a notion of continuity to the motion of time-varying sets. Roughly speaking, the following lemma shows that the set S in [START_REF] Miranda-Villatoro | Robust output regulation of strongly linear passive systems with multivalued maximally monotone controls[END_REF] "moves in a Lispchitz continuous way".

Additionally, [START_REF] Chen | Synchronization with on-off coupling: role of time scales in network dynamics[END_REF] implies that Θ † S(η 2 ) ⊆ Θ † S(η 1 ) + ∥ Θ † (η 1η 2 )∥ ∞ B and Θ † S(η 1 ) ⊆ Θ † S(η

2 ) + ∥ Θ † (η 1 - η 2 )∥ ∞ B. Consequently, d H ( Θ † S(η 1 ), Θ † S(η 2 )) ≤ ∥ Θ † (η 1 -η 2 )∥ ∞ ≤ ∥ Θ † (η 1 -η 2 )∥ 2 ≤ ∥ Θ † ∥∥η 1 -η 2 ∥ . (65)
Finally, the combination of ( 64) and (65) leads us to (63). The proof is thus complete.

Finally, taking the lim sup on both sides of(69) leads us to the desired result.

Figure 1 :

 1 Figure 1: Ideal electrical circuit realizing the coupling (6) between systems ν i and ν j .

Figure 2 :

 2 Figure 2: Time evolution of the graph of the normal cone to the time-varying set S(t) = [-r(t), r(t)].

Theorem 6 .

 6 Let Assumptions 1, 2, 4 and 5 hold. Then, there exist absolutely continuous functions x : [0, T ] → R n and r : [0, T ] → R |E| that satisfy (13) almost everywhere in [0, T ].

Figure 3 :

 3 Figure 3: Heterogeneous network of systems. At each time the input to the i-th vertex is computed from the outputs of its neighbors. The parameters of each vertex system are different, even if they belong to the same class.

Figure 4 :

 4 Figure 4: (Upper left) Time trajectories of states of vertex systems ν i . (Upper right) Time trajectories of coupling inputs u i for each vertex. (Bottom) Time trajectories of output signals of each vertex system. The coupling is given by[START_REF] Dörfler | Synchronization in complex networks of phase oscillators: a survey[END_REF].

Figure 7 :

 7 Figure 7: (Upper left) Time trajectories of states of vertex systems ν i . (Upper right) Time trajectories of coupling inputs u i for each vertex. (Bottom) Time trajectories of output signals for each vertex.

Figure 8 :

 8 Figure 8: Time evolution of the sum of squares signal e SOS (t) := 1 2 |V| i=1

Proposition 13 .

 13 For any ρ ∈ R |E| + , the set S(ρ) in (33) is nonempty with nonempty interior whenever ρ ∈ int R |E| + , and S(ρ) = Θ † S(ρ) + null Θ . (58)Moreover, i) β ∈ R n ∈ S(ρ) if, and only if, β ⊥ ∈ Θ † S(ρ). ii) N( S(ρ); β ⊥ ) = N( Θ † S(ρ); β ⊥ ) ∩ (null Θ) ⊥ .

Lemma 14 .

 14 For the set S defined in[START_REF] Miranda-Villatoro | Robust output regulation of strongly linear passive systems with multivalued maximally monotone controls[END_REF] there exists a positive constant L S > 0 such that for any two η 1 , η 2 ∈ R |E| and any β ⊥ ∈ (null Θ) ⊥ , the following holds| dist(β ⊥ , S(η 1 ))dist(β ⊥ , S(η 2 ))| ≤ L S ∥η 1η 2 ∥ . (63)Proof. It follows from Proposition 13 and (62) that,| dist(β ⊥ , S(η 1 ))dist(β ⊥ , S(η 2 ))| = | dist(β ⊥ , Θ † S(η 1 ))dist(β ⊥ , Θ † S(η 2 ))|≤ d H ( Θ † S(η 1 ), Θ † S(η 2 )) .

Lemma 15 .

 15 Let {β k } k∈N ⊂ R n , {ρ k } k∈N ⊂ R |E| be two convergent sequences such that β k → β and ρ k → ρ as k → ∞. Then, for any q ∈ R n , lim sup k→∞ σ ∂ dist(β k , Θ † S(ρ k )); q ≤ σ ∂ dist(β, Θ † S(ρ)); q .(66)Proof. Let p ∈ ∂ dist(β k , ΘS(ρ k )), it thus follows from Lemma 14 that, for anys ∈ R n , ⟨p, sβ k ⟩ ≤ dist(s, Θ † k ))dist(β k , Θ † S(ρ k )) ≤ dist(s, Θ † S(ρ))dist(β k , Θ † S(ρ)) + 2L S ∥ρ k -ρ∥ .Hence, for any ε > 0 there exists k * ∈ N such that, for all k ≥ k * ,∂ dist(β k , Θ † S(ρ k )) ⊂ ∂ dist(β k , Θ † S(ρ)) + εB . (67)For any nonempty, closed, convex set Θ † S(ρ) ⊂ R n the distance function dist( • , Θ † S(ρ)) is a proper, lower semicontinuous, convex function [7, Corollary 12.12 and Proposition 12.14]. Hence, ∂ dist( • , Θ † S(ρ)) is an upper-semicontinuous set-valued map. It thus follows from (67) that, for any ε > 0 there is k* ∈ N such that for all k ≥ k * , ∂ dist(β k , Θ † S(ρ k )) ⊂ ∂ dist(β, Θ † S(ρ k )) + εB .(68)Thus, for any ε > 0, the support function satisfies,σ(∂ dist(β k , Θ † S(ρ k )); q) = sup p∈∂ dist(β k , Θ † S(ρ k ))⟨p, q⟩ ≤ sup p∈∂ dist(β, Θ † S(ρ)) ⟨p, q⟩ + ε sup b∈B ⟨b, q⟩ = σ(∂ dist(β, Θ † S(ρ)); q) + ε∥q∥ .

The double arrow is used to emphasize that M maps points to sets. See for instance[START_REF] Rockafellar | Variational Analysis[END_REF] Chapter 5].

This work is supported by the French National Research Agency in the framework of the "Investissements d'avenir" program "ANR-15-IDEX-02" and the LabEx PERSYVAL "ANR-11-LABX-0025-01".


Chinnaraji Annamalai
email: anna@iitkgp.ac.in

A Novel Mass-Energy Equation from Lorentz Factor and Energy of Motion

Keywords: classical mechanics, electromagnetic energy, light waves, speed of light

The relativistic momentum is concerned with the motion of a body whose velocity approaches the speed of light. Also, the velocity of a moving object with mass is less than the speed of light because light has no mass and the speed of light is the upper limit for the speeds of all objects in the universe. In this article, a novel mass-energy equation is derived using the Lorentz factor and the kinetic energy of classical mechanics.

Introduction

Quanta of light are the smallest discrete packets of electromagnetic energy. For instance, sunlight is a kind of electromagnetic energy. A quantum of light that has no mass is known as a photon. We know that relativistic mechanics [START_REF] Annamalai | Einstein's Mass-Energy Equivalence is Not Applicable to Photon Energy[END_REF][START_REF] Annamalai | A Mathematical Approach to the Momentum Equations of Massless Photon and Particle with Relativistic Mass[END_REF][START_REF] Annamalai | Mass-Energy Equivalence: Light Energy[END_REF][START_REF] Annamalai | E=mc^2 : Mass-Energy Equivalence[END_REF][START_REF] Annamalai | Lorentz Factor and Time Dilation on the Special Theory of Relativity[END_REF] is concerned with the motion of bodies whose speeds approach the speed of light. In the following sections of this article, let us understand the relationship between the kinetic energy of classical mechanics and the Lorentz factor [START_REF] Annamalai | Lorentz Factor and Time Dilation on the Special Theory of Relativity[END_REF].

Kinetic Energy

Kinetic energy is a fundamental concept in classical mechanics that quantifies the mechanical energy performed by an object due to its motion.

The equation of kinetic energy is given below:

𝐸 = 1 2 𝑚𝑣 2 , ( 1
)
where m is the mass and v is the velocity.

By differentiating the equation of kinetic energy with respect to time t, we obtain

𝑑𝐸 𝑑𝑡 = 1 2 (𝑚2𝑣 𝑑𝑣 𝑑𝑡 )
By simplifying the above differential equation, we receive 𝑑𝐸 = 𝑚𝑣𝑑𝑣 = 𝑣𝑑(𝑚𝑣), (2) where the mass (m) of an object in motion is constant.

Momentum-Velocity Relation

In mechanics, a variable-mass object is an object whose mass varies with time [START_REF] Annamalai | Energy-Momentum Equivalence[END_REF][START_REF] Annamalai | Momentum-Velocity Equivalence[END_REF][START_REF] Annamalai | Momentum-Velocity Equivalence[END_REF][START_REF] Annamalai | Momentum-Velocity Equivalence[END_REF][START_REF] Annamalai | Energy-Work Equivalence[END_REF]. A rocket, a moving object that burns completely at regular time intervals, or any moving object that loses its mass at regular time intervals (displacement per second), is an example of a variable-mass object. In this article, the author introduces definitions of energy and momentum with the constant-mass and variable-mass of moving objects over time.

Definition for constant-mass:

The momentum (P) of a moving object with its constant-mass (m) is directly proportional to its velocity (v).

𝑃 𝛼 𝑣 ⇒ 𝑃 = 𝑚𝑣, where 𝑚 is constant.

Here, 𝑑𝑃 𝑑𝑡 = 𝐹 = 𝑑(𝑚𝑣) 𝑑𝑡 = 𝑚 𝑑𝑣 𝑑𝑡 .
Definition for variable-mass: Suppose a moving object loses its mass (m) at every unit of time and increases its velocity (v) proportionately with its mass (m). Then, the momentum (P) of the moving object is directly proportional to its velocity (v).

𝑃 𝛼 𝑣 ⇒ 𝑃 = 𝑚𝑣 = 𝑥 𝑖 𝑣 𝑖 at the unit 𝑡 𝑖 of time for 𝑖 = 1, 2, 3, 4, ⋯.

Then, 𝑑𝑃 𝑑𝑡 = 𝐹 = 𝑣 𝑖 𝑑𝑥 𝑖 𝑑𝑡 + 𝑥 𝑖 𝑑𝑣 𝑖 𝑑𝑡 = 𝑑(𝑚𝑣) 𝑑𝑡 .
Simply, we can write the above equation as follows:

𝑑𝑃 𝑑𝑡 = 𝐹 = 𝑣 𝑑𝑚 𝑑𝑡 + 𝑚 𝑑𝑣 𝑑𝑡 = 𝑑(𝑚𝑣) 𝑑𝑡 = 𝑚 𝑑𝑣 𝑑𝑡 . ( 3
)

Derivation of the Lorentz Factor

The Lorentz factor is the factor by which time, length, and mass change for an object moving at a speed close to the speed of light. By the following experiment, the author of this article derives the equation of the Lorentz factor [START_REF] Annamalai | Energy-Work Equivalence[END_REF][START_REF] Annamalai | Computation of Mass-Energy Equation from Lorentz Factor and Kinetic Energy[END_REF].

Figure 1. depicts the relativistic speed of the metal tube (AD). Let's consider that the metal tube shown in the Figure 1 is at a stationary position and the light takes time t to travel from the point A to the point D. The length of the metal tube is measured as follows: AB = BC = h = ct, where c is the speed of light and t is the time taken to move from the point A to the point D.

Suppose the metal tube with mass m and velocity v moves from the point A to the point B horizontally and the light travels from the point A along with the metal tube. When the light reaches the point C at time 𝑡′, the metal tube moves to point B at time 𝑡′. 𝐴𝐵 = 𝑣𝑡 ′ and 𝐴𝐶 = 𝑐𝑡 ′ According to the definition of the Pythagoras theorem, we can derive the Lorentz factor as follows:

𝐴𝐶 2 = 𝐴𝐵 2 + 𝐵𝐶 2 (4)

By substituting 𝐴𝐵 = 𝑣𝑡 ′ and 𝐴𝐶 = 𝑐𝑡 ′ in the equation ( 4), we obtain

(𝑐𝑡′) 2 = (𝑣𝑡′) 2 + (𝑐𝑡) 2 𝑡′ 2 (1 - 𝑉 2 𝐶 2 ) = 𝑡 2
By simplifying the above equation, we get

𝑡 ′ = 𝑡 √ 1 - 𝑣 2 𝑐 2 𝑡′ 𝑡 = 1 √ 1 - 𝑣 2 𝑐 2 𝛾 = 1 √ 1 - 𝑣 2 𝑐 2 (5)
Hence, 𝛾 is the Lorentz factor.

Novel Mass-Energy Equation

The equation of relativistic mass is m=𝑚 0 𝛾, where the Lorentz factor 𝑖𝑠 𝛾 = 1

√ 1 - 𝑣 2 𝑐 2 . Now, m = 𝑚 0 √ 1 - 𝑣 2 𝑐 2 (6)
Squaring both sides of the equation ( 6):

𝑚 2 = 𝑚 0 2 1 - 𝑣 2 𝑐 2 ⟹ 𝑚 2 (𝑐 2 -𝑣 2 ) = 𝑚 0 2 𝑐 2
By simplifying the above equation, we obtain

𝑚 2 𝑐 2 -𝑚 2 𝑣 2 = 𝑚 0 2 𝑐 2 (7)
By differentiating the equation ( 7) with respect to time, we get

𝑐 2 2𝑚 𝑑𝑚 𝑑𝑡 -[2𝑚𝑣 2 𝑑𝑚 𝑑𝑡 + 2𝑚 2 𝑣 𝑑𝑣 𝑑𝑡 ] = 0
By simplifying the above equation, we receive

𝑐 2 𝑑𝑚 = 𝑣𝑑(𝑚𝑣) ⇒ 𝑐𝑑(𝑐𝑚) = 𝑣𝑑(𝑚𝑣) (8)
Substituting the equation ( 8) in the differential equation ( 2) of kinetic energy:

𝑑𝐸 = 𝑐𝑑(𝑐𝑚) = 𝑣𝑑(𝑚𝑣), ( 9
)
where 𝑐 is a numeric constant and also one of the values of 𝑣.

Therefore, 𝑣 = 𝑐. [START_REF] Annamalai | Energy-Work Equivalence[END_REF] By substituting the equation [START_REF] Annamalai | Energy-Work Equivalence[END_REF] in the equation ( 1), we conclude that

𝐸 = 1 2 𝑚𝑐 2
, where 𝐸 is the energy of motion.

Let us understand the definition of momentum [START_REF] Annamalai | Energy-Momentum Equivalence[END_REF][START_REF] Annamalai | Momentum-Velocity Equivalence[END_REF][START_REF] Annamalai | Momentum-Velocity Equivalence[END_REF][START_REF] Annamalai | Momentum-Velocity Equivalence[END_REF][START_REF] Annamalai | Energy-Work Equivalence[END_REF]: Momentum (P) is the product of mass (m) of a moving object and its velocity (v). i.e., P = mv.

Here, v = c denotes the speed of light that is the upper limit of velocity for all moving objects in the universe and light has no mass. Therefore, momentum (P) = mc is wrong. i.e., the velocity (v) of an object with its mass (m) must be less than the speed of light (c).

𝑇ℎ𝑒𝑛, 𝐸 = 𝑃 𝑐 2 = 𝑚𝑣 𝑐 2 ≤ 1 2 𝑚𝑐 2 . Here, 𝐸 = 1 2 𝑚𝑐 2 is the upper limit of energy of motion.

Conclusion

In this article, a novel mass-energy equation has been derived using the Lorentz factor and the equation of kinetic energy in classical mechanics. This idea can enable the researchers for further involvements in the scientific research and development.

Electronic copy available at: https://ssrn.com/abstract=4629897


Nicolas Larrousse

: Huma-Num

Olivier Baude

Stéphane Pouyllau

Hélène Jouguet

Huma-Num, 10 years of building an infrastructure in SSH: a revolution in practices or an evolution?

Happy birthday Huma-Num! After 10 years of building an infrastructure dedicated to SSH in France: it's time to look back and see how far we have come.

Birth

At the beginning, around 2010, a time when digital services were only in their infancy, two existing infrastructures handled their use: the TGE Adonis, which stands for "Très Grand Equipement", and the Research Infrastructure "Corpus".

The first one was dedicated to providing digital services and the second one was designed to support SSH communities in digitizing corpora. In 2013, Huma-Num was created as a national infrastructure by the ministry of research and the CNRS by merging these two structures. The initial idea was to develop technological services for digital research data in close collaboration with research communities.

But what exactly is an infrastructure for SSH? This was a concept that had to be invented.

From Childhood to Adolescence

At the time of its creation, the most urgent task was to try to save data that were at best on a hard disk and sometimes on a floppy disk. We quickly discovered that preserving raw data was not sufficient -what was the point of preserving data if it was impossible to understand them? We had already gained some experience of this issue through long term preservation projects i but, as SSH are made up of diverse disciplines, we needed the support of experts. This is where Huma-Num consortia, the second part of Huma-Num inherited from the RI "Corpus", took over. A Huma-Num consortium is a group of people funded by Huma-Num, often from different institutions, sometimes from different disciplines, working on the same scientific objects, methods and themes. They provide expertise about their specific scientific objects in order to define good practices, training and tools ii .

It is an illustration of the fact that an infrastructure cannot be reduced to providing technological services, but that the expertise of communities and their involvement in the construction are key to building an infrastructure.

As for the technology, in 10 years we were well positioned to monitor the rapid evolution of needs. SSH research projects nowadays use more sophisticated tools (e.g. Geographical Information System, Rstudio, Deep Learning software, etc.), new types of data (e.g. Social media, LIDAR, 3D objects, etc.) and above all the volume of data has evolved very rapidly.

Toward maturity

Huma-Num has reached a certain level of maturity from a technological point of view even if there is still room for improvement as the needs and technologies are evolving rapidly: from a quantitative point of view, Huma-Num hosts more than 1 petabyte of data, about 400 virtual machines and 800 web sites. More generally, we are able to provide services for each step in the research data lifecycle. In particular, two services have been developed by Huma-Num to process, store and display research data while making them FAIR: NAKALA iii to share data and ISIDORE iv to disseminate them. Furthermore, Huma-Num's consortia are constantly changing and evolving to take into account new needs, for example supporting a standard such as the Music Encoding Initiative for the "Musica consortium" or providing a toolbox for historicized data for the "Paris Time Machine consortium".

At this point, we were able to focus on addressing one of the main issues we face: the quality of data, both of the digital object itself and its metadata. In order to address this issue, we can of course rely on the support provided by Huma-Num's consortia, but we need to have other complementary means to be able to support research projects throughout the whole data life cycle. Two years ago, we set up a dedicated service internally, currently made up of three people but which is expected to grow: they provide documentation, training and personalized support for research projects. This is of course not enough to take into account all needs and Huma-Num has established partnerships with regional SSH infrastructures, the "Maisons des Sciences de l'Homme v ", in order to provide additional support to local projects and thus fluidify the relationship with Huma-Num by clearing the way.

In the same vein, as a national infrastructure, it is important to establish contacts abroad to remain well-informed and also to be proactive in promoting services and standards developed by France. For the last 7 years, we have been establishing collaboration at the international level notably by participating in the construction of European Infrastructures such as DARIAH vi and participating in European infrastructure projects. We also have developed more targeted cooperation with Canada, for instance on innovative tools.

What's next?

We now have the necessary building blocks to develop new services and partnerships. The main one, which started this year, aims to establish a collaboration to link data and publication with other national infrastructures dedicated to publications, OpenEdition vii and METOPES viii . This ambitious project, called COMMONS, will last 8 years for a total cost of around 30 M€.

During the past 10 years, about 10,000 users have used the infrastructure, so the following question can be asked: Has Huma-Num influenced the ways of doing research in SSH or has it simply accompanied the general evolution of Digital Humanities? The Huma-Num lab was created in part to get an idea of the answers to these questions by investigating the evolution of current research practices in SSH through the prism of digital technology. This paper will present in greater detail what we learned the last 10 years the perspectives for the coming years.