# A Leonov Function for Almost Global Synchronization Conditions in Acyclic 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. A Leonov Function for Almost Global Synchronization Conditions in Acyclic Networks of Heterogeneous Kuramoto Oscillators. IFAC 2023 - 22nd IFAC World Congress, Jul 2023, Yokohama, Japan. hal-04176145

HAL Id: hal-04176145

https://inria.hal.science/hal-04176145

Submitted on 2 Aug 2023

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
# A Leonov Function for Almost Global Synchronization Conditions in Acyclic Networks of Heterogeneous Kuramoto Oscillators

Angel Mercado-Uribe* Jesús Mendoza-Ávila**

Denis Efimov** Johannes Schiffer*,**

* Brandenburgische Technische Universität Cottbus - Senftenberg, Germany (e-mail: {mercadou,schiffer}@b-tu.de).

** Inria, Univ. Lille, CNRS UMR 9189 - CRISTAL, F-59000 Lille, France (e-mail: {jesus.mendoza-avila,denis.efimov@inria.fr}).

*** Fraunhofer Research Institution for Energy Infrastructures and Geothermal Systems (IEG), 03046 Cottbus, Germany

**Abstract:** Sufficient conditions for almost global synchronization in acyclic networks of Kuramoto oscillators with heterogeneous coupling strengths and natural frequencies are presented. The result is established by employing the recently developed Leonov function framework for systems whose dynamics are periodic for all state variables. The synchronization property is accompanied by necessary and sufficient conditions to guarantee the existence of equilibria. The implications of these conditions on the network topology, the oscillator's coupling strengths and natural frequencies are discussed. Finally, the results are illustrated via a numerical example.

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

## 1. INTRODUCTION

The analysis of synchronization phenomena in networks of coupled oscillators is a long-standing fascinating, and pervasive topic in many different areas and disciplines, such as neuroscience, biology, chemical, physics, mathematics and engineering (Acebrón et al., 2005; Dörfler and Bullo, 2014) as well as, more recently, computational communities (Kassabov et al., 2022). The most prominent oscillator representation is the Kuramoto model (Kuramoto, 1975; Acebrón et al., 2005). This model consists of first-order oscillators with intrinsic natural frequencies coupled via weighted sinusoidal functions of their phase differences. The Kuramoto oscillator dynamics are sufficiently complex to exhibit a highly nontrivial and rich behavior. As a result, the Kuramoto model has become an essential benchmark for studying synchronization (Acebrón et al., 2005; Dörfler and Bullo, 2014; Rodrigues et al., 2016).

One of the main objectives in the analysis of Kuramoto oscillators is to identify conditions such that the oscillators reach frequency synchronization. Regrettably, the nonlinear terms imply that the system has infinite equilibria and significantly complicate the analysis. Consequently, deriving synchronization conditions for heterogeneous coupling weights and natural frequencies is a complex task, particularly in a *global* setting. So, simplified parametrizations and *local* synchronization are studied in the literature, see (Dörfler and Bullo, 2014; Rodrigues et al., 2016).

Under the common assumption of homogeneous oscillators, i.e., with identical natural frequencies, (Jadbabaie and others., 2004; Monzon and Paganini, 2005) show that the oscillators achieve global synchronization is obtained for any positive value of coupling strengths. Naturally, the analysis is more complex for heterogeneous (non-identical) Kuramoto oscillators. Hence, most available results are restricted to *local* synchronization, see (Dörfler and Bullo, 2011; Weiss et al., 2019; Zhu and Hill, 2020) and the surveys (Dörfler and Bullo, 2014; Rodrigues et al., 2016). For non-uniform natural frequencies but identical coupling strengths and graphs with diameter less than or equal to 2, i.e., the maximum shortest path from one node to another is less than 3, an analytic bound to reach synchronization is presented in (Gushchin et al., 2015b). These results are extended to non-identical coupling strengths in (Gushchin et al., 2015b) and arbitrary topologies but identical coupling strengths in (Gushchin et al., 2015a).

