复旦大学：《网络科学导论 Introduction to Network Science》教学参考文献_Synchronization in complex oscillator networks and smart grids

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 oscil￾lators 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 describ￾ing 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 synchroni￾zation 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 net￾work 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 net￾work 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 un￾directed 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 dy￾namics 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 synchronization￾enforcing 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 all￾aligned 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 elec￾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 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 net￾works. 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 be￾havior (5), populations of flashing fireflies (6), and crowd syn￾chrony 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 stimu￾lation (15), chemical oscillations (16), biological locomotion (17), rhythmic applause (18), and countless other synchronization phe￾nomena (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 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 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 cas￾cading failures and possibly result in widespread blackouts. In the face of the complexity of future smart grids and the inte￾gration challenges posed by renewable energy sources, a deeper understanding of synchronization is increasingly important. Despite the vast scientific interest, the search for sharp, con￾cise, 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 pro￾posed to quantify this tradeoff (21, 23–25, 30, 35–36). The cou￾pling 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 con￾ditions 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 thresh￾old 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 graph￾theoretical 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 ref￾erence 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 synchroni￾zation 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≤γ 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 dy￾namics 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 y0 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 con￾dition [2] can be quickly evaluated by numerically solving the linear system ω = Lθ. This linear system is possibly of high di￾mension, 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 lineari￾zation arguments. Energy Landscape Perspective. Condition [2] can also be un￾derstood in terms of an appealing energy landscape interpreta￾tion. 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 corre￾sponds to a Hookean potential. Condition [2] is then restated as follows: E(θ) features a phase-cohesive minimum with inter￾acting particles no further than γ apart if E0(θ) features a min￾imum 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 reformu￾lation 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 state￾ment 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 un￾modeled voltage dynamics. To address these uncertainties, condition [2] can be extended to a robust synchronization condition for variable parameters ωi ≤ ωi ≤ωi and 0 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 two￾particle 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 fre￾quencies. 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 ac￾curate 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 os￾cillators (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 syn￾chronization and robustness quickly in power networks under volatile operating conditions. Because real-world power net￾works 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 con￾sider 10 widely established Institute of Electrical and Electronics Engineers (IEEE) power network test cases provided by Zim￾merman 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 trans￾mission 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 fol￾lowing changes to the nominal network: (i) We assume fluctu￾ating 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 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% controlla￾ble 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) de￾scribed above, we numerically check for the existence of a syn￾chronous 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 char￾acterize 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 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 {223, 318}. At a nominal operating condition, the RTS 96 power net￾work 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 litera￾ture. 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 vari￾ous interesting network topologies and via Monte Carlo simu￾lations 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 characteriza￾tion 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 con￾tinuous 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) Genera￾tor 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

