Construction of Lyapunov Certificates for Systems on Torus Using Trigonometric Polynomials

1. Introduction

Current stability methods are still improving for dynamical systems such as humanoid robots and unmanned aerial vehicles, which are becoming increasingly prominent in today’s technology. The techniques that have been developed are available for very special classes of systems or are based on heuristic methods. With the increased computer processing power, algorithmic optimization methods have recently been developed to stabilize dynamic systems.

Lyapunov stability theorems are extensively utilized in the design of stabilizing feedback for various applications, including path-following robots [1], stabilization of fixed-wing unmanned aerial vehicles (UAV) [2], circular motion stability in UAV [3], unmanned surface vehicles [4], and magnetic field-operated satellites using the sum of squares method (SOS) [5]. Using the global Lyapunov functions, a closed-loop control system can reach a desired equilibrium point under every initial condition. However, two difficulties arise here: First, the global Lyapunov method can only directly be applied to systems on an Euclidean space [6]. For an equilibrium point to be globally stable in smooth vector fields, the state space must be equivalent (homeomorphic) to an Euclidean space [7]. In contrast, state space of a rotational dynamical system is often a multidimensional torus (for systems described by angular dynamical variables only) or a multidimensional cylinder (for systems described by both angular and non-angular variables). For such systems, it is not possible to stabilize an equilibrium point globally with continuous state feedback. The dual-Lyapunov method was introduced in [8] to overcome the two difficulties encountered while applying the global Lyapunov theorem. Whereas the Lyapunov method uses a continuous function decreasing along the singular solutions as a certificate of stability, the dual-Lyapunov method uses an increasing density function (Lyapunov density) along the solution sets as a certificate of almost global stability. Another difference between the two methods is that the dual-Lyapunov theorem guarantees convergence for almost every initial condition, while the global Lyapunov stability theorem guarantees that the solutions converge to an equilibrium point for every initial condition. [6,8] state that stability certification with the dual-Lyapunov theorem may be more advantageous than the classical Lyapunov method for systems expressed in angular variables. In [9], it was shown that the dual-Lyapunov method is more suitable for design with linear matrix inequalities since it requires an inequality that depends linearly on the design quantities, unlike the conventional Lyapunov method.

Since in rotational systems, the non-angular variables (e.g., angular velocity) depend on the angular variables (e.g., angle), fixing the angular variables at a value means convergence of the non-angular variables. More precisely, since angular state variables define a bounded manifold such as a torus, and since every continuous function (e.g., polynomials) on a bounded manifold is integrable, it may be possible to propose a systematic dual-Lyapunov design method for angular state variables in a polynomial approach. In this study, both the conventional local Lyapunov functions and the (global) Lyapunov density have been used to certify local stability and almost global stability, respectively, of an equilibrium or of a limit cycle. The polynomial approximation method based on the sum of squares (SOS) and semidefinite programming (SDP) [10], widespread in the literature in the last 20 years, the trigonometric polynomials [11], which has primary applications in signal processing, is used. SDP formulations of positive trigonometric polynomials have been utilized in signal processing applications such as FIR (finite impulse response) filters [12], beamformers [13], adapted wavelets [14], compaction filters [15] etc. In control engineering, trigonometric polynomials have been used in stability analysis of discrete-time systems [16].

In this work, we present the certifications of the almost global stability of dynamical systems on hypertorus in Theorem 1, in Section 3 by using Gram matrices and trace parametrizations of trigonometric polynomials. We use this theory to certify some physical phenomenons of coupled oscillators such as the almost global phase synchronization of uniform all-to-all coupled three oscillators in Example 1, the almost global phase-locking in three nonuniform coupled oscillators in Example 3, the almost global partial-phase-cohesion in nonuniform Kuramoto model with a tree-topology in Example 4. Using Gram matrices of trigonometric polynomials, we improve a technique for the certification of local Lyapunov stability of systems on hypertorus by using Putinar’s Positivstellensatz Lemma 2 in Theorem 2 in Section 3. We employ Theorem 2 to certify the local Phase Synchronization of uniform all-to-all coupled Oscillators for three and four oscillators in Example 2.

2. Preliminaries

2.1. Trigonometric Polynomials

A trigonometric polynomial $R\left(\mathbf{\theta }\right):{\mathbb{T}}^d={[0,2\mathbf{\pi })}^d\to \mathbb{R}$ is defined by

$$R\left(\mathbf{\theta }\right)=\sum _{\mathbf{k}=-{\mathbf{n}}_{\mathbf{r}}}^{{\mathbf{n}}_{\mathbf{r}}}{r}_{\mathbf{k}}{\mathrm{e}}^{i〈\mathbf{k},\mathbf{\theta }〉}=\sum _{k_1=-n_r\left(1\right)}^{n_r\left(1\right)}\dots \sum _{k_d=-n_r\left(d\right)}^{n_r\left(d\right)}{r}_{\mathbf{k}}{\mathrm{e}}^{i〈\mathbf{k},\mathbf{\theta }〉},$$

(1)

where ${\mathbf{n}}_{\mathbf{r}}=\left(n_r\left(1\right),n_r\left(2\right)\dots ,n_r\left(d\right)\right)\in {\mathbb{Z}}^d$ is the degree of the trigonometric polynomial, $\mathbf{k}=(k_1,k_2,\cdots ,k_d)\in {\mathbb{Z}}^d$; $\mathbf{\theta }=({\mathbf{\theta }}_1,{\mathbf{\theta }}_2,\cdots ,{\mathbf{\theta }}_d)\in {\mathbb{T}}^d$; and $r_{\mathbf{k}}$ are complex coefficients of $R\left(\mathbf{\theta }\right)$. Observe the difference in how the elements of ${\mathbb{Z}}^d$ are denoted depending on whether they represent degrees or indices. We make this distinction to make tracking of indices throughout the paper easier. Here, $〈\mathbf{k},\mathbf{\theta }〉={\sum }_{j=1}^d{k}_j{\theta }_j$ is an inner product. For real-valued a trigonometric polynomial $R\left(\mathbf{\theta }\right)$, the coefficients $r_{\mathbf{k}}$ must satisfy the symmetry (Hermitian) condition $r_{-\mathbf{k}}=r_{\mathbf{k}}^{*}$ where $x^{*}$ denotes the complex conjugate of x,11]. The inequality $-{\mathbf{n}}_{\mathbf{r}}\le \mathbf{k}\le {\mathbf{n}}_{\mathbf{r}}$ is understood elementwise and is equivalently denoted by $\left|\mathbf{k}\right|\le {\mathbf{n}}_{\mathbf{r}}$ and $\mathbf{k}\in [-{\mathbf{n}}_{\mathbf{r}},{\mathbf{n}}_{\mathbf{r}}]$. We denote the vectors $(0,\dots ,0),(1,\dots ,1)\in {\mathbb{Z}}^d$ by 0 and 1 respectively; and $max\{\mathbf{k},\mathbf{j}\}=\left(max\{k_1,j_1\},\dots ,max\{k_d,j_d\}\right)$.

The Gram matrix $G_R$ associated with $R\left(\mathbf{\theta }\right)$ is defined by

$$R\left(\mathbf{\theta }\right)={\left(z{\left(\mathbf{\theta }\right)}^{*}\right)}^T{G}_{R}z\left(\mathbf{\theta }\right),$$

(2)

where $G_R$ is a square matrix of size ${\prod }_{i=1}^d\left(n_r\left(i\right)+1\right)$, $z\left(\mathbf{\theta }\right)$ is a vector composed of monomials from the trigonometric basis on ${\mathbb{T}}^d$, and $z{\left(\mathbf{\theta }\right)}^{*}$ is its conjugate transpose. In particular, $z\left(\mathbf{\theta }\right)={\left[1,{\mathrm{e}}^{i{\mathbf{\theta }}_d},\cdots ,{\mathrm{e}}^{i{n}_r\left(d\right){\mathbf{\theta }}_d}\right]}^{\top }\otimes \cdots \otimes {\left[1,{\mathrm{e}}^{i{\mathbf{\theta }}_1},\cdots ,{\mathrm{e}}^{i{n}_r\left(1\right){\mathbf{\theta }}_1}\right]}^{\top }$.

If the Gram matrix $G_R$ is positive semidefinite, then $R\left(\mathbf{\theta }\right)\ge 0$,17]. Even more can be said in this case. Any positive semidefinite matrix can be decomposed as $U^{*}U$ for some upper triangular matrix U, known as its Cholesky decomposition. This implies that $R\left(\mathbf{\theta }\right)={\parallel Uz\left(\mathbf{\theta }\right)\parallel }^2$ for some upper triangular matrix U, hence can be expressed as a sum of squares (SOS) which are defined to have a general form

$$R\left(\mathbf{\theta }\right)=\sum _{\ell =1}^m{\left|P_{\ell }\left(\mathbf{\theta }\right)\right|}^2,$$

(3)

where each $P_{\ell }$ is a trigonometric polynomial.

Toeplitz matrices occur frequently in the study of trigonometric polynomials. An $n\times n$ matrix T is said to be the $k^{th}$ elementary Toeplitz matrix of size n if it is a $(0,1)-$matrix with ones on the $k^{th}$ diagonal, meaning $T_{(i,j)}=1$ if and only if $j-i=k$. We will denote the $k^{th}$ elementary Toeplitz matrix of size n by $T_k^n$. Clearly, $T_k^n$ is a zero matrix for any $k\ge n$. Define $T_{\mathbf{k}}^{\mathbf{n}}=T_{k_d}^{n\left(d\right)}\otimes T_{k_{d-1}}^{n(d-1)}\otimes \cdots \otimes T_{k_1}^{n\left(1\right)}$ where ⊗ denotes the Kronecker product. If $\mathbf{k}\notin (-\mathbf{n},\mathbf{n})$, then $k_i\ge n\left(i\right)$ for some i. In that case $T_{k_i}^{n\left(i\right)}$, hence $T_{\mathbf{k}}^{\mathbf{n}}$, is the zero matrix. Toeplitz matrices are utilized to relate the coefficients $r_{\mathbf{k}}$ of $R\left(\mathbf{\theta }\right)$ to the Gram matrix $G_R$ in (2) by a property known in the literature as trace parameterization [11], which states for any $\mathbf{k}\in [-{\mathbf{n}}_r,{\mathbf{n}}_r]$

$$r_{\mathbf{k}}=\mathrm{Tr}[T_{\mathbf{k}}^{\mathbf{n}}\mathbf{r}+\mathbf{1}{G}_R].$$

(4)

Positive semidefinite matrices can parametrize SOS trigonometric polynomials [11,18]. Equivalently, $R\left(\mathbf{\theta }\right)$ is SOS if and only if there exits positive semidefinite Gram matrix $G_R$ in (2) satisfying (4) (Theorem 3.15, [11]).

Given the trigonometric polynomial $R\left(\mathbf{\theta }\right)$ defined in (1), we can define the derivative with respect to ${\theta }_l$ as follows

