# The stochastic nature of power-grid frequency in South Africa

Leonardo Rydin Gorjão  Department of Computer Science, OsloMet – Oslo
Metropolitan University, N-0130 Oslo, Norway  Faculty of Science and
Technology, Norwegian University of Life Sciences, 1432 Ås, Norway  Jacques
Maritz  Department of Engineering Sciences, University of the Free State,
Bloemfontein 9301, South Africa

###### Abstract

In this work, we explore two mechanisms that explain non-Gaussian behaviour of
power-grid frequency recordings in the South African grid. We make use of a
Fokker–Planck approach to power-grid frequency that yields a direct relation
between common model parameters such as inertia, damping, and noise amplitude
and non-parametric estimations of the same directly from power-grid frequency
recordings. We propose two explanations for the non-Gaussian leptokurtic
distributions in South Africa: The first based on multiplicative noise in
power-grid frequency recordings, which we observe in South Africa; The second
based on the well-known scheduled and unscheduled load shedding and rolling
blackouts that beset South Africa. For the first we derive an analytic
expression of the effects of multiplicative noise that permits the estimation
of all statistical moments – and discuss drawbacks in comparison with the
data; For the second we employ a simple numerical analysis with a modular
power grid of South Africa. Both options help understand the statistics of
power-grid frequency in South Africa – particularly the presence of heavy
tails.

##  I Introduction

Power-grid systems power the world. The developed countries rely on the
availability of electrical power to the extent that no less than a 24-hour
service is conceivable. Yet this is a reality only in parts of the world, so
much that this is one of the key elements of the Sustainable Development Goals
(SDG) – SDG 7: “Ensure access to affordable, reliable, sustainable and modern
energy for all” [1]. One curious example of a thus far poorly managed power
system is the case of South Africa, a country still dominantly supplied via
coal-fired power plants [2]. At least 12% of the population does not have
access to power and roughly 10% cannot adequately afford electricity [3, 4,
5], particularly in rural areas [6]. Nevertheless, the power system has
continuously evolved and has successfully attracted renewable-energy centred
investments [7, 2]. Yet these developments stand on a frail and deteriorating
infrastructure. Since 2007 the maintenance and development of the South
African grid has been hindered by mismanagement [8, 9]. Mismanagement has
manifested itself in recent years with its inability to deliver power
according to the country’s demand. This lead to the implementation of ‘rolling
blackouts’ load shedding events across the country [10]. Load shedding has
marked deleterious societal effects [11, 12, 13]. In 2021, the citizens and
industries of South Africa were afflicted by a lack of power and periodic load
shedding for over $48$ days of the year. Additional to load shedding, the
utility experiences the overloading of transformers and substations due to
illegal connections, overloading, and vandalism. The latter adds to unplanned
outages (also known as non-technical losses) for parts of the network and
accounts for $10\%$ of the non-technical losses [14]. Energy losses
experienced by the utility are due to meter tampering, incorrect billing, and
cable theft [15]. The overloading of electrical infrastructure leads to
mandatory disconnection of municipalities from the distribution network.

Lack of power supply or the need for abrupt load-shedding leads to visible
changes in power-grid frequency, the key signature of a power-grid system [16,
17]. Alternating current power grids operate at a nominal frequency – $50$ Hz
in South Africa – which means all conventional power-generating units contain
a mechanical rotating generator tuned to this frequency. If the power
generation falls short of the power demand, frequency sags. If the power
generator surpasses demand, frequency rises. Power grids are complex systems
with inherently complex stochastic dynamics [18, 19]. Power generation, both
conventional but particularly renewable, is manifestly volatile [20, 21, 22,
23] and non-Gaussian behaviour permeates the statistics of most power system
[19, 24, 25]. Similarly, the demand side can be equated to a complex
stochastic process [26]. Thus, we take a complex-system approach to modelling
and analysing power-grid systems from a complex system perspective [27, 28].

In this work, we will examine three power-grid frequency recordings from 2021,
from Continental Europe, the Nordic Grid, and the South African grid, with a
focus on the latter. We investigate these recordings and contrast their
statistical and stochastic characteristics by employing a phenomenological
Fokker–Planck equation of an aggregated swing-equation model with stochastic
noise [29, 30, 31, 32]. We seek an explanation for the non-Gaussian
distribution of power-grid frequency observed in South Africa, wherein we
offer two explanations for the heavy tails of the distributions [18, 19, 24].
Firstly we identify the potential presence of multiplicative noise in the
recordings and find an analytic approximation for its effect on the
statistical moments of Langevin processes with multiplicative noise. We
contrast this with two other recordings of power-grid frequency in South
Africa and discuss the limitations of such an approach. Lastly, we simulate
power-grid frequency trajectories in a model grid of South Africa and mimic
the effect of load shedding. Both approaches help us explain the non-Gaussian
statistics of the power-grid frequency recordings of South Africa.

This work is organised as follows. In Sec. II we shortly introduced the
historical context situating the South African power grid, the ongoing load-
shedding events, and we showcase the direct effects at the University of the
Free State in Bloemfontein. In Sec. III we introduce both a low-level
dynamical model of power-grids affected by noise as well as Fokker–Planck
approximation of the power-grid frequency statistics. From this basis, we
present our theoretical and numerical results. In Sec. III we discuss the
outcomes and limitations of our data examination as extend some thoughts on
future work and the necessity to examine previously overlooked power grids.

##  II Background

###  II.1 Load-shedding in South Africa

In this work, we focus on power-grid frequency recordings from the Continental
European power grid, the Nordic Grid, and the South African power grid. It is
the latter that interests us the most for this examination. The South African
power grid is undergoing an unprecedented transformation.

The power sector in South Africa is dominated by the state-owned entity ESKOM,
which both imports and exports energy from surrounding states with some input
from nuclear and renewable energy (photovoltaic (PV), wind and concentrated
solar power (CSP) [33]. Unfortunately, the South African power system has been
experiencing an energy crisis since 2007 with electrical generation lagging
the electrical demand, ultimately leading to implemented rolling blackouts in
an effort to stabilise the national grid.

Due to extended periods of lack of maintenance and mismanagement resulted in
an unpredictable and unreliable power system in South Africa, that during 2021
experienced load-shedding for $1169$ hours with $1775$ GWh (with the upper
limit of $2521$ GWh) energy shed. The latter amounts to 13.3$\%$ of all hours
of 2021 being in the state of, mostly, Stage 2 load-shedding (for details see
Tab. 1 in Sec. III.5) [10]. The power system remained coal powered during
2021, but with $18.6\%$ contributed by renewable energy (of which $6.7\%$ were
variable sources) [2]. Concernedly, $3.2$ TWh of energy was generated by using
diesel.

With critical low reserve margins and increasing load-shedding, some consumers
are considering the implementation of self-sufficient microgrids with some
renewable input. The latter is preceded by the digitisation of the existing
electrical network topology (i.e., smart grids) that requires initial capital
input [34]. South African university campuses are naturally forced to evolve
according to the constraints of the national utility, the ability to swiftly
react to demand reduction signals being the most critical, see Fig. 1. In
cases of severe loss of supply, self-sufficient microgrids are utilised that
synchronise with available renewable supply to ensure supply continuity on
campus. The latter typically being hybrid PV-Diesel based microgrids.

During load-shedding events (i.e. demand overpowering generation), the utility
reduces load in stages of 1000 MW. A total of 6 stages of load-shedding are in
place in South Africa at this moment, with varying duration, usually designed
to stay in place for $i\times 6$ hours, with $i$ the stage. An overview of the
last three years of total load shed energy can be found below in Tab. 1 (and
in Ref. [10].). One should note that 1) load-shedding upper limits are never
actually achieved; 2) load-shedding hours are uncertain due to additional
switching time; 3) additional outages exist after switching. Unfortunately,
access to data relating to all these events is hard, and thus we focus
primarily on modelling uncertain load-shedding events with an equivalent
proportion of hours affected vs total hours in 2021. Therein we consider the
occasional events of total blackout, which in our case is represented by the
load of a region vanishing for a period of time.

