# Conditions for Global Synchronization in Networks of Heterogeneous Kuramoto Oscillators

José Ángel Mercado Uribe, Jesús Mendoza-Avila, Denis Efimov, Johannes Schiffer

► To cite this version:

José Ángel Mercado Uribe, Jesús Mendoza-Avila, Denis Efimov, Johannes Schiffer. Conditions for Global Synchronization in Networks of Heterogeneous Kuramoto Oscillators. Automatica, In press, 10.1016/j.automatica.2025.112180. hal-04864992

HAL Id: hal-04864992

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

Submitted on 5 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Conditions for Global Synchronization in Networks of Heterogeneous Kuramoto Oscillators

Angel Mercado-Uribe<sup>a</sup>, Jesus Mendoza-Avila<sup>b</sup>, Denis Efimov<sup>b</sup>, Johannes Schiffer<sup>a,c</sup>

## Abstract

The Kuramoto model is essential for studying synchronization. In this work, we present sufficient conditions for global synchronization in networks of heterogeneous Kuramoto oscillators in the absence of homoclinic and heteroclinic cycles. The result is established by constructing a suitable Leonov function candidate for the Kuramoto model, which provides sufficient conditions for almost global synchronization in networks with acyclic and meshed topologies. The synchronization property is accompanied by necessary and sufficient conditions to guarantee the existence of equilibria, which are satisfied if the conditions for synchronization hold. The implications of the main conditions and their relationship with the network topology and parameters are discussed. Finally, the results are illustrated via a numerical example.

**Key words:** Coupled oscillators, Kuramoto model, Synchronization, Complex networks, Nonlinear analysis.

## 1 Introduction

Synchronization phenomena can be found in many environments like nature, science, society, and technology [28, 29]. At a high level, this phenomenon can be seen as the coordination between multiple individual agents. Some classic examples are: pendulum clocks in the same wall, flashing fireflies, neurons, planets and synchronous generators [1, 28].

For this reason, the analysis of synchronization in networks of coupled oscillators has taken a crucial role over the years. Nowadays, it is one of the main pervasive topics in many disciplines, like neuroscience, biology, chemistry, physics, mathematics, engineering [1, 8], and more recently, computer science [17].

The behavior of a network of coupled oscillators is described by two main features [8, 14]: a population of oscillators that describe the coupling function and the interconnection topology, which is given by a graph. In the literature, many models of coupled oscillators can

be found. However, among all proposed models for describing the synchronization phenomena, the most successful, well-known and prominent oscillator representation is the Kuramoto model [1, 18]. This model considers a collection of coupled phase oscillators of first-order as population. Each agent is assumed to have an intrinsic natural frequency coupled by a weighted sinusoidal function of their phase differences.

Even though, at first glance, the coupled Kuramoto oscillator dynamics may seem to be relatively simple, they are complex enough due to their highly nontrivial and nonlinear behavior. In contrast to other models in the literature like the one proposed in [36], the Kuramoto model is sufficiently simple to work mathematically on it and still has a rich behavior to analyze the synchronization phenomenon. For this reason, the Kuramoto model has become a crucial benchmark for studying synchronization, and it can be found in many scientific disciplines [1, 8, 29].

The Kuramoto oscillators are usually studied in two configurations: they are said to be homogeneous if all natural frequencies are identical and heterogeneous if at least one of the natural frequencies is different. Likewise, the synchronization analysis can be associated with the phases or frequencies of each oscillator. The first case

<footnote><sup>a</sup>Brandenburg University of Technology Cottbus - Senftenberg, Germany.</footnote>

<footnote><sup>b</sup>Inria, Univ. Lille, CNRS UMR 9189 - CRISTAL, F-59000 Lille, France.</footnote>

<footnote><sup>c</sup>Fraunhofer Research Institution for Energy Infrastructures and Geothermal Systems (IEG), 03046 Cottbus, Germany</footnote>

<footnote>Email addresses: mercadou@b-tu.de (Angel Mercado-Uribe), jesus.mendoza-avila@inria.fr (Jesus Mendoza-Avila), denis.efimov@inria.fr (Denis Efimov), schiffer@b-tu.de (Johannes Schiffer).</footnote>
implies that all phases of the oscillators become iden-
tical over time, and the second one is obtained when
all frequencies tend to the same value [8]. One of the
primary analyses of coupled Kuramoto oscillators is fo-
cused on identifying conditions to reach global or local
phase and frequency synchronization. However, due to
the nonlinear periodic couplings between the oscillators,
the analysis requires investigations of a system with in-
finite equilibria and this is generally highly complicated.
Consequently, this usually forces us to perform the study
under some simplifications or considering only local syn-
chronization properties, see [8, 29].

For the case of the homogeneous oscillators, the global
phase and frequency synchronization can be obtained if a
locally exponentially stable phase synchronization man-
ifold exists. Likewise, this manifold is ensured for any
positive value of coupling strengths [16, 25]. In spite of
the analysis of homogeneous oscillators being relatively
simple, the heterogeneous case is by far more complex.
First, in this case, phase synchronization is not possi-
ble [8]. Second, the frequency synchronization analysis
implies a highly nonlinear system that depends on mul-
tiple factors, like the number of oscillators, maximum
difference of natural frequencies, and steady-state phase
difference. For this reason, most results in the literature
consider the local synchronization problem for the het-
erogeneous case, which includes the estimation of the re-
gion of attraction, see [7, 35, 37] and the surveys [8, 29].
Likewise, for global results, [14] considers graphs with di-
ameter less than or equal to 2, i.e., the maximum short-
est path from one node to another is less than 3, and
provide sufficient conditions for synchronization. After,
these results are extended for arbitrary topology in [13].
However, the last result is only locally obtained.

Most of the works in the literature for the analysis of
synchronization of the Kuramoto model are based on
Lyapunov theory. However, the periodic non-linearities
associated with the Kuramoto model coupling function
generate multiple equilibrium points. This complicates
the global analysis, and the synchronization results are
usually local, especially for the heterogeneous case. For
obtaining global results, the main problem is finding con-
ditions to guarantee the boundedness of solutions, see
[14]. Recently, the cell structure approach has been de-
veloped in the literature. This theory allows us to find
conditions for the global boundedness of trajectories in
periodic systems. This approach relaxes the conditions
of the classical Lyapunov theory by considering the pe-
riodicity of the system. The first result of this theory
was developed for systems with only one periodic state
in [19, 27]; this approach is known as the Leonov func-
tion theory. Some years later, this innovative result was
extended to multiple periodic variables in [10, 30]. Using
this approach, some robustness properties can be veri-
fied for state periodic systems as, for example, the con-
ditions of Input-to-State Stability (ISS) in [11] and in-
tegral ISS (iISS) in [23]. Additionally, control synthesis

of state periodic systems has been proposed in [21] for
boundedness of solutions through the concept of control
Leonov function, and for ISS- (iISS-) estabilization in
[22] through the concept of ISS- (iISS-)control Leonov
function.

The present work is focused on finding sufficient condi-
tions for global synchronization in networks of hetero-
geneous Kuramoto oscillators to complement our previ-
ous work [24]. As a consequence of the multiple equilib-
ria in the Kuramoto system, to manage the objective,
the recently developed Leonov function framework [30]
is employed instead of standard Lyapunov techniques.
In summary, the main contributions are:

* Characterization of the equilibrium sets and, for the acyclic case, also their local stability properties.

* Sufficient conditions for *global* boundedness of the Kuramoto system trajectories for the case of meshed and acyclic network topologies are presented using a suitable Leonov function.

* Building upon the results mentioned above, global synchronization of all oscillators is proven by using a LaSalle argument.

* Based on two corollaries to our main result, we show how the derived conditions can be associated with the main features of the network.

* The results obtained in this work are illustrated and discussed via a numerical example based on [14, Example 1].

Hence, the analysis in the present paper extends our
previous work [24] in several respects: first, we consider
meshed in addition to acyclic network topologies in the
analysis; second, besides necessary conditions on the ex-
istence of equilibria, we now also provide a characteri-
zation of their local stability properties for the acyclic
case; third, we provide a more detailed proof of our re-
sults; fourth, we extensively discuss the (physical) im-
plications of the proposed synchronization conditions as
well as the system behavior in the presence of hetero-
clinic cycles. Finally, we illustrate our findings via sev-
eral numerical examples.

The remainder of the paper is structured as follows. Some essential preliminaries are presented in Section 2. The detailed model is given in Section 3. The main result is provided in Section 4. A physical interpretation of the derived conditions is discussed in Section 5. A simulation example is shown in Section 6. Finally, conclusions and future results are discussed in Section 7.

## 2 Preliminaries

### 2.1 Notation

