Construction of Lyapunov Certificates for Systems on Torus Using Trigonometric Polynomials

Amethod for the certification of local and almost global stability of invariant sets of ordinary differential equations on torus is proposed. The method is based on using trigonometric polynomials as candidates of local Lyapunov functions or (global) Lyapunov densities (so-called dual Lyapunov functions) and converting the required local Lyapunov or dual Lyapunov conditions into semidefinite programming. We provide several examples with three and four oscillators to showcase the adaptability of our results, which can be used to certify partial or complete phase-locking/phase-cohesion and to validate the presence of an attracting/repelling set in various scenarios.

Keywords:

almost global synchronization

;

trigonometric polynomials

;

stability certificates

Subject:

Engineering - Control and Systems Engineering

This preprints belongs to the Topic

Modeling, Stability, and Control of Dynamic Systems and Their Applications

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 (θ) : T^{d} = {[0, 2 π)}^{d} \to R

is defined by

R (θ) = \sum_{k = - n_{r}}^{n_{r}} r_{k} e^{i 〈k, θ〉} = \sum_{k_{1} = - n_{r} (1)}^{n_{r} (1)} \dots \sum_{k_{d} = - n_{r} (d)}^{n_{r} (d)} r_{k} e^{i 〈k, θ〉},

(1)

where

n_{r} = (n_{r} (1), n_{r} (2) \dots, n_{r} (d)) \in Z^{d}

is the degree of the trigonometric polynomial,

k = (k_{1}, k_{2}, \dots, k_{d}) \in Z^{d}

;

θ = (θ_{1}, θ_{2}, \dots, θ_{d}) \in T^{d}

; and

r_{k}

are complex coefficients of

R (θ)

. Observe the difference in how the elements of

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,

〈k, θ〉 = \sum_{j = 1}^{d} k_{j} θ_{j}

is an inner product. For real-valued a trigonometric polynomial

R (θ)

, the coefficients

r_{k}

must satisfy the symmetry (Hermitian) condition

r_{- k} = r_{k}^{*}

where

x^{*}

denotes the complex conjugate of x,11]. The inequality

- n_{r} \leq k \leq n_{r}

is understood elementwise and is equivalently denoted by

| k | \leq n_{r}

and

k \in [- n_{r}, n_{r}]

. We denote the vectors

(0, \dots, 0), (1, \dots, 1) \in Z^{d}

by 0 and 1 respectively; and

max {k, j} = (max {k_{1}, j_{1}}, \dots, max {k_{d}, j_{d}})

The Gram matrix

G_{R}

associated with

R (θ)

is defined by

R (θ) = {(z {(θ)}^{*})}^{T} G_{R} z (θ),

(2)

where

G_{R}

is a square matrix of size

\prod_{i = 1}^{d} (n_{r} (i) + 1)

z (θ)

is a vector composed of monomials from the trigonometric basis on

T^{d}

, and

z {(θ)}^{*}

is its conjugate transpose. In particular,

z (θ) = {[1, e^{i θ_{d}}, \dots, e^{i n_{r} (d) θ_{d}}]}^{⊤} \otimes \dots \otimes {[1, e^{i θ_{1}}, \dots, e^{i n_{r} (1) θ_{1}}]}^{⊤}

If the Gram matrix

G_{R}

is positive semidefinite, then

