
==== Front
Nat Commun
Nat Commun
Nature Communications
2041-1723
Nature Publishing Group UK London

33674557
21290
10.1038/s41467-021-21290-5
Article
Asymmetry underlies stability in power grids
Molnar Ferenc 13
http://orcid.org/0000-0002-2147-0242
Nishikawa Takashi t-nishikawa@northwestern.edu

12
http://orcid.org/0000-0003-1794-4828
Motter Adilson E. 12
1 grid.16753.36 0000 0001 2299 3507 Department of Physics and Astronomy, Northwestern University, Evanston, IL USA
2 https://ror.org/000e0be47 grid.16753.36 0000 0001 2299 3507 Northwestern Institute on Complex Systems, Northwestern University, Evanston, IL USA
3 Present Address: SimpleRose Inc, 1017 Olive Street, Suite 800, Saint Louis, MO 63101 USA
5 3 2021
5 3 2021
2021
12 145726 8 2020
15 1 2021
© The Author(s) 2021, corrected publication 2024
2021
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Behavioral homogeneity is often critical for the functioning of network systems of interacting entities. In power grids, whose stable operation requires generator frequencies to be synchronized—and thus homogeneous—across the network, previous work suggests that the stability of synchronous states can be improved by making the generators homogeneous. Here, we show that a substantial additional improvement is possible by instead making the generators suitably heterogeneous. We develop a general method for attributing this counterintuitive effect to converse symmetry breaking, a recently established phenomenon in which the system must be asymmetric to maintain a stable symmetric state. These findings constitute the first demonstration of converse symmetry breaking in real-world systems, and our method promises to enable identification of this phenomenon in other networks whose functions rely on behavioral homogeneity.

Converse symmetry breaking is a counterintuitive phenomenon in which the system must have an asymmetry to stabilize a symmetric state. Molnar et al. demonstrate this effect in real power-grid networks and show that synchronous operation can be improved by inhomogeneities across power generators.

Subject terms

Energy grids and networks
Complex networks
Nonlinear phenomena
https://doi.org/10.13039/100006133 DOE | Advanced Research Projects Agency - Energy (Advanced Research Projects Agency - Energy - U.S. Department of Energy) DE-AR0000702 DE-AR0000702 DE-AR0000702 Molnar Ferenc Nishikawa Takashi Motter Adilson E. Northwestern University’s Finite Earth Initiative (supported by Leslie and Mac McQuown)Northwestern University’s Finite Earth Initiative (supported by Leslie and Mac McQuown)Northwestern University’s Finite Earth Initiative (supported by Leslie and Mac McQuown)issue-copyright-statement© Springer Nature Limited 2021
==== Body
pmcIntroduction

In an alternating current power grid, the generators provide electrical power that oscillates in time as sinusoidal waves. As these waves are superimposed before reaching the consumers, they need to be synchronized to the same frequency; otherwise, time-dependent cancellation between these waves would cause the delivered power to fluctuate, which can lead to equipment malfunction and damage1. Maintaining frequency synchronization is challenging because the system is complex in various ways, with every generator responding differently to the continual influence of disturbances and varying conditions2. Adding to the challenge is the increase in perturbations resulting from the ongoing integration of energy from intermittent sources3, the emergence of grid-connected microgrids4, and the expansion of an increasingly open electricity market5. Furthermore, the inherent heterogeneities in the parameters of system components and in the structure of the interaction network are perceived as obstacles to achieving synchronization. Consistent with the view that heterogeneities may generally inhibit frequency homogeneity, an earlier study showed that homogenizing the (otherwise heterogeneous) values of generator parameters can lead to stronger stability of synchronous states than in the original system6. An outstanding question remains, however, as to whether there is a heterogeneous parameter assignment (different from the nominal one) that would enable even stronger stability for synchronous states than the best homogeneous parameter assignment. Though motivated by its significance for power grids, this question is broadly relevant for improving the stability of homogeneous dynamics in complex network systems in general, including consensus dynamics in networks of human or robotic agents7,8, coordinated spiking of neurons in the brain9,10, and synchronization in communication networks11,12.

To gain insights into the potential role of heterogeneity in enhancing stability, it is instructive to first consider the case of damped harmonic oscillators. For a single oscillator, the optimal stability corresponds to the fastest convergence to the stable equilibrium and is achieved when the oscillator is critically damped: underdamping would lead to lingering oscillations around the equilibrium, and overdamping would lead to slowed convergence due to excess dragging. This optimization is exploited in door closers (devices that passively close doors in a controlled manner), which are designed to be critically damped for the door to close fast without slamming. When multiple damped oscillators are coupled, the damping giving rise to optimal stability will be influenced by the network interactions. More important, we can show that the optimal stability in such a network requires different oscillators to have different damping (even when their other parameters are all identical and they are positioned identically in the network), as illustrated in Fig. 1.Fig. 1 Stabilizing effect of heterogeneity in a mass-spring system.

a System consisting of a linear chain of three unit masses connected by two identical springs. The masses are constrained to move horizontally, and their dynamics are governed by the equation shown, where xi is the displacement of mass i relative to its equilibrium and bi is its damping coefficient. b Total potential energy of the springs vs. time for three different damping scenarios. The optimal damping (red), corresponding to the fastest energy decay, is achieved for b1≈2.47, b2≈3.17, b3≈1.47 (or equivalently, b1≈1.47, b2≈3.17, b3≈2.47), despite the fact that masses 1 and 3 are otherwise identical and identically coupled. Overdamping leads to a slower monotonic decay, while underdamping results in a slower oscillatory decay, as shown in blue by varying b1 and b3 by a factor of 5. In all cases, the initial conditions are (x1, x2, x3) = (1, 0, −1) and (x°1,x°2,x°3)=(0,0,0).