We denote the set of real and integer numbers by $\mathbb{R}$ and $\mathbb{Z}$, respectively. Let $\mathbb{R}_+ = \{s \in \mathbb{R} | s \ge 0\}$ and
$\mathbb{R}_{++} = \mathbb{R}_{+} \setminus \{0\}$ as well as $\mathbb{Z}_{+} = \mathbb{Z} \cap \mathbb{R}_{+}$ and $\mathbb{Z}_{++} = \mathbb{Z} \cap \mathbb{R}_{++}$. A real (integer) vector is denoted by $v \in \mathbb{R}^k$ ($v \in \mathbb{Z}^k$) with $k \in \mathbb{Z}_{++}$ components, its infinity norm is represented by $\|v\|_{\infty} = \max(|v_1|, ..., |v_k|)$. The identity matrix is denoted by $\mathbf{I}_k \in \mathbb{R}^{k \times k}$. Let $M \in \mathbb{R}^{n \times m}$ be a matrix. Then, its induced infinity norm is defined as $\|M\|_{\infty} = \max_{\{x \in \mathbb{R} \mid \|x\|_{\infty}=1\}} (\|Mx\|_{\infty})$. Likewise, the rank, range and kernel of the matrix $M$ are denoted by rank($M$), ran($M$) and ker($M$), respectively. For a square matrix $A \in \mathbb{R}^{n \times n}$, $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ denote its smallest and largest eigenvalues. Let $\mathbf{1}_k$ and $\mathbf{0}_k$ be the $k$-dimensional vectors of ones and zeros, respectively. For $x \in \mathbb{R}^k$ define $\sin(x) = [\sin(x_1), ..., \sin(x_k)]^\top$, $\cos(x) = [\cos(x_1), ..., \cos(x_k)]^\top$ and diag($x$) as a $k \times k$ diagonal matrix with diagonal entries $x_i$, $i = 1, ..., k$. The set $\mathbb{S}$ denotes the unit circle, an angle is a point $\theta \in \mathbb{S}$ and $v_\theta \in \mathbb{S}^k$ represents a vector in the $k^{th}$ sphere.

## 2.2 An auxiliary lemma

To analyze local stability properties of equilibria of the Kuramoto dynamics, in the sequel we make use of the following lemma.

**Lemma 1** Let $M \in \mathbb{R}^{n \times n}$ and $D \in \mathbb{R}^{n \times n}$ be a symmetric positive definite matrix and a diagonal matrix, respectively. Then, the matrices $D$ and $MD$ (or $DM$) have the same number of positive, negative and zero eigenvalues.

*Proof.* The fact that the principal square root of a positive definite matrix is real and positive definite [15, Theorem 7.2.6], implies that $\sqrt{M}$ may be chosen real. Next, by taking $\sqrt{M}$ as basis matrix, we see that the matrix products $MD$ and $DM$ are both similar to $\sqrt{MD}\sqrt{M}$ and, hence, both matrices have the same eigenvalues [5, Proposition 4.4.4]. Since $\sqrt{M}$ is real positive definite, the matrices $D$ and $\sqrt{MD}\sqrt{M}$ are congruent and by Sylvester's law of inertia [6, Theorem 9.13], the matrix $\sqrt{MD}\sqrt{M}$ has the same number of positive, negative and zero eigenvalues as the matrix $D$. □

## 2.3 Global boundedness of solutions of state periodic systems via Leonov functions

We recall some results on the global boundedness of solutions of state periodic systems via the notion of Leonov functions [10, 30]. Consider the system

$$ \dot{x} = f(x(t)), \quad \forall t \ge 0, \qquad (1) $$

where $x(t) \in \mathbb{R}^n$ are the states and $f(x) : \mathbb{R}^n \to \mathbb{R}^n$ is a $C^1$ map with $f(\mathbf{0}_n) = \mathbf{0}_n$. We assume that (1) admits uniquely defined solutions and denote the solution of (1) at time $t$ by $X(t, x)$, which fulfills $X(0, x) = x \in \mathbb{R}^n$.

In this work, we focus on systems, whose dynamics are periodic with respect to all state variables, which is formalized in the assumption below.

**Assumption 2** The vector field $f$ in (1) is $2\pi$-periodic, i.e., $f(x + 2\pi j) = f(x)$, where $j \in \mathbb{Z}^n$ is a multi-index vector, for all $x \in \mathbb{R}^n$.

Thus, the system (1) can be embedded into a manifold $\mathbb{S}^n$ by a simple projection of the variables $x$ on the set $\mathbb{S}^n$, where $\mathbb{S} = [0, 2\pi)$. By Assumption 2, the dynamics (1) have several equilibria. For the subsequent exposition, we introduce two auxiliary sets. Let $\mathcal{W} \subset \mathbb{R}^n$ be a compact set and $\mathcal{U} \subset \mathbb{R}^n$, such that for some $\underline{c} \in (0, 2\pi)$ and $\bar{c} \in (\pi, 2\pi)$ with $\underline{c} < \bar{c}$,

$$ \begin{gathered} \{x \in \mathbb{R}^n | \underline{c} \le \|x\|_{\infty} \le \bar{c}\} \subseteq \mathcal{W}, \\ \mathcal{U} = \bigcup_{r \in \mathbb{Z}_+} \mathcal{U}_r, \qquad (2) \\ \mathcal{U}_r = \{x \in \mathbb{R}^n | \|x\|_{\infty} = 2r\pi, f(x) = \mathbf{0}_n\}. \end{gathered} $$

Then, the set $\mathcal{U}$ includes all equilibria of the system (1), which are obtained by shifting the one at the origin. In general, there are also other equilibria of (1) that are not included in $\mathcal{U}$.

Now, we can introduce an appropriate definition of a strong Leonov function for the system (1), under Assumption 2, as follows, see [30]:

**Definition 3** A $C^1$ function $V : \mathbb{R}^n \to \mathbb{R}$ is a strong Leonov function for the system (1) if

$$ \inf_{x \in \mathcal{W}} V(x) > 0, \quad \sup_{x \in \mathcal{U}} V(x) \le 0, \qquad (3) $$

and there exists a continuous function $\lambda : \mathbb{R} \to \mathbb{R}$, satisfying $\lambda(0) = 0$ and $\lambda(s)s > 0$ for all $s \ne 0$, such that

$$ \frac{\partial V(x)}{\partial x} f(x) \le -\lambda(V(x)) \quad \forall x \in \mathbb{R}^n. \qquad (4) $$

In order to provide the preliminary theory, we define the sets

$$ \begin{gathered} \Omega = \{x \in \mathbb{R}^n | V(x) \le 0\}, \\ \mathcal{D}_s \supseteq \mathcal{D} = \{x \in \mathbb{R}^n | f(x) = \mathbf{0}_n\}, \end{gathered} \qquad (5) $$

where $\mathcal{D}_s$ contains all invariant sets for (1), which may also comprise unbounded continua of equilibria. For the subsequent analysis, we need the assumption below, whose formulation follows [2] and considers the solutions of the system (1) on the manifold $\mathbb{S}^n$.

**Assumption 4** The set $\mathcal{D}_s \cap \mathbb{S}^n$ admits a finite decomposition without cycles $\mathcal{D}_s \cap \mathbb{S}^n = \bigcup_{i=1}^k \mathcal{D}_{s_i}$ for some non-empty disjoint compact sets $\mathcal{D}_{s_i}$, which form a filtration ordering of $\mathcal{D}_s$, as detailed in [2, Definitions 3 and 4].
In general, systems with multiple invariant sets may exhibit solutions that travel between different equilibria, e.g., homoclinic or heteroclinic cycles. The presence of such cycles affects any stability criteria, which are based on Lyapunov-like functions [2]. The absence of such particular solutions is extremely hard to verify analytically and this can only be done numerically in most of the cases [4]. Therefore, Assumption 4 is standard in most works on multistable systems, see [2] and, for completeness, also needs to be added to the result from [30].

In the Kuramoto system, for example, the presence of homoclinic and heteroclinic cycles would imply that there are trajectories connecting distinct equilibrium points or an equilibrium point with itself and, as a consequence, the oscillators' frequencies cannot synchronize as it is shown in Section 6.

By using Definition 3, the following theorem is stated.

**Theorem 5** [30] If for the system (1) under Assumptions 2 and 4, there exists a strong Leonov function, then for all $x_0 \in \mathbb{R}^n$ the corresponding trajectories $X(t, x_0)$ are bounded for all $t \ge 0$.

Furthermore, the following corollary can be stated, which is instrumental for performing the global synchronization analysis of networks of heterogeneous Kuramoto oscillators in the present paper.

**Corollary 6** [30] If for the system (1) under Assumptions 2 and 4, there exists a Leonov function with

$$V(0) \le 0, \inf_{x \in U \setminus U_0} V(x) \le 0, \inf_{x \in W} V(x) > 0, \quad (6)$$

and satisfying the following dissipation inequality

$$\frac{\partial V(x)}{\partial x} f(x) < 0 \quad \forall x \in \mathbb{R}^n \setminus \{\mathcal{D} \cup \Omega\}, \quad (7)$$

then for all $x_0 \in \mathbb{R}^n$, the corresponding trajectories $X(t, x_0)$ are bounded for all $t \ge 0$.

# 3 Modeling networks of first-order Kuramoto oscillators

## 3.1 Model of coupled Kuramoto oscillators

We consider a network of first-order heterogeneous Kuramoto oscillators with $N > 1$ states, see [8]. The network topology is described by an undirected and connected graph $G = (\mathcal{N}, \mathcal{E})$, where the set of network nodes is denoted by $\mathcal{N} = \{1, 2, \dots, N\}$ and the set of edges by $\mathcal{E} = \{e_1, \dots, e_m\}$. By associating an arbitrary ordering to the edges, we introduce the node-edge incidence matrix $\mathcal{B} \in \mathbb{R}^{N \times m}$, whose entries are defined as $b_{i\ell} = 1$ if

node $i$ is the source of the $\ell$-th edge $e_\ell$, $b_{i\ell} = -1$ if $i$ is the sink of $e_\ell$ and $b_{i\ell} = 0$ otherwise. Thus, $\mathbf{1}_N^\top \mathcal{B} = \mathbf{0}_m^\top$.

