Collective nonlinear dynamics and self-organization in decentralized power grids

The ongoing transition to renewable energy supply comes with a restructuring of power grids, changing their effective interaction topologies, more and more strongly decentralizing them and substantially modifying their input, output, and response characteristics. All of these changes imply that power grids become increasingly affected by collective, nonlinear dynamic phenomena, structurally and dynamically more distributed and less predictable in space and time, more heterogeneous in its building blocks, and as a consequence less centrally controllable. Here cornerstone aspects of data-driven and mathematical modeling of collective dynamical phenomena emerging in real and model power grid networks by combining theories from nonlinear dynamics, stochastic processes and statistical physics, anomalous statistics, optimization, and graph theory are reviewed. The mathematical background required for adequate modeling and analysis approaches is introduced, an overview of power system models is given, and a range of collective dynamical phenomena are focused on, including synchronization and phase locking, flow (re)routing, Braess’s paradox, geometric frustration, and spreading and localization of perturbations and cascading failures, as well as the nonequilibrium dynamics of power grids, where fluctuations play a pivotal role.

Figure 1

Transmission elements. (a) Simple series admittance. (b) Unified transmission element consisting of an ideal transformer and a $\pi$-equivalent line model.

Figure 2

Schematic of a round-rotor synchronous machine (left image) and a salient pole machine (right image).

Figure 3

Equivalent circuit diagram of (a) a round-rotor synchronous machine and (b) a salient pole machine.

Figure 4

Load-frequency control scheme according to ENTSO-E ([242]).

Figure 5

Topology of a high voltage power grid. The visualization has been produced by Paul Cuffe using the methods described by [46] based on data of the NESTA archive. From [36].

Figure 6

Nose curve. (a) Voltage at the load node as a function of transmitted power for three different power factors: $\tan (\vartheta )=+0.3,0,-0.3$ (from left to right). (b) Voltage at the load node for $\tan (\vartheta )=0$ in the initial state (right) and after half of the connectivity is lost, i.e., $B_{ns}$ is halved (left). Parameters are $S_{\mathrm{base}}=100\mathrm{MW}$, $V_s=1\mathrm{pu}$, and $B_{ns}=10\mathrm{pu}$. Solid (dashed) lines indicate stable (unstable) solutions of the load flow equations.

Figure 7

Factors limiting real power transmission in ac power grids for two different characteristic cases (blue, short line with $l=100\mathrm{km}$; dashed red, long line with $l=200\mathrm{km}$). We consider a single transmission line connecting one power plant to a substation. Shown are the current $|I|$, the voltage magnitude $|V_n|$, and the voltage angle ${\delta }_n$ at the load node as a function of the real-power load $|P_n|$. The dotted line shows the limits $I_{\mathrm{th}}$ and $V_{\lim }$ for the current and the voltage.

Figure 8

The tilted washboard potential [Eq. (98)] as a mechanical analog to the dynamics of a single generator in the second-order model for (a) $P_{\mathrm{mech}}=0.1\mathrm{K}$, (b) $P_{\mathrm{mech}}=0.5\mathrm{K}$, (c) $P_{\mathrm{mech}}=0.9\mathrm{K}$, (d) $P_{\mathrm{mech}}=1.1\mathrm{K}$. The potential ${\mathcal{V}}_{\mathrm{eff}}$ is given in units of K. Green disks, stable fixed points; red crosses, saddle points.

Figure 9

Basin of attraction (red) of the stable fixed point of Eq. (96) in the $\delta$-$\dot{\delta }$ plane. Parameters are $2H/{\omega }_R=1$, $D=0.3$, $K=8$, and (a) $P_{\mathrm{mech}}=1$, (b) $P_{\mathrm{mech}}=2$, (c) $P_{\mathrm{mech}}=4$, and (d) $P_{\mathrm{mech}}=9$. The white dots are the stable and unstable fix points. The black line gives the limit cycle.

Figure 10

Stability map of a single generator coupled to an infinite grid in the second-order model. Three different parameter regions exist: a globally stable fixed point, a globally attractive limit cycle, and a region of coexistence of a fixed point and a limit cycle. Bifurcations between these regions are shown as colored lines. The red dashed line is the approximation (102) for low friction. See [187] and [143] for a more detailed analysis.

Figure 11

Survivability region (red) of the stable fixed point of Eq. (96), with the survival region given by the area between the lines $\dot{\delta }=\pm 5\mathrm{rad}/\mathrm{s}$. The basin of attraction of Fig. 9 is shaded in orange. The parameters are $2H/{\omega }_R=1$, $D=0.3$, $K=8$, and (a) $P=1$, (b) $P=2$, (c) $P=4$, and (d) $P=9$. The white dots are the stable and unstable fixed points. The black line gives the limit cycle.

