Synchronization in complex oscillator networks and smart grids Florian dorflera b, 1. Michael chertkoy and francesco bulloa aCenter for Control, Dynamical Systems, and Compu University of California, Santa Barbara, CA 93106; andCenter for Nonlinear Studies and Theory Division, Los Alamos National Laboratory, Los Alamos, NM 87545 Edited by Steven H. Strogatz, Cornell University, Ithaca, NY, and accepted by the Editorial Board November 14, 2012(received for review July 16, 2012) The emergence of synchronization in a network of coupled oscil- oscillators Vi with Newtonian dynamics, inertia coefficients Mi lators is a fascinating topic in various scientific disciplines. A widely and viscous damping D. The remaining oscillators v2 feature adopted model of a coupled oscillator network is characterized by first-order dynamics with time constants D. A perfect electrical a population of heterogeneous phase oscillators, a graph describ- analog of the coupled oscillator model [1] is given by the classic ing the interaction among them, and diffusive and sinusoidal structure-preserving power network model (4), our er coupling. It is known that a strongly coupled and sufficiently application of interest. Here, the first- and second-orde homogeneous network synchronizes, but the exact threshold from namics correspond to loads and generators, respectively, incoherence to synchrony is unknown. Here, we present a unique, right-hand sides depict the power injections @; and the concise, and closed-form condition for synchronization of the fully flows ai sin(0-e1) along transmission lines nonlinear, nonequilibrium, and dynamic network, Our synchro The rich dynamic behavior of the coupled oscillator model zation condition can be stated elegantly in terms of the network 1] arises from a competition between each oscillator's tendency topology and parameters or equivalently in terms of an intuitive, to align with its natural frequency o; and the synchronization linear, and static auxiliary system. Our results significantly improve enforcing coupling ay sin(; -e) with its neighbors. If all natural upon the existing conditions advocated thus far, they are provably frequencies o are identical, the coupled oscillator dynamics[ 1] exact for various interesting network topologies and parameters; collapse to a trivial phase-synchronized equilibrium, where all and biology as well as in engineered oscillator networks, such as aligned equilibrium. Moreover, even if the coupled oscillator electrical power networks We illustrate the validity, the accuracy, model [1] synchronizes, the motion of its center of mass still and the practical applicability of our results in complex network carries the fiux of angular rotation, respectively, the flux of elec- scenarios and in smart grid applications. trical power from generators to loads in a power grid. Despite all these complications, the main result of this article is that, for a broad range of network topologies and parameters, an elegant and easily verified criterion characterizes synchronization of the the scientific interest in the synchronization of coupled nonlinear and nonequilibrium dynamic oscillator network [1] oscillators can be traced back to Christiaan Huygens' seminal Review of Synchronization in Oscillator Networks work on"an odd kind sympathy"between coupled pendulum The coupled oscillator model [1] unifies various models in the locks(1),and it continues to fascinate the scientifc community literature, including dynamic models of electrical power net to date(2, 3). A mechanical analog of a coupled oscillator net- works. Modeling of electrical power networks is discussed in SI work is shown in Fig. 1A and consists of a group of particles Text in detail. For v2= 0, the coupled oscillator model [1] without colliding. Each particle is characterized by a phase angle havior(5), populations of flashing fireflies(6), and crowd syn- e, and has a preferred natural rotation frequency o. Pairs of chrony on London's Millennium bridge(7), as well as in Huygen,'s interacting particles i and j are coupled through an elastic spring pendulum clocks(8). For V1=9, the coupled oscillator model with stiffness air. Intuitively, a weakly pled oscillator net- (1)reduces to the celebrated Kuramoto model(9), which appears work with strongly, heterogeneous natural frequencies a does in coupled Josephson junctions(10), particle coordination(11), not display any coherent behavior, whereas a strongly network with sufficiently homogeneous natural freque spin glass models(12, 13), neuroscience(14), deep brain stimu lation(15), chemical oscillations(16), biological locomotion(17) enable to synchronization. These two qualitatively distinct rhythmic applause(18), and countless other synchronization phe regimes are illustrated in Fig. 1 B and C. nomena(19-21). Finally, coupled oscillator models of the form Formally, the interaction n such phase oscillators is shown in[1]are canonical models of coupled limit cycle oscillators modeled by a connected graph G(, E, A)with nodes=1 (22)and serve as prototypical examples in complex networks n, edges&cv x v, and positive weights ag >0 for each un- studies(23-25) directed edge i, k E &. For pairs of noninteracting oscillators The coupled oscillator dynamics [1] feature i and j, the coupling weight ay is 0. We assume that the node set effect of the coupling described by the graph G(V, E, A)and the is partitioned as V=VIU V2, and we consider the following desynchronizing effect of the dissimilar natural general coupled oscillator model M+D=o-∑qsin(-9),i∈n Author cont and F.B. designed research: F.D. performed research; FD analyzed data; and F D. M C, and F.B. wrote the paper. [11 The authors deare no conflict of interest. D=0 (B-B),i∈v mc时mm pnas. org/cgil/doi/10.1073/pnas. 1212134110 PNAs| February5.2013|vo.110|no.6|2005-2010
Synchronization in complex oscillator networks and smart grids Florian Dörflera,b,1, Michael Chertkovb , and Francesco Bulloa a Center for Control, Dynamical Systems, and Computation, University of California, Santa Barbara, CA 93106; and b Center for Nonlinear Studies and Theory Division, Los Alamos National Laboratory, Los Alamos, NM 87545 Edited by Steven H. Strogatz, Cornell University, Ithaca, NY, and accepted by the Editorial Board November 14, 2012 (received for review July 16, 2012) The emergence of synchronization in a network of coupled oscillators is a fascinating topic in various scientific disciplines. A widely adopted model of a coupled oscillator network is characterized by a population of heterogeneous phase oscillators, a graph describing the interaction among them, and diffusive and sinusoidal coupling. It is known that a strongly coupled and sufficiently homogeneous network synchronizes, but the exact threshold from incoherence to synchrony is unknown. Here, we present a unique, concise, and closed-form condition for synchronization of the fully nonlinear, nonequilibrium, and dynamic network. Our synchronization condition can be stated elegantly in terms of the network topology and parameters or equivalently in terms of an intuitive, linear, and static auxiliary system. Our results significantly improve upon the existing conditions advocated thus far, they are provably exact for various interesting network topologies and parameters; they are statistically correct for almost all networks; and they can be applied equally to synchronization phenomena arising in physics and biology as well as in engineered oscillator networks, such as electrical power networks. We illustrate the validity, the accuracy, and the practical applicability of our results in complex network scenarios and in smart grid applications. nonlinear dynamics | power grids The scientific interest in the synchronization of coupled oscillators can be traced back to Christiaan Huygens’ seminal work on “an odd kind sympathy” between coupled pendulum clocks (1), and it continues to fascinate the scientific community to date (2, 3). A mechanical analog of a coupled oscillator network is shown in Fig. 1A and consists of a group of particles constrained to rotate around a circle and assumed to move without colliding. Each particle is characterized by a phase angle θi and has a preferred natural rotation frequency ωi. Pairs of interacting particles i and j are coupled through an elastic spring with stiffness aij. Intuitively, a weakly coupled oscillator network with strongly heterogeneous natural frequencies ωi does not display any coherent behavior, whereas a strongly coupled network with sufficiently homogeneous natural frequencies is amenable to synchronization. These two qualitatively distinct regimes are illustrated in Fig. 1 B and C. Formally, the interaction among n such phase oscillators is modeled by a connected graph G(V, E, A) with nodes V = {1, ..., n}, edges E ⊂ V × V, and positive weights aij > 0 for each undirected edge {i, k} ∈ E. For pairs of noninteracting oscillators i and j, the coupling weight aij is 0. We assume that the node set is partitioned as V = V1 ∪ V2, and we consider the following general coupled oscillator model: Miθ €i + Diθ _ i = ωi − Xn j=1 aij sin θi −θj ; i ∈V1 Diθ _ i = ωi − Xn j=1 aij sin θi −θj ; i ∈V2: [1] The coupled oscillator model [1] consists of the second-order oscillators V1 with Newtonian dynamics, inertia coefficients Mi, and viscous damping Di. The remaining oscillators V2 feature first-order dynamics with time constants Di. A perfect electrical analog of the coupled oscillator model [1] is given by the classic structure-preserving power network model (4), our enabling application of interest. Here, the first- and second-order dynamics correspond to loads and generators, respectively, and the right-hand sides depict the power injections ωi and the power flows aij sin(θi − θj ) along transmission lines. The rich dynamic behavior of the coupled oscillator model [1] arises from a competition between each oscillator’s tendency to align with its natural frequency ωi and the synchronizationenforcing coupling aij sin(θi − θj ) with its neighbors. If all natural frequencies ωi are identical, the coupled oscillator dynamics [1] collapse to a trivial phase-synchronized equilibrium, where all angles θi are aligned. The dissimilar natural frequencies ωi, on the other hand, drive the oscillator network away from this allaligned equilibrium. Moreover, even if the coupled oscillator model [1] synchronizes, the motion of its center of mass still carries the flux of angular rotation, respectively, the flux of electrical power from generators to loads in a power grid. Despite all these complications, the main result of this article is that, for a broad range of network topologies and parameters, an elegant and easily verified criterion characterizes synchronization of the nonlinear and nonequilibrium dynamic oscillator network [1]. Review of Synchronization in Oscillator Networks The coupled oscillator model [1] unifies various models in the literature, including dynamic models of electrical power networks. Modeling of electrical power networks is discussed in SI Text in detail. For V2 = , the coupled oscillator model [1] appears in synchronization phenomena in animal flocking behavior (5), populations of flashing fireflies (6), and crowd synchrony on London’s Millennium bridge (7), as well as in Huygen’s pendulum clocks (8). For V1 = , the coupled oscillator model (1) reduces to the celebrated Kuramoto model (9), which appears in coupled Josephson junctions (10), particle coordination (11), spin glass models (12, 13), neuroscience (14), deep brain stimulation (15), chemical oscillations (16), biological locomotion (17), rhythmic applause (18), and countless other synchronization phenomena (19–21). Finally, coupled oscillator models of the form shown in [1] are canonical models of coupled limit cycle oscillators (22) and serve as prototypical examples in complex networks studies (23–25). The coupled oscillator dynamics [1] feature the synchronizing effect of the coupling described by the graph G(V, E, A) and the desynchronizing effect of the dissimilar natural frequencies ωi. Author contributions: F.D., M.C., and F.B. designed research; F.D. performed research; F.D. analyzed data; and F.D., M.C., and F.B. wrote the paper. The authors declare no conflict of interest. This article is a PNAS Direct Submission. S.H.S. is a guest editor invited by the Editorial Board. 1 To whom correspondence should be addressed. E-mail: dorfler@engineering.ucsb.edu. This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10. 1073/pnas.1212134110/-/DCSupplemental. www.pnas.org/cgi/doi/10.1073/pnas.1212134110 PNAS | February 5, 2013 | vol. 110 | no. 6 | 2005–2010 APPLIED MATHEMATICS
distance smaller than some angle y E [0, r/l, that is, 10,-0lsy for every edge{i,j}∈E Based on a previously unexplored analysis approach to the synchronization problem, we propose the following synchroni zation condition for the coupled oscillator model [1 C≈ Synchronization Condition. The coupled oscillator model [1]has a unique and stable solution 0* with synchronized frequencies and cohesive phases 0--01sy<x/2 for every connected pair of con nected oscillators fi,j∈εif ‖ olle s sin() [2 Fig. 1. Mechanical analog of a lator network (A)and its cly coupl Here, L' is the pseudoinverse of the network Laplacian matrix L thts a w, all parameters in simulations and lille, o=maxyijlecki-xil is the worst-case dissimilarity forx We establish the broad applicability of the proposed conditio [2 to various classes of networks via analytical and statistical The complex network com ity asks questions of the form ethods in the next section. Before that, we provide some equiv- "What are the conditions on the coupling and the dissimilarity alent formulations for the proposed condition [2] to develop such that a synchronizing behavior emerges?" Similar questions deeper intuition and obtain insightful conclusion in large-scale electrical power systems. Because synchronization Complex Network Interpretation. Surprisingly, topological or spec is pervasive in the operation of an interconnected power grid, tral connectivity measures, such as nodal degree or algebraic parameters and topology, the current load profile, and power advocated connectivity measures(23-25, 30,33-35)turn out to 27), when is it optimal(28), when is it stable(29, 30), and how This statement can be seen by introducing the matrix U of or- robust is it (31 local loss of synchrony can trigger cas- thonormal eigenvectors of the network placlan matr the face of the complexity of future smart grids and the inte- spectral viewpoint, condition [2] can be equivalently written as cading failures and possibly result in widespread blackouts. In corresponding eigenvalues 0= A1< B2 <. ..s im From gration challenges posed by renewable energy sources, a deeper ‖ U diag(0,1/2,…,1/xn)·(Ura)ls≤sn(y).B3 Despite the vast scientific interest, the search for sharp, co cise, and closed-form synchronization conditions for coupled In words, the natural frequencies o are projected on the network oscillator models of the form shown in [1] has been in vain so far. modes U, weighted by the inverse Laplacian eigenvalues, and Loosely speaking, synchronization occurs when the coupling ll e o evaluates the worst-case dissimilarity of this weighted dominates the dissimilarity. Various conditions have been pro- projection. A sufficient condition for the inequality [3] to be true posed to quantify this tradeoff (21, 23-25, 30, 35-36). The cou- is the algebraic connectivity condition A2> olle ao sin(y). Like pling is typically quantified by the nodal degree or the algebraic wise, a necessary condition for inequality 3 is 2-deg(G) connectivity of the graph G, and the dissimilarity is quantified an2lolls sin(y), where deg(G) is the maximum nodal degree by the magnitude or the spread of the natural frequencies in the graph G(, E, A).Clearly, compared with 3], this sufficient Sometimes, these conditions can be evaluated only numerically condition and this necessary condition feature only one of n-1 because they depend on the network state (33, 34)or arise from nonzero Laplacian eigenvalues and are overly conservative. a nontrivial linearization process, such as the Master stability function formalism(23, 24). To date, exact synchronization con- Kuramoto Oscillator Perspective. Notice that in the limit y=x/2, ditions are known only for simple coupling topologies( 17, 21, 37, condition [2] suggests that there exists a stable synchronized 68). For arbitrary topologies, only sufficient conditions are known solution if (25, 30, 33-35), as well as numerical investigations for random networks(39-41) Simulation studies indicate that the known ‖Lols<1 sufficient conditions are very conservative estimates on the thresh- old from incoherence to synchrony. Literally, every review article For classic Kuramoto oscillators coupled in a compl parameters(19-21, 23-25). In this article, we present a conti c duces to the condition K> maxj e(L. -ol, known for the ssic Kuramoto model (21) nd sharp synchronization condition that features elegant graph Power Network Perspective equilibrium equations of the scillator model [1]. given For the coupled oscillator model [1] and its applications, the equations, and they are often approximated. he AC power flow Synchronization Condition to as following notions of synchronization are appropriate. First, a @=Ei=ai(0i,)(31-34), known as the DC power flow solution has synchronized frequencies if all frequencies ]i are equations In vector notation, the DC power flow equations read identical to a common constant value @sync. If a synchronized as o= LO and their solution satis solution exists, it is known that the synchronization frequency is LToll2ao. According to condition[2 the worst phase distance Ok/Ck- Dk and that by working in a rotating ref- L oll o obtained by the dC power flow equations needs to be erence frame, one may assume @ync =0. Second, a solution has less than or equal to sin(r) such that the solution to the AC cohesive phases if every pair of connected oscillators has a phase power flow equations satisfies maxjjlee[0i -,lsy. Hence,our 2006iwww.pnas.org/cgi/doi/10.1073/pnas.1212134110
The complex network community asks questions of the form “What are the conditions on the coupling and the dissimilarity such that a synchronizing behavior emerges?” Similar questions also appear in all the aforementioned applications, for instance, in large-scale electrical power systems. Because synchronization is pervasive in the operation of an interconnected power grid, a central question is “Under which conditions on the network parameters and topology, the current load profile, and power generation does there exist a synchronous operating point (26, 27), when is it optimal (28), when is it stable (29, 30), and how robust is it (31–34)?” A local loss of synchrony can trigger cascading failures and possibly result in widespread blackouts. In the face of the complexity of future smart grids and the integration challenges posed by renewable energy sources, a deeper understanding of synchronization is increasingly important. Despite the vast scientific interest, the search for sharp, concise, and closed-form synchronization conditions for coupled oscillator models of the form shown in [1] has been in vain so far. Loosely speaking, synchronization occurs when the coupling dominates the dissimilarity. Various conditions have been proposed to quantify this tradeoff (21, 23–25, 30, 35–36). The coupling is typically quantified by the nodal degree or the algebraic connectivity of the graph G, and the dissimilarity is quantified by the magnitude or the spread of the natural frequencies ωi. Sometimes, these conditions can be evaluated only numerically because they depend on the network state (33, 34) or arise from a nontrivial linearization process, such as the Master stability function formalism (23, 24). To date, exact synchronization conditions are known only for simple coupling topologies (17, 21, 37, 38). For arbitrary topologies, only sufficient conditions are known (25, 30, 33–35), as well as numerical investigations for random networks (39–41). Simulation studies indicate that the known sufficient conditions are very conservative estimates on the threshold from incoherence to synchrony. Literally, every review article on synchronization concludes emphasizing the quest for exact synchronization conditions for arbitrary network topologies and parameters (19–21, 23–25). In this article, we present a concise and sharp synchronization condition that features elegant graphtheoretical and physical interpretations. Synchronization Condition For the coupled oscillator model [1] and its applications, the following notions of synchronization are appropriate. First, a solution has synchronized frequencies if all frequencies θ _ i are identical to a common constant value ωsync. If a synchronized solution exists, it is known that the synchronization frequency is ωsync =Pn k=1ωk= Pn k=1Dk and that by working in a rotating reference frame, one may assume ωsync = 0. Second, a solution has cohesive phases if every pair of connected oscillators has a phase distance smaller than some angle γ ∈ [0, π/2[, that is, jθi − θj j ≤ γ for every edge {i, j} ∈ E. Based on a previously unexplored analysis approach to the synchronization problem, we propose the following synchronization condition for the coupled oscillator model [1]. Synchronization Condition. The coupled oscillator model [1] has a unique and stable solution θ* with synchronized frequencies and cohesive phases jθ* i − θ* j j≤γ < π=2 for every connected pair of connected oscillators {i, j} ∈ E if L† ω E;∞ ≤ sinðγÞ: [2] Here, L† is the pseudoinverse of the network Laplacian matrix L and kxkE;∞ = maxfi;jg∈ E jxi −xjj is the worst-case dissimilarity for x = (x1, ..., xn) over the edges E. We establish the broad applicability of the proposed condition [2] to various classes of networks via analytical and statistical methods in the next section. Before that, we provide some equivalent formulations for the proposed condition [2] to develop deeper intuition and obtain insightful conclusions. Complex Network Interpretation. Surprisingly, topological or spectral connectivity measures, such as nodal degree or algebraic connectivity, are not key to synchronization. In fact, these often advocated connectivity measures (23–25, 30, 33–35) turn out to be conservative estimates of the synchronization condition [2]. This statement can be seen by introducing the matrix U of orthonormal eigenvectors of the network Laplacian matrix L with corresponding eigenvalues 0 = λ1 < λ2 ≤ ... ≤ λn. From this spectral viewpoint, condition [2] can be equivalently written as U diag ð0; 1=λ2; ... ; 1=λnÞ· UTω E;∞ ≤ sinðγÞ: [3] In words, the natural frequencies ω are projected on the network modes U, weighted by the inverse Laplacian eigenvalues, and k · kE;∞ evaluates the worst-case dissimilarity of this weighted projection. A sufficient condition for the inequality [3] to be true is the algebraic connectivity condition λ2 ≥ kωkE;∞ · sinðγÞ. Likewise, a necessary condition for inequality [3] is 2 · degðGÞ ≥ λn ≥kωkE;∞ · sinðγÞ, where deg(G) is the maximum nodal degree in the graph G(V, E, A). Clearly, compared with [3], this sufficient condition and this necessary condition feature only one of n − 1 nonzero Laplacian eigenvalues and are overly conservative. Kuramoto Oscillator Perspective. Notice that in the limit γ → π/2, condition [2] suggests that there exists a stable synchronized solution if L† ω E;∞ < 1: [4] For classic Kuramoto oscillators coupled in a complete graph with uniform weights aij = K/n, the synchronization condition [4] reduces to the condition K > maxi;j∈f1;...;ngjωi −ωjj, known for the classic Kuramoto model (21). Power Network Perspective. In power systems engineering, the equilibrium equations of the coupled oscillator model [1], given by ωi =Pn j=1aij sinðθi −θjÞ, are referred to as the AC power flow equations, and they are often approximated by their linearization ωi =Pn j=1aijðθi − θjÞ (31–34), known as the DC power flow equations. In vector notation, the DC power flow equations read as ω = Lθ and their solution satisfies maxfi;jg∈E jθi −θjj= kL†ωkE;∞. According to condition [2], the worst phase distance kL†ωkE;∞ obtained by the DC power flow equations needs to be less than or equal to sin(γ) such that the solution to the AC power flow equations satisfies maxfi;jg∈E jθi −θjj ≤γ. Hence, our A B C Fig. 1. Mechanical analog of a coupled oscillator network (A) and its dynamics in a strongly coupled network (B) and a weakly coupled network (C). With the exception of the coupling weights aij, all parameters in simulations B and C are identical. 2006 | www.pnas.org/cgi/doi/10.1073/pnas.1212134110 Dörfler et al.
extends the common dC power flow approxima five, (vi)arbitrary cycles with symmetrical parameters, and(vii) ngles y< I to large anglesy E [0 combinations of networks each satisfying one of the conditions(Or- (vi), which are connected to another and share no co Auxiliary Linear Perspective. As detailed in the previous section, A detailed and rigorous mathematical derivation and state the key term L'a in condition [2] equals the phase differences ment of the above analytical result can be found in SI Text obtained by the linear Laplacian equation o= LA. This linear interpretation is not only insightful but practical, because con- weights ay are known only with a certain degree of accuracy, or dition [2] can be quickly evaluated by numerically solving the they may be variable within certain ranges. For instance, in power linear system @=Le. This linear system is possibly of high di- networks, these variations arise from uncertain demand or un- but it inherits the parity of the graph G(, E, A). modeled voltage dynamics. To address these uncertainties, Thus, condition [2] can be verified efficiently even for large-scale condition [2] can be extended to a robust Snco to verity onization that our derivation of condition [2] is not based on any lineari- zation arguments ≤可· In this case, it is necessary and sufficient condition(2 at the vertices of the parameter space o E(oi,oil Energy Landscape Condition[2] can also be un- and ai etaj, ail to guarantee condition [2] for all possible derstood in terms aling energy landscape interpreta- parameter variations. The detailed results are reported in Si Tex tion. The coupled nodel [1] is a system of particles that alm t of particular network topologies and parameters, we establish its correctness and predictive power for a broad range of networks. Extensive simulation studies lead to the conclusion that the pr E(O=∑1(1-cos(-9)-∑· posed synchronization condition [2] is statistically correct. To verify this hypothesis, we conducted Monte Carlo simulation studies over a wide range of natural frequencies @i, network sizes where the first term is a pairwise nonlinear attraction among coupling weights a, and different random graph models of varying 9 driving the particles away from the"all-aligned"state. Because network models with topologies constructed from Erdos-Re the energy function E(0)is difficult to study, it is natural to graphs, random geometric graphs, and Watts-Strogatzsmall world ponds to a Hookean potential. Condition e. first tem o(0)2 networks, and the natural frequencies and coupling weights are look for a minimum of its second-order approximation ∑aesn(-9)212-∑=1日, amped from uniform distributions. In total, we constructed corre- 1.2.10 samples of such nominal networks, each with a connected acting particles no further than y apart if Eo() features a min- results can be found in si Text and allow us to establish the fol- imum with interacting particles no further from each other than lowing probabilistic result with a confidence level of at least 99% nd accuracy of at least 99% Analytical and Statistical Results Statistical Result for Nominal Networks With a 99.97%o probability Our analysis approach to the synchronization problem is based for a nominal network, condition 21 the synchronization condition [2]. In particular, we analytically chronization condition [ 2] holds true for alm establish the synchronization condition [2] for the following and parameters of the considered nominal network i deed we also show the existence of possibly thin sets of ne Analytical Result. The synchronization 2 is necessary and sufficiently tight. We refer to SI Text for an explicit family of sufficient for() the sparsest(acyclic)and (i) the densest(complete carefully engineered and"degenerate"counterexamples. Overall, and uniformly weighted) network topol GO,E, A),(ii) the our analytical and statistical results validate the correctness best(phase synchronizing)and (iv) the worst (cut-set inducing) of the proposed condition(2] natural frequencies, (v)cyclic topologies of length strictly less than After having established the statistical correctness of condition 21, we now investigate its predictive power for arbitrary ne works. Because we analytically establish that condition[2] is exact E(0),E0() for sufficiently small pairwise phase cohesiveness 10i-8l<1,we now investigate the other extreme, 10i-0,1=x/2. To test the corresponding condition [4] in a low-dimensional parameter space, we consider a complex network of Kuramoto oscillators 1 where all coupling weights ai are either 0 or 1 and the coupl B1-82Irad gain K>0 serves as a control parameter. If L is the corresponding Fig.2. Energy function E(0) and its quadratic approximation Eo(e)for a two. unweighted Laplacian matrix, condition [4] reads as K> article system are shown as solid and dashed curves, respectively, for the KcriticalAL'ol Of course, the condition k>Critical is only able(blue), marginal (green), and unstable (red) cases. The circles and sufficient, and the critical coupling may be smaller than kcritie diamonds represent stable critical points of E(e) and Eo(e). To test the accuracy of the condition k> Critical, we numerically Dorfler et al PNAS I February 5, 2013 I vol. 110 I no6 2007
condition extends the common DC power flow approximation from infinitesimally small angles γ 1 to large angles γ ∈ [0, π/2[. Auxiliary Linear Perspective. As detailed in the previous section, the key term L† ω in condition [2] equals the phase differences obtained by the linear Laplacian equation ω = Lθ. This linear interpretation is not only insightful but practical, because condition [2] can be quickly evaluated by numerically solving the linear system ω = Lθ. This linear system is possibly of high dimension, but it inherits the sparsity of the graph G(V, E, A). Thus, condition [2] can be verified efficiently even for large-scale sparse networks. Despite this linear interpretation, we emphasize that our derivation of condition [2] is not based on any linearization arguments. Energy Landscape Perspective. Condition [2] can also be understood in terms of an appealing energy landscape interpretation. The coupled oscillator model [1] is a system of particles that aim to minimize the energy function EðθÞ= X fi;jg∈E aij 1 −cos θi −θj − Xn i= 1 ωi · θi; where the first term is a pairwise nonlinear attraction among the particles and the second term represents the external force driving the particles away from the “all-aligned” state. Because the energy function E(θ) is difficult to study, it is natural to look for a minimum of its second-order approximation E0ðθÞ= P fi;jg∈E aijðθi − θjÞ 2 =2 −Pn i=1ωi · θi, where the first term corresponds to a Hookean potential. Condition [2] is then restated as follows: E(θ) features a phase-cohesive minimum with interacting particles no further than γ apart if E0(θ) features a minimum with interacting particles no further from each other than sin(γ), as illustrated in Fig. 2. Analytical and Statistical Results Our analysis approach to the synchronization problem is based on algebraic graph theory. We propose an equivalent reformulation of the synchronization problem that reveals the crucial role of cycles and cut-sets in the graph, and ultimately leads to the synchronization condition [2]. In particular, we analytically establish the synchronization condition [2] for the following interesting cases. Analytical Result. The synchronization condition [2] is necessary and sufficient for (i) the sparsest (acyclic) and (ii) the densest (complete and uniformly weighted) network topologies G(V, E, A), (iii) the best (phase synchronizing) and (iv) the worst (cut-set inducing) natural frequencies, (v) cyclic topologies of length strictly less than five, (vi) arbitrary cycles with symmetrical parameters, and (vii) combinations of networks each satisfying one of the conditions (i)– (vi), which are connected to another and share no common cycles. A detailed and rigorous mathematical derivation and statement of the above analytical result can be found in SI Text. In many applications, the natural frequencies ωi and coupling weights aij are known only with a certain degree of accuracy, or they may be variable within certain ranges. For instance, in power networks, these variations arise from uncertain demand or unmodeled voltage dynamics. To address these uncertainties, condition [2] can be extended to a robust synchronization condition for variable parameters ωi ≤ ωi ≤ωi and 0< aij ≤ aij ≤ aij. In this case, it is necessary and sufficient to verify condition [2] at the vertices of the parameter space ωi ∈ f ωi ;ωig and aij ∈f aij ; aijg to guarantee condition [2] for all possible parameter variations. The detailed results are reported in SI Text. After having analytically established condition [2] for a variety of particular network topologies and parameters, we establish its correctness and predictive power for a broad range of networks. Extensive simulation studies lead to the conclusion that the proposed synchronization condition [2] is statistically correct. To verify this hypothesis, we conducted Monte Carlo simulation studies over a wide range of natural frequencies ωi , network sizes n, coupling weights aij, and different random graph models of varying degrees of sparsity and randomness. We select a set of nominal network models with topologies constructed from Erdös–Rényi graphs, random geometric graphs, and Watts–Strogatz small world networks, and the natural frequencies and coupling weights are sampled from uniform distributions. In total, we constructed 1.2·106 samples of such nominal networks, each with a connected graph G(V, E, A) and natural frequencies ω satisfying kL†ωkE;∞ ≤ sinðγÞ for some γ < π/2. The detailed construction and precise results can be found in SI Text and allow us to establish the following probabilistic result with a confidence level of at least 99% and accuracy of at least 99%. Statistical Result for Nominal Networks. With a 99.97% probability for a nominal network, condition [2] guarantees the existence of a unique and stable solution θ* with synchronized frequencies and cohesive phases jθ* i −θ* j j≤γ for every connected pair {i, j} ∈ E. From this statistical result, we deduce that the proposed synchronization condition [2] holds true for almost all topologies and parameters of the considered nominal network models. Indeed, we also show the existence of possibly thin sets of network topologies and parameters for which our condition [2] is not sufficiently tight. We refer to SI Text for an explicit family of carefully engineered and “degenerate” counterexamples. Overall, our analytical and statistical results validate the correctness of the proposed condition [2]. After having established the statistical correctness of condition [2], we now investigate its predictive power for arbitrary networks. Because we analytically establish that condition [2] is exact for sufficiently small pairwise phase cohesiveness jθi − θj j 1, we now investigate the other extreme, maxfi;jg∈E jθi −θjj= π=2. To test the corresponding condition [4] in a low-dimensional parameter space, we consider a complex network of Kuramoto oscillators θ_ i = ωi −K · Xn j= 1 aij sin θi − θj ; i ∈f1; ... ; ng; [5] where all coupling weights aij are either 0 or 1 and the coupling gain K > 0 serves as a control parameter. If L is the corresponding unweighted Laplacian matrix, condition [4] reads as K > Kcritical ≜ jjL†ωjjE;∞. Of course, the condition K > Kcritical is only sufficient, and the critical coupling may be smaller than Kcritical. To test the accuracy of the condition K > Kcritical, we numerically Fig. 2. Energy function E(θ) and its quadratic approximation E0(θ) for a twoparticle system are shown as solid and dashed curves, respectively, for the stable (blue), marginal (green), and unstable (red) cases. The circles and diamonds represent stable critical points of E(θ) and E0(θ). Dörfler et al. PNAS | February 5, 2013 | vol. 110 | no. 6 | 2007 APPLIED MATHEMATICS
Erdos-Renyi Graph Random Geometric Graph: Small World Network C 150751 -----1◆n=40 E F P挂+=10 aa。a! of the exact critical coupling k in a complex Kuramoto oscillator network. The subfigures show K normalized by Lolles for an pability p of connecting two nodes, for a rand with connect probability p. Each data point is the mean over 100 samples of the respective om graph model for values of sampl istribution supported on (-1, 1] and for the network sizes n E [10, found the smallest value of K leading to synchrony with phase Applications in Power Networks We envision that condition [2] can be applied to assess syn Fig 3 reports our findings for various network sizes, connected chronization and robustness quickly in power networks under uencies. We refer to S/ Text for the detailed simulation setup. First, works ac erating conditions.Because re random graph models, and sample distributions of the natural fre- volatile op carefully engineered systems with particular network note from Fig 3A, B, D, and E that condition (4] is extremely ac- topologies and parameters, we do not extrapolate the statistica curate for a sparse graph, that is, for small p andn, as expected from results from the previous section to power grids. Rather, we con- our analytical results. Second, for a dense graph with s 1, Fig 3A, sider 10 widely established Institute of Electrical and Electronics B.D. and e confirms the results known for classic Kuramoto os Engineers (IEEE) power network test cases provided by Zim cillators(21): For a bipolar distribution, condition [4] is exact, and merman et al. (42)and Grigg et al. (43) for a uniform distribution, a small critical coup is obtained Under nominal operating conditions, the power generation is to meet the forecast demand while obeying the AC Watts-Strogatz small world network, that is, it has almost constant ow laws and respecting the thermal limits of each trans- curacy for various values of n and p. Fourth and finally, observe ne. Thermal limit constraints are precisely equivalent that condition [4] is always within a constant factor of the exact to phase cohesiveness requirements. To test the synchronization critical coupling, whereas other proposed conditions(23-25, 30, 33- condition [2] in a volatile smart grid scenario, we make the fol- 5)on the nodal degree or on the algebraic connectivity scale poorly lowing changes to the nominal network: (i) We assume fluctu with respect to network size n ating demand and randomize 50% of all loads to deviate from the Table 1. Evaluation of condition [2] for 10 IEEE test cases under volatile operating conditions Randomized test case(1,000 instances) Correctness* Accuracy, rad Chow 9 bus system 4.121810-5 0.12889 IEEE 14 bus system Always true 2.799510-4 0.16622 IEEE RTS 24 Always true 1.708910-3 0.22309 IEEE 30 bus system Always true 0.16430 New England 39 0.1682 IEEE 57 bus system Always true 0.20295 IEEE RTS 96 Always true 0.24593 IEEE 118 bus system Always true 5.995910-4 IEEE 300 bus system Always true 5,261810m4 043204 Polish 2383 bus system(winter 1999) Always true 4.218310-3 0.25144 Imaxuilezle-0l-arcsin(It' alle, m)I 2008Iwww.pnas.org/cgi/doi/10.1073/pnas.1212134110
found the smallest value of K leading to synchrony with phase cohesiveness π/2. Fig. 3 reports our findings for various network sizes, connected random graph models, and sample distributions of the natural frequencies. We refer to SI Text for the detailed simulation setup. First, note from Fig. 3 A, B, D, and E that condition [4] is extremely accurate for a sparse graph, that is, for small p and n, as expected from our analytical results. Second, for a dense graph with p ≈ 1, Fig. 3 A, B, D, and E confirms the results known for classic Kuramoto oscillators (21): For a bipolar distribution, condition [4] is exact, and for a uniform distribution, a small critical coupling is obtained. Third, Fig. 3 C and F shows that condition [4] is scale-free for a Watts–Strogatz small world network, that is, it has almost constant accuracy for various values of n and p. Fourth and finally, observe that condition [4] is always within a constant factor of the exact critical coupling, whereas other proposed conditions (23–25, 30, 33– 35) on the nodal degree or on the algebraic connectivity scale poorly with respect to network size n. Applications in Power Networks We envision that condition [2] can be applied to assess synchronization and robustness quickly in power networks under volatile operating conditions. Because real-world power networks are carefully engineered systems with particular network topologies and parameters, we do not extrapolate the statistical results from the previous section to power grids. Rather, we consider 10 widely established Institute of Electrical and Electronics Engineers (IEEE) power network test cases provided by Zimmerman et al. (42) and Grigg et al. (43). Under nominal operating conditions, the power generation is optimized to meet the forecast demand while obeying the AC power flow laws and respecting the thermal limits of each transmission line. Thermal limit constraints are precisely equivalent to phase cohesiveness requirements. To test the synchronization condition [2] in a volatile smart grid scenario, we make the following changes to the nominal network: (i) We assume fluctuating demand and randomize 50% of all loads to deviate from the ABC DEF Fig. 3. Numerical evaluation of the exact critical coupling K in a complex Kuramoto oscillator network. The subfigures show K normalized by kL†ωkE;∞ for an Erdös–Rényi graph with probability p of connecting two nodes, for a random geometric graph with connectivity radius p, and for a Watts–Strogatz small world network with rewiring probability p. Each data point is the mean over 100 samples of the respective random graph model for values of ωi sampled from a bipolar or a uniform distribution supported on [−1, 1] and for the network sizes n ∈ {10, 20, 40, 80, 160}. Table 1. Evaluation of condition [2] for 10 IEEE test cases under volatile operating conditions Randomized test case (1,000 instances) Correctness* Accuracy† , rad Cohesive‡ , rad Chow 9 bus system Always true 4.1218·10−5 0.12889 IEEE 14 bus system Always true 2.7995·10−4 0.16622 IEEE RTS 24 Always true 1.7089·10−3 0.22309 IEEE 30 bus system Always true 2.6140·10−4 0.16430 New England 39 Always true 6.6355·10−5 0.16821 IEEE 57 bus system Always true 2.0630·10−2 0.20295 IEEE RTS 96 Always true 2.6076·10−3 0.24593 IEEE 118 bus system Always true 5.9959·10−4 0.23524 IEEE 300 bus system Always true 5.2618·10−4 0.43204 Polish 2383 bus system (winter 1999) Always true 4.2183·10−3 0.25144 Accuracy and phase cohesiveness results are averaged over 1,000 instances of randomized load and generation. *Correctness: kL†ωkE;∞ ≤sinðγÞ0maxfi;jg∈E jθ* i −θ* j j≤γ: † Accuracy: jmaxfi;jg∈E jθ* i −θ* j j−arcsinðkL†ωkE;∞Þj: ‡ Phase cohesiveness: maxfi;jg∈E jθ* i −θ* j j: 2008 | www.pnas.org/cgi/doi/10.1073/pnas.1212134110 Dörfler et al
the following two contingencies have taken place, and we char acterize the rem margin. First, we assume ge nerator 323 is disconnected, possibly due to maintenance or failure events. Second, we consider the following imbalanced power dispatch situation: The power demand at each load in the southeastern area deviates from the nominally forecasted demand by a uniform and positive amount, and the resulting power deficiency is compen- sated for by uniformly increasing the generation in the north western area. This imbalance can arise, for example, due to a shortfall in predicted load and renewable energy generation Correspondingly, power is exported from the northwestern to the southeastern area via the transmission lines(121, 325) and 1223, 318. At a nominal operating condition, the RTS 96 power net work is sufficiently robust to tolerate each single one of these two Fg. 4. Lustration of contingencies in the rTS 96 power network. Here. predicts that the thermal limit of the transmission line (121, 325)is power are exported from the northwestern area to the southeastern area, reached at an additional loading of 22.20%. Indeed, the dynamic and generator 323 is tripped. prediction. It can be observed that synchronization is lost for an additional loading of 22.33% and that the areas separate via the forecasted loads;(ii) we assume that the grid is penetrated by transmission line 1121, 325). This separation triggers a cascade of renewables with severely fluctuating power outputs, for example, events, such as an outage of the transmission line (223, 318), and wind or solar farms, and we randomize 33% of all generating the power network is en route to a blackout. We remark that, if units to deviate from the nominally scheduled generation; and generator 323 is not disconnected and there are no thermal limit (iii) following the paradigm of smart operation of smart grids(44), constraints, by increasing the loading, we observe the classic loss of the fluctuations can be mitigated by fast-ramping generation, such synchrony through a saddle-node bifurcation. This bifurcation can as fast-response energy storage, including batteries and flywheels, also be predicted accurately by our results(a detailed description and controllable loads, such as large-scale server farms or fleets of is provided in SI Texr) plug-in hybrid electrical vehicles. Here, we assume that the grid is In summary, the results in this section confirm the validity, the uipped with 10% fast-ramping generation and 10% controlla- applicability, and the accuracy of the synchronization condition ble loads, and that the power ce(caused by fluctuating [2] in complex power network scenarios. demand and generation) is uniformly dispatched among these adjustable power sources. For each of the 10 iEEE test cases, we Discussion and Conclusions construct 1,000 random realizations of scenarios ((iin de- In this article, we studied the synchronization phenomen scribed above, we numerically check for the existence of a syn- class of coupled oscillator models proposed in the scient chronous solution, and we compare the numerical solution with ture. We proposed a surprisingly simple condition that the results predicted by our condition [2]. Our findings are predicts synchronization as a function of the parameters and the reported in Table 1, and a detailed description of the simulation topology of the underlying network. Our result, with its physical and setup can be found in SI Tex graph theoretical interpretations, significantly improves upon the It can be observed that condition 2 predicts the correct phase existing tests in the literature on cohesiveness 1ei-el along all transmission lines i,J) EE with our synchronization condition is established analytically for vari- extremely high accuracy, even for large-scale networks featuring ous interesting network topolog 83 nodes. These conclusions can also be extended to power lations for a broad range of generic networks. We validated our network models with variable parameters, which account for theoretical results for complex Kuramoto oscillator networks as uncertainty in demand or unmodeled voltage dynamics. Further well as for smart grid applications. provided in SI Tex Our results pose as many questions as they answer. Among the As a final test, we validate the synchronization condition [2] in important theoretical problems to be addressed is a characteriza a stressed power grid case study. We consider the IEEE Reliability tion of the set of all network topologies and parameters for which Test System 96(RTS 96)(43)illustrated in Fig 4. We assume that our proposed synchronization condition Lolle <1 is not A B nuous load increase from 22. 19% to 22.24%.(A)Angles e(t) that lose t<t 1977 rad of the e(t)[rad smission line (121, 325) is (t”) reached. (B)Angles e(t)at t= t.(C) E (121, 325) and 1223, 318 are plotted t t dashed curves. (D and E) tor phase space (e(t)e(t))before t>t d after t, where the loss of e(t)rad common synchronization frequency can be observed Dorfler et al PNAS I February 5, 2013 I vol. 110 I no6I 2009
forecasted loads; (ii) we assume that the grid is penetrated by renewables with severely fluctuating power outputs, for example, wind or solar farms, and we randomize 33% of all generating units to deviate from the nominally scheduled generation; and (iii) following the paradigm of smart operation of smart grids (44), the fluctuations can be mitigated by fast-ramping generation, such as fast-response energy storage, including batteries and flywheels, and controllable loads, such as large-scale server farms or fleets of plug-in hybrid electrical vehicles. Here, we assume that the grid is equipped with 10% fast-ramping generation and 10% controllable loads, and that the power imbalance (caused by fluctuating demand and generation) is uniformly dispatched among these adjustable power sources. For each of the 10 IEEE test cases, we construct 1,000 random realizations of scenarios (i)–(iii) described above, we numerically check for the existence of a synchronous solution, and we compare the numerical solution with the results predicted by our condition [2]. Our findings are reported in Table 1, and a detailed description of the simulation setup can be found in SI Text. It can be observed that condition [2] predicts the correct phase cohesiveness jθi − θjj along all transmission lines {i, j} ∈ E with extremely high accuracy, even for large-scale networks featuring 2,383 nodes. These conclusions can also be extended to power network models with variable parameters, which account for uncertainty in demand or unmodeled voltage dynamics. Further details are provided in SI Text. As a final test, we validate the synchronization condition [2] in a stressed power grid case study. We consider the IEEE Reliability Test System 96 (RTS 96) (43) illustrated in Fig. 4. We assume that the following two contingencies have taken place, and we characterize the remaining safety margin. First, we assume generator 323 is disconnected, possibly due to maintenance or failure events. Second, we consider the following imbalanced power dispatch situation: The power demand at each load in the southeastern area deviates from the nominally forecasted demand by a uniform and positive amount, and the resulting power deficiency is compensated for by uniformly increasing the generation in the northwestern area. This imbalance can arise, for example, due to a shortfall in predicted load and renewable energy generation. Correspondingly, power is exported from the northwestern to the southeastern area via the transmission lines {121, 325} and {223, 318}. At a nominal operating condition, the RTS 96 power network is sufficiently robust to tolerate each single one of these two contingencies, but the safety margin is now minimal. When both contingencies are combined, our synchronization condition [2] predicts that the thermal limit of the transmission line {121, 325} is reached at an additional loading of 22.20%. Indeed, the dynamic simulation scenario shown in Fig. 5 validates the accuracy of this prediction. It can be observed that synchronization is lost for an additional loading of 22.33% and that the areas separate via the transmission line {121, 325}. This separation triggers a cascade of events, such as an outage of the transmission line {223, 318}, and the power network is en route to a blackout. We remark that, if generator 323 is not disconnected and there are no thermal limit constraints, by increasing the loading, we observe the classic loss of synchrony through a saddle-node bifurcation. This bifurcation can also be predicted accurately by our results (a detailed description is provided in SI Text). In summary, the results in this section confirm the validity, the applicability, and the accuracy of the synchronization condition [2] in complex power network scenarios. Discussion and Conclusions In this article, we studied the synchronization phenomenon for broad class of coupled oscillator models proposed in the scientific literature. We proposed a surprisingly simple condition that accurately predicts synchronization as a function of the parameters and the topology of the underlying network. Our result, with its physical and graph theoretical interpretations, significantly improves upon the existing tests in the literature on synchronization. The correctness of our synchronization condition is established analytically for various interesting network topologies and via Monte Carlo simulations for a broad range of generic networks. We validated our theoretical results for complex Kuramoto oscillator networks as well as for smart grid applications. Our results pose as many questions as they answer. Among the important theoretical problems to be addressed is a characterization of the set of all network topologies and parameters for which our proposed synchronization condition kL†ωkE;∞ <1 is not Fig. 4. Illustration of contingencies in the RTS 96 power network. Here, square nodes are generators and round nodes are loads, large amounts of power are exported from the northwestern area to the southeastern area, and generator 323 is tripped. A B C D E Fig. 5. RTS 96 dynamics for a continuous load increase from 22.19% to 22.24%. (A) Angles θ(t) that lose synchrony at t* = 18.94 s, when the thermal limit γ* = 0.1977 rad of the transmission line {121, 325} is reached. (B) Angles θ(t) at t = t*. (C) Angular distances and the thermal limits γ* and γ**, where the lines {121, 325} and {223, 318} are plotted as dashed curves. (D and E) Generator phase space ðθðtÞ; θ _ðtÞÞ before and after t*, where the loss of a common synchronization frequency can be observed. Dörfler et al. PNAS | February 5, 2013 | vol. 110 | no. 6 | 2009 APPLIED MATHEMATICS