###  II.2 The visible consequences at the University of the Free State

The University of the Free State (UFS) Bloemfontein campus has a remarkable
resource in the form of a campus smart grid. The latter enables granular
control and monitoring of intensive loads based on optimised scenarios or in-
situ demand response requests from the utility. Dual metering on both sides of
the utility feed enables monitoring of the UFS load, during instances of load-
shedding, to determine the reduction needed to comply with the requirements
set out by the utility. Maintaining high levels of control granularity ensures
that the institution can act as a packetised energy manager [35], hence
responding to the exact requirements of the power grid operator and generating
the opportunity of ancillary frequency support to the power system using the
available generating resources on campus. Fig. 1 illustrates the concept of
acting as a packetised energy manager, ultimately responding to the
requirements of the utility to reduce demand for a specified duration. Demand
response to utility curtailment notices is driven by the ability to reduce
demand for a negotiated time window while remaining connected to the power
grid, hence frequency support is inherent, especially in the case of UFS
campus power systems that are largely supported by photovoltaic renewable
generators. De-loaded photovoltaic plants can provide frequency support via
the control of plant active power [36]. Load-shedding becomes the critical
element when designing rural microgrids, with the emphasis moving from
economic sense to power continuity, ultimately, converging in microgrids that
couple with existing PV generation and require islanded inertia and sensitive
frequency control.

Figure 1: Demand-side management (as response to utility load-shedding) at the
University of the Free State (UFS) Campus in Bloemfontein, South Africa, in
comparison with a typical day. This is one solution among many that public and
private institutions have had to implement to be able to operate
uninterruptedly. Displayed are the 19th and 20th of May 2022, the first of
this a day with campus demand response totalling $\sim\\!\\!12\,$GWh, the
second a typical day without the call of demand response

##  III Statistics of power-system frequency recordings

###  III.1 The swing equation

An operating AC power-grid system is a set of highly synchronised oscillators
operating at a rotational frequency multiple of $50$ Hz (or $60$ Hz). In a
reduced format, we can consider a general model for power grids as coupled
inertial oscillators in a network. A synchronous machine $j$, a generator or a
load, is described by its rotor-angle $\theta_{j}(t)$ and its angular
frequency $\omega_{j}(t)=\dot{\theta}_{i}(t)$. It obeys the equations of
motion [17, 37, 38, 39]

$\displaystyle\dot{\theta}_{j}$ $\displaystyle=\omega_{j},$ (1) $\displaystyle
M_{j}\dot{\omega}_{j}$
$\displaystyle=-D_{j}\omega_{j}+P_{j}^{\mathrm{m}}-P_{j}+\xi_{j},$

where $M_{j}$ is the inertial mass of the $j$-th rotating machine, $D_{j}$ its
damping factor, $P_{j}^{\mathrm{m}}$ comprises the self-generated mechanical
power (or if negative it acts as a load), and $P_{j}$ is the exchanged power
with the other oscillators in the network. The innocently looking exchange
power $P_{j}$ embodies the very complex structure of power-grid systems, i.e.,
$P_{j}$ is given by

$P_{j}=\sum_{\ell=1}^{N}\\!E_{j}E_{\ell}\left[B_{j,\ell}\sin(\delta_{j}\\!-\\!\delta_{\ell})\\!+\\!G_{j,\ell}\cos(\theta_{j}\\!-\\!\theta_{\ell})\right],$
(2)

where $E_{j}$ is the transient voltage at $j$-th machine, which is our reduced
swing equation model (1) is static. The parameters $G_{j,\ell}$ and
$B_{j,\ell}$ denote the real and imaginary parts of the nodal admittance
matrix $Y_{j,\ell}$ and encode the network structure. The conductance
$G_{j,\ell}$ comprises all the resistive terms of the network which we will
not deal with, i.e., we will consider the case of a purely lossless system
$G_{j,\ell}=0,\forall j,\ell$. Conductance effects can be disregarded as they
are necessarily balanced in a functional power system where the generated
power is higher than the consumed one to compensate for these losses. These
models are also known as second-order/inertial Kuramoto models [30, 31, 37,
40].

The network topology is fully embedded in the susceptance matrix $B_{j,\ell}$
which, apart from shunt susceptances, is equivalent to the graph Laplacian of
the network [41]. It should be noted that the swing equation (1) or other more
complex models only have a descriptive character of what happens in
functioning power-grid systems. For example, in our case, we will simulate
power-grid frequency trajectories in South Africa using a very reduced model
of the real-world South African power-grid system, wherein we consider a
reduced number of nodes in a coupled graph of interacting generators/loads.
These coarse-grained models are illustrative at best – yet they capture with
surprising accuracy the dynamics of the rotor angle $\theta_{j}$ and angular
frequency $\omega_{j}$ of the generators/loads in real systems.

Lastly, we have included a noise term $\xi_{j}$ at each node. In principle, we
could ask why there is no noise term in the rotor-angle equation. One can
similarly include a noise element in the angles, which will manifest itself in
some format on the frequency. Yet here we arrive at the core of what
stochastic noise $\xi_{j}$ is. On the one hand, we know we can broadly model
complex power grids with simplified models as the swing equation (1). On the
other hand, thousands of power-generating units and millions of consumers are
‘condensed’ into a single machine/load, a node in our network. One manifest
effect is that real-world recording of voltage angles, power-grid frequency,
power flows, or any other physically relevant measures, are noisy. These are
not solely dynamic processes, these are stochastic processes, noisy processes,
whence we can learn plenty from the system by examining the statistical
properties of various real-world recordings [42].

Figure 2: Three exemplary power-grid frequency trajectories (left) from
Continental Europe (green), Nordic Grid (orange), and South African grid
(purple), between the 2nd and 3rd of January 2021. The recordings are
vertically displaced for clarity. The probability density functions for the
whole of 2021 are shown in the centre and right panels, the latter in a
vertical logarithmic scale. All three power grids exhibit non-Gaussian
phenomena, particularly showing larger than normal kurtosis ($\kappa>3$) and
some skewness of the distributions. The power-grid frequency data from
Continental Europe and the Nordic Grid have a temporal resolution (sampling
rate) of $1$ s. The data from the South African grid has $15$ minutes
resolution ($900$ s). All figures produced with Matplotlib, NumPy, SciPy, and
pandas [43, 44, 45, 46].

###  III.2 Fokker–Planck description of power-grid frequency recordings