R (θ) \geq 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 (θ) = {∥ U z (θ) ∥}^{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 (θ) = \sum_{ℓ = 1}^{m} {| P_{ℓ} (θ) |}^{2},

(3)

where each

P_{ℓ}

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^{t h}

elementary Toeplitz matrix of size n if it is a

(0, 1) -

matrix with ones on the

k^{t h}

diagonal, meaning

T_{(i, j)} = 1

if and only if

j - i = k

. We will denote the

k^{t h}

elementary Toeplitz matrix of size n by

T_{k}^{n}

. Clearly,

T_{k}^{n}

is a zero matrix for any

k \geq n

. Define

T_{k}^{n} = T_{k_{d}}^{n (d)} \otimes T_{k_{d - 1}}^{n (d - 1)} \otimes \dots \otimes T_{k_{1}}^{n (1)}

where ⊗ denotes the Kronecker product. If

k \notin (- n, n)

, then

k_{i} \geq n (i)

for some i. In that case

T_{k_{i}}^{n (i)}

, hence

T_{k}^{n}

, is the zero matrix. Toeplitz matrices are utilized to relate the coefficients

r_{k}

R (θ)

to the Gram matrix

G_{R}

in (2) by a property known in the literature as trace parameterization [11], which states for any

k \in [- n_{r}, n_{r}]

r_{k} = Tr [T_{k}^{n} r + 1 G_{R}] .

(4)

Positive semidefinite matrices can parametrize SOS trigonometric polynomials [11,18]. Equivalently,

R (θ)

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 (θ)

defined in (1), we can define the derivative with respect to

θ_{l}

as follows

D_{l} R (θ) : = \frac{d R}{d θ_{l}} (θ) = \sum_{k = - n_{r}}^{n_{r}} d r_{k} e^{i 〈k, θ〉}

(5)

where

d r_{k} = i \times k_{l} \times r_{k}

. Thus, due to (4) the relationship between Gram matrix

G_{R}

R (θ)

and Gram matrix

G_{D_{l} R}

D_{l} R (θ)

is given

Tr [T_{k}^{n_{r} + 1} G_{D_{l} R}] = i k_{l} Tr [T_{k}^{n_{r} + 1} G_{R}],

(6)

where

T_{k}^{n_{r} + 1}

are the Toeplitz matrices corresponding to

z (θ)

in (2).

Consider another trigonometric polynomial

P (θ)

defined as

P (θ) = \sum_{k = - n_{p}}^{n_{p}} p_{k} e^{i 〈k, θ〉},

(7)

where

p_{k}

are the coefficients of

P (θ)

. The product of the polynomials

P (θ)

and

R (θ)

, denoted as

P R (θ)

, is also a trigonometric polynomial which can expanded as

P R (θ) = \sum_{k = - n_{pr}}^{n_{pr}} p r_{k} e^{i 〈k, θ〉},

(8)

where

n_{pr} = n_{p} + n_{r}

, and the coefficients

p r_{k}

are given by

p r_{k} = \sum_{j = - n_{pr}}^{n_{pr}} p_{k - j} r_{j} .

(9)

The Gram matrix

G_{P R}

associated with the product

P R (θ)

can be expressed in terms of Gram matrices

G_{P}

and

G_{R}

P (θ)

and

R (θ)

, respectively. Specifically, we have

Tr [T_{k}^{n_{pr} + 1} G_{P R}] = \sum_{j = - n_{pr}}^{n_{pr}} Tr [T_{k - j}^{n_{p} + 1} G_{P}] \times Tr [T_{j}^{n_{r} + 1} G_{R}] .

(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 (0) = 1 G_{R} 1^{⊤} = 0

2.2. A Class of Coupled Oscillators

In this article, we consider a class of coupled oscillators defined on the hypertorus

T^{d}

\dot{θ} = F (θ) = {(F_{1} (θ), \dots, F_{d} (θ))}^{⊤},

(11)

where each component function

F_{l} : T^{d} \to R

has a Fourier series expansion

F_{l} (θ) = \sum_{k = - n_{f}}^{n_{f}} f_{k}^{(l)} e^{i 〈k, θ〉},

(12)

with Fourier coefficients

f_{k}^{(l)}

. 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

n_{f} = {max}_{ℓ = 1}^{d} n_{f} (ℓ) \times 1

and adding zero fourier coefficients wherever necessary. Every smooth vector field on a compact manifold is complete [19], meaning that the flow

Φ_{t} : T^{d} \to T^{d}

is defined for each

t \in R

. Thus, the solutions of (11) exist for all time, irrespective of the initial condition. A set

A \subseteq T^{d}

is said to be forward-invariant (backward-invariant, respectively) under the flow of (11) if

Φ_{t} (A) \subseteq A

for all

t \geq 0

(

t \leq 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 $R ⊊ T^{d}$ such that all solution curves originating from $R$ converge to an equilibrium point in $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} (θ_{i} (t) - θ_{j} (t)) = 0

; phase-locked if

{lim}_{t \to \infty} (θ_{i} (t) - θ_{j} (t)) = θ_{i j}

for some

θ_{i j} \in T

; and frequency synchronized if

{lim}_{t \to \infty} (\frac{d θ_{i} (t)}{d t} - \frac{d θ_{j} (t)}{d t}) = 0

. The system (11) is said to exhibit

almost global phase synchronization (phase-locking, frequency synchronization, respectively): if the set of initial conditions on $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 $R ⊊ T^{d}$ such that the set of initial conditions on $T^{d}$ for which a pair of oscillators is not phase synchronized (not phase-locked, not frequency synchronized, respectively) is a subset of $R^{c}$ .

Under suitable conditions, for example, in the presence of phase-shift symmetry

F (θ + ϵ \cdot 1) = F (θ)

, the system (11) on

T^{d}

can be reduced to a phase-difference system on

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

Θ_{0} \neq 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 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

ρ : A^{c} \to R

and satisfies

C1.: $ρ (θ) > 0$ for almost every $θ \in A^{c}$ .
C2.: ρ is integrable away from A, that is, if $A_{ϵ} = {θ \in T^{d} : {inf}_{ψ \in A} d i s t (θ, ψ) < ϵ}$

$\int_{{(A_{ϵ})}^{c}} ρ (θ) μ (d θ) < \infty \forall ϵ > 0 .$
C3.: $\nabla \cdot (ρ F) (θ) > 0$ for almost every $θ \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

ρ

, whose details can be found in [8]. It is known that a density

ρ

satisfying the evolution constraint inevitably has a singularity in A (Lemma 3, [21]). Considering this and taking

ρ (θ) = V {(θ)}^{- 1}

where

V (θ)

is a trigonometric polynomial, we obtain the following corollary.

Corollary 1.

Assume there exist nonzero trigonometric polynomials

V (θ) \geq 0

and

W (θ) \geq 0

such that

V^{- 1} (0)

is forward-invariant under the flow of (11); and for almost every

θ \in T^{d}

W (θ) = - \nabla V (θ) \cdot F (θ) + V (θ) \times (\nabla \cdot F) (θ),

(13)

then almost all solutions of (11) converge to

V^{- 1} (0)

t \to \infty

. Thus

ρ (θ) = 1 / V (θ)

is a Lyapunov density for (11).

Proof.

Consider

ρ (θ) = 1 / V (θ)

where

V (θ) \geq 0

is a nonzero trigonometric polynomial and take

A = V^{- 1} (0)

. A is the set of all singularities of

ρ (θ)

, 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 (θ) \geq 0

. The set of singularities of

ρ

is nonempty due to Remark 3 [21], thus A is nonempty. Also,

ρ

is positive on

A^{c}

and bounded, hence integrable, over

{(A_{ϵ})}^{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} (0)

. □

2.2.1. Generalized Kuramoto Models with Coupling Topology

Consider the diffusively coupled system of d phase oscillators given by

{\dot{φ}}_{j} = ω_{j} + K \sum_{k = 1}^{d} a_{k j} g (φ_{k} - φ_{j}), φ_{i} \in T, i = 1, \dots, d,

(14)

where

ω_{i}

is the natural frequency of the oscillator

φ_{i}

;

A = (a_{k j})

is a

(0, 1) -

matrix defining the coupling topology, meaning

a_{k j} = 1

if and only if there is an edge from the oscillator

φ_{k}

φ_{j}

(in physical terms,

φ_{k}

influences

φ_{j}

); K is the coupling constant; and g is a coupling function expressed in terms of its first L harmonics

g (x) = \sum_{k = 1}^{L} α_{k} sin (k x + β_{k}),

(15)

where

| β_{k} | < π / 2

for all k. The system is called all-to-all coupled if

a_{k j} = 1

for all

k, j \in {1, \dots, d}

. System of coupled oscillators where all oscillators have identical natural frequencies

ω_{i} = ω_{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 (θ + ϵ \cdot 1) = F (θ)

; 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 (x) = sin (x)

, all-to-all coupled oscillators take the form of the well-known Kuramoto model

{\dot{φ}}_{j} = ω_{j} + K \sum_{k = 1}^{d} sin (φ_{k} - φ_{j}), φ_{i} \in T, i = 1, \dots, d,

(16)

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

When oscillators are nonuniform, that is

ω_{i} \neq ω_{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

ρ (θ) = 1 / V (θ)

in the former case, and a local Lyapunov function

V (θ)

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{θ} = F (θ)

where each component function

F_{l} : T^{d} \to R

has Fourier series expansion (12), if for a nonzero Hermitian matrix

G_{V} ⪰ 0

of size

\prod_{i = 1}^{d} (n_{v} (i) + 1)

having zero-sum, there exist a nonzero Hermitian matrix

G_{W} ⪰ 0

of size

\prod_{i = 1}^{d} (n_{f} (i) + n_{v} (i) + 1)

such that

Tr [T_{k}^{n_{f} + n_{v} + 1} G_{W}] = \sum_{l = 1}^{d} \sum_{j = - (n_{f} + n_{v})}^{n_{f} + n_{v}} i j_{l} (- Tr [T_{j}^{n_{v} + 1} G_{V}] \times f_{k - j}^{(l)} + f_{j}^{(l)} \times Tr [T_{k - j}^{n_{v} + 1} G_{V}]),

(17)

for all

| k | \leq n_{f} + n_{v}

, then the dynamical system satisfies Equation (13) for

V (θ)

and

W (θ)

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 (0) = 0

and (13) is satisfied. The condition

V (0) = 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{matrix} W (θ) & = & \sum_{l = 1}^{d} (- F_{l} (θ) \times \frac{\partial V (θ)}{\partial θ_{l}} + V (θ) \times \frac{\partial F_{l} (θ)}{\partial θ_{l}}), \end{matrix}

(18)

must be satisfied. If

G_{V}

G_{W}

and

G_{F_{l}}

denote the Gram matrices of

V (θ)

W (θ)

and

F_{l} (θ)

respectively, the equation can be rewritten in terms of Gram matrices as

\begin{matrix} G_{W} & = & \sum_{l = 1}^{d} (- G_{F_{l} \times D_{l} V} + G_{V \times D_{l} F_{l}}) . \end{matrix}

(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} {(G_{V})}_{(j, k)} = 0

. Clearly, the size of the matrix

G_{V}

determines the size of

G_{W}

. Assuming that the size of

G_{V}

\prod_{i = 1}^{d} (n_{v} (i) + 1)

and the size of each

G_{F_{l}}

\prod_{i = 1}^{d} (n_{f} (i) + 1)

, the sizes of both

G_{D_{l} V \times F_{l}}

and

G_{D_{l} F_{l} \times V}

\prod_{i = 1}^{d} (n_{f} (i) + n_{v} (i) + 1)

. Thus,

G_{W}

must be of the size

\prod_{i = 1}^{d} (n_{f} (i) + n_{v} (i) + 1)

. Thus for

k \in [- (n_{f} + n_{v}), n_{f} + n_{v}]

, by virtue of trace parametrization (4), the

k^{t h}

coefficient of both sides of (18) are equal if

\begin{matrix} Tr [T_{k}^{n_{f} + n_{v} + 1} G_{W}] & = & \sum_{l = 1}^{d} (- Tr [T_{k}^{n_{f} + n_{v} + 1} G_{F_{l} \times D_{l} V}] + Tr [T_{k}^{n_{f} + n_{v} + 1} G_{V \times D_{l} F_{l}}]) \end{matrix}

(20)

\begin{matrix} = & \sum_{l = 1}^{d} \sum_{j = - (n_{f} + n_{v})}^{n_{f} + n_{v}} (- Tr [T_{j}^{n_{v} + 1} G_{D_{l} V}] \times Tr [T_{k - j}^{n_{f} + 1} G_{F_{l}}] + Tr [T_{j}^{n_{f} + 1} G_{D_{l} F_{l}}] \times Tr [T_{k - j}^{n_{v} + 1} G_{V}]), \end{matrix}

(21)

\begin{matrix} = & \sum_{l = 1}^{d} \sum_{j = - (n_{f} + n_{v})}^{n_{f} + n_{v}} i j_{l} (- Tr [T_{j}^{n_{v} + 1} G_{V}] \times Tr [T_{k - j}^{n_{f} + 1} G_{F_{l}}] + Tr [T_{j}^{n_{f} + 1} G_{F_{l}}] \times Tr [T_{k - j}^{n_{v} + 1} G_{V}]), \end{matrix}

(22)

\begin{matrix} = & \sum_{l = 1}^{d} \sum_{j = - (n_{f} + n_{v})}^{n_{f} + n_{v}} i j_{l} (- Tr [T_{j}^{n_{v} + 1} G_{V}] \times f_{k - j}^{(l)} + f_{j}^{(l)} \times Tr [T_{k - j}^{n_{v} + 1} G_{V}]) . \end{matrix}

(23)

□

Remark 1.

If the condition `the sum of all entries of

G_{V}

is zero’ is dropped in Theorem 1, then

V (0) \neq 0

. However due to [21], the trigonometric polynomial

V (θ) = 1 / ρ (θ)

has a zero. Thus,

V^{- 1} (0)

is a nonempty compact set. If

V^{- 1} (0) \subseteq F^{- 1} (0)

, then the set

V^{- 1} (0)

is clearly invariant under the flow of (11) and almost all solutions of (11) converge to

V^{- 1} (0)

t \to \infty

due to Corollary 1. In particular when

V^{- 1} (0)

is a singleton, say

V^{- 1} (0) = {θ_{0}}

where

θ_{0}

is an equilibrium point of the system (11), then almost all solutions of (11) converge to

θ_{0}

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 (θ) \neg \equiv 0

such that

\nabla V (θ) \cdot F (θ)

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 (θ)

satisfying

W (θ) = - \nabla V (θ) \cdot F (θ) .

(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 (θ)

and

W (θ)

if their corresponding Gram matrices

G_{V}

and

G_{W}

satisfy

Tr [T_{k}^{n_{f} + n_{v} + 1} G_{W}] = - \sum_{l = 1}^{d} \sum_{j = - (n_{f} + n_{v})}^{n_{f} + n_{v}} i j_{l} Tr [T_{j}^{n_{v} + 1} G_{V}] \times f_{k - j}^{(l)},

(25)

for all

k \in [- (n_{f} + n_{v}), n_{f} + n_{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

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_{ℓ}

be trigonometric polynomials (say, of degree

n_{h}

) with complex coefficients for

ℓ = 1, \dots, L

and

D

be the domain

D = θ \in {[- π, π]}^{d} : H_{ℓ} (θ) \geq 0, ℓ = 1, \dots, L .

(26)

If a trigonometric polynomial

R (θ)

is positive on

D

(that is,

R (θ) > 0

for all

θ \in D

), then there exist sum-of-squares polynomials

S_{ℓ} (θ)

for

ℓ = 0, \dots, L

such that

R (θ) = S_{0} (θ) + \sum_{ℓ = 1}^{L} S_{ℓ} (θ) H_{ℓ} (θ) .

(27)

Without loss of generality, it may be assumed that all trigonometric polynomials

S_{ℓ}

share the same degree

n_{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 (θ)

has degree

n_{r}

and is positive on

D

defined as in (26), then there exist nonzero Hermitian matrices

G_{S_{ℓ}} ⪰ 0

of sizes

\prod_{i = 1}^{d} (n_{s} (i) + 1)

for

ℓ = 0, \dots, L

such that for any

| k | \leq max {n_{s} + n_{h}, n_{r}}

r_{k} = Tr [T_{k}^{n_{s} + 1} G_{S_{0}}] + \sum_{ℓ = 1}^{L} \sum_{j = - (n_{h} + n_{s})}^{n_{h} + n_{s}} Tr [T_{k} - j_{s}^{n} + 1 G_{S_{ℓ}}] \times Tr [T_{k}^{n_{h} + 1} G_{H_{ℓ}}] .

(28)

Remark 2.

A vector

θ \in T^{d}

is said to have diameter a, which is denoted by

Diam (θ) = a

, if

{max}_{i, j \in {1, \dots, d}} | θ_{i} - θ_{j} | = a

. The set

Arc (a) = {θ \in T^{d} : Diag (θ) < a}

which is a 2d-gonal hypercylinder with “diagonal" joining

- π \cdot 1

and

π \cdot 1

as its axis. Refer to Figure 1 (left) where

a = π / 2

d = 3

, and the region is a hexagonal cylinder. The terminology

Arc (a)

is used because all components

θ_{j}

of any given

θ \in Arc (a)

fall on an arc of length a in

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

Arc (π / 2)

is forward-invariant. They subsequently find a Lyapunov function in this domain. The conditions ensure that all solutions starting in

Arc (π / 2)

get phase synchronized.

Theorem 2.

Given a dynamical system

\dot{θ} = F (θ)

satisfying

F (0) = 0

,where each component function

F_{l} : T^{d} \to R

has Fourier series expansion (12);

H_{ℓ}

be trigonometric polynomials with Gram-matrices

G_{H_{ℓ}}

for

ℓ = 1, \dots, L

; and

D = {θ \in {[- π, π]}^{d} : H_{ℓ} (θ) \geq 0, ℓ = 1, \dots, L} .

If there exist zero-sum Hermitian matrices

G_{S_{ℓ}^{v}} ⪰ 0

of sizes

\prod_{i = 1}^{d} (n_{s^{v}} (i) + 1)

; Hermitian matrices

G_{S_{ℓ}^{w}} ⪰ 0

of sizes

\prod_{i = 1}^{d} (n_{f} (i) + n_{s^{v}} (i) + 1)

for each

ℓ \in {1, \dots, L}

satisfying

\begin{matrix} Tr [T_{k}^{n_{f} + n_{s^{v}} + 1} G_{S_{0}}^{w}] + \sum_{ℓ = 1}^{L} \sum_{| m | \leq n_{f} + n_{s^{v}} + n_{h}} Tr [T_{k - m}^{n_{f} + n_{s^{v}} + 1} G_{S_{ℓ}^{w}}] \times Tr [T_{m}^{n} h + 1 G_{H_{ℓ}}] = - \sum_{l = 1}^{d} \sum_{| j | \leq n_{f} + n_{s^{v}} + n_{h}} i \times \\ j_{l} f_{k - j}^{(l)} (Tr [T_{j}^{n_{s^{v}} + 1} G_{S_{0}}^{v}] + \sum_{ℓ = 1}^{L} \sum_{| m | \leq n_{s^{v}} + n_{h}} Tr [T_{j - m}^{n_{s^{v}} + 1} G_{S_{ℓ}^{v}}] \times Tr [T_{m}^{n} h + 1 G_{H_{ℓ}}]), \end{matrix}

(29)

for each

| k | \leq n_{f} + n_{s^{v}} + n_{h}

; and a forward-invariant region

R \subseteq D

containing the origin, then the solutions of (11) with initial conditions in

R

converge to the origin as

t \to \infty

. Hence, the system (11) is locally stable on

R

Proof.

We are interested in finding trigonometric polynomials

V (θ)

and

W (θ)

which are positive on

D

and satisfy (24) for all

θ \in D

. Note that each

H_{ℓ} (θ)

is a trigonometric polynomial of degree

n_{h}

. Thus, by Putinar’s Positivstellensatz V and W must have the forms

\begin{matrix} V (θ) & = & S_{0}^{v} (θ) + \sum_{ℓ = 1}^{L} H_{ℓ} (θ) S_{ℓ}^{v} (θ), \end{matrix}

(30)

\begin{matrix} W (θ) & = & S_{0}^{w} (θ) + \sum_{ℓ = 1}^{L} H_{ℓ} (θ) S_{ℓ}^{w} (θ), \end{matrix}

(31)

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

\begin{matrix} Tr [T_{k}^{n_{v} + 1} G_{V}] & = & Tr [T_{k}^{n_{s^{v}} + 1} G_{S_{0}}^{v}] + \sum_{ℓ = 1}^{L} \sum_{| m | \leq n_{s^{v}} + n_{h}} Tr [T_{k - m}^{n_{s^{v}} + 1} G_{S_{ℓ}^{v}}] \times Tr [T_{m}^{n} h + 1 G_{H_{ℓ}}], \end{matrix}

(32)

\begin{matrix} Tr [T_{k}^{n_{w} + 1} G_{W}] & = & Tr [T_{k}^{n_{s^{w}} + 1} G_{S_{0}}^{w}] + \sum_{ℓ = 1}^{L} \sum_{| m | \leq n_{s^{w}} + n_{h}} Tr [T_{k - m}^{n_{s^{w}} + 1} G_{S_{ℓ}^{w}}] \times Tr [T_{m}^{n} h + 1 G_{H_{ℓ}}] . \end{matrix}

(33)

Clearly,

n_{v} = n_{s^{v}} + n_{h}

and

n_{w} = n_{s^{w}} + n_{h}

. Note that

G_{S_{ℓ}^{v}}

are zero-sum matrices, hence

S_{ℓ}^{v} (0) = 0

, for all

ℓ \in {1, \dots, L}

. Thus,

V (0) = 0

. Also (32) occurs as the left side of (25). For the right side, substitute the first equality in Equation (25) to obtain

{RHS}_{e q . ()} = - \sum_{l = 1}^{d} \sum_{| j | \leq n_{f} + n_{s^{v}} + n_{h}} i j_{l} f_{k - j}^{(l)} (Tr [T_{j}^{n_{s^{v}} + 1} G_{S_{0}}^{v}] + \sum_{ℓ = 1}^{L} \sum_{| m | \leq n_{s^{v}} + n_{h}} Tr [T_{j - m}^{n_{s^{v}} + 1} G_{S_{ℓ}^{v}}] \times Tr [T_{m}^{n} h + 1 G_{H_{ℓ}}]) .

(34)

Also note that

n_{s^{w}}

must be equal to

n_{f} + n_{s^{v}}

for the equality (25) to be satisfied. Now, since the equality (29) implies that

V (θ)

is a Lyapunov function in

D \supseteq R

. Along with the fact that

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} (c) \subseteq D

. Then, the system is locally stable to the origin for all initial conditions from the set

V^{- 1} ([0, c])

since the boundary points

V^{- 1} (c)

form a level set and the flow is inwards, by definition. □

Corollary 2.

Let

\dot{θ} = (F_{1} (θ), \dots, F_{d} (θ))

be the all-to-all connected and uniform generalized Kuramoto model with the coupling function (15) with

α_{k} \geq 0

for all k. Suppose

{f_{k}^{(l)} : | k | \leq n_{f}}

are fourier coefficients of

F_{l}

l \in {1, \dots, d}

;

F (0) = 0

; and

r \in N

. Let

\begin{matrix} S_{\pm r}^{(p, q)} & = & \{\pm (0, \dots, 0, \underset{p^{t h} element}{r}, 0, \dots, 0, \underset{q^{t h} element}{- r}, 0, . . . . 0)\} \subseteq Z^{d} . \end{matrix}

If there exist zero-sum Hermitian matrices

G_{S_{0}^{v}}, G_{S_{(p, q)}^{v}} ⪰ 0

of sizes

\prod_{i = 1}^{d} (n_{s}^{v} (i) + 1)

; Hermitian matrices

G_{S_{0}^{w}}, G_{S_{(p, q)}^{w}} ⪰ 0

of sizes

\prod_{i = 1}^{d} (n_{f} (i) + n_{s}^{v} (i) + 1)

; satisfing

\begin{matrix} Tr [T_{k}^{n_{f} + n_{s^{v}} + 1} G_{S_{0}}^{w}] + \frac{1}{2} \sum_{\begin{matrix} p, q \in {1, \dots, d} \\ p \neq q \end{matrix}} \sum_{m \in S_{\pm r}^{(p, q)}} Tr [T_{k} - m_{f}^{n} + n_{s^{v}} + 1 G_{S_{(p, q)}^{w}}] = - \sum_{l = 1}^{d} \sum_{| j | \leq n_{f} + n_{s^{v}} + r \cdot 1} i j_{l} f_{k - j}^{(l)} \times \\ (Tr [T_{j}^{n_{s^{v}} + 1} G_{S_{0}}^{v}] + \frac{1}{2} \sum_{\begin{matrix} p, q \in {1, \dots, d} \\ p \neq q \end{matrix}} \sum_{m \in S_{\pm r}^{(p, q)}} Tr [T_{j} - m_{s^{v}}^{n} + 1 G_{S_{(p, q)}^{v}}]), \end{matrix}

(35)

for all

| k | \leq n_{f} + n_{s^{v}} + r \cdot 1

, then all solutions of the system (14) originating from

Arc (min \{\frac{π}{2 r}, \frac{2}{L} (\frac{π}{2} - {max}_{k} | β_{k} |)\})

converge to the origin.

Proof.

We are interested in the case when (24) is satisfied for all

θ \in Arc (π / 2 r)

for some

r \in N

. Let

D_{r} = Arc (π / 2 r)

, then

\begin{matrix} D_{r} & = & \{θ \in {[- π, π]}^{d} : H_{r}^{(p, q)} (θ) = cos ((θ_{p} - θ_{q}) r) \geq 0 \forall p, q \in {1, \dots, d}, p \neq q\} . \end{matrix}

(36)

Note that

H_{r}^{(p, q)} (θ) = \frac{1}{2} e^{i (r θ_{p} - r θ_{q})} + \frac{1}{2} e^{- i (r θ_{p} - r θ_{q})}

, and hence is a trigonometric polynomial of degree

n_{h} = r \cdot 1 \in 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)

Tr [T_{m}^{(r + 1) \cdot 1} G_{H_{r}^{(p, q)}}] = \{\begin{matrix} \begin{matrix} 1 / 2, & if m_{p} = - m_{q} = \pm r \\ 0, & otherwise \end{matrix} \end{matrix} .

(37)

Thus, the expression is nonzero only on the set

S_{\pm r}^{(p, q)}

given in the hypothesis. We substitute the expression in (29) and observe that the index in

\sum_{ℓ = 1}^{L}

changes to

\sum_{\begin{matrix} p, q \in {1, \dots, d} \\ p \neq q \end{matrix}}

. Also, the summands of

\sum_{m}

on both sides of (29) are nonzero (and equal to

1 / 2

) only when

m \in S_{\pm r}^{(p, q)}

due to (37). Hence, the result will follow once we find a forward invariant set in

D_{r}

Consider the set

τ_{θ} (i, j) = {t \in R_{+} : θ_{i} (t) \equiv θ_{j} (t) (mod 2 π)}

. Due to the analyticity of

θ_{k} (t)

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

R_{+}

and natural frequencies

ω_{i} = ω_{j}

; or the set is countable. We claim that the set

Arc (\frac{2}{L} (\frac{π}{2} - {max}_{k} | β_{k} |))

is forward-invariant under the flow. Let

θ (0) \in Arc (\frac{2}{L} (\frac{π}{2} - {max}_{k} | β_{k} |))

, then

θ_{k} (0)

’s can be plotted on a nonunique arc of length

\frac{2}{L} (\frac{π}{2} - {max}_{k} | β_{k} |)

which is less than

π

. Let

m (0), M (0) \in {1, \dots, d}

be the indices of

θ_{k} (0)

occurring first and last, respectively, when we travel clockwise on the arc. Note that the choices

m (0), M (0)

are irrespective of the choice of arc. We extend this definition to

m (t), M (t)

corresponding to

θ (t)

. We impose the positivity convention,

θ_{M (t)} (t) - θ_{m (t)} (t) \in [0, π]

. This means that

Diam (θ (t)) = θ_{M (t)} (t) - θ_{m (t)} (t)

Suppose

θ_{k} (0)

are distinct, then there exists an interval

[0, τ]

on which

θ_{k} (t)

’s do not intersect and hence, due to the dichotomy, each

τ_{θ} (i, j)

is countable. This implies that functions

m (t) \equiv m (0)

and

M (t) \equiv M (0)

on the interval.

\begin{matrix} \frac{d}{d t} Diam (θ (t)) & = & {\dot{θ}}_{M (t)} (t) - {\dot{θ}}_{m (t)} (t) \\ = & K (\sum_{j = 1}^{d} \{g (θ_{j} - θ_{M (t)}) - g (θ_{j} - θ_{m (t)})\}) \\ = & - 2 K \sum_{k = 1}^{L} α_{k} [\sum_{j \notin {m (t), M (t)}} sin (k \frac{θ_{M (t)} - θ_{m (t)}}{2}) cos (k (θ_{j} - \frac{θ_{M (t)} + θ_{m (t)}}{2}) + β_{k}) \\ + sin (k (θ_{M (t)} - θ_{m (t)})) cos (β_{k})] . \end{matrix}

θ (t) \in Arc (\frac{2}{L} (\frac{π}{2} - {max}_{k} | β_{k} |))

, then for each

k \in {1, \dots, L}

\begin{matrix} |k (θ_{j} - \frac{θ_{M (t)} + θ_{m (t)}}{2}) + β_{k}| & \leq & L (|\frac{θ_{j} - θ_{M (t)}}{2}| + |\frac{θ_{j} - θ_{m (t)}}{2}|) + max_{k} | β_{k} | \\ = & L (\frac{θ_{M (t)} - θ_{j}}{2} + \frac{θ_{j} - θ_{m (t)}}{2}) + max_{k} | β_{k} | \\ \leq & L (\frac{1}{L} (\frac{π}{2} - max_{k} | β_{k} |)) + max_{k} | β_{k} | = \frac{π}{2}, \end{matrix}

thus, all the cos terms appearing in the expression are nonnegative. Furthermore,

θ_{M (t)} - θ_{m (t)} < π / L

and all the sin terms appearing in the expression are nonegative. Thus,

Diam (θ (t))

is nonincreasing in

[0, τ]

. The set

\{t \in R_{+} : M or m is discontinuous at t\} \subseteq ⋃_{i, j} τ_{θ} (i, j),

hence is countable and can be rewritten as

{t_{1} < t_{2} < \dots}

. Then

M (t)

and

m (t)

are constant in each of the intervals

[t_{i}, t_{i + 1}]

for

i \in N \cup {0}

(where

t_{0} = 0

). We can repeat the process above to show that

Diam (θ (t))

is continuous and piecewise decreasing on each of these intervals. Thus if any vector

θ (0)

having distinct entries is contained in

Arc (\frac{2}{L} (\frac{π}{2} - {max}_{k} | β_{k} |))

, then so is its solution trajectory.

When

θ (0)

is such that no two oscillators are equal for all time, the same arguments as above apply since each

τ_{θ} (i, j)

is countable. If some oscillators coincide, there might be several choices for

M (t)

and

m (t)

at given time instants. We can choose the smallest index from

{1, \dots, d}

for uniqueness and repeat the process by taking

\{t \in R_{+} : M or m is discontinuous at t\} \subseteq ⋃_{\begin{matrix} i, j : θ_{i} (t) \neg \equiv θ_{j} (t) \forall t \end{matrix}} τ_{θ} (i, j) .

Thus,

Arc (\frac{2}{L} (\frac{π}{2} - {max}_{k} | β_{k} |))

is forward invariant and so is

Arc (a)

for any

a < \frac{2}{L} (\frac{π}{2} - {max}_{k} | β_{k} |)

. We thus have forward invariance and the existence of a local Lyapunov function on the set

Arc (\frac{2}{L} (\frac{π}{2} - {max}_{k} | β_{k} |)) \cap Arc (\frac{π}{2 r})

, 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{matrix} d (2 d^{2} - 2 d + 1) \times (\prod_{i = 1}^{d} (2 n_{f} (i) + 2 n_{s^{v}} (i) + 2 r + 3) - \prod_{i = 1}^{d} (2 n_{f} (i) + 2 n_{s^{v}} (i) + 2 r + 1)) . \end{matrix}

Corollary 3.

Given a dynamical system

\dot{θ} = F (θ)

satisfying

F (0) = 0

,where each component function

F_{l} : T^{d} \to R

has Fourier series expansion (12); and

ε \in (0, d)

. Let

\begin{matrix} D_{ε} & = & {θ \in {[- π, π]}^{d} : H_{ε} (θ) = ε - d + \sum_{i = 1}^{d} cos (θ_{i}) \geq 0}, and \\ B & = & \{\pm b^{(1)}, \pm b^{(2)}, \dots, \pm b^{(d)} : b_{j}^{(i)} = δ_{i j} \forall i, j \in {1, \dots, d}\} \subseteq Z^{d} . \end{matrix}

If there exist zero-sum Hermitian matrices

G_{S_{0}^{v}}, G_{S_{1}^{v}} ⪰ 0

of sizes

\prod_{i = 1}^{d} (n_{s^{v}} (i) + 1)

; Hermitian matrices

G_{S_{0}^{w}}, G_{S_{1}^{w}} ⪰ 0

of sizes

\prod_{i = 1}^{d} (n_{f} (i) + n_{s^{v}} (i) + 1)

satisfying

\begin{matrix} Tr [T_{k}^{n_{f} + n_{s^{v}} + 1} G_{S_{0}}^{w}] + (ε - d) Tr [T_{k}^{n_{f} + n_{s^{v}} + 1} G_{S_{1}}^{w}] + \frac{1}{2} \sum_{m \in B} Tr [T_{k - m}^{n_{f} + n_{s^{v}} + 1} G_{S_{1}^{w}}] = - \sum_{l = 1}^{d} \sum_{| j | \leq n_{f} + n_{s^{v}} + 1} i \times \\ j_{l} f_{k - j}^{(l)} (Tr [T_{j}^{n_{s^{v}} + 1} G_{S_{0}}^{v}] + (ε - d) Tr [T_{j}^{n_{s^{v}} + 1} G_{S_{1}}^{v}] + \frac{1}{2} \sum_{m \in B} Tr [T_{j - m}^{n_{s^{v}} + 1} G_{S_{1}^{v}}]), \end{matrix}

for each

| k | \leq n_{f} + n_{s^{v}} + 1

; and a forward-invariant region

R \subseteq D

containing the origin, then the solutions of (11) with initial conditions in

R

converge to the origin as

t \to \infty

. Hence, the system (11) is locally stable on

R

Proof.

The region defined by

D_{ε}

is a d-dimensional ellipsoid centered at the origin, which contracts as we increase

ε

; refer to Figure 1 (right). By definition of

H_{ε}

, we can observe that

n_{h} = 1

. Let

G_{H_{ε}}

be its Gram-matrix, then we have

\begin{matrix} Tr [T_{j} G_{H_{ε}}] = \{\begin{matrix} \begin{matrix} ε - d, & j = 0, \\ 1 / 2, & j \in B \\ 0, & otherwise \end{matrix} \end{matrix} . \end{matrix}

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

Remark 4.

Since the domain

D_{ε}

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

φ_{j} (t) \mapsto - ω \frac{t}{K} + φ_{j} (\frac{t}{k})

, reduces to

{\dot{φ}}_{j} = \sum_{k = 1}^{d} g (φ_{k} - φ_{j}), φ_{i} \in T, i = 1, \dots, 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 (x) = sin (x + β_{1}) - α_{2} sin (2 x),

(39)

where

0 ⩽ β_{1} ⩽ π / 2

and

0 ⩽ α_{2} ⩽ 0.5

. Suppose there are three oscillators (

d = 3

). Taking phase-difference variables

θ_{1} : = φ_{1} - φ_{3}

θ_{2} : = φ_{2} - φ_{3}

, we obtain the reduced phase-difference system of Kuramoto oscillators generated by (39) as

\dot{θ_{1}} = f^{(1)} (θ_{1}, θ_{2}) : = g (θ_{2} - θ_{1}) + g (- θ_{1}) - g (θ_{1}) - g (θ_{2}),

(40)

{\dot{θ}}_{2} = f^{(2)} (θ_{1}, θ_{2}) : = g (θ_{1} - θ_{2}) + g (- θ_{2}) - g (θ_{1}) - g (θ_{2}),

(41)

where

θ_{1}, θ_{2} \in 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

{(f^{(1)})}_{k_{1} + 3, k_{2} + 3} = f_{(k_{1}, k_{2})}^{(1)}

and

{(f^{(2)})}_{k_{1} + 3, k_{2} + 3} = f_{(k_{1}, k_{2})}^{(2)}

for

- 2 \leq k_{1}, k_{2} \leq 2

as follows

\begin{matrix} (f^{(1)}) & = & (\begin{matrix} 0 & 0 & α_{2} i & 0 & 0.5 α_{2} i \\ 0 & 0 & - cos β_{1} i & - 0.5 i e^{i β_{1}} & 0 \\ 0.5 α_{2} i & - 0.5 i e^{- i β_{1}} & 0 & 0.5 i e^{i β_{1}} & - 0.5 i α_{2} \\ 0 & 0.5 i e^{- i β_{1}} & i cos β_{1} & 0 & 0 \\ - 0.5 i α_{2} & 0 & - i α_{2} & 0 & 0 \end{matrix}), \end{matrix}

(42)

\begin{matrix} (f^{(2)}) & = & (\begin{matrix} 0 & 0 & 0.5 i α_{2} & 0 & - 0.5 i α_{2} \\ 0 & 0 & - 0.5 i e^{- i β_{1}} & 0.5 i e^{- i β_{1}} & 0 \\ i α_{2} & - i cos β_{1} & 0 & i cos β_{1} & - i α_{2} \\ 0 & - 0.5 i e^{i β_{1}} & 0.5 i e^{i β_{1}} & 0 & 0 \\ 0.5 i α_{2} & 0 & - 0.5 i α_{2} & 0 & 0 \end{matrix}) . \end{matrix}

(43)

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

\begin{matrix} f^{(1)} (θ_{1}, θ_{2}) & = \sum_{k_{1} = - 2}^{2} \sum_{k_{2} = - 2}^{2} f_{(k_{1}, k_{2})}^{(1)} e^{i (k_{1} θ_{1} + k_{2} θ_{2})}, \\ f^{(2)} (θ_{1}, θ_{2}) & = \sum_{k_{1} = - 2}^{2} \sum_{k_{2} = - 2}^{2} f_{(k_{1}, k_{2})}^{(2)} e^{i (k_{1} θ_{1} + k_{2} θ_{2})} . \end{matrix}

(44)

For

α_{2} = \frac{5}{9} \times 0.5 \approx 0.278

and

β_{1} = \frac{5}{9} \times \frac{π}{2} \approx 0.785

, employing Theorem 1 with SOS approach we obtain positive semidefinite Hermitian matrices

G_{V} \in M_{25} (C)

and

G_{W} \in M_{49} (C)

. Using the standard basis of trigonometric polynomials

z (θ)

defined as

z (θ) = [\begin{matrix} 1 & e^{i θ_{2}} & \dots & e^{4 i θ_{2}} \end{matrix}] \otimes [\begin{matrix} 1 & e^{i θ_{1}} & \dots & e^{4 i θ_{1}} \end{matrix}]

, we obtain

V (θ) = {(z {(θ)}^{*})}^{T} G_{V} z (θ)

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

(β_{1}, α_{2})

of (44). Due to Theorem 1, almost all solutions of the phase-difference system (44) converge the origin for these parameters

(α_{2}, β_{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

α_{2} = 0.5 cos β_{1}

, where

β_{1} \in [0, π / 2]

α_{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

(β_{1}, α_{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 (x) = sin (x + 1.2) - 0.16 sin (2 x)

. To certify the local synchronization of the system (38), we employ Theorem 2. Note that the coefficient of

s i n (2 x)

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

D_{1} = \{θ \in {[- π, π]}^{2} : | θ_{1} - θ_{2} | < \frac{π}{2}\} .

Of course, the equalities reduce to ones obtained in Corollary 2 for this domain. The region

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

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

D = \{θ \in {[- π, π]}^{2} : | θ_{1} - θ_{2} | < \frac{π}{2}, | θ_{1} | < \frac{π}{2}, | θ_{2} | < \frac{π}{2}\} \subseteq D_{1},

is forward-invariant but is not easy to prove mathematically. Computations show that

V^{- 1} ([0, 150]) \subseteq D

. By the definition of the Lyapunov function, we know that all vectors on the boundary satisfy

V (θ) = 150

and point inwards to the region. Hence,

V^{- 1} ([0, 150])

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} ([0, 150])

implies that the original system exhibits local phase synchronization for initial conditions from the set

V = {φ \in {[- π, π]}^{3} : V (φ_{1} - φ_{2}, φ_{1} - φ_{3}) \leq 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

φ \mapsto φ + ε \cdot 1

and does not change in the direction of the main diagonal (vector joining

- π \cdot 1

π \cdot 1

). Define an equivalence relation ∼ in

T^{3}

φ^{(1)} \sim φ^{(2)}

if and only if

φ^{(1)} - φ^{(2)} = ε \cdot 1

for some

ε \in R

. Thus, the phase-difference vector field on the quotient space

T^{3} / \sim

is obtained by projecting the original vector field to any plane perpendicular to the main diagonal. The invariance of

V^{- 1} ([0, 150])

under the phase-difference system is equivalent to the invariance of

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 (x) = sin (x + 1.2) - 0.16 sin (2 x)

exhibits local phase synchronization in the region

V

We similarly performed computations for system (38) of four oscillators with coupling function

g (x) = sin (x + π / 6) - 0.3889 sin (2 x)

. Considering the corresponding phase-difference system obtained by substituting

θ_{i} = φ_{i} - φ_{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 {[- π, π]}^{3} : | θ_{i} - θ_{j} | < \frac{π}{2} \forall i, j \in {1, 2, 3}\}

. The function is plotted in Figure 4, and the phase-difference system is locally stable in the region

V_{*}^{- 1} ([0, 60])

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{φ}}_{j} = ω_{j} + K \sum_{k = 1}^{d} g (φ_{k} - φ_{j}), φ_{i} \in T, i = 1, \dots, 3,

(45)

where

g (x) = sin (x + 0.785) - 0.278 sin (2 x)

ω_{1} = - 1.024

ω_{2} = 0.5

ω_{3} = 0

and

K = 2

. The corresponding phase-difference system is given by

\dot{θ_{1}} = f^{(1)} (θ_{1}, θ_{2}) : = - 1.024 + 2 (g (θ_{2} - θ_{1}) + g (- θ_{1}) - g (θ_{1}) - g (θ_{2})),

(46)

{\dot{θ}}_{2} = f^{(2)} (θ_{1}, θ_{2}) : = 0.5 + 2 (g (θ_{1} - θ_{2}) + g (- θ_{2}) - g (θ_{1}) - g (θ_{2})),

(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

{(f^{(1)})}_{3, 3} = - 1.024

and

{(f^{(2)})}_{3, 3} = 5.359

, respectively. The system has equibrium points

(θ_{1}, θ_{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 (θ)

. As discussed in the introduction, such a point must exist since the density

ρ = 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

(θ_{1} (t) - θ_{3} (t)) \to - 0.849402

and

(θ_{2} (t) - θ_{3} (t)) \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 < π / 2

and

T_{0} > 0

such that

| θ_{i} (t) - θ_{j} (t) | \leq c

for all

t > T_{0}

. A coupled system is said to be partially-phase-cohesive,28], for

O ⊊ {1, \dots, d}

if there exists

c > 0

and

T_{0} > 0

such that

{max}_{i, j \in O} | θ_{i} (t) - θ_{j} (t) | \leq c

for all

t \geq 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

ω_{1} = 1

ω_{2} = 3

and

ω_{3} = 1.5

, then using

θ_{1} = φ_{1} - φ_{3}

and

θ_{2} = φ_{2} - φ_{3}

, the corresponding phase difference is obtained as

\begin{matrix} {\dot{θ}}_{1} & = & - 0.5 - sin θ_{1} - 0.5 sin θ_{2} = F_{1} (θ) \end{matrix}

(48)

\begin{matrix} {\dot{θ}}_{2} & = & 1.5 - sin θ_{2} - 0.5 sin θ_{1} = F_{2} (θ) . \end{matrix}

(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{θ} = F (θ)

and

\dot{θ} = - F (θ)

respectively, satisfying (17). These give rise to functions

V_{F} (θ)

and

V_{B} (θ)

, expressions in the Appendix, such that their reciprocals are Lyapunov densities for forward and backward flow, respectively, of the system

\dot{θ} = F (θ)

. This implies that almost all solutions reach

V_{B}^{- 1} (0)

in backward time (hence it is a repellor) and reach

V_{A}^{- 1} (0)

in the forward time (hence it is an attractor), refer Figure (Figure 4) (right). Notice that almost all trajectories get trapped in the region

{θ \in {[- π, π]}^{3} : - 1.2 < θ_{1} < 0}

. This implies that oscillators

φ_{1}

and

φ_{3}

eventually satisfy

| φ_{1} - φ_{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} (0)

is the line

θ_{i} = θ_{0}

, then almost all trajectories of the phase-difference system satisfy

θ_{i} (t) \to θ_{0}

, which eventually implies that

(φ_{i} (t) - φ_{3} (t)) \to θ_{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{φ}}_{j} = ω_{j} + \sum_{\begin{matrix} i, k = 1 \\ i, k \neq j \end{matrix}}^{d} g (φ_{i} + φ_{k} - 2 φ_{j}), φ_{i} \in T, i = 1, \dots, 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.

Appendix A

The expression of

V_{F} (θ)

and

V_{B} (θ)

in Example 4 is given by

\begin{matrix} V_{F} (θ) & = & - 0.00108589 sin (t_{1} - 2 t_{2}) - 0.591088 sin (t_{1} - t_{2}) + 0.00894074 sin (2 (t_{1} - t_{2})) + 0.07452 sin (2 t_{1} - t_{2}) \\ + 0.215346 sin (t_{1} + t_{2}) + 0.00185937 sin (2 (t_{1} + t_{2})) - 0.00278604 sin (2 t_{1} + t_{2}) - 0.0138828 sin (t_{1} + 2 t_{2}) \\ + 0.0494588 cos (t_{1} - 2 t_{2}) - 0.226064 cos (t_{1} - t_{2}) - 0.0107154 cos (2 (t_{1} - t_{2})) + 0.133784 cos (2 t_{1} - t_{2}) \\ + 0.290884 cos (t_{1} + t_{2}) + 0.000185757 cos (2 (t_{1} + t_{2})) - 0.0583589 cos (2 t_{1} + t_{2}) - 0.421386 sin (2 t_{1}) \\ + 0.00646918 cos (t_{1} + 2 t_{2}) + 1.30953 sin (t_{1}) + 0.0666633 cos (2 t_{1}) - 1.54687 cos (t_{1}) \\ - 0.74803 sin (t_{2}) - 0.232342 cos (t_{2}) - 0.036627 cos (2 t_{2}) + 0.0472706 sin (t_{2}) cos (t_{2}) + 1.60567 \\ V_{B} (θ) & = & 0.591015 sin (t_{1} - t_{2}) - 0.00893979 sin (2 (t_{1} - t_{2})) + 0.0745105 sin (2 t_{1} - t_{2}) - 0.215319 sin (t_{1} + t_{2}) \\ - 0.00185914 sin (2 (t_{1} + t_{2})) - 0.00278548 sin (2 t_{1} + t_{2}) - 0.0138811 sin (t_{1} + 2 t_{2}) - 0.00108587 sin (t_{1} - 2 t_{2}) \\ - 0.226037 cos (t_{1} - t_{2}) - 0.0107141 cos (2 (t_{1} - t_{2})) - 0.133769 cos (2 t_{1} - t_{2}) + 0.290849 cos (t_{1} + t_{2}) \\ + 0.000185739 cos (2 (t_{1} + t_{2})) + 0.0583518 cos (2 t_{1} + t_{2}) - 0.00646834 cos (t_{1} + 2 t_{2}) - 0.0494526 cos (t_{1} - 2 t_{2}) \\ + 1.30938 sin (t_{1}) + 0.421337 sin (2 t_{1}) + 1.54668 cos (t_{1}) + 0.0666534 cos (2 t_{1}) - 0.747938 sin (t_{2}) \\ + 0.232313 cos (t_{2}) - 0.0366224 cos (2 t_{2}) - 0.02863235 sin (2 t_{2}) + 1.60547 \end{matrix}

References

Tedrake, R.; Manchester, I.R.; Tobenkin, M.; Roberts, J.W. LQR-trees: Feedback motion planning via sums-of-squares verification. The International Journal of Robotics Research 2010, 29, 1038–1052. [CrossRef]
Moore, J.; Tedrake, R. Control synthesis and verification for a perching UAV using LQR-Trees. 2012 IEEE 51st IEEE Conference on Decision and Control (CDC). IEEE, 2012, pp. 3707–3714. [CrossRef]
Yoon, S.; Park, S.; Kim, Y. Circular motion guidance law for coordinated standoff tracking of a moving target. IEEE Transactions on Aerospace and Electronic Systems 2013, 49, 2440–2462.
Liu, Z.; Zhang, Y.; Yu, X.; Yuan, C. Unmanned surface vehicles: An overview of developments and challenges. Annual Reviews in Control 2016, 41, 71–93. [CrossRef]
Misra, R.; Wisniewski, R.; Karabacak, Ö. Sum-of-Squares based computation of a Lyapunov function for proving stability of a satellite with electromagnetic actuation. IFAC-PapersOnLine 2020, 53, 7380–7385. [CrossRef]
Angeli, D. An almost global notion of input-to-state stability. IEEE Transactions on Automatic Control 2004, 49, 866–874. [CrossRef]
Koditschek, D.E.; others. The application of total energy as a Lyapunov function for mechanical control systems. Dynamics and control of multibody systems 1989, 97, 131–157. [CrossRef]
Rantzer, A. A dual to Lyapunov’s stability theorem. Systems & Control Letters 2001, 42, 161–168. [CrossRef]
Prajna, S.; Parrilo, P.A.; Rantzer, A. Nonlinear control synthesis by convex optimization. IEEE Transactions on Automatic Control 2004, 49, 310–314. [CrossRef]
Henrion, D.; Garulli, A. Positive polynomials in control; Vol. 312, Springer Science & Business Media, 2005.
Dumitrescu, B. Positive trigonometric polynomials and signal processing applications; Vol. 103, Springer, 2007.
Alkire, B.; Vandenberghe, L. Convex optimization problems involving finite autocorrelation sequences. Mathematical Programming 2002, 93, 331–359. [CrossRef]
Davidson, T.N.; Luo, Z.Q.; Sturm, J.F. Linear matrix inequality formulation of spectral mask constraints with applications to FIR filter design. IEEE Transactions on Signal Processing 2002, 50, 2702–2715. [CrossRef]
Zhang, J.H.; Janschek, K.; Böhme, J.; Zeng, Y.J. Multi-resolution dyadic wavelet denoising approach for extraction of visual evoked potentials in the brain. IEE Proceedings-Vision, Image and Signal Processing 2004, 151, 180–186. [CrossRef]
Niemistö, R.; Dumitrescu, B.; Tăbuş, I. SDP design procedures for near-optimum IIR compaction filters. Signal processing 2002, 82, 911–924. [CrossRef]
Henrion, D.; Vyhlídal, T. Positive trigonometric polynomials for strong stability of difference equations. IFAC Proceedings Volumes 2011, 44, 296–301. [CrossRef]
Mayer, T. Polyhedral faces in Gram spectrahedra of binary forms. Linear Algebra and its Applications 2021, 608, 133–157. [CrossRef]
McLean, J.W.; Woerdeman, H.J. Spectral factorizations and sums of squares representations via semidefinite programming. SIAM journal on matrix analysis and applications 2002, 23, 646–655. [CrossRef]
Lee, J. Introduction to Smooth Manifolds; Graduate Texts in Mathematics, Springer New York, 2012.
Kudeyt, M.; Kıvılcım, A.; Köksal-Ersöz, E.; İlhan, F.; Karabacak, Ö. Certification of almost global phase synchronization of all-to-all coupled phase oscillators. Chaos, Solitons & Fractals 2023, 174, 113838. [CrossRef]
Karabacak, Ö.; Wisniewski, R.; Leth, J. On the almost global stability of invariant sets. 2018 European Control Conference (ECC). IEEE, 2018, pp. 1648–1653. [CrossRef]
Ha, S.Y.; Ryoo, S.Y. Asymptotic phase-locking dynamics and critical coupling strength for the Kuramoto model. Communications in Mathematical Physics 2020, 377, 811–857. [CrossRef]
Chopra, N.; Spong, M.W. On Exponential Synchronization of Kuramoto Oscillators. IEEE Transactions on Automatic Control 2009, 54, 353–357. doi:10.1109/TAC.2008.2007884. [CrossRef]
Ha, S.Y.; Ha, T.; Kim, J.H. On the complete synchronization of the Kuramoto phase model. Physica D: Nonlinear Phenomena 2010, 239, 1692–1700. [CrossRef]
Ashwin, P.; Burylko, O.; Maistrenko, Y. Bifurcation to heteroclinic cycles and sensitivity in three and four coupled phase oscillators. Physica D: Nonlinear Phenomena 2008, 237, 454–466. [CrossRef]
Hansel, D.; Mato, G.; Meunier, C. Clustering and slow switching in globally coupled phase oscillators. Physical Review E 1993, 48, 3470. [CrossRef]
Canale, E.; Monzón, P. Almost global synchronization of symmetric Kuramoto coupled oscillators. Systems Structure and Control 2008, 8, 167–190. [CrossRef]
Qin, Y.; Kawano, Y.; Portoles, O.; Cao, M. Partial Phase Cohesiveness in Networks of Networks of Kuramoto Oscillators. IEEE Transactions on Automatic Control 2021, 66, 6100–6107. Conference Name: IEEE Transactions on Automatic Control, doi:10.1109/TAC.2021.3062005. [CrossRef]

Figure 1. (Left)

Arc (a)

for

a = π / 2

and

d = 3

is a solid hexagonal cylinder as discussed in Remark 2; and (Right)

D_{ε}

appearing in the hypothesis of Corollary 3 for

d = 3

and

ε = 0.6

Figure 1. (Left)

Arc (a)

for

a = π / 2

and

d = 3

is a solid hexagonal cylinder as discussed in Remark 2; and (Right)

D_{ε}

appearing in the hypothesis of Corollary 3 for

d = 3

and

ε = 0.6

Figure 2. Parameters

(β_{1}, α_{2})

for certifying almost global stability of (44) obtain by employing Theorem 1 and Gram matrices. Each blue star depicts a pair of parameter

(β_{1}, α_{2})

for which (44) is almost global phase synchronized. For each such pair, two positive semidefinite matrices

G_{V}

and

G_{W}

satisfying (17) exist and hence

ρ (θ) = 1 / V (θ)

is a Lyapunov density due to Corollary 1 and Theorem 1.

Figure 2. Parameters

(β_{1}, α_{2})

for certifying almost global stability of (44) obtain by employing Theorem 1 and Gram matrices. Each blue star depicts a pair of parameter

(β_{1}, α_{2})

for which (44) is almost global phase synchronized. For each such pair, two positive semidefinite matrices

G_{V}

and

G_{W}

satisfying (17) exist and hence

ρ (θ) = 1 / V (θ)

is a Lyapunov density due to Corollary 1 and Theorem 1.

Figure 3. (Left) flow and equilibrium points (red dots) of the phase-difference system (40)–(41); and several regions discussed in Example 2. The reduced vector field and the regions are obtained from their three-dimensional counterparts (right) on quotient-ing by the equivalence relation ∼ defined by

φ^{(1)} \sim φ^{(2)}

if and only if

φ^{(1)} - φ^{(2)} = ε \cdot 1

for some

ε \in R

φ^{(1)} \sim φ^{(2)}

if and only if

φ^{(1)} - φ^{(2)} = ε \cdot 1

for some

ε \in R

Figure 4. (Left) The plot of Lyapunov function

V_{*}

over

A r c (π / 2)

in Example 2 implying that the all-to-all coupled uniform oscillators (38) for coupling (39) with

β_{1} = π / 6

and

α_{2} = 0.3889

exhibit local synchronization. (Right) Verification that almost all solutions emerge from the repeller

V_{B}^{- 1} (0)

and reach the attractor

V_{F}^{- 1} (0)

for the phase-difference system in Example 4.

Figure 4. (Left) The plot of Lyapunov function

V_{*}

over

A r c (π / 2)

in Example 2 implying that the all-to-all coupled uniform oscillators (38) for coupling (39) with

β_{1} = π / 6

and

α_{2} = 0.3889

exhibit local synchronization. (Right) Verification that almost all solutions emerge from the repeller

V_{B}^{- 1} (0)

and reach the attractor

V_{F}^{- 1} (0)

for the phase-difference system in Example 4.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

Downloads

Views

Comments

Subscription

Notify me about updates to this article or when a peer-reviewed version is published.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Construction of Lyapunov Certificates for Systems on Torus Using Trigonometric Polynomials

Abstract

Keywords:

Subject:

1. Introduction

2. Preliminaries

2.1. Trigonometric Polynomials

2.2. A Class of Coupled Oscillators

2.2.1. Generalized Kuramoto Models with Coupling Topology

3. Trigonometric Polynomials in Stability Analysis

4. Examples

5. Conclusions

Appendix A

References

MDPI Initiatives

Important Links

Subscribe