Hybrid coupling rules for leaderless heterogeneous oscillators: uniform global asymptotic and finite-time synchronization

We investigate the engineering scenario where the objective is to synchronize heterogeneous oscillators in a distributed fashion. The internal dynamics of each oscillator are general enough to capture their time-varying natural frequency as well as physical couplings and unknown bounded terms. A communication layer is set in place to allow the oscillators to exchange synchronizing coupling actions through a tree-like leaderless network. In particular, we present a class of hybrid coupling rules depending only on local information to ensure uniform global practical or asymptotic synchronization, which is impossible to obtain by using the Kuramoto model customarily used in the literature. We further show that the synchronization set can be made uniformly globally prescribed finite-time stable by selecting the coupling function to be discontinuous at the origin. Novel mathematical tools on non-pathological functions and set-valued Lie derivatives are developed to carry out the stability analysis. The effectiveness of the approach is illustrated in simulations where we apply our synchronizing hybrid coupling rules to models of power grids previously used in the literature.


Introduction
The Kuramoto model (Kuramoto (1975)) is used in various research fields to describe and analyze the dynamics of a broad family of systems with oscillatory behavior (Acebron et al. (2005)) including neuroscience (Aokii (2015); Tass (2003); Cumin and Unsworth (2007)), chemistry (Forrester (2015)), power networks (Dörfler and Bullo (2012)) and natural sciences (Leonard et al. (2012)), to cite a few (see also (Strogatz (2003))).The many application areas where Kuramoto dynamics emerged from physical considerations motivated a detailed analysis of the synchronization properties of the model, first for the all-to-all connection case (Aeyels and Rogge (2004)), as originally described by Kuramoto, then for a general interconnection layout (Jadbabaie et al. (2004)), with a focus on the derivation of the least conservative lower bound for a stabilizing coupling gain (Jafarpour and Bullo (2019); Chopra and Spong (2009); Dörfler and Bullo (2011)).
Given its simple and accurate description of natural synchronization phenomena, the Kuramoto model has also inspired the design of distributed communication protocols in engineering applications where the coupling function among different agents can be arbitrarily assigned to achieve synchronization, as in the bio-inspired synchronization of moving particles in (Sepulchre et al. (2007)), the synchronized acquisition of oceanographic data from Autonomous Underwater Vehicles (Baldoni et al. (2007)), in clock synchronization (Kiss (2018)), in mobile sensors networks modeled as particles with coupled oscillator dynamics (Paley et al. (2007)), in monotone coupled oscillators (Mauroy and Sepulchre (2012)) or in other engineering applications surveyed in (Dörfler and Bullo (2014)).
While the sinusoidal coupling of Kuramoto models provides a powerful tool to obtain synchronization in coupled networks of oscillators, it also introduces some undesirable properties for engineering applications.For example, when the network comprises oscillators with the same natural frequency, it is now well-known that a system of Kuramoto oscillators admits, in addition to stable equilibria coinciding with the synchronization set, equilibria that are unstable (see, e.g., (Strogatz (2000); Sepulchre et al. (2007))).The downside of this result is that the closer a solution is initialized to an unstable equilibrium, the longer it will take for phase synchronization to arise: we talk of non-uniform convergence (Sepulchre et al. (2007)).Although non-uniform synchronization may naturally characterize certain physical (Oud (2006)) and biological systems, in general, it is not a desirable property for engineering applications.Indeed, the lack of uniformity may induce arbitrarily slow convergence to the attractor set and poor robustness properties (Miller and Pachter (1997)).Secondly, it may occur in the Kuramoto model that the angular phase mismatch between adjacent oscillators remains constant and different from zero indefinitely: in this case we talk of phase locking (Aeyels and Rogge (2004)), which hampers the capability to reach asymptotic collective synchronization.Thirdly, in critical applications, finite-time stability, instead of only asymptotic synchronization, may be a mandatory requirement (Polyakov (2011)).
In this work, we investigate the engineering scenario where the goal is to synthesize local coupling rules to synchronize a set of heterogeneous oscillators.We assume the model of the oscillators to be general enough to capture not only their (time-varying) natural frequency but also physical coupling actions and other unknown bounded terms, thus being able to represent, among many possibilities, networks of Kuramoto oscillators with heterogeneous time-varying natural frequencies.Furthermore, without loss of generality, we introduce suitable resets of the oscillators' phase coordinates, so that they are unwrapped to evolve in a compact set, which includes [−π, π] consistently with their angular nature.Consequently, we define hybrid 2π-unwinding mechanisms to ensure the forward completeness of the oscillating solutions.
To achieve uniform global phase synchronization, thereby overcoming the limitations of Kuramoto models, we equip the oscillators with a leaderless tree-like communication network to locally exchange coupling actions based on local information.This approach has been already exploited in the context of DC microgrids as in, e.g., (Cucuzzella et al. (2018)), or (Giraldo et al. (2019)), for a network of Kuramoto oscillators equipped with a leader.The selection of a tree-like graph, which can always be derived in a distributed way by using the algorithms surveyed in (Pandurangan et al. (2018)), is also not new while addressing a problem of distributed cooperative control: see (Mayhew et al. (2012b)) in the context of hybrid dynamical systems, or (Bai et al. (2011)) and (Alagoz et al. (2012)) for continuous-time networked systems and power grids, respectively.To define the coupling actions, we present novel hybrid coupling rules for which a Lyapunov-based analysis ensures uniform global (practical or asymptotic) phase synchronization.This result overcomes both the lack of uniform convergence and the phase-locking issues characterizing the Kuramoto model (Sepulchre et al. (2007)).Interestingly, we can design the coupling rules in such a way that the network of oscillators behaves like the original Kuramoto models when the oscillators are near phase synchronization.Furthermore, due to the mild properties that we require for our hybrid coupling function, discontinuous selections are allowed, like in (Coraggio et al. (2020)).When the discontinuity is at the origin, we prove finite-time stability properties.In particular, exact synchronization can be reached in a prescribed finite-time (Song et al. (2017)), and convergence is thus independent of the initial conditions.Compared to the related works in (Mauroy and Sepulchre (2012)) and (Wu and Li (2018)), the finite-time stability property we ensure is global and the convergence time can be arbitrarily prescribed, respectively.We resort for this purpose to non-smooth Lyapunov theory, in particular non-pathological Lyapunov functions and setvalued Lie derivatives (Bacciotti and Ceragioli (2003)), for which we provide new results and novel proof techniques that are of independent interest.Due to the possible presence of discontinuities in the coupling function, the stability analysis is carried out by focusing on the regularization of the dynamics, as typically done in the hybrid formalism of (Goebel et al., 2012, Ch. 4).Finally, simulations are provided to illustrate the theoretical guarantees and demonstrate the potential strength of our hybrid theoretical tools to address both first and second-order oscillators modeling generators in power grids considered in (Dörfler and Bullo (2012)).
The recent submission (Bosso et al. (2021a)) (see also (Bosso et al. (2021b))) also uses hybrid tools to obtain uniform global synchronization guarantees in a Kuramoto setting but in a different context, namely for second-order oscillators (where the ω i 's are states rather than external inputs) and, most importantly, for a network with a leader, which significantly changes the setting compared to the leaderless scenario investigated in this work, where no oscillator is insensible to the coupling actions from its neighbours.With respect to the preliminary version of this work in (Bertollo et al. ( 2020)), we include the next novel elements: relaxed requirements on the coupling function, time-varying, phasedependent, (possibly) non-identical natural frequencies, generalizing the two-agents theorems of (Bertollo et al. (2020)) to the case of n oscillators in addition to establishing a set of new stability results missing in (Bertollo et al. (2020)) (finitetime, practical properties and other ancillary results).
The rest of the paper is organized as follows.Notation and background material are given in Section 2. The local hybrid coupling rules and oscillators network model are derived in Section 3. In Section 4, we introduce the regularized version of the dynamics presented in Section 3. In Section 5, we present Lyapunov-based analysis tools establishing the asymptotic properties of our model, while prescribed finite-time results are given in Section 6. Numerical illustrations are provided in Section 7, while most of the technical aspects of our proofs requiring non-smooth analysis concepts are gathered in Section 8.A few proofs of minor importance are relegated to the Appendix.
Background on graph theory.We denote an unweighted undirected graph as u = ( , u ), where is the set of vertices, or nodes, and u ⊆ × is the set of edges, or arcs, composed by unordered pairs of nodes.If a pair (i, j) of nodes belongs to u , we say that those nodes are adjacent and that j is a neighbour of i and vice versa.Given two nodes x and y of an undirected graph u , we define as path from x to y a set of vertices starting with x and ending with y, such that consecutive vertices are adjacent.If there is a path between any couple of nodes, the graph is called connected, otherwise it is called disconnected.We define as subgraph of u a graph s = ( s , s ), where s ⊂ and s ⊂ u .An induced subgraph of u that is maximal, subject to be connected, is called a connected component of u .A cycle is a connected graph where every vertex has exactly two neighbours.An acyclic graph is a graph for which no subgraph is a cycle.A connected acyclic graph is called a tree.
We denote an unweighted directed graph as = ( , ), where ⊆ × is composed of ordered pairs, therefore arcs have a specific direction.An arc going from node i to node j is denoted by (i, j) ∈ .If a directed graph is obtained choosing an arbitrary direction for the edges of an undirected graph u , we call it an oriented graph, and we say that is obtained from an orientation of u .If (i, j) ∈ , we say that i belongs to the set of in-neighbors j of j, while j belongs to the set of out-neighbors i of i.The union of i and i gives the more generic set of neighbors i := i ∪ i of node i, containing all the nodes connected to it, in any direction.With B ∈ n×m we denote the incidence matrix of graph such that each column [B] , ∈ {1, . . ., m}, is associated to an edge (i, j) ∈ , and all entries of [B] are zero except for b i = −1 (the tail of edge ) and b j = 1 (the head of edge ), namely [B] = e j − e i .

Flow dynamics
Consider a networked system of n heterogeneous oscillators.To achieve synchronization, the oscillators locally exchange coupling actions through the unweighted undirected tree1 u := ( , u ) made of n nodes and thus m = n − 1 edges, n ∈ >1 .We assign an arbitrary orientation to u , which leads to the oriented tree = ( , ).In this scenario, the oscillator phase corresponding to node i, with i ∈ , is denoted θ i and has the next flow dynamics where ω i (θ , t) is a possibly an unknown term modeling the dynamics of the i-th oscillator, which can capture physical coupling actions, its time-varying natural frequency, and any other unknown bounded dynamics affecting the oscillator; see Section 7 for a numerical example.We assume that ω i is locally bounded, measurable in t, piecewise continuous in θ and such that ω i (θ , t) ∈ Ω := [ω m , ω M ] for any time t ≥ 0 and (θ , q) ∈ C, with ω m ≤ ω M ∈ , namely Ω is a compact interval of values2 .Since (1) possibly has a discontinuous right-hand side, the notion of solution should be carefully defined, and we postpone this discussion to Section 4 (where we also prove the existence of solutions) to avoid overloading the exposition.For now it suffices to say that a function θ is a solution of (1) if it is absolutely continuous (i.e., it coincides with the integral of its derivative) and satisfies (1) almost everywhere.
Phase θ i in (1) evolves in the set [−π − δ, π + δ], with δ ∈ (0, π), which thus covers the unit circle corresponding to phases taking values in [−π, π].Parameter δ > 0 inflates the set of angles [−π, π] to rule out Zeno solutions as explained in the following, see Section 3.3.Thus, δ is a regularization parameter chosen to be the same for each oscillator.Variable q i j , with (i, j) ∈ , is a logic state taking values in {−1, 0, 1}, which is constant during flows.Its role is to unwind the difference between the two phases θ j and θ i through jumps.Indeed, since θ j and θ i are angles, to evaluate their mismatch, loosely speaking, we have to consider their minimum mismatch modulo 2π: q i j is introduced for this purpose as clarified in Section 3.2.The vectors θ and q collect all the states θ i , i ∈ , and q i j , (i, j) ∈ , respectively, as formalized in the following, together with the formal definition of the flow set C, namely a compact subset of the state-space where the solutions are allowed to evolve continuously.The gain κ ∈ >0 is associated with the intensity of each coupling action and it is the same for each interconnection.Finally, the coupling action between each pair of nodes (i, j) ∈ is defined as σ(θ j − θ i + 2q i j π), where σ is the function used to penalize the phase mismatch θ j − θ i + 2q i j π between phases θ j and θ i , and it satisfies the next property.
Item a) of Property 1 ensures that σ is an odd function and thus implies σ(0) = 0, while item b) of Property 1 guarantees that σ(s) can only be zero at s = 0. Notice that the sine function, customarily used in the classical Kuramoto model, satisfies item a) but not item b) of Property 1, which is fundamental to establish the global uniform stability result of this work.Examples of functions σ satisfying Property 1 are depicted in Figure 1, together with the sine function for the sake of comparison.We emphasize that the mild assumptions of Property 1 allow considering, among others, intuitive discontinuous selections such as the sign function of Figure 1, which leads to an interesting parallel between (1) and the ternary controllers considered in (De Persis and Frasca (2013)).Another possible example of σ enjoying Property 1 is σ(s) = sin(s) + u(s), where u is such that items a) and b) of Property 1 hold.Note that when u is negligible compared to sin in a neighborhood of the origin, the model behaves locally like the classical Kuramoto network.Also, Property 1 comes with no loss of generality as we consider the scenario where we have the freedom to design the coupling rules among the oscillators and thus σ.
Collecting in the vector σ(x) ∈ m all the coupling actions σ(θ j − θ i + 2q i j π), with (i, j) ∈ , using the same order as the columns of B, the flow dynamics in ( 1) is written as , and where q ∈ {−1, 0, 1} m is the vector stacking all the q i j 's for (i, j) ∈ , ordered as in σ(x).Thus, the overall state x := (θ , q) evolves in the compact state space defined as The flow set C in (2) will be selected as the closed complement of the jump set D introduced next.