The prevalent analysis technique in the literature is the Lyapunov theory. However, in a heterogeneous setting, a Lyapunov-based analysis usually only yields *local*, but no *global* synchronization results. The main technical challenge is establishing global boundedness of solutions, see also (Gushchin et al., 2015b). One approach to analyze global boundedness properties for state periodic systems is the cell structure approach, which relaxes the definiteness and periodicity requirements of standard Lyapunov-base methods by exploiting the periodicity of the system dynamics. It was first presented for systems with a scalar periodic state in (Leonov, 1974; Noldus, 1977).

<footnote>* This work is partially supported by the ANR-DFG project SyNPIID - Project number 446182476.</footnote>
This seminal result was extended to systems with multiple periodic variables in (Efimov and Schiffer, 2019; Schiffer and Efimov, 2023). Also, this result was employed to study Input-to-State Stability in (ISS) in (Efimov et al., 2017) and integral ISS (iISS) in (Mendoza-Ávila et al., 2022).

This work's main contributions lie in analyzing global synchronization properties in acyclic networks of heterogeneous Kuramoto oscillators. As demonstrated by the counterexample in (Gushchin et al., 2015b), even for an acyclic network topology, the existence of locally stable equilibria does not imply any kind of global stability. To reach this objective, instead of standard Lyapunov techniques, we employ the recently developed Leonov function framework (Efimov and Schiffer, 2019; Schiffer and Efimov, 2023). In summary, the main contributions of the paper are:

* Necessary and sufficient conditions to ensure the existence of a synchronized motion are provided.

* By constructing a suitable Leonov function candidate, sufficient conditions for global boundedness of trajectories of the coupled oscillators are presented.

* Under the assumption of local stability, almost global synchronization of all oscillators is proven by using a LaSalle argument.

The obtained conditions are discussed and illustrated via an example, which was used in (Gushchin et al., 2015b).

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

## 2. PRELIMINARIES

### 2.1 Notation

We denote the real and integer numbers sets by $\mathbb{R}$ and $\mathbb{Z}$. Let $\mathbb{R}_+ = \{s \in \mathbb{R} | s \ge 0\}$ as well as $\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|)$. $\mathbf{I}_k \in \mathbb{R}^{k \times k}$ represents the identity matrix. Let $M \in \mathbb{R}^{n \times m}$ be a matrix. Then, its induced infinity norm is defined as $\|M\|_{\infty} = \max_{\|x\|_{\infty}=1} (\|Mx\|_{\infty})$. For a matrix $A \in \mathbb{R}^{n \times n}$, $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ denote its smallest and largest eigenvalue. 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$ and $\cos(x) = [\cos(x_1), ..., \cos(x_k)]^\top$ and diag($x$) to be 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 Global boundedness of solutions of state periodic systems via Leonov functions

We recall some notions of Leonov functions based on (Efimov and Schiffer, 2019; Schiffer and Efimov, 2023). 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 $C^1$ map with $f(\mathbf{0}_n) = \mathbf{0}_n$.

This work focuses on systems whose dynamics are periodic for all state variables by considering the assumption below.

**Assumption 1.** 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$, for all $x \in \mathbb{R}^n$.

By Assumption 1, the dynamics (1) have several equilibria. So, we introduce two auxiliary sets. Let $\mathcal{W} \subset \mathbb{R}^n$ be an $n^{th}$ sphere, such that for some $c \in (0, 2\pi)$ and $\bar{c} \in (\pi, 2\pi)$,

$$ \{x \in \mathbb{R}^n | 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\}, $$

where the set $\mathcal{U}$ includes all equilibria by shifting the one at the origin. However, there are also other equilibria of (1) not included in $\mathcal{U}$.

Then, we can introduce an appropriated definition of a strong Leonov function for system (1), under Assumption 1, as follows, see (Schiffer and Efimov, 2023):

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

$$ \inf_{x \in \mathcal{W}} V(x) > 0, \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 \forall 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) $$

By using Definition 1, the following theorem is stated.

**Theorem 1.** (Schiffer and Efimov, 2023). If for the system (1) under Assumption 1, 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$.

By defining the sets

$$ \Omega = \{x \in \mathbb{R}^n | V(x) \le 0\}, \quad \mathcal{D} = \{x \in \mathbb{R}^n | f(x) = \mathbf{0}_n\}, \quad (5) $$

the following corollary can be stated, which is instrumental in performing the global synchronization analysis.

**Corollary 1.** (Schiffer and Efimov, 2023). If for the system (1) under Assumption 1, a Leonov function with

