# Rare desynchronization events in power grids: On data implementation and
dimensional reductions

Tim Ritmeester and Hildegard Meyer-Ortmanns  School of Science, Jacobs
University Bremen, Campus Ring 1, 28759 Bremen, Germany  t.ritmeester@jacobs-
university.de<EMAIL_ADDRESS>
###### Abstract

We discuss the frequency of desynchronization events in power grids for
realistic data input. We focus on the role of time correlations in the
fluctuating power production and propose a new method for implementing colored
noise that reproduces non-Gaussian data by means of cumulants of data
increment distributions. Our desynchronization events are caused by overloads.
We extend known and propose different methods of dimensional reduction to
considerably reduce the high-dimensional phase space and to predict the rare
desynchronization events with reasonable computational costs. The first method
splits the system into two areas, connected by heavily loaded lines, and
treats each area as a single node. The second method considers a separation of
the timescales of power fluctuations and phase angle dynamics and completely
disregards the latter. The fact that this separation turns out to be
justified, albeit only to exponential accuracy in the strength of
fluctuations, means that the number of rare events does not sensitively depend
on inertia or damping for realistic heterogeneous parameters and long
correlation times. Neither does the number of desynchronization events
automatically increase with non-Gaussian fluctuations in the power production
as one might have expected. On the other hand, the analytical expressions for
the average time to desynchronization depend sensitively on the finite
correlation time of the fluctuating power input.

[figure]style=plain, subcapbesideposition = center

  * November 2022

##  1 Introduction

Power grids are confronted with various sources of stochastic fluctuations.
Fluctuations refer to production, consumption, and prices on the energy
market, where supply and demand are traded to cure imbalances between forecast
and actual needs of energy supply. Variations in supply and demand have always
existed, but in view of decarbonization, the increase in the contribution of
renewable energies has dramatically increased fluctuations due to larger and
more localized variability of wind and solar energy [1]. As a result, power
grids become increasingly affected by collective, nonlinear phenomena, which
are less predictable in space and time, more heterogeneous and less well
controllable. For a recent review see [2].

Wind and solar energies are currently fast-growing energy sectors in Europe.
In addition, more renewables in the grid go along with stronger fluctuations,
less inertia and less damping. How fluctuations induced by renewables
challenge the stability and quality of electrical power grids was analyzed in
[3], addressing the regime in which desynchronization events are not rare.
Here, feed-in fluctuations were modelled with realistic features of non-
vanishing temporal correlations, the Kolmogorov power spectrum and
intermittent increments. Less inertia and damping are expected to further
challenge the grid stability. The impact of increased fluctuations, reduced
inertia and reduced damping on the grid stability was studied in [4]. The
fluctuations were implemented as Gaussian white noise and assumed to be small.
The analysis was based on employing Kramer’s escape rate. Accordingly, inertia
has a minor effect on the escape rate, entering the prefactor of the
exponential, while the damping occurs in the exponent with a stronger effect.
The escape routes were shown to typically proceed in the vicinity of saddle
fixed points with a low potential barrier, caused by a single overloaded link.
As shown in [5], a stable fixed point of the swing equations can be only lost
via an inverse saddle node bifurcation. As long as the phase differences along
all transmission lines are smaller or equal $\pi/2$, which is the case under
normal grid operation [6, 7], the loss of a stable fixed point is equivalent
to an overload of one or more transmission lines. It is this case that we
consider later in this paper. However, if the condition on the phase
differences is violated, the bifurcation of a fixed point is generally not
associated with an overload.
Histograms of power increments of wind and solar data show clearly non-
Gaussian tails. Non-Gaussian fluctuations were considered in [8] to determine
rare events of power grid desynchronization, where the model was based on the
swing equations under colored noise. The computed rate of desynchronization
showed that higher moments of noise enter at specific powers of the coupling.
Their effect is either increasing or decreasing the rate of desynchronization
events, depending on how the noise statistics match the statistics of the
Fiedler mode, which is the slowest mode in the system. Their work makes use of
a WKB-approach to determine the average time to desynchronization. For earlier
work that uses the WKB approach for classical stochastic systems see, for
example, [9, 10]. We will also employ the WKB-approach, but use different
versions of dimensional reduction for obtaining our predictions.
When strong fluctuations in the production propagate through the grid, non-
Gaussian tails in the distribution of local grid-frequency data have been
observed, which themselves may induce large-scale blackouts in the worst case.
Here the authors of [11] focused on the role of time correlations in the
fluctuations. According to their results, noise disturbances with finite
skewness and positive kurtosis propagate through the entire grid when the
noise correlation time is larger than the network’s intrinsic time scales.
Otherwise, for rapidly fluctuating noise, the disturbances decay over short
distances. We will also focus on the role of time correlations in the non-
Gaussian noise, but pursue their impact on the average time to
desynchronization events.
The paper is organized as follows. In Sec. 2 we present the model in terms of
the swing equations with stochastic power input, described as non-Gaussian
colored noise to reproduce in particular real data of wind power increment
distributions. As network we choose a coarse-grained version of the Brazilian
grid with prototypical values of inertia and damping, while the power
production is artificially increased to simulate a heavily loaded grid at risk
of desynchronization. Examples of typical histograms of power increments are
shown, as measured from wind and solar data. Sec. 3 is devoted to the data
implementation. We propose a new approach, the so-called Fourier-method, and
compare it with a compound Poisson process, previously used, for example, in
[8]. The Fourier-method is analytically more convenient, and more accurately
reproduces renewable generation data. Both approaches can lead to different
predictions of the desynchronization times, depending on how the data are
implemented on the time scale of one integration step. In Sec. 4 we derive
various versions of dimensional reduction of the equations of motion, first
applied to the swing equations, next to Hamilton’s equations of motion that
show up in the WKB-approach. We reduce the phase space of Hamilton’s equations
for the full grid to twelve dimensions (in the synchronized subgraph
approximation), and to two dimensions (in the overload approximation). In Sec.
5 we explain the methods of finding the solutions. Sec. 6 contains a
comparison of the results for the average desynchronization times for
different versions of data implementation and approximation schemes. We
demonstrate the accuracy of these approximations by comparing their results to
those obtained for the full system. We also discuss the role of the skewness
of the power increment distributions. Though derived differently, our results
confirm that the number of rare desynchronization events may increase or
decrease due to more renewables, depending on where the fluctuations enter the
grid. Moreover, we find that neither less inertia nor less damping endangers
the grid stability per se with respect to desynchronization. The role of time
correlations in the fluctuations is clearly non-negligible, correlation times
are related to the time scales of inertia and damping. Sec. 7 contains our
summary and conclusions.

##  2 The model

After introducing the model (Sec. 2.1) and discussing the role of time
correlations in Sec. 2.2, we give some details about the Brazilian grid as a
grid that shares prototypical features in view of realistic data for
production, consumption and inertia (Sec. 2.3), followed by examples of real
data sets (Sec. 2.4).

###  2.1 Swing equation under realistic data input

If one neglects resistive power losses, assumes constant voltage magnitudes
and phase velocities much smaller than the grid frequency
($|\dot{\phi}|\ll\omega$), the dynamics of the voltage phase angles
$\bm{\phi}$ are modeled by the swing equations [5, 6, 12]:

$\displaystyle\frac{\operatorname{\text{d}\\!}\phi_{i}}{\operatorname{\text{d}\\!}t}$
$\displaystyle=v_{i}$ $\displaystyle
M_{i}\frac{\operatorname{\text{d}\\!}v_{i}}{\text{d}t}+\gamma_{i}v_{i}$
$\displaystyle=P_{i}+p_{i}(t)+\sum_{j}K_{ij}\sin(\phi_{j}-\phi_{i}),$ (1)

Here $M_{i}$ and $\gamma_{i}$ are respectively the inertia and the damping
constants of the oscillators, $K_{ij}$ is the susceptance of the power line
between nodes $i$ and $j$ (equal to $0$ if nodes $i$ and $j$ are not
adjacent). $K_{ij}\sin(\phi_{j}-\phi_{i})$ equals the flow $f_{ij}$ of power
from node $j$ to node $i$. For a detailed derivation of why the swing
equations are representative for modeling the grid dynamics (without including
the voltage dynamics) we refer to the work of [5, 6, 12]. The constants
$P_{i}$ are deterministic power injections with $\sum_{i}P_{i}=0$, and
$p_{i}(t)$ are unforeseen fluctuations in power input. Fluctuations in power
input are both non-Gaussian and colored [13, 14, 15].

A convenient noise model which incorporates both of these properties is the
Langevin-like dynamics (following [8, 16]):

$\displaystyle\frac{\text{d}p_{i}(t)}{\text{d}t}$ $\displaystyle=-\alpha
p_{i}(t)+\sigma_{i}\cdot\xi_{i}(t),$ (2)

where $\xi_{i}(t)$ is stationary (generally non-Gaussian) white noise with
$\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$,
and where $\sigma_{i}$ controls the strength of the noise. As the solution to
Eq. 2 is
$p_{i}(t)=\sigma_{i}\int_{-\infty}^{t}\operatorname{\text{d}\\!}t^{\prime}\,\exp(-\alpha(t-t^{\prime}))\xi_{i}(t^{\prime})$,
Eq. 2 leads to exponentially decaying time correlations $\frac{\langle
p_{i}(t)p_{i}(t^{\prime})\rangle}{\langle
p_{i}(t)p_{i}(t)\rangle}=\exp\big{(}-\alpha(t-t^{\prime})\big{)}$, so
$1/\alpha$ sets the correlation time of the fluctuations. This exponential
decay may not really reflect realistic time correlations in feed-in power
fluctuations, but provides a tractable approach to pursue the impact of a
finite (long) correlation time in contrast to a vanishing one as for white
noise. Furthermore this model for the noise allows also to reproduce the
measured increment distributions of power input as generated by real wind
fluctuations, though their power spectrum is not reproduced. The authors of
[3] found that in the regime where desynchronization events are common, non-
Gaussianity has only a small effect on desynchronization times; our approach
will clarify the more important role of non-Gaussian effects in the regime
where such events are rare.

###  2.2 Role of time correlations in power fluctuations

As our focus is on the impact of realistic parameter choices on the grid
dynamics, the role of the correlation time of fluctuations (here of the power
production, with time correlations parameterized by $1/\alpha)$ should be
considered. In the per-unit system power is dimensionless, and $M_{i}/P_{i}$
and $\gamma_{i}/P_{i}$ have units of $\text{time}^{2}$ and time, respectively.
The influence of the correlation time $1/\alpha$ on the dynamics of the swing
equations can be made explicit if we consider a rescaled time
$\tilde{t}\equiv\alpha t$, and define:

$\displaystyle\tilde{\phi}_{i}(\tilde{t})$
$\displaystyle\equiv\phi_{i}(\tilde{t}/\alpha)=\phi_{i}(t)\,,$ (3)
$\displaystyle\tilde{p}_{i}(\tilde{t})$ $\displaystyle\equiv
p_{i}(\tilde{t}/\alpha)=p_{i}(t)\,,$ (4)
$\displaystyle\tilde{\xi}_{i}(\tilde{t})$
$\displaystyle\equiv\alpha^{-1/2}\xi_{i}(\tilde{t}/\alpha)=\alpha^{-1/2}\xi_{i}(t)\,,$
(5)

then using
$\frac{\text{d}\tilde{\phi}_{i}(\tilde{t})}{\text{d}\tilde{t}}=\frac{1}{\alpha}\frac{\text{d}\phi_{i}(t)}{\text{d}t}$
and
$\frac{\text{d}^{2}\tilde{\phi}_{i}(\tilde{t})}{\text{d}\tilde{t}^{2}}=\frac{1}{\alpha^{2}}\frac{\text{d}^{2}\phi_{i}(t)}{\text{d}t^{2}}$,
the swing equations (Eq. 1) become:

$\displaystyle(\alpha^{2}{M}_{i})\cdot\frac{\text{d}^{2}\tilde{\phi}_{i}}{\text{d}\tilde{t}^{2}}+(\alpha{\gamma}_{i})\cdot\frac{\text{d}\tilde{\phi}_{i}}{\text{d}\tilde{t}}$
$\displaystyle=P_{i}+\tilde{p}_{i}(\tilde{t})+\sum_{j}K_{ij}\sin(\tilde{\phi}_{j}-\tilde{\phi}_{i})\,.$
(6)

Furthermore, using that
$\frac{\text{d}\tilde{p}_{i}(\tilde{t})}{\text{d}\tilde{t}}=\frac{1}{\alpha}\frac{\text{d}{p}_{i}(t)}{\text{d}t}$
and
$\delta(\tilde{t}_{1}-\tilde{t}_{2})=\delta(\alpha({t}_{1}-{t}_{2}))=\frac{1}{\alpha}\delta({t}_{1}-{t}_{2})$
turns Eq. 2 into:

$\displaystyle\frac{\text{d}\tilde{p}_{i}}{\text{d}\tilde{t}}=-\tilde{p}_{i}+(\alpha^{-1/2}\sigma_{i})\cdot\tilde{\xi}_{i}(\tilde{t})$
(7)

with
$\langle\langle\tilde{\xi}_{i}(\tilde{t})\tilde{\xi}_{i}(\tilde{t}^{\prime})\rangle\rangle=\delta(\tilde{t}-\tilde{t^{\prime}})$
and
$\langle\langle\tilde{\xi}_{i}(\tilde{t}_{1})\dots\tilde{\xi}_{i}(\tilde{t}_{n})\rangle\rangle=\big{(}\alpha^{\frac{n}{2}-1}{\Gamma}_{n}^{\xi_{i}}\big{)}\cdot\delta(\tilde{t}_{1}-\tilde{t}_{2})\dots\delta(\tilde{t}_{1}-\tilde{t}_{n})$,
where $\langle\langle...\rangle\rangle$ denote the cumulants (Sec. 3, [16,
17]), here of $\tilde{\xi}_{i}$. Therefore, setting the correlation time to
$1$ by rescaling time according to $t\rightarrow\alpha t$ sends
$M_{i}\rightarrow\alpha^{2}M_{i}$, $\gamma_{i}\rightarrow\alpha\gamma_{i}$,
$\sigma_{i}^{2}\rightarrow\sigma_{i}^{2}/\alpha$ and
${\Gamma}_{n}^{\xi_{i}}\rightarrow\alpha^{\frac{n}{2}-1}{\Gamma}_{n}^{\xi_{i}}$,
with the combinations $\alpha^{2}M_{i}$ and $\alpha\gamma_{i}$ invariant.
Increasing the timescale $1/\alpha$ of the noise is thus equivalent to
decreasing the inertia $\bm{M}$ and damping $\bm{\gamma}$. The right-hand side
of Eq. 2 is, in the per unit system, dimensionless, and therefore invariant
under a change in $1/\alpha$. The rescaled strength of noise increases for
increasing $1/\alpha$, while the higher order cumulants of non-Gaussian noise
get increasingly suppressed. This means, for non-Gaussian noise, the noise
terms become closer to Gaussian if $1/\alpha>1\,\text{s}$, and vice-versa. In
Sec. 6 this will be manifest in the results for the average time to
desynchronization.
A further discussion comparing $M_{i}$ and $\gamma_{i}$ to the timescale of
the noise is given in [11], where it is argued that typical timescales of
fluctuations are so long ($1/\alpha=100-300\;\text{s}$) that the impact of
$\bm{M}$ and $\bm{\gamma}$ on the propagation of non-Gaussian disturbances can
typically be neglected. The results of [11] are for typical fluctuations; our
results are derived for rare large fluctuations, thus complementing their
results.

###  2.3 The Brazilian grid as a prototypical case

As a realistic test case for the stability of the fixed point under stochastic
fluctuation we use the coarse-grained data on the Brazilian power network
presented in [18, 19, 20]. The network has $6$ nodes, realistic and spatially
heterogeneous parameters, and contains both synchronous generators/loads and
inductive generators/loads. The network topology is shown in Fig. 1. The
networks of (b) and (c) have the values of $P_{i}$, $\gamma_{i}$ and $M_{i}$
of nodes 3 and 4 swapped as compared to (a). The networks contain two nodes
with high power generation/consumption; nodes $4$ and $6$ in (a) and nodes $3$
and $6$ in (b) and (c). In (a) and (b) we artificially add stochastic
fluctuations $\bm{\xi}$ to the power input at these two nodes (crossed in Fig.
1) to simulate a fraction of the input being fluctuating rather than constant.
In (c) we add power fluctuations to all nodes. The dotted lines are the most
heavily loaded lines (with the smallest $\frac{K_{ij}-f_{ij}}{K_{ij}}$), one
in (a), two in (b) and (c). If the load on the grid is increased, the system
undergoes a bifurcation in which the system splits into two areas (with blue
and pink nodes, respectively) on either end of the line(s), which remain
synchronized internally (Sec. 4.1). It should be noticed that a splitting into
two areas (the blue and pink ones in Fig. 1) within which the nodes remain
internally synchronized are well known in power engineering and called a
system split. They occur also in much larger grids, such as the European one,
where such a split happened in January 2021 [7], see Fig. 2. We artificially
increase the load (by uniformly scaling all power inputs) on the network to
model a situation in which the network is operating close to capacity, with a
realistic safety margin to overload ($\frac{K_{ij}-f_{ij}}{K_{ij}}\approx
10\%$ for the most loaded line [21], where $f_{ij}$ is the power flow from $i$
to $j$, see Sec. 4.1). This network provides a generic test case in the sense
that approximate ratios between power input, inertia and damping are typical
for power grids; see [6, 12]. It is small enough that we can computationally
analyze it in detail, but large enough to provide interesting results.

(a)

(b)

(c)