Jump dynamics
We introduce jump rules to constrain each phase θ i to take values in [−π − δ, π + δ] as well as to guarantee that the argument θ j −θ i +2q i j π of σ in (1) belongs to dom σ = [−π− δ, π + δ] when flowing.To guarantee the latter property, define, for any (i, j) ∈ , the jump set and the associated difference inclusion where the entries of G i j : X ⇒ {−1, 0, 1} m are given by with (u, v), (i, j) ∈ .Set D i j in (4a) enforces a jump when θ j − θ i + 2q i j π is not in dom σ for (i, j) ∈ .Across a jump, according to (4b), only q i j changes in such a way that |θ j − θ i + 2q i j π| < π + δ after a jump as formalized in the next lemma whose proof is given in Appendix A to avoid breaking the flow of the exposition.
Figure 2: Projection of the flow and jump sets on (θ i , θ j ) for each value of q i j .
Lemma 1.For any (i, j) ∈ and x A second jump rule is introduced for when one of the oscillators i ∈ reaches |θ i | = π+δ.In this case, a jump of 2π is enforced so that the phase then belongs to (−π − δ, π + δ) while remaining the same modulo 2π.We define for this purpose where the entries of g i,θ : X → [−π − δ, π + δ] n and g i,q : X → {−1, 0, 1} m are defined as with j ∈ and (u, v) ∈ .The set D i , i ∈ , is defined as and In view of (5d), the jump rule (5a) is allowed when both |θ i | = π + δ and x is not in the interior of D uv for any (u, v) ∈ , where a jump may occur according to (4).Note that each function g i is continuous on its (not connected) domain because D i does not contain points with θ i = 0 for any i ∈ .
Finally, switching/jumping ruled by (5) unwinds the phase θ i without changing the phase mismatches between neighbours, defined as (θ j −θ i +2q i j ), as shown in the next lemma, whose proof is given in Appendix A. (6)