$$ V(0) \le 0, \inf_{x \in \mathcal{U} \setminus \mathcal{U}_0} V(x) \le 0, \inf_{x \in \mathcal{W}} V(x) > 0 \qquad (6) $$

exists and satisfies 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\}, \qquad (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 heterogeneous Kuramoto oscillators

We consider a network of first-order heterogeneous Kuramoto oscillators with $N > 1$ states, see (Dörfler and Bullo, 2014). The network topology is described by an undirected and connected graph $\mathcal{G} = (\mathcal{N}, \mathcal{E})$, where the set of network nodes is denoted by $\mathcal{N} = \{1, 2, ..., N\}$ and the set of edges by $\mathcal{E} = \{e_1, ..., 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}_{\ge 0} \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 \ne i, a_{ik} \ne 0\}$.

Regarding the topology, we make the assumption below.

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

Considering Assumption 2 and defining $m = N - 1$, the matrix of coupling strengths can be defined as

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

Likewise, as a result of Assumption 2, the node-edge incidence matrix $\mathcal{B}$ has full column rank. Consequently, its (Moore-Penrose) pseudo-inverse $\mathcal{B}^\dagger$ is given by

$$ \mathcal{B}^\dagger = (\mathcal{B}^\top \mathcal{B})^{-1} \mathcal{B}^\top \quad (10) $$

where $\mathcal{B}^\dagger$ is a left-inverse of $\mathcal{B}$, i.e., $\mathcal{B}^\dagger \mathcal{B} = \mathbf{I}_m$.

To write the network dynamics (8) compactly, it is convenient to introduce the potential function $U : \mathbb{R}^N \to \mathbb{R}$,

$$ U(\theta) = -\mathbf{1}_m^\top \mathbf{A} \cos (\mathcal{B}^\top \theta) . \quad (11) $$

The potential $U(\theta)$ and its gradient,

$$ \nabla U(\theta) = \mathcal{B} \mathbf{A} \sin(\mathcal{B}^\top \theta), \quad (12) $$

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

$$ U(\theta + \alpha \mathbf{1}_N) = U(\theta), \quad \nabla U(\theta + \alpha \mathbf{1}_N) = \nabla U(\theta), \quad (13) $$

because 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) = \mathbf{0}_N^\top. \quad (14) $$

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

$$ \dot{\theta} = \omega - \nabla U(\theta). \quad (15) $$

Similarly to (Schiffer et al., 2014, 2019), we introduce the notion below.

**Definition 2.** The system (15) admits a synchronized motion 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^s(t)), \quad (16) $$

where $\omega^* \in \mathbb{R}$ and $\theta_0^s \in \mathbb{R}^n$, such that

$$ |\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 (14) and (16), we obtain $\mathbf{1}_N^\top \mathcal{B} \mathbf{1}_N \omega^* = \mathbf{1}_N^\top \omega$. Thus, the synchronization frequency $\omega^*$ is determined by the average of the natural frequencies $\omega$ (Dörfler et al., 2013), 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 (17) $$

## 3.2 Change of coordinates

Inspired by (Arcak, 2007; Trip et al., 2016), we perform the following change of coordinates

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

It follows from (11) that

$$ U(\eta) = -\mathbf{1}_m^\top \mathbf{A} \cos(\eta), \quad \nabla U(\eta) = \mathbf{A} \sin(\eta), \quad (19) $$

and thus from (12) that

$$ \nabla U(\theta) = \mathcal{B} \nabla U(\eta). \quad (20) $$

Hence, in the new coordinates, the dynamics (15) become

$$ \dot{\eta} = \bar{\omega} - \mathcal{B}^\top \mathcal{B} \nabla U(\eta) = \bar{\omega} - \mathcal{B}^\top \mathcal{B} \mathbf{A} \sin(\eta), \quad (21) $$

where $\bar{\omega} = \mathcal{B}^\top \omega$ represents the differences in natural frequencies between any pair of coupled oscillators. We recall that by Assumption 2, the incidence matrix $\mathcal{B}$ has full column rank and hence $\mathcal{B}^\top \mathcal{B}$ is positive definite.

Clearly, any synchronized motion of the system (15) corresponds to an equilibrium of (21), i.e.

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