The dynamics of the oscillator at node $i \in \mathcal{N}$ is given by

$$\dot{\theta}_i = \omega_i - \sum_{k \in \mathcal{N}_i} a_{ik} \sin(\theta_i - \theta_j), \quad (8)$$

where $\theta_i : \mathbb{R}_+ \to \mathbb{R}$ is the phase of the $i$-th oscillator, $\omega_i \in \mathbb{R}$ is its natural frequency and $a_{ik} = a_{ki} = a_\ell \in \mathbb{R}_{>0}$ is the coupling strength of the $\ell$-th edge connecting nodes $i$ and $k$. The set of neighbors of a node $i \in \mathcal{N}$ is denoted by $\mathcal{N}_i := \{k \mid k \in \mathcal{N}, k \neq i, a_{ik} \neq 0\}$.

To write the network dynamics (8) compactly, we define the matrix of coupling strengths as follows

$$A = \operatorname{diag} \left( \left[ a_1, \dots, a_m \right]^\top \right) > 0, \quad A \in \mathbb{R}^{m \times m}. \quad (9)$$

Also, we introduce the potential function $U_\theta : \mathbb{R}^N \to \mathbb{R}$,

$$U_\theta(\theta) = -\mathbf{1}_m^\top A \cos(\mathcal{B}^\top \theta), \quad (10)$$

and its gradient,

$$\nabla U_\theta(\theta) = \mathcal{B} A \sin(\mathcal{B}^\top \theta), \quad (11)$$

possessing the following symmetry properties for any $\alpha \in \mathbb{R}$,

$$U_\theta(\theta + \alpha \mathbf{1}_N) = U_\theta(\theta), \quad \nabla U_\theta(\theta + \alpha \mathbf{1}_N) = \nabla U_\theta(\theta). \quad (12)$$

These symmetry properties arise from the fact that the coupling terms in (8) only depend upon angle differences. Furthermore, since $\mathbf{1}_N^\top \mathcal{B} = \mathbf{0}_m^\top$,

$$\mathbf{1}_N^\top \nabla U_\theta(\theta) = \mathbf{0}_N^\top. \quad (13)$$

Then, the model (8) can be compactly written in gradient-form, i.e.,

$$\dot{\theta} = \omega - \nabla U_\theta(\theta). \quad (14)$$

Similarly to [31, 32, 8], we introduce the notions below to characterize desired invariant solutions of the system (14).

**Definition 7** The system (14) is frequency synchronized if it has an invariant solution for all $t \ge 0$ of the form

$$\theta^s(t) = \omega^s t + \theta_0^s, \quad \omega^s = \mathbf{1}_N \omega^* = \omega - \nabla U_\theta(\theta^s(t)), \quad (15)$$

where $\omega^* \in \mathbb{R}$ and $\theta_0^s \in \mathbb{R}^N$. The system (14) admits a synchronized motion if it is frequency synchronized with

$$|\theta_{0,i}^s - \theta_{0,k}^s| < \frac{\pi}{2} \quad \forall i \in \mathcal{N}, \quad \forall k \in \mathcal{N}_i.$$
From (13) and (15), we obtain $\mathbf{1}_N^\top \mathbf{1}_N \omega^* = \mathbf{1}_N^\top \omega$. Thus, it follows that the synchronization frequency $\omega^*$ is uniquely determined by the average of the natural frequencies $\omega$ [9], i.e.,

$$ \omega^* = \frac{\mathbf{1}_N^\top \omega}{\mathbf{1}_N^\top \mathbf{1}_N} = \frac{1}{N} \sum_{i=1}^{N} \omega_i. \quad (16) $$

## 3.2 Change of coordinates

Inspired by [3, 31, 33], we perform the following change of coordinates

$$ \eta = \mathcal{B}^\top \theta \in \mathbb{R}^m. \quad (17) $$

It follows from (10) that

$$ U_{\eta}(\eta) = -\mathbf{1}_{m}^{\top} A \cos(\eta), \quad \nabla U_{\eta}(\eta) = A \sin(\eta), \quad (18) $$

and thus from (11), we obtain

$$ \nabla U_{\theta}(\theta) = \mathcal{B} \nabla U_{\eta}(\eta). \quad (19) $$

Hence, in the new coordinates $\eta$, the dynamics (14) become

$$ \dot{\eta} = \mathcal{B}^{\top} \omega - \mathcal{B}^{\top} \mathcal{B} \nabla U_{\eta}(\eta) = \mathcal{B}^{\top} \omega - \mathcal{B}^{\top} \mathcal{B} A \sin(\eta), \quad (20) $$

where $\mathcal{B}^\top \omega$ is the differences in natural frequencies between any pair of coupled oscillators.

Clearly, any synchronized solution of the system (14) corresponds to an equilibrium of (20), i.e.,

$$ \eta^s = \mathcal{B}^{\top} (\theta_0^s + \omega^s t + \alpha \mathbf{1}_N) = \mathcal{B}^{\top} (\theta_0^s + (\omega^s t + \alpha) \mathbf{1}_N) = \mathcal{B}^{\top} \theta_0^s. $$

Thus, we focus our analysis in the sequel on the dynamics (20). Thereby, the main objective of this work is to identify suitable coupling strengths $A$ and natural frequency differences $\bar{\omega} = \mathcal{B}^\top \omega$, such that the system (20) possesses a globally attractive equilibrium set, which implies that the trajectories of the system (8) synchronize.

## 3.3 Equilibrium points

The equilibrium set of the dynamics (20) is given by

$$ \mathcal{D} = \{\eta \in \mathbb{R}^m \mid \mathcal{B}^\top \omega = \mathcal{B}^\top \mathcal{B} A \sin(\eta)\}. \quad (21) $$

Since the magnitude of the sinusoidal function is bounded by one, solutions to (21) may not always exist. Therefore, we analyze the non-emptiness of the equilibrium set $\mathcal{D}$ in (21). This set being non-empty implies that any equilibrium point $\eta^*$ of the system (20), has to satisfy the following equation

$$ \mathcal{B}^{\top} [\omega - \mathcal{B} A \sin(\eta^*)] = \mathbf{0}_m. \quad (22) $$

To solve the previous equation, we introduce the pseudoinverse (or generalized inverse) of $\mathcal{B}$ as $\mathcal{B}^\dagger$, which satisfies $\mathcal{B}\mathcal{B}^\dagger\mathcal{B} = \mathcal{B}$, $\mathcal{B}^\dagger\mathcal{B}\mathcal{B}^\dagger = \mathcal{B}^\dagger$, $(\mathcal{B}\mathcal{B}^\dagger)^\top = \mathcal{B}\mathcal{B}^\dagger$ and $(\mathcal{B}^\dagger\mathcal{B})^\top = \mathcal{B}^\dagger\mathcal{B}$. Using the pseudoinverse $\mathcal{B}^\dagger$, the following proposition can be stated.

**Proposition 8** Consider the system (20). The set $\mathcal{D}$ in (21) is non-empty if and only if there exists a vector $v_{\mathcal{B}_0} \in \ker(\mathcal{B})$, such that

$$ \|A^{-1}(\mathcal{B}^\dagger \omega + v_{\mathcal{B}_0})\|_{\infty} \le 1. \quad (23) $$

Furthermore, if and only if (23) can be satisfied with strict inequality, then the equilibria set $\mathcal{D}$ is formed by $\eta_i^*$ (mod $2\pi$), where $\eta_i^*$, $i = 1, \dots, m$, are permutations of the two vectors

$$ \eta^* = \arcsin(A^{-1}(\mathcal{B}^\dagger \omega + v_{\mathcal{B}_0})), \quad \hat{\eta}^* = \pi \mathbf{1}_m - \eta^*, \quad (24) $$

for all $v_{\mathcal{B}_0} \in \ker(\mathcal{B})$ that satisfy (23) with strict inequality.

*Proof.* Since $\mathcal{B}$ is the incidence matrix of a connected graph, we have that $\ker(\mathcal{B}) = \ker(\mathcal{B}^\top) = N - 1$ with $\ker(\mathcal{B}^\top) = \operatorname{span}\{\mathbf{0}_N, \mathbf{1}_N\}$ [12]. Hence, from equation (22), the equilibria of the system (20) can be computed as

$$ \omega - \mathcal{B} A \sin(\eta^*) = v_0, \quad v_0 \in \operatorname{span}\{\mathbf{0}_N, \mathbf{1}_N\}. \quad (25) $$

Multiplying by the matrix $\mathcal{B}\mathcal{B}^\dagger$ from the left, we obtain

$$ \mathcal{B}\mathcal{B}^\dagger\omega - \mathcal{B}\mathcal{B}^\dagger\mathcal{B}A\sin(\eta^*) = \mathcal{B}\mathcal{B}^\dagger v_0. \quad (26) $$

Since $\ker(\mathcal{B}^\top) = \ker(\mathcal{B}^\top)$ as well as $\operatorname{ran}(\mathcal{B}^\top) = \operatorname{ran}(\mathcal{B}^\top)$ and using the property $(\mathcal{B}\mathcal{B}^\top)^\top = \mathcal{B}\mathcal{B}^\top$ together with the fact that $\operatorname{ran}(\mathcal{B}^\top) \cap \ker(\mathcal{B}) = \mathbf{0}_m$, we conclude that $\ker(\mathcal{B}\mathcal{B}^\top) = \ker(\mathcal{B}^\top)$. Therefore, although $\mathcal{B}\mathcal{B}^\top$ is not a full-rank matrix, all possible solutions of (22) are equivalent to those of (26). Furthermore, due to the fact that $v_0 \in \ker(\mathcal{B}^\top) = \operatorname{span}\{\mathbf{0}_N, \mathbf{1}_N\}$, (26) is equivalent to