sufficiently tight. We conjecture that this set is"thin"in an ap- generator dynamics. We envision that our synchronization con- propriate parameter space. Our results suggest that an exact dition enable emerging smart grid ap lications. su ch as p condition for synchronization of any arbitrary network is of the flow optimization subject to stability constraints, distance to fail- form ILTall2oo <c, and we conjecture that the constant c is always ure metric, and the design of control strategies to avoid cascading strictly positive, upper-bounded, and close to 1. However, nportant question not addressed in the present article the region of attraction of a synchronized solution. We Laboratory was that the latter depends on the gap in the presented condition. of the Nation On the application side, the results contained in this nal Labora. ntract DE need to be extended to more detailed power network mode A25396. This ted by National Science Foundation Grants IiS-0904501 and reactive (1893)Oeuvres Completes de Christiaan Huygens (Martinus Nijhoff, The 24. Boccaletti S, Latora V, Moreno Y, chave structure ogical Time(Springer, New York, NY), 2nd Ed del system stability 26 Pai MA (1999)Existence of solutions for the network/load 5. Ha SY, Jeong E, Kang M(2010)Emergent behaviour of a generalized Viscek-type 27. Dobson I (1992)Observations on the geometry of saddle node bifurcation and vol 23:3139-3156. Theory Appl39: 240-243. for synchrony in the firefly pteroptyx 28. cae. J Math Bio 29: 571-58 7. Strogatz SH, Abrams DM, McRobie A, Eckhardt B, ott E(2005)Theoretical mechanics: 29. Hill DI, Chen G(2006)Power systems as dynamic networks. IEEE International Sym- he Millennium Bridge. Nature 438(7064): 43-44. 8. Bennett M, Schatz MF, Rockwood H, wiesenfeld K (2002) Huygens's clocks Proc Math. 30. Dorfler F, Bullo F(2012) Synchronization and transient stability in power networks Phys AM J Contr 0:1616-164 9. Kuramoto oupled non-linear oscillators. 31. llic M(1992)Network theoretic conditions for existence and uniqueness of stead Diego,CApp2821-282810.11091scAs.1992230611 11. Paley DA Leonard NE, Sepulchre R, Grunbaum D, Parrish JK (2007) Osillator models circ sast 29 27903-31982) Steady-state security regions of power systems. / Trans 34. Wu FF, Kumagais(1980)Limits on Power injections for Power Flow Equations to Have Bolle D, Coolen ACC Perez-Vicente C(2001)Coupled dy- Research Laboratory, College of Engineering, University amics of fast spins and slow exchang ions in the XY spin glass. J Phys Math en34:3957-3984 ee N, Barahona M(2004)On the sta e ons. Phys Rev Lett 68(7): 1073-1076. Availableathttp:fieeexplore.ieeeorg/stamp/stamp.jsp?tp=&arnumber=1383983 ic bipolar pop- 15. Tass PA (2003)A model of desynchro ulation networks. Phys Rev E Stat Nonlin Soft Matter Phys 80(6 Pt 2): 066120. eural subpopulations. Bio/ Cybern 89(2): 81-88 JL (2002) Emerging coherence in a population of chemical 3):1676-1678 17.Ko rmentrout GB(1988)Coupled oscillators and the design of central pattern artie graph. SIAM J Appl Dyn Syst 8: 417- 39. Gomez-Gardenes J, Moreno Y, Arenas A(2007)Paths to synchronization on comple 18. Neda Z, Ravasz E, Vicsek T, Brechet Y, Barabasi AL (2000)Physics of the rhythmic ap- ause. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 61(6 Pt B): 6987- 40. Nishikawa T, Motter AE, Lai YC in oscillator synchronize? Phys Rev Lett 91(1): 014101 19. Strogatz SH (2000)From Kuramoto to Crawford: Exploring the onset of synchron 2004)Synchronization of Kuramoto oscillators in scale-free CP, Ritort F, Spigler R(2005)The Kuramoto model: A 42. Zimmerman RD, Muri an D (2011)MATPOWER: Steady-state op ple paradigm for synchronization phenomena. Rev Mod Phys 77: 137-185. ations, planning and analysis tools for power systems res and education. IEEE 21. Dorfler F. Bullo F(2011)On the critical coupling for Kuramoto oscillators. SIAM J Appl Dyn Syst10:10701099. 43. Grigg C et al. ( 1999)The IEEE Reliability Test System-1996 A report prepared by the 22. Hoppensteadt FC Izhikevich EM (1997)Weakly Connected Neural Networks(Springe Methods syst14:1010-1020. 23. Arenas A Diaz-Guilera A, Kurths J, Moreno Y, Zhou C(2008) Synchronization in com- 44. Varaiya PP, Wu FF, Bialek J(2011)Smart operation of smart grid: Risk-limiting rks. Phys Rep 469: 93-153. dispatch. Proc IEEE 99: 40-57. 2010Iwww.pnas.org/cgi/doi/10.1073/pnas.1212134110 orfler et al

sufficiently tight. We conjecture that this set is “thin” in an ap￾propriate parameter space. Our results suggest that an exact condition for synchronization of any arbitrary network is of the form kL†ωkE;∞ < c, and we conjecture that the constant c is always strictly positive, upper-bounded, and close to 1. However, another important question not addressed in the present article concerns the region of attraction of a synchronized solution. We conjecture that the latter depends on the gap in the presented condition. On the application side, the results contained in this paper need to be extended to more detailed power network models, in￾cluding voltage dynamics, reactive power flows, and higher order generator dynamics. We envision that our synchronization con￾ditions enable emerging smart grid applications, such as power flow optimization subject to stability constraints, distance to fail￾ure metric, and the design of control strategies to avoid cascading failures. ACKNOWLEDGMENTS. Research at Los Alamos National Laboratory was carried out under the auspices of the National Nuclear Security Adminis￾tration of the US Department of Energy at Los Alamos National Labora￾tory under Contract DE C52-06NA25396. This material is based, in part, on work supported by National Science Foundation Grants IIS-0904501 and CPS-1135819. 1. Huygens C (1893) Oeuvres Complètes de Christiaan Huygens (Martinus Nijhoff, The Hague, The Netherlands). 2. Strogatz SH (2003) SYNC: The Emerging Science of Spontaneous Order (Hyperion, New York, NY). 3. Winfree AT (2001) The Geometry of Biological Time (Springer, New York, NY), 2nd Ed. 4. Bergen AR, Hill DJ (1981) A structure preserving model for power system stability analysis. IEEE Trans Power Apparatus Syst 100:25–35. 5. Ha SY, Jeong E, Kang MJ (2010) Emergent behaviour of a generalized Viscek-type flocking model. Nonlinearity 23:3139–3156. 6. Ermentrout GB (1991) An adaptive model for synchrony in the firefly pteroptyx malaccae. J Math Biol 29:571–585. 7. Strogatz SH, Abrams DM, McRobie A, Eckhardt B, Ott E (2005) Theoretical mechanics: Crowd synchrony on the Millennium Bridge. Nature 438(7064):43–44. 8. Bennett M, Schatz MF, Rockwood H, Wiesenfeld K (2002) Huygens’s clocks. Proc Math. Phys Eng Sci 458:563–579. 9. Kuramoto Y (1975) Self-entrainment of a population of coupled non-linear oscillators. International Symposium on Mathematical Problems in Theoretical Physics, Lecture Notes in Physics, ed Araki H (Springer), Vol. 39, pp 420–422. 10. Wiesenfeld K, Colet P, Strogatz SH (1998) Frequency locking in Josephson arrays: Connection with the Kuramoto model. Phys Rev E Stat Phys Plasmas Fluids Relat In￾terdiscip Topics 57:1563–1569. 11. Paley DA, Leonard NE, Sepulchre R, Grunbaum D, Parrish JK (2007) Oscillator models and collective motion. IEEE Contr Syst Mag 27:89–105. 12. Jongen G, Anemüller J, Bollé D, Coolen ACC, Perez-Vicente C (2001) Coupled dy￾namics of fast spins and slow exchange interactions in the XY spin glass. J Phys Math Gen 34:3957–3984. 13. Daido H (1992) Quasientrainment and slow relaxation in a population of oscillators with random and frustrated interactions. Phys Rev Lett 68(7):1073–1076. 14. Varela F, Lachaux JP, Rodriguez E, Martinerie J (2001) The brainweb: Phase syn￾chronization and large-scale integration. Nat Rev Neurosci 2(4):229–239. 15. Tass PA (2003) A model of desynchronizing deep brain stimulation with a demand￾controlled coordinated reset of neural subpopulations. Biol Cybern 89(2):81–88. 16. Kiss IZ, Zhai Y, Hudson JL (2002) Emerging coherence in a population of chemical oscillators. Science 296(5573):1676–1678. 17. Kopell N, Ermentrout GB (1988) Coupled oscillators and the design of central pattern generators. Math Biosci 90:87–109. 18. Néda Z, Ravasz E, Vicsek T, Brechet Y, Barabási AL (2000) Physics of the rhythmic ap￾plause. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 61(6 Pt B):6987– 6992. 19. Strogatz SH (2000) From Kuramoto to Crawford: Exploring the onset of synchroni￾zation in populations of coupled oscillators. Physica D 143:1–20. 20. Acebrón JA, Bonilla LL, Vicente CJP, Ritort F, Spigler R (2005) The Kuramoto model: A simple paradigm for synchronization phenomena. Rev Mod Phys 77:137–185. 21. Dörfler F, Bullo F (2011) On the critical coupling for Kuramoto oscillators. SIAM J Appl Dyn Syst 10:1070–1099. 22. Hoppensteadt FC, Izhikevich EM (1997) Weakly Connected Neural Networks (Springer, New York, NY), Vol. 126. 23. Arenas A, Díaz-Guilera A, Kurths J, Moreno Y, Zhou C (2008) Synchronization in com￾plex networks. Phys Rep 469:93–153. 24. Boccaletti S, Latora V, Moreno Y, Chavez M, Hwang DU (2006) Complex networks: Structure and dynamics. Phys Rep 424:175–308. 25. Dörfler F, Bullo F (2012) Exploring synchronization in complex oscillator networks. IEEE Conference on Decision and Control. Available at http://arxiv.org/pdf/1209.1335. pdf. Accessed September 6, 2012. 26. Lesieutre BC, Sauer PW, Pai MA (1999) Existence of solutions for the network/load equations in power systems. IEEE Trans Circ Syst I Fundam Theory Appl 46:1003–1011. 27. Dobson I (1992) Observations on the geometry of saddle node bifurcation and voltage collapse in electrical power systems. IEEE Trans Circ Syst I Fundam Theory Appl 39:240–243. 28. Lavaei J, Tse D, Zhang B (2012) Geometry of power flows in tree networks. IEEE Power and Energy Society General Meeting, 10.1109/PESGM.2012.6344803. 29. Hill DJ, Chen G (2006) Power systems as dynamic networks. IEEE International Sym￾posium on Circuits and Systems, pp 722–725, 10.1109/ISCAS.2006.1692687. 30. Dörfler F, Bullo F (2012) Synchronization and transient stability in power networks and non-uniform Kuramoto oscillators. SIAM J Contr Optim 50:1616–1642. 31. Ilic M (1992) Network theoretic conditions for existence and uniqueness of steady state solutions to electric power circuits. IEEE International Symposium on Circuits and Systems, San Diego, CA pp 2821–2828, 10.1109/ISCAS.1992.230611. 32. Araposthatis A, Sastry S, Varaiya P (1981) Analysis of power-flow equation. Int J Electr Power Energy Syst 3:115–126. 33. Wu F, Kumagai S (1982) Steady-state security regions of power systems. IEEE Trans Circ Syst 29:703–711. 34. Wu FF, Kumagai S (1980) Limits on Power Injections for Power Flow Equations to Have Secure Solutions (Electronics Research Laboratory, College of Engineering, University of California, Berkeley, CA). 35. Jadbabaie A, Motee N, Barahona M (2004) On the stability of the Kuramoto model of coupled nonlinear oscillators. American Control Conference, Boston, MA, pp 4296–4301. Available at http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1383983 &isnumber=30144. Accessed June 30, 2004. 36. Buzna L, Lozano S, Díaz-Guilera A (2009) Synchronization in symmetric bipolar pop￾ulation networks. Phys Rev E Stat Nonlin Soft Matter Phys 80(6 Pt 2):066120. 37. Strogatz SH, Mirollo RE (1988) Phase-locking and critical phenomena in lattices of coupled nonlinear oscillators with random intrinsic frequencies. Physica D 31:143–168. 38. Verwoerd M, Mason O (2009) On computing the critical coupling coefficient for the Kuramoto model on a complete bipartite graph. SIAM J Appl Dyn Syst 8:417–453. 39. Gómez-Gardeñes J, Moreno Y, Arenas A (2007) Paths to synchronization on complex networks. Phys Rev Lett 98(3):034101. 40. Nishikawa T, Motter AE, Lai YC, Hoppensteadt FC (2003) Heterogeneity in oscillator networks: are smaller worlds easier to synchronize? Phys Rev Lett 91(1):014101. 41. Moreno Y, Pacheco AF (2004) Synchronization of Kuramoto oscillators in scale-free networks. Europhys Lett 68(4):603–609. 42. Zimmerman RD, Murillo-Sánchez CE, Gan D (2011) MATPOWER: Steady-state oper￾ations, planning, and analysis tools for power systems research and education. IEEE Trans Power Syst 26:12–19. 43. Grigg C, et al. (1999) The IEEE Reliability Test System—1996. A report prepared by the Reliability Test System Task Force of the Application of Probability Methods Sub￾committee. IEEE Trans Power Syst 14:1010–1020. 44. Varaiya PP, Wu FF, Bialek JW (2011) Smart operation of smart grid: Risk-limiting dispatch. Proc IEEE 99:40–57. 2010 | www.pnas.org/cgi/doi/10.1073/pnas.1212134110 Dörfler et al.

Supporting Information Dorfler et al.10.1073/pna5.1212134110 sI Text matrix LER"Xn is defined by L=diag(Eisai)_)-AIf ng information is organized as follows. is defined component-wise as Bx=l if node k is the sink node SI Mathematical Models and Synchronization Notions provides of edge and as Bk=-l if node k is the source node of edge G a description of the considered coupled oscillator model, in- all other elements are 0. For xER, the vector Bx has com luding a detailed modeling of a mechanical analog and a few ponents x-x for any oriented edge from to i, that is, bmaps ower network models. Furthermore, we state our definition of node variables x; and r to incremental edge variables xi-x. If ynchronization and compare various synchronization conditions diag(aihiilee)is the diagonal matrix of nonzero edge weights, oposed for oscillator networl L=B diagHaihuee)B For a vector xER", the incremental SI Mathematical Analysis of Synchronization provides a rigorous norm lilIl, o max leeki -xl used in the main text can be ex- chronization conditions proposed in the main text. Throughout our is connected, Ker(BT)=Ker(L)=span(In), all n-1 remaining ei nalysis, we provide various examples illustrating certain theoretical genvalues of L are real and strictly positive, and the second smallest cepts and results, and we also compare our results with existing eigenvalue A2()is called the algebraic connectivity. The orthogonal results in the synchronization and power network literature vector spaces Ker(B)and Ker(B)=Im(B )are spanned by SI Robust Synchronization in Presence of Uncertainty extends vectors associated with cycles and cut-sets in the graph(e.g our synchronization condition to the case where the network ref. 1, section 4; ref 2). In the following, we refer to Ker(B)and parameters can vary within prescribed upper and lower bounds. Im(B) as the cycle space and the cut-set space, respectively This parameter-varying approach can account for modeling uI certainties or unmodeled dynamics. Laplacian Inverses. Because the Laplacian matrix L is singular SI Statistical Synchronization Assessment provides a detailed ac we will frequently use its Moore-Penrose pseudo inverse LT.If ount of our Monte Carlo simulation studies and the complex UER"Xn is an orthonormal matrix of eigenvectors of L.the Kuramoto network studies. Throughout this section, we also recall singular value decomposition of L is L=U diag((0,h2 the basics of probability estimation by Monte Carlo methods that and its Moore-Penrose pseudo inverse L is given by Lf=U allow us to establish a statistical synchronization result in a math- diag((0,1/A2,.,1/n y)U. We will frequently use the identit LL'=L-L=ln-flnxn, which follows directly from the singu Finally, the subsection on Sy value dec sition. We also define the effective resistance be- Networks describes the detailed simulation setup for the random- tween nodes i and j by Ri=LI+Li-2Li. We refer to the study ized Institute of Electrical and Electronics Engineers(IEEE)test by Dorfler and Bullo(3)for further information on Laplacian ystems, provides the simulation data used for the dynamic Ieee inverses and on the resistance distance teliability Test System 96(RTS 96) power network simulations, illustrates a dynamic bifurcation scenario in the RTS 96 power SI Mathematical Models and Synchronization Notions etwork, and describes extensions of the results in the main text In this section, we introduce the mathematical model of coupled to variable load demands and load voltages phase oscillators considered in this article, present some syn- SI Preliminaries and notation chronization notions, and give a detailed account of the literature on synchronization of coupled phase oscillators. Vectors and Functions. Let ln and On be the n-dimensional vector of unit and zero entries, and let 1, be the orthogonal comple General Coupled oscillator Model Consider a weighted undirected canonical basis vector of R", that is, the ith entry of e! is 1 and all partitioned node set V=V1UV2, and edge set E induced by the other entries are 0. Let Inxn=ln. I be the(n xn)matrix of unit adjacency matrix AER"X. We assume that the graph G has no ntries. Given an n-tuple (x1, .. x), letxER be the associated self-loops i, il (i.e, ai=0 for all iEv). Associated with this vector. For an ordered index set I of cardinality I and a Id array graph, consider the following model of v1 >0 second-order Exihier, we define diag((xihier)ERIx to be the associated Newtonian and [v2120 first-order kinematic phase oscillators diagonal matrix. For xER, we define the vector-valued func- tions sin(x)=(sin(xu) n(n))and arcsin(x)=(arcsin(x M+D=0-∑asin(-9),i∈v arcsin(nn )) where the arcsin function is defined for the brand (x/2, x/2]. For a set xcR and a matrix AER, let Ax= [s1 D ai sin((4-9),i∈v2, eometry on n-Torus. The set s denotes the unit circle, an angle is a point eES, and an arc is a connected subset of s. The where 0; es and eieR are the phase and frequency of oscil- esic distance between two angles 01, 02ES is the minimum lator iEV, @i ER and Di>0 are the natural frequency and of the counterclockwise and clockwise arc lengths connecting damping coefficient of oscillator iEV, and Mi>0 is the inertial e, and 02. With slight abuse of notation, let 101-02l denote the constant of oscillator iEVi. The coupled oscillator model (SI geodesic distance between two angles 01, 02ES. Finally, the n-tonts, evolves on T" I and features an important symmetry, namel the product set T=S x.xS, is the direct sum of n unit circles. the rotational invariance of the angular variable 6. The interest ng dynamics of the algebraic Graph Theory. Given an undirected, connected, and led oscillator model [sl] arise from a competition between each oscillator,'s tendency to align with weighted graph G(V, E, A)induced by the symmetrical, reduc- its natural frequency o and the synchronization-enforcing cou- ible, and nonnegative adjacency matrix AER, the Laplacian pling ai sin(0i-8i) with its neighbors Dorfleretal.www.pnas.org/cgicontent/short/1212134110 1 of 1

Supporting Information Dörfler et al. 10.1073/pnas.1212134110 SI Text The supporting information is organized as follows. In SI Preliminaries and Notation, we introduce some notation and recall some preliminaries. SI Mathematical Models and Synchronization Notions provides a description of the considered coupled oscillator model, in￾cluding a detailed modeling of a mechanical analog and a few power network models. Furthermore, we state our definition of synchronization and compare various synchronization conditions proposed for oscillator networks. SI Mathematical Analysis of Synchronization provides a rigorous mathematical analysis of synchronization, which leads to the syn￾chronization conditions proposed in the main text. Throughout our analysis, we provide various examples illustrating certain theoretical concepts and results, and we also compare our results with existing results in the synchronization and power network literature. SI Robust Synchronization in Presence of Uncertainty extends our synchronization condition to the case where the network parameters can vary within prescribed upper and lower bounds. This parameter-varying approach can account for modeling un￾certainties or unmodeled dynamics. SI Statistical Synchronization Assessment provides a detailed ac￾count of our Monte Carlo simulation studies and the complex Kuramoto network studies. Throughout this section, we also recall the basics of probability estimation by Monte Carlo methods that allow us to establish a statistical synchronization result in a math￾ematically rigorous way. Finally, the subsection on Synchronization Assessment for Power Networks describes the detailed simulation setup for the random￾ized Institute of Electrical and Electronics Engineers (IEEE) test systems, provides the simulation data used for the dynamic IEEE Reliability Test System 96 (RTS 96) power network simulations, illustrates a dynamic bifurcation scenario in the RTS 96 power network, and describes extensions of the results in the main text to variable load demands and load voltages. SI Preliminaries and Notation Vectors and Functions. Let 1n and 0n be the n-dimensional vector of unit and zero entries, and let 1⊥ n be the orthogonal comple￾ment of 1n in Rn, that is, 1⊥ n ≜ fx ∈ Rn : x ⊥1ng. Let en i be the ith canonical basis vector of Rn, that is, the ith entry of en i is 1 and all other entries are 0. Let 1n× n =1n · 1T n be the ðn × nÞ matrix of unit entries. Given an n-tuple ðx1; ... ; xnÞ, let x∈ Rn be the associated vector. For an ordered index set I of cardinality jIj and a 1D array fxigi∈I , we define diagðfxigi∈I Þ∈ RjIj × jIj to be the associated diagonal matrix. For x ∈ Rn, we define the vector-valued func￾tions sinðxÞ =ðsinðx1Þ; ... ;sinðxnÞÞ and arcsinðxÞ=ðarcsinðx1Þ; ... ; arcsinðxnÞÞ, where the arcsin function is defined for the branch ½−π=2; π=2. For a set X ⊂ Rn and a matrix A∈ Rm × n, let AX = fy∈ Rm : y = Ax; x∈ Xg. Geometry on n-Torus. The set S1 denotes the unit circle, an angle is a point θ ∈S1 , and an arc is a connected subset of S1 . The geodesic distance between two angles θ1,θ2 ∈S1 is the minimum of the counterclockwise and clockwise arc lengths connecting θ1 and θ2. With slight abuse of notation, let jθ1 −θ2j denote the geodesic distance between two angles θ1; θ2 ∈ S1 . Finally, the n-torus, the product set Tn = S1 × ⋯ × S1 , is the direct sum of n unit circles. Algebraic Graph Theory. Given an undirected, connected, and weighted graph GðV; E;AÞ induced by the symmetrical, irreduc￾ible, and nonnegative adjacency matrix A∈ Rn × n, the Laplacian matrix L∈ Rn × n is defined by L= diag
f Pn j=1aijg n i=1 −A. If a number ℓ∈f1; ... ; jEjg and an arbitrary direction are assigned to each edge fi; jg ∈E, the (oriented) incidence matrix B∈ Rn × jEj is defined component-wise as Bkℓ = 1 if node k is the sink node of edge ℓ and as Bkℓ = − 1 if node k is the source node of edge ℓ; all other elements are 0. For x ∈ Rn, the vector BTx has com￾ponents xi −xj for any oriented edge from j to i, that is, BT maps node variables xi and xj to incremental edge variables xi − xj. If diagðfaijgfi;jg∈E Þ is the diagonal matrix of nonzero edge weights, L=B diagðfaijgfi;jg∈E ÞBT. For a vector x ∈ Rn, the incremental norm kxkE;∞ ≜ maxfi;jg∈E jxi −xjj used in the main text can be ex￾pressed via the incidence matrix B as kxkE;∞ = kBTxk∞. If the graph is connected, KerðBTÞ= KerðLÞ=spanð1nÞ, all n− 1 remaining ei￾genvalues of L are real and strictly positive, and the second smallest eigenvalue λ2ðLÞ is called the algebraic connectivity. The orthogonal vector spaces KerðBÞ and KerðBÞ ⊥ = ImðBTÞ are spanned by vectors associated with cycles and cut-sets in the graph (e.g., ref. 1, section 4; ref. 2). In the following, we refer to KerðBÞ and ImðBTÞ as the cycle space and the cut-set space, respectively. Laplacian Inverses. Because the Laplacian matrix L is singular, we will frequently use its Moore–Penrose pseudo inverse L†. If U ∈ Rn × n is an orthonormal matrix of eigenvectors of L, the singular value decomposition of L is L= U diagðf0; λ2; ... ; λngÞUT and its Moore–Penrose pseudo inverse L† is given by L† = U diagðf0; 1=λ2; ... ; 1=λngÞUT. We will frequently use the identity L ·L† = L† · L=In − 1 n 1n × n, which follows directly from the singular value decomposition. We also define the effective resistance be￾tween nodes i and j by Rij = L† ii + L† jj −2L† ij. We refer to the study by Dörfler and Bullo (3) for further information on Laplacian inverses and on the resistance distance. SI Mathematical Models and Synchronization Notions In this section, we introduce the mathematical model of coupled phase oscillators considered in this article, present some syn￾chronization notions, and give a detailed account of the literature on synchronization of coupled phase oscillators. General Coupled Oscillator Model. Consider a weighted, undirected, and connected graph GðV; E;AÞ with n nodes V =f1; ... ; ng, partitioned node set V =V1 ∪ V2, and edge set E induced by the adjacency matrix A∈ Rn × n. We assume that the graph G has no self-loops fi; ig (i.e., aii =0 for all i ∈V). Associated with this graph, consider the following model of jV1j≥ 0 second-order Newtonian and jV2j≥0 first-order kinematic phase oscillators: 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; [S1] where θi ∈S1 and θ _ i ∈ R1 are the phase and frequency of oscil￾lator i∈ V, ωi ∈ R1 and Di >0 are the natural frequency and damping coefficient of oscillator i ∈V, and Mi > 0 is the inertial constant of oscillator i ∈V1. The coupled oscillator model [S1] evolves on Tn × RjV1j and features an important symmetry, namely, the rotational invariance of the angular variable θ. The interest￾ing dynamics of the coupled oscillator model [S1] arise from a competition between each oscillator’s tendency to align with its natural frequency ωi and the synchronization-enforcing cou￾pling aijsinðθi − θjÞ with its neighbors. Dörfler et al. www.pnas.org/cgi/content/short/1212134110 1 of 17

As discussed in the main text, the coupled oscillator model admittance matrix YECx(augmented with the generator [Si] unifies various models proposed in the literature. For ex- transient reactances). If the network is lossless and the voltage ample, for the parameters Vi=0 and D;=l for all iE V2, it levels VAl at all nodes iEVIUV2 are constant, the maximum reduces to the celebrated Kuramoto model (4, 5) real-power transfer between any two nodes i,jEV,UV2 is Vil- Vl- J(Yi), where J(Yi) denotes the susceptance of the (G-),i∈{1,…,n} [S2] transmission line i,jEE. With this notation, the swing dynamics of generator i are given by We refer to specific reviews(6-10)for various theoretical results on the Kuramoto model [S2] and further synchronization ap- M+D=Pm-∑asin(-9) [S4 plications in natural sciences, technology, and social networks Here, we present a detailed modeling of the spring oscillator net- where 0, es and 0 ER are the generator rotor angle and work used as a mechanical analog in the main text and present a frequency; 0 ES for jEv2 are the voltage phase angles at the few power network models, which can be described by the coupled load buses; and Pmi>0, Mi>0, and Di>0 are the mechanical er input from the prime mover, the generator inertia constant, Mechanical Spring Network. Consider the spring network illus- For the load buses D2, we consider the following three load trated in Fig. SI consisting of a group of n particles constrained to models illustrated in Fig. S2. tate around a circle with unit radius. For simplicity, we assume that the particles are allowed to move freely on the circle and 1)Pv buses with frequency-dependent loads: All load buses are exchange their order without collisions. PV buses, that is, the active power demand Pui and the voltage magnitude Vil are specified for each bus. The real power drawn frequency BiER, and its inertial and damping coefficients are by load i consists of a constant term Pii>0 and a frequency Mi>0 and Di>0, respectively. The external forces and torques dependent term D ei with D>0, as illustrated in Fig. S2A. The cting on each particle are()a viscous damping force D; e; opposing esulting real-power balance equation is the direction of motion, (i)a nonconservative force O ER along the direction of motion depicting a preferred natural rotation fre Dn+P1=-∑ ai sin(-9),i∈ quency, and (inl) an elastic restoring torque between interacting particles i andj coupled by an ideal elastic spring with stiffness ai >0 and zero rest length. The topology of the spring network is described The dynamics [S4 and S5] are known as the structure-preserving by the weighted, undirected, and connected graph G(V, E,A power network model (12)and equal the coupled oscillator To compute the elastic torque between the particles, we param- model [Sl] for o=Pmi,i∈v,ando=-Pi,l∈以 etrize the position of each particle i by the unit vector pi= 2)PV buses with constant power loads: All load buses are Pv os(0:), sin(O)]E S cR. The elastic Hookean energy stored buses. each load features a constant real demand Pli>0 the springs is the function E: T-R given up to an additive and the load damping in[S5]is neglected, that is, Di=0 in[S: onstant by The corresponding circuit-theoretical model is shown in Fig S2B. If the angular distances 0i (0)-0i (0)<x/2 are bounded E(6)= for each transmission line iiEE(this condition will be pre- cisely established in the next section), the resulting differen G)) tial-algebraic system has the same local stability properties as the dynamics [S4 and S5 F; see ref. 13. Hence, all our resul lso apply locally to the structure-preserving power network ai (1-cos(0-Bi)) model [S4 and S5] with zero load damping D =0 for iEV2 3)Constant current and constant admittance loads: If each load where we used the trigonometric identity cos(a-B)=cos a cosB+ iEV2 is modeled as a constant current demand li and an(in- ductive)admittance Yi shunt to ground as illustrated in Fig. S2C sin a sin B in the last equality. Hence, we obtain the restoring tor- the linear current-balance equations are I=Y, where IEC que acting on particle i as and VEC are the vectors of nodal current injections and T()=-E()=-∑ ai sin(-9) the form [S3] known as the(lossless) network-redt stem model (14, 15). We refer to the studie Therefore, the network of spring-interconnected particles depicted in Fig. SI obeys the dynamics and Bullo (3)and Sauer and Pai (11) for a detailed of the network-reduced model M+D=-∑qin(-9),i∈{1…ns The above model [s4 and s5 is valid for an Ac grid with a synchronous generator and load models 1-3. We remark that chronous motor loads also assume the form S4 with In conclusion, the spring network in Fig. SI is a mechanical Pmi <0(16), and a first-principle modeling of a DC power source connected to an AC grid via a droop-controlled inverter also re- Power Network Model. The coupled oscillator model [SI] also in- sults in[S5)(further details are provided in ref. 17 er al analog of the coupled oscillator model [SI]with V2=0 cludes a variety of power network models We briefiy present different Remark 1(Voltage Dynamics). To conclude this modeling section, we ower network models compatible with the coupled oscillator model want to state a word of caution regarding the load models 1 and refer to work by Sauer and 7) tailed derivation from a Consider a connected power network with generators Vi and the assumption of constant voltage magnitudes is load buses V2. The network is described by the symmetrical nodal because voltage magnitudes are controlled at the generators and the Dorfleretal.www.pnas.org/cgicontent/short/1212134110 2 of 1

As discussed in the main text, the coupled oscillator model [S1] unifies various models proposed in the literature. For ex￾ample, for the parameters V1 = ™ and Di =1 for all i ∈V2, it reduces to the celebrated Kuramoto model (4, 5) _ θi = ωi − Xn j=1 aij sin
θi −θj ; i ∈f1; ... ; ng: [S2] We refer to specific reviews (6–10) for various theoretical results on the Kuramoto model [S2] and further synchronization ap￾plications in natural sciences, technology, and social networks. Here, we present a detailed modeling of the spring oscillator net￾work used as a mechanical analog in the main text and present a few power network models, which can be described by the coupled oscillator model [S1]. Mechanical Spring Network. Consider the spring network illus￾trated in Fig. S1 consisting of a group of n particles constrained to rotate around a circle with unit radius. For simplicity, we assume that the particles are allowed to move freely on the circle and exchange their order without collisions. Each particle is characterized by its phase angle θi ∈ S1 and frequency θ _ i ∈ R, and its inertial and damping coefficients are Mi >0 and Di >0, respectively. The external forces and torques acting on each particle are (i) a viscous damping force Diθ _ i opposing the direction of motion, (ii) a nonconservative force ωi ∈ R along the direction of motion depicting a preferred natural rotation fre￾quency, and (iii) an elastic restoring torque between interacting particles i and j coupled by an ideal elastic spring with stiffness aij >0 and zero rest length. The topology of the spring network is described by the weighted, undirected, and connected graph GðV; E;AÞ. To compute the elastic torque between the particles, we param￾etrize the position of each particle i by the unit vector pi = ½cosðθiÞ;sinðθiÞT ∈ S1 ⊂ R2 . The elastic Hookean energy stored in the springs is the function E : Tn → R given up to an additive constant by EðθÞ= X fi;jg∈E aij 2   pi−pj  2 2 = X fi;jg∈E aij
1 −cosðθiÞcos
θj − sinðθiÞsin
θj = X fi;jg∈E aij
1 −cos
θi −θj ; where we used the trigonometric identity cosðα−βÞ=cos α cos β + sin α sin β in the last equality. Hence, we obtain the restoring tor￾que acting on particle i as TiðθÞ = − ∂ ∂θi EðθÞ= − Xn j=1 aij sin
θi − θj : Therefore, the network of spring-interconnected particles depicted in Fig. S1 obeys the dynamics Miθ € i + Diθ _ i = ωi − Xn j=1 aij sin
θi −θj ; i∈ f1; ... ; ng: [S3] In conclusion, the spring network in Fig. S1 is a mechanical analog of the coupled oscillator model [S1] with V2 = ™. Power Network Model. The coupled oscillator model [S1] also in￾cludes a variety of power networkmodels.We briefly present different power network models compatible with the coupled oscillator model [S1] and refer to work by Sauer and Pai (ref. 11, chapter 7) for a de￾tailed derivation from a higher order first-principle model. Consider a connected power network with generators V1 and load buses V2. The network is described by the symmetrical nodal admittance matrix Y ∈ ℂn × n (augmented with the generator transient reactances). If the network is lossless and the voltage levels jVij at all nodes i∈ V1 ∪ V2 are constant, the maximum real-power transfer between any two nodes i; j∈V1 ∪ V2 is aij = jVij · jVjj · IðYijÞ, where IðYijÞ denotes the susceptance of the transmission line fi; jg∈ E. With this notation, the swing dynamics of generator i are given by Mi €θi + Diθ _ i =Pm;i − Xn j=1 aij sin
θi − θj ; i∈ V1; [S4] where θi ∈ S1 and θ _ i ∈ R1 are the generator rotor angle and frequency; θj ∈ S1 for j∈ V2 are the voltage phase angles at the load buses; and Pm;i >0, Mi >0, and Di >0 are the mechanical power input from the prime mover, the generator inertia constant, and the damping coefficient, respectively. For the load buses V2, we consider the following three load models illustrated in Fig. S2. 1) PV buses with frequency-dependent loads: All load buses are PV buses, that is, the active power demand Pl;i and the voltage magnitude jVij are specified for each bus. The real power drawn by load i consists of a constant term Pl;i > 0 and a frequency￾dependent term Diθ _ i with Di > 0, as illustrated in Fig. S2A. The resulting real-power balance equation is Diθ _ i +Pl;i = − Xn j=1 aij sin
θi −θj ; i∈ V2: [S5] The dynamics [S4 and S5] are known as the structure-preserving power network model (12) and equal the coupled oscillator model [S1] for ωi =Pm;i, i ∈V1, and ωi = −Pl;i, i∈ V2. 2) PV buses with constant power loads: All load buses are PV buses, each load features a constant real-power demand Pl;i >0, and the load damping in [S5] is neglected, that is, Di = 0 in [S5]. The corresponding circuit-theoretical model is shown in Fig. S2B. If the angular distances jθiðtÞ−θjðtÞj<π=2 are bounded for each transmission line fi; jg∈ E (this condition will be pre￾cisely established in the next section), the resulting differen￾tial-algebraic system has the same local stability properties as the dynamics [S4 and S5]; see ref. 13. Hence, all our results also apply locally to the structure-preserving power network model [S4 and S5] with zero load damping Di =0 for i∈ V2. 3) Constant current and constant admittance loads: If each load i ∈V2 is modeled as a constant current demand Ii and an (in￾ductive) admittance Yi;shunt to ground as illustrated in Fig. S2C, the linear current-balance equations are I = YV, where I ∈ ℂn and V ∈ℂn are the vectors of nodal current injections and voltages. After elimination of the bus variables Vi, i∈ V2, through Kron reduction [S3], the resulting dynamics assume the form [S3] known as the (lossless) network-reduced power system model (14, 15). We refer to the studies by Dörfler and Bullo (3) and Sauer and Pai (11) for a detailed derivation of the network-reduced model. The above model [S4 and S5] is valid for an AC grid with a synchronous generator and load models 1–3. We remark that synchronous motor loads also assume the form [S4] with Pm;i <0 (16), and a first-principle modeling of a DC power source connected to an AC grid via a droop-controlled inverter also re￾sults in [S5] (further details are provided in ref. 17). Remark 1 (Voltage Dynamics). To conclude this modeling section, we want to state a word of caution regarding the load models. The PV load models [S1 and S2] assume constant voltage magnitudes jVij at the loads. Under normal operating conditions, the assumption of constant voltage magnitudes is well justified because voltage magnitudes are controlled at the generators and the Dörfler et al. www.pnas.org/cgi/content/short/1212134110 2 of 17