Returning to our noise-disturbed swing equation (1) we can now focus on
writing a Fokker–Planck equation for the probability density
$\rho(\theta,\omega,t)$ for a single node [29, 30, 32]. Let us consider the
noise is delta correlated
$\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=2B\delta_{i,j}\delta(t-t^{\prime})$,
with $B$ the diffusion constant. Now considering precisely the previous
argument that our models are descriptive of $N\gg 1~{}~{}(N\to\infty)$, we can
take the ‘mean-field’ approximation yielding [29, 30]

$\frac{\partial\rho}{\partial
t}=-\omega\frac{\partial\rho}{\partial\theta}+\frac{1}{m}\frac{\partial}{\partial\omega}\left[\left(-d\omega+P_{j}^{\mathrm{m}}-P_{j}\right)\rho\right]+\frac{B}{m^{2}}\frac{\partial^{2}}{\partial\omega^{2}}\rho.$
(3)

We assume that $M_{j}=m$ and $D_{j}=m$. We can separate the probability
density $\rho(\theta,\omega,t)=\phi(\theta,t)\eta(\omega,t)$ and focus solely
on the statistics of the angular frequency $\eta(\omega,t)$, following Acebrón
& Spigler [29]. This results in

$\frac{\partial\eta}{\partial
t}=-\frac{d}{m}\frac{\partial}{\partial\omega}\omega\eta+\frac{B}{m^{2}}\frac{\partial^{2}\eta}{\partial\omega^{2}}.$
(4)

The static solution at $t\to\infty$ (and having appropriate decay conditions
of $\eta(\omega)=0$ at $\omega\to\pm\infty$) is the well known Gaussian
distribution

$\eta(\omega)=\sqrt{\frac{md}{2\pi B}}e^{-\frac{md}{2B}\omega^{2}}.$ (5)

There are various important details to notice here. First of all, we have
started with uncorrelated Gaussian white noise, which is the most natural
assumption as it is the simplest noise model. One can easily argue that power-
grid frequencies are not Markovian processes and that these are likely driven
by correlated noises [27, 47, 48, 49]. Non-Markovianity is ubiquitous across
energy systems [24, 50, 51], yet we do not necessarily need to address it
directly while modelling. In this work, we are interested in understanding the
statistics of power-grid frequency recordings in South Africa, for which we
have access to recordings with a time sampling of $15$ minutes. This time
resolution is much beyond the intrinsic Einstein–Markov length of the process,
i.e., the time scale at which a non-Markovian process can safely be modelled
in a Markovian setting [28]. For an inhomogeneous spatially extended parabolic
partial differential equation of coupled swing equations, see Ref. [52].

Figure 3: Drift and diffusion estimated for three power-grid frequency
recordings from Continental Europe, the Nordic Grid, and the South African
grid. (Left panel) The drift $D_{1}$ in the Nordic Grid and South Africa is
linear with a negative slope, indicating a linear response, as expected in
(4). The case of Continental Europe is less trivial, the shaded green area
shows the area wherein there is no control active – the ‘deadband’. (Right
panel) the diffusion $D_{2}$ in Continental Europe and Nordic Grid can be said
to be constant. In South Africa, we observe a functional dependence on angular
frequency. Data covers the year 2021 [53].

Secondly, we can already suspect that our Gaussian description of the
distribution of the frequencies (5) might not be sufficient. Schäfer et al.
[19] report on more complex statistics of power-grid frequency in Continental
Europe and the Nordic Grid, among others, and provide a descriptive analysis
based on a convolution of statistics which results in power-grid frequency
distribution far beyond Gaussianity. Various subsequent works offer
explanations for the aforementioned phenomena [54, 55, 24]. The most
straightforward explanation of the non-Gaussianity of power-grid frequency
statistics in Continental Europe and the Nordic Grid is due to changes in
dispatch. At fixed intervals, usually every hour, with some minor in-between
blocks of 15 minutes, different power plants, big or small, sell their power
generation capacities on an electricity exchange market. This is a fairly
complex problem of bidding and offering, which, for our purpose, boils down to
one known and very clear effect: at every trading block (1 hour or 15
minutes), power-grid frequency sees a large overshoot or a large sag. This
effect comes about because the generation needs to adapt to an ever-changing
consumption, that is, the generation adapts to the predicted consumption for
the following 1-hour block. If there is more generation than consumption,
power-grid frequency goes up from its nominal value of $50$ Hz and vice versa,
lower generation implies a sag of the frequency.

In Fig. 2 (left) we display a snippet of the three frequency recordings
Continental Europe, the Nordic Grid, and the South African grid, between the
2nd and 3rd of January 2021. We can see that these recordings look different
from one another. In the centre panel, we display the probability density
$\rho(f)$ of the three recordings, covering all of 2021. Similarly, we exhibit
the same thing in the right panel on a vertical logarithmic scale. Here we can
notice the evident heavy tails of all three recordings. A purely Gaussian
distribution is shown in black dashed lines and is an inverter parabola in a
vertical logarithmic scale.

To best understand the distribution of power-grid frequency in the various
grid, we can directly examine the skewness $s$ and kurtosis $\kappa$, given
respectively by

$\displaystyle
s_{X}=\frac{\operatorname{E}\left[(X-\mu)^{3}\right]}{\left(\operatorname{E}\left[(X-\mu)^{2}\right]\right)^{3/2}},\quad\kappa_{X}=\frac{\operatorname{E}\left[(X-\mu)^{4}\right]}{\left(\operatorname{E}\left[(X-\mu)^{2}\right]\right)^{2}}.$
(6)

A Gaussian distribution has skewness $s=0$ and a kurtosis $\kappa=3$. That is,
it is purely symmetric and does not have any heavy (or light) tails. Being
facetious, Gaussian distributions are very normal. What we observe is that
none of the recordings are symmetric ($s\neq 0$) and they all show heavy tails
($\kappa>3$). The aforementioned changes in the dispatch/market activity in
Continental Europe and the Nordic Grid explain this behaviour. But what about
South Africa?

###  III.3 Understanding linear response and the effect of noise

We have introduced a Fokker–Planck equation describing the probability density
$\rho(\omega)$, wherein we can directly evaluate the effects of the inertia
$m$ in the system and the presence of noise, mediated by its amplitude $B$. We
can now test, using our frequency recordings, how much our simplistic
Fokker–Planck model (4) can describe the data. In a more general setting, we
write a Fokker–Planck equation as [56, 57, 28]

$\frac{\partial\eta}{\partial
t}=\frac{\partial}{\partial\omega}D_{1}(\omega)\eta+\frac{1}{2}\frac{\partial^{2}}{\partial\omega^{2}}D_{2}(\omega)\eta,$
(7)

where $D_{1}(\omega)$ is known as the drift and $D_{2}(\omega)$ the diffusion.
In general, $D_{n}(\omega)$ are the Kramers–Moyal coefficients. If we compare
this with our Fokker–Planck equation, we see that

$D_{1}(\omega)=-\frac{d}{m}\omega,\quad D_{2}(\omega)=D_{2}=\frac{B}{m^{2}}.$
(8)

First, we notice the linear response of the drift $D_{1}(\omega)$, which is a
linear function of $\omega$. Second, we note that given out choice of noise
$\xi$ mediated by $B$, the diffusion $D_{2}(\omega)=D_{2}$ is not a function
of the angular frequency.

