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 appropriate 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, including voltage dynamics, reactive power flows, and higher order generator dynamics. We envision that our synchronization conditions enable emerging smart grid applications, such as power flow optimization subject to stability constraints, distance to failure 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 Administration of the US Department of Energy at Los Alamos National Laboratory 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 Interdiscip 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 dynamics 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 synchronization and large-scale integration. Nat Rev Neurosci 2(4):229–239. 15. Tass PA (2003) A model of desynchronizing deep brain stimulation with a demandcontrolled 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 applause. 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 synchronization 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 complex 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 Symposium 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 population 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 operations, 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 Subcommittee. 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, including 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 synchronization 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 uncertainties or unmodeled dynamics. SI Statistical Synchronization Assessment provides a detailed account 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 mathematically rigorous way. Finally, the subsection on Synchronization Assessment for Power Networks describes the detailed simulation setup for the randomized 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 complement 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 functions 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, irreducible, 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 components 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 expressed via the incidence matrix B as kxkE;∞ = kBTxk∞. If the graph is connected, KerðBTÞ= KerðLÞ=spanð1nÞ, all n− 1 remaining eigenvalues 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 between 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 synchronization 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 oscillator 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 interesting 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 coupling 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 example, 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 applications in natural sciences, technology, and social networks. Here, we present a detailed modeling of the spring oscillator network 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 illustrated 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 frequency, 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 parametrize 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 torque 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 includes 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 detailed 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 frequencydependent 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 precisely established in the next section), the resulting differential-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 (inductive) 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 results 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