Synchronization of
Kuramoto Oscillators:
Inverse Taylor
Expansions^{†}^{†}thanks: This work was supported in part by the
U.S. Department of Energy (DOE) Solar Energy Technologies Office
under Contract No. DEEE00001583.
Abstract
Synchronization in networks of coupled oscillators is a widely studied topic with extensive scientific and engineering applications. In this paper, we study the frequency synchronization problem for networks of Kuramoto oscillators with arbitrary topology and heterogeneous edge weights. We propose a novel equivalent transcription for the equilibrium synchronization equation. Using this transcription, we develop a power series expansion to compute the synchronized solution of the Kuramoto model as well as a sufficient condition for the strong convergence of this series expansion. Truncating the power series provides (i) an efficient approximation scheme for computing the synchronized solution, and (ii) a simpletocheck, statisticallycorrect hierarchy of increasingly accurate synchronization tests. This hierarchy of tests provides a theoretical foundation for and generalizes the bestknown approximate synchronization test in the literature. Our numerical experiments illustrate the accuracy and the computational efficiency of the truncated series approximation compared to existing iterative methods and existing synchronization tests.
Key words. Kuramoto oscillators, frequency synchronization, synchronization manifold, Taylor series, power networks
AMS subject classifications. 34D06, 34C15, 93D20, 37C25, 37M20
1 Introduction
Collective synchronization is an interesting behaviour which lies at the heart of various natural phenomena. The celebrated Kuramoto model [22] is one of the simplest models for studying synchronization in a network of coupled oscillators. Kuramoto model has been successfully used to model the synchronization behaviour of a wide range of physical, chemical, and biological system [2]. Examples include the power grids [15, 10], automated vehicle coordination [21, 31], pacemakers in heart [42], clock synchronization [32], and neural networks [13]; see also [6, Chapter 13] for additional examples. One of the most interesting types of synchronization is frequency synchronization, where all oscillators reach the same rotational frequency with possibly different phases. It is wellknown that the Kuramoto model can exhibit a transition from incoherence to frequency synchronization. For many applications, such as power networks, it is important to have an accurate estimate of this transition to synchronization. This is essential, particularly as the grid is being pushed closer to its maximum capacity due to increases in the load demand and penetration of renewable energy units. Finding sharp conditions to determine when this transition happens continues to be a challenging as well as a critical problem.
Literature Review
The problem of finding conditions for existence of a stable synchronized solution for the Kuramoto model of coupled oscillators has been studied extensively in the literature. For complete graphs with homogeneous weights, the order parameter is used to implicitly determine the exact critical coupling needed for a synchronized solution [3, 25, 39]. For acyclic graphs with heterogeneous weights, a necessary and sufficient condition is developed for synchronization of the heterogeneous Kuramoto model [12]. In addition, Lyapunov analysis applied to the complete graph is used in [8] to give a sufficient condition and in [9] to give an explicit necessary and sufficient condition for existence of a synchronized solution. However for general topology graphs, such a complete characterization of frequency synchronization does not exist. For general graphs with heterogeneous weights, several necessary conditions and sufficient conditions for existence of stable synchronized solution have been reported in the literature. [36] requires sufficiently large nodal degrees relative to the natural frequencies, [4] uses the cutset in the graph, and [18] states that the algebraic connectivity must be sufficiently large compared to the difference in natural frequencies of connected oscillators. Recently, a novel cutset projection operator has been introduced to rigorously prove a simpletocheck, sufficient condition for synchronization of Kuramoto model [19]. Using numerous simulations, it is shown that this new sufficient condition scales better to large networks [19].
Despite these deep results in the literature, the existing synchronization conditions usually provide conservative estimates for the synchronization threshold. In an effort to come closer to finding the exact synchronization threshold, [12] and [17] introduce a statistically accurate approximate test for synchronization that depends on the network parameters and topology derived from the linearized Kuramoto map and the converging power series expansion of the phase angles of the Kuramoto oscillators, respectively.
If existence of a frequency synchronized solution can be guaranteed, then the next step is computing the synchronized solutions. A common method to approximate the solution is to linearize the equations. This will result in studying the equations of the form , where is the Laplacian matrix of the network [23, 33, 35]. The angles can be approximately solved very efficiently, even for extremely large, sparse graphs [40, Theorem 3.1]. However, when phase differences of the oscillators are large, this linear approximation is not very accurate.
In order to compute the synchronization manifold of the nonlinear Kuramoto equations, one can employ iterative numerical algorithms such as Newton–Raphson or Gauss–Seidel [34, 37, 14]. Unfortunately, these algorithms do not guarantee convergence to the synchronized solutions and failure of these algorithms could be due to numerical instability, an initialization issue, or nonexistence of the solution. Another approach is to use numerical polynomial homotopy continuation (NPHC). It is guaranteed that NPHC will find all stable and unstable manifolds of the Kuramoto model, but this method is not computationally tractable for large networks; [24] uses NPHC to study the homogeneous Kuramoto model for particular graph topologies with up to nodes. Finally, [41] gives an approximate analytical solution for stable synchronization manifolds using the order parameter. However, this approximation scheme is only applicable to the uniformweight Kuramoto model with alltoall connections.
A wide range of methods for finding synchronized solutions of the Kuramoto model stem from the power network literature, where different techniques are used to find the solutions of the AC power flow equations. Here, we only review two of these approaches. The first approach is called Holomorphic Embedding LoadFlow Method (HELM) and has been proposed to find all the solutions of power flow equations [38]. While HELM is based on advanced results and concepts from complex analysis, its numerical implementation is recursive and straightforward [38, 30]. However, HELM is reported to be much slower than the Newton–Raphson methods [30]. The second approach is the optimization approach, whereas an optimal power flow problem (OPF) is used to solve for the AC power flow equations. The OPF problems have been studied extensively in the power network literature, e.g., see [27, 28, 26]. Thus, one can use the numerical algorithms for the optimization problem to find the synchronized solution of the Kuramoto model. Unfortunately, due to the nonconvex nature of the OPFs, these algorithms usually result in an approximation of the synchronized solution.
Contribution
The contributions of this paper are both theoretical and computational. From a theoretical viewpoint, first, we review important properties of the Kuramoto model of coupled oscillators and, as a minor contribution, we provide a rigorous proof for the following wellknown folk theorem: frequency synchronization is equivalent with the existence of a stable synchronization manifold (see [25] and [39] for statement of this result without proof). Second, by introducing the notions of edge vectors and flow vectors in graphs, we propose four equivalent transcriptions for the synchronization manifold of the Kuramoto model: node, flow, constrained edge, and unconstrained edge balance equations. While the first three formulations have already been studied in the literature (see [12] and [19]), the unconstrained edge balance equations provide a novel important characterization of the synchronization manifold. Our main technical results are (1) a sufficient condition for existence of a unique solution for unconstrained edge balance equations and (2) a recursive expression for each term of the Taylor series expansion for this solution of the unconstrained edge balance equations. Additionally, we prove that, if our simpletocheck sufficient condition is satisfied, then the Taylor series expansion for the solution converges strongly. We also provide an algorithm to symbolically compute all terms of the expansion. Third and final, using the onetoone correspondence between solutions of the unconstrained edge balance equations and synchronized solutions of the Kuramoto model, we propose a power series expansion for the synchronized solutions of the Kuramoto model and an estimate on the region of convergence of the power series.
From a computational viewpoint, first, we propose a method to approximate the synchronization manifold of the Kuramoto model using the truncated power series. We present several numerical experiments using IEEE test cases and random graphs to illustrate (1) the accuracy of the truncated series and (2) the computational efficiency of the new methods for computing the synchronization manifold. We show that the seventh order approximate method has low absolute error when applied to IEEE test cases with weakly coupled oscillators. The truncated series, up to the seventh order, have comparable computational efficiency to Newton–Raphson when solving for solutions with static graph topology and multiple natural frequencies, or power injections. Second, based on our novel power series approach, we propose a hierarchy of approximate tests for synchronization of the Kuramoto model; our approach provides a theoretical basis for and generalize the stateoftheart approximate synchronization test in the literature [12]. With numerical analysis, we verify the accuracy of our family of approximate tests for several random graphs and numerous IEEE test cases. In each of these cases, we show that our new approximate tests are a significant improvement compared to the bestknown approximate condition given in [12].
Finally, we compare this paper with our preliminary conference article [17]. In short, this paper presents a substantially more complete and comprehensive treatment of the power series approach to synchronization of Kuramoto oscillators. Specifically, while [17] presents a power series expansion for nodal phase angles, this paper develops a novel power series expansion for the flows in the network. Using the Banach FixedPoint Theorem, we provide an estimate on the domain of convergence of the power series which is substantially larger than the estimates given in [17]. Moreover, our numerical analysis shows that the hierarchy of approximate synchronization tests obtained by truncating this power series is more accurate than the estimate tests proposed in [17].
Paper organization
In Section LABEL:sec:prelim, we give preliminaries and notation used in the paper. In Section LABEL:sec:kuramotoLABEL:sec:auxiliary we review the Kuramoto model, frequency synchronization, and give several equivalent formulations of the algebraic Kuramoto equation. Sections LABEL:sec:cutflowbalancesolution and LABEL:sec:syncOfKuramoto contain the paper’s main theoretical results and a family of approximate synchronization tests. Finally, Section LABEL:sec:numerical contains numerical experiments analyzing the approximate synchronization tests and efficiency of computation methods for the synchronization manifold.
2 Preliminaries and notation
Vectors and functions
Let , , and denote the set of nonnegative integers, the dimensional real Euclidean space, and the dimensional complex Euclidean space, respectively. For , let denote the double factorial. For and , the real polydisk with center and radius is
Similarly, for and , the complex polydisk with center and radius is
Let and be dimensional column vectors of ones and zeros respectively. For , let and be the diagonal matrix with , for every . For with , let , where
For every , we denote the torus by . For every , the clockwise rotation of by the angle is the function defined by
Using the rotation function, one defines an equivalence relation on the torus as follows: For every two points , we say if there exists such that . For every , the equivalence class of is denoted by . The quotient space of under the equivalence relation is denoted by .
Algebraic graph theory
Let be a weighted undirected connected graph with the node set and the edge set with elements. We assume that has no selfloops and the weights of the edges are described by the nonnegative, symmetric adjacency matrix . The Laplacian matrix of the graph is . Define the diagonal edge weight matrix by . It is known that the Laplacian is . Since is singular, we use the Moore–Penrose pseudoinverse which has the following properties: , , , and . In addition, for a connected graph . The weighted cutset projection matrix is the oblique projection onto parallel to given by
The weighted cutset projection matrix is idempotent, and and are its eigenvalues with algebraic (and geometric) multiplicity and , respectively. Additional properties of are in [19, Theorem 5]. Similarly, the weighted cycle projection matrix is the oblique projection onto parallel to given by
Analytic functions and power series
A multiindex is a member of . For every , we define . For , the formal expression
(1) 
where , for every is called a formal power series around point . The power series converges strongly at point if all rearrangement of the terms of the series converges. For every , the domain of convergence of (LABEL:eq:5) around is defined as the set of all points such that the power series converges strongly at point . While for , one can show that the domain of convergence is an open interval around , for the domain of convergence of a power series is not necessarily an open polydisk around . An open set is a Reinhardt domain if, for every and every , we have . The Reinhardt domains can be considered as the generalization of the disks on the complex plane to higher dimensions.
3 The heterogeneous Kuramoto model
The Kuramoto model is a system of oscillators, where each oscillator has a natural frequency and its state is represented by a phase angle . The interconnection of these oscillators are described using a weighted undirected connected graph , with nodes , edges , and positive weights . The dynamics for the heterogeneous Kuramoto model is given by:
(2) 
In matrix language, one can write this differential equations as:
(3) 
where is the phase vector, is the natural frequency vector, and is the incidence matrix for the graph . One can show that if is a solution for the Kuramoto model (LABEL:eq:kuramoto_model) then, for every , the curve is also a solution of (LABEL:eq:kuramoto_model). Therefore, for the rest of this paper, we consider the state space of the Kuramoto model (LABEL:eq:kuramoto_model) to be .
Definition 1 (Frequency synchronization).
A solution of the coupled oscillator model (LABEL:eq:kuramoto_model) achieves frequency synchronization if there exists a frequency such that
By summing all the equations in (LABEL:eq:2), one can show that if a solution of (LABEL:eq:kuramoto_model) achieves frequency synchronization then . Therefore, without loss of generality, we can assume that in the Kuramoto model (LABEL:eq:kuramoto_model), we have and .
Definition 2 (Synchronization manifold).
Let be a solution of the algebraic equation
(4) 
Then is called a synchronization manifold for the Kuramoto model (LABEL:eq:kuramoto_model).
The following theorem reduces the problem of local frequency synchronization in the Kuramoto model (LABEL:eq:kuramoto_model) to the existence of a solution for the algebraic equations (LABEL:eq:synchronization_manifold).
Theorem 3 (Characterization of frequency synchronization).
For the heterogeneous Kuramoto model (LABEL:eq:kuramoto_model) on graph , the following statements are equivalent:

there exists an open set such that every solution of the Kuramoto model (LABEL:eq:kuramoto_model) starting in set achieves frequency synchronization;

there exists a locally asymptotically stable synchronization manifold for (LABEL:eq:kuramoto_model).
Additionally, if any of equivalent conditions LABEL:p1:frequency or LABEL:p2:equilibrium holds, then, for every , we have .
Proof.
Regarding , if the solution achieves frequency synchronization, then for all . Consider a sequence of natural numbers and the corresponding sequence in . Since is a compact metric space, it is sequentially compact [29, Theorem 28.2]. This means that there is a subsequence such that is convergent. Then exists and . Therefore is a synchronization manifold because it is a solution for equation (LABEL:eq:synchronization_manifold) and is locally asymptotically stable since all solutions starting in reach .
Regarding , by the definition of local asymptotic stability, there exists some such that the open set is defined to be where is the synchronization manifold. Then for solutions starting in , so is also a frequency synchronized solution for (LABEL:eq:kuramoto_model).
The last statement follows from the proofs of and .
In many application, such as power networks, not only is it important to study the frequency synchronization of the Kuramoto oscillators but also it is essential to bound the position of the synchronization manifold due to some security constraints for the grid. An important class of security constraints are thermal constraints which are usually expressed as bounds on the geodesic distances , for . The geodesic distance is defined as the minimum of the clockwise and counterclockwise arc lengths between the phase angles . Let be an undirected weighted connected graph with edge set and let We define the cohesive subset by
For every , we define the embedded cohesive subset by:
where . Note that, in general, we have . We refer to [19] for additional properties of embedded cohesive subset. In particular, it is shown that is diffeomorphic with , for every [19, Theorem 8]. Using this result, in the rest of this paper we identify the set with .
4 Equivalent transcriptions of the equilibrium manifold
Consider an undirected graph with vertex set and edge set with . We start by introducing three vector spaces defined by :