In this paper, we first demonstrate that an analogous effect occurs in power-grid networks: heterogeneity in generator parameters can robustly enhance both the linear and the nonlinear stability of synchronous states in power grids from North America and Europe. Since these systems have heterogeneity in the network structure in addition to the tunable generator parameters, one possibility is that the effect arises entirely from compensation: stability reduction due to one heterogeneity is compensated by another heterogeneity, leading to a stability enhancement when the latter heterogeneity is added. An alternative, which we validate here, involves the recently established phenomenon of converse symmetry breaking (CSB)13, in which the stability of a symmetric state requires the system’s symmetry to be broken. Owing to its counterintuitive nature, this phenomenon had not been recognized until it was recently predicted and experimentally confirmed13,14 for synchronization in oscillator networks (a class of network dynamics widely studied in the literature15–17). Despite its conceptual generality and potential to underlie symmetric states of many systems, this phenomenon has not yet been observed outside laboratory settings. The symmetry relevant here is node-permutation symmetry, since in a synchronized state the states of different nodes are equal and can be permuted without altering the state of the system. For power grids, CSB would translate to a stability enhancement mechanism in which maintaining the stability of synchronous (and thus symmetric) states requires the generator parameters to be heterogeneous (thus making the system asymmetric). By systematically removing all the other system heterogeneities and isolating the effect of the generator heterogeneity, we establish that CSB is responsible for a significant portion of the stability improvement observed in the power grids we consider. This offers insights into mechanisms underlying the parameter heterogeneity that arises when the generators are tuned to damp oscillations18,19 (e.g., by adjusting devices called power system stabilizers). Our results are of particular relevance given that CSB has thus far not been observed in any real-world system outside controlled laboratory conditions, let alone power-grid networks.

Results

Power-grid dynamics and stability

To describe the dynamics of n generators in a power-grid network, we represent each generator node as a constant voltage source behind a reactance (the so-called classical model) and their interactions through intermediate non-generator nodes as effective impedances (a process known as Kron reduction)20. We assume that the system is operating near a synchronous state in which the voltage frequencies of the n generators are all equal to a constant reference frequency ωs, and we examine whether the homogeneous state is stable against dynamical perturbations. Such perturbations, whether they are small or large, may come for instance from sudden changes in generation and/or demand due to shifting weather condition at wind or solar farms, variations in power consumption, switching on/off connections to microgrids, etc. The short-term dynamics (of the order of one second or less) are then governed by the so-called swing equation20,21:1 δi¨+βiδi°=ai−∑k≠iciksinδi−δk−γik,

where δi is the phase angle variable for generator i (representing the generator’s internal electrical angle, relative to a reference frame rotating at the reference frequency ωs); βi≡Di/(2Hi) is an effective damping parameter (corresponding to bi in the mass-spring system of Fig. 1), with constant Di capturing both mechanical and electrical damping and constant Hi representing the generator’s inertia; ai is a parameter representing the net power driving the generator (i.e., the mechanical power provided to the generator, minus the power demanded by the network, including loss due to damping); and cik and γik are respectively the coupling strength and phase shift characterizing the electrical interactions between the generators. The parameters in Eq. (1) for a given system are determined by computing the active and reactive power flows between network nodes from system data and using them to calculate the complex-valued effective interaction (and thus its magnitude cik and angle γik) between every pair of generators. In real power grids, stable system operation is ensured by a hierarchy of controllers that adjust generator power outputs and thus the parameters in Eq. (1). Here, however, these parameters can be regarded as constants, since the lowest level of control (known as the primary control) is modeled as a damping-like effect captured by the βi term in Eq. (1), while the upper-level controls (known as the secondary and tertiary controls) act on time scales much longer than that of the short-term generator dynamics described by the model. In addition, fluctuations in power generation and demand on the time scales of minutes or longer (which can come, e.g., from renewable energy sources) do not affect the short-term dynamics. Equation (1) has recently been studied extensively in the network dynamics community3,6,22–25.

We first analyze the stability of the synchronous state against small perturbations. The synchronous state corresponds to a fixed point of Eq. (1) given by δi=δi* and δ°i=0, which represents frequency synchronization because δ°i is the frequency relative to the reference ωs. The Jacobian matrix of Eq. (1) at this point can be written as2 J=OI−P−B,

where O and I denote the n × n null and identity matrices, respectively; P = (Pik) is the n × n matrix defined by3 Pik=−cikcos(δi*−δk*−γik),i≠k,−∑k′≠iPik′,i=k,

which expresses the effective interactions between the generators; and B is the n × n diagonal matrix with βi as its diagonal elements. We note that, while the form of the Jacobian matrix for coupled damped harmonic oscillators is the same as in Eq. (2), power grids are different in that they can have P ≠ PT because cik ≠ cki in general and because γik appears in Eq. (3). The stability under noiseless conditions is determined by the Lyapunov exponent defined as λmax≡maxi≥2Re(λi), where λi are the eigenvalues of J. The identically zero eigenvalue, which comes from the zero row-sum property of P and is denoted here by λ1, is excluded because it is associated with the invariance of the equation under uniform shift of phases. If λmax<0, then the synchronous state is asymptotically stable, and smaller λmax implies stronger stability (this is known as small-signal stability analysis in power system engineering). Since real power-grid dynamics are noisy due to power generation/demand fluctuations and various other disturbances occurring on short time scales, λmax needs to be sufficiently negative to keep the system close to the synchronous state. Indeed, a previous study14 showed that, for broad classes of noise dynamics, there is a (negative) threshold value of λmax for such stability: the system is stable if and only if λmax is below the threshold. This stability threshold depends on the noise intensity level. For impulse-like disturbances, the intensity level corresponds to the maximum deviation of δi that can be induced by a single disturbance, such as a sudden loss of a generator or a spike in power demand. For continual disturbances, the intensity level can be quantified by the variances of the fluctuating power generation and demand, which can be modeled by adding a randomly varying term to the parameter ai. Since the stability threshold is generally lower for higher noise levels, the lower the value of λmax for a given power grid, the more intense disturbances and fluctuations the system can endure without losing stability. Incidentally, the optimal damping in the mass-spring system of Fig. 1 is given precisely by minimizing λmax for that system.