Figure 12

Desynchronization of the swing equation due to a fluctuating power input. (a),(b) When the input power $P_{\mathrm{mech}}$ fluctuates, the generator can lose synchrony to the grid after an escape time $\tau$. (c) Kramers’s escape rate theory predicts the escape process for Gaussian white noise. The theoretical prediction [Eq. (107) (black lines)] accurately predicts the mean escape times $\overline{\tau }$ obtained from direct numerical simulations averaged over 2000 sample processes (symbols). Results are shown for an extremely loaded connection with $K=1$ and $P=0.95$ and a damping coefficient of $D=8\times {10}^{-5}{\mathrm{s}}^{-1}\times {\omega }_R$ (intermediate damping). From [200].

Figure 13

Multiplicity of sync states: scaling of the number of stable fixed points in a grid with two cycles. The grid consists of two rings with $n+1$ nodes sharing one edge such that $\mathtt{N}=2n$. Edge weights are assumed to be homogeneous and loads vanish ($P_i=0$). The analytic scaling result $0.1576(n^2+2n)$ (solid line) matches the numerically exact results (circles) well. Adapted from [142].

Figure 14

Nonlinear stability of different attractors in a small sample network. (a) Structure of the network with two classes of nodes with injections $P_i=\pm 1$. (b) The basin bifurcation diagram for the system with the synchronized, partially synchronized, and synchronized regimes. From [80].

Figure 15

Topological properties determine the single-node basin stability and survivability. Left panel: survivability (Surv) and basin stability (BS) for nodes in a synthetic power grid network ensemble, showing pronounced differences in the behaviors of various node types. Right panel: node classification according to [168]. Edges that are not in any cycles are trees. The nodes in trees are distinguished as roots (part of a cycle and a tree), sprouts (degree 1 nodes next to a root), leaf (degree 1 nodes not next to roots), and inner tree nodes. Sprouts are further distinguished according to their neighbor degree as dense ($k_{av}>5$) and sparse ($k_{av}\le 5$). Adapted from [168].

Figure 16

Main results of [95]. (a)–(d) Slices of the phase space in the dimensions belonging to a node $k$. Different colors show different asymptotic states. Peach is the sync state, blue is a solitary state, and green is a solitary rotating against its power infeed. (e) Average single-node basin stability for the entire ensemble of networks, for the different types of asymptotic states, as well as the number of desynchronized oscillators following a perturbation (solid purple line). Adapted from [95].

Figure 17

Impact of a link failure in a square lattice with uniform edge weights. (a) Normalized change of the nodal potentials ${\alpha }_n$ for a single failing link located in the center of the network. Both the size of the nodes and the color code represent ${\alpha }_n$. (b) Normalized change of the link flows ${\mathrm{\Delta }F}_e$ for the same topology. Arrows and color represent the direction and strength of flow changes, respectively. From [230].

Figure 18

Left panel: German transmission grid topology (220 and 380 kV) ([144]). Right panel: Double-logarithmic plot of transmitted power change $⟨|\mathrm{\Delta }{F}_{mn}^{(\ell )}|⟩(r)$ as a function of distance $r$ when a transmission line $(\ell )$ between 880 different pairs of nodes $i$ and $j$ was added. We averaged this quantity over all edges with the same distance $r$ to the added line $(\ell )$ and performed an ensemble average over $R=100$ realizations with randomly distributed generators and loads. The power-law expected for a two-dimensional regular grid $\mathrm{\Delta }F\sim r^{-2}$ is shown for comparison (black line). The colored data belong to special subsets of added edges, as discussed in the text. Error bars, 95% confidence level. Adapted from [107].

Figure 19

Localization of the impact of line failures as seen with different distance measures. (a) Flow changes after the outage of a single line (in red) in the IEEE 30-bus test grid. (b) Geodesic distance to the failing link. (c) Rerouting distance to the failing link accurately predicts the spatial distribution of the flow changes. (d),(e) Flow changes $|\mathrm{\Delta }{F}_e|$ vs distance to the failing link $\ell$ on a logarithmic scale. One clearly observes that the rerouting distance used in (e) yields a higher correlation than the ordinary geodesic distance shown in (d). Adapted from [230].

Figure 20