the node space is ; elements of this space are called node vectors;

the edge space is ; elements of this space are called edge vectors; and

the flow vector space is ; elements of belonging to this space are called by flow vectors.
It is easy to see that an edge vector is a flow vector if and only if there exists a node vector such that .
Next, we introduce four different balance equations on an undirected graph with incidence matrix , weight matrix , cutset projection , and cycle projection . Given a node vector , define the shorthand flow vector . The node balance equations in the unknown node vector is
(5) 
The flow balance equations in the unknown flow vector is
(6) 
The constrained edge balance equations in the unknown edge vector is
(7) 
The unconstrained edge balance equations in the unknown edge vector is
(8) 
We now present equivalent characterizations for synchronization manifold of the Kuramoto model (LABEL:eq:kuramoto_model).
Theorem 4 (Characterization of synchronization manifold).
Consider an undirected connected graph with incidence matrix , weight matrix , cutset projection , and cycle projection . Given a node vector , define the shorthand . Pick an angle . Then the following statements are equivalent:

there exists a unique locally exponentially stable synchronization manifold for the Kuramoto model (LABEL:eq:kuramoto_model) in ;

the node balance equations (LABEL:eq:nodal_sync) have a unique solution in ;

the flow balance equations (LABEL:eq:edge_sync) have a unique solution with ;