So, we focus our analysis on the dynamics (21). Therefore, the main objective lies in identifying suitable coupling strengths $\mathbf{A}$ and natural frequency differences $\bar{\omega}$, such that (21) possesses an almost globally stable equilibrium set.

## 3.3 Equilibrium points

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

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

Since the sinusoidal function is bounded, the existence of equilibria is not always possible. Therefore, we provide a necessary and sufficient condition for the existence of equilibria of the system (21). This result is a direct corollary of the derived for second-order microgrid systems in (Schiffer et al., 2019, Proposition 4). Also, related results have been reported using the pseudo-inverse of the network Laplacian in (Dörfler et al., 2013; Simpson-Porco et al., 2013) and (Manik et al., 2017, Corollary 2).

**Proposition 1.** Consider the system (21) with Assumption 2. The system (21) possesses equilibria, i.e., equation (22) has a solution, if and only if

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

Then, also all equilibria are isolated. Furthermore, if and only if (23) is satisfied with strict inequality, then the equilibria are given 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}^\top \mathcal{B})^{-1} \bar{\omega}), \quad \hat{\eta}^* = \pi \mathbf{1}_m - \eta^*. \quad (24) $$

*Proof.* Any equilibrium $\eta^* \in \mathbb{R}^m$ of (21) belongs to (22). By using algebra, we obtain the following system of equations

$$ A^{-1}(\mathcal{B}^\top \mathcal{B})^{-1} \bar{\omega} = \mathbf{sin}(\eta^*). \quad (25) $$

Clearly, (25) has solutions iff (23) is satisfied and each equation has the solutions given by (24). $\square$

## 3.4 Error dynamics

Now, we consider the error dynamics with respect to a equilibrium point. So, the following assumption is stated.

**Assumption 3.** Condition (23) is satisfied with strict inequality. The equilibrium point $\eta^*$ in (24), which satisfies $\|\eta^*\|_{\infty} < \frac{\pi}{2}$ (mod $2\pi$), is locally exponentially stable. All other equilibria, which satisfy $\|\eta^*\|_{\infty} \in (\frac{\pi}{2}, \frac{3\pi}{2})$ (mod $2\pi$) are unstable and the Jacobian of the dynamics (21) evaluated at any unstable equilibrium point has at least one eigenvalue with positive real part.
With Assumption 3 and by defining the variable $\tilde{\eta} = \eta - \eta^*$, the error dynamics of (21) results in

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

with

$$ \zeta(\tilde{\eta}) = \nabla U(\tilde{\eta} + \eta^*) - \nabla U(\eta^*) = A[\mathbf{sin}(\tilde{\eta} + \eta^*) - \mathbf{sin}(\eta^*)]. $$

## 4. SUFFICIENT CONDITIONS FOR ALMOST GLOBAL STABILITY

In this section, we derive sufficient conditions for almost global stability of (26). To establish the result, we propose the following candidate Leonov function

$$ V(\tilde{\eta}, \eta^*) = U(\tilde{\eta} + \eta^*) - U(\eta^*) - \nabla U^\top(\eta^*) \tilde{\eta} + 0.5\alpha \zeta^\top(\tilde{\eta}) A^{-1} \zeta(\tilde{\eta}) = \sum_{j=1}^{m} a_j V(\tilde{\eta}_j, \eta_j^*), \quad (27) $$

where $\alpha \in [0, 1)$ is a design parameter and

$$ V(\tilde{\eta}_j, \eta_j^*) = -\cos(\tilde{\eta}_j + \eta_j^*) + \cos(\eta_j^*) - \tilde{\eta}_j \sin(\eta_j^*) + 0.5\alpha (\sin(\tilde{\eta}_j + \eta_j^*) - \sin(\eta_j^*))^2. \quad (28) $$

### 4.1 Preliminary lemmata

A crucial point in the application of the Leonov function framework, see Section 2.2, is to establish that the Leonov function candidate $V$ in (27) is strictly positive on the set $W$ defined in (2). To this end, two lemmas are stated.

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

$$ \begin{aligned} S_1 &= \{x \in \mathbb{R} \mid x = 2r\pi \quad \forall r \in \mathbb{Z}\}, \\ S_2 &= \{x \in \mathbb{R} \mid x = \pi - 2x^* + 2r\pi \quad \forall r \in \mathbb{Z}\}. \end{aligned} \quad (29) $$

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