The importance of community structures in the spreading of failures. (a) Community structure of the Scandinavian power grid obtained by spectral clustering. (b) Magnitude of the LODFs for a failing link at the border between southern Norway and Sweden (red line). (c) Magnitude of the line outage distribution factors $|{\mathrm{LODF}}_{j\ell }|$ as a function of the unweighted edge distance of the transmission lines $j,\ell$. The cases in which $j$ and $\ell$ are in the same or different communities are compared. The line gives the median, while the shaded area gives the 25%–75% quantile. The grid data were taken from [99].

Figure 21

Collective effects can lead to a complete reversal of the flow changes in an $N-2$ failure. The color code of the lines indicates the magnitude of flow, with red indicating failing links. (a) Flows in the initial unperturbed grid. (b),(c) Flow changes after the individual failure of two links ($\ell$ and $k$, respectively, in red). In both cases, the flow in the top left link ($e$; boldface font) is greater than in the unperturbed grid: $\mathrm{\Delta }{F}_l^{(\ell )}\approx 0.3$, $\mathrm{\Delta }{F}_l^{(k)}\approx 0.1$. (d) Flows after the simultaneous failure of both links $\ell$ and $k$. The flow in the top left link is smaller in magnitude than in the unperturbed grid and even changes its direction: $\mathrm{\Delta }{F}_e^{(\ell ,k)}\approx -0.7$. From [111].

Figure 22

(a) Real power flows $F_{ab}$ in MW in a test grid based on the Austrian grid topology. The two marked links are investigated in detail. (b) The ratio of the flow and the redundant capacity $F_{a\to b}^{(0)}/K_{a\to b}^{\mathrm{red}}$ provides a powerful indicator for critical links. (c),(d) Illustration of rerouting pathways after the outage of links 1 and 2, respectively. Several rerouting pathways exist for link 1 such that the redundant capacity $K^{\mathrm{red}}$ is rather high, and a line failure is uncritical. In contrast, a failure of link 2 will cause secondary outages. The min cut for the rerouting problem consists of a single link (marked with an arrow). This bottleneck has a free capacity of 864 MW such that a rerouting is impossible. Grid data are based on the work of [103], with a slight adaptation of power injections. Austria was artificially islanded for illustrative purposes.

Figure 23

Fluctuations of power generation and consumption can induce overloads of transmission lines. Left panel: nominal line loads $|{\mu }_{L,\ell }|=|⟨L_{\ell }⟩|$ in the SciGRID model of the German power grid. The nominal operation is determined from one-hour averages of the renewable feed-in and load. The feed-in of the dispatchable power plants and the power flows are then computed via a linear OPF; see Sec. 3b. Right panel: probability of a transmission line overload $\mathcal{P}(|L_{\ell }|\ge 1)$ when the true power injections fluctuate around the nominal values. Vulnerable lines are identified using small values of the indicator $(1-|{\mu }_{L,\ell }|{)}^2/{\sigma }_{\ell }^2$. From [162].

Figure 24

Influence graphs are used to describe the propagation of cascading failures and to identify critical infrastructures. Left panel: power flows in an elementary test grid consisting of six nodes. Right panel: influence graph summarizing the likelihood that the failure of one line will induce the failure of another line. The nodes of this influence graph correspond to the lines in the original grid as identified by their end points. Adapted from [98].

Figure 25

Statistical analysis of cascading failures. The coreness of a link determines its vulnerability to fail (upper panels) or to be disconnected (lower panels) in a cascade of failures. Shown is the fraction of lines that (c) fail or (e) get disconnected in at least one cascade. (d) Probability of failure $⟨{\mathcal{P}}_{\ell }^{(p)}⟩$. (f) Probability of disconnection $⟨{\mathcal{P}}_{\ell }^{(c)}⟩$ averaged over all cascades and operating conditions. Simulations were carried out for the three North American grids under various operating conditions; see the text for details. Adapted from [269].

Figure 26

Statistics of power outages in North American power grids (1984–1998). The empirical probability of a power outage as a function of the outage size on a logarithmic scale. The data suggest a power-law distribution with exponent $1.3\le \beta \le 2.0$. From [49].

Figure 27

Transient effects can increase the vulnerability of a grid to cascading failures. (a) An elementary test grid with two effective generators ($P=+1.5{\mathrm{s}}^{-2}$; squares) and three effective consumers ($P=-1{\mathrm{s}}^{-2}$; circles). All lines have $K_{ij}=10/6$ and $F_{ij}^{\max }=1$. At time $t=1\mathrm{s}$, the line $(2,4)$ fails. (b) In a full dynamical study, line ratings are exceeded transiently, leading to secondary line outages. Shown is the number of failed lines as a a function of time. (c) In a quasistatic picture, the grid immediately relaxes back to a stable steady state that respects the line ratings. Shown are the steady-state line flows $|F_{\ell }|$ before and after the initial failure. (d) Time evolution of the line flows $|F_{\ell }(t)|$ taking into account both the continuous dynamics and the outage of lines. From [202].