the constrained edge balance equations (LABEL:eq:aux_sync) have a unique solution with ;

the unconstrained edge balance equations (LABEL:eq:cutflow) have a unique solution with .
Moreover, if one of the above equivalent conditions hold, then
Proof.
The implications and are easy to show.
Regarding , if is the unique solution to the flow balance equations (LABEL:eq:edge_sync), then satisfies and is a solution for the edge balance equations (LABEL:eq:aux_sync). Now, we show that is the unique solution for the constrained edge balance equations (LABEL:eq:aux_sync) such that . Suppose that is another solution of the constrained edge balance equations (LABEL:eq:aux_sync) satisfying . Then, by the constrained edge balance equations (LABEL:eq:aux_sync), there exists such that and . This implies that and . Therefore, and satisfies the flow balance equations (LABEL:eq:edge_sync). However, this is in contradiction with the facts that and that is the unique solution of the flow balance equations (LABEL:eq:edge_sync).
Regarding , if is a solution of constrained edge balance equations (LABEL:eq:aux_sync) satisfying , then
(9)  
(10) 
Because , the inclusion (LABEL:eq:second_aux) implies that
(11) 
By adding equations (LABEL:eq:first_aux) and (LABEL:eq:third_aux), we obtain . This means that satisfies unconstrained edge balance equations (LABEL:eq:cutflow).
Regarding , if solves the unconstrained edge balance equations (LABEL:eq:cutflow), then
(12) 
Leftmultiplying both sides of equations (LABEL:eq:fourth) by and using the facts that , and , we obtain
Leftmultiplying both side of the equations (LABEL:eq:fourth) by we obtain
This last equality implies that . Thus, there exists a vector such that . First, note that and, by multiplying both side of this equation by , we obtain
Moreover, . Thus, we have and . This implies that is a synchronization manifold for the Kuramoto model (LABEL:eq:kuramoto_model) in . The uniqueness follows from [19, Theorem 10, statement (ii)].
5 Solvability of the unconstrained edge balance equations
The unconstrained edge balance equations (LABEL:eq:cutflow) allow us to focus on a single analytic map whose inverse can be used in computing the synchronization solutions of Kuramoto model. In this section, we study the solvability of these equations and find their inverse on a suitable domain. We start with relaxing the condition and complexifing the equations (LABEL:eq:cutflow). This extension will allow us to use the theory of several complex variables to find the Taylor series expansion for the inversion of the complexified equations and prove the strong convergence of the Taylor series. We then restrict back to real domain and use the constrained to find the solutions of the unconstrained edge balance equations (LABEL:eq:cutflow). We start with some useful definitions. Given an undirected graph with cutset projection and cycle projection , define the complex edge balance map by
and the real edge balance map by
With this notation, the unconstrained edge balance equations (LABEL:eq:cutflow) read , together with the constraints .
Next, we define the scalar function by:
The graph of function on the interval is shown in Figure (LABEL:fig:functionh).
Since is continuous and strictly monotonicallydecreasing, its inverse exists and is denoted by . Although we do not have an analytical form for , it is simple to compute numerically.
We are now ready to provide an estimate on the image of the maps and and to present a power series expansion for the inverse maps and on suitable domains.
Theorem 5 (Properties of the complex edge balance map).
Consider an undirected connected graph with cutset projection and cycle projection . Select such that and define by
Then the following statements holds:

there exists a unique such that ; that is unconstrained edge balance equations have a unique solution;

there exists a holomorphic map such that
that is the edge balance map is invertible on ;

the power series
converges strongly to , where, for every , the term is a homogeneous polynomial of order in defined iteratively by:
Proof.
Regarding statement LABEL:p1:existenceG, define the map by
The map appears from bringing all terms of the unconstrained edge equation (LABEL:eq:cutflow) to the left hand side and adding to both sides.
First, we show that . For , we compute
Moreover, for , we have . These equalities imply that
(13) 
where, for the last inequality, we used the fact that . By the definition of , we have
Noting the fact that , we obtain
Now, by replacing the above equation into inequality (LABEL:eq:important), we have
Thus, by the Banach FixedPoint Theorem, there exists a unique fixed point for . By construction, this fixed element satisfies
This completes the proof of statement LABEL:p1:existenceG.
Regarding statement LABEL:p3:inverseG, by statement LABEL:p1:existenceG, for every such that , there exists a unique such that . This implies that has a unique inverse which satisfies the equalities in statement LABEL:p3:inverseG. Now we show that is holomorphic on . Note that, for every , the derivative of the map at point is given by:
We first show that is invertible. Suppose that, there exists such that . This means that and . The first equality implies that and the second inequality implies that . Therefore, there exists such that . Thus, we get
(14) 
Moreover, is a diagonal matrix with positive diagonal elements. Therefore, equations (LABEL:eq:laplacian) implies that and as a result . This proves that the derivative is invertible. Now, by the Inverse Function Theorem [1, Theorem 2.5.2], the maps and are locally holomorphic and therefore they are holomorphic on their domains. This completes the proof of statement LABEL:p3:inverseG.
Regarding statement LABEL:p4:power_seriesG, we first find the formal power series representation for . Suppose that is the formal power series for . Then we have
By replacing the power series for and using the power series expansion of , we obtain
(15) 
By equating the same order terms on the both side of equation (LABEL:eq:series_equality) and using the fact that , we obtain that and , for every . Simple bookkeeping shows that the recursive formula in statement LABEL:p4:power_seriesG holds for the odd terms in the power series.
Finally, we prove that the formal power series converges on the domain . Note that statement LABEL:p3:inverseG implies that the map is holomorphic and and that the set is a Reinhardt domain. Therefore, [16, Theorem 2.4.5] implies that the power series converges strongly on the domain .
Remark 1 (Properties of the real edge balance map).
A similar result as Theorem LABEL:thm:propertyGcomplex holds for the real edge balance map by replacing the complex variables by their real counterparts. The proof is straightforward by restricting the results in Theorem LABEL:thm:propertyGcomplex to the real Euclidean space.
6 Inverse Taylor expansion for Kuramoto model
In this section we study the synchronization of the Kuramoto model (LABEL:eq:kuramoto_model) by applying the results on the unconstrained edge balance equations (LABEL:eq:cutflow) from Theorem LABEL:thm:propertyGcomplex in the previous section.
Theorem 6 (Inverse Taylor expansion).
Consider the Kuramoto model (LABEL:eq:kuramoto_model) with undirected connected graph , weighted cutset projection , and weighted cycle projection . Given frequencies satisfying
(T0) 
define by
Then the following statements hold:

there exists a unique locally stable synchronization manifold in ; and

the power series
(16) converges strongly to where, for every , the term is a homogeneous polynomial of order in defined iteratively as in Theorem LABEL:thm:propertyGcomplexLABEL:p4:power_seriesG.
This theorem is an immediate application of Theorem LABEL:thm:existence_uniqueness_S on the equivalent transcriptions and of Theorem LABEL:thm:propertyGcomplex on the properties of the maps and .
Proof (Proof of Theorem LABEL:thm:suff_condition+power_series).
Regarding statement LABEL:p1:ex, we use (the real version of) Theorem LABEL:thm:propertyGcomplexLABEL:p1:existenceG with . Since