Overall model
In view of Sections 3.1-B, the overall hybrid model is given by ẋ where f is defined in (2), and using (4a) and (5d), with X defined in (3).The set-valued jump map G is defined in terms of its graph, which is given by with g i and G ext i j as per (4b), ( 5a)-(5c).Figure 2 shows three projections of the state space X on the plane (θ i , θ j ) for some (i, j) ∈ , which corresponds to a union of three squares, one for each value of q i j .
Remark 1.Since we envision engineering applications, each phase θ i with i ∈ may be reconstructed from the angular measurements provided by sensors.Due to the wide variety of outputs provided by commercial sensors, a relevant task is to extrapolate a continuous measurement from a sensor that may return values whose wrapping around 2π is unknown; see, for example, (Reigosa et al. (2018)) and (Anandan and George ( 2017)).In this scenario, we can implement an algorithm to extract a continuous measurement of the phase satisfying (2).In particular, following a rationale similar to that proposed in (Mayhew et al., 2012a, Figure 1) for a setting with sampled measurements, we may continuously update an estimate θ i,e of θ i .Indeed, for each sensor output θ i,so , we may extract the lifted measurement as the closest one to θ i,e when performing 2π-wraps θ + i,e = θ i,so + 2π argmin h∈{−1,0,1} |θ i,so − θ i,e + 2hπ|.This rule parallels the selection of ( 15) and (27b) of (Mayhew et al. (2012a)) for the simpler case of 1 and scalar angular measurements.