$$ \begin{aligned} \mathbf{0}_m &= \mathcal{B}\mathcal{B}^\dagger\omega - \mathcal{B}\mathcal{B}^\dagger\mathcal{B}A\sin(\eta^*) \\ &= \mathcal{B}[\mathcal{B}^\dagger\omega - A\sin(\eta^*)] \\ &= \mathcal{B}[\mathcal{B}^\dagger\omega - A\sin(\eta^*) + v_{\mathcal{B}_0} - v_{\mathcal{B}_0}], \end{aligned} $$

where we have used the identity $\mathcal{B}\mathcal{B}^\dagger\mathcal{B} = \mathcal{B}$ and $v_{\mathcal{B}_0} \in \ker(\mathcal{B})$. Due to $\mathcal{B}v_{\mathcal{B}_0} = \mathbf{0}_m$, we obtain

$$ \mathcal{B}[\mathcal{B}^\dagger\omega - A\sin(\eta^*) + v_{\mathcal{B}_0}] = \mathbf{0}_m. \quad (27) $$

Solving the previous equation, we find that the solutions for $\eta^*$ are given by

$$ \sin(\eta^*) = A^{-1}(\mathcal{B}^\dagger \omega + v_{\mathcal{B}_0}), $$
and these correspond to (24).

From (23) and (24), the system (20) may exhibit a con-
tinuum of equilibria, since inequality (23) may, in gen-
eral, be satisfied for a continuous range of $v_{B_0} \in \ker(\mathcal{B})$.
A similar result is reported in [9], where the pseudoin-
verse of the Laplacian matrix is used instead of $\mathcal{B}^\dagger$.

**Remark 9** The dimension of $\eta$ in (17) depends on the range of $\mathcal{B}$. Therefore, isolated equilibria of the system (20) can only be obtained for the acyclic case. In order to remove this degree of freedom in the meshed case, a commonly used approach is to use one phase angle $\theta_i$ of the system (14) as a reference [20, 32]. However, in the present setting, it seems technically more convenient to work with the coordinates $\eta$ in (17), since this simplifies the construction of a Leonov function that satisfies the conditions of Corollary 6.

Based on Proposition 8, we consider the following natu-
ral assumption.

**Assumption 10** The condition (23) is satisfied with strict inequality for some $v_{B_0} \in \ker(\mathcal{B})$.

Assumption 10 ensures the existence of equilibria for the
system (20), meaning that the set $\mathcal{D}$ in (21) is non-empty.

3.4 Error dynamics

With Assumptions 10 and defining the variable $\tilde{\eta} = \eta - \eta^*$, where $\eta^*$ represents an equilibrium point in $\mathcal{D}$, the error dynamics of (20) results in

$$
\dot{\tilde{\eta}} = - \mathcal{B}^{\top} \mathcal{B} \zeta(\tilde{\eta}), \quad (28)
$$

with

$$
\begin{align}
\zeta(\tilde{\eta}) &= \nabla U_{\eta}(\tilde{\eta} + \eta^*) - \nabla U_{\eta}(\eta^*) \nonumber \\
&= A [\mathbf{sin}(\tilde{\eta} + \eta^*) - \mathbf{sin}(\eta^*)]. \tag{29}
\end{align}
$$

In terms of $\tilde{\eta}$, the equilibria set $\mathcal{D}$ can be written as

$$
\mathcal{D} = \{ \tilde{\eta} \in \mathbb{R}^m \mid \mathcal{B}^{\top} \mathcal{B} \zeta(\tilde{\eta}) = \mathbf{0}_m \}. \quad (30)
$$

4 Sufficient Conditions for Global Synchroniza-
tion

In this section, we derive sufficient conditions for global
attractiveness of the set $\mathcal{D}$ in (30), which implies that
asymptotically all phase differences remain constant and
hence, their frequencies are synchronized. For the case
of acyclic network topologies, we additionally present
conditions under which for all initial conditions, except
a set of measure zero, the solutions of the system (28)
converge to a stable equilibrium point.

To establish the result, we invoke Corollary 6, which
is instrumental in guaranteeing global boundedness of
solutions. To this end, we propose the following Leonov
function candidate, with a design parameter $\alpha \in [0, 1)$,

$$
\begin{equation}
\begin{split}
V(\tilde{\eta}, \eta^*) ={}& U_{\eta}(\tilde{\eta} + \eta^*) - U_{\eta}(\eta^*) - \nabla U_{\eta}^{\top}(\eta^*) \tilde{\eta} + \\
& \frac{\alpha}{2} \zeta^{\top}(\tilde{\eta}) A^{-1} \zeta(\tilde{\eta}).
\end{split}
\tag{31}
\end{equation}
$$

Using (18) and the fact that $A$ is diagonal, the function $V(\tilde{\eta}, \eta^*)$ in (31) can be expanded as

$$
\begin{align*}
V(\tilde{\eta}, \eta^*) &= \mathbb{1}_m^\top A (\cos(\eta^*) - \cos(\eta)) - (A \sin(\eta^*))^\top \tilde{\eta} + \\
&\quad \frac{\alpha}{2} \zeta^\top(\tilde{\eta}) A^{-1} \zeta(\tilde{\eta}) \\
&= \mathbb{1}_m^\top A [\cos(\eta^*) - \cos(\eta)] - (A \sin(\eta^*))^\top \tilde{\eta} + \\
&\quad \frac{\alpha}{2} \mathbb{1}_m^\top A (\sin(\tilde{\eta} + \eta^*) - \sin(\eta^*))^2,
\end{align*}
$$

with the power $(\cdot)^2$ applied element-wise. The previous
equality can be rewritten as follows

$$
V(\tilde{\eta}, \eta^*) = \sum_{j=1}^{m} a_j \mathcal{V}(\tilde{\eta}_j, \eta_j^*),
$$

where

$$
\begin{equation}
\begin{split}
\mathcal{V}(\tilde{\eta}_j, \eta_j^*) = & - \cos(\tilde{\eta}_j + \eta_j^*) + \cos(\eta_j^*) - \tilde{\eta}_j \sin(\eta_j^*) + \\
& \frac{\alpha}{2} (\sin(\tilde{\eta}_j + \eta_j^*) - \sin(\eta_j^*))^2.
\end{split}
\tag{32}
\end{equation}
$$

4.1 Preliminary lemmata

A crucial point in the application of the Leonov func-
tion framework, see Section 2.3, is to establish that the
Leonov function candidate $V$ in (31) is strictly positive
on the set $\mathcal{W}$ defined in (2). In order to ensure this prop-
erty for the dynamics (28), two preliminary lemmata are
stated.

**Lemma 11** Let $x \in \mathbb{R}$ be a variable and $x^* \in \mathbb{R}$ be a constant. The function $\mathcal{V}(x, x^*)$ in (32) has critical points in the following two sets

$$
\begin{align}
S_1 &= \{x \in \mathbb{R} \mid x = 2r\pi \quad \forall r \in \mathbb{Z}\}, \nonumber \\
S_2 &= \{x \in \mathbb{R} \mid x = \pi - 2x^* + 2r\pi \quad \forall r \in \mathbb{Z}\}. \tag{33}
\end{align}
$$

Moreover, if $|x^*| < \frac{\pi}{2}$, the critical points in $S_1$ are min-
ima and those in $S_2$ are maxima.

*Proof.* In order to calculate the maximum and minimum points of $\mathcal{V}(x, x^*)$ in (32), we use its partial derivative with respect to $x$, i.e.,

$$
\frac{\partial}{\partial x} \mathcal{V}(x, x^*) = (1 + \alpha \cos(x + x^*)) (\sin(x + x^*) - \sin(x^*)).
$$
The first component of the product on the right-hand side cannot be zero since $\alpha \in [0, 1)$. Therefore, the critical points of $V_F$ are obtained, when the second term is zero, which is zero on the set $S = S_1 \cup S_2$. To determine the extremal properties of the critical points, we compute the second partial derivative of $V$, i.e.,

$$ \frac{\partial^2}{\partial x^2} V(x, x^*) = (1 + \alpha \cos(x + x^*)) \cos(x + x^*) - \alpha \sin(x + x^*) (\sin(x + x^*) - \sin(x^*)). $$

By evaluating it on the set $S$, we have

$$ \left. \frac{\partial^2}{\partial x^2} V(x, x^*) \right|_{x \in S} = \left( (1 + \alpha \cos(x + x^*)) \cos(x + x^*) \right) \bigg|_{x \in S}. $$

Since $\alpha \in [0, 1)$, we see that if $|x^*| < \frac{\pi}{2}$, then

$$ \left. \frac{\partial^2}{\partial x^2} V(x) \right|_{x \in S_1} > 0, \quad \left. \frac{\partial^2}{\partial x^2} V(x) \right|_{x \in S_2} < 0, $$

which implies that the points in $S_1$ are minima and the points in $S_2$ are maxima. $\square$