We are now interested in obtaining the drift $D_{1}$ and diffusion $D_{2}$
directly from the frequency data. We can make use of (Nadaraya–Watson) non-
parametric estimators [58, 59] to extract the functional form of $D_{1}$ and
$D_{2}$ from the data [60, 61, 62, 63, 53, 64]. For a time series $x(t)$, like
power-grid frequency recordings, the Kramers–Moyal coefficients can be
estimated using a kernel density estimation [61]

$\displaystyle D_{m}(x)$ $\displaystyle\sim\frac{1}{m!}\frac{1}{\Delta
t}\langle(x(t+\Delta t)-x(t))^{m}|x(t)=x\rangle$ (9)
$\displaystyle\sim\frac{1}{m!}\frac{1}{\Delta
t}\frac{1}{N}\sum_{i=1}^{N-1}K(x-x_{i})(x_{i+1}-x_{i})^{m},$

where $K(x)$ is a normalised kernel with a given bandwidth $h$ (similar to
what the number of bins is in a histogram), such that $K(x)=1/hK(x/h)$. We use
the commonly employed Epanechnikov kernel $K(x)=3/4(1-x^{2})$, with support
$|x|<1$ [65]. Furthermore, we apply some finite-time corrections to improve
the quality of the estimations by accounting for the fact that the recordings
have a finite-time sampling (cf. Refs. [60, 64]).

In Fig. 3 we display the drift $D_{1}(f/\sigma(f))$ over the normalised power-
grid frequency $f$ (left) and the diffusion $D_{2}(f/\sigma(f))$, also
normalised (right). In the insets, the non-normalised versions of each
aforementioned Kramers–Moyal coefficient can be found.

Let us start by pointing out what immediately agrees with our understanding of
power-grid frequency. Both the Nordic Grid and the South African grid show a
clear linear response of the drift. This response almost overlaps, which
points to an identical ratio of the damping and inertia $d/m$, which is also
expected, as these go hand-in-hand. When we examine the Continental European
grid, we find a much stranger picture. Continental Europe has a band of
frequency where no control is active, i.e., there is no active control to
maintain the frequency at its nominal value, marked in the shaded area on the
left panel, known as the ‘deadband’ [66]. In the drift we even seem to see and
effect of ‘repulsion’, i.e., the drift has a positive slope, pushing the
frequency off from the nominal value. Moreover, at frequency deviations larger
than $|\sigma(f)|>1.5$, it also seems like there is no control or linear
response. This agrees to some extent with the control mechanism in Continental
Europe denoted as ‘secondary control’, which is a slow mechanism to restore
large frequency deviations to a nominal value. Being a slow mechanism, it
cannot be captured in our stochastic approach and functions as a slow
dynamical effect. We should note that the Nordic Grid and South African
systems similarly have their own ‘deadband’, as well as primary and secondary
control mechanisms, yet these are not detectable using our analysis.

Figure 4: Probability of occurrence of sags or overshoots in relation to the
recorded power-grid frequency in South Africa in 2021. (Left panel) Two-
dimensional probability density of the recorded frequency $f(t)$ and the
minimal $f_{\mathrm{MIN}}$ and maximal $f_{\mathrm{MAX}}$ recorded frequency.
(Right panel) Two-dimensional probability density of the recorded frequency
$f(t)$ and the difference $f_{\mathrm{MIN}}-f_{\mathrm{MAX}}$. When the
frequency is far from the nominal value of $50$ Hz, the minimal and maximal
frequencies come together, pointing to a strong control at play when the
system is far from equilibrium.

When we examine the diffusion $D_{2}(\omega)$ of the three recordings, seen on
the right panel of Fig. 3, we observe an almost constant value for both
Continental Europe and the Nordic Grid. The case of the South African grid is
a bit different. The diffusion is state-dependent, i.e., the value of
$D_{2}(\omega)$ is roughly proportional to $\omega^{2}$ and not a constant.
Thus, at first glance, it seems that when the South African power-grid
frequency is not at its nominal value of $50$ Hz but either slightly above or
slightly below, the amplitude of the noise is larger. We see this contrasts
with our initial simplistic assumption when writing (4), wherein we started
with simple state-independent or additive Gaussian white noise. This does not
seem to be the case in South Africa. Thus, we return to (7) and write a more
general static solution to the Fokker–Planck equation as

$\eta_{D_{2}}(\omega)\approx\sqrt{\frac{d}{2\pi
mD_{2}(\omega)}}e^{-\frac{d}{2mD_{2}(\omega)}\omega^{2}},$ (10)

wherein we have that the variance of $\eta(\omega)$ depends on the angular
frequency. This is the static solution to our swing equation with noise (1)
with multiplicative noise (or state-dependent noise). We note that this is
only approximate since this density is not always normalisable. In fact, with
$D_{1}(\omega)$ a linear function of $\omega$, the underlying equation of
motion is identical to an Ornstein–Uhlenbeck process, from which we know its
closed-form solutions [56, 57, 28]. Unfortunately for the simple case of
multiplicative/state-dependent noise $D_{2}(\omega)=b+c\omega^{2}$, with
$b,c>0$ two constants, that resemble our case in South Africa, one cannot
evaluate the expected values $\operatorname{E}[X^{n}]$ explicitely. We can
nevertheless approximate the density function for small $c>0$. Let us consider
the infinite Taylor expansion of $\eta(\omega)$ for small $c$, denoted
$\eta_{c}(\omega)$, given by

$\displaystyle\eta_{c}(\omega)$
$\displaystyle=\eta(\omega)\left[1+\frac{c}{2b}\omega^{2}\left(\frac{d\omega^{2}}{mb}-1\right)\right.$
(11)
$\displaystyle\qquad\qquad\quad+\left.\frac{1}{2}\frac{c^{2}}{4b^{2}}\omega^{4}\left(\frac{d^{2}\omega^{4}}{m^{2}b^{2}}-6\frac{d\omega^{2}}{mb}+3\right)+\dots\right]$
$\displaystyle=\eta(\omega)\sum\limits_{n=0}^{\infty}\frac{z^{2n}}{n!}\operatorname{He}_{2n}(z)\left(\\!\sqrt{\frac{mc}{2d}}\right)^{2n}\\!\\!\\!\\!\\!\\!,~{}~{}\mathrm{with}~{}z^{2}=\frac{d\omega^{2}}{mb},$

with $\operatorname{He}_{2n}(z)$ the Hermite polynomials. From the generating
function of the Hermite polynomials

$e^{\frac{1}{2}zt-t^{2}}=\sum\limits_{n=0}^{\infty}\frac{\operatorname{He}_{n}(z)}{n!}t^{n}$
(12)

we get

$e^{-\frac{t^{2}}{2}}\cosh(zt)=\sum\limits_{n=0}^{\infty}\frac{\operatorname{He}_{2n}(z)}{(2n)!}t^{2n},$
(13)

from which we draw our approximate expression by substituting
$t=\sqrt{\frac{mc}{2d}}z=\sqrt{\frac{c}{2b}}\omega$

$\eta_{c}(\omega)=\eta(\omega)e^{-\frac{c\omega^{2}}{4b}}\cosh\left(\frac{\sqrt{dc}}{\sqrt{2m}b}\omega^{2}\right).$
(14)

