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., [5,17,19,21,22,47] 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 [36]; biological and artificial neural networks in neuroscience and computer science [32]; gene regulatory networks in biology [26]; distributed resource allocation in operations research [57]; and opinion formation of individuals in sociology [9]; 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 [44,41,45,52].In such contexts, linear, and more generally, smooth approaches, (such as the master-stability function and contraction theory, see e.g., [4,50,41]), 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, [34,28,39].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 [12,16,28,30,34,39,55].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 [28,34,39,55] that under suitable assumptions (such as connectivity of the network, 1 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, [12,30] showed, under similar assumptions as in the diffusive case, that a better performance is obtained with high-gain time-varying interconnections.In [30] 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 [12] 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, [16] established interesting properties of the network that are not shared by their smooth counterparts.Indeed, [16] 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 [16] 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., [31,33,51].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 [14,15,24,27,42,49,54,56]. 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 [43,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 [30], adaptable dead-zone coupling [12], constrained differential inclusions [40], and perturbed Moreau's sweeping processes [35].Intuitively, it can be seen as a nonsmooth generalization of the funnel coupling in [30] and the dead-zone coupling in [12].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 [48,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 [30] 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 [10,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∥ = ∥x − proj(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∥ . ( Let φ : R l → R ∪ {+∞} be a proper, convex, lower semicontinuous function.The convex subdifferential of φ at x is defined as the set-valued map One of the most important properties of convex subdifferentials, that is central in the developments that follow, concerns its maximal monotonicity.A set-valued map 1 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 ) 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) The proximal map to φ is related to the convex subdifferential of φ via the following expression, see e.g., [7,Theorem 16.34], prox(φ; x) = (Id +∂φ) The right-hand side of (3) is also called the resolvent of ∂φ at x, see e.g., [7,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 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 Thus, it follows from (3) that proj(U; 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., [23], 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, 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 Throughout the paper it is assumed that each function ) 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, 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|} 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|} 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., [34,39].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.[16,33,51].
Motivated by the robustness properties of the set-valued hard-thresholds controllers studied in [33], it is proposed to use the following set-valued and time-dependent coupling law 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 where the sets C (i,j) are constant and satisfy the following standing assumption.
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 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), Note that with such initial condition, r (i,j) (t) > 0 for all t ∈ [0, +∞) and 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]) 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) (10a) The substitution of (10b) into the complementarity conditions (10c)-(10d) yields It is clear from (11) 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 j − y i ≤ r (i,j) . Hence, 1.If y j − y i = −r (i,j) , then it follows from (11b) and (10a) that I D2 = 0 and u (i,j) = −I D1 ≤ 0.
Putting all cases together we arrive at the expression 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 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.

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 (6) is written in compact form as Finally, the set S(r(t)) is such that 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 (13) 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 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 (13), see e.g., [11].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 |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 (15).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 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 We are now ready to state the main result concerning the existence of solutions of (35).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 [44], 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), where Thus, the derivative of the positive definite function 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 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 (13) 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, 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 (5).
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.
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 (5).
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 (6), 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., [1,6].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 [2], 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 where x(t) is the state of the network (13) at time t ≥ 0. The discretized output dynamics is The new variable ỹk+1 in (18) 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 [3,37].
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 where L = (Θ ⊺ ⊗ I m )Λ and Λ = ( C B) . 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, where 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 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 and the substitution of ( 22) back into (21b) gives us an explicit expression for v k+1 as from where (19) follows.Also note that (22) implies that wk+1 and (18a) leads us to If the map F and ζ k , are uniformly bounded then there is a finite M > 0 such that 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 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 [44], incremental dissipativeness [45], or QUAD-ness of the vector field [18], 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 (18) 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 [38] 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 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, where δ = 0.5.In this setting, the input to each vertex is computed in a centralized fashion via (19).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  (19).
Table 1: Nominal parameters for the three classes of systems for the network in Figure 3.
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 i − y j + Ru (i,j) ) .
Hence, y i − y j ∈ (Id +N(S (i,j) ; • ))(y i − y j + Ru (i,j) ), and from (4) we retrieve the explicit expression for u (i,j) as 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 (26) for the scalar case.
Figure 6: Time evolution of the graph of the coupling law in (26).In this case as t ↑ ∞, R t ↓ 0 and the graph of the coupling map approaches the graph of the normal cone 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 i − y 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 (12), see Theorem 6 and the remark after it.
Corollary 11.Under the assumptions of Theorem 6, the regularized coupling where u k+1 (i,j) is given by and the regularization parameters 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 Hence, ∥R k (i,j) u k+1 (i,j) ∥ = dist(y k i − y k j ; S(r k+1 (i,j) )).Consequently, 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 (28).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 (29) 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 (19) and (28) under the same sampling rate.As expected, there is a trade-off between centrality and since the centralized scheme (19) leads to a lower error compared to the fully distributed scheme (28).j∼i ∥y i (t k ) − y j (t k )∥ 2 for coupling schemes (19)(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 (13).The proof is an adaptation and an extension of the proof in [25] 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 where R satisfies Assumption 4. Thus, ( 13) is transformed into where where S(r(t is closed and convex.Hence, (31) transforms into where 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 η = Πη.

Construction of approximate solutions
Clearly, that well-posedness of ( 13) is equivalent to that of (35).Thus, we continue by discretizing in time the dynamics (35) using a semi-implicit Euler method.
where h := t k+1 − t k > 0 is the time-step, and It thus follows from the maximal monotonicity of the subdifferential map and Minty's theorem [7,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 [7,Theorem 3.14]), gives an explicit expression for z k+1 ⊥ satisfying (36a) as 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 (37) 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 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.
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.
To that end, let t ∈ [0, T ] \ {t 0 , t 1 , . . ., t N } and consider the following inequality Let us focus on the first term on the right-hand side of the inequality.It follows from ( 36)- (37) and the definition of projection that 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 and the substitution of ( 42) into (41) yields the estimate where L F denotes the Lipschitz constant of F and ∥ ζ(t k , z k )∥ ≤ Mζ ∥z k ∥ + M ζ follows from Assumption 1.
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 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 , 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 (33) "moves in a Lispchitz continuous way".
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 :
Figure 1: Ideal electrical circuit realizing the coupling (6) between systems ν i and ν j .

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

Figure 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 :
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(19).

Figure 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 :
Figure 8: Time evolution of the sum of squares error signal e SOS (t) := 1 2 |V| i=1