*Proof.* This proof is excluded due to page limitations. □

Lemma 1 shows that if $\|\eta^*\|_{\infty} < \frac{\pi}{2}$, the Leonov function $V$ in (27) is rising near to the origin and has a maximum before $\|\eta\|_{\infty} = \pi$ for any nonzero $\eta^* \in \mathbb{R}^m$. Since $\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 that any $\eta_j^*$, $j \in \mathcal{N}$, can take, fix $x = \pm\pi$, let $x^* \in \mathbb{R}$ and consider the function

$$ \mathcal{F}_V(x^*) = \mathcal{V}(\pm\pi, x^*) = 2\cos(x^*) + 2\alpha \sin^2(x^*) - (\pm\pi) \sin(x^*) \quad (30) $$

along with the following lemma.

**Lemma 2.** Let $\alpha \in [0, 1)$. 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.* This proof is excluded due to page limitations. □

Since the function $\mathcal{F}_V$ in (30) is analytically complicated, it is in general not possible to derive an 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(\bar{x}^*(\alpha)) = 0$.</figcaption></figure>

## 4.2 Global boundedness of trajectories

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

**Theorem 2.** Consider the system (26) with Assumptions 2 and 3. Let $\alpha \in [0, 1)$ be such that

$$ B^\top B + \frac{1}{2}\alpha [\operatorname{diag}(\cos(\phi)) B^\top B + B^\top B \operatorname{diag}(\cos(\phi))] > 0, \quad (31) $$

for all $\phi \in \mathbb{S}^m$. Then, all solutions of the system (26) are bounded if

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

*Proof.* The claim is established by using Corollary 1 with the candidate Leonov function $V$ in (27). 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(\tilde{\eta} + \eta^*)$ as well as $\zeta(\tilde{\eta} + \eta^*)$ are bounded functions in $\tilde{\eta}$, the first two conditions in (6) are satisfied. Lemma 1 together with Lemma 2 ensure that (32) holds, so the last condition in (6) is also verified. Hence, all requirements in (6) are satisfied.

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

$$ \dot{V} = -\zeta^\top (B^\top B + \frac{1}{2}\alpha [\operatorname{diag}(\cos(\eta)) B^\top B + B^\top B \operatorname{diag}(\cos(\eta))]) \zeta, \quad (33) $$

where we recall that $\eta = \tilde{\eta} + \eta^*$. Furthermore, from (26) together with condition (31), we see that $\dot{V}$ is negative outside the equilibrium set. Hence, (7) is also satisfied. Therefore, by invoking Corollary 1, we conclude that all trajectories of the system (26) are bounded. □

**Remark 1.** Since $B^\top B > 0$, there always is a sufficiently small $\alpha$, such that condition (31) is satisfied.

## 4.3 Main result: Almost global stability

Having established conditions for global boundedness of solutions, we are now ready to present our main result. For this purpose, we introduce the set

$$ \mathcal{D}_0 = \{\tilde{\eta} \in \mathbb{R}^m \mid \tilde{\eta} = 2r\pi \mathbf{1}_m \quad \forall r \in \mathbb{Z}\}. \quad (34) $$

**Theorem 3.** Consider the system (26) with Assumptions 2 and 3. Suppose the conditions (31) and (32) are satisfied. Then, the set $\mathcal{D}_0$ in (34) is almost globally asymptotically stable, i.e., for all initial conditions, except a set of measure zero, the solutions of the system (26) converge to $\mathcal{D}_0$.

*Proof.* Under Assumptions 2 and 3, by Theorem 2, if (31) and (32) are satisfied, all trajectories of the system (26) are bounded. Therefore, we must prove that almost all trajectories converge to the set $\mathcal{D}_0$.

For this purpose, we see from (33) that if (31) is satisfied, the derivative of the Leonov function $V$ is negative semidefinite. Thus, by invoking LaSalle's invariance principle (van der Schaft, 2000), the trajectories of the system (26) asymptotically converge to the set

$$ \mathcal{D}_0 \subset \{\tilde{\eta} \in \mathbb{R}^m \mid \zeta(\tilde{\eta}) = \mathbf{0}_m\}. $$