Enhancing stability with generator heterogeneity

We now study λmax=λmax(β) as a function of β≡(β1,…,βn) for a selection of power grids whose dynamics can be described by Eq. (1) with the parameter values based on data. Using the same model, it was previously shown6 that, under the constraint that all βi’s have the same value, λmax is minimized when β = β=, where β=≡(β=,…,β=) and β=≡2α2, with α2 denoting the smallest nonidentically zero eigenvalue of matrix P. The eigenvalue α2 is associated with the least stable eigenmode, and we assume that it is real and positive (as confirmed in all systems we consider). It was further shown that, at this homogeneous optimal point β=, the function λmax(β) is non-differentiable (which precludes the use of a standard derivative test), but its one-sided derivative along any given straight-line direction is positive, i.e., the directional derivative D𝛃′λmax(β=) is positive in the direction of any n-dimensional vector β′. Thus, moving away from β= along any straight line would necessarily increase λmax from the local minimum value λmax(β=)=−α2 and hence only reduce the stability of the synchronous state.

Despite the apparent impossibility of improving on λmax(β=) locally, we first show that there can be curved paths starting at β= along which λmax can be further minimized with heterogeneous βi. Indeed, Fig. 2a illustrates using a 3-generator system that such curved paths exist and can connect β= to the (unique) global minimum, which we denote by β≠ as its components are all different. The corresponding optimal λmax(β≠)≈−9.41 represents more than 8% improvement over λmax(β=)≈−8.69. In general, if a curved path starts at β=, and if λmax decreases monotonically along that path, then it cannot be oriented in an arbitrary direction in the β-space. We show that it needs to be tangent to a system-specific plane (or hyperplane of co-dimension one for n > 3), denoted here by L and defined by the equation ∑i=1nuiviβi=0, where ui and vi are the ith component of the left and right eigenvectors, respectively, associated with the eigenvalue α2. This result, illustrated by the three example paths in Fig. 2a, follows from the derivation of a formula for λmax and the full analytical characterization of the stability landscape near β= (both presented in Supplementary Note 1).Fig. 2 Enhancing the stability of the 3-generator system along curved paths.

a Examples of curved paths (red, green, and blue) in the β-space from β= to β≠ along which the Lyapunov exponent λmax decreases monotonically. The numerically generated curves confirm our theoretical prediction that these paths are all tangent to the plane L at the point β=. The droplines indicate that these curves are outside L. The contour levels of λmax are shown on L. Also shown is the plane M, which is defined as the plane perpendicular to L that contains β= and β≠. b Contour levels of λmax on the plane M, which contains the green path from a. The orange lines in a and b indicate the intersection between planes L and M. The red dashed curves (and the green path) trace cusp surfaces associated with degeneracies of the real parts of the eigenvalues of J that determine λmax. Three types of degeneracy are indicated: a real eigenvalue equal to the real parts of a pair of complex conjugate eigenvalues (a, a ± bi), two different complex conjugate eigenvalues with equal real parts (a ± b1i, a ± b2i), and two equal real eigenvalues (a, a). For details on the system, see “Methods”.

The curved paths of decreasing λmax are part of the complex structure of the stability landscape. These paths generally lie on a cusp surface, defined by the property that, at any point on the surface, λmax is non-differentiable and locally minimum along any direction transverse to the surface. The three paths shown in Fig. 2a all lie on the same cusp surface, which contains both β= and β≠. The intersection between this cusp surface and the plane M (the one perpendicular to L) is the green path of monotonically decreasing λmax shown in Fig. 2. In fact, there are infinitely many different paths of decreasing λmax on this cusp surface. Of these paths, the red and blue paths shown in Fig. 2a share the additional property of being an intersection between pairs of cusp surfaces. Each of these paths consists of at least two parts that are intersections between different pairs of cusp surfaces, which explains the kinks observed in Fig. 2 as points at which the curve switches from one intersecting surface to another. In larger systems, we find that their higher-dimensional β-spaces are sectioned by many entangled cusp hypersurfaces associated with spectral degeneracies (as illustrated in Supplementary Fig. 2 using the four larger systems we will introduce below). Their intersections, which themselves form cusp hypersurfaces of lower dimensions, are expected to contain curved paths of monotonically decreasing λmax. The existence of kinks and cusp surfaces in the stability landscape, which makes numerical search for global optima challenging, is not unique to power grids nor phase oscillator networks. It is a consequence of a much more general mathematical observation that the largest real part of the eigenvalues of a matrix (known as the spectral abscissa), such as λmax we consider here, is a non-smooth, non-convex, and non-Lipschitz function of the matrix elements26.

Stabilizing heterogeneity in real power grids

Having established that heterogeneous β≠ can improve stability over the homogeneous β= for a small example system, we now show that this result extends to much larger, real-world power grids. Specifically, we study the 48-generator NPCC portion of the North American power grid and the 69-generator German portion of the European power grid. Assessing the stability against small perturbations based on Eqs. (1)–(3) has the advantage of reducing the complexity of these systems to a single matrix P and its eigenvalues. For each system, we identify a local minimum β≠ that has heterogeneous βi (and thus is distinct from β=) and achieves the lowest λmax over 200 independent runs of simulated annealing. We find simulated annealing to be more effective than other methods in locating a minimum on a non-differentiable landscape27. The resulting stability improvement over β= is substantial: λmax(β=)−λmax(β≠)=0.66 for the NPCC network and λmax(β=)−λmax(β≠)=0.42 for the German network. The optimized βi assignment in β≠ exhibits substantial heterogeneity across each network and also across the corresponding geographical area (Fig. 3).Fig. 3 Heterogeneity of optimized generator parameters βi for two power grids.