Figure 1: a) Coarse-grained network based on the Brazilian power grid [18, 19,
20] (Sec. 2.3). The dotted lines are the most heavily loaded lines, one in
(a), two in (b) and (c); after a desynchronization event the system splits
into two areas (with blue and pink nodes, respectively). The two crossed nodes
in (a) and (b) contain large power generation/load, to which we add stochastic
fluctuations. In c) stochastic fluctuations (with the same variance) are added
to all nodes. The networks in (b) and (c) differ from a) as the values of
$P_{i}$, $\gamma_{i}$ and $M_{i}$ for the $3$rd and $4$th nodes are swapped.

As mentioned before, the power inputs are all artificially scaled to produce a
grid closer to overload, to model a heavily loaded grid. The resulting set of
parameters is given in Appendix A. After swapping the values of $P_{i}$,
$M_{i}$ and $\gamma_{i}$ for the third and fourth nodes in Fig. 1(b) and (c),
the power input is again rescaled to a safety margin of $9.1\%$. Fluctuations
are artificially added at node $4$ and node $6$ for Fig. 1(a), at nodes $3$
and $4$ for (b), and at all nodes in (c). For (a), the line connecting node
$4$ and node $5$ (with $K_{45}=13.624$) is the one closest to overload, as
indicated in Fig. 1(a) with a dotted line. In (b) and (c) the most heavily
loaded lines are those connecting node $5$ to nodes $1$ and $2$, with
$K_{15}=4.082$ and $K_{25}=4.444$. For heavily loaded grids, the safety
margin, $(K_{ij}-f_{ij})/K_{ij}$ for the most heavily loaded line, is
generally on the order of $5\%-30\%$ [6, 21]. Note that inertia $\bm{M}$ and
damping $\bm{\gamma}$ are roughly proportional to the power input at the
respective nodes, but about an order of magnitude smaller than $\bm{P}$.

(a)

(b)

Figure 2: A system split that occurred in Europe in January 2021 [7]. The two
areas, indicated in (a), lost synchronization between each other due to an
overload of the red and green lines connecting the two areas indicated in (b),
with the phase angles on either end of the lines reaching a difference of
$\frac{\pi}{2}$. Figures included (with permission) from [7].

The calculations that we performed to determine the parameters of the
Brazilian grid follow [12], they are not specific to the Brazilian grid. In
general, also for different grids like IEEE test-cases [18, 19, 20], the
(order of magnitude of the) ratios between $\bm{P}$ and $\bm{M}$ and
$\bm{\gamma}$ will be similar to the ones used here. What is furthermore
common to power grids, including the Brazilian grid, is that most nodes have
only little generation/load, and a few nodes have a very high generation or
very high load, rather than power generation/load being spread relatively
evenly over the grid. With increasing renewables, however, more inductive and
less synchronous machines will be attached, leading to an associated decrease
in the values of $M_{i}$ and $\gamma_{i}$ as compared to the total power input
at that node.

###  2.4 Wind versus solar data

Our noise terms $\xi_{i}$ will be constructed such that the resulting
increments of power input $p_{i}(t+\tau)-p_{i}(t)$ for a given time interval
$\tau$ match a realistic power input distribution. For the power input
increments we will consider both Gaussian distributions and distributions
taken from data on wind and solar power generation [13, 22] (representing a
month of power generation in Bremen and Oldenburg, respectively), with data
resolution $\tau=1\,\text{s}$. We will vary the strength of the fluctuations
by rescaling the power generation. Histograms for the distribution of power
increments from [22] are given in Fig. 3. For estimating desynchronization
times we consider only wind data and will approximate histogram (a) with
constructed increments $p_{i}(t+\tau)-p_{i}(t)$. For a comparison of the data
implementation procedure we consider also solar data. Moments are
straightforward to calculate from the data, cumulants can then be calculated
from the moments according to e.g. [23]. The first $8$ cumulants of these
distributions are shown in Table 1.

Table 1: The $n-th$ cumulants of the increments $p(t+\tau)-p(t)$ for a Gaussian distribution and for the wind and solar power generation data corresponding to Fig. 3, all standardized for comparison (i.e. rescaled to have zero mean and unit standard deviation). The values show that higher order cumulants can become quite large. Distribution / n | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ | $8$
---|---|---|---|---|---|---|---|---
Gaussian | $0$ | $1$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$
Wind | $0$ | $1$ | $4.198$ | $141.5$ | $2.826\cdot 10^{3}$ | $1.113\cdot 10^{5}$ | $3.847\cdot 10^{6}$ | $1.836\cdot 10^{8}$
Solar | $0$ | $1$ | $-2.497$ | $124.6$ | $-898.7$ | $3.992\cdot 10^{4}$ | $-4.310\cdot 10^{5}$ | $2.135\cdot 10^{7}$

The authors of [11] explicitly analyze the time correlations $\frac{\langle
p(t)p(t^{\prime})\rangle}{\langle p(t)^{2}\rangle}$ for the wind and solar
data of [22] and find that correlations typically decay on a timescale of
$2-5$ minutes. Accordingly, we here set $1/\alpha=100\,\text{s}$, and also
consider $1/\alpha=300\,\text{s}$ in Sec. 6.

(a)

(b)

Figure 3: Histogram of $1$-second power increments meant to be approximated by
$p(t+\tau)-p(t)$, with $10^{3}$ bins, for (a) wind data and (b) solar data
(for the brightest hour of the day) from [22], showing non-Gaussian tails.
Standardized for comparison (i.e. rescaled such that both figures have zero
mean and unit standard deviation). These histograms should be reproduced by
our data implementation.

##  3 Data Implementation

Real data sets of power production are often presented by the corresponding
distribution of power increments, usually showing significantly broader tails
than Gaussian distributions. In the past, fluctuations in production and
consumption were often modelled as Gaussian white noise. To overcome this
approximation, we develop a new version of implementing colored noise (Sec.
3.1). We compare this approach with an implementation as a compound Poisson
process (Sec. 3.2) and point to the importance of time correlations in the
fluctuations.

After all, we want to approximate data in terms of fluctuating increments. To
reproduce the measured increments, we should construct a corresponding
stochastic process $\xi_{i}(t)$ for which the resulting distribution of
$p_{i}(t+\tau)-p_{i}(t)$ accurately reproduces the distribution of data
$\text{pr}(\text{increment data})$. We can either quantify $\xi_{i}(t)$
directly, or, equivalently, quantify it by the distribution of its increments
$Z_{\Delta t,i}\equiv\int_{t}^{t+\Delta
t}\xi_{i}(t^{\prime})\text{d}t^{\prime}$ for any given $\Delta t$.

Let us start with the first option. The noise term $\xi_{i}(t)$ is in general
singular, such that we cannot define an instantaneous probability distribution
of $\xi_{i}$ [16]. Instead we can specify $\xi_{i}$ by its cumulants
$\Gamma^{\xi_{i}}_{n}$ [16, 23], defined by:

$\displaystyle\langle\langle\xi_{i}(t_{1})\cdot\dots\cdot\xi_{i}(t_{n})\rangle\rangle=\Gamma^{\xi_{i}}_{n}\delta(t_{1}-t_{2})\cdot\dots\cdot\delta(t_{1}-t_{n}),$
(8)

where we will always rescale the noise term such that $\Gamma^{\xi_{i}}_{1}=0$
and $\Gamma^{\xi_{i}}_{2}=1$. As defined in Sec. 2.2, the double angular
brackets $\langle\langle\dots\rangle\rangle$ denote cumulants rather than
moments, e.g. $\langle\langle x^{2}\rangle\rangle=\langle x^{2}\rangle-\langle
x\rangle^{2}$. Here we will always assume that all cumulants are finite, and
that the cumulant generating function
$C^{\xi_{i}}(J)=\sum_{n=1}^{\infty}\Gamma^{\xi_{i}}_{n}J^{n}/n!$ exists for
all real $J$. For a more detailed description of white noise processes we
refer to the discussions in [16, 17]. The cumulants $\Gamma^{\xi_{i}}_{n}$ are
constructed such that increments of the power input $p_{i}(t+\tau)-p_{i}(t)$
for a given time interval $\tau$ match a realistic power input distribution.
The cumulants $\Gamma^{\xi_{i}}_{n}$ can be related to the cumulants of
$p_{i}(t+\tau)-p_{i}(t)$ by considering that Eq. 2 gives
$p_{i}(t)=\sigma_{i}\int_{-\infty}^{t}\operatorname{\text{d}\\!}t^{\prime}\,\exp(-\alpha(t-t^{\prime}))\xi_{i}(t^{\prime})$,
that the cumulants of a sum are the sum of the cumulants [16], and that
multiplication by a constant multiplies the $n$-th cumulant by that constant
to the power $n$. Using Eq. 8 we then get that
$\Gamma_{n}^{p_{i}(t+\tau)-p_{i}(t)}$ are related to $\Gamma^{\xi_{i}}_{n}$
according to:

$\displaystyle\Gamma_{n}^{p_{i}(t+\tau)-p_{i}(t)}$
$\displaystyle=\Gamma^{\xi_{i}}_{n}\cdot\sigma_{i}^{n}\int_{-\infty}^{t+\tau}\operatorname{\text{d}\\!}t^{\prime}\,\Big{[}\exp(-\alpha(t+\tau-t^{\prime}))-\mathbbm{1}(t^{\prime}\leq
t)\exp(-\alpha(t-t^{\prime}))\Big{]}^{n}$ (9)
$\displaystyle=\Gamma^{\xi_{i}}_{n}\cdot\sigma_{i}^{n}\Big{(}\int_{-\infty}^{t}\operatorname{\text{d}\\!}t^{\prime}\,\exp(-n\alpha(t-t^{\prime}))\big{[}\exp(-\alpha\tau)-1\big{]}^{n}$
(10)
$\displaystyle\qquad+\int_{t}^{t+\tau}\operatorname{\text{d}\\!}t^{\prime}\,\exp(-n\alpha(t+\tau-t^{\prime}))\Big{)}$
$\displaystyle=\Gamma^{\xi_{i}}_{n}\cdot\frac{\sigma_{i}^{n}}{n\alpha}\Big{(}\big{[}\exp(-\alpha\tau)-1\big{]}^{n}+\big{[}1-\exp(-n\alpha\tau)\big{]}\Big{)}.$
(11)

Note that $\Gamma_{n}^{p_{i}(t+\tau)-p_{i}(t)}$ are constructed directly from
the dynamics of $p(t)$ according to Eq. 2 including the parameter $\alpha$ and
therefore are the quantities that should be primarily identified with the
data. To get the statistics of $p_{i}(t+\tau)-p_{i}(t)$ to approximate the
power increments found from the data, i.e.
$\Gamma_{n}^{p_{i}(t+\tau)-p_{i}(t)}\stackrel{{\scriptstyle\approx}}{{\rightarrow}}\Gamma_{n}^{\text{data}}$
where $\stackrel{{\scriptstyle\approx}}{{\rightarrow}}$ means ’should
approximate’, we can then reverse the relation of Eq. 11 to get that
$\Gamma_{n}^{\xi_{i}}$ should be chosen such that:

$\displaystyle\sigma_{i}^{n}\Gamma_{n}^{\xi_{i}}$
$\displaystyle\stackrel{{\scriptstyle\approx}}{{\rightarrow}}\Gamma_{n}^{\text{data}}\cdot
n\alpha\cdot\Big{(}\big{[}\exp(-\alpha\tau)-1\big{]}^{n}+\big{[}1-\exp(-n\alpha\tau)\big{]}\Big{)}^{-1}\,.$
(12)

The alternative characterization of $\xi_{i}(t^{\prime})$ is by the
distribution of its increments $Z_{\Delta t,i}\equiv\int_{t}^{t+\Delta
t}\xi_{i}(t^{\prime})\text{d}t^{\prime}$ for any given $\Delta t$. From here
on we always consider $\Delta t$ to be measured in units of seconds, even if
the units are not explicitly written out, to avoid cluttering the expressions.
The cumulants of $Z_{\Delta t,i}$ can be expressed in terms of
$\Gamma_{n}^{\xi_{i}}$. They are found from Eq. 8 by using that the cumulants
of a sum equal the sum of the cumulants:

$\displaystyle\Gamma_{n}^{Z_{\Delta t,i}}=\Delta t\,\Gamma_{n}^{\xi_{i}}.$
(13)

Which of the two representations of the noise (in terms of the cumulants of
$\xi_{i}$ or $Z_{\Delta t,i}$) is used for replacing
$\Gamma_{n}^{p_{i}(t+\tau)-p_{i}(t)}$ is then a matter of convenience that
depends on the application.

###  3.1 Numerical implementation and generation of increments $Z_{\Delta t}$

The swing equations are implemented numerically by Eulers method [24, 25] with
step-size $\delta t=10^{-3}\,\text{s}$. The power input fluctuations
$p_{i}(t)$ are updated accordingly:

$\displaystyle p_{i}(t+\delta t)-p_{i}(t)$
$\displaystyle=-\alpha_{i}p_{i}(t)\cdot\delta t+\sigma_{i}Z_{\delta
t,_{i}}+\text{h.o.}\,$ (14) $\displaystyle Z_{\delta t,i}$
$\displaystyle=\int_{t}^{t+\delta t}\text{d}t^{\prime}\xi_{i}(t^{\prime}).$
(15)

Euler’s method corresponds to dropping the higher order terms (h.o.). At each
time-step the random variables $Z_{\delta t,i}$ should be drawn from
$\text{pr}\big{(}Z_{\delta t,i}\big{)}=\text{pr}\big{(}\int_{t}^{t+\delta
t}\text{d}t^{\prime}\xi_{i}(t^{\prime})\big{)}$. Although
$\Gamma_{n}^{Z_{\delta t,i}}=\delta t\cdot\Gamma_{n}^{\xi_{i}}$ themselves are
determined from $p_{i}(t+\tau)-p_{i}(t)$ according to Eq. 12, numerically
constructing a distribution from its cumulants is in general not tractable for
strongly non-Gaussian distributions [26, 27]. Here we will give our solution
to this problem, and compare it to the alternative approach of [8] in Sec.
3.2. We may expand the cumulants $\Gamma_{n}^{Z_{\Delta t,i}}$ (for $n\geq 2$)
in $\alpha\Delta t$:

$\displaystyle\sigma^{n}\Gamma_{n}^{Z_{\Delta t,i}}$
$\displaystyle=\Gamma_{n}^{p_{i}(t+\Delta t)-p_{i}(t)}\cdot n\alpha\Delta
t\cdot\Big{(}\big{[}\exp(-\alpha\Delta
t)-1\big{]}^{n}+\big{[}1-\exp(-n\alpha\Delta t)\big{]}\Big{)}^{-1}$
$\displaystyle=\Gamma_{n}^{p_{i}(t+\Delta t)-p_{i}(t)}\cdot[1+O([\alpha\Delta
t])].$ (16)

Since $\sigma_{i}^{n}\Gamma_{n}^{Z_{\Delta t,i}}$ are the cumulants of
$\sigma_{i}Z_{\Delta t,i}$, and since $\alpha\tau\approx 10^{-2}$ is small,
setting $\Delta t=\tau$ in Eq. 16 means that the difference of using
$\text{pr}(\text{increment data})$ for constructing either
$\operatorname{\text{pr}}\big{(}p_{i}(t+\tau)-p_{i}(t)\big{)}$ or directly
$\text{pr}(\sigma_{i}Z_{\tau,i})$ is of $O([\alpha\tau])$. Thus, neglecting
the $O(\alpha\tau)$ term:

$\displaystyle\operatorname{\text{pr}}(\sigma_{i}Z_{\tau,i})\approx\operatorname{\text{pr}}\big{(}p_{i}(t+\tau)-p_{i}(t)\big{)}\stackrel{{\scriptstyle\approx}}{{\rightarrow}}\operatorname{\text{pr}}(\text{increment
data})\,,$ (17)

such that the distribution $Z_{\tau,i}$ can simply be approximated by the
distribution of the increments found from the data, rescaled by $\sigma_{i}$.
From here on we will therefore set:

$\displaystyle\operatorname{\text{pr}}(\sigma_{i}Z_{\tau,i})=\operatorname{\text{pr}}(\text{increment
data})\,.$ (18)

What is left is using this distribution of $Z_{\tau,i}$ to construct the
distribution of $Z_{\delta t,i}$ for the finer time resolution $\delta t$.
111Formally, random variables $Z_{\tau}$ are only increments of a (time-
homogeneous) white noise process if the distribution of $Z_{\tau}$ is
infinitely divisible; see [28]. Here we assume this is the case (to a
sufficiently good approximation). If $\operatorname{\text{pr}}(Z_{\tau}$) is
not infinitely divisible, then for $\delta t\rightarrow 0$ the distribution of
$Z_{\delta t}$ calculated by the procedure shown in this section may give
complex probabilities rather than real, positive numbers. We have not found
this to be the case for $\delta t=10^{-3}-10^{-4}\,\text{s}$ and
discretization up to several thousand bins.  Using the fact that the Fourier
transform of the distribution of a sum of independent random variables equals
the product of the Fourier transforms of the distributions of the random
variables [16], we get:

$\displaystyle\mathcal{F}\big{(}\operatorname{\text{pr}}(Z_{\tau})\big{)}[\omega]$
$\displaystyle=\Big{[}\mathcal{F}\big{(}\operatorname{\text{pr}}(Z_{\delta
t})\big{)}[\omega]\Big{]}^{\tau/\delta t}\,.$ (19)

For Gaussian $\text{pr}(Z_{\tau})=N(0,1)$, this is solved by
$\operatorname{\text{pr}}(Z_{\delta t})=N(0,\delta t/\tau)$, in which case
(for $\text{pr}(Z_{\tau})$ given) the standard Euler method for Gaussian SDEs
is obtained [17]. In general, Eq. 19 should be solved for
$\operatorname{\text{pr}}(Z_{\delta t})$ numerically.
Once $\operatorname{\text{pr}}(Z_{\tau})$ is discretized (i.e. represented as
a histogram),
$\mathcal{F}\big{(}\operatorname{\text{pr}}(Z_{\tau})\big{)}[\omega]$ is
easily found by a discrete Fourier transform (for which we used the Fast
Fourier Transform library implemented in Scipy [29]). Similarly, once
$\mathcal{F}\big{(}\operatorname{\text{pr}}(Z_{\delta t})\big{)}[\omega]$ is
found, $\operatorname{\text{pr}}(Z_{\delta t})$ is easily obtained by the
inverse discrete Fourier transform. The complex logarithm is multi-valued,
however, such that:

$\displaystyle\mathcal{F}\big{(}\operatorname{\text{pr}}(Z_{\delta
t})\big{)}[\omega]$ $\displaystyle=\exp\Big{(}\frac{\delta
t}{\tau}\Big{[}\operatorname{\text{Log}}\mathcal{F}\big{(}\operatorname{\text{pr}}(Z_{\tau})\big{)}[\omega]+n(\omega)\cdot
2\pi i\Big{]}\Big{)},\quad n(\omega)\in\mathbb{Z}$ (20)

where $\operatorname{\text{Log}}$ denotes the principal branch of the
logarithm (i.e. with imaginary part in $(-\pi,\pi]$). The Fourier transform
$\mathcal{F}\big{(}\operatorname{\text{pr}}(Z_{\delta t})\big{)}[\omega]$ is
continuous in $\omega$ and equals zero for $\omega=0$, so $n(\omega)$ should
be chosen accordingly. Discretizing $\text{pr}(Z_{\tau})$ with resolution
$\Delta Z_{\tau}$ requires frequencies up to $\omega=\pm\frac{\pi}{\Delta
Z_{\tau}}$; even when discretizing the $\text{pr}(Z_{\tau})$ taken from the
data by up to $10^{4}$ bins,
$\text{Im}\Big{(}\operatorname{\text{Log}}\mathcal{F}\big{(}\operatorname{\text{pr}}(Z_{\tau})\big{)}[\omega]+n(\omega)\cdot
2\pi i\Big{)}\in(-\pi,\pi]$, such that we can simply set $n(\omega)=0$ (see
Fig. 4)
222$\text{Im}\,\log\big{(}\mathcal{F}(\operatorname{\text{pr}}(Z_{\delta
t}))[w]\big{)}$ can be shown to be equal to $\frac{\delta
t}{\tau}\times\big{(}\Gamma^{Z_{\tau}}_{1}\,\omega-\Gamma^{Z_{\tau}}_{3}\,\omega^{3}/3!+\Gamma^{Z_{\tau}}_{5}\,\omega^{5}/5!\dots\big{)}$.
A translation of the distribution $\text{pr}(Z_{\tau})$ along the $x$-axis
changes $\Gamma^{Z_{\tau}}_{1}$ and hence changes the linear term of
$\text{Im}\,\big{(}\mathcal{F}(\operatorname{\text{pr}}(Z_{\delta
t}))[w]\big{)}$. We found that in order for setting $n(\omega)=0$ to be valid
it is important to ensure that the mode of the discretized distribution
$\operatorname{\text{pr}}(Z_{\tau})$ equals $0$, by shifting the $x$-axis if
this is not the case.. Once $\operatorname{\text{pr}}(Z_{\delta t,i})$ is
found, random increments $Z_{\delta t,i}$ can be generated with standard
sampling methods; we used ’rv_histogram’ implemented in Scipy [29]. Examples
of the resulting time-series are shown in Fig. 5 for different distributions.
We term our method, based on Eq. 19, the Fourier-method in the following.

|
---|---

(a)

(b)

Figure 4:
$\text{Im}\,\log\big{(}\mathcal{F}(\operatorname{\text{pr}}(Z_{\tau}))[w]\big{)}$
for the standardized wind data (a) and solar data (b) from [22], discretized
in $10^{4}$ and $3\cdot 10^{3}$ bins respectively. The values remain in
$(-\pi,\pi]$, such that within the considered range of frequencies the
principal branch of the logarithm satisfies the requirement that the Fourier
transform is continuous.

|  |
---|---|---

(a)

(b)

(c)

Figure 5: Example of time-series generated from Eq. 2 by the procedure of Sec.
3.1 with $\alpha=10^{-2}\,\text{s}$, for various (standardized) distributions
of $p(t+\tau)-p(t)$. a) Gaussian distribution. b) Wind power generation [22].
c) Solar power generation during the sunniest hour of the day [22].

###  3.2 Approximation by a compound Poisson process and its small
$\tau$-approximation

An alternative method for reconstructing non-Gaussian noise from data was used
in [8], without explicitly mentioning the involved approximations, by
representing the noise by a compound Poisson process. A continuous time Markov
process, in this case assumed for $p(t)$, can be characterized by the rescaled
moments of its transition distribution, the Kramers-Moyal coefficients
$M^{(n)}(p)\equiv\lim_{\Delta t\rightarrow 0}\frac{1}{\Delta
t}\langle\big{(}p(t+\Delta t)-p(t)\big{)}^{n}|p(t)=p\rangle$ [17, 30]. The
procedure of [8] effectively approximates the Kramers-Moyal coefficients for
$n\geq 2$ (which, for the process of Eq. 2, are independent of $p$) by their
values for finite $\Delta t$, with $\Delta t=\tau$ equal to the data
resolution. The procedure of [8] approximates the correct Kramers-Moyal
coefficients arbitrarily well for data resolution $\tau\rightarrow 0$ (see
[30], or the derivations that follow in this section). Here we will show,
taking Eq. 16 as a starting point, that this compound Poisson (’CP’)
approximation is an approximation in addition to that of Eq. 16 (i.e. the
approximation entering also the Fourier-method), and explicitly derive the
error terms for the CP-approximation. It turns out that for finite $\Delta t$
the Kramers-Moyal coefficients are better approximated by cumulants of the
data than by its moments; the Fourier-Method corresponds to the former, while
the CP approximation corresponds to the latter.
We will find that although the procedure of [8] is justified for non-Gaussian
noise in the limit of infinite measurement resolution $\tau\rightarrow 0$, it
is uncontrolled in the sense that it is not clear $\grave{a}$ priori how small
$\tau$ should be for this mapping to be a good approximation (see also [31]).
For $\tau=1\,\text{s}$, the approximation is relatively accurate for the wind
generation data from [22], while it gives significant errors for the solar
generation data from [22], as we discuss below in Table 2, and fails for
Gaussian fluctuations (for which the authors of [8] do not use the compound
Poisson process), in contrast to our Fourier-method that applies also to
Gaussian noise.
Under mild assumptions [17], the formal time-derivative of any continuous-time
Markov process, and in particular $\xi_{i}(t)$ (entering in $\dot{p}_{i}(t)$
according to Eq. 2), can be decomposed into a sum of Gaussian noise and a
compound Poisson process; see [17, 28, 30] for derivations 333Gaussian noise
itself can be obtained as the limit of a sequence of compound Poisson
processes [16, 28], but this will not be useful for our purposes.. The
Gaussian noise generates continuous Brownian motion of $p_{i}(t)$, and the
compound Poisson process generates discontinuous ’jumps’. A compound Poisson
process consists of a series of delta functions, randomly distributed over
time with the average number of pulses per second given by a rate $\rho$. Each
of the delta pulses has a prefactor $z$ randomly drawn from a size
distribution $\lambda(z)$, which determines the distribution of jump sizes of
$p_{i}(t)$. If the set of possible jump sizes $z$ is discrete, the compound
Poisson process may alternatively be seen as a sum of Poisson processes, each
with its own rate and its own jump size $z$ [8, 16], which makes it
’compound’.

The presence of jumps can be visualized in Fig. 5: For Fig. 5(a) the noise
$\xi_{i}$ is Gaussian, and $p_{i}(t)$ follows a continuous (but non-
differentiable 444That is, the formal derivative $\dot{p}_{i}(t)$ does not
have a finite value for any $t$ [16, 17].) path, while for non-Gaussian noise
(Figs. 5(b-c) this continuous movement is randomly interrupted by sudden
’jumps’ with various jump sizes. Discontinuous jumps in $p_{i}(t)$ correspond
to delta functions in its derivative, i.e. in $\xi_{i}(t)$. If the
contribution of the Brownian motion part of $p_{i}(t)$ is assumed to be
negligible compared to that of the discontinuous jumps, we may approximate
$\xi_{i}(t)$ by a compound Poisson process.
For a given time interval $\Delta t$ (not necessarily equal to the discrete
time-steps in the numerical simulation), the random increments
$Z^{\text{cp}}_{\Delta t}\equiv\int_{t}^{t+\Delta
t}\operatorname{\text{d}\\!}t^{\prime}\xi^{\text{cp}}(t^{\prime})$ of a
compound Poisson process $\xi^{\text{cp}}(t)$ can be generated as follows
[16]:

  * •

The total number of pulses in the time interval $[t,t+\Delta t)$ is generated
as a Poisson random variable $x_{\Delta t}$ with parameter $\rho\Delta t$.

  * •

The contribution of each of the pulses is assigned a random size $z$,
independently drawn from a distribution $\lambda(z)$. The increment
$Z^{\text{cp}}_{\Delta t}$ is then the sum of the sizes $z$ of the
contribution of each of the pulses.

Analytically, instead of subsequently generating the compound Poisson process,
an alternative representation of $Z^{\text{cp}}_{\Delta t}$ is more
convenient: we can discretize $z$ according to $z\in\\{\dots,-\Delta
z,0,\Delta z,\dots\\}$ with $\Delta z\rightarrow 0$ and denote the number of
pulses in $[t,t+\Delta t)$ with contribution smaller than $z$ by $x_{\Delta
t}(z)$. Then the increments can be written as $Z^{\text{cp}}_{\Delta
t}=\sum_{z}z\cdot\Delta x_{\Delta t}(z)$, where $\Delta x_{\Delta t}(z)\equiv
x_{\Delta t}(z+\Delta z)-x_{\Delta t}(z)$ give the number of pulses with size
in $[z,z+\Delta z)$, which are independent Poisson random variables with
parameter $\rho\Delta t\lambda(z)\cdot\Delta z$. If the statistics of the
compound Poisson process $\xi^{\text{cp}}(t)$ should approximate those of
$\xi(t)$, $\rho$ and $\lambda(z)$ should be chosen such that the increments
$Z^{\text{cp}}_{\Delta t}$ have the same distribution as $Z_{\Delta
t,i}\equiv\int_{t}^{t+\Delta
t}\operatorname{\text{d}\\!}t^{\prime}\xi_{i}(t^{\prime})$. The random
increments $Z^{\text{cp}}_{\Delta t}$ have the cumulant generating function:

$\displaystyle\log\langle\exp(J\cdot Z^{\text{cp}}_{\Delta
t})\rangle_{\operatorname{\text{pr}}(Z^{\text{cp}}_{\Delta t})}$
$\displaystyle=\sum_{z}\log\langle\exp(J\cdot z\cdot\Delta x_{\Delta
t}(z))\rangle_{\operatorname{\text{pr}}(\Delta x_{\Delta t}(z))}$
$\displaystyle=\sum_{z}\log\Big{(}\sum_{n=0}^{\infty}\exp(Jzn)\operatorname{\text{pr}}(\Delta
x_{\Delta t}(z)=n)\Big{)}$
$\displaystyle=\sum_{z}\log\Big{(}\exp\big{(}-\rho\Delta t\lambda(z)\Delta
z\big{)}\cdot\sum_{n=0}^{\infty}\frac{\big{[}\exp(Jz)\rho\Delta
t\lambda(z)\Delta z\big{]}^{n}}{n!}\Big{)}$ $\displaystyle=\rho{\Delta
t}\cdot\sum_{z}\lambda(z)(\exp(J\cdot z)-1)\Delta z$
$\displaystyle=\rho{\Delta t}\cdot\big{(}\langle\exp(J\cdot
z)\rangle_{\lambda(z)}-1\big{)}\,,$ (21)

where in the first line we used that the cumulant generating function of a sum
is the sum of the cumulant generating functions of the individual variables
[16, 17], and in the second to fourth lines we evaluated the expectation value
over $\Delta x_{\Delta t}(z)\sim\text{Poiss}\big{(}\rho\Delta
t\lambda(z)\Delta z\big{)}$.
The left-hand side is the cumulant generating function of
$\text{pr}(Z^{\text{cp}}_{\Delta t})$, while the right-hand side is (up to a
constant) proportional to the moment generating function of $\lambda(z)$.
Thus:

$\displaystyle\rho\Delta t\mu^{\lambda}_{n}=\Gamma^{Z^{\text{cp}}_{\Delta
t}}_{n}=\Delta t\Gamma^{\xi^{\text{cp}}}_{n}\,,$ (22)

where $\mu^{\lambda}_{n}$ are the moments of $\lambda(z)$. To construct the
compound Poisson process such that the random variables $Z^{\text{cp}}_{\tau}$
follow approximately the same distribution as $Z_{\tau,i}$, we find from
setting $\Delta t=\tau$ in Eq. 22 that $\rho$ and $\lambda(z)$ should be such
that:

$\displaystyle\rho\mu^{\lambda}_{n}$
$\displaystyle\stackrel{{\scriptstyle\approx}}{{\rightarrow}}\Gamma^{Z_{\tau,i}}_{n}/\tau\,.$
(23)

Now it would seem like we have not come any closer to approximating $\xi$,
since in general finding a distribution $\lambda(z)$ with these moments is not
straightforward [26, 27]. However, for small $\tau$ the cumulants and moments
of $Z_{\tau,i}$ approach each other:

$\displaystyle\langle\exp(J\cdot
Z_{\tau,i})\rangle_{\operatorname{\text{pr}}(Z_{\tau_{i}})}-1\,$
$\displaystyle=\exp\log\big{(}\langle\exp(J\cdot
Z_{\tau,i})\rangle_{\operatorname{\text{pr}}(Z_{\tau,i})}\big{)}-1\,$
$\displaystyle=\log\langle\exp(J\cdot
Z_{\tau,i})\rangle_{\operatorname{\text{pr}}(Z_{\tau,i})}+O\big{(}\big{[}\log\langle\exp(J\cdot
Z_{\tau,i})\rangle_{\operatorname{\text{pr}}(Z_{\tau,i})}\big{]}^{2}\big{)}$
$\displaystyle=\log\langle\exp(J\cdot
Z_{\tau,i})\rangle_{\operatorname{\text{pr}}(Z_{\tau,i})}+O\big{(}\tau^{2})\,,$
(24)

where we expanded the exponential in its Taylor series and used that by Eq. 13
the cumulant generating function $\log\langle\exp(J\cdot
Z_{\tau,i})\rangle_{\operatorname{\text{pr}}(Z_{\tau,i})}=\Gamma^{Z_{\tau_{i}}}_{1}\frac{J}{1!}+\Gamma^{Z_{\tau_{i}}}_{2}\frac{J^{2}}{2!}+\dots\sim\tau$.
If the distribution of $Z_{\tau,i}$ is neither Gaussian nor a delta function,
all of its cumulants are non-zero [17]. Expanding both the moment generating
function (the left-hand side of Eq. 24) and the cumulant generating function
(on the right-hand side of Eq. 24) in $J$, we then find that for all orders
$\mu_{n}^{Z_{\tau,i}}=\Gamma_{n}^{Z_{\tau,i}}[1+O(\tau)]$. Neglecting the
$O\big{(}\tau)$ term, Eq. 23 is then solved by setting $\rho=1/\tau$ and
$\lambda(z)=\operatorname{\text{pr}}(Z_{\tau,i}=z)$. Therefore, consistently,
we may replace $\mu^{\lambda}_{n}$ with $\mu_{n}^{\text{data}}$ read off from
the data. After all, this means that $\Gamma^{\xi^{\text{cp}}}_{n}$ can be
replaced by $\mu_{n}^{\text{data}}$ which we will use in a comparison between
both methods in Sec. 6.

The identification $\lambda(z)=\operatorname{\text{pr}}(Z_{\tau,i}=z)$ may be
interpreted as follows: $\xi_{i}$ can always be composed into the sum of
Gaussian noise and a compound Poisson process with a certain
($\tau$-independent) rate. If for small $\tau$ the increments $Z_{\tau,i}$ are
dominated by the contribution of the compound Poisson process, then in the
limit $\tau\rightarrow 0$ the probability of having two or more pulses in the
time interval $[t,t+\tau)$ becomes negligible, such that the distribution of
pulse contributions $\lambda(z)$ becomes (up to an irrelevant delta function
at $z=0$ corresponding to having no pulses at all) equal to the increment
distribution $\operatorname{\text{pr}}(Z_{\tau,i}=z)$.

The case of $\xi_{i}$ being Gaussian noise is an exception: In Eq. 24 the
leading order contribution to $\mu_{n}^{Z_{\tau,i}}$ is the
$O\big{(}\tau^{2})$ term (all cumulants of order $m>2$ are zero), which can no
longer be neglected. While for the non-Gaussian case the approximation is
valid for $\tau\rightarrow 0$, it is not obvious how small $\tau$ should be
for the approximation to be accurate. 555This is in contrast with the
approximation of Eq. 16, where $\alpha\tau$ is unitless, and the approximation
(of replacing $\text{pr}(\sigma_{i}Z_{\tau,i})$ by
$\text{pr}(p_{i}(t+\tau)-p_{i}(t))$) can thus be seen to be justified whenever
$\tau\ll 1/\alpha$.

##  4 Dimensional Reduction of the equations of motion

The different versions of dimensional reduction of the equations of motion
given in this section are based on the fact that an instability of the
synchronized state is caused by overloaded lines. To apply these reductions,
we first have to identify transmission lines that are endangered by possible
overloads (Sec. 4.1). In Sec. 4.2 we apply the so-called synchronized subgraph
approximation to the swing equations (Eqs.1-2,3), for which desynchronization
times are calculated by stochastic sampling with many iterations. As an
alternative analytical approach, applicable to very long desynchronization
times, we apply the WKB-approach to the differential Chapman-Kolmogorov
equation that is presented in Sec. 4.3 and derived in Appendix B. Inserting
the WKB-ansatz into this equation, we derive in Sec. 4.4 Hamilton’s equations
of motion for the full system, including reductions to 12- and two-dimensional
phase spaces, respectively. For the 12-dimensional case we later make use of
the iterative action minimization method (Sec. 5.2) [32], for the two-
dimensional case we derive the analytical form of the action. The WKB-approach
avoids the need for stochastic sampling in the regime where desynchronization
events are very rare.

###  4.1 Partition of the network: Identifying endangered transmission lines

In the absence of fluctuations $p_{i}(t)$, fixed points of the swing equations
satisfy (besides $v_{i}=0$) for all $i$:

$P_{i}+\sum_{j}K_{ij}\sin(\phi_{j}-\phi_{i})=0\,,$ (25)

Deterministically, the system is in stable operation as long as it resides at
a fixed point which is linearly stable. Linear stability of a fixed point
$\bm{\phi}^{\ast}$ is determined by the eigenvalues of the Laplacian matrix
$L(\bm{\phi}^{\ast})$ [5] with elements:

$L_{ik}(\bm{\phi}^{\ast})\equiv\frac{\partial\sum_{k}K_{ik}\sin(\phi_{k}-\phi_{i})}{\partial\phi_{j}}\biggr{\rvert}_{\bm{\phi}^{\ast}}\;.$
(26)

The Laplacian always has a zero eigenvalue corresponding to a global
translation of the phase angles. Assuming the network has a single connected
component, stability of the fixed point is then determined by the remaining
eigenvalues. The fixed point is stable if its second highest eigenvalue (with
zero being the highest eigenvalue), the Fiedler value (with the corresponding
eigenvector called the Fiedler vector), is negative. The authors of [5] show
that a sufficient (but not necessary) condition for linear stability of the
fixed point is $|\phi^{\ast}_{i}-\phi^{\ast}_{j}|<\pi/2$ for any two nodes $i$
and $j$ connected by a line. In power engineering this condition is used as a
standard criterion for stable synchronized operation [6, 7], and it is this
type of fixed points that we will focus on in this work. Linear stability can
then be interpreted as the lines not being overloaded:

$|\phi_{i}^{\ast}-\phi_{j}^{\ast}|<\pi/2\quad\leftrightarrow\quad|f_{ij}|<K_{ij}\,,$
(27)

where $f_{ij}=K_{ij}\sin(\phi_{i}^{\ast}-\phi_{j}^{\ast})$ is the power flow
from node $i$ to node $j$. If the load on some parts of the grid is increased
to the point where $|\phi^{\ast}_{i}-\phi^{\ast}_{j}|\rightarrow\pi/2$ (i.e.,
$|f_{ij}|\rightarrow K_{ij}$) for some lines $(ij)$, the system undergoes a
saddle-node bifurcation, in which the stable fixed point collides with a
saddle-node (for which $|\phi^{s}_{i}-\phi^{s}_{j}|>\pi/2$) and annihilates
[5]. In [5] it is shown that for any bifurcation due to an overload, the
system splits into two areas (which we will call $A_{1}\equiv\\{i\text{ in
area }1\\}$ and $A_{2}\equiv\\{j\text{ in area }2\\}$), for which all lines
connecting the two areas are overloaded: $|f_{ij}|=K_{ij}$ for all lines
$(ij)$ s.t. $i\in A_{1}$ and $j\in A_{2}$. At the bifurcation, the two areas
desynchronize, while the nodes remain synchronized within each area; in
particular, the Fiedler vector becomes
$(\sqrt{\frac{|A_{2}|/|A_{1}|}{|A_{1}|+|A_{2}|}}:i\in
A_{1},-\sqrt{\frac{|A_{1}|/|A_{2}|}{|A_{1}|+|A_{2}|}}:j\in A_{2})$.

###  4.2 Swing equations in the synchronized subgraph approximation

Although the stable fixed point $\bm{\phi}^{\ast}$ is linearly stable against
fluctuations $\bm{p}(t)$, strong fluctuations may still kick the system out of
the basin of attraction of $\bm{\phi}^{\ast}$. As mentioned in Sec. 4.1, at
overload the network splits into two areas, each of which remains internally
synchronized. A natural reduction is then to model each of the two areas as a
single oscillator, and approximate the rate of desynchronization between the
two areas by the rate of desynchronization between the two reduced
oscillators. Such a reduction was formulated at the level of Hamilton’s
equations (see Sec. 4.3) in [8] and dubbed the ’synchronized subgraph
approximation’. Here instead we formulate this reduction more generally, and
directly at the level of the swing equations: We specify the two areas of the
grid which are at risk for a system split (shown for the Brazilian grid in
Fig. 1), choose area $1$ to be the area in which the total power production
$\sum_{i\in A_{1}}P_{i}=-\sum_{j\in A_{2}}P_{j}$ is positive (where the sets
$A_{1}$ and $A_{2}$ are as in Sec. 4.1); then we write
$\phi_{1}=\phi_{i}-(\phi_{i}^{SN}-\frac{\pi}{4})$ for all $i\in A_{1}$ and
$\phi_{2}=\phi_{j}-(\phi_{j}^{SN}+\frac{\pi}{4})$ for all $j\in A_{2}$ (where
$\bm{\phi}^{SN}$ are the phase angles at bifurcation), and sum the swing
equations over all $i\in A_{1}$ and $j\in A_{2}$, respectively, to obtain:

$\displaystyle M_{1}\dot{v}_{1}+\gamma_{1}v_{1}$
$\displaystyle=p_{1}(t)+P+K\sin(\phi_{2}-\phi_{1})$ (28) $\displaystyle
M_{2}\dot{v}_{2}+\gamma_{2}v_{2}$
$\displaystyle=p_{2}(t)-P+K\sin(\phi_{1}-\phi_{2})$ (29)
$\displaystyle\dot{p}_{1}$ $\displaystyle=-\alpha
p_{1}+\sigma_{1}\cdot\xi_{1}$ (30) $\displaystyle\dot{p}_{2}$
$\displaystyle=-\alpha p_{2}+\sigma_{2}\cdot\xi_{2}\,,$ (31)

where:

$\displaystyle P$ $\displaystyle\equiv\sum_{i\in A_{1}}P_{i}=-\sum_{j\in
A_{1}}P_{j}\qquad K$ $\displaystyle\equiv\sum_{i\in A_{1},j\in
A_{2}}K_{ij}\qquad\gamma_{k}$ $\displaystyle\equiv\sum_{i\in A_{k}}\gamma_{i}$
$\displaystyle M_{k}$ $\displaystyle\equiv\sum_{i\in A_{k}}M_{i}\qquad\qquad$
$\displaystyle\sigma_{k}\xi_{k}\equiv\sum_{i\in
A_{k}}\sigma_{i}\xi_{i}\qquad\sigma_{k}^{2}$ $\displaystyle\equiv\sum_{i\in
A_{k}}\sigma_{i}^{2}$ (32)

with $k\in\\{1,2\\}$. This thus gives a reduction of the swing equations in
which each of the two areas involved in a system split is represented by a
single oscillator. This reduced representation is a generalization of the
reduction of [8] in the sense that it applies to any network topology (as long
as the bifurcation is due to an overload), and is written directly in terms of
the swing equations. Furthermore, our description of the noise terms in terms
of their cumulants is essential for quantifying the noise terms within this
dimensional reduction. The fact that the cumulants of a sum are the sums of
the cumulants gives:

$\displaystyle\Gamma_{n}^{\xi_{k}}$ $\displaystyle=\sum_{i\in
A_{k}}\frac{\sigma^{n}_{i}}{\sigma_{k}^{n}}\Gamma_{n}^{\xi_{i}}\qquad\quad
C^{\xi_{k}}(\sigma_{k}J)=\sum_{i\in A_{k}}C^{\xi_{i}}(\sigma_{i}J)$ (33)

for $k\in\\{1,2\\}$.
Lastly, note that the power flow between the two nodes equals $P$. The
bifurcation therefore occurs at $P=K$. The stable and saddle fixed points are
found to be (up to a global phase shift):

$\displaystyle\phi^{\ast}_{1}$ $\displaystyle=-\phi_{2}^{\ast}=$
$\displaystyle\frac{1}{2}\arcsin(\frac{P}{K})\,,$ (34)
$\displaystyle\phi^{s}_{1}$ $\displaystyle=-\phi_{2}^{s}=$
$\displaystyle\frac{\pi}{2}-\frac{1}{2}\arcsin(\frac{P}{K})\,,$ (35)

respectively.

###  4.3 The WKB approach

To some extent, the problem becomes analytically tractable in the limit in
which typical fluctuations are small compared to the distance to overload
($\sigma\ll K-P$). To this end, we first have to derive an equation for the
time evolution of the probability to find the system in a state with variables
$\bm{\phi},\bm{v},\bm{p}$ at time $t$, subject to the swings equations (Eqs.
1-2) and for specified initial conditions.

####  4.3.1 The differential Chapman-Kolmogorov equation

Under mild assumptions, further specified in [17], for any continuous time
Markov process, here for the swing equations with fluctuations according to
Eq. 2, the time evolution of the probability distribution
$\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)$ to find phases
$\bm{\phi}$, frequencies $\bm{v}$, power input $\bm{p}$, at time $t$ is given
by an integro-differential equation, called the differential Chapman-
Kolmogorov equation [17]. In Appendix B we derive this evolution equation from
Eqs. 1-2 to be given as:

$\displaystyle\frac{\partial\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)}{\partial
t}=\sum_{i}\Big{[}$
$\displaystyle-\frac{\partial}{\partial\phi_{i}}[v_{i}\cdot\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)]$
$\displaystyle-\frac{\partial}{\partial
v_{i}}\big{[}\frac{1}{M_{i}}\big{(}-\gamma_{i}v_{i}+\overline{P}_{i}+p_{i}+\sum_{j}K_{ij}\sin(\phi_{j}-\phi_{i})\big{)}\cdot\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)\big{]}$
$\displaystyle-\frac{\partial}{\partial
p_{i}}\big{[}-\alpha_{i}p_{i}\cdot\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)\big{]}$
$\displaystyle+\int\operatorname{\text{d}\\!}z_{i}w_{i}(z_{i})\cdot\big{(}\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p}-\bm{\mathbbm{1}}_{i}\sigma_{i}z_{i},t)-\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)\big{)}\Big{]}\,,$
(36)

with the transition rates $w_{i}(z_{i})$ (uniquely) specified by their moment-
generating function
$\langle\exp(Jz_{i})\rangle_{w_{i}(z_{i})}=1+C^{\xi_{i}}(J)\equiv
1+\Gamma_{1}^{\xi_{i}}\frac{J}{1!}+\Gamma_{2}^{\xi_{i}}\frac{J^{2}}{2!}+\dots$.
In combination with the WKB ansatz, $w_{i}(z)$ leads to the occurrence of the
cumulant generating function $C^{\xi_{i}}$ in Hamilton’s equation of motion
below. For a data implementation according to the Fourier-method, the
cumulants are given by
$\Gamma_{n}^{\xi_{i}}=\Gamma_{n}^{\text{data}}/(\tau\sigma_{i}^{n})$ (Eqs. 13,
18), and $w_{i}(z_{i})$ can in principle be numerically calculated from its
Fourier transform, see Eq. 84. For Gaussian noise,
$w_{i}(z_{i})=\delta(z_{i})+\frac{1}{2}\frac{\text{d}^{2}\delta(z_{i})}{\text{d}z_{i}^{2}}$
turns the differential Chapman-Kolmogorov equation into a Fokker-Planck
equation. If the power input is generated as a compound Poisson process,
$w_{i}(z_{i})=\rho\lambda_{i}(z_{i})+(1-\rho)\delta(z_{i})$. In general,
$w_{i}(z_{i})$ is a weighted sum of the expression for Gaussian noise and the
expression for the compound Poisson process [17]. For further details see
Appendix B.

####  4.3.2 Hamilton’s equations for the full system

If typical fluctuations are small compared to the distance to overload
($\sigma\ll K-P$), fluctuations that lead to desynchronization become very
rare, and the rate of escape can be found by the WKB-approximation. The WKB-
approximation assumes that escape from the basin of attraction of the stable
fixed point becomes dominated by a single trajectory, along which escape is
overwhelmingly more likely than along any other trajectory. This ’optimal
path’ is found by solving Hamilton’s equations [8, 9, 33, 34] subject to
boundary conditions. Initial values are at the stable fixed point, while the
final values are at the edge of the basin of attraction, in this case we fix
it at a saddle point.

The question remains which saddle point to choose: The system undergoes a
saddle-node bifurcation whenever the lines between any two areas of the grid
are overloaded. Once the lines that are in danger of being overloaded are
specified, we choose as final values the saddle-node which emerges along with
the stable fixed point at this bifurcation. Just like the synchronized
subgraph approximation, the WKB method thus requires the lines which are in
danger of being overloaded to be specified. Here we choose the partition as in
Fig. 1. In general the saddle-node which gives the lowest average time to
desynchronization should be chosen. Following [9, 34] we assume that in this
regime the quasi-stationary distribution follows the WKB-form:

$\displaystyle\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)$
$\displaystyle\propto\exp(-S(\bm{\phi},\bm{v},\bm{p},t))\,,$ (37)

with $S\rightarrow\infty$ as the noise strength goes to zero. Often a large
factor in front of $S$ in Eq. 37 is explicitly split off. The WKB-approach
applies to fluctuations which are rare. For example, when WKB is applied to
large populations in the context of evolutionary game theory [9, 34, 35],
large fluctuations are suppressed if the system size $N$ is large. Typical
fluctuations in population sizes are of $O(\sqrt{N})$, rare ones are of
$O(N)$. In this case a factor $N$ is split off.

For our case, if the strength of noise is parameterized by $\sigma\ll 1$,
typical fluctuations in $\bm{\phi},\bm{v},\bm{p}$ are of $O(\sigma)$, while
the WKB approach applies to rare ones of $O(1)$. For Gaussian noise, the WKB-
approach is formally derived in the limit $1/\sigma^{2}\rightarrow\infty$ in
[36, 37], for which the action is proportional to $1/\sigma^{2}$, which may in
principle be split off from the action. For the non-Gaussian case the precise
dependence of $S$ on $\sigma$ is less obvious. Furthermore, in our case
desynchronization events are rare either because the noise strength is too
small to overcome the distance to the boundary of the basin of attraction of
the synchronized state, or because the load of the system is so small that the
distance to the boundary of the basin of attraction is very large, even
compared to a larger noise strength. Inserting the WKB ansatz into the
differential Chapman-Kolmogorov equation, and defining auxiliary momenta
$\bm{\lambda^{\phi}}\equiv\frac{\partial S}{\partial{\bm{\phi}}}$ ,
$\bm{\lambda^{v}}\equiv\frac{\partial S}{\partial{\bm{v}}}$ ,
$\bm{\lambda^{p}}\equiv\frac{\partial S}{\partial{\bm{p}}}$, gives:

$\displaystyle-\frac{\partial S}{\partial t}=\sum_{i}$
$\displaystyle\Big{[}v_{i}\cdot\lambda^{\phi}_{i}-\alpha_{i}p_{i}\cdot\lambda^{p}_{i}+\frac{1}{M_{i}}\big{(}-\gamma_{i}v_{i}+P_{i}+p_{i}+\sum_{j}K_{ij}\sin(\phi_{j}-\phi_{i})\big{)}\cdot\lambda^{v}_{i}\Big{]}$
(38)
$\displaystyle+\int\text{d}z_{i}w_{i}(z_{i})\cdot\Big{[}\exp\big{(}-S(\bm{\phi},\bm{v},\bm{p}-\mathbbm{1}_{i}\sigma_{i}z_{i},t)+S(\bm{\phi},\bm{v},\bm{p})\big{)}-1\Big{]}+\text{constant}\,,$
(39)

where we are free to set the constant to zero by absorbing it in a
$t$-dependent prefactor in Eq. 37. Expanding the action in the exponential to
first order in $\sigma_{i}$ and using
$\int\text{d}z_{i}w_{i}(z_{i})\exp(\sigma_{i}\lambda_{i}^{p}z_{i})=1+C^{\xi_{i}}(\sigma_{i}\lambda_{i}^{p})$
(Eq. 39) gives the Hamilton-Jacobi equation:

$\displaystyle-\frac{\partial S}{\partial t}=H$ (40)

with the Hamiltonian:

${}H=\sum_{i}\Big{[}v_{i}\lambda^{\phi}_{i}-\alpha_{i}p_{i}\lambda^{p}_{i}+C^{\xi_{i}}(\sigma_{i}\lambda^{p}_{i})+\big{(}-\gamma_{i}v_{i}+P_{i}+p_{i}+\sum_{j}K_{ij}\sin(\phi_{j}-\phi_{i})\big{)}\cdot\frac{\lambda^{v}_{i}}{M_{i}}\Big{]}.$
(41)

Note that in contrast to [8], our Hamiltonian contains the cumulant generating
function $C^{\xi_{i}}$ of the fluctuations $\xi_{i}$. The cumulants are
related to the power increment fluctuations (generated by the Fourier-method)
according to Eq. 11 and quantify non-Gaussian effects: for Gaussian noise all
cumulants $\Gamma_{n}^{\xi_{i}}$ with $n>2$ are zero. The expansion of the
cumulant generating function
$C^{\xi_{i}}(\sigma_{i}\lambda^{p}_{i})=\sigma_{i}\Gamma_{1}^{\xi_{i}}\cdot\frac{\lambda_{i}^{p}}{1!}+\sigma^{2}_{i}\Gamma_{2}^{\xi_{i}}\cdot\frac{(\lambda_{i}^{p})^{2}}{2!}+\dots$
can be alternatively obtained if the Kramers-Moyal expansion [16] is applied
to the differential Chapman-Kolmogorov equation (i.e. the Taylor expansion of
$\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p}-\bm{\mathbbm{1}}_{i}\sigma_{i}z_{i},t)$
in $\sigma_{i}z_{i}$).