Thus, to show that $\mathcal{D}_0$ is almost globally asymptotically stable, note that Assumption 3 implies that the Jacobian of the dynamics (26) evaluated at any equilibrium point outside of the set $\mathcal{D}_0$ has at least one eigenvalue
with positive real part. Therefore, invoking (Monzon and Potrie, 2006, Proposition 11), we conclude that the region of attraction of any unstable equilibrium point has zero Lebesgue measure. Hence, for almost all initial conditions, the solutions of the system (26) converge to $\mathcal{D}_0$. $\square$

## 4.4 Interpretation of our global synchronization conditions

Next, we simplify our main synchronization conditions (31) and (32) in order to illustrate their implications. Naturally, these simplifications result in more conservative conditions. However, they will help us to link the synchronization conditions with the underlying graph topological properties of the network. For their presentation, let us introduce the Laplacian matrix of the graph $\mathcal{G}$ given by, see (Godsil and Royle, 2001),

$$ \mathcal{L} = \mathcal{B}\mathcal{B}^{\top}. \quad (35) $$

By Assumption 2, $\mathcal{G}$ is connected, then $\mathcal{L} \ge 0$ with a simple zero eigenvalue $\lambda_1 = 0$, (Godsil and Royle, 2001). Furthermore, we denote by $\lambda_2(\mathcal{L})$ the first nonzero eigenvalue of $\mathcal{L}$, also known as the algebraic connectivity of the graph, (Godsil and Royle, 2001).

Considering the previous relationship between $\mathcal{L}$ and $\mathcal{B}^{\top}\mathcal{B}$, we have the following two corollaries to Theorem 2.

**Corollary 2.** The bound (32) is satisfied if

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

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

*Sketch of the proof* Under the condition (32), the equality $\bar{\omega} = \mathcal{B}^{\top}\mathcal{B}\text{Asin}(\eta^*)$, which contains all equilibrium points, has a solution if and only if

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

Using algebra and bounds for norms, we can verify that if (36) holds, the condition (37) is satisfied. $\square$

**Corollary 3.** The inequality (31) is satisfied if

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

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

*Proof.* This proof is excluded due to page limitations. $\square$

Concerning condition (36), we observe that by recalling the change of coordinates (18), it follows that $0 \le \bar{x}(\alpha) \le 0.9 \le 1.57 = 0.5\pi$ (see Figure 1) can be interpreted as the maximal allowable stationary phase angle difference (mod $2\pi$). Likewise, $\omega_{\Delta_{\max}}$ represents the worst case difference in the oscillator's natural frequencies. Thus, since the sine function is locally proportional to its argument, (36) implies that the minimum network coupling strength $\lambda_{\min}(A)$ has to scale inversely proportional to the maximal allowable stationary phase angle difference $\bar{x}(\alpha)$ (mod $2\pi$). That is, in order to ensure synchronization to a solution with smaller stationary phase angle differences, larger coupling strengths are required. In addition, $\lambda_{\min}(A)$ has to dominate over the heterogeneity in natural frequencies $\omega_{\Delta_{\max}}$. Moreover, 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)$.

Condition (38) depends on the network topology, encompassed in the Laplacian matrix $\mathcal{L}$ in (35) and restricts

the heterogeneity of the network topology in terms of the distance of the eigenvalues of $\mathcal{L}$. These topological properties restrict the magnitude of $\alpha$ and, thus, the maximal allowable stationary phase angle difference, see Figure 1.

By comparing (37) with (23), we see that (36) ensures the existence of equilibria. In particular, by choosing $\alpha = 0$ and considering $\bar{x}(\alpha = 0) = 0.5667$, we find that (38) is always satisfied and (37) represents a constrained version of (23). This fact implies that all stationary solutions of the system (26), whose phase angle differences are less than

$$ |\theta_i^* - \theta_k^*| = |\eta_\ell^*| < 0.5667 < \frac{1}{5.5437}\pi \quad \forall i \in \mathcal{N}, \forall k \in \mathcal{N}_i, $$

are almost globally asymptotically stable. The range of admissible stationary phase angles can be increased by maximizing $\alpha$ subject to (31). Thereby, the maximum admissible value for $|\eta_\ell^*|$ with our conditions is $\bar{x}(\alpha) = 0.9 \approx \pi/3.5$, see Figure 1.