Regularized hybrid dynamics
Model ( 7) is a time-varying hybrid system with a possibly discontinuous right-hand side, due to the mild properties of σ, see Property 1.Hence, solutions may be understood in the generalized sense of (Goebel et al. ( 2012)).We consider for this purpose the regularization of ( 7), so that stability properties for the regularized system carry over to the nominal and generalized solutions of (7).In particular, following (Goebel et al., 2012, Page 79) we consider where C, D, and G coincide with those in (7), and the setvalued map F regularizes f in (2) as with the sets is the stacking (with the same ordering as in σ) of the setvalued maps σ i j defined as for all (i, j) ∈ .Since the jump set, flow set, and jump map of hybrid system (8) coincide with those of (7), and for any x ∈ X and ω ∈ Ω, f (x, ω) ∈ F (x), we study the stability properties of solutions of (7) by concentrating on the regularized dynamics (8).In addition to clarifying the nature of solutions of (7), which may, among other things, present sliding behavior (see Section 6), the advantage of using (8) instead of ( 7) is that (8) satisfies the so-called hybrid basic conditions (Goebel et al., 2012, As. 6.5), which ensure its wellposedness (Goebel et al., 2012, Thm. 6.30).
Proof: Sets C and D, as defined in (4a), (5d), (7b), (7c) are closed, as required by (Goebel et al., 2012, As. 6.5 (A1)).On the other hand, F is the Krasovskii regularization of a function, which satisfies the HBC in view of its locally boundedness on C and (Goebel et al., 2012, Lemma 5.16) as shown in (Goebel et al., 2012, Ex. 6.6), thus (Goebel et al., 2012, As. 6.5 (A2)) is satisfied.Lastly, each g i and G ext i j has a closed graph, and so due to (7d) the graph of G is closed as well.As consequence, according to (Goebel et al., 2012, Lemma 5.10), G is outer semicontinuous and it is also locally bounded relative to D, thereby satisfying (Goebel et al., 2012, As. 6.5 (A3)).
Among other useful properties, Lemma 3 guarantees intrinsic robustness of the stability property established later in Sections 5 and 6, see (Goebel et al., 2012, Ch. 7).To conclude this section, we note that all maximal solutions to (8) are complete3 and exhibit a (uniform) average dwell-time property, thereby excluding Zeno phenomena.We emphasize that, through Lemmas 1 and 2, the parameter δ plays a key role in establishing that no complete discrete solution exists.In particular, the fact that δ = 0 and δ = π is key for being able to exclude Zeno solutions.
Proposition 1.All solutions to (8) enjoy a uniform average dwell-time property.Namely, there exist τ D ∈ >0 and J 0 ∈ ≥0 such that, for any solution x to (8) and for any pair of hybrid times such that t + j ≥ s + r with (s, r), (t, j) ∈ dom x, Proof: We first recall that, in view of Lemma 2, for any i ∈ , x ∈ D i and x + = g i (x) while the other θ j , j = i ∈ remain unchanged across such a jump and so does θ j − θ i + 2q i j π for all (i, j) ∈ .We also recall that, from Lemma 1, for any (i, j) ∈ , x ∈ D i j and while the θ u and the other θ h −θ k +2q hk π remain unchanged for any u ∈ and (h, k) = (i, j) ∈ .From uniform global boundedness of the right hand-side F (x) of the flow dynamics (a consequence of the local boundedness of F and of the boundedness of X ), all solutions satisfy a global Lipschitz property with respect to the flowing time and ( 11) and ( 12) imply a uniform average dwell time on the jumps from n i=1 D i and (i, j)∈ D i j , respectively.Finally, the uniform average dwell time property of solutions jumping from D derives directly from Lemma 2 and (4b)-(4c).Indeed, (4b)-(4c) imply that jumping from D i j does not affect the triggering condition in D u or D hk , for any u ∈ and (i, j) = (h, k) ∈ .In a similar way, Lemma 2 implies that jumping from D i does not affect the triggering condition in D j or D i j or D ji , for any j = i ∈ and (i, j) ∈ or ( j, i) ∈ .Thus, a uniform average dwell time on jumps from D stems from the global Lipschitz property of the solutions with respect to the flowing time, along with ( 11) and ( 12).
We now check that maximal solutions of (8) are complete by proving that the conditions of (Goebel et al., 2012, Prop. 6.10) (Goebel et al., 2012, Prop. 6.10 (VC)) holds for any ξ ∈ C \ D. On the other hand, the state space X in (3) is bounded, thus item (b) in (Goebel et al., 2012, Prop. 6.10) is excluded.To rule out item (c) in (Goebel et al., 2012, Prop. 6.10), from Lemmas 1 and 2, we have G(D) ⊂ X = C ∪ D. Hence, we can apply (Goebel et al., 2012, Prop. 6.10) to conclude that all maximal solutions are complete, thus obtaining t-completeness of solutions in view of their uniform average dwell time property established above.