The stationary solution $S_{stat}(\bm{\phi},\bm{v},\bm{p})$ to the Hamilton-
Jacobi equation can be equivalently represented in terms of Hamilton’s
equations:

$\displaystyle
S_{stat}(\bm{\phi},\bm{v},\bm{p})=\int_{t_{0}}^{t_{1}}\operatorname{\text{d}\\!}t\sum_{i}\big{[}\lambda^{\phi}_{i}\frac{\operatorname{\text{d}\\!}\phi_{i}}{\operatorname{\text{d}\\!}t}+\lambda^{v}_{i}\frac{\operatorname{\text{d}\\!}v_{i}}{\operatorname{\text{d}\\!}t}+\lambda^{p}_{i}\frac{\operatorname{\text{d}\\!}p_{i}}{\operatorname{\text{d}\\!}t}]\,,$
(42)

with $t_{1}\gg t_{0}$ and where $\bm{\phi}(t)$, $\bm{v}(t)$, $\bm{p}(t)$,
$\bm{\lambda}^{\bm{\phi}}(t)$, $\bm{\lambda}^{\bm{v}}(t)$ and
$\bm{\lambda}^{\bm{p}}(t)$ are solutions of Hamilton’s equations:

$\displaystyle\frac{\operatorname{\text{d}\\!}\phi_{i}}{\operatorname{\text{d}\\!}t}$
$\displaystyle=$ $\displaystyle\frac{\partial H}{\partial\lambda^{\phi}_{i}}$
$\displaystyle=v_{i}\,,$
$\displaystyle\frac{\operatorname{\text{d}\\!}v_{i}}{\operatorname{\text{d}\\!}t}$
$\displaystyle=$ $\displaystyle\frac{\partial H}{\partial\lambda^{v}_{i}}$
$\displaystyle=\Big{(}-\gamma_{i}v_{i}+P_{i}+p_{i}+\sum_{j}K_{ij}\sin(\phi_{j}-\phi_{i})\Big{)}/M_{i}\,,$
$\displaystyle\frac{\operatorname{\text{d}\\!}p_{i}}{\operatorname{\text{d}\\!}t}$
$\displaystyle=$ $\displaystyle\frac{\partial H}{\partial\lambda^{p}_{i}}$
$\displaystyle=-\alpha
p_{i}(t)+\sigma_{i}C^{{}^{\prime}\xi_{i}}(\sigma_{i}\lambda_{i}^{p})\,,$
$\displaystyle\frac{\operatorname{\text{d}\\!}\lambda^{\phi}_{i}}{\operatorname{\text{d}\\!}t}$
$\displaystyle=$ $\displaystyle-\frac{\partial H}{\partial\phi_{i}}$
$\displaystyle=-\sum_{j}\big{(}K_{ij}\cos(\phi_{j}-\phi_{i})\cdot[\lambda^{v}_{j}/M_{j}-\lambda^{v}_{i}/M_{i}]\big{)}\,,$
$\displaystyle\frac{\operatorname{\text{d}\\!}\lambda^{v}_{i}}{\operatorname{\text{d}\\!}t}$
$\displaystyle=$ $\displaystyle-\frac{\partial H}{\partial v_{i}}$
$\displaystyle=-\lambda^{\phi}_{i}+\gamma_{i}\lambda^{v}_{i}/M_{i}\,,$
$\displaystyle\frac{\operatorname{\text{d}\\!}\lambda^{p}_{i}}{\operatorname{\text{d}\\!}t}$
$\displaystyle=$ $\displaystyle-\frac{\partial H}{\partial p_{i}}$
$\displaystyle=\alpha\lambda^{p}_{i}-\lambda^{v}_{i}/M_{i}\,,$ (43)

with all the parameters as in the swing equations. Here
$C^{{}^{\prime}\xi_{i}}$ denotes the derivative of the cumulant generating
function of $\xi_{i}$. While in our Hamilton’s equations the higher order
cumulants naturally occur, it is higher order moments in the approach of [8].
Furthermore note that Hamilton’s equations for $\bm{\phi}$, $\bm{v}$ and
$\bm{p}$ are the same as the swing equations with the noise terms
$\sigma_{i}\xi_{i}$ replaced by
$C^{{}^{\prime}\xi_{i}}(\sigma_{i}\lambda_{i}^{p})$, with the momenta
$\lambda_{i}^{p}(t)$ determined by the remaining equations.
Hamilton’s equations are subject to the boundary conditions that the system is
at the stable fixed point $\bm{\phi}^{\ast}$ (and other variables equal to
$0$) at time $t_{0}$ and at $(\bm{\phi},\bm{v},\bm{p})$ (and other variables
equal to $0$) at time $t_{1}$. The mean time of escape via the saddle
$\bm{\phi}^{s}$ is then estimated as [8, 9, 34]

$\displaystyle\log\langle T\rangle=\text{constant}+S(\bm{\phi}^{s},0,0)\,,$
(44)

with the constant independent of the strength of the fluctuations (though it
does depend on other parameters).
Let us consider the special case of Gaussian noise. For Gaussian noise,
$C_{i}^{\prime}(\lambda_{i}^{p})=\sigma^{2}_{i}\lambda_{i}^{p}$, moreover,
$S\propto 1/(\bm{\sigma}\cdot\bm{\sigma})$ [33, 36, 37]. The proportionality
of the action to the inverse of the noise strength can still be seen from
Hamilton’s equations by rescaling the noise strength by some constant $c$,
i.e. $\sigma^{2}_{i}\rightarrow{\sigma^{\prime}}_{i}^{2}=c\cdot\sigma^{2}_{i}$
at all nodes. For Gaussian noise, Hamilton’s equations are linear in the
momenta. If $\lambda^{\phi}_{i}$, $\lambda^{v}_{i}$, $\lambda^{p}_{i}$ solve
Hamilton’s equations for noise strengths $\sigma^{2}_{i}$, then
${\lambda^{\prime}}^{\phi}_{i}=\lambda^{\phi}_{i}/c$,
${\lambda^{\prime}}^{v}_{i}=\lambda^{v}_{i}/c$,
${\lambda^{\prime}}^{p}_{i}=\lambda^{p}_{i}/c$ (and $\bm{\phi}$, $\bm{v}$,
$\bm{p}$ unchanged) solve Hamilton’s equations for noise strengths
${\sigma^{\prime}}^{2}_{i}$, thereby changing the action:

$\displaystyle S\rightarrow S^{\prime}$
$\displaystyle=\int_{t_{0}}^{t_{1}}\operatorname{\text{d}\\!}t\sum_{i}\big{[}{\lambda^{\prime}}^{\phi}_{i}\frac{\operatorname{\text{d}\\!}\phi_{i}}{\operatorname{\text{d}\\!}t}+{\lambda^{\prime}}^{v}_{i}\frac{\operatorname{\text{d}\\!}v_{i}}{\operatorname{\text{d}\\!}t}+{\lambda^{\prime}}^{p}_{i}\frac{\operatorname{\text{d}\\!}p_{i}}{\operatorname{\text{d}\\!}t}\big{]}$
(45)
$\displaystyle=\frac{1}{c}\int_{t_{0}}^{t_{1}}\operatorname{\text{d}\\!}t\sum_{i}\big{[}\lambda^{\phi}_{i}\frac{\operatorname{\text{d}\\!}\phi_{i}}{\operatorname{\text{d}\\!}t}+\lambda^{v}_{i}\frac{\operatorname{\text{d}\\!}v_{i}}{\operatorname{\text{d}\\!}t}+\lambda^{p}_{i}\frac{\operatorname{\text{d}\\!}p_{i}}{\operatorname{\text{d}\\!}t}\big{]}$
(46) $\displaystyle=S/c\,.$ (47)

Thus, the action is proportional to the inverse of the noise strength.

####  4.3.3 Hamilton’s equations in reduced phase spaces

When the WKB approach is applied to the subgraph approximation as introduced
for the swing equations in Sec. 4.2, Hamilton’s equations of motion are given
by Eqs. 43 with $i\in\\{1,2\\}$, labeling the two disjoint areas in which the
system splits upon an overload. The parameters
$P_{i},\gamma_{i},M_{i},\sigma_{i}$ are given again by Eqs. 32 (with
$P_{1}=-P_{2}=P$ and $K_{12}=K$). The terms representing the presence of noise
are given by $\sum_{i\in
A_{k}}C^{\prime\xi_{i}}(\sigma_{i}\lambda_{k}^{p})=C^{\prime\xi_{k}}(\sigma_{k}\lambda_{k}^{p})$
with $k\in\\{1,2\\}$ for the two areas, where we use again that the sum of
cumulants is the cumulant of sums (Eq. 33). (The same kind of commutativity
does not hold for the moments.) With the three variables
$(\bm{\phi},\bm{v},\bm{p})$, each being two-component vectors, and the
associated auxiliary momenta, the phase space becomes 12-dimensional.
Numerical integration of Hamilton’s equations of motion by means of the
iterative action minimization method (IAMM) [32] turns out to be challenging
and time-consuming due to the existence of multiple time-scales in the
evolution of the momenta. We describe this method in Sec. 5.2.
A further reduction of phase space, disregarding the dynamics of the phase
angles and their time derivatives altogether, is suggested by the fact that
operators of the power grid attribute system splits to lines being overloaded,
rather than due to other transient behaviour of the phase angles [6, 7]. Of
course, in isolation this argument does not justify to drop the phase
dynamics, but we will further justify this approximation in Sec. 5.3. Besides
$P=(P_{1}-P_{2})/2$, the difference in power input in the two regions involved
in the system split is quantified by:

$\displaystyle p(t)\equiv(p_{1}(t)-p_{2}(t))/2\,,$ (48)

with time-evolution:

$\displaystyle\dot{p}$ $\displaystyle=-\alpha p+\sigma\xi\,,$ (49)
$\displaystyle\sigma\xi$
$\displaystyle\equiv(\sigma_{1}\xi_{1}-\sigma_{2}\xi_{2})/2$ (50)
$\displaystyle\sigma^{2}$
$\displaystyle\equiv\frac{\sigma_{1}^{2}+\sigma_{2}^{2}}{4}\,,$ (51)

with areas $1$ and $2$ as in the subgraph approximation of Sec. 4.2. An
overload is hit when $P+p(t)=K$. The associated reduction in Hamilton’s
equations is obtained by dropping the dynamics of
$\bm{\phi},\bm{v},\bm{\lambda^{\phi}},\bm{\lambda^{v}}$:

$\displaystyle\frac{\mathrm{d}p}{\mathrm{d}t}$ $\displaystyle=-\alpha p+\sigma
C^{{}^{\prime}\xi}(\sigma\lambda^{p})\,,$ (52)
$\displaystyle\frac{\mathrm{d}\lambda^{p}}{\mathrm{d}t}$
$\displaystyle=\alpha\lambda^{p}\,,$ (53) $\displaystyle S$
$\displaystyle=\int_{-\infty}^{{t_{o}}}\mathrm{dt}\lambda^{p}\dot{p}\,,$ (54)

where $C^{{}^{\prime}\xi}(\sigma\lambda)$ is the derivative of the cumulant
generating function of $\xi$, and with boundary conditions $p({t_{o}})=K-P$,
and $p(-\infty)=\lambda^{p}(-\infty)=0$, where $t_{o}$ denotes the time at
which the overload occurs. The system of equations is invariant under shifting
of the time axis, and hence the action $S$ will be independent of ${t_{o}}$.
Hamilton’s equations are then solved for $t<{t_{o}}$ (for $t>{t_{o}}$ the
variables are assumed to follow their deterministic trajectories), by:

$\displaystyle p(t)$
$\displaystyle=\int_{-\infty}^{t}\mathrm{d}t^{\prime}\exp(-\alpha(t-t^{\prime}))\sigma
C^{{}^{\prime}\xi}(\sigma\lambda^{p}(t^{\prime}))\,,$ (55)
$\displaystyle\lambda^{p}(t)$
$\displaystyle=\lambda^{p}({t_{o}})\exp(\alpha(t-{t_{o}}))\,,$ (56)

where $\lambda^{p}({t_{o}})$ is found from Eq. 52 by expanding the cumulant
generating function:

$\displaystyle p({t})$
$\displaystyle=\int_{-\infty}^{{t}}\mathrm{d}t^{\prime}\exp\big{(}-\alpha(t-t^{\prime})\big{)}\sum_{n=2}^{\infty}\sigma^{n}\Gamma^{\xi}_{n}\frac{\lambda^{p}(t^{\prime})^{n-1}}{(n-1)!}$
$\displaystyle=\sum_{n=2}^{\infty}\sigma^{n}\Gamma^{\xi}_{n}\int_{-\infty}^{{t}}\mathrm{d}t^{\prime}\exp\big{(}\alpha(t^{\prime}-t)\big{)}\frac{\lambda^{p}({t})^{n-1}\exp\big{(}\alpha\cdot(n-1)\cdot(t^{\prime}-t)\big{)}}{(n-1)!}$
$\displaystyle=\sum_{n=2}^{\infty}\sigma^{n}\Gamma^{\xi}_{n}\frac{\lambda^{p}({t})^{n-1}}{\alpha
n!}\,.$ (57)

Setting $p({t_{o}})=K-P$, this equation can be solved (numerically) for
$\lambda^{p}({t_{o}})$. Integrating the action of Eq. 54 term by term (where
it is convenient to substitute $\lambda^{p}$ as integration variable and use
Eqs. 52-53 and 57), for Gaussian noise the result is:

$\displaystyle S$
$\displaystyle=\alpha\Big{[}\frac{K-P}{\sigma}\Big{]}^{2}\,,$ (58)

while for non-Gaussian noise the action (containing Eq. 58 as a special case)
gives

$\displaystyle S=\sum_{n=2}^{\infty}\frac{[\lambda^{p}({t_{o}})]^{n}}{\alpha
n!}\frac{n-1}{n}\sigma^{n}\Gamma^{\xi}_{n}.$ (59)

Thus, for the two-dimensional phase space $(p,\lambda^{p})$ we have an
analytical expression for the action up to the fact that Eq. 57 is solved
numerically.

####  4.3.4 Expansion of the action in the distance to overload

The action of Eq. 59 can be expanded in $\frac{K-P}{\sigma}$ by inserting
$\lambda^{p}(t_{o})=c_{1}\Big{[}\frac{K-P}{\sigma}\Big{]}+c_{2}\Big{[}\frac{K-P}{\sigma}\Big{]}^{2}+\dots$
(where constants $c_{1},c_{2},\dots$ are to be determined from Eq. 57). The
lowest two orders are:

$\displaystyle
S=\frac{1}{\alpha}\Big{[}\Big{(}\alpha\frac{K-P}{\sigma}\Big{)}^{2}-\frac{4}{9}\Gamma^{\xi}_{3}\Big{(}\alpha\frac{K-P}{\sigma}\Big{)}^{3}+\text{h.o.}\Big{]}\,,$
(60)

with in general the lowest order contribution of $\Gamma^{\xi}_{n}$ being at
order $\frac{1}{\alpha}\Big{(}\alpha\frac{K-P}{\sigma}\Big{)}^{n}$, such that
when desynchronization events become very rare (large $\frac{K-P}{\sigma}$)
the average desynchronization time $\langle T\rangle$ becomes increasingly
sensitive to higher order cumulants (which require increasingly large amounts
of data to accurately estimate them [38]), corresponding to non-Gaussian
effects. This explains why the authors of [3], studying the regime where
desynchronization are common, found that for their systems non-Gaussianity had
only a small effect on desynchronization times. A large correlation time
(small $\alpha$) suppresses the higher order terms, so that
$\alpha\frac{K-P}{\sigma}$ is a good expansion parameter. Eq. 60 can be used
to determine how the contribution of the skewness $\Gamma^{\xi}_{3}$ affects
the desynchronization time. In terms of the stochastic fluctuations
$\\{p_{i}(t)\\}$ entering in the full network, $\Gamma^{\xi}_{3}$ is found
from Eqs. 32, 50 and 16 as:

$\displaystyle\Gamma^{\xi}_{3}=\frac{1}{2^{3}\sigma^{3}\tau}\Big{[}\sum_{i\in
A_{1}}\Gamma_{3}^{p_{i}(t+\tau)-p_{i}(t)}-\sum_{j\in
A_{2}}\Gamma_{3}^{p_{j}(t+\tau)-p_{j}(t)}\Big{]}\cdot[1+O(\alpha\tau)]\,,$
(61)

where $A_{1}$ is the area with positive power production $\sum_{i\in
A_{1}}P_{i}$ and $A_{2}$ is the area with negative power production
$\sum_{j\in A_{2}}P_{j}=-\sum_{i\in A_{1}}P_{i}$. Eqs. 60-61 therefore show
that whether a third cumulant (skewness) in the noise input increases or
decreases the time to desynchronization depends on where in the network it
enters: if it enters in the area with positive power input (here area $A_{1}$)
a positive (negative) skewness decreases (increases) the average time to
desynchronization, while if it enters in the area with negative power input
(here area $A_{2}$) a positive (negative) skewness increases (decreases) the
average time to desynchronization.

The authors of [8], assuming the parameters are spatially homogeneous, made an
expansion similar to Eq. 60. Their expansion is formally in $\frac{K-P}{P}$,
parameterizing the distance to the bifurcation point, and the main difference
in the assumptions on their network parameters is that they set the
fluctuation timescale to $1/\alpha=1\,\text{s}$ (while here we set it, in
accordance with [11], two orders of magnitude larger) 666Formally, for
spatially homogeneous systems their expansion is valid in the limit where
$\frac{K-P}{P}$ is much smaller than all other relevant quantities, but for
the realistic values we consider here we have
$\frac{M\alpha^{2}}{|P|}\ll\frac{K-P}{P}$ and
$\frac{\gamma\alpha}{|P|}\ll\frac{K-P}{P}$, invalidating their expansion in
our case.. Thus, our approximations may not be accurate for the cases studied
in [8] and vice-versa. Nevertheless, the expressions of [8] (for spatially
homogeneous parameters) can be directly compared to the expansion of Eq. 60 if
we consider both as an expansion in $K$ around $K=P$, that is, an expansion in
the distance to the bifurcation point, as summarized in Table 2.
Qualitatively, their predicted influence of the third cumulant is the same as
we found above.