Figure 28

Braess’s paradox: loss of the synchronous fixed point due to grid extensions. (a)–(c) Topology of the network. The vertices generate or consume the power $P_i=\pm P$, and the transmission lines have a capacity $K$. We consider the cases where (b) the capacity of one line is doubles and (c) a new line is added. (d)–(f) Dynamics of the grid given by the equations of motion (111) starting from the initial state ${\theta }_i(0)={\dot{\theta }}_i(0)=0$. The initial grid relaxes to a synchronized fixed point, while synchronization becomes impossible after the grid extensions. The parameters are $K=1.03P$, $D=P$, and $J=1$. Adapted from [260].

Figure 29

Non-Gaussian frequency deviations with a reduced number of extreme events. (a) Distribution of frequencies around a set frequency (50 Hz) in the Continental European grid (yellow) deviates from the best Gaussian distribution (parabolic blue solid line). In particular, the probability (per Hz) of large events of more than 0.1 Hz frequency deviations is increased by more than a factor of 10 compared to the Gaussian fit. (b) Counts of “extreme events,” i.e., the number of frequency deviations larger than 0.1 Hz resolved for each minute, within full hours, accumulated for 2017 and 2011, respectively. Compiled using data from [186] and the methods of [201].

Figure 30

Arrival time $t^{*}=(t-t_0)/\tau$ of a short pulse disturbance in units of $\tau =J/D$ vs the geometric distance $r^{*}=r/a$ in units of the transmission line length $a$. The green rhombi represent the numerical results for (a) a square grid and (b) German transmission grid topology: $KJ/D^2{\omega }_R={10}^5$, a realistic value for high voltage transmission grids. Also shown are the ballistic equation (164) with the fitted velocity $v$ (red line) and the analytic velocity equation (165) (red dashed line). The black squares indicate node-resolved arrival times $t_i^{*}$. Unfilled squares indicate arrival times, averaged over all nodes at the same distance (c) for a Cayley tree with branching $b=2$, and (d) for the German transmission grid. In (c) and (d) $KJ/D^2{\omega }_R=10$. Diffusive equation (166) fitted (pink line) and with the analytical equation (167) (pink dashed line). Also shown are the ballistic equation (164) fitted (red line) and with analytic velocity (red dashed line). $\sigma =P/P_c=0.1$, with critical power $P_c$ for the respective grids. Adapted from [232].

Figure 31

Normalized Fiedler vector intensity $|{\varphi }_{2,i}|/{\max }_n|{\varphi }_{2,n}|$ relative to its maximum (shown in black) for the German transmission grid (left panel) featuring maximal intensity at the south and north borders and for a random grid with strong localization at a single node (black; at the bottom right of the right panel). Adapted [234].

Figure 32

(a) Variance of increments $⟨{\widehat{D}}_{\mathrm{\Delta }t}{\omega }_i^2⟩$ for a chain of $\mathsf{N}=50$ oscillators. The straight lines have a slope $m=-1/\xi$ with Eq. (169). (b) Eigenvalues of Laplacian ${\mathrm{\Lambda }}_n$ for various inertia $J$. (c) Fits of slope $m$ (orange) converge to an analytical $m=-1/\xi$ [Eq. (169)] (blue) with increasing inertia $J$. Vertical dashed line, $J_c$; error bars, $2\sigma$ confidence bound. From [89].

Figure 33

Kurtosis $k$ of frequency increment PDF in a chain of $\mathsf{N}=100$ nodes. The frequency increment distribution $p({\widehat{D}}_{\mathrm{\Delta }t}{\omega }_i)$ deforms only slowly toward an almost Gaussian distribution. Insets (from left to right): $i\in \{2,20,50\}$, $\theta =0.01\mathrm{s}$. From [89].

Figure 34

Localization, distributed resonances, and bulk oscillations. In a small network example (a), the joint implications on response strength (b) of driving one node (labeled 0) are illustrated as a function of driving frequencies [(c)–(e)] and at four selected response nodes (labeled 1–4). Whereas at high driving frequencies responses are localized (algebraically in frequency, exponentially in internode distance), they are irregular across frequencies and among nodes in a resonance regime of intermediate driving frequencies. At low frequencies, globally homogeneous bulk oscillations emerge where the entire grid follows the driving signal (even if at only one node). All response strengths are quantified in terms of the amplitude of the frequency response relative to their low frequency limit as $\omega \to 0$. From X. Zhang; cf. [275].