Lemma 11 shows that if $\|\eta^*\|_{\infty} < \frac{\pi}{2}$, the Leonov function $V$ in (31) is rising near to the origin and has a maximum before $\|\eta\|_{\infty} = \pi$ for any nonzero $\eta^* \in \mathbb{R}^m$. For the existence of $\bar{c} \in (\pi, 2\pi)$, see (2), in the present case, it is thus necessary that $V(\tilde{\eta}, \eta^*) > 0$ for all $\tilde{\eta} \in \{\tilde{\eta} \in \mathbb{R}^m \mid \|\tilde{\eta}\|_{\infty} = \pi\}$. Hence, to determine the maximum value, fix $x = \pm\pi$, let $x^* \in \mathbb{R}$ and consider the function

$$ \begin{aligned} \mathcal{F}_V(x^*) &= \mathcal{V}(\pm\pi, x^*) \\ &= 2\cos(x^*) + 2\alpha \sin^2(x^*) \mp \pi \sin(x^*) \end{aligned} \quad (34) $$

along with the following lemma.

**Lemma 12** Let $\alpha \in [0, 1)$ be a constant and let $x^* \in \mathbb{R}$. Then, there exists $\bar{x}^*(\alpha) \ge 0.5669$, such that $\mathcal{F}_V(x^*) > 0$ for all $x^* \in (-\bar{x}^*(\alpha), \bar{x}^*(\alpha))$.

*Proof.* At first, consider the reduced function

$$ \mathcal{F}_V(x^*) = 2 \cos(x^*) \mp \pi \sin(x^*), $$

which satisfies $\mathcal{F}_V(0) = 2 > 0$. Therefore, this function is positive until the nearest zero crossing point $\bar{x}^*$. To determine this point, we consider

$$ \mathcal{F}_V(\bar{x}^*) = 2 \cos(\bar{x}^*) \mp \pi \sin(\bar{x}^*) = 0, $$

which is satisfied for $\tan(\bar{x}^*) = \pm \frac{2}{\pi}$ and, hence, $\bar{x}^* = \pm \arctan(\frac{2}{\pi}) = \pm 0.5669$. Thus, $\mathcal{F}_V$ takes only positive values for $x^* \in (-0.5669, 0.5669)$. Note that the term, which is associated to $\alpha$ in (34) is always non-negative, i.e., $2\alpha \sin^2(x^*) \ge 0$ for all $x^* \in \mathbb{R}$. Consequently, the limits of $x^*$ such that $\mathcal{F}_V(x^*) > 0$ for $x^* \ge \bar{x}^*$ can be increased by enlarging $\alpha$. $\square$

Since the function $\mathcal{F}_V$ in (34) is highly nonlinear, there is no compact analytic expression for $\bar{x}^*(\alpha)$. However, using a numerical solver, we can obtain the graph in Figure 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1: $\mathcal{F}_V(x^*(\alpha)) = 0$.</figcaption></figure>

## 4.2 Global boundedness of trajectories

By using the Leonov function framework, the following theorem is formulated.

**Proposition 13** Consider the system (28) with Assumptions 4 and 10. Let $\alpha \in [0, 1)$ be such that

$$ x^{\top} (\mathcal{B}^{\top} \mathcal{B} + \frac{\alpha}{2} \mathcal{M}(\phi)) x > 0, \quad (35) $$

for all $\phi \in \mathbb{S}^m$ and $x \in \mathbb{R}^m \setminus \mathcal{D}$, with $\mathcal{D}$ given in (21) and

$$ \mathcal{M}(\phi) = [\operatorname{diag}(\cos(\phi)) \mathcal{B}^{\top} \mathcal{B} + \mathcal{B}^{\top} \mathcal{B} \operatorname{diag}(\cos(\phi))] $$

Then, all solutions of the system (28) are bounded if there exists $\eta^*$ in (24) that satisfies

$$ \|\eta^*\|_{\infty} < \bar{x}^*(\alpha). \quad (36) $$

*Proof.* The claim is established by using Corollary 6 with the candidate Leonov function $V$ in (31). We begin by evaluating the requirements in (6). Clearly, $V(0) = 0$ and, since $V$ contains a linear term in $\tilde{\eta}$ and both $U_\eta(\tilde{\eta} + \eta^*)$ as well as $\zeta(\tilde{\eta} + \eta^*)$ are bounded functions in $\tilde{\eta}$, the first two conditions in (6) are satisfied. Lemmas 11 and 12 ensure that if (36) holds, the last condition in (6) is also verified. Consequently, all requirements in (6) are satisfied.

The derivative of $V$ along solutions of the system (28) is

$$ \dot{V} = [\nabla U_{\eta}^{\top}(\tilde{\eta} + \eta^*) - \nabla U_{\eta}^{\top}(\eta^*)] \dot{\tilde{\eta}} + \frac{\alpha}{2} [\zeta^{\top}(\tilde{\eta}) A^{-1} \dot{\zeta}(\tilde{\eta}) + \dot{\zeta}^{\top}(\tilde{\eta}) A^{-1} \zeta(\tilde{\eta})]. $$

Recalling from (29)

$$ \begin{aligned} \zeta(\tilde{\eta}) &= \nabla U_{\eta}(\tilde{\eta} + \eta^*) - \nabla U_{\eta}(\eta^*) \\ &= A[\sin(\tilde{\eta} + \eta^*) - \sin(\eta^*)], \end{aligned} $$

we have $\dot{\zeta}(\tilde{\eta}) = A\cos(\tilde{\eta} + \eta^*)\dot{\tilde{\eta}}$ and the derivative of $V$
can be written as

$$
\begin{align}
\dot{V} &= \zeta^{\top}(\tilde{\eta})\dot{\tilde{\eta}} + \frac{\alpha}{2} [\zeta^{\top}(\tilde{\eta}) \operatorname{diag}(\cos(\tilde{\eta} + \eta^*))\dot{\tilde{\eta}} + \nonumber \\
&\qquad \dot{\tilde{\eta}}^{\top} \operatorname{diag}(\cos(\tilde{\eta} + \eta^*))\zeta(\tilde{\eta})], \\
&= -\zeta^{\top}(\tilde{\eta}) \left( \mathcal{B}^{\top}\mathcal{B} + \frac{\alpha}{2} [\operatorname{diag}(\cos(\tilde{\eta} + \eta^*))\mathcal{B}^{\top}\mathcal{B} + \right. \nonumber \\
&\qquad \left. \mathcal{B}^{\top}\mathcal{B} \operatorname{diag}(\cos(\tilde{\eta} + \eta^*))] \right) \zeta(\tilde{\eta}). \tag{37}
\end{align}
$$

From (28) together with condition (35), we see that $\dot{V}$ is negative outside the equilibrium set $\mathcal{D}$ for any $\tilde{\eta} \in \mathbb{R}^m$. Hence, all conditions of Corollary 6 are satisfied.

By Proposition 8, in general, the set $\mathcal{D}$ may contain an unbounded continuum of equilibria. Nevertheless, this setting is already addressed in the proof of Corollary 6 given in [30]. Thus, with Assumptions 4, we conclude that all trajectories of the system (28) are bounded. $\square$

**Remark 14** Due to the fact that $\mathcal{B}^{\top}\mathcal{B} \ge 0$ and $x^{\top}\mathcal{B}^{\top}\mathcal{B}x = 0$ only if $x \in \mathcal{D}$, condition (35) is always satisfied for $\alpha = 0$. Additionally, for the acyclic case, we have $\mathcal{B}^{\top}\mathcal{B} > 0$. Thus, there always exists a sufficiently small $\alpha \in (0, 1)$, such that condition (35) is satisfied.

### 4.3 Global attractivity for the case of meshed network topologies

Having established conditions for global boundedness of solutions, we can state the global attractivity of the set $\mathcal{D}$, which means that the system (14) reaches frequency synchronization, see Definition 7.

**Theorem 15** Consider the system (28) with Assumptions 4 and 10. Suppose that the conditions (35) and (36) are satisfied. Then, the set $\mathcal{D}$ in (30) is globally attractive.

*Proof.* Under Assumptions 4 and 10, it follows that all trajectories of the system (28) are bounded if the conditions (35) and (36) are satisfied by invoking Proposition 13. Therefore, it remains to prove that all trajectories converge to the set $\mathcal{D}$.

For this purpose, we see from (37) that if (35) is satisfied, the derivative of the Leonov function $V$ is negative semidefinite. Thus, by invoking LaSalle's invariance principle [34], the trajectories of the system (28) converge to the set $\mathcal{D}$.

### 4.4 Almost global stability for the case of acyclic network topologies

In the acyclic case, additionally to the global attractiveness of the set $\mathcal{D}$, it is possible to conclude almost global stability of a subset of $\mathcal{D}$, which contains all asymptotically stable equilibria. In order to present this result, we consider the following assumption.

**Assumption 16** The graph $\mathcal{G}$ describing the topology of the Kuramoto oscillator network (8) is a connected tree.

Under Assumption 16, the pseudoinverse can be calculated as $\mathcal{B}^{\dagger} = (\mathcal{B}^{\top}\mathcal{B})^{-1}\mathcal{B}^{\top}$ and $\ker(\mathcal{B}) = \mathbf{0}_{n-1}$. Hence, the condition (23) results in

$$ \|A^{-1}(\mathcal{B}^{\top}\mathcal{B})^{-1}\mathcal{B}^{\top}\omega\|_{\infty} \le 1. \quad (38) $$

Then, also all equilibria are isolated. Furthermore, if and only if (38) is satisfied with strict inequality, the solution in (24) becomes