With this closed-form approximation of the probability density $\eta(\omega)$
for small $c$ we can now find the expected values $\operatorname{E}[X^{n}]$,
from which we obtain a kurtosis $\kappa_{X}^{c}$

$\kappa_{X}^{c}=3.2140,$ (15)

with $d/m=0.6723$, $b=0.0023$, and $c=0.0467$ estimated from the drift and
diffusion. This is close to our empirical value of $\kappa=3.33$. We should
note that this is a crude approximation, but it helps understand the effect of
multiplicative noise on generating leptokurtic distributions (i.e.,
distributions with kurtosis $\kappa>3$).

We will now contrast our analysis with the aid of two other power-grid
frequency recordings from South Africa.

###  III.4 Counter-intuitive frequency deviations

We now address an interesting – and quite a counter-intuitive effect – in
power-grid frequency recordings in South Africa. Alongside recordings of the
frequency in South Africa with a sampling of $15$ minutes, we have access to
the smallest and largest recorded frequency measurements within each
$15$-minute block. Therein we have three time series we can examine: The
frequency $f(t)$ we have thus far examined, the recorded maximal frequency
$f_{\mathrm{MAX}}(t)$, and the recorded minimal frequency
$f_{\mathrm{MIN}}(t)$.

In Fig. 4 we display the two-dimensional probability density functions of
$\rho(f,f_{\mathrm{MAX}})$ and $\rho(f,f_{\mathrm{MIN}})$ (left) and the
difference $f_{\mathrm{MAX}}-f_{\mathrm{MIN}}$ (right). If we focus on the
left panel, we can see that an increase of the frequency $f$ similarly
increases both the $f_{\mathrm{MIN}}$ and the $f_{\mathrm{MAX}}$. This is
somewhat counter-intuitive. If the frequency $f$ is above its nominal value of
$50$ Hz within a period of 15 minutes we expect that the maximal frequency
$f_{\mathrm{MAX}}$ follows this effect, but not necessarily the minimal
frequency $f_{\mathrm{MIN}}$. In the end, one would hope that the frequency is
brought down back to $50$ Hz in the span of $15$ minutes. This follows in
reverse if the frequency drops below $50$ Hz.

One potential explanation for this would assume that we do not have a
stationary process, but instead that the average (on the centre of the
distribution) is moving with the frequency. This supports the idea that what
large overshoots and nadirs of the frequency we observe in the South African
grid follow a similar phenomenon to those in Continental Europe and the Nordic
grid: These are dynamical effects of the market. This being the case, we would
expect that the minimal and maximal frequency would shift in tandem. If there
is too much generation in the system the minimal frequency is likely to also
have increased. Similarly, too much consumption brings down the maximal
frequency.

What remains to be explained is shown in the right panel of Fig. 4. Here we
see the difference between $f_{\mathrm{MAX}}$ and $f_{\mathrm{MIN}}$. We see
that when the system is out of equilibrium these two measured frequencies come
together. This points to a strong control mechanism starting up as the
frequency deviates from its nominal value. Yet this is in clear contradiction
with our previous result wherein we showed there is multiplicative noise in
the system. If indeed there is multiplicative noise at frequencies far from
nominal the distribution must widen and, statistically speaking,
$f_{\mathrm{MAX}}$ and $f_{\mathrm{MIN}}$ should grow further apart, not come
together.

At the present state, we cannot offer an explanation for this. The obvious
remark points to the low time sampling of the frequency data, which inevitably
leads to finite-time effects when estimating the drift and diffusion [60, 62,
64]. Possibly at higher sampling rates of the power-grid frequency this effect
subsides.

Altogether, we have offered a first explanation for the non-Gaussian
distribution of the power-grid frequency recordings in South Africa – with its
various caveats. Yet, as we have seen for Continental Europe and the Nordic
Grid, noise in power-grid frequency does not seem to be multiplicative/state-
dependent, but constant across the frequencies. Moreover, we have also
discussed that these recordings similarly have high kurtosis and that this
arises primarily due to changes in dispatch and market effects. In South
Africa, similar large deviations of power-grid frequency – which lead to high
kurtosis of the distribution – can occur due to scheduled and unscheduled
load-shedding and rolling blackout events. We now explore this avenue of
effects on a coarse-grain modular power-grid model of South Africa with
numerical solvers of the swing equation with noise (1).

###  III.5 The effects of rolling blackouts and load shedding

Figure 5: Model power-grid network of the South Africa Grid based on PyPSA-ZA
[67], with 26 nodes and 41 lines. The classification into generator and load
nodes is based on various available data [67, 68, 69] and a power-flow
balancing obtain with pandapower [70].

In order to simulate power-grid voltage frequency dynamics, we make use of a
coarse-grain model of the South African grid introduced by Hörsch & Calitz
[67]. The model comprises $26$ geographically distributed nodes, located
roughly at the cultural centres of the $9$ provinces of South Africa, coupled
via $41$ aggregated high-voltage transmission lines. A representation of the
model is given in Fig. 5. Each node in this power-grid model is considered as
modelled via the swing equation (1) diffusively coupled following (2). The
total capacity at each node is assessed via the freely available data from the
lists of power plants in South Africa found in Wikipedia [68]. The capacity of
each power plant is added to the nearest geographic node of the power-grid
model. To account for the consumption, the estimation of the population in
2021 in each region nearest to a node is also summed [69]. The total
electricity demand for 2021 was $219\,423$ GWh, following ESKOM’s 2021
Integrated Report [33], from which we estimate the total average power
consumption and split it according to each node’s population. Inertia $M_{i}$
and damping $D_{i}$ for each node $i$ is estimated following Pagnier & Jacquod
[71]. The transmission line parameters, i.e., the susceptance, are given in
Ref. [67].

To simulate the frequency (and rotor-angle) dynamics we first solve the
alternating current (AC) power-flow problem using pandapower [70]. It,
firstly, ensures the physical power-flow problem is solvable, but also
adequately balances generation and consumption at each node, resulting in the
generated power $P_{j}^{\mathrm{m}}$ at each node. Nodes with negative
generation, i.e., consumption, act as loads in the system. The sketch of which
nodes act as generator and which act and load is given in Fig. 5. For the
numerical simulations we assume a lossless system, i.e., $G_{j,\ell}=0$, a
fairly common assumption in high-voltage transmission grids. We consider as
well a small Gaussian white noise added to each node to account for some
stochasticity in the problem.

The target is now to model load-shedding events to showcase how these can
influence power-grid frequency dynamics. As we have little access to the
actual operation, implementation, phasing, and duration of the load-shedding
events, we focus instead on matching the few statistics we do have access: the
number of hours under different load-shedding stages in South Africa. For
reference, we showcase in Tab. 1 the years 2019, 2020, and 2021.

year |  St. 1 [GWh] |  St. 2 [GWh] |  St. 3 [GWh] |  St. 4 [GWh] |  St. 5 [GWh] |  St. 6 [GWh] |  # hours
---|---|---|---|---|---|---|---
2019 | 43 | 618 | 93 | 568 | – | 30 | 530
2020 | 133 | 1192 | 141 | 332 | – | – | 859
2021 | 79 | 1848 | 210 | 384 | – | – | 1169
Table 1: Total estimated load shed energy as reported by the Council for
Scientific and Industrial Research – Energy Centre of South Africa (see p. 171
in Ref. [10]) at different load-shedding stages (St. 1 to St. 6). Shown as
well are the total number of hours under load-shedding in South Africa in
2019, 2020, and 2021.