However, there are also significant differences in the two expansions; in
particular, the lowest order contribution to the action found in [8] is at
order $(K-P)^{3/2}$, and is proportional to $\gamma$ and $\alpha^{2}$. In
contrast, the lowest order contribution to the action of Eq. 59 is at order
$(K-P)^{2}$, is proportional to $\alpha$ and completely independent of
$\gamma$. For Eq. 59 the order $(K-P)^{2}$ is exact for the Gaussian case,
while in the results of [8] all orders are relevant also for the Gaussian
case. Lastly, in Eq. 59 the lowest order at which the $n$-th cumulant enters
the action is $O\big{(}(K-P)^{n}\big{)}$, while for [8] this is at
$O\big{(}(K-P)^{n-1/2}\big{)}$. These differences illustrate the sensitivity
to different analytical approaches and the importance of the correlation time
$1/\alpha$ of the stochastic fluctuations (which a simple white noise process
for $p$, corresponding to $1/\alpha\rightarrow 0$, would not include). Further
work should therefore extend this investigation of the interplay between non-
Gaussianity and time-correlations also to more realistic power spectra [3].

| Eq. 59 | [8]
---|---|---
| Lowest order contribution
---
$\frac{\alpha}{\sigma^{2}}(K-P)^{2}$ | $\propto\frac{\gamma\alpha^{2}}{\sqrt{P}\sigma^{2}}(K-P)^{3/2}$
| Lowest order at which
---
$\Gamma^{\xi}_{n}$ contributes
$O\Big{(}\Big{[}\alpha\frac{K-P}{\sigma}\Big{]}^{n}\Big{)}$ | $O\Big{(}\Big{[}\frac{K-P}{P}\Big{]}^{n-1/2}\Big{)}$
| Lowest order impact
---
of positive skewness
| Increase of $\langle T\rangle$ in area
---
with positive generation
and vice-versa
| Increase of $\langle T\rangle$ in the area
---
with positive generation
and vice-versa
Table 2: Comparing the action of Eq. 59 to the expressions of [8].

##  5 Methods

###  5.1 Sampling of the swing equations with the Fourier-method

The pedestrian way to calculate the average desynchronization time is to
numerically integrate the swing equations Eqs. 1 with fluctuation dynamics Eq.
2, choosing $Z_{\Delta t,i}$ of Eq. 14 with probability $\text{pr}(Z_{\Delta
t,i})$, determined by the Fourier-method, that is, calculating
$\text{pr}(Z_{\Delta t,i})$ as inverse Fourier transformation of Eq. 20. For
the subgraph approximation of the swing equations, the sampling refers to the
same method, but applied to Eqs. 28-32.

Here, we calculate the average desynchronization time, for given parameters,
by calculating an average over $6000$ samples. The CPU time necessary for this
computation is roughly proportional to the average desynchronization time. If
the average desynchronization time is $\sim 100$ minutes, calculating the
average desynchronization time for the given set of parameters (i.e. for each
of the points shown in the figures) takes $\approx 200$ hours of CPU time on
an Intel E5-2640 processor. Doing the same computation for an average
desynchronization time of $\sim 1$ year would then cost $\approx 10^{6}$ hours
of CPU time, per set of parameters, which would clearly be infeasible. To
capture really rare blackouts, one should resort to the WKB-approach.
Nevertheless, in a range where both methods apply and are feasible, the WKB
approach should be supplemented by sampling to validate the results of the
various approximations within the WKB approach, to calculate the constant in
$\log\langle T\rangle$ (which is not predicted by WKB), to find how rare
events should be for the WKB approach to provide an accurate approximation,
and, if of interest, to capture also the regime where desynchronization events
are not rare, where WKB does not apply.

###  5.2 The iterative action minimization method

To solve Hamilton’s equations of motion for the synchronized subgraph
approximation, we have to deal with the 12-dimensional phase space. Here we
use the iterative action minimization method [32] to numerically find the
solutions. The method requires a discretization of Hamilton’s equations. We
take a time interval $[T_{0},T_{1}]$, and discretize time into $n=200$ points
along a grid $\\{t_{i}:i=-n/2+1,-n/2+2,\dots,n/2\\}$ with:

$t_{i}=T_{0}+(T_{1}-T_{0})\frac{\sum_{j=1}^{i}[\cosh(a\cdot
j)^{2}/(\cosh(a\cdot j)^{2}+b-1)]}{\sum_{j=1}^{n}[\cosh(a\cdot
j)^{2}/(\cosh(a\cdot j)^{2}+b-1)]}\,.$ (62)

For our concrete case we used $T_{1}=-T_{0}=1075$, $a=0.055$, $b=10$.
At grid points $i=5,\dots,n-4$, we approximate the derivatives that occur in
Hamilton’s equations by $8$-th order central differences [39]:

$\displaystyle\frac{\mathrm{d}f(t)}{\mathrm{d}t}\biggr{\rvert}_{t_{i}}$
$\displaystyle\approx\sum_{j=i-4}^{i+4}c_{ij}f(t_{j})\,,$ (63) $\displaystyle
c_{ij}$ $\displaystyle\equiv\begin{cases}\sum_{l\in(i-4,\dots,i+4)\setminus
i}\frac{1}{t_{i}-t_{l}}&\text{if}\quad j=i\ \,,\vspace{0.5cm}\\\
\frac{\prod_{k\in(i-4,\dots,i+4)\setminus
i,j}(t_{i}-t_{k})}{\prod_{k\in(i-4,\dots,i+4)\setminus
j}(t_{j}-t_{k})}&\text{else}\,.\end{cases}$ (64)

For $i=1,2,3,4,n-3,n-2,n-1,n$ we replace Hamilton’s equations by the
constraints that all the variables equal the boundary conditions, e.g.,
$\bm{\phi}(t_{1})=\bm{\phi^{\ast}}$ and $\bm{\phi}(t_{n})=\bm{\phi^{s}}$ for
the phase angles. Here $\bm{\phi^{\ast}}$ and $\bm{\phi^{s}}$ are only
determined up to a global phase shift, which we set according to Eqs. 34-35.
Together, this discretizes Hamilton’s equations into a set of $2\cdot 6\cdot
n$ equations, which can be solved with an off-the-shelf solver. We used the
Levenberg-Marquardt algorithm implemented in Scipy [40, 41]. The method
requires an initial guess; we find that for the parameters of the Brazilian
network of Sec. 2.3 and $1/\alpha=100\text{s}$ the algorithm does not converge
without a very accurate guess.

To overcome this, we first artificially multiply all $\gamma_{i}$ values by a
factor thousand. As a guess we then use the near-bifurcation solution given in
[8] (this guess needs spatially homogeneous parameters, for which we insert
the average values of the actual parameters), and use the solution as a new
guess for slightly lower values of $\gamma$, until the actual values are
reached. As a further slight improvement to the accuracy of the solution, we
afterwards increase the discretization to $n=800$, using the solution for
$n=200$ (linearly interpolated to the extra gridpoints) as an initial guess.
Decreasing the damping $\gamma$ too quickly leads to oscillations in the
numerical solution for the momenta (and eventually for all the variables). If
the values are lowered too fast, oscillations in the solution will increase
and prevent convergence. An example is shown in Fig. 6. The numerical solution
shows that, as a function of $t$, momentum builds up over time, and then
-after the overload is hit- very quickly drops to $0$. The difference in these
two timescales creates problems with the discretization, as we would need to
discretize time in a very large number of steps to capture both scales
accurately. This leads to oscillations in the numerical solution of the
momentum near the overload. The damping should be decreased in small ($\approx
50$) steps, and Hamilton’s equations solved numerically take up to several
hours for low damping each time. This makes this approach infeasible for
exploring a large space of parameters. Nevertheless, as we will now see, the
optimal path that is obtained with this approach is still useful as a
justification for the overload approximation of Sec. 4.3.3.

(a)

(b)

Figure 6: Momentum $\lambda^{p}$ found numerically by the IAMM (zoom-in in b),
compared to that predicted by the overload approximation. The numerical
solution shows oscillations near the overload, presumably due to the rapid
drop to zero after the overload. Here $t_{o}$ denotes the time at which the
phase difference reaches $\frac{\pi}{2}$ (Hamilton’s equations are invariant
under a time translation and hence $t_{o}$ is not fixed by the equations).

###  5.3 Justification for the overload approximation

The overload approximation amounts to a restriction of the dynamics to that of
the power fluctuations. Combined with the WKB-method this led to Hamilton’s
equations of motion in two-dimensional phase space according to Eqs. 52-53.
Looking purely at the mathematical expressions, what would suggest this
reduction? If we start with the subgraph approximation and subtract Eq. 29
from Eq. 28, we obtain:

$\displaystyle(M_{1}\dot{v}_{1}-M_{2}\dot{v}_{2})/2+(\gamma_{1}v_{1}-\gamma_{2}v_{2})/2$
$\displaystyle=P+p(t)+K\sin(\phi_{2}-\phi_{1})\,,$ (65)

with $p(t)\equiv(p_{1}(t)-p_{2}(t))/2$. The scaling of Eq. 6 shows that
$\dot{v}\sim\alpha^{2}$ and $v\sim\alpha$; more precisely, since $M_{i}$ and
$\gamma_{i}$ have units of $\text{time}^{2}$ and time, respectively, only
their values relative to the timescale $1/\alpha$ matter. The fact that here
$\frac{M_{i}\alpha^{2}}{P}\sim 10^{-6}$ and $\frac{\gamma_{i}\alpha}{P}\sim
10^{-3}$ (Sec. 2.3) motivates the assumption that the left-hand side of Eq. 65
is negligible along the optimal path in the WKB-method. This amounts to a time
scale separation between the phase dynamics and the correlation time between
power increment fluctuations. In the WKB-approach this means that the dynamics
gets independent of inertia and damping. The result for inertia is in
agreement with [11] and [42], according to which inertia has little impact on
the voltage angle fluctuations for long correlation times.

The phase difference $\phi_{2}-\phi_{1}$ then simply follows an equation that
determines the fixed point of the swing equations, satisfying
$p(t)+P+K\sin(\phi_{2}-\phi_{1})=0$. For $P+p(t)<K$ (i.e. no overload) there
are the two known fixed points (as in Eqs. 34-35), the stable fixed point and
the saddle-node. The stable fixed point itself agrees with the solution of the
static power flow equations, here not further considered. 777The fixed point
solution of Eq.1 is obtained for the active power and voltage phases. It
corresponds to the solution of the steady load flow if the reactive power and
the dynamics of the voltage magnitude are neglected. On the engineering side,
much effort was invested in determining a probabilistic load flow [43]. In
contrast to these approaches, we are not interested in the full solution to
the steady load flow under diverse sources of fluctuations, but only in the
fluctuating power flowing between the two areas connected by the overloaded
lines: In the WKB-approach, the probability for a desynchronization event is
then easily calculated with much less effort than required in the
probabilistic load flow approaches. For $P+p(t)\rightarrow K$ the fixed points
merge into one. If the line gets overloaded at some $t=t_{o}$, the system
resides at the stable fixed point for all $t<t_{o}$ and resides at the saddle
fixed point for $t>t_{o}$ (with the two fixed points merged at $t=t_{o}$),
then the phase difference follows the trajectory:

$\displaystyle\phi_{1}(t)-\phi_{2}(t)$
$\displaystyle=\begin{cases}\arcsin(\frac{P+p(t)}{K})\quad&\text{for}\quad
t<{t_{o}}\\\ \pi/{2}\quad&\text{for}\quad t={t_{o}}\\\
\pi-\arcsin(\frac{P+p(t)}{K})\quad&\text{for}\quad t>{t_{o}}\,.\end{cases}$
(66)

Fig. 7 shows that the optimal path computed numerically from Hamilton’s
equations for the synchronized subgraph approximation (Sec. 5.2) is
indistinguishable from the path predicted by Eq. 66, such that indeed
neglecting the left-hand side of Eq. 65 is justified at the level of the
action.

(a)  Figure 7: The phase difference between the two areas involved in the
system split (here reduced to two oscillators), along the optimal path to
desynchronization (obtained by a numerical solution of Hamilton’s equations,
Sec. 4.3.3) for the Brazilian power network (Sec. 2.3) with Gaussian noise
input and distance to overload $\frac{K-P}{K}=9.1\%$ (the path is independent
of $\sigma$), as compared to the fixed point of the swing equations with power
input $P+p(t)$ (Eq. 66). Here $t_{o}$ denotes the time at which the phase
difference reaches $\frac{\pi}{2}$ (Hamilton’s equations are invariant under a
time translation and hence $t_{o}$ is not fixed by the equations).

The authors of [11, 44] show that, for typical fluctuations, dropping the
phase dynamics is always valid in the regime of long correlation time. In
contrast, we argue that for rare events this time-scale separation, here
leading to the overload approximation, should only be seen as an approximation
at the level of the action (Eqs. 44, 59) and the optimal path (Fig. 7) rather
than of the full desynchronization time (the full desynchronization time would
include the prefactor of the exponential).888This is unlike the synchronized
subgraph approximation (whose partition of the grid is suggested by the values
of the slowest mode of the system, the Fiedler mode), which approximates both
small and large fluctuations. Both approximations turn out to also apply away
from the bifurcation point. This means it gives the order of magnitude of the
average time to desynchronization only up to exponential accuracy, even for
very large correlation time of the power input.

This can be seen as follows. If we consider the dynamics of $p$ according to
Eq. 49 in isolation, for Gaussian noise $p$ follows an Ornstein-Uhlenbeck
process [45]. The results of [45] then give that in the limit of large
$(K-P)/\sigma$ the expected hitting time for $p\rightarrow K-P$ follows:

$\displaystyle\langle
t_{o}\rangle=\frac{\sqrt{\pi}}{\alpha^{3/2}}\frac{\sigma}{(K-P)}\exp\Big{(}\alpha\Big{[}\frac{K-P}{\sigma}\Big{]}^{2}\Big{)}\,,$
(67)

where the exponent agrees with Eq. 58. This suggests the approximation
$\langle t_{o}\rangle\approx\langle T\rangle$ for the expected
desynchronization time. However, in escape problems with Gaussian noise the
prefactor (corresponding to the constant in Eq. 44) is known to depend
sensitively on the dynamics near the boundary of the basin of attraction [46].
For the swing equations, the dynamics of $p$ alone cannot be expected to give
a good approximation to the full system near the overload: For
$p(t)\rightarrow(K-P)$, Eq. 66 predicts $(v_{1}-v_{2})\rightarrow\infty$, such
that close to overload neglecting the left-hand side of Eq. 65 is not
justified. The prefactor in Eq. 67 is valid for an isolated Ornstein-Uhlenbeck
process, but as a prediction for the desynchronization time of the swing
equations it should be taken with care; in Sec. 6 we show that this analytical
formula for the overload hitting time does not give an accurate approximation
for the prefactor of the desynchronization time $\langle T\rangle$, when
compared to the values of $\langle T\rangle$ calculated with the sampling of
the full swing equations. Still, it is in the right ballpark.

##  6 Results

###  6.1 Comparison of the Fourier method and the compound Poisson process for
data implementation

In the following we assume that the approximation used for the Fourier method
according to Eq. 19 gives an accurate approximation of the data. For the
construction of a compound Poisson process, the relative error on the
cumulants of the increments induced by setting $\rho=1/\tau$ and
$\lambda(z)=\operatorname{\text{pr}}(Z_{\tau,i}=z)$ is then quantified by
$\Gamma_{n}^{Z^{\text{cp}}_{\tau}}/\Gamma_{n}^{Z_{\tau,i}}-1=\mu_{n}^{\lambda}/\Gamma_{n}^{Z_{\tau,i}}-1=\mu_{n}^{Z_{\tau,i}}/\Gamma_{n}^{Z_{\tau,i}}-1=\mu_{n}^{\text{data}}/\Gamma_{n}^{\text{data}}-1$,
where the equality signs follow from Eq. 22,
$\lambda(z)=\operatorname{\text{pr}}(Z_{\tau,i}=z)$, and Eq. 18, respectively,
which can be evaluated from the data on the increment distribution. Just to
recapitulate: $\Gamma_{n}^{Z^{\text{cp}}_{\tau}}$ refer to the n-th cumulants
of the power increments of the compound Poisson process,
$\Gamma_{n}^{Z_{\tau,i}}$ to the cumulants of the power increments over the
interval $\tau$, identified with those of the data $\Gamma_{n}^{\text{data}}$,
$\mu_{n}^{\lambda}$ to the n-th moments of $\lambda_{i}(z)=w_{i}(z)$,
$\mu_{n}^{Z_{\tau,i}}$ and $\mu_{n}^{\text{data}}$ accordingly. This means
that we can take the difference between moments and cumulants of the data as
indicating the error in the reproduction of histograms by a compound Poisson
process, if we assume Eq. 18 gives an accurate approximation of the data. The
errors on the first few cumulant ratios are shown in Table 3.