Synchronization set and its stability property
To analyze the synchronization properties of system (8), consider the set Because the network is a tree, for any x ∈ , the phases θ i and θ j coincide modulo 2π not only for any (i, j) ∈ but also for any i ∈ and j ∈ \ {i}.In other words, when x ∈ , all the oscillators are synchronized even if they do not share a direct link.We therefore call the synchronization set.Our main result below establishes a practical asymptotic stability result for , as a function of the coupling gain κ appearing in the flow map (1).The "practical" tuning of κ depends on the following two parameters: Parameter λ ensures a detectability property of the distance |x| in (13) from the norm | σ|, for any σ ∈ Σ(x).In particular, we have from the results in (Godsil and Royle (2001)).and τ i 0, such that w = lim Lemma 4. Since is a tree, λ := λ min (B B) > 0.
Proof: Consider a tree graph with incidence matrix B ∈ n×n−1 .By definition has only one bipartite connected component.From (Godsil and Royle, 2001, Thm. 8.2.1), rank(B) = n − 1 and thus dim(null(B)) = 0 by way of the fundamental theorem of linear algebra.Therefore, B y = 0 for any y ∈ n−1 \{0 n−1 } which implies y B B y = |B y| 2 > 0 for any y ∈ n−1 \ {0 n−1 }.Consequently all the eigenvalues of B B are strictly positive, thus completing the proof.
Remark 2. The smallest eigenvalue λ of B B and its positivity established in Lemma 4 play a fundamental role on the speed of convergence of the closed-loop solutions to the synchronization set.As is a tree, positivity of λ is ensured by Lemma 4. In more general cases with not being a tree, the leaderless context considered in this paper, where the synchronized motion emerges from the network, poses significant obstructions to achieving global results.A simple insightful example of a cyclic graph is discussed in Section 5.2, which provides a clear illustration of the motivation behind requiring that is a tree.We emphasize that a similar obstruction is experienced in prior work (Mayhew et al. (2012b)) where, in a different context, a similar assumption on the network is required.
We are now ready to state the main result of this paper, corresponding to a practical bound on the distance of x from that is uniform in κ.We state the bound in our main theorem below, whose proof is given in Section 8.2, and then illustrate its relevance on a number of corollaries given next.
Theorem 1.Given set in (13), there exists a class function β • and a class gain γ • , both of them independent of κ, such that, for any κ > 0, all solutions x of (8) satisfy for all (t, j) ∈ dom x and with c := max s∈dom σ σ(s).
The bound (15) in Theorem 1 is the sum of two terms: β • captures the phases tendency to synchronize, while function γ • depends on the mismatch among the (possibly) nonidentical, time-varying natural frequencies of the oscillators, which hampers asymptotic phase synchronization in general.Therefore, Theorem 1 provides an insightful bound (15) illustrating the trend of the continuous-time evolution of the hybrid solutions to (8).Notice that β • and γ • can be constructed by following similar steps as the ones in (Sontag and Wang, 1995, Lemma 2.14), noting that the resulting bound is often subject to some conservatism.On the other hand, because β • and γ • are independent of κ and λ, (15) still provides valuable quantitative information.Indeed, in view of (15), increasing κ speeds up the transient and reduces the asymptotic phase disagreement caused by the non-identical time-varying natural frequencies.Equation ( 15) also highlights the impact of the algebraic connectivity λ of (Mesbah and Egerstedt, 2010, pages 23-24) on the phase synchronization, by giving information on the scalability of our algorithm.Recall that λ is influenced by several parameters of the undirected graph, such as the maximum degree and the number of nodes (Rad et al. (2011)).This continuous-time focus in (15) in Theorem 1 is motivated by Proposition 1.It is also of interest to establish a bound similar to (15) while measuring the elapsed time in terms of t + j and not only in terms of t, as usually done when defining bounds for solutions to hybrid systems, which allows us to ensure stronger stability properties, in particular, uniformity and robustness (Goebel et al., 2012, Chp.7).Hence, combining Theorem 1 with Proposition 1, we obtain the following second main result.
Theorem 2. For each value κ > 0, there exists a class function β such that all solutions to (8) satisfy for all (t, j) ∈ dom x, with γ • as in Theorem 1.
Theorem 2 implies that the oscillator phases uniformly converge to any desired neighborhood of by taking κ sufficiently large, thus the practical nature of the result.We also immediately conclude from Theorem 2 and Lemma 3 that the stability property in ( 16) is robust in the sense of item (a) of (Goebel et al., 2012, Def. 7.18), according to (Goebel et al., 2012, Thm. 7.21).
We may draw an important additional conclusion from Theorem 2 corresponding to a global practical bound stemming from the fact that the function γ • in ( 16) is independent of κ.
Lastly, in the case of uniform frequencies ω(θ , t) = 1 n ω(t), for all t ≥ 0, with ω(t) ∈ Ω, we have B ω(t) = 0, for all t ≥ 0.Then, we can exploit the fact that the term |ω| at the right-hand side of ( 16) stems from upper bounding |B ω| 1 as in ( 14), which allows obtaining the following asymptotic property of .

Cyclic graphs and their potential issues
Before proceeding with the technical derivations needed to prove Theorem 1, we devote some attention to the issues pointed out in Remark 2 about the need for the graph to be a tree, similar to (Mayhew et al. (2012b)).
More generally, when the graph is not a tree, the kernel of matrix B contains additional elements besides the zero vector.Consequently, we can have Bσ(x) = 0 n even when x / ∈ , and λ = 0 in ( 14).As a result, is not globally attractive.
Function V enjoys useful relations with the distance of x from the synchronization set , as formalized next.
To prove Theorem 1, it is also fundamental to formalize the relation between the distance of x from the set and Σ(x) in ( 8), as done in the next lemma.Lemma 6.There exists a class ∞ function η such that, for each x ∈ X , η(|x| ) ≤ | σ| 2 , ∀ σ ∈ Σ(x).
Proof: From items a) and b) of Property 1, |σ(s)| ≥ α(|s|) for any s ∈ dom σ.Thus, in view of ( 9), |ς| ≥ α(|s|), for any s ∈ dom σ and ς ∈ σ(s).We recall that, for any x ∈ X , Σ(x) is the stacking of all the set-valued maps σ i j defined in ( 10), (i, j) ∈ .Hence, by definition of , for any x ∈ X \ , there exists at least one element θi j = 0, with (i, j) ∈ , thus ) is a suitable lower bound for | σ|, for any x ∈ X .Since θi j is a function of the states and max (i, j)∈ α(•) is positive definite and radially unbounded, as X is compact, then (Goebel et al., 2012, Page 54) implies that there exists η ∈ ∞ such that η(|x| ) ≤ | σ| 2 holds for each x ∈ X and for all σ ∈ Σ(x), thus concluding the proof.
Function V is locally Lipschitz due to the properties of σ and characterizing its variation when evaluated along the solutions of the hybrid inclusion (8) requires using tools from non-smooth analysis.To avoid breaking the flow of the exposition, we postpone to Section 8.2 those technical derivations and summarize the corresponding conclusions in the next proposition, a key result for proving Theorem 1.