$$D_{l}R\left(\mathbf{\theta }\right):={\displaystyle \frac{\mathrm{d}R}{\mathrm{d}{\mathbf{\theta }}_l}}\left(\mathbf{\theta }\right)=\sum _{\mathbf{k}=-{\mathbf{n}}_{\mathbf{r}}}^{{\mathbf{n}}_{\mathbf{r}}}d{r}_{\mathbf{k}}{\mathrm{e}}^{i〈\mathbf{k},\mathbf{\theta }〉}$$

(5)

where $d{r}_{\mathbf{k}}=i\times k_l\times r_{\mathbf{k}}$. Thus, due to (4) the relationship between Gram matrix $G_R$ of $R\left(\mathbf{\theta }\right)$ and Gram matrix $G_{D_{l}R}$ of $D_{l}R\left(\mathbf{\theta }\right)$ is given

$$\mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{r}}+\mathbf{1}}{G}_{D_{l}R}\right]=i{k}_l\mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{r}}+\mathbf{1}}{G}_R\right],$$

(6)

where $T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{r}}+\mathbf{1}}$ are the Toeplitz matrices corresponding to $z\left(\theta \right)$ in (2).

Consider another trigonometric polynomial $P\left(\mathbf{\theta }\right)$ defined as

$$P\left(\mathbf{\theta }\right)=\sum _{\mathbf{k}=-{\mathbf{n}}_{\mathbf{p}}}^{{\mathbf{n}}_{\mathbf{p}}}{p}_{\mathbf{k}}{\mathrm{e}}^{i〈\mathbf{k},\mathbf{\theta }〉},$$

(7)

where $p_{\mathbf{k}}$ are the coefficients of $P\left(\mathbf{\theta }\right)$. The product of the polynomials $P\left(\mathbf{\theta }\right)$ and $R\left(\mathbf{\theta }\right)$, denoted as $PR\left(\mathbf{\theta }\right)$, is also a trigonometric polynomial which can expanded as

$$PR\left(\mathbf{\theta }\right)=\sum _{\mathbf{k}=-{\mathbf{n}}_{\mathbf{pr}}}^{{\mathbf{n}}_{\mathbf{pr}}}p{r}_{\mathbf{k}}{\mathrm{e}}^{i〈\mathbf{k},\mathbf{\theta }〉},$$

(8)

where ${\mathbf{n}}_{\mathbf{pr}}={\mathbf{n}}_{\mathbf{p}}+{\mathbf{n}}_{\mathbf{r}}$, and the coefficients $p{r}_{\mathbf{k}}$ are given by

$$p{r}_{\mathbf{k}}=\sum _{\mathbf{j}=-{\mathbf{n}}_{\mathbf{pr}}}^{{\mathbf{n}}_{\mathbf{pr}}}{p}_{\mathbf{k}-\mathbf{j}}{r}_{\mathbf{j}}.$$

(9)

The Gram matrix $G_{PR}$ associated with the product $PR\left(\mathbf{\theta }\right)$ can be expressed in terms of Gram matrices $G_P$ and $G_R$ of $P\left(\mathbf{\theta }\right)$ and $R\left(\mathbf{\theta }\right)$, respectively. Specifically, we have

$$\mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{pr}}+\mathbf{1}}{G}_{PR}\right]=\sum _{\mathbf{j}=-{\mathbf{n}}_{\mathbf{pr}}}^{{\mathbf{n}}_{\mathbf{pr}}}\mathrm{Tr}\left[T_{\mathbf{k}-\mathbf{j}}^{{\mathbf{n}}_{\mathbf{p}}+\mathbf{1}}{G}_P\right]\times \mathrm{Tr}\left[T_{\mathbf{j}}^{{\mathbf{n}}_{\mathbf{r}}+\mathbf{1}}{G}_R\right].$$

(10)

We end this section by defining a zero-sum matrix or a matrix having zero-sum. As the name suggests, all elements of such a matrix sum up to 0. If the Gram-matrix $G_R$ in (2) has zero-sum, then it is clear that $R\left(0\right)=1G_R{1}^{\top }=0$.

2.2. A Class of Coupled Oscillators

In this article, we consider a class of coupled oscillators defined on the hypertorus ${\mathbb{T}}^d$ as

$$\dot{\mathbf{\theta }}=F\left(\mathbf{\theta }\right)={\left(F_1\left(\mathbf{\theta }\right),\dots ,F_d\left(\mathbf{\theta }\right)\right)}^{\top },$$

(11)

where each component function $F_l:{\mathbb{T}}^d\to \mathbb{R}$ has a Fourier series expansion

$$F_l\left(\mathbf{\theta }\right)=\sum _{\mathbf{k}=-{\mathbf{n}}_{\mathbf{f}}}^{{\mathbf{n}}_{\mathbf{f}}}{f}_{\mathbf{k}}^{\left(l\right)}{e}^{i〈\mathbf{k},\mathbf{\theta }〉},$$

(12)

with Fourier coefficients $f_{\mathbf{k}}^{\left(l\right)}$. Without loss of generality, we can assume that all component functions can be expressed using the same number of harmonics. This can be done by taking ${\mathbf{n}}_{\mathbf{f}}={max}_{\ell =1}^d{n}_f(\ell )\times \mathbf{1}$ and adding zero fourier coefficients wherever necessary. Every smooth vector field on a compact manifold is complete [19], meaning that the flow ${\Phi }_t:{\mathbb{T}}^d\to {\mathbb{T}}^d$ is defined for each $t\in \mathbb{R}$. Thus, the solutions of (11) exist for all time, irrespective of the initial condition. A set $A\subseteq {\mathbb{T}}^d$ is said to be forward-invariant (backward-invariant, respectively) under the flow of (11) if ${\Phi }_t\left(A\right)\subseteq A$ for all $t\ge 0$ ($t\le 0$, respectively). A is forward-invariant if and only if $A^c$ is backward-invariant.

Global (asymptotic) stability has been a well-explored issue in systems theory for the past several decades and is characterized by all solution curves converging to the equilibrium point. Due to topological restrictions, systems on a hypertorus cannot be globally asymptotically stable. However, in some cases, the system (11) might exhibit

• almost global stability: if the set of divergent solution curves (with respect to an equilibrium point) has zero measure; and
• local stability: if there exists a region $\mathfrak{R}⊊{\mathbb{T}}^d$ such that all solution curves originating from $\mathfrak{R}$ converge to an equilibrium point in $\mathfrak{R}$.

Another interesting notion for coupled oscillators is that of synchronization. We say that the oscillators i and j are phase synchronized if ${lim}_{t\to \infty }\left({\theta }_i\left(t\right)-{\theta }_j\left(t\right)\right)=0$; phase-locked if ${lim}_{t\to \infty }\left({\theta }_i\left(t\right)-{\theta }_j\left(t\right)\right)={\theta }_{ij}$ for some ${\theta }_{ij}\in \mathbb{T}$; and frequency synchronized if ${lim}_{t\to \infty }\left({\displaystyle \frac{\mathrm{d}{\theta }_{\mathrm{i}}\left(\mathrm{t}\right)}{\mathrm{d}\mathrm{t}}}-{\displaystyle \frac{\mathrm{d}{\theta }_{\mathrm{j}}\left(\mathrm{t}\right)}{\mathrm{d}\mathrm{t}}}\right)=0$. The system (11) is said to exhibit

• almost global phase synchronization (phase-locking, frequency synchronization, respectively): if the set of initial conditions on ${\mathbb{T}}^d$ for which a pair of oscillators is not phase synchronized (not phase-locked, not frequency synchronized, respectively) has zero measure; and
• local phase synchronization (phase-locking, frequency synchronization, respectively): if there exists a region $\mathfrak{R}⊊{\mathbb{T}}^d$ such that the set of initial conditions on ${\mathbb{T}}^d$ for which a pair of oscillators is not phase synchronized (not phase-locked, not frequency synchronized, respectively) is a subset of ${\mathfrak{R}}^c$.

Under suitable conditions, for example, in the presence of phase-shift symmetry $F(\mathbf{\theta }+\mathbf{ϵ}·1)=F\left(\mathbf{\theta }\right)$, the system (11) on ${\mathbb{T}}^d$ can be reduced to a phase-difference system on ${\mathbb{T}}^{d-1}$, refer [20]. In that case, the system (11) exhibits almost global phase synchronization (phase-locking, respectively) if and only if its corresponding phase-difference system exhibits almost global stability to the origin (to an equilibrium point ${\Theta }_0\ne 0$, respectively). This relationship has also been briefly touched upon in Example 2 for the local counterpart of the statement.