To model the load-shedding events, we consider an abrupt change in demand by
randomly selecting a node and a random subset of its neighbours. The total
load-shedding is split uniformly across these nodes and the power-flow balance
is recalculated. To mimic the system close to reality, we simulate a similar
ratio of load-shedding events covering $\sim\\!\\!13\%$ of time in a
simulation lasting 336 hours. The change in power leads to large frequency
excursions – one of the causes for the elevated kurtosis of the frequency
statistics. In Fig. 6 we display the statistics of the modelled system as
measured at the node located in Bloemfontein (see Fig. 5). We can observe a
direct impact on the kurtosis as the model system adjusts to the change in the
power flows due to the sudden events of load shedding. This leads to large
frequency excursions and thus elevated kurtosis. We note that in our
simulations the system never loses power balance, simply readjusted generation
directly as demand drops. The dynamical simulations simply readjust the angle
and frequency dynamic variables to these changes.

One can easily argue that most likely load shedding leads to a temporarily
unbalanced system – which would imply the dynamical system would attain a new
fixed point not at $f=50\,$Hz (or $\omega=0$). This, arguably, would explain
the data better – yet numerically simulating unbalanced systems results in the
loss of stable fixed points and the uncoupling of some of the nodes in the
system. Thus, we pursue a simpler approach and consider an abrupt change of
flows, which inherently maintain power balance (i.e., the system converges
back to $\omega_{i}(t)=0$ after some time) but affects the angles
$\theta_{i}$, leading to excursions of the angular frequency. Overall, as
expected, the statistics of the frequency show an increased kurtosis
$\kappa_{\mathrm{ZA}}^{\mathrm{sim}}=3.76$, which, in our limited simulations,
presents itself above the empirical kurtosis $\kappa_{\mathrm{ZA}}=3.33$. One
of the main reasons for such an elevated kurtosis in the numerical simulation
is that this leads to large angular frequency excursions in both positive and
negative directions.

Figure 6: Numerical simulation of power-grid frequency obeying the swing-
equation model (1) with 26 machine/nodes and 41 transmission lines/edges. The
simulated system is affected by load-shedding events 39 out of 283 hours
($\sim 13.8\%$ of the time). We can observe a direct effect on the kurtosis
$\kappa$ as large excursions of the frequency take place as the power flow
readjust to a rapid change of the load. The data is truncated at $\pm~{}6$
standard deviations as various control mechanisms would enter into frequency
control that cannot be easily simulated with (1).

##  IV Conclusion

In this work we have examined power-grid frequency recordings from the South
African grid for the year of 2021, having a 15-minute temporal resolution. We
contrasted these recordings with similar recordings from Continental Europe
and the Nordic Grid, both having a much higher temporal resolution of 1
second. At this much coarser time resolution, we observe very similar
statistics of power-grid frequency in South Africa as in Europe. All
distributions have heavy tails (kurtosis $\kappa>3$) and are slightly skewed,
with South Africa having a kurtosis of $\kappa_{\mathrm{ZA}}=3.33$.

Using a simplistic Fokker–Planck equation to describe the statistics of power-
grid frequency, we examine the drift and diffusion of all three recordings via
a non-parametric estimation of the Kramers–Moyal coefficients. Overall, we
find an identical response mechanism in the recordings, which is given by the
ratio of overall damping and inertia in the system. That is to say, for all
recordings we estimate an identical linear-restoring-force mechanics of the
frequency pushing the system to the nominal value of $50$ Hz. Examining the
diffusion in the recordings we come across what resembles
multiplicative/state-dependent noise in the South African recordings. We
derive an approximate analytical expression to treat multiplicative noise and
show that with such an expression one can estimate the statistical moments
$\operatorname{E}[X^{n}]$. From this expression, we obtain our first
explanation for the large kurtosis, resulting in
$\kappa_{\mathrm{ZA}}^{c}=3.21$.

This, however, disagrees with a subsequent analysis of two additional South
Africa power-grid time series we include in this work. Having access to both
the minimal and maximal measure power-grid frequency in each 15-minute block
of 2021, we find a counter-intuitive relation of their statistics when the
frequency is far from its nominal value. Away from equilibrium, the minimal
and maximal power-grid frequency recordings become closer to each other – this
disagree with our Fokker–Plank approach and contrasts with the presence of
multiplicative noise, which should drive the maximal and minimal frequency
apart, not together.

The narrowing of the minimal and maximal frequency, along with other details,
points to a simpler explanation for the non-Gaussian statistics of power-grid
frequency in South Africa: The detrimental effect of unscheduled load
shedding. Interestingly, large frequency deviations in Europe are explained by
various dispatch changes and electricity market activities (cf. Ref. [19,
54]), which, from a point-of-view solely of power-grid frequency, means large
deviations are induced by strong changes in generation and consumption – much
alike load shedding events. To this effect, we simulate a model power-grid
system of South Africa composed of 26 nodes and 41 transmission lines. We
consider a losses system of inertial masses coupled in the network obeying a
swing equation, wherein generation and consumption are drawn from freely
available data of the power grid and the population is South Africa. In order
to mimic the effects of load-shedding we simulate the system as affected by
drops in the load for roughly $13\%$ of the time. We consider the effects at a
central location in South Africa and show that the simulation lead as well to
an increased kurtosis of $\kappa_{\mathrm{ZA}}^{\mathrm{sim}}=3.76$.

Overall we provide a explanation for the statistics of power-grid frequency in
South Africa – which shows similar characteristics to more studied power grids
in Europe – yet this grid has contended (and still contends) with various
load-shedding events, which invariantly affects frequency stability. We
provide two explanations for the presence of elevated kurtosis of the
frequency statistics in South Africa. We hope this work can bring into focus
hitherto overlooked power grids around the world, particularly those that have
been or are afflicted by serious infrastructural, economical, or societal
problems.

###### Acknowledgements.

We thank Pedro G. Lind and Philipp C. Böttcher for their early suggestions and
remarks on the manuscript.