$$ \eta^* = \arcsin(A^{-1}(\mathcal{B}^{\top}\mathcal{B})^{-1}\mathcal{B}^{\top}\omega), \quad \hat{\eta}^* = \pi \mathbf{1}_m - \eta^*. \quad (39) $$

Likewise, the proposition below characterizes the (local) stability properties of the equilibria of the system (20).

**Proposition 17** Consider the system (20) under Assumptions 10 and 16. The equilibrium point $\eta^*$ given in (39) is locally exponentially stable. All other equilibria are unstable and the Jacobian of the dynamics (20) evaluated at any unstable equilibrium point has at least one eigenvalue with positive real part.

*Proof.* Consider the linearization of the system (20) around a point $\tilde{\eta} \in \mathbb{R}^m$, whose Jacobian is given by

$$ -\mathcal{B}^{\top}\mathcal{B} \left( \left. \frac{\partial}{\partial \eta} \mathbf{A}\sin(\eta) \right|_{\eta=\bar{\eta}} \right) = -\mathcal{B}^{\top}\mathcal{B}\mathbf{A}\operatorname{diag}(\cos(\bar{\eta})). \quad (40) $$

From Assumption 16, $\mathcal{B}^{\top}\mathcal{B}$ is a positive definite matrix and $A$ is a diagonal positive definite matrix. Therefore, if $\|\bar{\eta}\|_{\infty} < \frac{\pi}{2}$, as $\bar{\eta} = \eta^*$ in (39), then all eigenvalues of $\mathbf{A}\operatorname{diag}(\cos(\eta^*))$ are positive. Otherwise, if $\bar{\eta} = \eta^*$ in (39), then $\|\bar{\eta}\|_{\infty} \in (\frac{\pi}{2}, \frac{3\pi}{2})$ and all eigenvalues of $\mathbf{A}\operatorname{diag}(\cos(\eta^*))$ are negative. Thus, by invoking Lemma 1, we conclude that $\eta^*$ given in (39) is an exponentially stable equilibrium point, all other equilibria, which are obtained from permutations of $\eta^*$ and $\hat{\eta}^*$, are unstable and their Jacobian has at least one eigenvalue with positive real part. $\square$

In light of Proposition 17, we use $\eta^*$ in (39) to define the shifted variable $\tilde{\eta} = \eta - \eta^*$, such that the stable invariant set is given by

$$ \mathcal{D}_0 = \{\tilde{\eta} \in \mathbb{R}^m \mid \tilde{\eta} = 2r\pi \ \forall \ r \in \mathbb{Z}^m\}, \quad (41) $$

and the following corollary to Theorem 15 can be stated.

**Corollary 18** Consider the system (28) with Assumptions 4, 10 and 16. Suppose that the conditions (35) and (36) are satisfied. Then, the set $\mathcal{D}_0$ is almost globally asymptotically stable.

*Proof.* From Theorem 15, all trajectories converge to the set $\mathcal{D}$ if the conditions (35) and (36) hold. Therefore, it
remains to prove that almost all trajectories converge to the set $D_0$.

To show that $D_0$ is almost globally asymptotically stable, we use Proposition 17. This shows that the Jacobian of the dynamics (28) evaluated at any equilibrium point outside of the set $D_0$ has at least one eigenvalue with positive real part. By invoking [26, Proposition 11], we conclude that the region of attraction of any unstable equilibrium point has zero Lebesgue measure. Hence, for all initial conditions except a set of measure zero, the solutions of the system (28) converge asymptotically to $D_0$. □

**Remark 19** In the meshed case, due to a possible presence of a continuum of equilibria in $D$, separating the stable and unstable equilibria sets may be difficult. Furthermore, as a consequence of such a continuum of equilibria, the unstable invariant set may have an attraction domain with non-zero Lebesgue measure.

# 5 Interpretation of the global synchronization conditions

To illustrate the implications of the synchronization conditions, we simplify the conditions (35) and (36), focusing on the acyclic case. Although, this simplification may come at the expense of potentially more conservative conditions, they are going to help us to link the synchronization conditions with the underlying graph topological properties of the network. For their presentation, we introduce the Laplacian matrix of the graph $G$ given by, see [12],

$$ \mathcal{L} = BB^\top. \quad (42) $$

The graph $G$ is connected by assumption, hence, $\mathcal{L} \ge 0$ with a simple zero eigenvalue $\lambda_1(\mathcal{L}) = 0$ [12]. Furthermore, we denote by $\lambda_2(\mathcal{L})$ the first nonzero eigenvalue of $\mathcal{L}$, also known as the algebraic connectivity of the graph [12].

## 5.1 Maximal allowable stationary phase angle difference

From (36), we have $\|\eta^*\|_{\infty} < 0.5669$ for $\alpha = 0$. Recalling the change of coordinates (17), this restriction can be seen as the maximal allowable stationary phase angle difference (mod $2\pi$). This condition implies that the maximum distance for the general case with a meshed network topology is given by $\|\eta^*\|_{\infty} < 0.5669$ with the Leonov function (31). However, for the acyclic case, we can extend the maximal allowable stationary phase angle difference through $0.5669 \le \bar{x}^*(\alpha) \le 0.9 < 0.5\pi$ (see Figure 1) by using the parameter $\alpha$ in the Leonov function $V$ in (31).

As it can be seen in Figure 1, the value of the maximal allowable stationary phase angle difference given by $\bar{x}(\alpha)$

is increasing with respect to $\alpha \in [0, 1)$. Consequently, the condition (35) determines this maximum distance. Considering the relationship between $\mathcal{L}$ and $B^\top B$, we have the following corollary to Proposition 13 for the acyclic case.

**Corollary 20** Under Assumption 16, the inequality (35) is satisfied if

$$ \alpha < \sqrt{\lambda_2(\mathcal{L})/\lambda_{\max}(\mathcal{L})}, \quad (43) $$

with $\mathcal{L}$ given in (42).

*Proof.* We recall that Assumption 16 implies that $B^\top B$ is real positive definite. Let the symmetric matrix $\bar{B}^{\frac{1}{2}}$ be the unique real positive square root of $B^\top B$, i.e., $B^\top B = \bar{B}^{\frac{1}{2}}\bar{B}^{\frac{1}{2}}$ [15, Theorem 7.2.6]. Then, by performing a congruence transformation with $\bar{B}^{-\frac{1}{2}}$ and invoking Sylvester's law of inertia [6, Theorem 9.13] the matrix in (35) is positive definite iff

$$ \mathbf{I}_m + \frac{1}{2}\alpha (M(\phi) + M^\top(\phi)) > 0, \quad (44) $$

where we introduced the short-hand

$$ M(\phi) = \bar{B}^{-\frac{1}{2}} \mathrm{diag}(\cos(\phi)) \bar{B}^{\frac{1}{2}}. $$

By using the 2-norm properties, we have that

$$ \|M(\phi)\|_2 \le \|\bar{B}^{-\frac{1}{2}}\|_2 \|\bar{B}^{\frac{1}{2}}\|_2 = \sqrt{\frac{\lambda_{\max}(\bar{B}^\top B)}{\lambda_{\min}(\bar{B}^\top B)}}. $$

All nonzero eigenvalues of $B^\top B$ and $\mathcal{L}$ in (42) are identical [5, Proposition 4.4.9]. Thus, the inequality (35) is satisfied if

$$ 1 - \alpha \sqrt{\frac{\lambda_{\max}(\bar{B}^\top B)}{\lambda_{\min}(\bar{B}^\top B)}} = 1 - \alpha \sqrt{\frac{\lambda_{\max}(\mathcal{L})}{\lambda_2(\mathcal{L})}} > 0, $$

which is equivalent to (43). □

From Corollary 20, we see that condition (43) depends on the network topology, encompassed in the Laplacian matrix $\mathcal{L}$ in (42) and imposes a restriction on the heterogeneity of the network topology in terms of the distance of the eigenvalues of $\mathcal{L}$. These topological properties in turn restrict the magnitude of $\alpha$ and, thus, the maximal allowable stationary phase angle difference (mod $2\pi$).

## 5.2 Minimum coupling strength

Once the maximal allowable stationary phase angle difference has been determined, the minimum coupling
strength can be calculated. From (36) together with (24), we can see that the minimum coupling strength depends on the stationary states, the pseudoinverse $B^\dagger$ and the natural frequencies. Likewise, the inequality (36) depends on the already determined $\alpha$ from (35) or (43). By using the relationship between the eigenvalues of $\mathcal{L}$ and $B^\top B$, we have the following corollary to Proposition 13 for the acyclic case.

**Corollary 21** Under Assumption 16, the bound (36) is satisfied if

$$ \lambda_2(\mathcal{L})\lambda_{\min}(A) > \sqrt{m} \frac{\omega_{\Delta_{\max}}}{\sin(\bar{x}(\alpha))}, \quad (45) $$

with $\omega_{\Delta_{\max}} = |\omega_{\max} - \omega_{\min}|$.

*Proof.* We recall that the equilibrium point has to satisfy

$$ B^\top \omega = B^\top B A \sin(\eta^*) . $$

Under Assumption 16 and condition (36), the previous equality has a solution if and only if

$$ \|A^{-1}(B^\top B)^{-1}B^\top \omega\|_{\infty} < \sin(\bar{x}(\alpha)). \quad (46) $$

Since $A$ is a diagonal matrix and it is the first term in the matrix product on the left hand-side of (46), we have