In this article, we are interested in obtaining a certificate for almost global stability and a certificate for local stability of the sysem (11). Lyapunov method, which relies on finding an energy-like function that decreases along solution trajectories of the system, has been a tool to establish the global stability of systems. We will employ the Lyapunov method locally to establish the local stability of (11). The Lyapunov method has a sibling, the dual-Lyapunov method proposed by Rantzer (Theorem 1, [8]), which is a well-known tool to certify the almost global stability of such systems. Unlike the Lyapunov method, the dual-Lyapunov method relies on finding a `Lyapunov density’ satisfying certain density propagation inequalities. The following lemma can be seen as a modified version of Rantzer’s theorem for the system (11).

(Theorem 4.2, [21])).Lemma 1 (Adapted from  Let $A\subseteq {\mathbb{T}}^d$ be a compact set such that $A^c$ is backward-invariant under the flow of (11). If there exists a continuously differentiable function $\rho :A^c\to \mathbb{R}$ and satisfies

C1.

$\rho \left(\mathbf{\theta }\right)>0$ for almost every $\mathbf{\theta }\in A^c$.

C2.

ρ is integrable away from A, that is, if $A_{ϵ}=\{\mathbf{\theta }\in {\mathbb{T}}^d:{inf}_{\mathbf{\psi }\in A}dist(\mathbf{\theta },\mathbf{\psi })<\mathbf{ϵ}\}$

$${\int }_{{\left(A_{ϵ}\right)}^c}\rho \left(\mathbf{\theta }\right)\mathbf{\mu }\left(d\mathbf{\theta }\right)<\infty \forall \mathbf{ϵ}>0.$$

C3.

$\nabla ·\left(\rho F\right)\left(\mathbf{\theta }\right)>0$ for almost every $\mathbf{\theta }\in A^c$.

Then almost all solutions of (11) converge to A as $t\to \infty$.

Note that the flow of (11) is defined for all time; hence, the convergence of solution curves can be discussed. The second condition is a constraint on the evolution of density $\rho$, whose details can be found in [8]. It is known that a density $\rho$ satisfying the evolution constraint inevitably has a singularity in A (Lemma 3, [21]). Considering this and taking $\rho \left(\mathbf{\theta }\right)=V{\left(\mathbf{\theta }\right)}^{-1}$ where $V\left(\mathbf{\theta }\right)$ is a trigonometric polynomial, we obtain the following corollary.

Corollary 1.

Assume there exist nonzero trigonometric polynomials $V\left(\mathbf{\theta }\right)\ge 0$ and $W\left(\mathbf{\theta }\right)\ge 0$ such that $V^{-1}\left(0\right)$ is forward-invariant under the flow of (11); and for almost every $\mathbf{\theta }\in {\mathbb{T}}^d$

$$W\left(\mathbf{\theta }\right)=-\nabla V\left(\mathbf{\theta }\right)·F\left(\mathbf{\theta }\right)+V\left(\mathbf{\theta }\right)\times (\nabla ·F)\left(\mathbf{\theta }\right),$$

(13)

then almost all solutions of (11) converge to $V^{-1}\left(0\right)$ as $t\to \infty$. Thus $\rho \left(\mathbf{\theta }\right)=1/V\left(\mathbf{\theta }\right)$ is a Lyapunov density for (11).

Proof. 

Consider $\rho \left(\mathbf{\theta }\right)=1/V\left(\mathbf{\theta }\right)$ where $V\left(\mathbf{\theta }\right)\ge 0$ is a nonzero trigonometric polynomial and take $A=V^{-1}\left(0\right)$. A is the set of all singularities of $\rho \left(\mathbf{\theta }\right)$, is a closed subset of a compact space (hence compact), and has zero measure. Thus, A is a forward-invariant compact subset under the flow of (11). This implies $A^c$ is a backward-invariant subset under the flow of (11). On substitution, condition C3 reduces to (13) for some nonzero trigonometric polynomial $W\left(\mathbf{\theta }\right)\ge 0$. The set of singularities of $\rho$ is nonempty due to Remark 3 [21], thus A is nonempty. Also, $\rho$ is positive on $A^c$ and bounded, hence integrable, over ${\left(A_{ϵ}\right)}^c$ for any $ϵ$. Thus, conditions C1 and C2 is satisfied. Hence, by Lemma 1, almost all solutions of (11) converge to the set $A=V^{-1}\left(0\right)$. □

2.2.1. Generalized Kuramoto Models with Coupling Topology

Consider the diffusively coupled system of d phase oscillators given by

$${\dot{\phi }}_j={\omega }_j+K\sum _{k=1}^d{a}_{kj}g({\phi }_k-{\phi }_j),{\phi }_i\in \mathbb{T},i=1,\cdots ,d,$$

(14)

where ${\omega }_i$ is the natural frequency of the oscillator ${\phi }_i$ ; $A=\left(a_{kj}\right)$ is a $(0,1)-$matrix defining the coupling topology, meaning $a_{kj}=1$ if and only if there is an edge from the oscillator ${\phi }_k$ to ${\phi }_j$ (in physical terms, ${\phi }_k$ influences ${\phi }_j$); K is the coupling constant; and g is a coupling function expressed in terms of its first L harmonics

$$g\left(x\right)=\sum _{k=1}^L{\alpha }_{k}sin(kx+{\beta }_k),$$

(15)

where $\left|{\beta }_k\right|<\pi /2$ for all k. The system is called all-to-all coupled if $a_{kj}=1$ for all $k,j\in \{1,\dots ,d\}$. System of coupled oscillators where all oscillators have identical natural frequencies ${\omega }_i={\omega }_j$ for all $i,j\in \{1,\dots ,d\}$ are called uniform coupled oscillators. It is easy to observe that all oscillators should necessarily share a common natural frequency for phase synchronization for an all-to-all connected system. Note that (14) has phase shift symmetry $F(\mathbf{\theta }+\mathbf{ϵ}·1)=F\left(\mathbf{\theta }\right)$; hence we reiterate that local (or almost global) synchronization can be established by proving local (or almost global) stability of its corresponding phase-difference system, refer [20] and Example 2. When $g\left(x\right)=sin\left(x\right)$, all-to-all coupled oscillators take the form of the well-known Kuramoto model

$${\dot{\phi }}_j={\omega }_j+K\sum _{k=1}^{d}sin({\phi }_k-{\phi }_j),{\phi }_i\in \mathbb{T},i=1,\cdots ,d,$$

(16)

hence we refer to phase oscillator (14) as a generalized Kuramoto model with coupling topology.

When oscillators are nonuniform, that is ${\omega }_i\ne {\omega }_j$ for some pair $i,j\in \{1,\dots ,d\}$ in (14), phase synchronization is not possible. It is known that nonuniform oscillators (14) exhibit frequency synchronization for sufficiently large coupling K. Frequency synchronization is necessary for phase-locking but not sufficient. However, for the Kuramoto model (16), both are equivalent Remark 2.1 [22].

Having discussed all the prerequisites required for the paper, we present the main results in the next section.

3. Trigonometric Polynomials in Stability Analysis

In this section, we will obtain a certificate for almost global stability and a certificate for local stability of the sysem (11). The main idea is to find a Lyapunov density of the form $\rho \left(\mathbf{\theta }\right)=1/V\left(\mathbf{\theta }\right)$ in the former case, and a local Lyapunov function $V\left(\mathbf{\theta }\right)$ in the latter case, for some nonnegative trigonometric polynomial V. We will observe that the existence of a Lyapunov density (function, respectively) is established if the coefficients of V are zeros of some multivariate polynomial. The main problem in such a constraint is that obtaining zeroes of multivariate polynomials is extremely hard. However, trace parameterization turns the expression into a Semi-Definite Programming (SDP) problem involving matrices, the feasibility of which can be checked using solvers such as Sedumi in MATLAB. We now present our first result.

Theorem 1.

Given a dynamical system $\dot{\theta }=F\left(\theta \right)$ where each component function $F_l:{\mathbb{T}}^d\to \mathbb{R}$ has Fourier series expansion (12), if for a nonzero Hermitian matrix $G_V⪰0$ of size ${\prod }_{i=1}^d\left(n_v\left(i\right)+1\right)$ having zero-sum, there exist a nonzero Hermitian matrix $G_W⪰0$ of size ${\prod }_{i=1}^d\left(n_f\left(i\right)+n_v\left(i\right)+1\right)$ such that

$$\mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_W\right]=\sum _{l=1}^d\sum _{\mathbf{j}=-({\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}})}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}}i{j}_l\left(-\mathrm{Tr}\left[T_{\mathbf{j}}^{{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_V\right]\times f_{\mathbf{k}-\mathbf{j}}^{\left(l\right)}+f_{\mathbf{j}}^{\left(l\right)}\times \mathrm{Tr}\left[T_{\mathbf{k}-\mathbf{j}}^{{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_V\right]\right),$$

(17)

for all $\left|\mathbf{k}\right|\le {\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}$, then the dynamical system satisfies Equation (13) for $V\left(\mathbf{\theta }\right)$ and $W\left(\mathbf{\theta }\right)$ corresponding to $G_V$ and $G_W$ as per (2), respectively. Subsequently, almost all solutions of (11) converge to the origin as $t\to \infty$.

Proof. 

The system is almost globally stable if there exist positive definite trigonometric polynomials V and W such that $V\left(\mathbf{0}\right)=0$ and (13) is satisfied. The condition $V\left(\mathbf{0}\right)=0$ is equivalent to the sum of all entries of the Gram matrix $G_V$ being 0. Expanding the expression of gradient of V and divergence of F appearing in (13), we obtain that

$$\begin{array}{ccc}\hfill W\left(\mathbf{\theta }\right)& =& \sum _{l=1}^d\left(-F_l\left(\mathbf{\theta }\right)\times {\displaystyle \frac{\partial V\left(\mathbf{\theta }\right)}{\partial {\mathbf{\theta }}_l}}+V\left(\mathbf{\theta }\right)\times {\displaystyle \frac{\partial F_l\left(\mathbf{\theta }\right)}{\partial {\mathbf{\theta }}_l}}\right),\hfill \end{array}$$

(18)

must be satisfied. If $G_V$, $G_W$ and $G_{F_l}$ denote the Gram matrices of $V\left(\mathbf{\theta }\right)$, $W\left(\mathbf{\theta }\right)$ and $F_l\left(\mathbf{\theta }\right)$ respectively, the equation can be rewritten in terms of Gram matrices as

$$\begin{array}{ccc}\hfill G_W& =& \sum _{l=1}^d\left(-G_{F_l\times D_{l}V}+G_{V\times D_l{F}_l}\right).\hfill \end{array}$$

(19)

Thus, almost global stability is established if there exist positive definite matrices $G_V$ and $G_W$ satisfying the above equality and ${\sum }_{j,k}{\left(G_V\right)}_{(j,k)}=0$. Clearly, the size of the matrix $G_V$ determines the size of $G_W$. Assuming that the size of $G_V$ is ${\prod }_{i=1}^d\left(n_v\left(i\right)+1\right)$ and the size of each $G_{F_l}$ is ${\prod }_{i=1}^d\left(n_f\left(i\right)+1\right)$, the sizes of both $G_{D_{l}V\times F_l}$ and $G_{D_l{F}_l\times V}$ is ${\prod }_{i=1}^d\left(n_f\left(i\right)+n_v\left(i\right)+1\right)$. Thus, $G_W$ must be of the size ${\prod }_{i=1}^d\left(n_f\left(i\right)+n_v\left(i\right)+1\right)$. Thus for $\mathbf{k}\in [-({\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}),{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}]$, by virtue of trace parametrization (4), the ${\mathbf{k}}^{th}$ coefficient of both sides of (18) are equal if

$$\begin{array}{ccc}\hfill \mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_W\right]& =& \sum _{l=1}^d\left(-\mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_{F_l\times D_{l}V}\right]+\mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_{V\times D_l{F}_l}\right]\right)\hfill \end{array}$$

(20)

$$\begin{array}{ccc}& =& \sum _{l=1}^d\sum _{\mathbf{j}=-({\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}})}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}}\left(-\mathrm{Tr}\left[T_{\mathbf{j}}^{{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_{D_{l}V}\right]\times \mathrm{Tr}\left[T_{\mathbf{k}-\mathbf{j}}^{{\mathbf{n}}_{\mathbf{f}}+\mathbf{1}}{G}_{F_l}\right]+\mathrm{Tr}\left[T_{\mathbf{j}}^{{\mathbf{n}}_{\mathbf{f}}+\mathbf{1}}{G}_{D_l{F}_l}\right]\times \mathrm{Tr}\left[T_{\mathbf{k}-\mathbf{j}}^{{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_V\right]\right),\hfill \end{array}$$

(21)

$$\begin{array}{ccc}& =& \sum _{l=1}^d\sum _{\mathbf{j}=-({\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}})}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}}i{j}_l\left(-\mathrm{Tr}\left[T_{\mathbf{j}}^{{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_V\right]\times \mathrm{Tr}\left[T_{\mathbf{k}-\mathbf{j}}^{{\mathbf{n}}_{\mathbf{f}}+\mathbf{1}}{G}_{F_l}\right]+\mathrm{Tr}\left[T_{\mathbf{j}}^{{\mathbf{n}}_{\mathbf{f}}+\mathbf{1}}{G}_{F_l}\right]\times \mathrm{Tr}\left[T_{\mathbf{k}-\mathbf{j}}^{{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_V\right]\right),\hfill \end{array}$$

(22)

$$\begin{array}{ccc}& =& \sum _{l=1}^d\sum _{\mathbf{j}=-({\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}})}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}}i{j}_l\left(-\mathrm{Tr}\left[T_{\mathbf{j}}^{{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_V\right]\times f_{\mathbf{k}-\mathbf{j}}^{\left(l\right)}+f_{\mathbf{j}}^{\left(l\right)}\times \mathrm{Tr}\left[T_{\mathbf{k}-\mathbf{j}}^{{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_V\right]\right).\hfill \end{array}$$

(23)

□

Remark 1.

If the condition `the sum of all entries of $G_V$ is zero’ is dropped in Theorem 1, then $V\left(\mathbf{0}\right)\ne 0$. However due to [21], the trigonometric polynomial $V\left(\mathbf{\theta }\right)=1/\mathbf{\rho }\left(\mathbf{\theta }\right)$ has a zero. Thus, $V^{-1}\left(0\right)$ is a nonempty compact set. If $V^{-1}\left(0\right)\subseteq F^{-1}\left(0\right)$, then the set $V^{-1}\left(0\right)$ is clearly invariant under the flow of (11) and almost all solutions of (11) converge to $V^{-1}\left(0\right)$ as $t\to \infty$ due to Corollary 1. In particular when $V^{-1}\left(0\right)$ is a singleton, say $V^{-1}\left(0\right)=\left\{{\mathbf{\theta }}_0\right\}$ where ${\mathbf{\theta }}_0$ is an equilibrium point of the system (11), then almost all solutions of (11) converge to ${\mathbf{\theta }}_0$ as $t\to \infty$.