a Portion of the North American power grid corresponding to the former Northeast Power Coordinating Council (NPCC) region. b German portion of the European power grid. In both panels, a color-coded circle represents a generator or an aggregate of generators (see “Methods” for the aggregation procedure used), with the color indicating the corresponding optimal βi in the vector β≠. The arrows above the color bars indicate β=, the optimal uniform value of βi. The λmax values for the uniform and non-uniform optimal βi (in the vectors β= and β≠, respectively) are indicated at the bottom of each panel. The radius of the circle is proportional to the real power output of the generator in megawatts (MW). Small green dots indicate non-generator nodes. Details on these systems, including data sources, can be found in “Methods”.

To validate the prevalence of such stability-enhancing heterogeneity, we also analyze λmax as a function of system stress level (to be precisely defined below) for four different systems, including the two used in Fig. 3 (see “Methods” for detailed descriptions of the systems and data sources). To increase or decrease the level of stress in a given system, we scale the power output of all generators and the power demand at all nodes by a common constant factor. We then re-compute the power flows across the entire network and the parameters of Eq. (1). The system stress level is then defined to be the common scaling factor used in this procedure. Thus, a stress level of 1 for a given system corresponds to the original demand level in the corresponding dataset. For each stress level, we estimate λmax at β = β≠ from 200 independent simulated annealing runs. Over the entire range of stress levels considered, we consistently observe a smaller λmax for β≠ compared to β=, the optimal homogeneous βi assignment, and to β0, the original βi assignment in the dataset (Fig. 4a).Fig. 4 Improving the stability of power grids with heterogeneity in β.

The columns correspond to the four systems we consider. a Improved Lyapunov exponent λmax as functions of the system stress level for the heterogeneous optimum β≠ (blue), the homogeneous optimum β= (red), and the original parameter β0 (black). The cases shown in Fig. 3 are indicated in the second and the last plot. b Change of λmax under perturbations of size ε applied to β≠ (blue). We show λmax as a function of ε, where solid and dashed curves indicate the average and the maximum, respectively, over perturbations in 1000 random directions. Note that the maximum corresponds to the worst case scenario. For comparison, we also show the average of λmax when β= is perturbed (red). c Fraction f of trajectories that converge to synchronous states before a given cutoff time tc for β≠ (blue), β= (red), and β0 (black). Note that f for β0 remains zero for all tc < 10 seconds in all cases.

To test the robustness of the identified optimal λmax against uncertainties in the βi values, we study how λmax changes under perturbations along random directions in the β-space in the vicinity of β≠ and (for comparison) in the vicinity of β=. For the stress level of 1 and for each random direction, we compute λmax as a function of the perturbation size ε, measured in 2-norm. The resulting statistics from 1000 random directions indicate that, for each system, there is a sizable neighborhood of the optimum β≠ in which λmax is significantly lower than at β=, representing a stability improvement against small perturbations (Fig. 4b).

To show that the improvement is also observed for stability against large perturbations, we define a generalized notion of attraction basin as a set of initial conditions whose corresponding trajectories satisfy a criterion for convergence to synchronous states (a variation of the so-called basin stability28). Here, the convergence criterion we use is that the instantaneous frequency enters into a narrow band around ωs (within ± 0.3 Hz) and remains inside the band until tmax=10 seconds. This criterion is similar to what is typically used for transient stability analysis in power system engineering. It also captures a variety of synchronous states, including not only those corresponding to fixed points of Eq. (1) (with constant phase angle differences), but also those corresponding to time-dependent solutions of Eq. (1). To account for large perturbations, we consider initial conditions with arbitrary phase angles and frequencies within 1 Hz of the nominal frequency (60 Hz for the New England and NPCC systems; 50 Hz for the U.K. and German systems). Each initial condition can be regarded as resulting from a large impulse-like disturbance, such as a disconnection of a significant portion of the grid or a system-wide demand surge. The size of the basin can then be quantified using the fraction f of the corresponding trajectories that converge before a given cutoff time tc, i.e., the fraction of those that satisfy ∣δ°i(t)∣/(2π)≤0.3 Hz for all t∈[tc,tmax]. For each tc, the fraction f is estimated using 1000 initial conditions sampled randomly and uniformly from all states satisfying the criteria described above. As shown in Fig. 4c, we find that the estimated f is significantly larger for β≠ than for β= (which in turn is much larger than for β0). This indicates that the likelihood for the system to return to stable operation after a large disturbance is higher for the heterogeneous optimal βi than for the homogeneous optimal ones. We also observe that larger systems tend to exhibit larger increase in the size of the asymptotic basins (i.e., in the value of f for tc → ∞).

Isolating converse symmetry breaking

Since real power systems generally have heterogeneity in ai, cik, and γik, the stability improvement enabled by the βi heterogeneity (and the associated system asymmetry) could in principle be a compensation for heterogeneity in the network structure, power demand and generation, or other component parameters (and the associated system asymmetries). To illustrate that no such compensation is needed and that CSB can be responsible for stability improvement, we use an example system consisting of four generators connected to each other and to one load (see Supplementary Fig. 1 for a system diagram). This system is symmetric with respect to the permutation of generators 2 and 3 if β2 = β3, and this symmetry is reflected in the property that P2j = P3j for all j in the corresponding interaction matrix P (Fig. 5a). The minimum λmax possible for this symmetric system is λmax≈−2.40, which can be decreased further by >20% to λmax≈−2.97 if the β2 = β3 constraint is lifted (Fig. 5b–d). This demonstrates CSB for this system under a range of noise levels: breaking the system’s symmetry under the permutation of generators 2 and 3 is required for λmax to cross the stability threshold and make the (symmetric) synchronous state stable. We note that the observation of CSB depends on the system’s symmetry. While CSB is observed in this 4-generator system (with a two-generator permutation symmetry), we do not observe CSB in a variant of the system with the four-generator permutation symmetry. We also note that, while the optimal βi assignment does not share the two-generator permutation of the system, the two-dimensional stability landscape does, and it features a pair of equally optimal assignments related to each other through the symmetry (Fig. 5b). It is instructive to compare this result with the mass-spring system in Fig. 1, where similar breaking of a permutation symmetry (between masses 1 and 3) for a symmetric landscape (where optimal b1 and b3 are necessarily different but can be swapped) is shown to underlie optimal damping.Fig. 5 Illustrating CSB in a 4-generator example system.