## References

  * United Nations General Assembly [2015] United Nations General Assembly, Transforming our World: the 2030 Agenda for Sustainable Development (2015), a/RES/70/1.
  * Thopil  _et al._ [2018] M. S. Thopil, R. C. Bansal, L. Zhang, and G. Sharma, A review of grid connected distributed generation using renewable energy sources in South Africa, Energy Strategy Reviews 21, 88 (2018).
  * Bohlmann and Inglesi-Lotz [2018] J. A. Bohlmann and R. Inglesi-Lotz, Analysing the South African residential sector’s energy profile, Renewable and Sustainable Energy Reviews 96, 240 (2018).
  * Lawrence [2020] A. Lawrence, Energy decentralization in South Africa: Why past failure points to future success, Renewable and Sustainable Energy Reviews 120, 109659 (2020).
  * Ayamolowo  _et al._ [2022] O. J. Ayamolowo, P. Manditereza, and K. Kusakana, South Africa power reforms: The path to a dominant renewable energy-sourced grid, Energy Reports 8, 1208 (2022).
  * Meyer and Overen [2021] E. L. Meyer and O. K. Overen, Towards a sustainable rural electrification scheme in South Africa: Analysis of the status quo, Energy Reports 7, 4273 (2021).
  * Jain and Jain [2017] S. Jain and P. K. Jain, The rise of renewable energy implementation in South Africa, Energy Procedia 143, 721 (2017).
  * Baker and Phillips [2019] L. Baker and J. Phillips, Tensions in the transition: The politics of electricity distribution in South Africa, Environment and Planning C: Politics and Space 37, 177 (2019).
  * Buraimoh  _et al._ [2020] E. Buraimoh, A. A. Adebiyi, O. J. Ayamolowo, and I. E. Davidson, South Africa electricity supply system: The past, present and the future, in  _2020 IEEE PES/IAS PowerAfrica_ (2020) pp. 1–5.
  * Pierce and Ferreira [2022] W. Pierce and B. A. Ferreira, Statistics of utility-scale power generation in South Africa in 2021 (2022), presentation at the CSIR Energy Centre, April 2022.
  * Gehringer  _et al._ [2018] C. Gehringer, H. Rode, and M. Schomaker, The effect of electrical load shedding on pediatric hospital admissions in South Africa, Epidemiology 29, 841 (2018).
  * Laher  _et al._ [2019] A. Laher, B. van Aardt, A. Craythorne, M. van Welie, D. Malinga, and S. Madi, ‘Getting out of the dark’: Implications of load shedding on healthcare in South Africa and strategies to enhance preparedness, South African Medical Journal 109, 899 (2019).
  * Aidoo and Briggs [2019] K. Aidoo and R. C. Briggs, Underpowered: Rolling blackouts in Africa disproportionately hurt the poor, African Studies Review 62, 112 (2019).
  * Louw [2019] Q. E. Louw, The impact of non-technical losses: A South African perspective compared to global trends (2019).
  * Khonjelwayo and Nthakheni [2021] B. Khonjelwayo and T. Nthakheni, Determining the causes of electricity losses and the role of management in curbing them: A case study of city of Tshwane metropolitan municipality, South Africa, Journal of Energy in Southern Africa 32, 45 (2021).
  * Kundur [1994] P. S. Kundur,  _Power System Stability and Control_ , 1st ed. (McGraw-Hill, New York, 1994).
  * Machowski  _et al._ [2020] J. Machowski, Z. Lubosny, J. W. Bialek, and J. R. Bumby,  _Power System Dynamics: Stability and Control_ , 3rd ed. (John Wiley & Sons, New York, 2020).
  * Pagani and Aiello [2013] G. A. Pagani and M. Aiello, The power grid as a complex network: A survey, Physica A: Statistical Mechanics and its Applications 392, 2688 (2013).
  * Schäfer  _et al._ [2018] B. Schäfer, C. Beck, K. Aihara, D. Witthaut, and M. Timme, Non-Gaussian power grid frequency fluctuations characterized by Lévy-stable laws and superstatistics, Nature Energy 3, 119 (2018).
  * Anvari  _et al._ [2016] 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 18, 063027 (2016).
  * Schmietendorf  _et al._ [2017] 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 90, 222 (2017).
  * Wolff  _et al._ [2019] M. F. Wolff, K. Schmietendorf, P. G. Lind, O. Kamps, J. Peinke, and P. Maass, Heterogeneities in electricity grids strongly enhance non-Gaussian features of frequency fluctuations under stochastic power input, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 103149 (2019).
  * Qin  _et al._ [2022] H. Qin, Y. Cao, H. Zhang, Z. Gao, D. Yang, and H. Ma, Quantitative analysis of frequency dynamics in real power system, in  _2022 7th Asia Conference on Power and Electrical Engineering (ACPEE)_ (2022) pp. 1091–1095.
  * Rydin Gorjão  _et al._ [2021a] L. Rydin Gorjão, B. Schäfer, D. Witthaut, and C. Beck, Spatio-temporal complexity of power-grid frequency fluctuations, New Journal of Physics 23, 073016 (2021a).
  * Tyloo [2022] M. Tyloo, Layered complex networks as fluctuation amplifiers, Journal of Physics: Complexity 3, 03LT01 (2022).
  * Han  _et al._ [2022a] C. Han, D. Witthaut, L. Rydin Gorjão, and Ph. C. Böttcher, Collective effects and synchronization of demand in real-time demand response, Journal of Physics: Complexity 3, 025002 (2022a).
  * Friedrich  _et al._ [2011] R. Friedrich, J. Peinke, M. Sahimi, and M. R. R. Tabar, Approaching complexity by stochastic methods: From biological systems to turbulence, Physics Reports 506, 87 (2011).
  * Tabar [2019] M. R. R. Tabar,  _Analysis and Data-Based Reconstruction of Complex Nonlinear Dynamical Systems_ (Springer International Publishing, 2019).
  * Acebrón and Spigler [1998] J. A. Acebrón and R. Spigler, Adaptive frequency model for phase-frequency synchronization in large populations of globally coupled nonlinear oscillators, Physical Review Letters 81, 2229 (1998).
  * Acebrón  _et al._ [2000] J. A. Acebrón, L. L. Bonilla, and R. Spigler, Synchronization in populations of globally coupled oscillators with inertial effects, Physical Review E 62, 3437 (2000).
  * Filatrella  _et al._ [2008] G. Filatrella, A. H. Nielsen, and N. F. Pedersen, Analysis of a power grid using a Kuramoto-like model, The European Physical Journal B-Condensed Matter and Complex Systems 61, 485 (2008).
  * Rodrigues  _et al._ [2016] F. A. Rodrigues, T. K. D. Peron, P. Ji, and J. Kurths, The Kuramoto model in complex networks, Physics Reports 610, 1 (2016).
  * ESKOM [2021] ESKOM, ESKOM: Integrated Report 2021 (2021), for the Financial, Sustainability, and the Integrated reports, see https://www.eskom.co.za/investors/integrated-results/.
  * Motjoadi  _et al._ [2020] V. Motjoadi, P. Bokoro, and M. Onibonoje, A review of microgrid-based approach to rural electrification in South Africa: Architecture and policy framework, Energies 13, 2193 (2020).
  * Duffaut Espinosa  _et al._ [2018] L. A. Duffaut Espinosa, M. Almassalkhi, P. Hines, and J. Frolik, System properties of packetized energy management for aggregated diverse resources, in  _2018 Power Systems Computation Conference (PSCC)_ (2018) pp. 1–7.
  * Liao  _et al._ [2015] S. Liao, J. Xu, Y. Sun, W. Gao, L. Xu, L. Huang, X. Li, J. Gu, and J. Dong, WAMS-based frequency regulation strategy for photovoltaic system in isolated power systems, in  _2015 IEEE Power & Energy Society General Meeting_ (2015) pp. 1–5.
  * Schmietendorf  _et al._ [2014] K. Schmietendorf, J. Peinke, R. Friedrich, and O. Kamps, Self-organized synchronization and voltage stability in networks of synchronous machines, The European Physical Journal Special Topics 223, 2577 (2014).
  * Sharafutdinov  _et al._ [2018] K. Sharafutdinov, L. Rydin Gorjão, M. Matthiae, T. Faulwasser, and D. Witthaut, Rotor-angle versus voltage instability in the third-order model for synchronous generators, Chaos: An Interdisciplinary Journal of Nonlinear Science 28, 033117 (2018).
  * Böttcher  _et al._ [2022] P. C. Böttcher, D. Witthaut, and L. Rydin Gorjão, Dynamic stability of electric power grids: Tracking the interplay of the network structure, transmission losses, and voltage dynamics, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 053117 (2022).
  * Witthaut  _et al._ [2022] 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 94, 015005 (2022).
  * Newman [2018] M. Newman,  _Networks: An Introduction_ , 2nd ed. (Oxford University Press, 2018).
  * Rydin Gorjão  _et al._ [2020a] L. Rydin Gorjão, R. Jumar, H. Maass, V. Hagenmeyer, G. C. Yalcin, J. Kruse, M. Timme, C. Beck, D. Witthaut, and B. Schäfer, Open database analysis of scaling and spatio-temporal properties of power grid frequencies, Nature Communications 11, 6362 (2020a).
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2D Graphics Environment, Computing in Science Engineering 9, 90 (2007).
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. Fernández del Río, M. Wiebe, P. Peterson, Pearu andGérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature 585, 357 (2020).
  * Virtanen  _et al._ [2020] 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, İ. 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, P. van Mulbregt, and SciPy 1.0 Contributors, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17, 261 (2020).
  * The pandas development team [2020] The pandas development team, pandas-dev/pandas: Pandas (2020), Zenodo 3509134.
  * [47] M. Tyloo, J. Hindes, and Ph. Jacquod, Finite-time correlations boost large voltage-angle fluctuations in electric power grids, ArXiv 2203.00590 (2022).
  * Salcedo-Sanz  _et al._ [2022] S. Salcedo-Sanz, D. Casillas-Pérez, J. Del Ser, C. Casanova-Mateo, L. Cuadra, M. Piles, and G. Camps-Valls, Persistence in complex systems, Physics Reports 957, 1 (2022).
  * Kraljic [2022] D. Kraljic, Towards realistic statistical models of the grid frequency, IEEE Transactions on Power Systems (2022), (Early Access).
  * Han  _et al._ [2022b] C. Han, H. Hilger, E. Mix, P. C. Böttcher, M. Reyers, C. Beck, D. Witthaut, and L. Rydin Gorjão, Complexity and persistence of price time series of the european electricity spot market, PRX Energy 1, 013002 (2022b).
  * Rydin Gorjão  _et al._ [2022] L. Rydin Gorjão, L. Vanfretti, D. Witthaut, C. Beck, and B. Schäfer, Phase and amplitude synchronization in power-grid frequency fluctuations in the Nordic Grid, IEEE Access 10, 18065 (2022).
  * Pagnier  _et al._ [2022] L. Pagnier, J. Fritzsch, Ph. Jacquod, and M. Chertkov, Toward model reduction for power system transients with physics-informed PDE, IEEE Access 10, 65118 (2022).
  * Rydin Gorjão and Meirinhos [2019] L. Rydin Gorjão and F. Meirinhos, kramersmoyal: Kramers–Moyal coefficients for stochastic processes, Journal of Open Source Software 4, 1693 (2019).
  * Rydin Gorjão  _et al._ [2020b] L. Rydin Gorjão, M. Anvari, H. Kantz, C. Beck, D. Witthaut, M. Timme, and B. Schäfer, Data-driven model of the power-grid frequency dynamics, IEEE Access 8, 43082 (2020b).
  * Farmer and Rix [2021] W. J. Farmer and A. J. Rix, Understanding the impact of network topology on frequency stability considering continuous spatial-temporal disturbances from wind generation, International Journal of Electrical Power & Energy Systems 129, 106776 (2021).
  * Risken and Frank [1996] H. Risken and T. Frank,  _The Fokker–Planck Equation_, 2nd ed. (Springer-Verlag, Berlin, Heidelberg, 1996).
  * Gardiner [2009] C. Gardiner,  _Stochastic Methods: A Handbook for the Natural and Social Sciences_ , 4th ed. (Springer-Verlag Berlin Heidelberg, 2009).
  * Nadaraya [1964] E. A. Nadaraya, On estimating regression, Theory of Probability & Its Applications 9, 141 (1964).
  * Watson [1964] G. S. Watson, Smooth regression analysis, Sankhyā: The Indian Journal of Statistics, Series A 26, 359 (1964).
  * Gottschall and Peinke [2008] J. Gottschall and J. Peinke, On the definition and handling of different drift and diffusion estimates, New Journal of Physics 10, 083034 (2008).
  * Lamouroux and Lehnertz [2009] D. Lamouroux and K. Lehnertz, Kernel-based regression of drift and diffusion coefficients of stochastic processes, Physics Letters A 373, 3507 (2009).
  * Anteneodo and Riera [2009] C. Anteneodo and R. Riera, Arbitrary-order corrections for finite-time drift and diffusion coefficients, Physical Review E 80, 031103 (2009).
  * Anteneodo and Queirós [2010] C. Anteneodo and S. M. D. Queirós, Low-sampling-rate Kramers–Moyal coefficients, Physical Review E 82, 041122 (2010).
  * Rydin Gorjão  _et al._ [2021b] L. Rydin Gorjão, D. Witthaut, K. Lehnertz, and P. G. Lind, Arbitrary-order finite-time corrections for the Kramers–Moyal operator, Entropy 23, 517 (2021b).
  * Epanechnikov [1967] V. A. Epanechnikov, Non-parametric estimation of a multivariate probability density, Theory of Probability & Its Applications 14, 153 (1967).
  * Vorobev  _et al._ [2019] P. Vorobev, D. M. Greenwood, J. H. Bell, J. W. Bialek, P. C. Taylor, and K. Turitsyn, Deadbands, droop, and inertia impact on power system frequency distribution, IEEE Transactions on Power Systems 34, 3098 (2019).
  * [67] J. Hörsch and J. Calitz, PyPSA-ZA: Investment and operation co-optimization of integrating wind and solar in South Africa at high spatial and temporal detail, ArXiv 1710.11199 (2017).
  * Wikipedia [2022] Wikipedia, List of power stations in South Africa (2022), last accessed 30th May 2022\.
  * [69] WorldPop: Spatial distribution of the population of South Africa in 2020 (unconstrained, from 100m resolution cell-grid), WorldPop – School of Geography and Environmental Science, University of Southampton; Department of Geography and Geosciences, University of Louisville; Departement de Geographie, Universite de Namur; Center for International Earth Science Information Network (CIESIN), Columbia University (2018). Global High Resolution Population Denominators Project - Funded by The Bill and Melinda Gates Foundation (OPP1134076).
  * Thurner  _et al._ [2018] L. Thurner, A. Scheidler, F. Schäfer, J.-H. Menke, J. Dollichon, F. Meier, S. Meinecke, and M. Braun, pandapower–An open-source python tool for convenient modeling, analysis, and optimization of electric power systems, IEEE Transactions on Power Systems 33, 6510 (2018).
  * Pagnier and Jacquod [2019] L. Pagnier and Ph. Jacquod, Inertia location and slow network modes determine disturbance propagation in large-scale power grids, PLOS ONE 14, 1 (2019).