Prescribed finite-time stability properties
A useful outcome of the mild regularity conditions that we require from σ in Property 1 is that defining σ to be discontinuous at the origin, as in the sign function represented in Figure 1, leads to desirable sliding-like behavior of the solutions in the attractor .This sliding property induces interesting advantages of the behavior of solutions, as compared to the general asymptotic and practical properties characterized in Section 5.
A first advantage is that, even with non-uniform natural frequencies, we prove uniform global asymptotic stability of for a large enough coupling gain κ, due to the well-known ability of sliding-mode mechanisms to dominate unknown additive bounded disturbances acting on the dynamics.A second advantage is that the Lyapunov decrease characterized in Proposition 2 can be associated with a guaranteed constant negative upper bound outside , which implies finite-time convergence.Finally, since this constant upper bound can be made arbitrarily negative by taking κ sufficiently large, we actually prove prescribed finite-time convergence (see (Song et al. (2017))) when using these special discontinuous functions σ, whose peculiar features are characterized in the next lemma.
Paralleling the structure of Proposition 2, the next proposition, whose proof is postponed to Section 8.3, is a key result for proving Theorem 3. Proposition 3. Consider system (8) and function V in (19)-( 20).If σ is discontinuous at the origin, then there exist µ ∈ >0 independent of ω in (14) and κ > 0 such that for each κ ≥ κ any solution x of (8) satisfies (denoting (ii) for all j ∈ {0, . . ., J − 1}, Exploiting Lemma 7 and Proposition 3, we can follow similar steps to those in the proof of Theorem 1 to show the following main result on uniform global asymptotic stability and prescribed finite-time stability of for ( 8).
Theorem 3. If σ is discontinuous at the origin, then set in (13) is prescribed finite-flowing-time stable for (8), i.e., for each T > 0 there exists κ > 0 such that for each κ ≥ κ : (i) there exists β ∈ such that all solutions x satisfy (18); (ii) all solutions x satisfy, x(t, j) ∈ for all (t, j) ∈ dom x with t ≥ T .

Proof:
We start showing that G(D∩ ) ⊂ .Indeed, we notice that D i j ∩ = for any (i, j) ∈ , and thus G(D i j ∩ ) ⊂ trivially holds.Moreover, from Lemma 2, it holds that G(D i ∩ ) ⊂ for any i ∈ .Hence, from (7b), we conclude that G(D ∩ ) ⊂ .To establish that is (strongly) forward invariant for (8), it is left to prove that solutions cannot leave while flowing.We proceed by contradiction and for this purpose suppose there exists a solution x bad to (8) such that x bad (0, 0) ∈ and x bad (t * , 0) / ∈ for some t * > 0 with (t * , 0) ∈ dom x bad .From continuity of flowing solutions between any two successive jumps and closedness of , there exists x bad (t, 0) ∈ for all t ∈ [0, t * ) and x bad (t * , 0) / ∈ .Hence, from ( 19) and ( 20) and positive definiteness of V , we have 0 = V (x bad (0, 0)) < V (x bad (t * , 0)).However, since the solution is flowing, integrating (24a) over the continuous time interval [0, t * ] we obtain V (x bad (t * , 0)) < V (x bad (0, 0)), which establishes a contradiction.Consequently, a solution cannot leave while flowing.We have proven that the set is (strongly) forward invariant, implying that if x(t, j) ∈ X \ then x(t , j ) ∈ X \ , for any t + j ≤ t + j, with (t , j ), (t, j) ∈ dom x.Let κ ≥ κ with κ defined in Proposition 3 and x be a solution to (8).Combining (24a) with the non-increase condition (24b) and the forward invariance of , we obtain by integration for any (t, j) whenever x(t, j) ∈ X \ , and thus for any x(t, j) ∈ X .Equation ( 26) can then be converted to a bound on |x(t, j)| using ( 21).Hence, we follow the same steps used in the proofs of Theorem 1 and 2 to obtain (18), thus concluding the proof of item (i) in Theorem 3.
In view of ( 25) and from the positive definiteness of V with respect to , we conclude that solutions to ( 8 Notice that phases synchronize at most at continuous time T = 1 κ 2v λµ 2 , in view of ( 26).Therefore, we may decrease T at will by selecting a larger µ as in Lemma 7 and/or by increasing the coupling gain κ.  27), and the red graph is the communication tree graph that we design for the hybrid coupling rules presented in Section 3.

Numerical illustration
In this section, we apply our control scheme to globally, uniformly, synchronize the phases of n = 10 strongly damped generators physically coupled over an all-to-all network connection (represented by the blue edges in Figure 3) over a set of nodes whose dynamics are approximated by nonuniform first-order Kuramoto oscillators given by (Dörfler and Bullo, 2012, eq. (2.8)).This fully connected dynamics can be accurately described by the terms ω i 's in (1) with the following selection: ).The physical all-to-all coupling among the oscillators (blue edges in Figure 3) is modeled by the sine functions sin(θ j − θ i + φ i j ) of the angular mismatch between the oscillators offsetted by the constant angle φ i j ∈ uni([0, atan(0.25)]).Furthermore, each physical coupling is scaled by the gain κ i j = κ ji ∈ uni([0.7,1.2]).This allows fully embedding in our time-varying generalized natural frequencies ω i of (1) the physical couplings of the oscillators.Each high-frequency disturbance d i : ≥0 → [0, 5] is defined as d i (s) = 0 if s ∈ [0, 5.2] ∪ [6.0, 11] and d i (s) =5 sin(50χ i s + φ i ) if s ∈ (5.2, 6).Thus ω i (θ , t) not only captures the time-varying natural frequency of the i-th oscillator but also the physical coupling actions and disturbances influencing its dynamics.The parameters have been selected as in (Dörfler and Bullo (2012)) to model realistic, strongly damped, generators.On the other hand, the synchronizing coupling actions are exchanged through our communication graph u = ( , u ), whose edges are depicted in red in Figure 3.These "cyber" coupling actions are represented by the functions σ's in (1), whose design is performed according to our solution of Section 3. Summarizing, the combination of the (blue) physical layer and the (red) "cyber" communication layer of Figure 3 generates a cyber-physical system whose dynamics is represented by ( 2), (7), with ω capturing the physical layer and σ capturing the hybrid feedback control action.We initialize the oscillators with q(0, 0) = 0 9 and the initial phases are chosen in such a way that the oscillators are equally spaced on the unit circle.Finally, we select δ = π 4 .The evolution of the phases, θ i s, and the angular errors between any two neighbours in , namely min h∈{−1,0,1} (|θ j − θ i + 2hπ|), are reported 5 in the top two rows of Figure 4 and in Figure 5, for different selections of σ, and κ = 576π 10 , which ensure finite-time synchronization due to the bound reported in Lemma 10 in Section 8.3.When no communication layer is considered (left plots), the oscillators do not synchronize.When the communication layer is implemented and σ is given as the ramp function, practical synchronization is achieved as established in Theorem 1 and shown in Figures 4 and 5. (Non-uniform) practical synchronization can also be achieved in the absence of a communication layer by selecting larger values for the κ i j 's.On the other hand, the sign function, which is discontinuous at 0, also leads to a finite-time synchronization property in agreement with Theorem 3, see Figures 4 and 5. Furthermore, the set of plots in the bottom row of Figure 4 shows that each phase θ i maps continuously the angular values identifying oscillator i on the unit circle, in agreement with Section 3.1 and Lemma 2.
The same exact hybrid controller dynamics is finally exploited in a more sophisticated context of non-strongly damped generators (rather than the strongly damped case, as considered above and in Figures 4 and 5).Following (Dörfler and Bullo ( 2012)), such behavior is modeled by a fully connected graph comprising the second-order (rather than first-order) heterogeneous oscillators in (Dörfler and Bullo, 2012, eq. (2.3)).Once again, this physical interconnection is well represented by the blue edges in Figure 3 and dynamics (1) with the following selection, generalizing (27 This dynamically generalized selection uses the same parameters as in the previous set of simulations, with the addition of the constant mass parameters m i ∈ uni([2, 12] 1 120π ) for each i ∈ , which is defined according to (Dörfler and Bullo, 2012, eq. 2.3).We apply our hybrid feedback control algorithm to this generalized scenario by augmenting once again the physical layer with a "cyber" communication layer represented by the red edges in Figure 3, inducing the stabilizing action of inputs σ in the hybrid interconnection ( 2), (7).We initialize q and θ as in the previous simulations and θi (0, 0) ∈ uni([−0.1,0.1]) for each i ∈ .Finally, we still select δ = π 4 in (3).The evolution of the phases is reported in Figure 6, for different selections of σ, and κ = 576π 10 .Similarly to what happens for the first-order oscillators of Figures 4 and 5, when no communication layer is considered, the second-order oscillators do not synchronize.When the generators are equipped with the communication network and σ is instead defined as the ramp function, practical synchronization is achieved, as predicted by Theorem 1 and as shown in Figure 6.On the other hand, considering the (discontinuous at 0) sign function to generate the synchronizing hybrid coupling actions again leads to a finite-time synchronization property, thus confirming Theorem 3.

Results on scalar non-pathological functions
The proofs of Propositions 2 and 3, which are instrumental for proving our main results of Sections 5 and 6, require exploiting results from non-smooth analysis, because V in ( 19) is not differentiable everywhere.A further complication emerges from the fact that, since σ may be discontinuous, the flow map in dynamics ( 8) is outer semicontinuous, but not inner semicontinuous.The lack of inner semicontinuity prevents us from exploiting the "almost everywhere" conditions of (Della Rossa et al. (2021)) and references therein.Instead, one could resort to conditions involving Clarke's generalized directional derivative and Clarke's generalized gradient, which can be defined as (see (Clarke, 1990, page 11)) where Ω u is the set (of Lebesgue measure zero) where V is not differentiable, and is any other set of Lebesgue measure zero.However, due to the peculiar dynamics considered here (resembling, for example, the undesirable conservativeness highlighted in (Della Rossa, 2020, Ex. 2.2)), Lyapunov decrease conditions based on Clarke's generalized gradient would be too conservative and impossible to prove.Due to the above motivation, in this section we prove Propositions 2 and 3 by exploiting the results of (Valadier (1989); Bacciotti and Ceragioli (2003)), whose proof is also reported in (Della Rossa, 2020, Lemma 2.23), establishing a link between the time derivative d dt V (φ(t)) of a Lyapunov-like function V evaluated along a generic solution φ of a continuous-time system, and the so-called set-valued Lie derivative (Bacciotti and Ceragioli (1999)) with ∂ V defined in (29).In the following, we characterize some features of the set-valued Lie derivative, useful for the next technical derivations.
Lemma 8. Consider a set S ⊂ n , F : dom F ⊂ n ⇒ n with S ⊂ dom (F ) and a locally Lipschitz V : dom V ⊆ n → such that S ⊂ dom (V ).Given any function ϕ Proof: Consider any x ∈ S. In view of (30), and by the fact 〈ϕ(x, f ), f 〉 as to be proven.
Whenever function V is non-pathological (according to the definition given next), (Della Rossa, 2020, Lemma 2.23) ensures that d dt V (φ(t)) ∈ V F (φ(t)) for almost all t in the domain of φ.We provide below the definition of nonpathological functions and we prove that function V in (19) enjoys that property.Our result below can be seen as a corollary of the fact that piecewise C 1 functions (in the sense of (Scholtes (2012))) are non-pathological.This result has been recently published in (Della Rossa et al., 2022, Lemma 4).The scalar case is a consequence of Proposition 5 in (Valadier (1989)).An alternative proof is reported here.Definition 1. (Valadier (1989)) A locally Lipschitz function W : dom W ⊆ → is non-pathological if, given any absolutely continuous function φ : ≥0 → dom W , we have that for almost every t ∈ ≥0 there exists a t ∈ satisfying 〈w, φ(t)〉 = a t , ∀w ∈ ∂ W (φ(t)). (32) In other words, for almost every t ∈ ≥0 , ∂ W (φ(t)) is a subset of an affine subspace orthogonal to φ(t).
Proposition 4. Any function W : dom W ⊆ → that is piecewise continuously differentiable is non-pathological.
For this paper, the interest of Proposition 4 stands in the fact that it implies that function V in ( 19) is non-pathological (Valadier (1989); Bacciotti and Ceragioli (1999)), being the sum of piecewise continuously differentiable scalar functions.
Remark 3.An alternative proof of Proposition 4, might be to first establish that piecewise C 1 functions from → can be represented as a max-min of C 1 functions from → (a similar result has been proven, for example, in (Xu et al., 2017, Thm. 1 and Prop. 1) with reference to piecewise affine functions), and then obtain Proposition 4 as a corollary of (Della Rossa, 2020;Valadier, 1989;Bacciotti and Ceragioli, 1999, Lemma 2.20), which establish non-pathological properties of max-min functions.
Proof of Proposition 4: Let W : dom W ⊆ → be piecewise continuously differentiable.Let φ : ≥0 → dom W be an absolutely continuous scalar function and suppose it is differentiable at t ∈ ≥0 .We split the analysis in three cases.a) φ(t) = 0.Then, for all w ∈ ∂ W (φ(t)), 〈w, φ(t)〉 = 〈w, 0〉 = 0. Thus, ( 32) is satisfied with a t = 0. Since φ is absolutely continuous, then the set Φ where it is not differentiable is of Lebesgue measure zero.We conclude that ( 32) is satisfied for almost all t ∈ ≥0 because we have arbitrarily selected t ∈ ≥0 \ Φ, as to be proven.

Proof of Proposition 2
The following lemma establishes geometric properties of V that are used, together with Lemma 5 and Proposition 4, to show the trajectory-based results of Proposition 2. Lemma 9. Consider system (8), function V in (19)-( 20) and c as in Theorem 1.There exists α 3 ∈ ∞ , independent of ω in (14) and of κ, such that Proof: We prove the two equations one by one.
Proof of (33b).We split the analysis in two cases.First, let i ∈ , x ∈ D i and x + = g i (x) as in Lemma 2. In view of the equality in (6) and the definition of V i j in (20), we have Second, let (i, j) ∈ , x ∈ D i j and x + ∈ G i j (x) as in Lemma 1.In view of item a) of Property 1, On the other hand, in view of (4a), (4b) and Lemma 1, On the other hand, from the definition of G i j in (4b), V uv (x + ) = V uv (x) for any (u, v) = (i, j) ∈ .Therefore V (x + ) − V (x) = V i j (x + ) − V i j (x) < 0, since the arguments of all the other elements of the summation in (19) do not change.
Based on Lemma 9, we can now prove Proposition 2.