a Effective interaction network connecting the generators and given by matrix P. Both the thickness and color of the arrow connecting node j to node i represent the interaction strength ∣Pij∣, with the thickness proportional to ∣Pij∣ and the color encoding ∣Pij∣ normalized by its maximum over all i and j with i ≠ j. b–d Dependence of λmax on β2 and β3, with the values of β1 and β4 set to the values indicated in a. In b, λmax is color-coded to visualize the full 2D landscape. In c, λmax is shown as a function of β2 along the white dashed line in b, corresponding to β2 = β3. It attains its minimum value ≈ −2.40 at β2 = β3 ≈ 4.80 (red circle). In d, λmax as a function of β2 along the black dashed line in b attains its minimum value ≈ −2.97 at (β2, β3) ≈ (4.27, 5.17) (white circle). Thus, the substantially improved optimal λmax is possible only when the permutation symmetry between generators 2 and 3 is broken. For details on the system, see “Methods”.

Having established that βi heterogeneity alone can enhance stability through CSB, we now introduce a systematic method to separate CSB from other mechanisms that involve interplays between multiple heterogeneities. For this purpose, we transform matrix P for each system used in Fig. 4, which does not have a pairwise node permutation symmetry, to a slightly different matrix P′ that does have the symmetry. More precisely, for a given pair of nodes i1 and i2, we define this symmetrized matrix P′ by Pi1j′=Pi2j′≡(Pi1j+Pi2j)/2 for all j ≠ i1 nor i2, making it symmetric under the permutation of nodes i1 and i2. To elucidate CSB for each system, we choose a node pair that simultaneously minimizes the difference between P and P′ and maximizes the amount of stability improvement observed for P′. For the four systems in Fig. 4, the stability of the symmetrized system can clearly be enhanced by allowing βi1≠βi2, and much of the enhancement is maintained as one interpolates from the symmetrized system back to the original system (Fig. 6). This indicates that a significant portion of the stability enhancement for the original system can be attributed to CSB.Fig. 6 Isolating the CSB effect in power grids.

Each column synthesizes results for the system indicated at the top. a 2D stability landscape λmax(βi1,βi2), where i1 and i2 are the nodes whose permutation holds the symmetrized matrix P′ invariant. In each panel, the red circle marks the optimal (βi1,βi2) on the diagonal βi1=βi2 (white dashed line), while the white circle marks the optimal when βi1≠βi2 is allowed. The other βi are fixed at the values in β≠ identified for a stress level of 1 in Fig. 4. b Stability λmax as a function of βi1 along the black dashed line connecting the red and white circles in a, respectively. c Input strength patterns of the nodes i1 and i2 for the original matrix P. The color of each arrow indicates ∣Pij∣, normalized by the largest value of ∣Pij∣ shown. For nodes i1 and i2, we show incoming links from the top six common neighbors in terms of the input strength. Also shown is the distance d from the original matrix P to its symmetrized version P′, in which the two nodes receive identical incoming (weighted) links, defined as d≡∥P−P′∥2/∥P∥2, where ∥ ⋅ ∥2 denotes the matrix 2-norm. d Change in the optimal λmax with (red) and without (blue) the constraint βi1=βi2, as we interpolate between the original matrix P and its symmetrized version P′. Node indexing for all four systems are described in “Methods”.

Discussion

Our demonstration that heterogeneity of generators can enhance the stability of synchronous states in a range of power grids suggests that there is large previously under-explored potential for tuning and upgrading current systems for better stability. Since larger conventional generators have larger inertia and thus larger impact on the stability of other generators, tuning of their parameters may be particularly beneficial. While we focused on the heterogeneity of a specific generator parameter here, further stability enhancement is likely to be possible by exploiting heterogeneity in other generator parameters and in the parameters of other network components as well as in the network topology. We suggest that such stability enhancement opportunities exist beyond power systems and may extend to any network whose function benefits from homogeneous dynamics and whose stability depends on tunable system parameters. For example, the results presented here suggest that CSB can potentially be observed for coupled oscillatory flows in microfluidic networks and for networks of coupled chemical reactors whose oscillatory node dynamics is close to a Hopf bifurcation. It is known13 that such systems can be parameterized so that their Jacobian matrices take a form that generalizes Eq. (2) and thus is conducive to the emergence of CSB. The approach we developed here to isolate CSB is versatile and can be applied broadly to systems for which different heterogeneities co-occur. Determining how prevalent CSB is and how it depends on the properties of the system (e.g., the network size, link distribution, and node dynamics) are important questions for future research.

It is instructive to interpret our results and contrast them with past approaches in network optimization. In seeking the best approach, one may form two complementary hypotheses. One hypothesis, invoked in the past, was that the stability of the desired homogeneous states would be optimal when the system is homogeneous; the approach would thus be to limit the optimization search to the low-dimensional parameter subspace corresponding to networks with identical parameter values for all nodes. The other hypothesis, validated here, is that optimal stability of the desired homogeneous states is generally obtained with heterogeneous parameter assignment, which implies that the search for this optimum requires exploring the high-dimensional parameter space without making a priori assumptions on how the parameters of different nodes are related. Recognizing this can lead to new control approaches designed to manipulate these parameters for further optimization of stability. We suggest that the fresh opportunities for network optimization and control revealed in this study apply to network systems in general and thus have the potential to inspire new discoveries in many different disciplines.

Methods

Power-grid datasets