$$ \|A^{-1}(B^\top B)^{-1}B^\top \omega\|_{\infty} \le \frac{\|(B^\top B)^{-1}\|_{\infty} \|B^\top \omega\|_{\infty}}{\lambda_{\min}(A)} . $$

Also, all nonzero eigenvalues of $B^\top B$ and $\mathcal{L}$ in (42) are identical [5, Proposition 4.4.9]. Thus, we obtain

$$ \|(B^\top B)^{-1}\|_{\infty} \le \sqrt{m} \|(B^\top B)^{-1}\|_{2} = \frac{\sqrt{m}}{\lambda_{\min}(B^\top B)} = \frac{\sqrt{m}}{\lambda_{2}(\mathcal{L})}, $$

which in combination with $\|B^\top \omega\|_{\infty} \le \omega_{\Delta_{\max}} = |\omega_{\max} - \omega_{\min}|$ yields

$$ \|A^{-1}(B^\top B)^{-1}B^\top \omega\|_{\infty} \le \frac{\sqrt{m\omega_{\Delta_{\max}}}}{\lambda_2(\mathcal{L})\lambda_{\min}(A)}. $$

Hence, (45) implies that condition (46) is satisfied. □

From condition (45), we can deduce that the minimum coupling strength depends on:

* $\omega_{\Delta_{\max}}$, which represents the worst case difference in the oscillator's natural frequencies. Here, we can note that the classical condition for homogeneous natural frequencies is recovered [16, 25], since then (45) is satisfied for any positive coupling strengths.

* $\lambda_2(\mathcal{L})$, which is the algebraic connectivity of the graph. We see that a lower intrinsic algebraic connectivity of the network, i.e., $\lambda_2(\mathcal{L})$, has to be compensated by larger coupling strengths, i.e., $\lambda_{\min}(A)$.

* $\sin(\bar{x}(\alpha))$, which is the sine of the maximal allowable stationary phase angle difference (mod $2\pi$). Thus, since the sine function is locally proportional to its argument, (45) implies that the minimum network coupling strength $\lambda_{\min}(A)$ has to scale inversely proportional to $\bar{x}(\alpha)$. This implies that smaller stationary angles $\eta^*$ require larger coupling strengths. This finding is consistent with Proposition 8.

We recall that the range of admissible stationary phase angles can be increased by maximizing $\alpha$ subject to (35). Thereby, the maximum admissible value for $|\eta^*|$ with our conditions is $\bar{x}(\alpha) = 0.9 \approx \pi/3.5$, see Fig. 1.

Additionally, by using the inequality (36), we can also design a coupling strength to obtain a final stationary phase angle difference (mod $2\pi$), which has to be less than the maximum one. However, to ensure synchronization to a solution with smaller stationary phase angle differences, larger coupling strengths are required.

# 6 Simulation Example

For illustrating our findings, we consider the example presented in [14, Example 1]. We recall that this example shows that the mere existence of equilibria for (14) does not ensure almost global synchronization. In this counterexample, three oscillators in an acyclic graph and the following parameters are considered:

$$ \mathcal{B} = \begin{bmatrix} 1_{1 \times 2} \\ -\mathbf{I}_2 \end{bmatrix}, \quad A = \mathbf{I}_2, \quad \omega = \begin{bmatrix} 1.9 \\ -0.951 \mathbf{1}_2 \end{bmatrix}. \quad (47) $$

For this set of parameters, the inequality (38) is satisfied and, hence, by Proposition 8, there exist equilibrium points for the system (20). However, we recall that this condition is not sufficient for ensuring global synchronization. Actually, with the parameters in (47), we can obtain the graph in Figure 2. In this figure, we can see that there exist cyclic trajectories between saddle equilibria, which are represented by cyan lines. As a consequence, we can observe unbounded invariant sets between the cycles. To generate this figure, we consider a random mesh of 10,000 initial conditions $\theta(0) \in S^3$ and we indeed find initial conditions, whose corresponding trajectories do not converge during the simulation (in red). Likewise, these initial conditions were shifted by $\pm 2\pi$ to obtain a wider view of the problem.

From Figure 2, we can see that the Kuramoto system with the parameters in (47) exhibits heteroclinic cycles and hence, Assumption 4 is not satisfied. We stress that verifying Assumption 4 is a very hard problem by itself.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Simulation example, $A = I_2$. Saddle points are denoted by $\star$, stable equilibria by $\star$ and unstable equilibria by $\star$.</figcaption></figure>

However, our numerical simulations show that these cycles can disappear with appropriate parameters, as is shown next. For this reason, although Assumption 4 cannot be guaranteed, we employ the conditions obtained in Proposition 13 and Theorem 15 to calculate conditions for the coupling strength matrix $A$. We recall that (35) and, thus, the upper bound for $\alpha$ merely depend on the graph topology. Consequently, at first, we calculate the maximum possible $\alpha$, such that condition (35) holds for all $\phi \in S^m$. For the current example with a precision<fnref n="1" /> of $10^{-3}$, this condition results in $\alpha < 0.866$. We observe that the obtained upper bound for $\alpha$ is very close to the maximum admissible value with our conditions, i.e., $\alpha = 1$ and $\bar{x}(1) = 0.9 \approx \pi/3.5$, see Fig. 1. By choosing $\alpha = 0.865$, we obtain $\bar{x}(\alpha) = 0.8159$ and condition (36) is satisfied if $\lambda_{\min}(A) > 1.3043$. We, hence, fix the coupling strength matrix to $A = 1.305I_2$.

For this set of parameters, we obtain the graph in Figure 3. We can see that despite the conditions in Proposition 13 and Theorem 15 being satisfied, there still exist cycles between some of the saddle equilibria of the system. Hence, Assumption 4 does not hold and the Leonov-based analysis on boundedness of solutions is inconclusive. Yet, we observe that compared to Figure 2 the size of the invariant sets in Figure 3 is strongly reduced.

Prompted by this observation, we slowly increase the magnitude of $A = aI_2$, while keeping $\alpha$ constant. By doing so, we find that the cycles disappear for $a \ge 1.3436$, as we can see in Figure 4. This corresponds to an increased magnitude of $a$ of just about 3% compared to the minimal value $a = 1.3043$ required for our proposed Leonov function to satisfy all conditions in Proposition 13 and Theorem 15. Thus, we conjecture that, while not being sufficient, for the considered example the derived conditions are close to being necessary for

avoiding the existence of cycles. Therefore, the obtained conditions in Theorem 13 are a good approximation to obtain almost global synchronization since nowadays, there exist no tools nor methods to formally exclude the existence of cycles.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Simulation example for $A = 1.305I_2$. Saddle points are denoted by $\star$, stable equilibria by $\star$ and unstable equilibria by $\star$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Simulation example for $A = 1.3436I_2$. Saddle points are denoted by $\star$, stable equilibria by $\star$ and unstable equilibria by $\star$.</figcaption></figure>

As a consequence of that our conditions are too close to the boundary of stability, we can use the more conservative bounds in (45). Using this inequality, we obtain $\lambda_{\min}(A) > 5.5312$, which guarantees the almost global stability for this example since there are no homoclinic orbits for $\lambda_{\min}(A) \ge 1.3436$. Therefore, despite condition (38) could be more conservative, this condition may guarantee the absence of homoclinic orbits in the system and the almost global stability of the set $D_0$. To illustrate this case, we consider the case of $n$ oscillators

<footnote n="1">¹ The calculations in [24] are done with a precision of $10^{-2}$.</footnote>
whose parameters satisfy $A = a_0 \mathbf{I}_{n-1}$ with $a_0 > 0$,

$$ \mathcal{B} = \begin{bmatrix} \mathbf{1}_{1 \times n-1} \\ -\mathbf{I}_{n-1} \end{bmatrix} \in \mathbb{R}^{n \times n-1}, \omega = \begin{bmatrix} n-1 - \frac{1}{10} \\ -\mathbf{1}_{n-1} + \frac{\mathbf{1}_{n-1}}{10(n-1)} \end{bmatrix}. $$

Using these parameters, we can obtain the gains in Table 1². In this table, we can see that the estimated gain convergence is lower for the all considered cases in this work. Thus, using equation (45) could be a good way to ensure the almost global stability of the equilibria set.

<table><caption>Table 1
Comparison of gains for different number of oscillators</caption><thead><tr><th>N</th><th>Eq (36)</th><th>Eq (45)</th><th>Convergence</th></tr></thead><tbody><tr><td>3</td><td>a<sub>0</sub> > 1.3043</td><td>a<sub>0</sub> ≥ 5.5312</td><td>a<sub>0</sub> ≈ 1.344</td></tr><tr><td>4</td><td>a<sub>0</sub> > 1.3609</td><td>a<sub>0</sub> ≥ 9.4285</td><td>a<sub>0</sub> ≈ 1.464</td></tr><tr><td>5</td><td>a<sub>0</sub> > 1.4202</td><td>a<sub>0</sub> ≥ 14.2025</td><td>a<sub>0</sub> ≈ 1.660</td></tr></tbody></table>

## 7 Conclusions

We have presented sufficient conditions for global frequency synchronization in meshed and acyclic networks of heterogeneons Kuramoto oscillators. The main result has been proven by using the recently proposed Leonov function approach for systems with state periodic dynamics. Compared to standard Lyapunov functions, the requirements for Leonov functions are more relaxed. By exploiting these features, we have constructed a suitable Leonov function candidate for the particular case of Kuramoto oscillators. To this end, we also have reformulated the problem of synchronization via a change of coordinates into a standard stability analysis problem. We have then provided lemmata, by which we can ensure that all conditions in Corollary 6 are satisfied. This step allows us to establish global boundedness of solutions.

