Asymmetry underlies stability in power grids

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.


Introduction
In an alternating current power grid, the generators provide electrical power that oscillates in time as sinusoidal waves.Because 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 damage 1 .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 conditions 2 .Adding to the challenge is the increase in perturbations resulting from the ongoing integration of energy from intermittent sources 3 , the emergence of grid-connected microgrids 4 , and the expansion of an increasingly open electricity market 5 .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 system 6 .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 agents 7,8 , coordinated spiking of neurons in the brain 9,10 , and synchronization in communication networks 11,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 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 x i is the displacement of mass i relative to its equilibrium and b i 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 b 1 ≈ 2.47, b 2 ≈ 3.17, b 3 ≈ 1.47 (or equivalently, b 1 ≈ 1.47, b 2 ≈ 3.17, b 3 ≈ 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 b 1 and b 3 by a factor of 5.In all cases, the initial conditions are (x 1 , x 2 , x 3 ) = (1, 0, −1) and ( ẋ1 , ẋ2 , ẋ3 ) = (0, 0, 0).illustrated in Fig. 1.
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.Due to its counterintuitive nature, this phenomenon had not been recognized until it was recently predicted and experimentally confirmed 13,14 for synchronization in oscillator networks (a class of network dynamics widely studied in the literature [15][16][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 syn-chronized 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 oscillations 18,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 powergrid 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 equation 20,21 : 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 ); is an effective damping parameter (corresponding to b i in the mass-spring system of Fig. 1), with constant D i capturing both mechanical and electrical damping and constant H i representing the generator's inertia; a i 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 net-work, including loss due to damping); and c ik 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 c ik 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 community 3,6,[22][23][24][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 as where O and I denote the n × n null and identity matrices, respectively; P = (P ik ) is the n × n matrix defined by 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 ̸ = P T because c ik ̸ = c ki in general and because γ ik appears in Eq. (3).The stability under noiseless conditions is determined by the Lyapunov exponent defined as λ max ≡ max i≥2 Re(λ 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 study 14 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 a i .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 shown 6 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 nondifferentiable (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 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 n i=1 u i v i β i = 0, where u i and v i 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).
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 elements 26 .
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 48generator 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 landscape 27 .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).
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 recompute 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).
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 1,000 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 stability 28 ).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 t max = 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 t c , i.e., the fraction of those that satisfy | δi (t)|/(2π) ≤ 0.3 Hz for all t ∈ [t c , t max ].For each t c , the fraction f is estimated using 1,000 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 t c → ∞).We show λ max as a function of ε, where solid and dashed curves indicate the average and the maximum, respectively, over perturbations in 1,000 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 t c for β ̸ = (blue), β = (red), and β 0 (black).Note that f for β 0 remains zero for all t c < 10 seconds in all cases.
Isolating converse symmetry breaking.Since real power systems generally have heterogeneity in a i , c ik , 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 P 2j = P 3j 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 more than 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 b 1 and b 3 are necessarily different but can be swapped) is shown to underlie optimal damping.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 i 1 and i 2 , we define this symmetrized matrix P ′ by P ′ i 1 j = P ′ i 2 j ≡ (P i 1 j + P i 2 j )/2 for all j ̸ = i 1 nor i 2 , making it symmetric under the permutation of nodes i 1 and i 2 .To elucidate CSB for each system, we chose 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 .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.New England (10-gen) NPCC (48-gen) Germany (69-gen) U.K. (66-gen) normalized 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 (β i 1 , β i 2 ), where i 1 and i 2 are the nodes whose permutation holds the symmetrized matrix P ′ invariant.In each panel, the red circle marks the optimal (β i 1 , β i 2 ) on the diagonal β i 1 = β i 2 (white dashed line), while the white circle marks the optimal when β i 1 ̸ = β i 2 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 β i 1 along the black dashed line connecting the red and white circles in a, respectively.c Input strength patterns of the nodes i 1 and i 2 for the original matrix P. The color of each arrow indicates |P ij |, normalized by the largest value of |P ij | shown.For nodes i 1 and i 2 , 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 β i 1 = β i 2 , as we interpolate between the original matrix P and its symmetrized version P ′ .Node indexing is described for all four systems in Methods.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.

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 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 known 13 that such systems can be parameterized so that their Jacobian matrices take a form generalizing Eq. ( 2), and they are thus capable of exhibiting 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.Realizing 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 H i , D i , and x int,i used in our stability calculations.The parameters H i and D i are the inertia and damping constants, respectively, that define the effective damping parameter through the relation The parameter x int,i represents the internal reactance of generator i and is used in the calculation of the parameters a i , c ij , 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 toolbox 29 .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 MAT-POWER toolbox 32 , 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 grid 33 , we used the data file (data48em.m)available in the PST toolbox 29 .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 H i and x int,i for all generators (while it assumes D i = 0).We generated D i 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 Earth 35 .
• 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, H i , D i , and x int,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 calculation, 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 model 37 .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 Earth 35 .
• 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 toolbox 32 , is available from our GitHub repository 38 .
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 Zhukovbased aggregation method similar to that described in Ref. 33.In this method, the equivalent generator has inertia constant i H i and damping constant i D i , where the sums are taken over the generators i in the subset.The effective damping parameter of the equivalent generator 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.
in power-grid networks (2021).GitHub: https://github.com/fmolnar-notredame/csb_powergridDOI: 10.5281/zenodo.4437866distinct, and indexed so that The assumption that these eigenvalues are distinct will be crucial in the arguments that follow.When ε = 0, we have B = β = I, so and this implies that β 2 = (ν 2 + ν) = α j , or equivalently ν = −1 ± 1 − α j /α 2 /2, whenever ν is an eigenvalue of J.We thus index ν j ± (ε) so that are not relevant for determining the stability of synchronous states since these eigenvalues are associated with perturbation modes that do not affect synchronization.Thus, the Lyapunov exponent λ max determining the stability is the largest real component among the remaining eigenvalues: For ε = 0, which corresponds to the point β = , it follows from the formula for ν j ± (0) in Eq. (S5) and α 2 < α j , ∀j ≥ 3, that Since eigenvalues are continuous functions of the matrix elements and γ is a continuous function, ν j ± (ε) changes with ε continuously.Thus, for ε ̸ = 0, the eigenvalue ν 2 + (ε) determines the maximum in Eq. (S6) for sufficiently small ε, and hence we have where we factored the characteristic polynomial of J into quadratic factors as (noting that d 1 = 0 always holds because ν = 0 is an eigenvalue of J associated with α 1 = 0).Equations (S8) and (S9) imply that λ max can be viewed as a function of c 2 and d 2 and thus defines a landscape over the (c 2 , d 2 )-plane, while the path γ(ε) in the β-space corresponds a path on the (c 2 , d 2 )-plane defined by the functions c 2 = c 2 (ε) and , we see that the condition for λ max to be decreasing along the path γ(ε) starting at β = in the β-space is equivalent to the condition that the corresponding path (c 2 (ε), d 2 (ε)) starting at (c 2 , d 2 ) = (1, 1/4) immediately enters the following triangular region of the (c 2 , d 2 )-plane: The latter condition implies that the derivatives of c 2 (ε) and d 2 (ε) at ε = 0 satisfy 2d ′ 2 (0) ≥ c ′ 2 (0).We may assume c ′ 2 (0) ≥ 0 without loss of generality.(If not, we simply need to make the change of variable, ε → −ε.)In Sec. 3, we show that, for a generic choice of the matrix P, we have d ′ 2 (0) = 0.This implies that, if λ max decreases along the path, we have c ′ 2 (0) = 0.In general, if c ′ 2 (0) = d ′ 2 (0) = 0, the nonlinearity of c 2 (ε) and d 2 (ε) determines whether λ max decreases along the path in either direction.We further show that the condition c ′ 2 (0) = 0 translates to the path γ being tangent to a specific hyperplane, which we denote by L, at the point β = .
2 Equations relating coefficients c i and d i to parameter ε We now derive equations that define the coefficients c i and d i implicitly as functions c i = c i (ε) and d i = d i (ε).We first rewrite the characteristic polynomial of J as where we used a similarity transformation based on the diagonalization of the matrix P), and D is the diagonal matrix with diagonal elements Along the path γ(ε) = (γ 1 (ε), . . ., γ n (ε)), the components of the matrix B can be expressed as where u iℓ and v iℓ are the ℓth component of the left and right eigenvectors of P associated with the eigenvalue α i , respectively.We now use the definition of the determinant, where the summation σ is taken over all possible permutations σ of indices, and sign(σ) is the sign of permutation σ (i.e., sign(σ) = 1 when σ is an even permutation and sign(σ) = −1 when σ is an odd permutation).With this definition, we can write the coefficient of the ν k term (for k = 1, . . ., 2n − 1; no constant term corresponding to k = 0, since ν = 0 is an eigenvalue) of the polynomial in Eq. (S11) as σ sign(σ) where the summation {k i } is taken over all possible combinations of k i = 0, 1, 2, i = 1, . . ., n.
The function χ is an indicator function defined by χ corresponding to the matrices D, B, and I in Eq. (S11), respectively (δ ij is the Kronecker delta function).Equations (S14) and (S15) provide an expression for the coefficients of det(J − νI) in terms of β, and thus of ε, for a given system and its matrix P.
We now derive a different expression for det(J − νI), this time in terms of the characteristic polynomial coefficients c i and d i , using Eq.(S9).Ignoring the dependence of c i and d i on ε for the moment and regarding them as independent variables, the coefficient of the ν k term can be written as where the summation is the same as in Eq. (S14) and Setting Eqs.(S14) and (S16) equal to each other, we obtain a set of nonlinear equations that must be satisfied by the variables c 1 , . . ., c n , d 1 , . . ., d n , and ε: or, in vector form, F(c 1 , . . ., c n , d 1 , . . ., d n , ε) = 0, (S19) where we have defined F ≡ (F 1 , . . ., F 2n−1 ) T , and the functions F k are given by (S20) This is the equation that implicitly defines the functions c i = c i (ε) and d i = d i (ε).Note that, when c i = 1, d i = α i , and ε = 0 (corresponding to the point β = ), we have implying i.e., Eq. (S18) is satisfied.The functions c i (ε) and d i (ε) defined through Eq. (S20) thus satisfy c i (0) = 1 and d i (0) = α i .
3 Characterizing descending paths on λ max -landscape Here, we will apply the IFT to show that the functions c i (ε) and d i (ε) are continuously differentiable for small ε (and thus define a smooth curve in the neighborhood of the point (1, α i ) in the (c i , d i )-plane), and we determine their first derivatives.The condition under which we can apply the IFT to Eq. (S18) at the point (c 1 , . . ., c n , d 1 , . . ., d n , ε) = (1, . . ., 1, α 1 , . . ., α n , 0) is that the (2n − 1) × (2n − 1) matrix is non-singular, where the elements of G are all evaluated at that point.We note that d 1 is excluded from the set of variables here because d 1 = 0 always holds.We also note that G (and whether it is singular or not) is completely determined by α 2 , . . ., α n , and hence by the matrix P (see examples in Sec. 4).For notational convenience, define x s = c s for s = 1, . . ., n and x s = d s−n+1 for s = n + 1, . . ., 2n − 1. Differentiating Eq. (S20), we find an expression for the (k, s)-element of G: where the summation is defined as in Eq. (S14); we denote ŝ ≡ s − n + 1; the values of a are given by Eq. (S21); and we have defined Since the eigenvalues of P are assumed to be distinct, we have We have numerically verified that this condition holds for all networks considered in the main text.We first seek to show that the matrix G is non-singular by proving that the sth and s ′ th columns of G are linearly independent for all distinct pairs s and s ′ .To do this, we first note that Eq. (S23) simplifies in the special cases k = 1, k = 2, k = 2n − 2, and k = 2n − 1, as follows.
For k = 1, to satisfy i k i = k, we must have k i = 1 for exactly one value of i and k i = 0 for all the others (recall that each k i is either 0, 1, or 2).For the case 1 ≤ s ≤ n, we can derive a simplified formula, but it is not needed below, so we will skip that case here.For the case n + 1 ≤ s ≤ 2n − 1, to have χ s, k, {k i } = 1 in Eq. (S24) we must have k t = 1 for some t ̸ = ŝ and k i = 0 for all i ̸ = t (including i = ŝ).Since there are n − 1 possibilities for t, there are that many nonzero terms in the summation in Eq. (S23), which reduces to For k = 2, the condition i k i = k implies that we either have (a) k i = 2 for exactly one value of i and k i = 0 for all the others, or (b) k i = 1 for two different values of i, and k i = 0 for all the others.For 1 ≤ s ≤ n, we have χ s, k, {k i } = 0 for the terms in Eq. (S23) corresponding to case (a), according to Eq. (S24).For nonzero terms in Eq. (S23) corresponding to case (b), we have k s = 1, k t = 1 with some t ̸ = s, and k i = 0 for all i ̸ = s, t.Putting the two cases together, we obtain A simplified expression for G 2s for n + 1 ≤ s ≤ 2n − 1 can also be obtained, but it is not needed for our purpose.
For k = 2n − 2, satisfying i k i = k requires that either (a) k i = 0 for exactly one value of i and k i = 2 for all the others, or (b) k i = 1 for two different values of i, and k i = 2 for all the others.For 1 ≤ s ≤ n, similarly to the case of k = 2, there is no nonzero term in Eq. (S23) corresponding to case (a), and for the nonzero terms in Eq. (S23) corresponding to case (b), we have k s = 1, k t = 1 with some t ̸ = s, and k i = 2 for all i ̸ = s, t.Putting these two cases together, we obtain For n + 1 ≤ s ≤ 2n − 1, case (b) does not correspond to any nonzero term in Eq. (S23) because k i ̸ = 0 for all i.Case (a), on the other hand, yields exactly one term with k ŝ = 0, leading to For k = 2n − 1, the only way to satisfy i k i = k is to have just one k i = 1 and all the other k i = 2.For 1 ≤ s ≤ n, we must have k s = 1 for a nonzero term, so only one term survives in Eq. (S23), and hence For n + 1 ≤ s ≤ 2n − 1, there is no nonzero term since there is no i for which k i = 0 (see Eq. (S24)), and thus Now we can use Eqs.(S25)-(S30) to show that any pair of columns of G are linearly independent.For 1 ≤ s ≤ n and n + 1 ≤ s ′ ≤ 2n − 1, the last two components (k = 2n − 2 and k = 2n − 1) of the sth and s ′ th column vectors form the two-dimensional vectors (n − 1, 1) T and (1, 0) T , respectively, which are linearly independent.This implies that the full (2n − 1)-dimensional vectors in the sth and s ′ th columns of G are also linearly independent.For 1 ≤ s < s ′ ≤ n, since the last components (k = 2n − 1) of the sth and s ′ th column vectors are both equal to one, it suffices to show that the second component (k = 2) is different in order to establish that they are linearly independent.From Eq. (S26), we have since α s ′ > α s and α i > 0, ∀i, and hence all terms in the summation are positive.Thus, we have G 2s ̸ = G 2s ′ , implying that the sth and s ′ th column vectors are linearly independent.For n + 1 ≤ s < s ′ ≤ 2n − 1, the argument is similar to the case of 1 ≤ s < s ′ ≤ n; it suffices to show that the first component (k = 1) is different, and Eq.(S25) yields where we have defined ŝ′ ≡ s ′ − n + 1 (and recall that ŝ = s − n + 1).Combining all of the above, we have that the sth and s ′ th column vectors of G are linearly independent for all pairs of distinct s and s ′ , which establishes that G is non-singular.
The IFT can now be applied to Eq. (S18) to conclude that the functions c i (ε) and d i (ε) are continuously differentiable.Furthermore, their derivatives satisfy the set of equations obtained by substituting these functions into Eq.(S18) and differentiating both sides with respect to ε: Using Eq. (S20), this can be written as where we have now written the dependence of a on ε explicitly.Setting ε = 0 and using Eq.(S21), we obtain which can be put in vector form as G(x−y) = 0 using the notations x = (x ′ 1 (0), . . ., x ′ 2n−1 (0)) T and y = (y ′ 1 (0), . . ., y ′ 2n−1 (0)) T .Thus, since G is non-singular, we have x = y, and hence x ′ i (0) = y ′ i (0).It then follows that Therefore, G is non-singular (regardless of the value of α 2 > 0).
Example 2: For n = 3, it is convenient to first write out the terms in the first sum of the first line in Eq. (S20) and differentiate them with respect to c i and d i (since the second sum does not contain c i or d i ).The terms from the first sum are Using this, we see that the matrix of partial derivatives is From this, we obtain det(G) = −α 3 (1 − 4α 3 ) 2 /64 = −α 3 (α 2 − α 3 ) 2 /(256α 3 2 ) and see that G is non-singular if and only if 0 < α 2 < α 3 .This condition is indeed satisfied by the 3-generator system considered in the main text (0 < α 2 ≈ 75.5 < α 3 ≈ 178.5).

Fig. 2 :
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 ± b 1 i, a ± b 2 i), and two equal real eigenvalues (a, a).For details on the system, see Methods.

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.

Fig. 5 :
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 |P ij |, with the thickness proportional to |P ij | and the color encoding |P ij | 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.