## 5. SIMULATION EXAMPLE

To illustrate our findings, we consider the example presented in (Gushchin et al., 2015b, Example 1). In this counterexample, three oscillators in an acyclic graph with $\epsilon$ being a small positive constant and the following parameters are considered:

$$ \mathcal{B} = \begin{bmatrix} 1 & 1 \\ -1 & 0 \\ 0 & -1 \end{bmatrix}, \quad A = \mathbf{I}_3, \quad \omega = \begin{bmatrix} 2 - \epsilon \\ -1 + 0.5\epsilon \\ -1 + 0.5\epsilon \end{bmatrix}, \quad \theta(0) = \begin{bmatrix} 0 \\ \frac{\pi}{2} \\ -\frac{\pi}{2} \end{bmatrix}, $$

The condition (23) results in $|1 - 0.5\epsilon| \le 1$. Therefore, for any $\epsilon \in [0, 4]$, there exist equilibrium points for the system (21). However, to ensure almost global synchronization, we need (32), which is, in general, more restrictive than (23). At first, we calculate the maximum possible $\alpha$, such that the condition (31) holds for all $\phi \in S^m$. For the current example with a precision of $10^{-2}$, this condition results in $\alpha < 0.86$. Choosing $\alpha = 0.85$, we obtain $\bar{x}(\alpha) = 0.79$. We recall that (31) and, thus, the upper bound for $\alpha$ merely depends on the graph topology. For the current graph, the obtained upper bound for $\alpha$ is very close to the maximum admissible value with our conditions, i.e., $\bar{x}(\alpha) = 0.9 \approx \pi/3.5$, see Figure 1. Thus, considering $\epsilon = 0.1$, the condition (32) is satisfied if $\lambda_{\min}(A) > 1.348$. But in (Gushchin et al., 2015b, Example 1), the coupling strength matrix is set to $A = \mathbf{I}_3$, which implies $\lambda_{\min}(A) = 1$. Hence, it is very far from our condition (32).

On the contrary, by considering the matrix $A = 1.35\mathbf{I}_3$ and $\epsilon = 0.1$, condition (32) holds. Consequently, the oscillators will globally synchronize from almost any initial condition. To illustrate this, by choosing the same equilibrium point, i.e., $\theta(0) = [0, \pi/2, -\pi/2]^\top$, we obtain the behavior presented in Figure 2. In Figure 2a, we can see that the frequencies converge to zero. Meanwhile, the trajectories $\theta$ are displayed in Figure 2b. We see that the components of $\theta$ converge to zero, which implies that all oscillators tend to the stable equilibrium point $\theta^* = [0.5204, -0.2602, -0.2602]^\top$, i.e., $\eta^* = [0.7806, 0.7806]^\top$ which can be calculated from (24). The infinity norm of this equilibrium point is equal to $\|\eta^*\|_{\infty} = 0.7806$ (mod $2\pi$), which according to (37) is less than the maximum allowable distance for $A = 1.35\mathbf{I}_3$.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Simulation example</figcaption></figure>

## 6. CONCLUSIONS

We have developed a global synchronization analysis for an acyclic network of heterogeneous Kuramoto oscillators. At first, we have determined the necessary and sufficient conditions for the existence of a synchronized motion. Subsequently, assuming knowledge of the local stability properties, we have used a change of coordinates and presented sufficient conditions under which all oscillators synchronize from all initial conditions except a set of Lebesgue measure zero. In contrast to previous works, the conditions presented in this paper do not depend on the magnitudes of initial phase angles.

The main result has been proven by using the multivariable Leonov function approach. In particular, we provide two lemmata to ensure that all conditions in Corollary 1 are satisfied. The application of the proposed conditions is illustrated via a numerical example obtained by redesigning the counterexample presented in (Gushchin et al., 2015b). The derived conditions can be associated with the main features of the network properties. 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. This parameter $\alpha$ and the coupling strengths constrain the maximum admissible difference between natural frequencies.

Current and future work is geared towards extending the analysis to generic, meshed network topologies as well as to other classes of oscillators.

## REFERENCES