The global attractiveness analysis is completed by one additional step, which consists of establishing necessary and sufficient conditions for the existence of equilibria. By combining them with the Leonov function analysis and the LaSalle's invariance principle, we have obtained sufficient conditions for global frequency synchronization of the network. Additionally, in the acyclic case, we derive the local stability/instability properties for all equilibria to show that proposed conditions ensure the almost global synchronized motion of the network.

In the acyclic case, we have also discussed, how the derived conditions can be associated with the main features of the network properties. Particularly, we have shown that the parameter $\alpha$, which provides an upper bound for the admissible steady-state angle differences,

is related to the network topology by means of the ratio of the minimum over the maximum eigenvalues of the graph Laplacian. The parameter $\alpha$ and the coupling strengths $A$, then constrain the maximum admissible difference between natural frequencies.

Finally, the results have been illustrated by considering the example presented in [14]. For this example, we have redesigned the matrix of coupling strengths $A$, such that almost global synchronization can be obtained. Via this example we have also examined the existence of cycles connecting different saddle equilibria of the Kuramoto system and these cycles produce an unbounded invariant set. We have found that for certain parameter configurations, such cycles may exist even if all conditions on the proposed Leonov function are satisfied. However, by only slightly tightening the parametric requirements obtained via the Leonov analysis, the cycles disappear. This may indicate that the synchronization conditions derived from a Leonov-based analysis are not overly conservative, but instead close to necessary to preclude the existence of cycles. We also emphasize that there no exists any formal method to find conditions to avoid the presence of cycles.

Current and future work is geared towards extending the analysis to other classes of oscillators, relaxing assumptions for the general case and deriving tools to formally exclude the existence of cycles.

## Acknowledgements

This work is partially supported by the ANR-DFG project SyNpID - Project number 446182476.

## References

[1] Juan A. Acebrón, L. L. Bonilla, Conrad J. Pérez Vicente, Félix Ritort, and Renato Spigler. The Kuramoto model: A simple paradigm for synchronization phenomena. *Reviews of Modern Physics*, 77:137–185, 2005.

[2] David Angeli and Denis Efimov. Characterizations of input-to-state stability for systems with multiple invariant sets. *IEEE Transactions on Automatic Control*, 60(12):3242–3256, 2015.

[3] Murat Arcak. Passivity as a design tool for group coordination. *IEEE Transactions on Automatic Control*, 52(8):1380–1390, 2007.

[4] Viktor Avrutin, Björn Schenke, and Laura Gardini. Calculation of homoclinic and heteroclinic orbits in 1d maps. *Communications in Nonlinear Science and Numerical Simulation*, 22(1):1201–1214, 2015.

[5] Dennis S. Bernstein. *Matrix Mathematics: Theory, Facts, and Formulas: Second Edition*. Princeton University Press, 2 edition, 2009.

<footnote n="2">² For $n = 4, 5$, the gains for convergence were estimated with a random mesh of 1000 initial conditions and increasing the gain until all trajectories converge to the equilibria set.</footnote>
[6] James B. Carrel. *Groups, Matrices, and Vector Spaces: A Group Theoretic Approach to Linear Algebra*. Springer, 2017.

[7] Florian Dörfler and Francesco Bullo. On the critical coupling for Kuramoto oscillators. *SIAM Journal on Applied Dynamical Systems*, 10(3):1070–1099, 2011.

[8] Florian Dörfler and Francesco Bullo. Synchronization in complex networks of phase oscillators: A survey. *Automatica*, 50(6):1539–1564, 2014.

[9] Florian Dörfler, Michael Chertkov, and Francesco Bullo. Synchronization in complex oscillator networks and smart grids. *Proceedings of the National Academy of Sciences*, 110(6):2005–2010, 2013.

[10] Denis Efimov and Johannes Schiffer. On boundedness of solutions of state periodic systems: A multivariable cell structure approach. *IEEE Transactions on Automatic Control*, 64(10):4094–4104, 2019.

[11] Denis Efimov, Johannes Schiffer, Nikita Barabanov, and Romeo Ortega. A relaxed characterization of ISS for periodic systems with multiple invariant sets. *European Journal of Control*, 37:1–7, 2017.

[12] Chris Godsil and Gordon F. Royle. *Algebraic graph theory*, volume 207. Springer Science & Business Media, 2001.

[13] Andrey Gushchin, Enrique Mallada, and Ao Tang. Synchronization of heterogeneous Kuramoto oscillators with arbitrary topology. In *American Control Conference*, pages 637–644, Chicago, IL, USA, 2015.

[14] Andrey Gushchin, Enrique Mallada, and Ao Tang. Synchronization of heterogeneous Kuramoto oscillators with graphs of diameter two. In *49th Annual Conference on Information Sciences and Systems*, pages 1–6, Baltimore, MD, USA, 2015.

[15] Roger A. Horn and Charles R. Johnson. *Matrix Analysis*. Cambridge University Press, 2 edition, 2013.

[16] A. Jadbabaie, N. Motee, and M. Barahona. On the stability of the Kuramoto model of coupled nonlinear oscillators. In *American Control Conference*, volume 5, pages 4296–4301, Boston, MA, USA, 2004.

[17] Martin Kassabov, Steven H. Strogatz, and Alex Townsend. A global synchronization theorem for oscillators on a random graph. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 32(9):093119, 2022.

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

[19] Gennady A. Leonov. On the boundedness of the trajectories of phase systems. *Siberian Mathematical Journal*, 16(3):491–495, 1974.

[20] Dhagash Mehta, Noah S. Daleo, Florian Dörfler, and Jonathan D. Hauenstein. Algebraic geometry of the Kuramoto model: Equilibria and stability analysis. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 25(5):053103, 2015.

[21] Jesús Mendoza-Avila, Denis Efimov, José Mercado Uribe, and Johannes Schiffer. Design of controls for boundedness of trajectories of multistable state periodic systems. In *62nd IEEE Conference on Decision and Control*, pages 6647–6652, Singapore, 2023.

[22] Jesús Mendoza-Avila, Angel Mercado-Uribe, Denis Efimov, and Johannes Schiffer. Design of controls for ISS and integral ISS stabilization of multistable state periodic systems. *IFAC-PapersOnLine*, 56(2):4472–4477, 2023.

[23] Jesús Mendoza-Ávila, Denis Efimov, Angel Mercado-Uribe, and Johannes Schiffer. On relaxed conditions of integral ISS for multistable periodic systems. In *61st IEEE Conference on Decision and Control*, Cancun, Mexico, 2022. IEEE.

[24] Angel Mercado-Uribe, Jesús Mendoza-Avila, Denis Efimov, and Johannes Schiffer. A leonov function for almost global synchronization conditions in acyclic networks of heterogeneous kuramoto oscillators. *IFAC-PapersOnLine*, 56(2):9505–9510, 2023.

[25] Pablo Monzon and Fernando Paganini. Global considerations on the Kuramoto model of sinusoidally coupled oscillators. In *44th IEEE Conference on Decision and Control*, pages 3923–3928, Seville, Spain, 2005. IEEE.

[26] Pablo Monzon and Rafael Potrie. Local and global aspects of almost global stability. In *45th IEEE Conference on Decision and Control*, pages 5120–5125, San Diego, CA, USA, 2006. IEE.

[27] Erik J. Noldus. New direct Lyapunov-type method for studying synchronization problems. *Automatica*, 13(2):139–151, 1977.

[28] Arkady Pikovsky, Michael Rosenblum, and Jürgen Kurths. *Synchronization: A Universal Concept in Nonlinear Sciences*. Cambridge University Press, 2003.

[29] Francisco A. Rodrigues, Thomas K. DM. Peron, Peng Ji, and Jürgen Kurths. The kuramoto model in complex networks. *Physics Reports*, 610:1–98, 2016.

[30] Johannes Schiffer and Denis Efimov. Strong and weak Leonov functions for global boundedness of state periodic systems. *IEEE Transactions on Automatic Control*, 68(12):7958–7965, 2023.

[31] Johannes Schiffer, Denis Efimov, and Romeo Ortega. Global synchronization analysis of droop-controlled microgrids—a multivariable cell structure approach. *Automatica*, 109:108550, 2019.

[32] Johannes Schiffer, Romeo Ortega, Alessandro Astolfi, Jörg Raisch, and Tevfik Sezi. Conditions for stability of droop-controlled inverter-based microgrids. *Automatica*, 50(10):2457–2469, 2014.

[33] Sebastian Trip, Mathias Bürger, and Claudio De
Persis. An internal model approach to (optimal) fre-
quency regulation in power grids with time-varying
voltages. Automatica, 64:240–253, 2016.

[34] Arjan van der Schaft. *L2-Gain and passivity techniques in nonlinear control*. Springer, Netherlands, 2000.

[35] George Weiss, Florian Dörfler, and Yoash Lebron. A stability theorem for networks containing synchronous generators. *Systems & Control Letters*, 134:104561, 2019.

[36] Arthur T. Winfree. *The Geometry of Biological Time*. Springer, USA, 1980.

[37] Lijun Zhu and David J Hill. Synchronization of Kuramoto oscillators: A regional stability framework. *IEEE Transactions on Automatic Control*, 65(12):5070–5082, 2020.