Here, we describe the sources of data for the six power-grid networks considered (the 3-generator system in Fig. 2; the New England, NPCC, U.K., and German systems in Figs. 3, 4, and 6; and the 4-generator system in Fig. 5). For each system, the data provide the net injected real power at all generator nodes, the power demand at all non-generator nodes, and the parameters of all power lines and transformers. These parameters are sufficient to determine all active and reactive power flows in the system using a standard power flow calculation. The data also provide the generators’ dynamic parameters Hi, Di, and xint,i used in our stability calculations. The parameters Hi and Di are the inertia and damping constants, respectively, that define the effective damping parameter through the relation βi = Di/(2Hi). The parameter xint,i represents the internal reactance of generator i and is used in the calculation of the parameters ai, cij, and γij. In each system, nodes are indexed as in the original data source (except for the German power grid; see below).

3-generator test system (3-gen). For this IEEE 3-generator, 9-node test system, which appeared in ref. 20, we used the data file (data3m9b.m) available in the PST toolbox29. This system represents the Western System Coordinating Council (WSCC), which was part of the region now called the Western Electricity Coordinating Council (WECC) in the North American power grid. The data file provides all necessary dynamical parameters for each generator.

New England test system (10-gen). For the IEEE 10-generator, 39-node test system, as described in refs. 30 and 31, we used the data file (case39.m) available in the MATPOWER toolbox32, with dynamic parameters added manually from ref. 30. This is a reduced model representing the New England portion of the Eastern Interconnection in the North American power grid, with one generator representing the connection to the rest of the grid.

NPCC power grid (48-gen). For the 48-generator, 140-node NPCC power grid33, we used the data file (data48em.m) available in the PST toolbox29. The system represents the former NPCC region of the Eastern Interconnection in the North American power grid and includes an equivalent generator/load node representing the rest of the Interconnection. The data file provides Hi and xint,i for all generators (while it assumes Di = 0). We generated Di randomly by sampling from the uniform distribution on the interval [1, 3] (in per unit on the system base, as specified by the data file). The geographic coordinates of the nodes used in Fig. 3a were extracted from ref. 34, and the coastline and boundary data used to draw the map were obtained from Natural Earth35.

U.K. power grid (66-gen). For the 66-generator, 29-node U.K. power grid, we used the data file (GBreducednetwork.m) available from ref. 36. The system represents a reduced model for the power grid of Great Britain. The dynamical parameters, Hi, Di, and xint,i, were generated randomly by sampling from the uniform distribution on the intervals, [1, 5], [1, 3], and [0.001, 0.101], respectively. The generated parameters values for each generator are in per unit on its own machine base, i.e., normalized by the reference values computed from the power base for the generator (chosen to be 1.5 times the maximum real power generation provided in the data file). For stability calculations, we converted these values to the corresponding values in per unit on a common system base.

German power grid (69-gen). For the 69-generator, 228-node German power grid, we created the data from the ENTSO-E 2009 Winter model37. The ENTSO-E model is a DC power flow model of the continental Europe and contains 1,486 nodes and 565 generators. We first created a dynamical model for the entire ENTSO-E network by solving the DC power flow and converting it to an AC power flow solution (assuming a 0.95 power factor at each node), and then generating dynamical parameters using the same method as for the U.K. grid. For any node with multiple generators attached, the net reactive power injection was distributed among these generators in proportion to their real power generation. From this full ENTSO-E model, we extracted the German portion by eliminating (using Kron reduction) all the nodes outside Germany (identified using the country label “D” representing Germany in the dataset). We re-indexed the extracted nodes consecutively, preserving the original ordering. The geographic coordinates of the nodes used in Fig. 3b were extracted from the PowerWorld data files available from ref. 37, and the coastline and boundary data used to draw the map were obtained from Natural Earth35.

4-generator example system. For the 4-generator, 5-node example system used in Fig. 5, we show a full system diagram in Supplementary Fig. 1, indicating the main parameters of the components. When the damping parameters of generators 2 and 3 are equal (i.e., β2 = β3), the system is symmetric under the permutation of these generators. MATLAB code for running simulations on this system, which includes the full set of parameters and uses the MATPOWER toolbox32, is available from our GitHub repository38.

Aggregation of generators and effective damping parameter βi

If a subset of generators are synchronized in the sense that δi − δj is constant in time for any two generators i and j in the subset, then they can be represented by a single equivalent generator using a Zhukov-based aggregation method similar to that described in ref. 33. In this method, the equivalent generator has inertia constant ∑iHi and damping constant ∑iDi, where the sums are taken over the generators i in the subset. The effective damping parameter of the equivalent generator is then ∑iDi/(2∑iHi)=D¯/(2H¯), where D¯ and H¯ are respectively the average of the inertia and damping constants of the generators in the subset. Thus, the aggregation does not introduce any artifactual heterogeneity.

Supplementary information

Supplementary Information

Supplementary information

The online version contains supplementary material available at 10.1038/s41467-021-21290-5.

Acknowledgements

The authors thank Alex Mercanti and Yuanzhao Zhang for insightful discussions. This research was supported by Northwestern University’s Finite Earth Initiative (supported by Leslie and Mac McQuown) and ARPA-E Award No. DE-AR0000702 and benefited from logistical support provided by Northwestern University’s Institute for Sustainability and Energy.

Author contributions

F.M., T.N., and A.E.M. designed the research, analyzed the results, and wrote the paper. F.M. performed the simulations. All authors approved the final manuscript.

Data availability

Data on all six systems we consider (described in “Methods”) and detailed data of the core results presented in the figures are available from our GitHub repository38.

Code availability

Essential code for reproducing the core results in all figures, as well as scripts for generating plain versions of the figures, is available from the GitHub repository38.

Competing interests

The authors declare no competing interests.

Peer review information Nature Communications thanks Benjamin Schäfer and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Change history

9/18/2024

A Correction to this paper has been published: 10.1038/s41467-024-52495-z
==== Refs
References