Table 3: The relative error $\Gamma_{n}^{Z^{\text{cp}}_{\tau}}/\Gamma_{n}^{Z_{\tau}}-1$ between cumulants of the increments $Z_{\tau}^{\text{cp}}$ and $Z_{\tau}$, induced by setting $\rho=1/\tau$ and $\lambda(z)=\operatorname{\text{pr}}(Z_{\tau}=z)$ in the compound Poisson process. Shown for various distributions of data for power increments, calculated as $\Gamma_{n}^{Z^{\text{cp}}_{\tau}}/\Gamma_{n}^{Z_{\tau}}-1=\mu_{n}^{\text{data}}/\Gamma_{n}^{\text{data}}-1$ (see Sec. 3.2). In particular, while the approximation is relatively accurate for the wind power generation data, it fails for Gaussian distributions, and gives significant errors on higher-order cumulants (e.g. $\sim 30\%$ for $n=12$) for the solar power generation data.  Distribution | $n=2$ | $4$ | $6$ | $8$ | $10$ | $12$
---|---|---|---|---|---|---
Gaussian (zero-mean) | $0.$ | $\infty$ | $\infty$ | $\infty$ | $\infty$ | $\infty$
Wind | $0.000$ | $0.021$ | $0.021$ | $0.025$ | $0.034$ | $0.053$
Solar | $0.000$ | $0.024$ | $0.049$ | $0.085$ | $0.152$ | $0.302$

The impact of these errors in the cumulants on the average desynchronization
time are shown in Fig. 8.

(a)

(b)

Figure 8: (a) Impact of the data implementation on the average
desynchronization time according to the Fourier (compound Poisson) method, red
(blue) color, respectively, for three values of $1/\alpha$: from top to
bottom, $1/\alpha=10^{2}\,\text{s},10^{1}\,\text{s},10^{0}\,\text{s}$
(crosses, squares, circles). The Gaussian case (green) is reached for
($1/\alpha\rightarrow\infty$). (b) Histograms of standardized increments
$Z_{\tau}$ for $5\cdot 10^{6}$ samples, $10^{3}$ bins (where the Fourier-
method samples from the exact distribution $\text{pr}(Z_{\tau}=z)$), colors as
in (a). For further explanations see the text.

As the difference is only clearly seen for very large times, we determine it
within the overload approximation, for which the WKB approximation becomes
fully analytically tractable. The implementation is either according to Sec.
3.1 (the Fourier-method), or by a compound Poisson process as in Sec. 3.2
(’CP’). We have chosen the noise close to a Gaussian distribution, as there
the difference between the Fourier and compound Poisson approaches is largest:

$\text{pr}(Z_{\tau}=z)=0.9\cdot
N[0,1](z)+0.05\cdot\sqrt{2}\text{exp}(-\sqrt{2}|z|)\,,$ (68)

where $N$ denotes the normal distribution, with the parameters given in square
brackets. In the overload approximation, $\log\langle T\rangle$ is predicted
only up to a constant, which for comparison we here set equal to $5.05$ (which
matches the results of sampling for Gaussian noise and
$1/\alpha=100\;\text{s}$) for all cases. Recall that within the overload
approximation noise enters via the cumulants $\Gamma^{\xi}_{n}$ in Eq. 59. In
the Fourier-method $\Gamma^{\xi}_{n}$ is approximated according to Eq. 13 with
$\Delta t=\tau$ and $\Gamma_{n}^{Z_{\tau,i}}$ identified with the cumulants of
the data, while in the CP-method the cumulants are approximated by Eq. 22 with
$\Delta t=\tau$ and $\Gamma_{n}^{Z_{\tau}^{cp}}$ approximated by the moments
of the data.

For the correlation time $1/\alpha\rightarrow\infty$, the difference to the
results for a Gaussian distribution disappears, while the difference gets
larger for small $1/\alpha$. In Fig. 8 we choose from top to bottom,
$1/\alpha=10^{2}\,\text{s},10^{1}\,\text{s},10^{0}\,\text{s}$ (crosses,
squares, circles). Note that for large desynchronization times the difference
between the two ways of including the noise gets larger. For example, the
Fourier-method predicts the average desynchronization time to be $3.55$ days
for $\Big{(}\frac{K-P}{\sigma}\Big{)}^{2}=760\,\text{s}$ and $17.6$ years for
$\Big{(}\frac{K-P}{\sigma}\Big{)}^{2}=1515\,\text{s}$, both for
$1/\alpha=100\;\text{s}$ (where differences are smallest), while the compound
Poisson approximation predicts an average desynchronization time of $3.10$
days and $10.6$ years, respectively. If we compare the increment data
$Z_{\tau}$, artificially generated according to the Fourier-method or to the
compound Poisson process, Fig. 8(b) shows the different histograms of
standardized increments $Z_{\tau}$ for $5\cdot 10^{6}$ samples, $10^{3}$ bins
(where the Fourier-method samples from the exact distribution
$\text{pr}(Z_{\tau}=z)$ without the need of performing a Fourier transform).

###  6.2 Comparison of the average desynchronization time for the different
dimensional reductions

In Figs. 9 (for Gaussian fluctuations)

(a)

(b)

Figure 9: Average desynchronization time $\langle T\rangle$ (in seconds) for
Gaussian fluctuations on a model of the Brazilian power network (2.3)
according to Fig. 1(a), calculated either by direct simulation of the swing
equations ($\delta t=10^{-3}\,\text{s}$, $\langle T\rangle$ averaged over 6000
samples) for the full network and the reduced network, by the WKB method (Sec.
4.3.3) for the reduced network or by the overload prediction of Sec. 4.3.3.
(a) for constant distance to overload $(K-P)/K=9.1\%$ and varying $\sigma$,
(b) for varying distances to overload $(K-P)/K=5.9-22.8\%$ and constant
$\sigma$. For further explanations see the text.

and 10 (for wind data) we compare the average time to desynchronization
between sampling of the full set of swing equations, sampling according to
their subgraph approximation, the WKB-overload approximation, and, in Fig.
9(a), the WKB-approximation in the subgraph approximation (12-dimensional
phase space). Within the overload approximation of Sec. 4.3.3, the average
desynchronization time only depends on $(K-P)/\sigma$ rather than on the
individual quantities, as can be found by dividing Eq. 49 by $\sigma$ or by
considering the action Eq. 58. Panels (a) and (b) in both figures differ in
how desynchronization events become rare: (a) for constant distance to
overload $(K-P)/K=9.1\%$ and varying $\sigma$ towards smaller values, (b) for
increasing distances to overload $(K-P)/K=5.9-22.8\%$ and constant $\sigma$.

(a)

(b)

Figure 10: Average desynchronization time $\langle T\rangle$ (in seconds) for
fluctuations according to (rescaled) wind data on a model of the Brazilian
power network (2.3) according to Fig.1(a), calculated either by direct
simulation of the swing equations for the full network and reduced network or
by the overload prediction of Sec. 4.3.3 (Eqs. 57-59 truncated after $8$
cumulants). Otherwise the same as Fig. 9.

WKB predicts $\log\langle T\rangle$ only up to a constant (in general
independent of $\sigma$, for the overload approximation also independent of
$K-P$), which is here chosen to match the results by direct simulation. In all
four panels we use the realization of the Brazilian grid according to Fig.
1(a) and see very good agreement for $(K-P)\gg\sigma$. Sampling within the
synchronized subgraph approximation is furthermore also an excellent
approximation for events which are not rare, i.e. where the WKB approximation
does not apply, and of the constant in $\log\langle T\rangle$ which is not
predicted by the WKB approximation. Here $K$ is the conductance and $P$ the
power flow through the most heavily loaded line, $K-P$ is the distance to
overload, $\sigma$ is the noise strength (as in Eq. 51), here chosen to be
equal for both nodes at which stochastic fluctuations enter.

In Fig. 11 we compare the average time to desynchronization for the sampling
of the full and reduced set of swing equations (with Gaussian noise) with the
full analytical solution for the overload approximation according to Eq. 67
(i.e. including the prefactor), with the latter giving a somewhat inaccurate
prediction. In contrast to the reduction by the synchronized subgraph
approximation, the overload approximation is thus only valid at the level of
the optimal path and the action, as already argued in Sec. 5.3.

(a)  Figure 11: Equivalent to Fig. 9, with the average escape time $\langle
T\rangle$ inaccurately (but not too badly) predicted by the full analytical
formula (Eq. 67, including the prefactor) for the overload approximation for
large $(K-P)/\sigma$ (Eq. 67) rather than by WKB.

(a)

(b)

Figure 12: Average desynchronization time $\langle T\rangle$ for Gaussian
fluctuations on the model of Fig. 1(b) and (c) (’uniform’),
$1/\alpha=100\,\text{s}$, (b) $1/\alpha=300\,\text{s}$. Otherwise the same as
Fig. 9(a).

To consider a more interesting situation in which the overload is not along an
edge to a leaf node, we repeat in Figs. 12 and 13 the comparison between the
various approximations for a realization of the Brazilian grid according to
Fig. 1(b) and (c), for which the network partitions into half along two lines.
Additionally, in the network of Fig. 1(c) we consider adding noise with the
same distribution to all different nodes of the network, rather than noise
entering only at two nodes. Figs. 12(a) and 13 show that the results obtained
for the network of Fig. 1(a) still hold with these adaptations. The analysis
of [11] finds that appropriate values for $1/\alpha$ are $2$-$5$ minutes; to
consider the higher end of this range, we consider $1/\alpha=300\,\text{s}$ in
Fig. 12(b) and find that the same results hold as for $1/\alpha=100\,\text{s}$
in Fig. 12(a).

(a)

(b)

Figure 13: Equivalent of Fig. 12 with fluctuations distributed according to
the (rescaled) wind data, $\delta t=10^{-4}\,\text{s}$, and each data point an
average over $1800$ samples. (a) Fluctuations entering in two nodes (Fig.
1(b)), (b) fluctuations entering at all nodes (Fig. 1(c)).

###  6.3 Skewness of the action: more or less desynchronization events for
non-Gaussian data

As we have derived in Sec. 4.3.4, a third cumulant (skewness) in the noise
input can increase or decrease the time to desynchronization. What it actually
does depends on where in the network it enters. As stated earlier, if it
enters in the area with positive power input a positive (negative) skewness
decreases (increases) the average time to desynchronization, while if it
enters in the area with negative power input a positive (negative) skewness
increases (decreases) the average time to desynchronization. A concrete
example of this impact is shown in Fig. 14.

Figure 14: Impact of the skewness on the average desynchronization time
$\langle T\rangle$, showing that a positive skewness in the noise entering in
an area with positive power production $P$ decreases $\langle T\rangle$ and
vice-versa. For further explanations see the text.

In the legend of Fig. 14 ’positive skewness’ refers to $\Gamma^{\xi}_{3}$ of
Eq. 61. The noise terms in areas $1$ (positive production) and $2$ (negative
production), $\xi_{1}$ and $\xi_{2}$ respectively, are constructed to have
tuneable $\Gamma^{\xi}_{3}$. Similar to [8], we construct compound Poisson
processes according to Sec. 3.2. For the case of positive skewness we choose
$\rho=101$ and a jump size distribution

$\lambda_{1}(z)=\frac{100}{101}\delta\Big{(}z+\frac{1}{\sqrt{100\cdot
1^{2}+1\cdot
100^{2}}}\Big{)}+\frac{1}{101}\delta\Big{(}z-\frac{100}{\sqrt{100\cdot
1^{2}+1\cdot 100^{2}}}\Big{)}$ (69)

for $\xi_{1}$, and

$\lambda_{2}(z)=\lambda_{1}(-z)$ (70)

for $\xi_{2}$. By Eqs. 22 and 13 this gives cumulants

$\displaystyle\Gamma_{n}^{\xi_{1}}$ $\displaystyle=$ $\displaystyle
100\cdot\Big{(}\frac{-1}{\sqrt{100\cdot 1^{2}+1\cdot
100^{2}}}\Big{)}^{n}+1\cdot\Big{(}\frac{100}{\sqrt{100\cdot 1^{2}+1\cdot
100^{2}}}\Big{)}^{n}$ $\displaystyle\Gamma_{n}^{\xi_{2}}$ $\displaystyle=$
$\displaystyle\Gamma_{n}^{\xi_{1}}\;\text{for even $n$}$
$\displaystyle\Gamma_{n}^{\xi_{2}}$ $\displaystyle=$
$\displaystyle-\Gamma_{n}^{\xi_{1}}\;\text{for odd $n$},$
$\displaystyle\Gamma^{\xi}_{n}$ $\displaystyle=$
$\displaystyle\frac{1}{2^{n-1}}\big{(}\Gamma_{n}^{\xi_{1}}-\Gamma_{n}^{\xi_{2}}).$
(71)

For the case of negative skewness we choose the same but with $\xi_{1}$ and
$\xi_{2}$ swapped. The case of ’no skewness’ corresponds to the same choice
but with $(1/\sqrt{2})\cdot(\xi_{1}+\xi_{2})$ entering (independently) at each
node. The three cases have identical even cumulants $\Gamma^{\xi}_{n}$, and
positive, negative, and zero odd cumulants respectively. Eq. 61 predicts the
case with $\Gamma^{\xi}_{3}>0$ ($\Gamma^{\xi}_{3}<0$) to decrease (increase)
$\langle T\rangle$ relative to the case with $\Gamma^{\xi}_{3}=0$, as is here
confirmed by calculating $\langle T\rangle$ from the overload prediction of
Eq. 59, truncated after the first $8$ cumulants (and otherwise the same
parameters as in Fig. 9).

##  7 Summary and Conclusions

We considered rare desynchronization events in a coarse-grained form of the
Brazilian grid with prototypical data assignments for inertia, damping, power
input and time correlations of the fluctuating production. As model we used
the swing equations with stochastic power input that is supposed to reproduce
histograms of measured power increment data with non-Gaussian features. We
proposed a so-called Fourier-method to reconstruct the measured histograms
that avoids inherent approximations in an alternative reconstruction by means
of a compound Poisson process. For solar data both approaches lead to clear
differences in the increment histograms and increment cumulants. Our approach
provides a description of power fluctuations in terms of the cumulants of its
increments (rather than its moments, as in the compound Poisson approach),
whose convenient properties in summations allow for an understanding of non-
Gaussian effects in heterogeneous power grids.
In view of large power grids, a dimensional reduction of the high-dimensional
phase space is rather appreciated. Here we extended the so-called synchronized
subgraph approximation [8] to the swing equations, arbitrary topologies and
heterogeneous parameters. Although the partition into (subgraph) areas of the
full grid is suggested by the values of the Fiedler vector at and close to the
bifurcation point, the partition into disjoint synchronized areas after an
overload remains applicable when a large fluctuation is necessary to induce
the desynchronization far off the bifurcation point.
For really rare events, so rare that a sampling of the swing equations for the
original or the reduced grid is unfeasible time-consuming, we applied the WKB-
approach. The WKB-approach amounts to solving classical Hamilton’s equations
of motion. In accordance with our different data implementation, we derive
Hamilton’s equations of motion directly in terms of (the derivative of) the
cumulant generating function of the stochastic power increment input; this
means for applying the WKB-approach, we have to calculate only the cumulants
of the power increment input data. Due to the sudden drop of the involved
auxiliary momenta when the desynchronization is due to an overload, the WKB-
approach is only feasible when the phase space is further reduced. Via the
subgraph approximation phase space gets reduced to twelve dimensions, and the
solutions have been found by the iterative action minimization method.
A further reduction to two dimensions was based on the so-called overload
approximation. This approximation applies when the optimal path between the
stable and saddle fixed point, at which the system escapes in the WKB-
approach, coincides with the time evolution of the phase difference between
the two oscillators (areas) as a function of the power input. The trajectory
of this phase difference follows the same equation that determines the fixed
points, that is, which is obtained via dropping damping and inertia terms.
Here an analytical solution can be found that gives the right order of
magnitude of the average desynchronization time and agrees with results of the
other approaches.
The action in the WKB approach is only determined up to a constant which is
independent of the noise strength. For small enough noise strength or a
sufficiently large distance to the bifurcation point, WKB is justified and all
four approximations agree well with the results of the stochastic swing
equations for the full grid.
For an increased contribution of renewables, an earlier expectation was to
observe more frequent blackout events due to non-Gaussian tails of the power
increments. Like [8] (but derived differently from [8]) we have shown that the
skewness of the power increment distribution can lead to either more or less
desynchronization events, depending on where in the grid the fluctuating input
enters. Neither does lower inertia due to more renewables necessarily lead to
more instabilities. As we have shown, for our realistic parameter values, the
order of magnitude of the average time to desynchronization is independent of
inertia and damping. Inertia and damping should not be considered in
isolation; what matters is their relation to the time correlations of the
fluctuations. Tuning the correlation time $1/\alpha$ amounts to tuning the
inertia and damping.

In general, the role of time correlations in the fluctuating power production
is not negligible. Our results have been obtained for
$1/\alpha=100-300\;\text{s}$, corresponding to measured correlations in [11].
It is $\alpha$ that determines the degree to which the results are impacted by
non-Gaussianity of the fluctuating input. Along with the strength of the noise
$\sigma$ it enters the action in the WKB-approach. Setting the correlation
time $1/\alpha=0\;\text{s}$ as for white noise is definitely not justified for
typical wind or solar data. Further work should therefore extend this
investigation of the interplay between non-Gaussianity and finite time-
correlations also to account for more realistic power spectra [3].

Due to the exponential sensitivity of the average time to desynchronization,
small differences in parameters may get strongly amplified in the very
prediction of how rare blackouts are. Therefore care is needed with respect to
results, for which parameters or data are homogenized for simplicity.

## Acknowledgment

We thank Jason Hindes (U.S. Naval Research Laboratory) and Bhumika Thakur
(Jacobs University Bremen) for useful discussions. Financial support of the
Bundesministerium für Bildung und Forschung (BMBF) (grant number 03EK3055D) is
gratefully acknowledged.