active power flow aisin(0i-0i)=Vil. l- 3(Y). sin(0i -ei) be- when the coupling dominates the dissimilarity. Various con- ween two nodes i, jEVIU v2 is primarily govemed by the angular ditions are proposed in the power systems and synchronization diference 0-B, and not by the voltage magnitudes Vil, Vil. The literature to quantify this tradeoff between coupling and dissimi- latter assumption is known as the "decoupling assumption"in the larity. The coupling is typically quantified by the algebraic con- power systems community. Whereas the model in (S4 and S5] is well- nectivity A2(L)(15, 19-23)or by the weighted nodal degree adopted for power sy ai (3, 15, 24-26), and the dissimilarity is quantified load voltage magnitudes ceases to hold in a heavily stressed gnd (near by either absolute norms loll, or incremental(relative)norms a bifurcation point), where additional dynamic phenomena can occu B ollo, where typically, PE 12, oo .Sometimes, these condition such as voltage collapse at the loads( 18). In short, the coupling weights can be evaluated only numerically because they are state-dependent "Likewise, if the shunt admitance loads in the load model(S3 (% 24)or arise from a nontrivial linearization process, such as the Ire not constant(.g, constant power loads can be transformed cise and accurate results are only known for specific topologies, woltage-dependent shunt admittances), the Kron reduction process may, such as complete graphs(9, 28), linear chains(29, 30), and com- be ill-posed or the admittance matri of the network-reduced model depends on the load voltages. In the latter case, the coupling weights aj plete bipartite graphs (31)with uniform weights. are again not constant but depend on the load voltage For arbitrary coupling topologies, only sufficient conditions are To account explicitly for such unmodeled voltage dynamics af- random networks(21, 32-34). To best of our t stigations for known(15, 19, 20, 24), as well as numerical inv sharpest and provably correct synchronization conditions for nization in the Presence of uncertainty bitrary top ref. 10, theorem 4.7). For arbitrary undirected, connected Synchronization Notions. The following subsets of the n-torus T" and weighted graphs G(V, E, A), simulation studies indicate that the are essential for the synchronization problem: ForyE[0, r/2 , let known sufficient conditions(15, 19, 20, 24)are conservative esti operty 10i - sy for (i,jEE. Also, let AG(r)be the interior review article on synchronization concludes with the open problem f△c(y) of finding sharp synchronization conditions(6, 7, 9, 22, 23, 35). ition 1. A solution (0, 0): R20-(T, RIV) to the coupled Mathematical Analysis of Synchronization illator model[Sl] is said to be synchronized if Bo E AG(r)and This section presents our analysis of the synchronization problem there exists sync ER such that e(r=0o+ r(mod 2r)and in the coupled oscillator model [S1] e(()=@syne vl for all 120 In other words, here, synchronized trajectories have the prop- Algebraic Approach to Synchronization. Here, we present a pre- erties of frequency synchronization and phase cohesiveness, that viously unexplored analysis approach that reduces the synchro is, all oscillators rotate with the same synchronization frequency nization problem to an equivalent algebraic problem that reveals ync, and all their phases belong to the set AG(r). For a power the crucial role of cycles and cut-sets in the graph topology.In network model [S4 and S5), the notion of phase cohesiveness is a first analysis step, we reduce the synchronization problem for equivalent to bounded flows aisin(B-i)lsa sin(r)for all the coupled oscillator model [Sl] to a simpler problem,namely transmission lines i,JE&. For the coupled oscillator model [S1], the explicit synchro- stability of a first-order model. It turns out that existence and nization frequency is given by async 4 2i=10:/Li= Di(a detailed pled oscillator model (Si can be entirely described by means of the first-order Kuramoto model [S2 ency obtain ayne=0(or, equivalently, oel) corresponding to Lemma 1(Synchronization Equivalence). Consider the coupled os- work applicat Given a point rES and an angle sE[ 0, 2x], let rot (r)ES be the rotation of r counterclockwise by the angle s. For i)[e] is a locally exponentially stable synchronization manifold of (n1,…,mn)∈T, define the equivalence class the Kuramoto model (s2 (1,…,n)={(rot(r1),…,rot(rn)∈T∈[0.2x]} ii)(0. 0n is a locally exponentially stable synchronization mani Clearly,f(r1,…,mn)∈△c(r),then[n,…,mn)c△c(7) If the equivalent statements (i)and (ii) are true, then locally Definition 2. Given eeAG() for some rE (0, x/2 the set near, their respective synchronization manifolds,the (l0 o1)CTXRIVl is a synchronization manifold of the coupled oscillator model [SI] and the Kuramoto model[S2]. Note that a synchronized solution takes value in a synchro mosely speaking, the topological conjugacy result means that oscillators [S2), the state space T, the set AG(x /2), and the the trajectories of the two plots in Fig. S4 can be continuously ynchronization manifold g associated with an angle array deformed to match each other while preserving parameterization of g=(0, 0)eT are illustrated in Fig. $3 time. Lemma I is illustrated in Fig. S4, and its proof can be found in the study by Dorfer and Bullo(ref 9, theorems 5.1 and 5.3). Existing Synchronization Conditions. The coupled oscillator dy By Lemma I, the local synchronization problem for the cou- amics[S1], and the Kuramoto dynamics[S2]for that matter, pled oscillator model [SI] reduces to the synchronization feature(i)the synchronizing effect of the coupling described by problem for the first-order Kuramoto model [S2). Henceforth, we the weighted edges of the graph G(V, E, A)and() the de- restrict ourselves to the Kuramoto model [S2 The following g effect of the dissimilar natural frequencies result is known in the synchronization literature(15, 20)as well oEI at the nodes. Loosely speaking, synchronization occurs as in power systems, where the saturation of a transmission line Dorfleretal.www.pnas.org/cgicontent/short/1212134110 3 of 1

active power flow aijsinðθi − θjÞ= jVij · jVjj · IðYijÞ·sinðθi − θjÞ be￾tween two nodes i; j∈V1∪ V2 is primarily governed by the angular difference θi −θj and not by the voltage magnitudes jVij; jVjj. The latter assumption is known as the “decoupling assumption” in the power systems community. Whereas the model in [S4 and S5] is well￾adopted for power systems stability studies, the assumption of constant load voltage magnitudes ceases to hold in a heavily stressed grid (near a bifurcation point), where additional dynamic phenomena can occur, such as voltage collapse at the loads (18). In short, the coupling weights aij are not necessarily constant. Likewise, if the shunt admittance loads in the load model [S3] are not constant (e.g., constant power loads can be transformed to voltage-dependent shunt admittances), the Kron reduction process may be ill-posed or the admittance matrix of the network-reduced model depends on the load voltages. In the latter case, the coupling weights aij are again not constant but depend on the load voltages. To account explicitly for such unmodeled voltage dynamics af￾fecting the coupling weights aij, we study the coupled oscillator model [S1] with interval-valued parameters in SI Robust Synchro￾nization in the Presence of Uncertainty. Synchronization Notions. The following subsets of the n-torus Tn are essential for the synchronization problem: For γ ∈ ½0; π=2½, let ΔGðγÞ ⊂Tn be the closed set of angle arrays ðθ1; ... ; θnÞ with the property jθi − θjj≤ γ for fi; jg∈ E. Also, let ΔGðγÞ be the interior of ΔGðγÞ. Definition 1. A solution ðθ; θ _Þ : R≥0 → ðTn; RjV1j Þ to the coupled oscillator model [S1] is said to be synchronized if θ0 ∈ ΔGðγÞ and there exists ωsync ∈ R1 such that θðtÞ= θ0 + ωsync1ntðmod 2πÞ and θ _ðtÞ =ωsync1jV1j for all t≥0. In other words, here, synchronized trajectories have the prop￾erties of frequency synchronization and phase cohesiveness, that is, all oscillators rotate with the same synchronization frequency ωsync, and all their phases belong to the set ΔGðγÞ. For a power network model [S4 and S5], the notion of phase cohesiveness is equivalent to bounded flows jaijsinðθi −θjÞj≤aijsinðγÞ for all transmission lines fi; jg ∈E. For the coupled oscillator model [S1], the explicit synchro￾nization frequency is given by ωsync ≜ Pn i=1ωi= Pn i=1Di (a detailed derivation is provided in ref. 9). By transforming to a rotating frame with frequency ωsync and by replacing ωi by ωi − Diωsync, we obtain ωsync =0 (or, equivalently, ω∈1⊥ n ) corresponding to balanced power injections P i∈V1 Pm;i =P i∈V2 Pl;i in power net￾work applications. Hence, without loss of generality, we assume that ω ∈1⊥ n such that ωsync =0. Given a point r ∈ S1 and an angle s∈½0; 2π, let rotsðrÞ ∈S1 be the rotation of r counterclockwise by the angle s. For ðr1; ... ; rnÞ ∈Tn, define the equivalence class ½ðr1; ... ; rnÞ= fðrotsðr1Þ; ... ;rotsðrnÞÞ∈ Tnjs∈½0; 2πg: Clearly, if ðr1; ... ; rnÞ∈ ΔGðγÞ, then ½ðr1; ... ; rnÞ ⊂ ΔGðγÞ. Definition 2. Given θ ∈ ΔGðγÞ for some γ ∈½0; π=2½, the set ð½θ; 0jV1jÞ⊂ Tn × RjV1j is a synchronization manifold of the coupled oscillator model [S1]. Note that a synchronized solution takes value in a synchroni￾zation manifold due to rotational symmetry. For two first-order oscillators [S2], the state space T2, the set ΔGðπ=2Þ, and the synchronization manifold ½θ* associated with an angle array θ* = ðθ* 1; θ* 2Þ∈ T2 are illustrated in Fig. S3. Existing Synchronization Conditions. The coupled oscillator dy￾namics [S1], and the Kuramoto dynamics [S2] for that matter, feature (i) the synchronizing effect of the coupling described by the weighted edges of the graph GðV; E;AÞ and (ii) the de￾synchronizing effect of the dissimilar natural frequencies ω∈1⊥ n at the nodes. Loosely speaking, synchronization occurs when the coupling dominates the dissimilarity. Various con￾ditions are proposed in the power systems and synchronization literature to quantify this tradeoff between coupling and dissimi￾larity. The coupling is typically quantified by the algebraic con￾nectivity λ2ðLÞ (15, 19–23) or by the weighted nodal degree degi ≜Pn j=1aij (3, 15, 24–26), and the dissimilarity is quantified by either absolute norms kωkp or incremental (relative) norms kBTωkp, where, typically, p∈ f2; ∞g. Sometimes, these conditions can be evaluated only numerically because they are state-dependent (19, 24) or arise from a nontrivial linearization process, such as the Master stability function formalism (22, 23, 27). In general, con￾cise and accurate results are only known for specific topologies, such as complete graphs (9, 28), linear chains (29, 30), and com￾plete bipartite graphs (31) with uniform weights. For arbitrary coupling topologies, only sufficient conditions are known (15, 19, 20, 24), as well as numerical investigations for random networks (21, 32–34). To best of our knowledge, the sharpest and provably correct synchronization conditions for ar￾bitrary topologies assume the form λ2ðLÞ>
P fi;jg∈E jωi −ωjj 2 1=2 (ref. 10, theorem 4.7). For arbitrary undirected, connected, and weighted graphs GðV; E;AÞ, simulation studies indicate that the known sufficient conditions (15, 19, 20, 24) are conservative esti￾mates on the threshold from incoherence to synchrony, and every review article on synchronization concludes with the open problem of finding sharp synchronization conditions (6, 7, 9, 22, 23, 35). Mathematical Analysis of Synchronization This section presents our analysis of the synchronization problem in the coupled oscillator model [S1]. Algebraic Approach to Synchronization. Here, we present a pre￾viously unexplored analysis approach that reduces the synchro￾nization problem to an equivalent algebraic problem that reveals the crucial role of cycles and cut-sets in the graph topology. In a first analysis step, we reduce the synchronization problem for the coupled oscillator model [S1] to a simpler problem, namely, stability of a first-order model. It turns out that existence and local exponential stability of synchronized solutions of the cou￾pled oscillator model [S1] can be entirely described by means of the first-order Kuramoto model [S2]. Lemma 1 (Synchronization Equivalence). Consider the coupled os￾cillator model [S1] and the Kuramoto model [S2]. The following statements are equivalent for any γ ∈½0; π=2½ and any synchro￾nization manifold ð½θ; 0jV1jÞ⊂ ΔGðγÞ × RjV1j : i) ½θ is a locally exponentially stable synchronization manifold of the Kuramoto model [S2]. ii) ð½θ; 0jV1jÞ is a locally exponentially stable synchronization mani￾fold of the coupled oscillator model [S1]. If the equivalent statements (i) and (ii) are true, then locally near their respective synchronization manifolds, the coupled oscillator model [S1] and the Kuramoto model [S2], together with the frequency dynamics d dtθ _= − M−1Dθ _, are topologically conjugate. Loosely speaking, the topological conjugacy result means that the trajectories of the two plots in Fig. S4 can be continuously deformed to match each other while preserving parameterization of time. Lemma 1 is illustrated in Fig. S4, and its proof can be found in the study by Dörfler and Bullo (ref. 9, theorems 5.1 and 5.3). By Lemma 1, the local synchronization problem for the cou￾pled oscillator model [S1] reduces to the synchronization problem for the first-order Kuramoto model [S2]. Henceforth, we restrict ourselves to the Kuramoto model [S2]. The following result is known in the synchronization literature (15, 20) as well as in power systems, where the saturation of a transmission line Dörfler et al. www.pnas.org/cgi/content/short/1212134110 3 of 17

corresponds to a singularity of the load flow Jacobian, resulting The following statements hold in a saddle-node bifurcation(12, 13, 18, 19, 24, 36-42) 1)Absolute boundedness: If there exists a synchronized solution emma 2(Stable Synchronization in(AGx/2). Consider the kuramoto model (S2] with a connected graph G(V, E, A), and let yE 0,x/2 deg sin(y)≥ oil for all i∈{1,…,n} The following statements hold. 1)Jacobian: The Jacobian of the Kuramoto model evaluated at tion BEAG(r) 2)Incremental boundedness: If there exists a synchronized solu- J(=-Bdg({a(-9)s) (cg+dg)sn)2-tra∈:. 2)Stability: If there exists an equilibrium point g"EAG(), it be- Proof: The first condition arises becat ∈ sin(n), sin(n) for BEAG(r), and the fixed-point equation [S6 has no solution if condition (S8 is not satisfied. 3)Uniqueness: This old is unique in△c(x/2) Alternatively, because @El#, multiplication of the fixed-point equation [S7 by the vector (e-e)El for i,JE& or Proof: Because we have that(a-∑=1 aik sin(日1-)= juivalently, subtraction of the ith and jth fixed-point equation -∑k=1 aik cos(2- aik sin( -Ox))=ai [S6) yields the following equation for all (i,jleE: cos(Oi-0i), the negative Jacobian of the right-hand side of the Kuramoto model (S2] equals the Laplacian matrix of the connected graph G(V, E, A), where aj=ai cos(0 a-ai=>(aik sin(0;-0k)-ajx sin(@-0x)).[S101 Equivalently, in compact notation, the Jacobian is given (0)=-Bdiag(aj cos(0:-Bi)ijlec)b. Th ent Again, [S10] has no solution in AG(r) if condition [S9 is not The Jacobian J(O)evaluated at an equilibrium point 0*EAG(y) sat is negative semidefinite with rank n-1. Its nullspace is In and In the following, we aim to find sufficient and sharp e rises from the rotational symmetry of the right-hand side of the under which the fixed-point equation [S7] admits Kuramoto model [S2](an illustration is provided in Fig. S3). BEAG(r). We resort to a rather straightforward solut Consequently, the equilibrium point B EAG(r) is locally(trans- By formally replacing each term sin(ei-0i)in the sally) exponentially stable. Moreover, the corresponding equ quations [S7] by an auxiliary scalar variable vi librium manifold [e]E AG(r) is locally exponentially stable. This point equation [S7] is equivalently written as completes the proof of statement 2. by f: T -e statement 3, we denote ght-hand side of [S2 =B diag((a,)yer)y f()=-∑4小in(-),i∈{1…,n y=sin( B'e [S12] where yEREl is a vector with elements yi. We will refer to In ref. 39(corollary 1), it is shown that f- o Is a one-to-one [Sil] as the auxiliary fixed-point equation and characterize its unction on AG(r/2)modulo rotational symmetry, that is, for properties in the following theorem fu e1EAG(/2)and z e 4g(/ 2), we have that f(en)=f()if and Theorem 1(Properties of the Fixed-Point Equations). Consider the only if[01]=J02]. This proves uniqueness of the equilibrium Kuramoto model (S2) with graph G(V, E, A)and oE1#, its fixed manifold in△c(x/2) By Lemma 2, the problem of finding a locally stable syn- point equation [S7], and the auxiliary fixed-point equation[Sil] ch eEAG()for some yE(0, r /2[. The fixed-point equations of 1) Exact solution: Every solution of the auxiliary fixed-point equa- tion [Sil] is of the form y=BL o+who [s13 wi- S6] where the homogeneous solution WhomER satisfies In a compact notation, the fixed quations [S6] are diag(le 2) nization condition: Let yE0, r/2[ The followi a=B diag(ag)wilee)sin(e) three statements are equivalent. [S71 () There exists a solution 0 AG(r) to the fixed-point equation The following conditions show that the nat uencies o have i) There exists a solution e∈△c()to to be absolutely and incrementally bounded has to be sufficiently large such that fixed points of [S6] exist. BLo+Whom=sin(B'e) [SI4] yE0, r/2 and define the weighted nodal degree deg, A tion [Si1] of the form [S13] satisfying the norm constraint for each node i∈{1,…,n} Iyll o ssin(r)and the cycle constraint arcsin()EIm(B) Dorfleretal.www.pnas.org/cgicontent/short/1212134110 4 of 1

corresponds to a singularity of the load flow Jacobian, resulting in a saddle-node bifurcation (12, 13, 18, 19, 24, 36–42). Lemma 2 (Stable Synchronization in (ΔGπ=2)). Consider the Kuramoto model [S2] with a connected graph GðV; E;AÞ, and let γ ∈ ½0; π=2½. The following statements hold: 1) Jacobian: The Jacobian of the Kuramoto model evaluated at θ ∈Tn is given by JðθÞ = −B diag aijcos
θi −θj  fi;jg∈E  BT: 2) Stability: If there exists an equilibrium point θ* ∈ ΔGðγÞ, it be￾longs to a locally exponentially stable equilibrium manifold ½θ*∈ ΔGðγÞ. 3) Uniqueness: This equilibrium manifold is unique in ΔGðπ=2Þ. Proof: Because we have that ∂ ∂θi ðωi −Pn k=1aik sinðθi −θkÞÞ= −Pn k=1aik cosðθi − θkÞ and ∂ ∂θj ðωi −Pn k=1aik sinðθi − θkÞÞ= aij cosðθi −θjÞ, the negative Jacobian of the right-hand side of the Kuramoto model [S2] equals the Laplacian matrix of the connected graph GðV; E;A~Þ, where ~aij =aijcosðθi − θjÞ. Equivalently, in compact notation, the Jacobian is given by JðθÞ= −B diagðfaijcosðθi−θjÞgfi;jg∈E ÞBT. This completes the proof of statement 1. The Jacobian JðθÞ evaluated at an equilibrium point θ* ∈ ΔGðγÞ is negative semidefinite with rank n−1. Its nullspace is 1n and arises from the rotational symmetry of the right-hand side of the Kuramoto model [S2] (an illustration is provided in Fig. S3). Consequently, the equilibrium point θ* ∈ ΔGðγÞ is locally (trans￾versally) exponentially stable. Moreover, the corresponding equi￾librium manifold ½θ*∈ ΔGðγÞ is locally exponentially stable. This completes the proof of statement 2. To prove statement 3, we denote the right-hand side of [S2] by f : Tn → Rn, where f is defined component-wise by fiðθÞ= ωi − Xn j=1 aijsin
θi − θj ; i ∈f1; ... ; ng: In ref. 39 (corollary 1), it is shown that f −ω is a one-to-one function on ΔGðπ=2Þ modulo rotational symmetry, that is, for θ1 ∈ ΔGðπ=2Þ and θ2 ∈ ΔGðπ=2Þ, we have that fðθ1Þ = fðθ2Þ if and only if ½θ1 =½θ2. This proves uniqueness of the equilibrium manifold in ΔGðπ=2Þ. ■ By Lemma 2, the problem of finding a locally stable syn￾chronization manifold reduces to that of finding a fixed point θ* ∈ ΔGðγÞ for some γ ∈ ½0; π=2½. The fixed-point equations of the Kuramoto model [S2] read as follows: ωi = Xn j=1 aijsin
θi − θj ; i ∈f1; ... ; ng: [S6] In a compact notation, the fixed-point equations [S6] are ω =B diag aij fi;jg∈E  sin
BTθ : [S7] The following conditions show that the natural frequencies ω have to be absolutely and incrementally bounded and the nodal degree has to be sufficiently large such that fixed points of [S6] exist. Lemma 3 (Necessary Synchronization Conditions). Consider the Kuramoto model [S2] with graph GðV; E; AÞ and ω∈ 1⊥ n . Let γ ∈½0; π=2½, and define the weighted nodal degree degi ≜Pn j=1aij for each node i∈ f1; ... ; ng. The following statements hold: 1) Absolute boundedness: If there exists a synchronized solution θ ∈ ΔGðγÞ, degi sinðγÞ ≥jωij for all i ∈f1; ... ; ng: [S8] 2) Incremental boundedness: If there exists a synchronized solu￾tion θ ∈ ΔGðγÞ,  degi + degj  sinðγÞ ≥jωi −ωjj for all fi; jg ∈E: [S9] Proof: The first condition arises because sinðθi − θjÞ∈ ½−sinðγÞ;sinðγÞ for θ ∈ ΔGðγÞ, and the fixed-point equation [S6] has no solution if condition [S8] is not satisfied. Alternatively, because ω∈1⊥ n , multiplication of the fixed-point equation [S7] by the vector ðen i −en j Þ ∈1⊥ n for fi; jg∈E or, equivalently, subtraction of the ith and jth fixed-point equation [S6] yields the following equation for all fi; jg ∈E: ωi −ωj = Xn k=1
aik sinðθi − θkÞ− ajk sin
θj − θk : [S10] Again, [S10] has no solution in ΔGðγÞ if condition [S9] is not satisfied. ■ In the following, we aim to find sufficient and sharp conditions under which the fixed-point equation [S7] admits a solution θ* ∈ ΔGðγÞ. We resort to a rather straightforward solution ansatz. By formally replacing each term sinðθi −θjÞ in the fixed-point equations [S7] by an auxiliary scalar variable ψij, the fixed￾point equation [S7] is equivalently written as ω =B diag aij fi;jg∈E  ψ; [S11] ψ =sin
BTθ ; [S12] where ψ ∈ RjEj is a vector with elements ψij. We will refer to [S11] as the auxiliary fixed-point equation and characterize its properties in the following theorem. Theorem 1 (Properties of the Fixed-Point Equations). Consider the Kuramoto model [S2] with graph GðV; E;AÞ and ω∈ 1⊥ n , its fixed point equation [S7], and the auxiliary fixed-point equation [S11]. The following statements hold: 1) Exact solution: Every solution of the auxiliary fixed-point equa￾tion [S11] is of the form ψ = BTL† ω +ψhom; [S13] where the homogeneous solution ψhom ∈ RjEj satisfies diagðfaijgfi;jg∈E Þ ψhom ∈ KerðBÞ. 2) Exact synchronization condition: Let γ ∈½0; π=2½. The following three statements are equivalent: (i) There exists a solution θ* ∈ ΔGðγÞ to the fixed-point equation [S7]. (ii) There exists a solution θ ∈ ΔGðγÞ to BTL† ω +ψhom = sin
BTθ [S14] for some ψhom ∈diagðf1=aijgfi; jg∈E ÞkerðBÞ. (iii) There exists a solution ψ ∈ RjEjto the auxiliary fixed-point equa￾tion [S11] of the form [S13] satisfying the norm constraint kψk∞ ≤sinðγÞ and the cycle constraint arcsinðψÞ∈ImðBTÞ. Dörfler et al. www.pnas.org/cgi/content/short/1212134110 4 of 17