1. Machowski J Lubosny Z Bialek JW Bumby JR Power System Dynamics: Stability and Control 2020 Hoboken John Wiley
Machowski, J., Lubosny, Z., Bialek, J. W. & Bumby, J. R. Power System Dynamics: Stability and Control (John Wiley, Hoboken, 2020).
2. Backhaus S Chertkov M Getting a grip on the electrical grid Phys. Today 2013 66 42 48 10.1063/PT.3.1979
Backhaus, S. & Chertkov, M. Getting a grip on the electrical grid. Phys. Today 66, 42–48 (2013).
3. Schäfer B Beck C Aihara K Witthaut D Timme M Non-Gaussian power grid frequency fluctuations characterized by Lévy-stable laws and superstatistics Nat. Energy 2018 3 119 126 10.1038/s41560-017-0058-z
Schäfer, B., Beck, C., Aihara, K., Witthaut, D. & Timme, M. Non-Gaussian power grid frequency fluctuations characterized by Lévy-stable laws and superstatistics. Nat. Energy 3, 119–126 (2018).
4. Olivares DE Trends in microgrid control IEEE T. Smart Grid 2014 5 1905 1919 10.1109/TSG.2013.2295514
Olivares, D. E. et al. Trends in microgrid control. IEEE T. Smart Grid 5, 1905–1919 (2014).
5. Griffin, J. M. & Puller, S. L. eds. Electricity Deregulation: Choices and Challenges (University of Chicago Press, Chicago, 2009).
6. Motter AE Myers SA Anghel M Nishikawa T Spontaneous synchrony in power-grid networks Nat. Phys. 2013 9 191 197 10.1038/nphys2535
Motter, A. E., Myers, S. A., Anghel, M. & Nishikawa, T. Spontaneous synchrony in power-grid networks. Nat. Phys. 9, 191–197 (2013).
7. Judd S Kearns M Vorobeychik Y Behavioral dynamics and influence in networked coloring and consensus Proc. Natl Acad. Sci. USA 2010 107 14978 14982 10.1073/pnas.1001280107 20696936
Judd, S., Kearns, M. & Vorobeychik, Y. Behavioral dynamics and influence in networked coloring and consensus. Proc. Natl Acad. Sci. USA 107, 14978–14982 (2010).20696936
8. Ren W Cao Y Distributed Coordination of Multi-agent Networks 2010 London Springer
Ren, W. & Cao, Y. Distributed Coordination of Multi-agent Networks (Springer, London, 2010).
9. Axmacher N Mormann F Fernández G Elger CE Fell J Memory formation by neuronal synchronization Brain Res. Rev. 2006 52 170 182 10.1016/j.brainresrev.2006.01.007 16545463
Axmacher, N., Mormann, F., Fernández, G., Elger, C. E. & Fell, J. Memory formation by neuronal synchronization. Brain Res. Rev. 52, 170–182 (2006).16545463
10. Penn Y Segal M Moses E Network synchronization in hippocampal neurons Proc. Natl Acad. Sci. USA 2016 113 3341 3346 10.1073/pnas.1515105113 26961000
Penn, Y., Segal, M. & Moses, E. Network synchronization in hippocampal neurons. Proc. Natl Acad. Sci. USA 113, 3341–3346 (2016).26961000
11. Bregni S Synchronization of Digital Telecommunications Networks 2002 West Sussex John Wiley & Sons
Bregni, S. Synchronization of Digital Telecommunications Networks (John Wiley & Sons, West Sussex, 2002).
12. Wu Y-C Chaudhari Q Serpedin E Clock synchronization of wireless sensor networks IEEE Signal. Proc. Mag. 2010 28 124 138 10.1109/MSP.2010.938757
Wu, Y.-C., Chaudhari, Q. & Serpedin, E. Clock synchronization of wireless sensor networks. IEEE Signal. Proc. Mag. 28, 124–138 (2010).
13. Nishikawa T Motter AE Symmetric states requiring system asymmetry Phys. Rev. Lett. 2016 117 114101 10.1103/PhysRevLett.117.114101 27661690
Nishikawa, T. & Motter, A. E. Symmetric states requiring system asymmetry. Phys. Rev. Lett. 117, 114101 (2016).27661690
14. Molnar F Nishikawa T Motter AE Network experiment demonstrates converse symmetry breaking Nat. Phys. 2020 16 351 356 10.1038/s41567-019-0742-y
Molnar, F., Nishikawa, T. & Motter, A. E. Network experiment demonstrates converse symmetry breaking. Nat. Phys. 16, 351–356 (2020).
15. Arenas A Díaz-Guilera A Kurths J Moreno Y Zhou C Synchronization in complex networks Phys. Rep. 2008 469 93 153 10.1016/j.physrep.2008.09.002
Arenas, A., Díaz-Guilera, A., Kurths, J., Moreno, Y. & Zhou, C. Synchronization in complex networks. Phys. Rep. 469, 93–153 (2008).
16. Kiss IZ Rusin CG Kori H Hudson JL Engineering complex dynamical structures: sequential patterns and desynchronization Science 2007 316 1886 1889 10.1126/science.1140858 17525302
Kiss, I. Z., Rusin, C. G., Kori, H. & Hudson, J. L. Engineering complex dynamical structures: sequential patterns and desynchronization. Science 316, 1886–1889 (2007).17525302
17. Matheny MH Exotic states in a simple network of nanoelectromechanical oscillators Science 2019 363 eaav7932 10.1126/science.aav7932 30846570
Matheny, M. H. et al. Exotic states in a simple network of nanoelectromechanical oscillators. Science 363, eaav7932 (2019).30846570
18. Rogers G Power System Oscillations 2012 New York Springer Science & Business Media
Rogers, G. Power System Oscillations (Springer Science & Business Media, New York, 2012).
19. Obaid ZA Cipcigan LM Muhssin MT Power system oscillations and control: classifications and PSSs’ design methods: a review Renew. Sust. Energ. Rev. 2017 79 839 849 10.1016/j.rser.2017.05.103
Obaid, Z. A., Cipcigan, L. M. & Muhssin, M. T. Power system oscillations and control: classifications and PSSs’ design methods: a review. Renew. Sust. Energ. Rev. 79, 839–849 (2017).
20. Anderson PM Fouad AA Power System Control and Stability 2003 Piscataway, NJ IEEE Press
Anderson, P. M. & Fouad, A. A. Power System Control and Stability (IEEE Press, Piscataway, NJ, 2003).
21. Nishikawa T Motter AE Comparative analysis of existing models for power-grid synchronization New J. Phys. 2015 17 015012 10.1088/1367-2630/17/1/015012
Nishikawa, T. & Motter, A. E. Comparative analysis of existing models for power-grid synchronization. New J. Phys. 17, 015012 (2015).
22. Susuki Y Mezić I Hikihara T Coherent swing instability of power grids J. Nonlinear Sci. 2011 21 403 439 10.1007/s00332-010-9087-5
Susuki, Y., Mezić, I. & Hikihara, T. Coherent swing instability of power grids. J. Nonlinear Sci. 21, 403–439 (2011).
23. Menck PJ Heitzig J Kurths J Schellnhuber HJ How dead ends undermine power grid stability Nat. Commun. 2014 5 3969 10.1038/ncomms4969 24910217
Menck, P. J., Heitzig, J., Kurths, J. & Schellnhuber, H. J. How dead ends undermine power grid stability. Nat. Commun. 5, 3969 (2014).24910217
24. Yang Y Motter AE Cascading failures as continuous phase-space transitions Phys. Rev. Lett. 2017 119 248302 10.1103/PhysRevLett.119.248302 29286707
Yang, Y. & Motter, A. E. Cascading failures as continuous phase-space transitions. Phys. Rev. Lett. 119, 248302 (2017).29286707
25. Schäfer B Witthaut D Timme M Latora V Dynamically induced cascading failures in power grids Nat. Commun. 2018 9 1 13 29317637
Schäfer, B., Witthaut, D., Timme, M. & Latora, V. Dynamically induced cascading failures in power grids. Nat. Commun. 9, 1–13 (2018).29317637
26. Burke JV Overton ML Variational analysis of non-Lipschitz spectral functions Math. Programming 2001 90 317 351
Burke, J. V. & Overton, M. L. Variational analysis of non-Lipschitz spectral functions. Math. Programming 90, 317–351 (2001).
27. Nishikawa T Molnar F Motter AE Stability landscape of power-grid synchronization IFAC-PapersOnLine 2015 48 1 6 10.1016/j.ifacol.2015.11.001
Nishikawa, T., Molnar, F. & Motter, A. E. Stability landscape of power-grid synchronization. IFAC-PapersOnLine 48, 1–6 (2015).
28. Menck PJ Heitzig J Marwan N Kurths J How basin stability complements the linear-stability paradigm Nat. Phys. 2013 9 89 92 10.1038/nphys2516
Menck, P. J., Heitzig, J., Marwan, N. & Kurths, J. How basin stability complements the linear-stability paradigm. Nat. Phys. 9, 89–92 (2013).
29. Chow JH Cheung KW A toolbox for power system dynamics and control engineering education and research IEEE Trans. Power Syst. 1992 7 1559 1564 10.1109/59.207380
Chow, J. H. & Cheung, K. W. A toolbox for power system dynamics and control engineering education and research. IEEE Trans. Power Syst. 7, 1559–1564 (1992).
30. Pai M Energy Function Analysis for Power System Stability 1989 Norwell Kluwer Academic Publishers
Pai, M. Energy Function Analysis for Power System Stability (Kluwer Academic Publishers, Norwell, 1989).
31. Athay T Podmore R Virmani S A practical method for the direct analysis of transient stability IEEE Trans. Power Appar. Syst. 1979 PAS-98 573 584 10.1109/TPAS.1979.319407
Athay, T., Podmore, R. & Virmani, S. A practical method for the direct analysis of transient stability. IEEE Trans. Power Appar. Syst. PAS-98, 573–584 (1979).
32. Zimmerman RD Murillo-Sánchez CE Thomas RJ MATPOWER: steady-state operations, planning and analysis tools for power systems research and education IEEE Trans. Power Syst 2011 26 12 19 10.1109/TPWRS.2010.2051168
Zimmerman, R. D., Murillo-Sánchez, C. E. & Thomas, R. J. MATPOWER: steady-state operations, planning and analysis tools for power systems research and education. IEEE Trans. Power Syst. 26, 12–19 (2011).
33. Chow JH Power System Coherency and Model Reduction 2013 New York Springer
Chow, J. H. Power System Coherency and Model Reduction (Springer, New York, 2013).
34. Qi, J., Sun, K. & Kang, W. Optimal PMU placement for power system dynamic state estimation by using empirical observability Gramian. 2015 IEEE Power & Energy Society General Meeting (Denver, CO, USA, 2015).
35. Patterson, T. & Kelso, N. V. Natural Earth: Free Vector and Raster Map Data, March 2018, http://www.naturalearthdata.com/ (2018).
36. Bukhsh, W. A. & McKinnon, K. Network Data of Real Transmission Networks, April 2013, https://www.maths.ed.ac.uk/optenergy/NetworkData/reducedGB (2013).
37. Hutcheon, N. & Bialek, J. W. Updated and validated power flow model of the main continental European transmission network. 2013 IEEE Grenoble Conference, Grenoble, France, p. 1–5. Data files available at http://www.powerworld.com/bialek (2013).
38. Molnar, F., Nishikawa, T. & Motter, A. E. Asymmetry underlies stability in power grids (this paper), GitHub repository: code and data for analyzing converse symmetry breaking in power-grid networks, 10.5281/zenodo.4437866 (2021).