Given a system, the existence of a Lyapunov function, which decreases along the system’s solution trajectories, is a certificate of global stability of the system. Thus, the system (11) would be globally stable if there exists a continuously differentiable positive semidefinite function $V\left(\mathbf{\theta }\right)\neg \equiv 0$ such that $\nabla V\left(\mathbf{\theta }\right)·F\left(\mathbf{\theta }\right)$ is negative semidefinite. Thus, to find a nonnegative trigonometric polynomial that serves as a Lyapunov function for (11), we need to find a corresponding nonnegative trigonometric polynomial $W\left(\mathbf{\theta }\right)$ satisfying

$$W\left(\mathbf{\theta }\right)=-\nabla V\left(\mathbf{\theta }\right)·F\left(\mathbf{\theta }\right).$$

(24)

The Equation (24) also has a corresponding gram matrix representation. The equation is obtained from (13) by removing the second summand. Reproducing the steps of the proof of Theorem 1, we obtain (17) without the second summand of each term inside the bracket. We thus obtain that the expression (24) is satisfied by $V\left(\mathbf{\theta }\right)$ and $W\left(\mathbf{\theta }\right)$ if their corresponding Gram matrices $G_V$ and $G_W$ satisfy

$$\mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_W\right]=-\sum _{l=1}^d\sum _{\mathbf{j}=-({\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}})}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}}i{j}_l\mathrm{Tr}\left[T_{\mathbf{j}}^{{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_V\right]\times f_{\mathbf{k}-\mathbf{j}}^{\left(l\right)},$$

(25)

for all $\mathbf{k}\in [-({\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}),{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{\mathbf{v}}]$. As stated before, the expression is irrelevant in establishing global stability as no system on a hypertorus is globally stable due to topological restrictions. This will, however, be used later to establish the local stability of systems on a hypertorus by considering the equality only on certain subsets of ${\mathbb{T}}^d$.

We now recall a celebrated result, which characterizes the trigonometric polynomials that are positive in certain domains. The domain comprises points in the state space where a given set of trigonometric polynomials is positive.

Lemma 2

(Putinar’s Positivstellensatz; adapted from [11] as per notations of this paper). Let $H_{\ell }$ be trigonometric polynomials (say, of degree ${\mathbf{n}}_{\mathbf{h}}$) with complex coefficients for $\ell =1,\dots ,\mathcal{L}$ and $\mathcal{D}$ be the domain

$$\mathcal{D}=\mathbf{\theta }\in {[-\mathbf{\pi },\mathbf{\pi }]}^d:H_{\ell }\left(\mathbf{\theta }\right)\ge 0,\ell =1,\dots ,\mathcal{L}.$$

(26)

If a trigonometric polynomial $R\left(\mathbf{\theta }\right)$ is positive on $\mathcal{D}$ (that is, $R\left(\mathbf{\theta }\right)>0$ for all $\mathbf{\theta }\in \mathcal{D}$), then there exist sum-of-squares polynomials $S_{\ell }\left(\mathbf{\theta }\right)$ for $\ell =0,\dots ,\mathcal{L}$ such that

$$\mathit{R}\left(\mathbf{\theta }\right)={\mathit{S}}_0\left(\mathbf{\theta }\right)+\sum _{\ell =1}^{\mathcal{L}}{S}_{\ell }\left(\mathbf{\theta }\right){\mathit{H}}_{\ell }\left(\mathbf{\theta }\right).$$

(27)

Without loss of generality, it may be assumed that all trigonometric polynomials $S_{\ell }$ share the same degree ${\mathbf{n}}_{\mathbf{s}}$. Also, [11] states that `the degrees of sum-of-squares polynomials may be arbitrarily high, at least theoretically.’ The lemma can be transformed in terms of Gram matrices by using trace parametrization (4) as follows:

Lemma 3

(adapted from [11] as per notations of this paper). If the trigonometric polynomial $R\left(\mathbf{\theta }\right)$ has degree ${\mathbf{n}}_{\mathbf{r}}$ and is positive on $\mathcal{D}$ defined as in (26), then there exist nonzero Hermitian matrices $G_{S_{\ell }}⪰0$ of sizes ${\prod }_{i=1}^d\left(n_s\left(i\right)+1\right)$ for $\ell =0,\dots ,\mathcal{L}$ such that for any $\left|\mathbf{k}\right|\le max\{{\mathbf{n}}_{\mathbf{s}}+{\mathbf{n}}_{\mathbf{h}},{\mathbf{n}}_{\mathbf{r}}\}$

$$r_{\mathbf{k}}=\mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{s}}+\mathbf{1}}{G}_{S_0}\right]+\sum _{\ell =1}^{\mathcal{L}}\sum _{\mathbf{j}=-({\mathbf{n}}_{\mathbf{h}}+{\mathbf{n}}_{\mathbf{s}})}^{{\mathbf{n}}_{\mathbf{h}}+{\mathbf{n}}_{\mathbf{s}}}\mathrm{Tr}[T_{\mathbf{k}}-{\mathbf{j}}_{\mathbf{s}}^{\mathbf{n}}+\mathbf{1}{G}_{S_{\ell }}]\times \mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{h}}+\mathbf{1}}{G}_{H_{\ell }}\right].$$

(28)

Remark 2.

A vector $\mathbf{\theta }\in {\mathbb{T}}^d$ is said to have diameter a, which is denoted by $\mathrm{Diam}\left(\mathbf{\theta }\right)=a$, if ${max}_{i,j\in \{1,\dots ,d\}}|{\theta }_i-{\theta }_j|=a$. The set $\mathrm{Arc}\left(a\right)=\{\mathbf{\theta }\in {\mathbb{T}}^d:\mathrm{Diag}\left(\mathbf{\theta }\right)<a\}$ which is a 2d-gonal hypercylinder with “diagonal" joining $-\pi ·\mathbf{1}$ and $\pi ·\mathbf{1}$ as its axis. Refer to Figure 1 (left) where $a=\pi /2$, $d=3$, and the region is a hexagonal cylinder. The terminology $\mathrm{Arc}\left(a\right)$ is used because all components ${\theta }_j$ of any given $\mathbf{\theta }\in \mathrm{Arc}\left(a\right)$ fall on an arc of length a in $\mathbb{T}$. In the literature [23,24], the authors study nonuniform all-to-all connected Kuramoto oscillators (16) and find conditions for frequency synchronization and local phase synchronization, respectively. They show that $\mathrm{Arc}(\pi /2)$ is forward-invariant. They subsequently find a Lyapunov function in this domain. The conditions ensure that all solutions starting in $\mathrm{Arc}(\pi /2)$ get phase synchronized.

Theorem 2.

Given a dynamical system $\dot{\theta }=F\left(\theta \right)$ satisfying $F\left(0\right)=0$,where each component function $F_l:{\mathbb{T}}^d\to \mathbb{R}$ has Fourier series expansion (12); $H_{\ell }$ be trigonometric polynomials with Gram-matrices $G_{H_{\ell }}$ for $\ell =1,\dots ,\mathcal{L}$; and

$$\mathcal{D}=\{\theta \in {[-\mathbf{\pi },\mathbf{\pi }]}^d:H_{\ell }\left(\theta \right)\ge 0,\ell =1,\dots ,\mathcal{L}\}.$$

If there exist zero-sum Hermitian matrices $G_{S_{\ell }^v}⪰0$ of sizes ${\prod }_{i=1}^d\left(n_{s^v}\left(i\right)+1\right)$; Hermitian matrices $G_{S_{\ell }^w}⪰0$ of sizes ${\prod }_{i=1}^d\left(n_f\left(i\right)+n_{s^v}\left(i\right)+1\right)$ for each $\ell \in \{1,\dots ,\mathcal{L}\}$ satisfying

$$\begin{array}{c}\hfill \mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_0}^w\right]+\sum _{\ell =1}^{\mathcal{L}}\sum _{\left|\mathbf{m}\right|\le {\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+{\mathbf{n}}_{\mathbf{h}}}\mathrm{Tr}\left[T_{\mathbf{k}-\mathbf{m}}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_{\ell }^w}\right]\times \mathrm{Tr}[T_{\mathbf{m}}^{\mathbf{n}}\mathbf{h}+\mathbf{1}{G}_{H_{\ell }}]=-\sum _{l=1}^d\sum _{\left|\mathbf{j}\right|\le {\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+{\mathbf{n}}_{\mathbf{h}}}i\times \\ \hfill j_l{f}_{\mathbf{k}-\mathbf{j}}^{\left(l\right)}\left(\mathrm{Tr}\left[T_{\mathbf{j}}^{{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_0}^v\right]+\sum _{\ell =1}^{\mathcal{L}}\sum _{\left|\mathbf{m}\right|\le {\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+{\mathbf{n}}_{\mathbf{h}}}\mathrm{Tr}\left[T_{\mathbf{j}-\mathbf{m}}^{{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_{\ell }^v}\right]\times \mathrm{Tr}[T_{\mathbf{m}}^{\mathbf{n}}\mathbf{h}+\mathbf{1}{G}_{H_{\ell }}]\right),\end{array}$$

(29)

for each $\left|\mathbf{k}\right|\le {\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+{\mathbf{n}}_{\mathbf{h}}$; and a forward-invariant region $\mathfrak{R}\subseteq \mathcal{D}$ containing the origin, then the solutions of (11) with initial conditions in $\mathfrak{R}$ converge to the origin as $t\to \infty$. Hence, the system (11) is locally stable on $\mathfrak{R}$.

Proof. 

We are interested in finding trigonometric polynomials $V\left(\theta \right)$ and $W\left(\theta \right)$ which are positive on $\mathcal{D}$ and satisfy (24) for all $\theta \in \mathcal{D}$. Note that each $H_{\ell }\left(\theta \right)$ is a trigonometric polynomial of degree ${\mathbf{n}}_{\mathbf{h}}$. Thus, by Putinar’s Positivstellensatz V and W must have the forms

$$\begin{array}{ccc}\hfill V\left(\theta \right)& =& S_0^v\left(\theta \right)+\sum _{\ell =1}^{\mathcal{L}}{H}_{\ell }\left(\theta \right)S_{\ell }^v\left(\theta \right),\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill W\left(\theta \right)& =& S_0^w\left(\theta \right)+\sum _{\ell =1}^{\mathcal{L}}{H}_{\ell }\left(\theta \right)S_{\ell }^w\left(\theta \right),\hfill \end{array}$$

(31)

and their corresponding Gram matrix representations can be obtained using Lemma 3 as

$$\begin{array}{ccc}\hfill \mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{v}}+\mathbf{1}}{G}_V\right]& =& \mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_0}^v\right]+\sum _{\ell =1}^{\mathcal{L}}\sum _{\left|\mathbf{m}\right|\le {\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+{\mathbf{n}}_{\mathbf{h}}}\mathrm{Tr}\left[T_{\mathbf{k}-\mathbf{m}}^{{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_{\ell }^v}\right]\times \mathrm{Tr}[T_{\mathbf{m}}^{\mathbf{n}}\mathbf{h}+\mathbf{1}{G}_{H_{\ell }}],\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill \mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{w}}+\mathbf{1}}{G}_W\right]& =& \mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{w}}}+\mathbf{1}}{G}_{S_0}^w\right]+\sum _{\ell =1}^{\mathcal{L}}\sum _{\left|\mathbf{m}\right|\le {\mathbf{n}}_{{\mathbf{s}}^{\mathbf{w}}}+{\mathbf{n}}_{\mathbf{h}}}\mathrm{Tr}\left[T_{\mathbf{k}-\mathbf{m}}^{{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{w}}}+\mathbf{1}}{G}_{S_{\ell }^w}\right]\times \mathrm{Tr}[T_{\mathbf{m}}^{\mathbf{n}}\mathbf{h}+\mathbf{1}{G}_{H_{\ell }}].\hfill \end{array}$$

(33)

Clearly, ${\mathbf{n}}_{\mathbf{v}}={\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+{\mathbf{n}}_{\mathbf{h}}$ and ${\mathbf{n}}_{\mathbf{w}}={\mathbf{n}}_{{\mathbf{s}}^{\mathbf{w}}}+{\mathbf{n}}_{\mathbf{h}}$. Note that $G_{S_{\ell }^v}$ are zero-sum matrices, hence $S_{\ell }^v\left(0\right)=0$, for all $\ell \in \{1,\dots ,\mathcal{L}\}$. Thus, $V\left(0\right)=0$. Also (32) occurs as the left side of (25). For the right side, substitute the first equality in Equation (25) to obtain

$${\mathrm{RHS}}_{eq.\left(\right)}=-\sum _{l=1}^d\sum _{\left|\mathbf{j}\right|\le {\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+{\mathbf{n}}_{\mathbf{h}}}i{j}_l{f}_{\mathbf{k}-\mathbf{j}}^{\left(l\right)}\left(\mathrm{Tr}\left[T_{\mathbf{j}}^{{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_0}^v\right]+\sum _{\ell =1}^{\mathcal{L}}\sum _{\left|\mathbf{m}\right|\le {\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+{\mathbf{n}}_{\mathbf{h}}}\mathrm{Tr}\left[T_{\mathbf{j}-\mathbf{m}}^{{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_{\ell }^v}\right]\times \mathrm{Tr}[T_{\mathbf{m}}^{\mathbf{n}}\mathbf{h}+\mathbf{1}{G}_{H_{\ell }}]\right).$$

(34)

Also note that ${\mathbf{n}}_{{\mathbf{s}}^{\mathbf{w}}}$ must be equal to ${\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}$ for the equality (25) to be satisfied. Now, since the equality (29) implies that $V\left(\theta \right)$ is a Lyapunov function in $\mathcal{D}\supseteq \mathfrak{R}$. Along with the fact that $\mathfrak{R}$ is a forward-invariant set, we have the result.

The existence of such a forward-invariant set can be argued as follows. Once the Lyapunov function is obtained, we find a value $c>0$ such that $V^{-1}\left(c\right)\subseteq \mathcal{D}$. Then, the system is locally stable to the origin for all initial conditions from the set $V^{-1}\left([0,c]\right)$ since the boundary points $V^{-1}\left(c\right)$ form a level set and the flow is inwards, by definition. □

Corollary 2.

Let $\dot{\theta }=\left(F_1\left(\theta \right),\dots ,F_d\left(\theta \right)\right)$ be the all-to-all connected and uniform generalized Kuramoto model with the coupling function (15) with ${\alpha }_k\ge 0$ for all k. Suppose $\{f_{\mathbf{k}}^{\left(l\right)}:\left|\mathbf{k}\right|\le {\mathbf{n}}_{\mathbf{f}}\}$ are fourier coefficients of $F_l$, $l\in \{1,\dots ,d\}$; $F\left(0\right)=0$; and $r\in \mathbb{N}$. Let

$$\begin{array}{ccc}\hfill {\mathfrak{S}}_{\pm r}^{(p,q)}& =& \left\{\pm (0,\dots ,0,\underset{p^{th}\mathrm{element}}{r},0,\dots ,0,\underset{q^{th}\mathrm{element}}{-r},0,....0)\right\}\subseteq {\mathbb{Z}}^d.\hfill \end{array}$$

If there exist zero-sum Hermitian matrices $G_{S_0^v},G_{S_{(p,q)}^v}⪰0$ of sizes ${\prod }_{i=1}^d\left(n_s^v\left(i\right)+1\right)$; Hermitian matrices $G_{S_0^w},G_{S_{(p,q)}^w}⪰0$ of sizes ${\prod }_{i=1}^d\left(n_f\left(i\right)+n_s^v\left(i\right)+1\right)$; satisfing

$$\begin{array}{c}\hfill \mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_0}^w\right]+{\displaystyle \frac{1}{2}}\sum _{\begin{array}{c}p,q\in \{1,\dots ,d\}\\ p\ne q\end{array}}\sum _{\mathbf{m}\in {\mathfrak{S}}_{\pm r}^{(p,q)}}\mathrm{Tr}[T_{\mathbf{k}}-{\mathbf{m}}_{\mathbf{f}}^{\mathbf{n}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}{G}_{S_{(p,q)}^w}]=-\sum _{l=1}^d\sum _{\left|\mathbf{j}\right|\le {\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+r·\mathbf{1}}i{j}_l{f}_{\mathbf{k}-\mathbf{j}}^{\left(l\right)}\times \\ \hfill \left(\mathrm{Tr}\left[T_{\mathbf{j}}^{{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_0}^v\right]+{\displaystyle \frac{1}{2}}\sum _{\begin{array}{c}p,q\in \{1,\dots ,d\}\\ p\ne q\end{array}}\sum _{\mathbf{m}\in {\mathfrak{S}}_{\pm r}^{(p,q)}}\mathrm{Tr}[T_{\mathbf{j}}-{\mathbf{m}}_{{\mathbf{s}}^{\mathbf{v}}}^{\mathbf{n}}+\mathbf{1}{G}_{S_{(p,q)}^v}]\right),\end{array}$$

(35)

for all $\left|\mathbf{k}\right|\le {\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+r·\mathbf{1}$, then all solutions of the system (14) originating from $\mathrm{Arc}\left(min\left\{{\displaystyle \frac{\pi }{2r}},{\displaystyle \frac{2}{L}}\left({\displaystyle \frac{\pi }{2}}-{max}_k\left|{\beta }_k\right|\right)\right\}\right)$ converge to the origin.

Proof. 

We are interested in the case when (24) is satisfied for all $\theta \in \mathrm{Arc}(\mathbf{\pi }/2r)$ for some $r\in \mathbb{N}$. Let ${\mathcal{D}}_r=\mathrm{Arc}(\pi /2r)$, then

$$\begin{array}{ccc}\hfill {\mathcal{D}}_r& =& \left\{\theta \in {[-\mathbf{\pi },\mathbf{\pi }]}^d:H_r^{(p,q)}\left(\theta \right)=cos\left(({\theta }_p-{\theta }_q)r\right)\ge 0\forall p,q\in \{1,\dots ,d\},p\ne q\right\}.\hfill \end{array}$$

(36)

Note that $H_r^{(p,q)}\left(\theta \right)={\displaystyle \frac{1}{2}}{\mathrm{e}}^{\mathrm{i}(\mathrm{r}{\theta }_{\mathrm{p}}-\mathrm{r}{\theta }_{\mathrm{q}})}+{\displaystyle \frac{1}{2}}{\mathrm{e}}^{-\mathrm{i}(\mathrm{r}{\theta }_{\mathrm{p}}-\mathrm{r}{\theta }_{\mathrm{q}})}$, and hence is a trigonometric polynomial of degree ${\mathbf{n}}_{\mathbf{h}}=r·\mathbf{1}\in {\mathbb{Z}}^d$. Its Gram matrix $G_{H_r^{(p,q)}}$ is a square matrix of size ${(r+1)}^d$ and, due to trace parameterization, satisfies (4)

$$\mathrm{Tr}\left[T_{\mathbf{m}}^{(r+1)·\mathbf{1}}{G}_{H_r^{(p,q)}}\right]=\left\{\begin{array}{c}\begin{array}{ccc}\hfill & 1/2,\hfill & \hfill \mathrm{if}{m}_p=-m_q=\pm r\\ \hfill & 0,\hfill & \hfill \mathrm{otherwise}\end{array}\hfill \end{array}\right..$$

(37)

Thus, the expression is nonzero only on the set ${\mathfrak{S}}_{\pm r}^{(p,q)}$ given in the hypothesis. We substitute the expression in (29) and observe that the index in ${\sum }_{\ell =1}^{\mathcal{L}}$ changes to ${\sum }_{\begin{array}{c}p,q\in \{1,\dots ,d\}\\ p\ne q\end{array}}$. Also, the summands of ${\sum }_{\mathbf{m}}$ on both sides of (29) are nonzero (and equal to $1/2$) only when $\mathbf{m}\in {\mathfrak{S}}_{\pm r}^{(p,q)}$ due to (37). Hence, the result will follow once we find a forward invariant set in ${\mathcal{D}}_r$.

Consider the set ${\tau }_{\theta }(i,j)=\{t\in {\mathbb{R}}_{+}:{\theta }_i\left(t\right)\equiv {\theta }_j\left(t\right)(mod2\pi )\}$. Due to the analyticity of ${\theta }_k\left(t\right)$ for each $k=1,\dots ,d$, following the proof of Lemma 2.1 [24], there is a dichotomy: the set either contains a limit point, in which case it is equal to ${\mathbb{R}}_{+}$ and natural frequencies ${\omega }_i={\omega }_j$; or the set is countable. We claim that the set $\mathrm{Arc}\left({\displaystyle \frac{2}{L}}\left({\displaystyle \frac{\pi }{2}}-{max}_k\left|{\beta }_k\right|\right)\right)$ is forward-invariant under the flow. Let $\theta \left(0\right)\in \mathrm{Arc}\left({\displaystyle \frac{2}{L}}\left({\displaystyle \frac{\mathbf{\pi }}{2}}-{max}_k\left|{\mathbf{\beta }}_k\right|\right)\right)$, then ${\theta }_k\left(0\right)$’s can be plotted on a nonunique arc of length ${\displaystyle \frac{2}{L}}\left({\displaystyle \frac{\pi }{2}}-{max}_k\left|{\beta }_k\right|\right)$ which is less than $\pi$. Let $m\left(0\right),M\left(0\right)\in \{1,\dots ,d\}$ be the indices of ${\theta }_k\left(0\right)$ occurring first and last, respectively, when we travel clockwise on the arc. Note that the choices $m\left(0\right),M\left(0\right)$ are irrespective of the choice of arc. We extend this definition to $m\left(t\right),M\left(t\right)$ corresponding to $\theta \left(t\right)$. We impose the positivity convention, ${\theta }_{M\left(t\right)}\left(t\right)-{\theta }_{m\left(t\right)}\left(t\right)\in [0,\pi ]$. This means that $\mathrm{Diam}\left(\theta \left(t\right)\right)={\theta }_{M\left(t\right)}\left(t\right)-{\theta }_{m\left(t\right)}\left(t\right)$.

Suppose ${\theta }_k\left(0\right)$ are distinct, then there exists an interval $[0,\tau ]$ on which ${\theta }_k\left(t\right)$’s do not intersect and hence, due to the dichotomy, each ${\tau }_{\theta }(i,j)$ is countable. This implies that functions $m\left(t\right)\equiv m\left(0\right)$ and $M\left(t\right)\equiv M\left(0\right)$ on the interval.

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\mathrm{Diam}\left(\theta \left(t\right)\right)& =& {\dot{\theta }}_{M\left(t\right)}\left(t\right)-{\dot{\theta }}_{m\left(t\right)}\left(t\right)\hfill \\ & =& K\left(\sum _{j=1}^d\left\{g\left({\theta }_j-{\theta }_{M\left(t\right)}\right)-g\left({\theta }_j-{\theta }_{m\left(t\right)}\right)\right\}\right)\hfill \\ & =& -2K\sum _{k=1}^L{\alpha }_k\left[\sum _{j\notin \left\{m\right(t),M(t\left)\right\}}sin\left(k{\displaystyle \frac{{\theta }_{M\left(t\right)}-{\theta }_{m\left(t\right)}}{2}}\right)cos\left(k\left({\theta }_j-{\displaystyle \frac{{\theta }_{M\left(t\right)}+{\theta }_{m\left(t\right)}}{2}}\right)+{\beta }_k\right)\right.\hfill \\ & & \left.+sin\left(k\left({\theta }_{M\left(t\right)}-{\theta }_{m\left(t\right)}\right)\right)cos\left({\beta }_k\right)\right].\hfill \end{array}$$

If $\theta \left(t\right)\in \mathrm{Arc}\left({\displaystyle \frac{2}{L}}\left({\displaystyle \frac{\mathbf{\pi }}{2}}-{max}_k\left|{\mathbf{\beta }}_k\right|\right)\right)$, then for each $k\in \{1,\dots ,L\}$

$$\begin{array}{ccc}\hfill \left|k\left({\theta }_j-{\displaystyle \frac{{\theta }_{M\left(t\right)}+{\theta }_{m\left(t\right)}}{2}}\right)+{\beta }_k\right|& \le & L\left(\left|{\displaystyle \frac{{\theta }_j-{\theta }_{M\left(t\right)}}{2}}\right|+\left|{\displaystyle \frac{{\theta }_j-{\theta }_{m\left(t\right)}}{2}}\right|\right)+\underset{k}{max}\left|{\beta }_k\right|\hfill \\ & =& L\left({\displaystyle \frac{{\theta }_{M\left(t\right)}-{\theta }_j}{2}}+{\displaystyle \frac{{\theta }_j-{\theta }_{m\left(t\right)}}{2}}\right)+\underset{k}{max}\left|{\beta }_k\right|\hfill \\ & \le & L\left({\displaystyle \frac{1}{L}}\left({\displaystyle \frac{\pi }{2}}-\underset{k}{max}\left|{\beta }_k\right|\right)\right)+\underset{k}{max}\left|{\beta }_k\right|={\displaystyle \frac{\pi }{2}},\hfill \end{array}$$

thus, all the cos terms appearing in the expression are nonnegative. Furthermore, ${\theta }_{M\left(t\right)}-{\theta }_{m\left(t\right)}<\pi /L$ and all the sin terms appearing in the expression are nonegative. Thus, $\mathrm{Diam}\left(\theta \right(t\left)\right)$ is nonincreasing in $[0,\tau ]$. The set

$$\left\{t\in {\mathbb{R}}_{+}:M\mathrm{or}m\mathrm{is}\mathrm{discontinuous}\mathrm{at}t\right\}\subseteq \bigcup _{i,j}{\tau }_{\theta }(i,j),$$

hence is countable and can be rewritten as $\{t_1<t_2<\dots \}$. Then $M\left(t\right)$ and $m\left(t\right)$ are constant in each of the intervals $[t_i,t_{i+1}]$ for $i\in \mathbb{N}\cup \left\{0\right\}$ (where $t_0=0$). We can repeat the process above to show that $\mathrm{Diam}\left(\theta \right(t\left)\right)$ is continuous and piecewise decreasing on each of these intervals. Thus if any vector $\theta \left(0\right)$ having distinct entries is contained in $\mathrm{Arc}\left({\displaystyle \frac{2}{L}}\left({\displaystyle \frac{\pi }{2}}-{max}_k\left|{\beta }_k\right|\right)\right)$, then so is its solution trajectory.

When $\theta \left(0\right)$ is such that no two oscillators are equal for all time, the same arguments as above apply since each ${\tau }_{\theta }(i,j)$ is countable. If some oscillators coincide, there might be several choices for $M\left(t\right)$ and $m\left(t\right)$ at given time instants. We can choose the smallest index from $\{1,\dots ,d\}$ for uniqueness and repeat the process by taking

$$\left\{t\in {\mathbb{R}}_{+}:M\mathrm{or}m\mathrm{is}\mathrm{discontinuous}\mathrm{at}t\right\}\subseteq \bigcup _{\begin{array}{c}i,j:{\theta }_i\left(t\right)\neg \equiv {\theta }_j\left(t\right)\forall t\end{array}}{\tau }_{\theta }(i,j).$$

Thus, $\mathrm{Arc}\left({\displaystyle \frac{2}{L}}\left({\displaystyle \frac{\pi }{2}}-{max}_k\left|{\beta }_k\right|\right)\right)$ is forward invariant and so is $\mathrm{Arc}\left(a\right)$ for any $a<{\displaystyle \frac{2}{L}}\left({\displaystyle \frac{\pi }{2}}-{max}_k\left|{\beta }_k\right|\right)$. We thus have forward invariance and the existence of a local Lyapunov function on the set $\mathrm{Arc}\left({\displaystyle \frac{2}{L}}\left({\displaystyle \frac{\pi }{2}}-{max}_k\left|{\beta }_k\right|\right)\right)\cap \mathrm{Arc}\left({\displaystyle \frac{\pi }{2r}}\right)$, hence the result follows. □

Remark 3.

Choosing a larger r in Corollary 2, even though it shrinks the region where the Lyapunov function is being searched for, increases the problem’s computational complexity. Finding a local Lyapunov function for $r+1$ instead of r increases the number of summands on the right-hand side of (35) by

$$\begin{array}{c}\hfill d(2d^2-2d+1)\times \left(\prod _{i=1}^d\left(2n_f\left(i\right)+2n_{s^v}\left(i\right)+2r+3\right)-\prod _{i=1}^d\left(2n_f\left(i\right)+2n_{s^v}\left(i\right)+2r+1\right)\right).\end{array}$$

Corollary 3.

Given a dynamical system $\dot{\theta }=F\left(\theta \right)$ satisfying $F\left(0\right)=0$,where each component function $F_l:{\mathbb{T}}^d\to \mathbb{R}$ has Fourier series expansion (12); and $\epsilon \in (0,d)$. Let

$$\begin{array}{ccc}\hfill {\mathcal{D}}_{\epsilon }& =& \{\theta \in {[-\mathbf{\pi },\mathbf{\pi }]}^d:H_{\mathbf{\epsilon }}\left(\theta \right)=\mathbf{\epsilon }-d+\sum _{i=1}^{d}cos\left({\theta }_i\right)\ge 0\},\mathrm{and}\hfill \\ \hfill \mathfrak{B}& =& \left\{\pm {\mathbf{b}}^{\left(\mathbf{1}\right)},\pm {\mathbf{b}}^{\left(\mathbf{2}\right)},\dots ,\pm {\mathbf{b}}^{\left(\mathbf{d}\right)}:b_j^{\left(i\right)}={\delta }_{ij}\forall i,j\in \{1,\dots ,d\}\right\}\subseteq {\mathbb{Z}}^d.\hfill \end{array}$$

If there exist zero-sum Hermitian matrices $G_{S_0^v},G_{S_1^v}⪰0$ of sizes ${\prod }_{i=1}^d\left(n_{s^v}\left(i\right)+1\right)$; Hermitian matrices $G_{S_0^w},G_{S_1^w}⪰0$ of sizes ${\prod }_{i=1}^d\left(n_f\left(i\right)+n_{s^v}\left(i\right)+1\right)$ satisfying

$$\begin{array}{c}\hfill \mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_0}^w\right]+(\epsilon -d)\mathrm{Tr}\left[T_{\mathbf{k}}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_1}^w\right]+{\displaystyle \frac{1}{2}}\sum _{\mathbf{m}\in \mathfrak{B}}\mathrm{Tr}\left[T_{\mathbf{k}-\mathbf{m}}^{{\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_1^w}\right]=-\sum _{l=1}^d\sum _{\left|\mathbf{j}\right|\le {\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}i\times \\ \hfill j_l{f}_{\mathbf{k}-\mathbf{j}}^{\left(l\right)}\left(\mathrm{Tr}\left[T_{\mathbf{j}}^{{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_0}^v\right]+(\epsilon -d)\mathrm{Tr}\left[T_{\mathbf{j}}^{{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_1}^v\right]+{\displaystyle \frac{1}{2}}\sum _{\mathbf{m}\in \mathfrak{B}}\mathrm{Tr}\left[T_{\mathbf{j}-\mathbf{m}}^{{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}}{G}_{S_1^v}\right]\right),\end{array}$$

for each $\left|\mathbf{k}\right|\le {\mathbf{n}}_{\mathbf{f}}+{\mathbf{n}}_{{\mathbf{s}}^{\mathbf{v}}}+\mathbf{1}$; and a forward-invariant region $\mathfrak{R}\subseteq \mathcal{D}$ containing the origin, then the solutions of (11) with initial conditions in $\mathfrak{R}$ converge to the origin as $t\to \infty$. Hence, the system (11) is locally stable on $\mathfrak{R}$.

Proof. 

The region defined by ${\mathcal{D}}_{\epsilon }$ is a d-dimensional ellipsoid centered at the origin, which contracts as we increase $\epsilon$; refer to Figure 1 (right). By definition of $H_{\epsilon }$, we can observe that ${\mathbf{n}}_{\mathbf{h}}=1$. Let $G_{H_{\epsilon }}$ be its Gram-matrix, then we have

$$\begin{array}{c}\hfill \mathrm{Tr}\left[T_{\mathbf{j}}{G}_{H_{\epsilon }}\right]=\left\{\begin{array}{c}\begin{array}{ccc}& \epsilon -d,\hfill & \hfill \mathbf{j}=\mathbf{0},\\ & 1/2,\hfill & \hfill \mathbf{j}\in \mathfrak{B}\\ & 0,\hfill & \hfill \mathrm{otherwise}\end{array}\hfill \end{array}\right..\end{array}$$

Hence, the result follows from Theorem 2 by rewriting the expression (29). □

Remark 4.

Since the domain ${\mathcal{D}}_{\epsilon }$ in Corollary 3 is constrained by the positivity of exactly one trigonometric polynomial, the computational power required to check the feasibility is sufficiently less. Also, checking for feasibility in a larger or smaller local region does not affect the computation time significantly since the same number of computations are required irrespective of the value of ε unlike for the case discussed in Remark 3.

4. Examples

We now present several examples of coupled systems with three and four oscillators, showcasing the usefulness of our results. Apart from almost global and local results for stability, synchrony, and phase-locking, the results can also be employed to validate the presence of an attracting/ repelling set, phase cohesion, and partial phase-locking of almost all solutions.

Example 1

(Almost global phase Synchronization of uniform all-to-all coupled Oscillators). Consider a system of all-to-all uniformly coupled oscillators. This system under the transformation ${\phi }_j\left(t\right)\mapsto -\omega {\displaystyle \frac{t}{K}}+{\phi }_j\left({\displaystyle \frac{t}{k}}\right)$, reduces to

$${\dot{\phi }}_j=\sum _{k=1}^{d}g({\phi }_k-{\phi }_j),{\phi }_i\in \mathbb{T},i=1,\cdots ,d.$$

(38)

We now exhibit the employment of Theorem 1 to establish the almost global phase synchronization of (38). For this, we consider a well-known coupling function g occurring in [25,26] defined by

$$g\left(x\right)=sin(x+{\beta }_1)-{\alpha }_{2}sin\left(2x\right),$$

(39)

where $0⩽{\beta }_1⩽\pi /2$ and $0⩽{\alpha }_2⩽0.5$. Suppose there are three oscillators ($d=3$). Taking phase-difference variables ${\theta }_1:={\phi }_1-{\phi }_3$, ${\theta }_2:={\phi }_2-{\phi }_3$, we obtain the reduced phase-difference system of Kuramoto oscillators generated by (39) as

$$\dot{{\theta }_1}=f^{\left(1\right)}({\theta }_1,{\theta }_2):=g({\theta }_2-{\theta }_1)+g(-{\theta }_1)-g\left({\theta }_1\right)-g\left({\theta }_2\right),$$

(40)

$${\dot{\theta }}_2=f^{\left(2\right)}({\theta }_1,{\theta }_2):=g({\theta }_1-{\theta }_2)+g(-{\theta }_2)-g\left({\theta }_1\right)-g\left({\theta }_2\right),$$

(41)

where ${\theta }_1,{\theta }_2\in \mathbb{T}$. Almost global phase synchronization of (38) is equivalent to almost global stability of the phase-difference system. The Fourier (exponential) coefficients of the component functions (40) and (41) can be written in the form of $5\times 5$ matrices such that ${\left(f^{\left(1\right)}\right)}_{k_1+3,k_2+3}=f_{(k_1,k_2)}^{\left(1\right)}$ and ${\left(f^{\left(2\right)}\right)}_{k_1+3,k_2+3}=f_{(k_1,k_2)}^{\left(2\right)}$ for $-2\le k_1,k_2\le 2$ as follows

$$\begin{array}{ccc}\hfill \left(f^{\left(1\right)}\right)& =& \left(\begin{array}{ccccc}0& 0& {\alpha }_{2}i& 0& 0.5{\alpha }_{2}i\\ 0& 0& -cos{\beta }_{1}i& -0.5i{e}^{i{\beta }_1}& 0\\ 0.5{\alpha }_{2}i& -0.5i{e}^{-i{\beta }_1}& 0& 0.5i{e}^{i{\beta }_1}& -0.5i{\alpha }_2\\ 0& 0.5i{e}^{-i{\beta }_1}& icos{\beta }_1& 0& 0\\ -0.5i{\alpha }_2& 0& -i{\alpha }_2& 0& 0\end{array}\right),\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill \left(f^{\left(2\right)}\right)& =& \left(\begin{array}{ccccc}0& 0& 0.5i{\alpha }_2& 0& -0.5i{\alpha }_2\\ 0& 0& -0.5i{e}^{-i{\beta }_1}& 0.5i{e}^{-i{\beta }_1}& 0\\ i{\alpha }_2& -icos{\beta }_1& 0& icos{\beta }_1& -i{\alpha }_2\\ 0& -0.5i{e}^{i{\beta }_1}& 0.5i{e}^{i{\beta }_1}& 0& 0\\ 0.5i{\alpha }_2& 0& -0.5i{\alpha }_2& 0& 0\end{array}\right).\hfill \end{array}$$

(43)

Using (42) and (43), Fourier expansion of the system is

$$\begin{array}{cc}\hfill f^{\left(1\right)}({\theta }_1,{\theta }_2)& =\sum _{k_1=-2}^2\sum _{k_2=-2}^2{f}_{(k_1,k_2)}^{\left(1\right)}{e}^{i(k_1{\theta }_1+k_2{\theta }_2)},\hfill \\ \hfill f^{\left(2\right)}({\theta }_1,{\theta }_2)& =\sum _{k_1=-2}^2\sum _{k_2=-2}^2{f}_{(k_1,k_2)}^{\left(2\right)}{e}^{i(k_1{\theta }_1+k_2{\theta }_2)}.\hfill \end{array}$$

(44)

For ${\alpha }_2={\displaystyle \frac{5}{9}}\times 0.5\approx 0.278$ and ${\beta }_1={\displaystyle \frac{5}{9}}\times {\displaystyle \frac{\pi }{2}}\approx 0.785$, employing Theorem 1 with SOS approach we obtain positive semidefinite Hermitian matrices $G_V\in M_{25}\left(\mathbb{C}\right)$ and $G_W\in M_{49}\left(\mathbb{C}\right)$. Using the standard basis of trigonometric polynomials $z\left(\theta \right)$ defined as $z\left(\theta \right)=\left[\begin{array}{cccc}1& e^{i{\theta }_2}& \dots & e^{4i{\theta }_2}\end{array}\right]\otimes \left[\begin{array}{cccc}1& e^{i{\theta }_1}& \dots & e^{4i{\theta }_1}\end{array}\right]$, we obtain $V\left(\theta \right)={\left(z{\left(\theta \right)}^{*}\right)}^T{G}_{V}z\left(\theta \right)$ our Mathematica code [LINK HERE]. The expression is omitted here due to its large size; nonetheless, one such expression will be presented later.

Employing Theorem 1 with SOS approach using our GitHub code, we can obtain positive semidefinite matrices $G_V$ and $G_W$ satisfying (17) for several parameters $({\beta }_1,{\alpha }_2)$ of (44). Due to Theorem 1, almost all solutions of the phase-difference system (44) converge the origin for these parameters $({\alpha }_2,{\beta }_1)$. These parameters values indicated in Figure 2 with blue stars give rise to complete phase synchronization of the coupled oscillators (40) and (41). The region examined in [25] is bounded by the red curve ${\alpha }_2=0.5cos{\beta }_1$, where ${\beta }_1\in [0,\pi /2]$, ${\alpha }_2\in [0,0.5]$.

Figure 2 demonstrates that our approach in Theorem 1 efficiently ensures almost global phase synchronization of the system (44) for parameters $({\beta }_1,{\alpha }_2)$ by producing Lyapnunov densities. Using our method, we can produce Lyapunov densities in the whole region where the almost global phase synchronization of the system (44) has already been established in [25].

The above method can also obtain Lyapunov density for identical oscillators where the connection topology is arbitrary, such as when the underlying topology is a spanning tree. It is known that such a topology gives rise to almost global phase synchronization when all subsystems have identical natural frequencies [27]. Note that the above results were proved without using the dual-Lyapunov method.

Example 2

(Local Phase Synchronization of uniform all-to-all coupled Oscillators for $d=3$ and $d=4$).

Consider the system (38) of three oscillators with coupling function $g\left(x\right)=sin(x+1.2)-0.16sin\left(2x\right)$. To certify the local synchronization of the system (38), we employ Theorem 2. Note that the coefficient of $sin\left(2x\right)$ in (39) is negative, so the forward-invariance of arcs is not easy to prove mathematically, hence Corollary 2 cannot be used directly. We aim to use the equalities implying the existence of a local Lyapunov function for the phase-difference system (40) and (41) in the domain

$${\mathcal{D}}_1=\left\{\theta \in {[-\mathbf{\pi },\mathbf{\pi }]}^2:|{\theta }_1-{\theta }_2|<{\displaystyle \frac{\mathbf{\pi }}{2}}\right\}.$$

Of course, the equalities reduce to ones obtained in Corollary 2 for this domain. The region ${\mathcal{D}}_1$ is plotted in Figure 3 (left). We refer to this figure for the following argument. Using the equalities, we use our code to obtain a local Lyapunov function for the phase-difference system (40) and (41), but it has been skipped due to space constraints. The set ${\mathcal{D}}_1$ is not invariant under the phase-difference flow, as seen in the figure. The subset (of the space of phase-differences) given by

$$\mathfrak{D}=\left\{\theta \in {[-\mathbf{\pi },\mathbf{\pi }]}^2:|{\theta }_1-{\theta }_2|<{\displaystyle \frac{\mathbf{\pi }}{2}},\left|{\theta }_1\right|<{\displaystyle \frac{\mathbf{\pi }}{2}},\left|{\theta }_2\right|<{\displaystyle \frac{\mathbf{\pi }}{2}}\right\}\subseteq {\mathcal{D}}_1,$$

is forward-invariant but is not easy to prove mathematically. Computations show that $V^{-1}\left([0,150]\right)\subseteq \mathfrak{D}$. By the definition of the Lyapunov function, we know that all vectors on the boundary satisfy $V\left(\theta \right)=150$ and point inwards to the region. Hence, $V^{-1}\left([0,150]\right)$ is a forward-invariant set. Thus, local stability of the phase-difference system follows.

We claim that the local stability of phase-difference system in $V^{-1}\left([0,150]\right)$ implies that the original system exhibits local phase synchronization for initial conditions from the set

$$\mathcal{V}=\{\phi \in {[-\mathbf{\pi },\mathbf{\pi }]}^3:V({\phi }_1-{\phi }_2,{\phi }_1-{\phi }_3)\le 150\}.$$

The system of phase-differences defined by (40) and (41) is obtained on realizing the phase-shift symmetry in the system (38). The phase-shift symmetry means that the vector field is invariant under the translation $\phi \mapsto \phi +\mathbf{\epsilon }·\mathbf{1}$ and does not change in the direction of the main diagonal (vector joining $-\pi ·\mathbf{1}$ to $\pi ·\mathbf{1}$). Define an equivalence relation ∼ in ${\mathbb{T}}^3$ as ${\phi }^{\left(1\right)}\sim {\phi }^{\left(2\right)}$ if and only if ${\phi }^{\left(1\right)}-{\phi }^{\left(2\right)}=\epsilon ·1$ for some $\epsilon \in \mathbb{R}$. Thus, the phase-difference vector field on the quotient space ${\mathbb{T}}^3/\sim$ is obtained by projecting the original vector field to any plane perpendicular to the main diagonal. The invariance of $V^{-1}\left([0,150]\right)$ under the phase-difference system is equivalent to the invariance of $\mathcal{V}$ under the original system. Additionally, it is straightforward that the stability of the phase-difference system implies synchronization of the original system. Thus, the system (38) with coupling function $g\left(x\right)=sin(x+1.2)-0.16sin\left(2x\right)$ exhibits local phase synchronization in the region $\mathcal{V}$.

We similarly performed computations for system (38) of four oscillators with coupling function $g\left(x\right)=sin(x+\pi /6)-0.3889sin\left(2x\right)$. Considering the corresponding phase-difference system obtained by substituting ${\theta }_i={\phi }_i-{\phi }_4$ for $i=1,2,3$, we obtained matrices $G_{S_0^V}$, $G_{S_{(1,2)}^V}$, $G_{S_{(1,3)}^V}$ and $G_{S_{(1,2)}^V}$ satisfying (35) for $r=1$. This implies the existence of a local Lyapunov function $V_{*}$ in $\left\{\theta \in {[-\mathbf{\pi },\mathbf{\pi }]}^3:|{\theta }_i-{\theta }_j|<{\displaystyle \frac{\mathbf{\pi }}{2}}\forall i,j\in \{1,2,3\}\right\}$. The function is plotted in Figure 4, and the phase-difference system is locally stable in the region $V_{*}^{-1}\left([0,60]\right)$.

Recall that Theorem 1 and Remark 1 can be used to certify almost global phase-locking of nonuniform oscillators if the singularity of the obtained Lyapunov density for its phase-difference system is a singleton. We demonstrate it with the help of the following example.

Example 3

(Almost global phase-locking in three nonuniform coupled oscillators).

$${\dot{\phi }}_j={\omega }_j+K\sum _{k=1}^{d}g({\phi }_k-{\phi }_j),{\phi }_i\in \mathbb{T},i=1,\cdots ,3,$$

(45)

where $g\left(x\right)=sin(x+0.785)-0.278sin\left(2x\right)$, ${\omega }_1=-1.024$, ${\omega }_2=0.5$, ${\omega }_3=0$ and $K=2$. The corresponding phase-difference system is given by

$$\dot{{\theta }_1}=f^{\left(1\right)}({\theta }_1,{\theta }_2):=-1.024+2(g({\theta }_2-{\theta }_1)+g(-{\theta }_1)-g\left({\theta }_1\right)-g\left({\theta }_2\right)),$$

(46)

$${\dot{\theta }}_2=f^{\left(2\right)}({\theta }_1,{\theta }_2):=0.5+2(g({\theta }_1-{\theta }_2)+g(-{\theta }_2)-g\left({\theta }_1\right)-g\left({\theta }_2\right)),$$

(47)

which differ from (40) and (41) just by a constant. The corresponding Gram-matrices can thus be obtained from (42) and (43) by changing ${\left(f^{\left(1\right)}\right)}_{3,3}=-1.024$ and ${\left(f^{\left(2\right)}\right)}_{3,3}=5.359$, respectively. The system has equibrium points $({\theta }_1,{\theta }_2)=(-0.849402,0.170344)$ and $(-2.09039,0.872725)$. Checking the set of equalities (17) without the zero-sum condition on $G_V$, the code returns a Hermitian matrix $G_V$ whose corresponding polynomial has a unique minimum approximately at $(-0.849148,0.170431)$ according to Mathematica 13.2, where it assumes the value $3.3\times {10}^{-9}$. Thus, it must imply that this point is a zero (approx.) of the function $V\left(\theta \right)$. As discussed in the introduction, such a point must exist since the density $\rho =1/V$ inevitably has a singularity [21]. Note that small errors arise due to approximations in numerical solutions. Thus, the phase-difference system exhibits almost global stability to the point $(-0.849402,0.170344)$. This implies that the original system (45) exhibits almost global phase-locking with $({\theta }_1\left(t\right)-{\theta }_3\left(t\right))\to -0.849402$ and $({\theta }_2\left(t\right)-{\theta }_3\left(t\right))\to 0.170344$ for almost all initial conditions.

In general systems, even more complex dynamics can be captured, such as the phase-cohesiveness of oscillators. An oscillator pair $(i,j)$ is phase-cohesive if there exists a constant $0<c<\pi /2$ and $T_0>0$ such that $|{\theta }_i\left(t\right)-{\theta }_j\left(t\right)|\le c$ for all $t>T_0$. A coupled system is said to be partially-phase-cohesive,28], for $\mathcal{O}⊊\{1,\dots ,d\}$ if there exists $c>0$ and $T_0>0$ such that ${max}_{i,j\in \mathcal{O}}|{\theta }_i\left(t\right)-{\theta }_j\left(t\right)|\le c$ for all $t\ge T_0$.

Example 4.(Almost global partial-phase-cohesion in nonuniform Kuramoto model with a tree-topology: characterizing limit cycles) Consider a coupled nonuniform Kuramoto model (14) of three oscillators with topology $1⟷3⟷2$, meaning $a_{12}=a_{21}=a_{13}=a_{31}=1$, and coupling $K=0.6$. Suppose ${\omega }_1=1$, ${\omega }_2=3$ and ${\omega }_3=1.5$, then using ${\theta }_1={\phi }_1-{\phi }_3$ and ${\theta }_2={\phi }_2-{\phi }_3$, the corresponding phase difference is obtained as

$$\begin{array}{ccc}\hfill {\dot{\theta }}_1& =& -0.5-sin{\theta }_1-0.5sin{\theta }_2=F_1\left(\theta \right)\hfill \end{array}$$

(48)

$$\begin{array}{ccc}\hfill {\dot{\theta }}_2& =& 1.5-sin{\theta }_2-0.5sin{\theta }_1=F_2\left(\theta \right).\hfill \end{array}$$

(49)

The system does not have an equilibrium point. The code returns $9\times 9$ matrices $G_{V_F}$ and $G_{V_B}$ for $\dot{\theta }=F\left(\theta \right)$ and $\dot{\theta }=-F\left(\theta \right)$ respectively, satisfying (17). These give rise to functions $V_F\left(\theta \right)$ and $V_B\left(\theta \right)$, expressions in the Appendix, such that their reciprocals are Lyapunov densities for forward and backward flow, respectively, of the system $\dot{\theta }=F\left(\theta \right)$. This implies that almost all solutions reach $V_B^{-1}\left(0\right)$ in backward time (hence it is a repellor) and reach $V_A^{-1}\left(0\right)$ in the forward time (hence it is an attractor), refer Figure (Figure 4) (right). Notice that almost all trajectories get trapped in the region $\{\theta \in {[-\mathbf{\pi },\mathbf{\pi }]}^3:-1.2<{\theta }_1<0\}$. This implies that oscillators ${\phi }_1$ and ${\phi }_3$ eventually satisfy $|{\phi }_1-{\phi }_3|<1.2$ for almost all initial conditions. Hence, the system is almost globally $(1,3)-$phase-cohesive. Note that the existence of a Lyapunov density for systems with no equilibrium points will certify the presence of a limit cycle.

The reader might have realized that if it so happens that $V_F^{-1}\left(0\right)$ is the line ${\theta }_i={\theta }_0$, then almost all trajectories of the phase-difference system satisfy ${\theta }_i\left(t\right)\to {\theta }_0$, which eventually implies that $({\phi }_i\left(t\right)-{\phi }_3\left(t\right))\to {\theta }_0$. Thus, the original system exhibits almost global (1,3)-phase synchronization. It should be noted if the system has multiple attractors, then the above process will inevitably fail. In such cases, if a point on the limit cycle is known, the code of local stability can be slightly modified to get a local characterization of the limit cycle. Gluing parts of the limit cycle obtained through this procedure could give a bigger picture, in theory, but it seems highly impractical. The methods in this article also extend to coupled oscillators, which depend on the collective interaction of more than two oscillators. For the interaction of three oscillators, these take the form

$${\dot{\phi }}_j={\omega }_j+\sum _{\begin{array}{c}i,k=1\\ i,k\ne j\end{array}}^{d}g({\phi }_i+{\phi }_k-2{\phi }_j),{\phi }_i\in \mathbb{T},i=1,\cdots ,d,$$

and are known to exhibit complex dynamics. We end this section by shedding light upon the fact that, due to phase-shift symmetry, all the methods and observations discussed via the preceding examples can be used as a framework to analyze the dynamics of coupled oscillators with higher-order coupling.

5. Conclusions

We develop a certification of almost global stability of the dynamical system on hypertorus using Gram matrix representation and trace parametrization of multivariate trigonometric polynomials in Theorem 1. We employ Theorem 1 to the system of all-to-all uniformly coupled oscillators of $N=3$ in Example 1. We obtain the parameters plane for certifying the almost global phase synchronization of the system in Figure 2. In Example 3, we obtain a Gram matrix to certify almost global phase-locking in three nonuniform coupled oscillators using Theorem 1. We also use Theorem 1 for almost global partial-phase-cohesion in a nonuniform Kuramoto model with a tree-topology in Example 4. Using Putinar’s Positivstellensatz Lemma 2, we construct a novel approach to certify the local Lyapunov stability on systems on hypertorus in Theorem 2. We apply this theory to certify local phase synchronization of uniform all-to-all coupled for three or four oscillators in Example 2 with local solutions Figure 3.