* Acebrón, J.A. et al. (2005). The Kuramoto model: A simple paradigm for synchronization phenomena. *Rev. Mod. Phys.*, 77.

* Arcak, M. (2007). Passivity as a design tool for group coordination. *IEEE Transactions on Automatic Control*.

* Dörfler, F. and Bullo, F. (2011). On the critical coupling for Kuramoto oscillators. *SIAM Journal on Applied Dynamical Systems*, 10(3).

* Dörfler, F. and Bullo, F. (2014). Synchronization in complex networks of phase oscillators: A survey. *Automatica*.

* Dörfler, F., Chertkov, M., and Bullo, F. (2013). Synchronization in complex oscillator networks and smart grids. *Proceedings of the National Academy of Sciences*.

* Efimov, D. and Schiffer, J. (2019). On boundedness of solutions of state periodic systems: A multivariable cell structure approach. *IEEE Transactions on Automatic Control*, 64(10).

* Efimov, D. et al. (2017). A relaxed characterization of ISS for periodic systems with multiple invariant sets. *European Journal of Control*, 37.

* Godsil, C. and Royle, G.F. (2001). *Algebraic graph theory*, volume 207. Springer Science & Business Media.

* Gushchin, A., Mallada, E., and Tang, A. (2015a). Synchronization of heterogeneous Kuramoto oscillators with arbitrary topology. In *2015 American Control Conference*.

* Gushchin, A., Mallada, E., and Tang, A. (2015b). Synchronization of heterogeneous Kuramoto oscillators with graphs of diameter two. In *2015 49th Annual Conference on Information Sciences and Systems (CISS)*.

* Jadbabaie, A. and others. (2004). On the stability of the Kuramoto model of coupled nonlinear oscillators. In *American Control Conference (ACC)*, volume 5.

* Kassabov, M. et al. (2022). A global synchronization theorem for oscillators on a random graph. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 32(9).

* Kuramoto, Y. (1975). Self-entrainment of a population of coupled non-linear oscillators. In H. Araki (ed.), *International Symposium on Mathematical Problems in Theoretical Physics*. Springer Berlin Heidelberg.

* Leonov, G.A. (1974). On the boundedness of the trajectories of phase systems. *Siberian Mathematical Journal*.

* Manik, D. et al. (2017). Cycle flows and multistability in oscillatory networks. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 27(8).

* Mendoza-Ávila, J. et al. (2022). On relaxed conditions of integral ISS for multistable periodic systems. In *61th IEEE Conference on Decision and Control (CDC)*.

* Monzon, P. and Paganini, F. (2005). Global considerations on the Kuramoto model of sinusoidally coupled oscillators. In *44th IEEE CDC*.

* Monzon, P. and Potrie, R. (2006). Local and global aspects of almost global stability. In *45th IEEE CDC*.

* Noldus, E.J. (1977). New direct Lyapunov-type method for studying synchronization problems. *Automatica*, 13(2).

* Rodrigues, F.A. et al. (2016). The Kuramoto model in complex networks. *Physics Reports*, 610.

* Schiffer, J. and Efimov, D. (2023). Strong and weak Leonov functions for global boundedness of state periodic systems. *IEEE Transactions on Automatic Control*.

* Schiffer, J., Efimov, D., and Ortega, R. (2019). Global synchronization analysis of droop-controlled microgrids—a multivariable cell structure approach. *Automatica*, 109.

* Schiffer, J. et al. (2014). Conditions for stability of droop-controlled inverter-based microgrids. *Automatica*.

* Simpson-Porco, J.W., Dörfler, F., and Bullo, F. (2013). Synchronization and power sharing for droop-controlled inverters in islanded microgrids. *Automatica*, 49(9).

* Trip, S., Bürger, M., and De Persis, C. (2016). An internal model approach to (optimal) frequency regulation in power grids with time-varying voltages. *Automatica*, 64.

* van der Schaft, A. (2000). *L2-Gain and passivity techniques in nonlinear control*. Springer, Netherlands.

* Weiss, G., Dörfler, F., and Lebron, Y. (2019). A stability theorem for networks containing synchronous generators. *Systems & Control Letters*, 134.

* Zhu, L. and Hill, D.J. (2020). Synchronization of Kuramoto oscillators: A regional stability framework. *IEEE Transactions on Automatic Control*, 65(12).