## References

  * [1] J. Kabouris and F. Kanellos, “Impacts of Large Scale Wind Penetration on Energy Supply Industry,” Energies, vol. 2, pp. 1031–1041, Nov. 2009.
  * [2] D. Witthaut, F. Hellmann, J. Kurths, S. Kettemann, H. Meyer-Ortmanns, and M. Timme, “Collective nonlinear dynamics and self-organization in decentralized power grids,” Reviews of Modern Physics, vol. 94, p. 015005, Feb. 2022.
  * [3] K. Schmietendorf, J. Peinke, and O. Kamps, “The impact of turbulent renewable energy production on power grid stability and quality,” The European Physical Journal B, vol. 90, p. 222, Nov. 2017.
  * [4] B. Schäfer, M. Matthiae, X. Zhang, M. Rohden, M. Timme, and D. Witthaut, “Escape routes, weak links, and desynchronization in fluctuation-driven networks,” Physical Review E, vol. 95, p. 060203, June 2017.
  * [5] D. Manik, D. Witthaut, B. Schäfer, M. Matthiae, A. Sorge, M. Rohden, E. Katifori, and M. Timme, “Supply networks: Instabilities without overload,” The European Physical Journal Special Topics, vol. 223, pp. 2527–2547, oct 2014.
  * [6] J. Machowski, J. W. Bialek, and J. R. Bumby, Power system dynamics: stability and control.  Chichester, U.K: Wiley, 2nd ed ed., 2008.  OCLC: ocn232130756.
  * [7] T. E. P. on the separation of the Continental Europe Synchronous Area of 08 January 2021, “Continental europe synchronous area separation on 08 january 2021,” jul 2021.  Official report on the European system split of 08-01-2021. Available at: https://www.entsoe.eu/news/2021/07/15/final-report-on-the-separation-of-the-continental-europe-power-system-on-8-january-2021/.
  * [8] J. Hindes, P. Jacquod, and I. B. Schwartz, “Network desynchronization by non-gaussian fluctuations,” Physical Review E, vol. 100, p. 052314, nov 2019.
  * [9] M. I. Dykman, E. Mori, J. Ross, and P. M. Hunt, “Large fluctuations and optimal paths in chemical kinetics,” The Journal of Chemical Physics, vol. 100, pp. 5735–5750, apr 1994.
  * [10] M. Dykman, I. Schwartz, and A. Landsman, “Disease Extinction in the Presence of Random Vaccination,” Physical Review Letters, vol. 101, p. 078101, Aug. 2008.
  * [11] M. Tyloo, J. Hindes, and P. Jacquod, “Finite-time Correlations Boost Large Voltage-Angle Fluctuations in Electric Power Grids,” Mar. 2022.  Number: arXiv:2203.00590 arXiv:2203.00590 [cond-mat, physics:nlin, physics:physics].
  * [12] T. Nishikawa and A. E. Motter, “Comparative analysis of existing models for power-grid synchronization,” New Journal of Physics, vol. 17, p. 015012, jan 2015.
  * [13] M. Anvari, G. Lohmann, M. Wächter, P. Milan, E. Lorenz, D. Heinemann, M. R. R. Tabar, and J. Peinke, “Short term fluctuations of wind and solar power systems,” New Journal of Physics, vol. 18, p. 063027, jun 2016.
  * [14] V. Lopes, T. Scholz, A. Estanqueiro, and A. Novais, “On the use of Markov chain models for the analysis of wind power time-series,” 2012 11th International Conference on Environment and Electrical Engineering, pp. 770–775, May 2012.
  * [15] H. Haehne, K. Schmietendorf, S. Tamrakar, J. Peinke, and S. Kettemann, “Propagation of wind-power-induced fluctuations in power grids,” Physical Review E, vol. 99, p. 050301, May 2019.
  * [16] N. G. v. Kampen, Stochastic processes in physics and chemistry.  North-Holland personal library, Amsterdam ; New York: North-Holland, rev. and enl. ed., 1992.
  * [17] C. W. Gardiner and C. W. Gardiner, Stochastic methods: a handbook for the natural and social sciences.  Springer series in synergetics, Berlin: Springer, 4th ed ed., 2009.
  * [18] A. B. Birchfield, T. Xu, K. M. Gegner, K. S. Shetye, and T. J. Overbye, “Grid structural characteristics as validation criteria for synthetic networks,” IEEE Transactions on Power Systems, vol. 32, no. 4, pp. 3258–3265, 2017\.  Data available at: https://electricgrids.engr.tamu.edu/electric-grid-test-cases/.
  * [19] T. Xu, A. B. Birchfield, K. S. Shetye, and T. J. Overbye, “Creation of synthetic electric grid models for transient stability studies,” Proc.10th Bulk Power Systems Dynamics and Control Symposium (IREP 2017), sep 2017\.
  * [20] T. Xu, A. B. Birchfield, and T. J. Overbye, “Modeling, tuning, and validating system dynamics in synthetic electric grids,” IEEE Transactions on Power Systems, vol. 33, no. 6, pp. 6501–6509, 2018.
  * [21] J. Machowski, P. Kacejko, S. Robak, P. Miller, and M. Wancerz, “Simplified angle and voltage stability criteria for power system planning based on the short-circuit power: SIMPLIFIED ANGLE AND VOLTAGE STABILITY CRITERIA,” International Transactions on Electrical Energy Systems, vol. 25, pp. 3096–3108, Nov. 2015.
  * [22] Available at http://www.uni-oldenburg.de/fileadmin/user_upload/physik/ag/twist/Forschung/Daten/AnvariEtAl2016_ExampleData_WindPowerAndIrradianceData.zip.
  * [23] P. J. Smith, “A Recursive Formulation of the Old Problem of Obtaining Moments from Cumulants and Vice Versa,” The American Statistician, vol. 49, pp. 217–218, May 1995.
  * [24] P. Protter and D. Talay, “The Euler Scheme for Levy Driven Stochastic Differential Equations,” The Annals of Probability, vol. 25, no. 1, pp. 393–423, 1997.  Publisher: Institute of Mathematical Statistics.
  * [25] R. Chen, J. Hu, Y. Peng, and ,School of Management, Fudan University, Guoshun Road 670, Shanghai 200433, “Simulation of levy-Driven models and its application in finance,” Numerical Algebra, Control & Optimization, vol. 2, no. 4, pp. 749–765, 2012.
  * [26] S. Blinnikov and R. Moessner, “Expansions for nearly Gaussian distributions,” Astronomy and Astrophysics Supplement Series, vol. 130, pp. 193–205, May 1998.
  * [27] J.-P. Chateau and D. Dufresne, “Gram-Charlier Processes and Applications to Option Pricing,” Journal of Probability and Statistics, vol. 2017, p. e8690491, Feb. 2017.  Publisher: Hindawi.
  * [28] F. W. Steutel and K. v. Harn, Infinite divisibility of probability distributions on the real line.  No. 259 in Monographs and textbooks in pure and applied mathematics, New York: Marcel Dekker, 2004.
  * [29] P. e. a. Virtanen, “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python,” Nature Methods, vol. 17, pp. 261–272, 2020.
  * [30] M. Anvari, M. R. R. Tabar, J. Peinke, and K. Lehnertz, “Disentangling the stochastic behavior of complex time series,” Scientific Reports, vol. 6, p. 35435, Dec. 2016.
  * [31] K. Lehnertz, L. Zabawa, and M. R. R. Tabar, “Characterizing abrupt transitions in stochastic dynamics,” New Journal of Physics, vol. 20, p. 113043, Nov. 2018.
  * [32] B. S. Lindley and I. B. Schwartz, “An iterative action minimizing method for computing optimal paths in stochastic dynamical systems,” Physica D: Nonlinear Phenomena, vol. 255, pp. 22–30, July 2013.
  * [33] H. Touchette, “The large deviation approach to statistical mechanics,” Physics Reports, vol. 478, pp. 1–69, July 2009.
  * [34] M. Assaf and B. Meerson, “WKB theory of large deviations in stochastic populations,” Journal of Physics A: Mathematical and Theoretical, vol. 50, p. 263001, jun 2017.
  * [35] S. R. Serrao, T. Ritmeester, and H. Meyer-Ortmanns, “Rare extinction events in cyclic predator–prey games,” Journal of Physics A: Mathematical and Theoretical, vol. 54, p. 235001, June 2021.
  * [36] M. I. Freĭdlin, A. D. Wentzell, and J. Szücs, Random perturbations of dynamical systems.  No. 260 in Grundlehren der mathematischen Wissenschaften, Heidelberg: Springer, 3d edition ed., 2012.
  * [37] A. Dembo and O. Zeitouni, Large Deviations Techniques and Applications, vol. 38 of Stochastic Modelling and Applied Probability.  Berlin, Heidelberg: Springer Berlin Heidelberg, 2010.
  * [38] L. H. Chan, K. Chen, C. Li, C. W. Wong, and C. Y. Yau, “On higher-order moment and cumulant estimation,” Journal of Statistical Computation and Simulation, vol. 90, pp. 747–771, Mar. 2020.
  * [39] a. k. Singh and B. Bhadauria, “Finite Difference Formulae for Unequal Sub-Intervals Using Lagrange’s Interpolation Formula,” Journal of Math. Analysis, vol. 3, pp. 815–827, jan 2009.
  * [40] J. J. Moré, “The Levenberg-Marquardt algorithm: Implementation and theory,” in Numerical Analysis (G. A. Watson, ed.), Lecture Notes in Mathematics, (Berlin, Heidelberg), pp. 105–116, Springer, 1978.
  * [41] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, I. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, and P. van Mulbregt, “SciPy 1.0: fundamental algorithms for scientific computing in Python,” Nature Methods, vol. 17, pp. 261–272, mar 2020\.  Documentation of the Levenberg-Marquardt algorithm available at https://docs.scipy.org/doc/scipy/reference/optimize.root-lm.html.
  * [42] M. Tyloo and P. Jacquod, “Primary Control Effort Under Fluctuating Power Generation in Realistic High-Voltage Power Networks,” IEEE Control Systems Letters, vol. 5, pp. 929–934, July 2021.
  * [43] P. Chen, Z. Chen, and B. Bak-Jensen, “Probabilistic load flow: A review,” in 2008 Third International Conference on Electric Utility Deregulation and Restructuring and Power Technologies, (Nanjing, China), pp. 1586–1591, IEEE, Apr. 2008.
  * [44] X. Zhang, S. Hallerberg, M. Matthiae, D. Witthaut, and M. Timme, “Fluctuation-induced distributed resonances in oscillatory networks,” Science Advances, vol. 5, p. eaav1027, July 2019.
  * [45] L. M. Ricciardi and S. Sato, “First-passage-time density and moments of the ornstein-uhlenbeck process,” Journal of Applied Probability, vol. 25, pp. 43–57, Mar. 1988.
  * [46] F. Bouchet and J. Reygner, “Generalisation of the Eyring–Kramers Transition Rate Formula to Irreversible Diffusion Processes,” Annales Henri Poincaré, vol. 17, pp. 3499–3532, Dec. 2016.
  * [47] UCTE, “Ucte operation handbook,” 2004.  Available at: https://eepublicdownloads.entsoe.eu/clean-documents/pre2015/publications/entsoe/Operation_Handbook/Policy_1_Appendix%20_final.pdf.

##  Appendix A Parameters for the reduced Brazilian grid

The data ([18]) contains the topology of the grid, the reactance for each line
(equal to $1/K_{ij}$ when resistance is neglected, as is the case for the
swing equations) and the generators and loads connected to each node, along
with the production/consumption at these nodes. The data contains a
transformer, with a node for each end. As no (relevant) data for the
transformer is given, and we cannot model the transformer within the swing
equations (without voltage dynamics), we treat this as a single node. For
synchronous generators and synchronous loads the constants of inertia $H_{i}$
are given, from which $M_{i}=\frac{2H_{i}}{\omega}|P_{i}|$ can be determined
(where $|P_{i}|$ is the production/consumption at the generator/load). For the
Brazilian grid, the angular frequency $\omega=2\pi\cdot 60\;\text{Hz}$. For
these generators and loads, we assume (as argued in [12]) that the damping
coefficient is dominated by a droop of $2\%$ (see [6] for a detailed
description of droop). The damping of inductive generators and loads is due to
a small self-regulating effect of $\approx 1\%/\text{Hz}$ [47]; for these
devices the inertia is set to zero. If there are multiple generators/loads
connected to a node, we sum their inertia, damping, power input. One of the
nodes does not have any inertia; at this node we set the inertia to a small
value ($0.01\,\text{s}^{2}$) to avoid numerical instability.

This gives the following sets of parameters, here shown at $9.1\%$ from
overload, in the per unit system (i.e. in units of a reference power that
itself is here not further specified, see [12, 6]):

$\displaystyle\bm{P}$
$\displaystyle=\begin{pmatrix}-1.38304&1.17679&1.56501&12.38545&-1.05426&-12.68997\end{pmatrix}\,,$
(72) $\displaystyle\bm{M}$
$\displaystyle=\begin{pmatrix}0.06201&0.04950&0.05716&0.29119&0.01&0.14902\end{pmatrix}\,\text{s}^{2}\,,$
(73) $\displaystyle\bm{\gamma}$
$\displaystyle=\begin{pmatrix}0.35049373&0.27670489&0.31923646&1.43640985&0.11831416&0.75322515\end{pmatrix}\,\text{s}\,,$
(74) $\displaystyle\bm{K}$
$\displaystyle=\begin{pmatrix}0.&0.&26.316&0.&4.082&0.\\\
0.&0.&13.158&0.&4.444&0.\\\ 26.316&13.158&0.&0.&0.&0.\\\
0.&0.&0.&0.&13.624&0.\\\ 4.082&4.444&0.&13.624&0.&17.544\\\
0.&0.&0.&0.&17.544&0.\end{pmatrix}\,.$ (75)

##  Appendix B Derivation of the differential Chapman-Kolmogorov equation

Under mild assumptions, further specified in [17], for any continuous time
Markov process, here for the swing equations with fluctuations according to
Eq. 2, the time evolution of
$\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)$ is given by an integro-
differential equation called the differential Chapman-Kolmogorov equation
[17]. For $\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)$ this equation
can be derived from Eqs. 1-2 by first considering the discrete evolution of
$\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)$:

$\displaystyle\frac{\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t+\Delta
t)-\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)}{\Delta t}$
$\displaystyle=\frac{1}{\Delta
t}\int\operatorname{\text{d}\\!}\bm{\phi}^{\prime}\operatorname{\text{d}\\!}\bm{v}^{\prime}\operatorname{\text{d}\\!}\bm{p}^{\prime}$
(76)
$\displaystyle\times\Big{[}\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t+\Delta
t|\bm{\phi}^{\prime},\bm{v}^{\prime},\bm{p}^{\prime},t)\operatorname{\text{pr}}(\bm{\phi}^{\prime},\bm{v}^{\prime},\bm{p}^{\prime},t)$
$\displaystyle-\operatorname{\text{pr}}(\bm{\phi}^{\prime},\bm{v}^{\prime},\bm{p}^{\prime},t+\Delta
t|\bm{\phi},\bm{v},\bm{p},t)\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)\Big{]}\,.$

For a given realization of the noise $\xi$, the transition probabilities
$\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t+\Delta
t|\bm{\phi}^{\prime},\bm{v}^{\prime},\bm{p}^{\prime},t)$ and
$\operatorname{\text{pr}}(\bm{\phi}^{\prime},\bm{v}^{\prime},\bm{p}^{\prime},t+\Delta
t|\bm{\phi},\bm{v},\bm{p},t)$ are simply delta functions enforcing the swing
equations (Eqs. 1-2). For small $\Delta t$, we get,

$\displaystyle\frac{1}{\Delta
t}\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t+\Delta
t|\bm{\phi}^{\prime},\bm{v}^{\prime},\bm{p}^{\prime},t)$
$\displaystyle=\frac{1}{\Delta
t}\prod_{i}\Big{[}\int\operatorname{\text{d}\\!}z_{i}\operatorname{\text{pr}}(Z_{\Delta
t,i}=z_{i})\times\delta\Big{(}(\phi_{i}-\phi_{i}^{\prime})-\Delta t\cdot
v_{i}\Big{)}$
$\displaystyle\times\delta\Big{(}(v_{i}-v_{i}^{\prime})-\frac{\Delta
t}{M_{i}}\cdot\big{[}-\gamma_{i}v_{i}^{\prime}+\overline{P}_{i}+p_{i}^{\prime}+\sum_{j}K_{ij}\sin(\phi_{j}-\phi_{i})\big{]}\Big{)}$
$\displaystyle\times\delta\Big{(}(p_{i}-p^{\prime}_{i})-[-\alpha_{i}p_{i}^{\prime}\Delta
t+\sigma_{i}z_{i}]\Big{)}\Big{]}\,,$ (77)

and equivalent expressions for
$\operatorname{\text{pr}}(\bm{\phi}^{\prime},\bm{v}^{\prime},\bm{p}^{\prime},t+\Delta
t|\bm{\phi},\bm{v},\bm{p},t)$. Expanding the delta functions to first order in
$\Delta t$ (with the derivative of the delta function defined such that
$\int\operatorname{\text{d}\\!}xf(x)\frac{\text{d}\delta(x)}{\text{d}x}=-\frac{\operatorname{\text{d}\\!}f(x)}{\operatorname{\text{d}\\!}x}\biggr{\rvert}_{0}$),
the limit $\Delta t\rightarrow 0$ leads from Eq. 76 to:

$\displaystyle\frac{\partial\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)}{\partial
t}=\sum_{i}\Big{[}$
$\displaystyle-\frac{\partial}{\partial\phi_{i}}[v_{i}\cdot\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)]$
$\displaystyle-\frac{\partial}{\partial
v_{i}}\big{[}\frac{1}{M_{i}}\big{(}-\gamma_{i}v_{i}+\overline{P}_{i}+p_{i}+\sum_{j}K_{ij}\sin(\phi_{j}-\phi_{i})\big{)}\cdot\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)\big{]}]$
$\displaystyle-\frac{\partial}{\partial
p_{i}}\big{[}-\alpha_{i}p_{i}\cdot\operatorname{\text{pr}}(\bm{\phi},\bm{v},\bm{p},t)\big{]}$
$\displaystyle+\lim_{\Delta t\rightarrow
0}\int\operatorname{\text{d}\\!}z_{i}\frac{\operatorname{\text{pr}}(Z_{\Delta
t,i}=z_{i})}{\Delta