Proof of Proposition 3
Paralleling Section 8.2, we establish via the next lemma geometric properties of V that can be used, together with Lemma 5 and Proposition 4, to show the trajectory-based result of Proposition 3. Proof: We only prove (38a), as (38b) follows from the same arguments as those proving (33b) in Proposition 2. For each x ∈ C \ and each f ∈ F (x), we may proceed as in (34) and then exploit from Lemma 7 that, sup V F (x) ≤ −κλµ 2 + cω. (39) By selecting κ ≥ κ = 2cω λµ 2 , (39) yields sup V F (x) ≤ − 1 2 κλµ 2 , which proves (38a).
Based on Lemma 10, we are now ready to prove Proposition 3.

Conclusions
We presented a cyber-physical hybrid model of leaderless heterogeneous first-order oscillators, where global uniform asymptotic and/or finite-time synchronization is obtained in a distributed way via hybrid coupling.More specifically we establish that the synchronization set for the proposed model enjoys uniform asymptotic practical stability property.Thanks to the mild requirements on the coupling function, the stability result was strengthened to a prescribed finite-time property when the coupling function is discontinuous at the origin.Finally, we proved a useful statement on scalar non-pathological functions, exploited here for the nonsmooth Lyapunov analysis in our main theorems.We believe that this work demonstrates the potential of hybrid systems theoretical tools to overcome the fundamental limitations of continuous-time systems.
Future extensions of this work include addressing graphs with cycles (not trees) and investigating the case with leaders as done for a second-order Kuramoto model in (Bosso et al. (2021a)).Additional challenges may include studying the converse problem of globally de-synchronizing the network (Franci et al. (2012)) via hybrid approaches.
Acknowledgement.The authors would like to thank Matteo Della Rossa and Francesca Maria Ceragioli for useful discussions about the results of Section 8.1, Elena Panteley for useful discussions and suggestions given during the preparation of the manuscript, and the reviewers of the submission, which allowed us, through their comments, to improve the overall quality of the work.

Figure 1 :
Figure 1: Examples of functions σ satisfying Property 1, together with the sine function (which does not satisfy Property 1).

Figure 3 :
Figure3: The networks of oscillators described in Section 7: the blue graph depicts physical couplings captured by ω i in (27), and the red graph is the communication tree graph that we design for the hybrid coupling rules presented in Section 3.

Figure 5 :
Figure 5: Evolution of the phase errors for κ = 576π 10 and different selections of σ and communication configurations (from top to bottom: no communication layer, ramp and sign functions).