# Microscopic Fluctuations in Power-Grid Frequency Recordings at the

Authors: Benjamin Schäfer; Leonardo Rydin Gorjão; G. Cigdem Yalcin; Ellen Förstner; Richard Jumar; Heiko Maass; Uwe Kühnapfel; Veit Hagenmeyer
Journal: Complexity
 (2023)

---

## Abstract

Complex systems, such as the power grid, are essential for our daily lives. Many
complex systems display multifractal behavior, correlated fluctuations and power
laws. Whether the power-grid frequency, an indicator of the balance of supply
and demand in the electricity grid, also displays such complexity remains a
mostly open question. Within the present article, we utilize highly resolved
measurements to quantify the properties of the power-grid frequency, making
three key contributions: First, we demonstrate the existence of power laws in
power-grid frequency measurements. Second, we show that below one second, the
dynamics may fundamentally change, including a suddenly increasing power
spectral density, emergence of multifractality and a change of correlation
behavior. Third, we provide a simplified stochastic model involving positively
correlated noise to reproduce the observed dynamics, possibly linked to
frequency-dependent loads. Finally, we stress the need for high-quality
measurements and discuss how we obtained the data analyzed here.


---

## Body

## 1. Introduction

### 1.1. Motivation

Complex systems, sometimes natural, sometimes man-made, surround us at all scales
as an integral part of everyday life, ranging from neural networks at small
scales to global pandemics or the World Wide Web at large scales [1]. To satisfy the most basic needs of
modern people in everyday life, we rely on power systems, which are an important
example of a complex system [2, 3]. This dependence ranges from the simple
cooling of food, vital medical care, long-distance transportation, instant
communication, and industrial automation to advanced scientific experiments,
many of which are impossible without electricity. Since reliability and
stability are essential for our power grid, analysis of power systems as complex
systems has become critically important [4].A distinguishing feature of complex systems is the emergence of new phenomena and
complex dynamics based on the various nonlinear interactions of a large
collection of components. One particularly interesting and universal property is
the observation of power laws [5]. Under
power laws, we understand that a function decays as fx∼x−β, for large x, with some constant β>0. Such power laws have been investigated in real-world data
sets from a wide variety of fields, including biological and technical systems
[6]. Empirically, they are often
discussed as the emergence of heavy or fat tails, e.g., web nodes with an
extremely large number of connections exist more often than intuitively
expected. From a theoretical point of view, the scale-invariance and
universality of power laws explain their popularity in many models [7]. Furthermore, the fractal dimension of a
system can be used to investigate dynamic processes in complex systems and is
indeed one of the most basic and important quantities to characterize a system
[8]. It captures how a system with
many degrees of freedom effectively only occupies a fraction of the
high-dimensional state space.In the present paper, we focus on power grids, i.e., electrical supply networks
and their dynamics as a complex system. While the network itself easily
constitutes a complex system [4], less
attention has been paid to individual local dynamics [9]. One key quantity for control and monitoring purposes of
the power system is the power-grid frequency. It follows the balance of
generation and consumption and allows insights into the general stability of the
grid as well as the impact of renewable generators [10–12]. So
far, only a few studies have investigated power laws or even multifractality in
power grids [13], partially because data
are not easily available. Consistently, frequency fluctuations and their
fractality on a short time scale are not fully understood.


### 1.2. Literature Review

The dynamics of the power-grid frequency are hard to describe analytically and
even numerically due to the nonlinear interactions of power generators and
consumers, the complex network topology, the large number of individuals as well
as companies, countries, and added stochastic effects, e.g., stemming from
renewable generation [9, 14]. This stresses the need to acquire
high-quality and trustworthy data on power systems as an approach to
quantitatively describe and understand the power grid as a complex system.
Unfortunately, reliable, high-quality continental, and global measurements from
grid frequency are only rarely made publicly available [15].To obtain high-quality measurements, initiatives such as grid radar [16] and GridEye/FNET [17] have been developed. From the academic side, the
electrical data recorder (EDR) has been developed as a device for
high-resolution time series acquisition in low-voltage distribution grids. It
provides UTC-time-stamped frequency estimates up to ten times per second and
allows long-term storage of full voltage and/or current waveforms for further
analysis [18–21]. The frequency is estimated using the well-known
zero-crossing technique that enables straightforward and fast processing and is
widely developed by scientists [22–25].


### 1.3. Paper Contribution and Organization

Within the present article, we make three key contributions: first, we observe
power-law decay in the power spectral density (Section2.1.1). Second, we show how the frequency displays
complex behavior on time scales below one second. In particular, we observe an
increasing spectral density, emerging multifractality and a switch from
correlated to anticorrelated time series. We further highlight differences in
short-term fluctuations among various power grids (Section 2.1.2). Third, we propose a simple stochastic process
that reproduces key properties of the empirical spectrum (Section 2.1.3). Next, we emphasize the need for
trustworthy, highly resolved measurements of the power-grid frequency and show
one way how these can be obtained (Section 2.2). We close with a discussion of our results (Section 3). Note that an earlier version of this
article was published as a preprint [26].


## 1.1. Motivation

Complex systems, sometimes natural, sometimes man-made, surround us at all scales
as an integral part of everyday life, ranging from neural networks at small
scales to global pandemics or the World Wide Web at large scales [1]. To satisfy the most basic needs of
modern people in everyday life, we rely on power systems, which are an important
example of a complex system [2, 3]. This dependence ranges from the simple
cooling of food, vital medical care, long-distance transportation, instant
communication, and industrial automation to advanced scientific experiments,
many of which are impossible without electricity. Since reliability and
stability are essential for our power grid, analysis of power systems as complex
systems has become critically important [4].A distinguishing feature of complex systems is the emergence of new phenomena and
complex dynamics based on the various nonlinear interactions of a large
collection of components. One particularly interesting and universal property is
the observation of power laws [5]. Under
power laws, we understand that a function decays as fx∼x−β, for large x, with some constant β>0. Such power laws have been investigated in real-world data
sets from a wide variety of fields, including biological and technical systems
[6]. Empirically, they are often
discussed as the emergence of heavy or fat tails, e.g., web nodes with an
extremely large number of connections exist more often than intuitively
expected. From a theoretical point of view, the scale-invariance and
universality of power laws explain their popularity in many models [7]. Furthermore, the fractal dimension of a
system can be used to investigate dynamic processes in complex systems and is
indeed one of the most basic and important quantities to characterize a system
[8]. It captures how a system with
many degrees of freedom effectively only occupies a fraction of the
high-dimensional state space.In the present paper, we focus on power grids, i.e., electrical supply networks
and their dynamics as a complex system. While the network itself easily
constitutes a complex system [4], less
attention has been paid to individual local dynamics [9]. One key quantity for control and monitoring purposes of
the power system is the power-grid frequency. It follows the balance of
generation and consumption and allows insights into the general stability of the
grid as well as the impact of renewable generators [10–12]. So
far, only a few studies have investigated power laws or even multifractality in
power grids [13], partially because data
are not easily available. Consistently, frequency fluctuations and their
fractality on a short time scale are not fully understood.


## 1.2. Literature Review

The dynamics of the power-grid frequency are hard to describe analytically and
even numerically due to the nonlinear interactions of power generators and
consumers, the complex network topology, the large number of individuals as well
as companies, countries, and added stochastic effects, e.g., stemming from
renewable generation [9, 14]. This stresses the need to acquire
high-quality and trustworthy data on power systems as an approach to
quantitatively describe and understand the power grid as a complex system.
Unfortunately, reliable, high-quality continental, and global measurements from
grid frequency are only rarely made publicly available [15].To obtain high-quality measurements, initiatives such as grid radar [16] and GridEye/FNET [17] have been developed. From the academic side, the
electrical data recorder (EDR) has been developed as a device for
high-resolution time series acquisition in low-voltage distribution grids. It
provides UTC-time-stamped frequency estimates up to ten times per second and
allows long-term storage of full voltage and/or current waveforms for further
analysis [18–21]. The frequency is estimated using the well-known
zero-crossing technique that enables straightforward and fast processing and is
widely developed by scientists [22–25].


## 1.3. Paper Contribution and Organization

Within the present article, we make three key contributions: first, we observe
power-law decay in the power spectral density (Section2.1.1). Second, we show how the frequency displays
complex behavior on time scales below one second. In particular, we observe an
increasing spectral density, emerging multifractality and a switch from
correlated to anticorrelated time series. We further highlight differences in
short-term fluctuations among various power grids (Section 2.1.2). Third, we propose a simple stochastic process
that reproduces key properties of the empirical spectrum (Section 2.1.3). Next, we emphasize the need for
trustworthy, highly resolved measurements of the power-grid frequency and show
one way how these can be obtained (Section 2.2). We close with a discussion of our results (Section 3). Note that an earlier version of this
article was published as a preprint [26].


## 2. Results

To analyze empirical complex systems, such as the power-grid frequency dynamics, we
require data. Some measurements with the one-second resolution are readily available
from transmission system operators (TSOs) for analysis and comparison with
independent measurements [15]. Meanwhile,
time series with the higher temporal resolution is not easy to obtain. Still, this
higher temporal resolution is critical to uncover key properties of the underlying
system, see e.g., [27, 28] and Figure 1:
representative frequency trajectories on 1 s and on 100 ms
(0.1 s) resolutions clearly display different behavior, with notable
fluctuations and high volatility on the short time scale.Figure 1
Snippets of power-grid frequency trajectories. Power-grid frequency
trajectories from 2019-07-27 06:00 in Continental Europe, for
20 seconds, including 1-second resolution and 100 ms
resolution data. The smaller 100 ms resolution recordings fluctuate
strongly around the 1 s resolution data, sometimes with very large
jumps. The frequency dynamics are revealed as fractal for this exemplary
trajectory: if we “zoom in” to higher temporal resolutions, we
observe that the data is self-affine, i.e., the data displays
self-similarity in its fluctuations. In the inset, we focus on a 2-second
span, where we can observe both small fluctuations (in the 1st second) and
large fluctuations (in the 2nd second).We address observable changes in the spectrum towards short time scales as well as
the correlation and fractal behavior of the time series. To understand these
observations and ensure high data quality, we outline how the recordings have been
obtained. We, therefore, describe the EDR measurement process and focus on its
reliability.


### 2.1. Emerging Phenomena on Short-Time Scales

The following observations are based on newly available 100 ms frequency
estimates and particularly focus on phenomena emerging at the newly-accessible
time scale below one second. The acquisition of the high-resolution frequency
data is described later in Subsection2.2.As one key quantity, we compute the power spectral density (PSD)Sf, which allows us to analyze the power distribution of
different oscillations in the components of the frequency recorded from the
system in a stochastic process [29]. As
with many processes observed in nature, we also expect the power spectrum of
self-affine signals to decay following a power law Sf∼f−β [30, 31]. PSD is a common tool that allows us to
perform further analysis on data originating from stochastic processes, such as
data obtained from electroencephalogram (EEG) [32], as well as for power-grid frequency signal analysis, which we
are currently considering here [29].
Furthermore, we may also estimate other properties such as the fractality, i.e.,
self-similarity on various scales, or the correlation of the time series, which
are closely connected to the spectrum [30, 33, 34].


#### 2.1.1. Power Law Analysis

Computing the PSD for the power-grid frequency recordings reveals new
qualitative dynamics on the short time scale at essentially all measurement
locations, Figure2. We utilize
several independent recordings (see Figure 2(a) for a map) from the Continental European grid (Figures
2(c) and 2(d)) and from three different synchronous areas
Figure 2(b). Let us briefly remark
on two observations we do not discuss in detail.Figure 2
Power-grid dynamics change qualitatively at the 1-second time scale.
We utilize recordings from various locations in Europe (a) to then
plot the spectral densitySf of the power-grid frequency (b)–(d). This
analysis includes independent synchronous areas (b) as well as
Continental Europe data obtained from grid radar (c) and from the
EDR (d). Note that we plot Sf and also report the corresponding time scales t=1/f. While all recordings initially display an
approximately decaying power law β, several recordings display a rising spectrum with
slope γ for time scales below t∼1s.
(a)(b)(c)(d)First, on the longer time scale of∼3…10 seconds, we observe pronounced interarea oscillations,
where different regions within one synchronous area oscillate against each
other. We notice the spatial dependency in the spectrum as different
locations show different magnitudes and positions (in terms of time) of
these peaks. This is explained by the fact that, e.g., north-south and
east-west oscillations have different frequency; see also [15, 27, 35] for a more
detailed discussion of this phenomenon.Second, there are various narrow and distinct peaks on the time scale below
one second in the spectrum; see, e.g., Figure2(d) Karlsruhe and Oldenburg. Curiously, not all
recordings from the Continental European area display such peaks, e.g.,
neither Lisbon nor Istanbul. Other synchronous areas again show mixed
results: recordings from Sweden only show less pronounced peaks, while the
peaks are again present in the Faroe Island recordings. These peaks could
arise from characteristic eigenfrequencies of the power grid when subject to
a disturbance; see Figure 3. While
the short circuit with an active power step of approximately 440MW only leads to a frequency change of the order 10−6 (0.15 mHz) in a realistic grid simulation (see
Figure 3(a)), the dynamics take place
on the second to subsecond scale with eigenfrequencies of about
1.5 Hz. Note the 3 full oscillations between 2 and 4 seconds
in Figures 3(a) and 3(c). Further investigations will be
necessary to fully explore the origin and relevance of these peaks. In the
current study, we instead focus on another key observation: The qualitative
change in the spectrum around 1 second and the emergence of complex
phenomena at time scales <1 second (greyed out).Figure 3
Frequency and power dynamics in a transmission grid display
characteristic subsecond dynamics. We display the local frequencies
at 380/220/110-kV station DAXLA with 3 powerstation-blocks (a), the
real power of one block (b), as well as the absolute voltage
amplitude (c), and the absolute current on the line that displayed a
short circuit to trigger the event (d). The transient dynamics is
triggered via a three-phase short circuit between 0.5 and
0.65 seconds on a 110 kV line between Daxlanden and
Leopoldshafen (connection point of the KIT campus North). Line and
node parameters use realistic values for the South-West German power
grid. Frequency and voltage are given in the per-unit (p.u.) system,
i.e., at 1 p.u., they are at the reference value.On the long time scale oft>1 second, all measurements display a power law-like decay of
the spectral density Sf∼f−β∼tβ, with f the frequency and t=1/f the period. For ease of interpretation, we plot Sf as a function of f but also provide the corresponding 1/t values; see Figure 2. This decay follows approximately f−β, where we determine the exponent −β as the slope of the dashed line in the double-logarithmic
plot for each measurement site separately. As expected, the decay constant β is positive, indicating a decrease in the spectral density
for shorter periods (large frequencies) of the recordings. This decreasing
spectral density, i.e., the power law, has been observed in numerous
systems, ranging from many man-made systems, such as the World Wide Web,
over transport networks to biological systems [5, 6]. Power laws
remain an active field of research, emerging both on microscopic and
macroscopic scales. It often remains open to which extent the individual
units or their collective dynamics contribute to the observed power laws
[36].Meanwhile, the spectral density starts to increase on short time scales oft<1 second. We again characterize the slope of the spectrum
via a linear fit in the double-logarithmic plot and thereby determine the
slope γ, which in this case is smaller than zero for all
recordings from the Continental European power grid (Figures 2(c) and 2(d)). Interestingly, this negative γ is not exclusive to the Continental European region; see
Figure 2(b), where we show the power
spectral density of five other power-grid recordings, namely from the power
grids on the Faroe Islands, Russia (St. Petersburg), and Nordic (Sweden).
The measurements taken on the Faroe Islands and in St. Petersburg used again
the EDR, which was also used to make the recordings in Figure 2(a). In addition, we also include
three recordings provided by the Swedish TSO in three undisclosed distinct
locations in Sweden for the month of January 2020. These locations are
reported simply as north Sweden, middle
Sweden, and south Sweden, as provided by the
Swedish TSO.Let us systematically investigate the different slopesβ and γ of the spectral decays observed in Figure 2 as a scatter plot in Figure 4, see also Supplementary Material for
numerical details. First, we note that all measurements consistently report
a decaying power law with exponent β∼2.5. However, the dynamics at short time scales, given by γ is more varied γ∼−2…2. We identify three regimes: the Swedish recordings display
a continuous decay of the spectrum. Next, we observe several recordings with
small absolute values of γ, implying the spectrum becomes almost flat. Finally, many
recordings, mostly from the Continental European grid, display negative γ values and have an increasing spectrum, pointing towards a
complex process emerging at these time scales.Figure 4
Spectral measurements on short time scales are either positively,
negatively or not correlated. We plot estimates of the power law
exponentsβ (long time scale) and γ (short time scale).Given this surprising increase of the spectrum on short time scales, let us
elaborate on why this is almost certainly not a problem of the measurement
system but a genuine effect of the underlying power system. Subsection2.2. deals with the EDR data
acquisition and analyzes it towards reliability and trust. Additionally, we
make the following observations: First, the EDR measurements within one grid
are consistent: We clearly observe negative γ values in Continental Europe at several independent
measurement points, namely in Istanbul, Oldenburg, Karlsruhe, and Lisbon.
Second, the negative γ values are further confirmed for Continental Europe by
measurements from grid radar (Figure 2(c)), relying on a completely independent hardware and software
solution. Third, the EDR measurements show γ≈0 in some other synchronous areas, such as in Russia (St.
Petersburg) and on the Faroe Islands.Interestingly, the time series provided by the Swedish TSO has a positiveγ, while the EDR-recorded time series in Stockholm (recorded
during a different time period) yields a low negative value of γ. As we do not know the reason for this discrepancy, we
list some potential causes: the measurements by the EDR and the Swedish TSO
were taken at different times and could hence measure different states of
the power system. Furthermore, the EDR was connected to a power plug in a
hotel, i.e., to the low-voltage distribution grid used for residential
areas. Meanwhile, the Swedish TSO has access to measurements at power plants
or directly at the high-voltage transmission grid. While the frequency is
typically a global property, the negative γ values could possibly arise from the complex dynamics of
components at the local distribution layer, such as inverters and nonlinear
loads. Finally, we do not have access to the raw data recorded by the
Swedish TSO and hence cannot infer anything about their measurement setup,
or filtering while we discuss the limits of our measurements in detail.


#### 2.1.2. Multifractality and Correlations

In our sample trajectories in Figure1, we already noted a certain multifractal structure to the
power-grid frequency recordings. One natural question arising is: what is
the exact nature of the fractal or multifractal structure that underlies
power-grid frequency?Power-grid systems are complex dynamical systems affected and impacted by
various endogenous and exogenous pressures from both within (control
mechanics, thermal effects, and oscillations) and from outside (changes in
generation and consumption, weather, and failures). Every impactful element
acts at various different time scales, and these time scales range from
years (ever-growing generation), to months (weather cycles), to weeks
(workweek-weekend cycles), to hours (market and dispatch changes), to
minutes (generation adjustments and automated control), to seconds (dynamic
interference, travelling waves, and occasional failures), to milliseconds
(electrodynamic effects) [10, 37].To best describe the power-grid frequency on very short time scales, we model
the time series as a combination of stochastic processes. To that end, we
focus on the incremental structure of our recordingsft, given by(1)Δfr=ft+r−ft,where r is our incremental lag. From a stochastic process
point-of-view it is usually denoted “increments.” These
increments carry information relating directly to the process’
fastest time scales, i.e., a process changes on a time scale of order r, as well as information relating to the existence of
memory and fractality in the recordings.One of the most common characteristics of the incremental time seriesΔfr is the “scaling” of the increments in
relation to the incremental lag r. This is commonly encapsulated in the structure function SqΔfr,r(2)SqΔfr,r=Δfrq∼rqhq,where hq is the generalized Hurst exponent, also known as the
Hölder exponent, which captures the change in the (noncentered and
nonnormalized) statistical moments, where q is the power of the statistical moment. hq is related to the scaling exponent τq, which is often used in turbulence analysis, as follows: τq=qhq−1 [38].We utilize the generalized Hurst exponent for two purposes: to investigate
correlations and to explore multifractality. For an uncorrelated Brownian
motion, we haveh2=0.5, while a correlated process has h2>0.5 and an anticorrelated process h2<0.5. Furthermore, a monofractal process is one where the
function hq is essentially constant, i.e., independent of q, while a multifractal process is characterized by a q-dependent generalized Hurst exponent hq.If the original time series yields aq-dependent hq, i.e., a multifractal process, this could either be due to
different long-range correlations or due to a broad probability
distribution. To ascertain how this multifractality emerges, we shuffle the
power-grid frequency recordings, i.e., randomize the position of each time
series entry. We can then calculate the generalized Hurst exponents hshuffledq for the shuffled data. Since shuffling the data removes
correlations, time series that are not distributional
multifractal will have hq=h2=1/2, i.e., they will be uncorrelated at all scales. If, on the
contrary, we find that hq is nonconstant, this is an indication of distributional
multifractality.Computing the generalized Hurst exponent reveals again a complex structure of
the time series on the short time scale; see Figure5. On a time scale of 0.2 to 1.0 seconds, the
generalized Hurst exponent is clearly dependent on q, indicating a multifractal process coherent with earlier
observations (Figure 1). This
multifractality persists even when shuffling, indicating a distributional
multifractality. For longer time scales, particularly above
4 seconds, the time series becomes monofractal (constant hq). Interestingly, the correlation changes over the time
scales: On short time scales of 0.4 second to 1 second, the
power-grid frequency is anticorrelated as indicated by a small value of h2<0.5. For these time scales, essentially, any disturbance on
the grid away from the reference frequency of 50Hz or 60Hz is followed by an opposite effect. While we do not expect
control actions to play a major role [10] on this time scale, the effect is similar, any disturbance
is counteracted, leading to a frequency time series close to the reference
value. A possible explanation for the emerging anticorrelation could be the
self-regulating effect of frequency-dependent loads [40], i.e., a frequency drop would be met by a reduced
load and thereby counteract the frequency change. Meanwhile, the time series
becomes positively correlated on longer time scales of several seconds and
more, as indicated by a large value of h2>0.5. This persistent behavior aligns with our expectations of
a power system subject to load and generation ramps on time scales of
seconds to minutes and an overall inertial system.Figure 5
Frequency recordings are multifractal and anticorrelated on short and
correlated on longer time scales. We plot the generalized Hurst
exponenthq of power-grid frequency recordings from the Lisbon
recording in Continental Europe at different time scales, from top
to bottom: 0.4∼1 s, 1∼4 s, and 4∼10 s. The smallest time scales show a
nonconstant exponent hq, which points to the existence of a large number
of degrees of freedom in the system. As the time scales become
larger, the system gets progressively more ordered, as fluctuations
are quenched by the strong synchrony in the system. At a time scale
of 4∼10 s the shuffled grid frequency recording
becomes indistinguishable from noise hq=0.5 [39].Let us briefly discuss the role of distributional multifractality before then
utilizing the insights about correlations in the time series to reproduce
the observed spectrum. We do not know for certain why and how the
multifractality arises in the frequency recordings discussed. However, we do
know that one mechanism giving rise to a multifractal time series is the
superposition of several processes at a shorter time scale, formalized under
the concept of superstatistics by Beck and Cohen [41]. Assuming the observed frequency time series arises
from such superstatistics, it would offer further directions of research. In
superstatistics, a system changes from one state to another while we treat
all observations as if they originated from the same state. In terms of the
energy system, these different states of the system could be related to
different generation mixes and demand patterns or due to drivers such as
superstatistical behavior in the wind [42]. Note that superstatistics has been discussed on power-grid
frequency recordings before [11], and
there are further hints that the dynamics of power-grid recordings get more
complex on shorter time scales. In fact, if we consider the temporal scales
wherein we observe the strongest evidence of multifractality, i.e., for t<4 s, this agrees well with previous results [27, 28], wherein power-grid frequency increments are shown to be
well described by a convolution of Gaussian statistics with a parametric
function of varying variance. All of this points to even more research is
necessary to understand the complexity of power-grid frequency time series
on this short time scale.


#### 2.1.3. Stochastic Process

Let us return to the observation of an increasing spectrum and try to explain
this via physical processes on a short time scale. We have established from
our generalized Hurst exponent analysis that the frequency is positively
correlated on a longer time scale of several seconds and negatively
correlated on shorter time scales. For now, we consider simple mono-fractal
processes and explore how these can reproduce the observed spectrum. Indeed,
a very simple stochastic process can already reproduce such a spectrum,
namely combining a fractional, positively correlated
Ornstein–Uhlenbeck process with added, negatively correlated
fluctuations; see Figure6. Before
discussing Figure 6 in detail, let us
motivate the usage of an Ornstein–Uhlenbeck process:Figure 6
Reproducing the spectrum using an Ornstein-Uhlenbeck process with
added correlated short-term fluctuations.Following earlier work [11, 43], we may model the frequency
dynamics of the lossless high-voltage grid on the short time scale as follows:(3)ddtθi=ωi,Middtωi=Pi+σiξi−giωi+∑j=1NKijsinθj−θi,where we assume that each node i in the grid behaves approximately like a synchronous
machine whose state is characterized by its voltage phase angle θi and angular velocity ωi=2πfi−fR, which is proportional to the deviation of the frequency
from the reference frequency of fR=50Hz or fR=60Hz. The dynamics of the network are then determined by the
following parameters: inertia Mi, mechanical power Pi (positive for effective generators and negative for
effective consumers), fluctuations ξi with amplitude σi, damping gi, and the coupling matrix Kij.As we typically do not have precise information about the full network
topology, all nodal parameters and states, we want to simplify this to an
equation for the bulk angular velocityω-≔∑i=1NMiωi/∑i=1NMi by assuming a homogeneous damping to inertia ratio, c=gi/Mi [44] and balanced
average power ∑i=1NPi=0. We obtain the aggregated swing equation as follows:(4)ddtω-=−cω-+σ¯ξ¯t,with aggregated fluctuation amplitude σ¯ and aggregated fluctuations ξ¯, which depend on the precise nature of the noise [11]. The bulk frequency dynamics given
in equation (4) describes
the average behavior in one synchronous area, such as the Continental
European grid, neglecting specific local properties. Crucially, we have
simplified the nonlinear network dynamics to a relatively simple
Ornstein–Uhlenbeck process. One main difference to a standard
Ornstein–Uhlenbeck process is that the aggregated fluctuations ξ¯ are typically not well described by ordinary Brownian
motion, but instead, the effective fluctuations acting on the system are
either heavy-tailed itself or arise from a superposition of several
stochastic processes, so-called superstatistics [11, 41],
mentioned earlier.Regardless of the precise nature of the fluctuationsξ¯, for an Ornstein–Uhlenbeck process, we would expect
a spectrum following a simple power law (dotted line in Figure 6). In order to obtain an explicit
expression for the power-spectral density, recall that the energy of a time
series is given by(5)E≔∫−∞∞xt2dt.Parseval’s theorem dictates that(6)∫−∞∞xt2dt=∫−∞∞x∧f2df,with x∧f the Fourier transform of xt, i.e.,(7)x∧f=∫−∞∞xte−i2πftdt,with i the imaginary unit. Using Parseval’s theorem, we
see that the integrand of the energy in equation (5) is simply the spectral
density Sf given as follows:(8)Sf≔x∧f2.In the following, we are interested in a solution to the aggregated swing
equation (4) in Fourier
space. Defining ω-∧ and ξ¯∧ as the Fourier transforms of the bulk angular velocity ω- and the aggregated noise amplitude ξ¯, respectively, and noting that c and σ¯ are constants, yields(9)ifω-∧+cω-∧=σ¯ξ¯∧,which we can express in relation to ω-∧ as follows:(10)ω-∧=σ¯ξ¯∧if+c.Inserting this into our expression of the power spectral density in equation
(8) results in the power
spectral density of ω- as follows:(11)Sf=σ¯2ξ¯∧2f2+c2.Sinceσ¯2 and c2 are constants, the spectral density dependence on the
frequency f is solely determined by the 1/f2-term and the noise contribution ξ¯. If the noise is homogeneous in the frequency, i.e., it is
white noise, then we obtain a power law Sf∼f−2, i.e., Sf decays with an exponent of −2, consistent with our observations.To reproduce the observed spectrum on short time scales, we consider the
addition of negatively correlated Brownian fluctuationsηt so that the full dynamics reads(12)ω-t=−c∫0tω-t′dt′+σ∫0tξ-t′dt′+ηt.Indeed, the experimental spectra can be reproduced by this synthetic process;
see Figure6. We observe that instead
of a continued decaying power law, both the synthetic process but also the
measurements (in this case from Lisbon) do change qualitatively at the
one-second time scale with a positive slope of the spectrum.In this section, we have demonstrated new characteristics of the power-grid
frequency dynamics on a time scale below one second. Our analysis relied on
the availability of high-resolution frequency readings. In the next section,
we discuss how such recordings were obtained and validated.


### 2.2. Obtaining High-Quality Recordings from the Power Grid

To explain and investigate the origin of our frequency measurements, we begin by
defining the power-grid frequency. We model the power-grid voltage curveut at a measurement location in the grid by(13)ut=Asinθt,wherein θt is the instantaneous phase at time t and A the amplitude. The instantaneous frequency of this signal is
then defined as follows:(14)finstt=12πdθtdt.Withfinstt we have the definition of the power-grid frequency.
Unfortunately, though, finstt is not directly observable, so we need a measurement process
to determine some repeated approximations of finstt from ut.In the following, we will lay out this process in Subsections2.2.1–2.2.3. Making use of already published developments
and we focus on accuracy in II B 4. We extend the works by addressing the impact
of signal disturbances on Sf in II B 5.


#### 2.2.1. Data Acquisition

Our measurement device, the electrical data recorder (EDR), was introduced
some years ago and is being developed continuously; see also [18–21]. In the following, we provide a general overview of
the data acquisition and frequency calculation. We then specify how to
derive frequency data with high temporal resolution.The EDR captures the voltageut at a given connection point with a sampling frequency of
25 kHz. Using the time information provided by a GPS satellite
receiver, captured data are assigned a standardized time and date
information (UTC timestamp). Via the time stamps, the acquisition and the
subsequent data processing inherit the precision of GPS as a local standard
of the second. This allows to directly overlay and compare data from
different locations without the need for manual alignment in postprocessing.
Based on the recording of the ut time series, we extract features, such as the signal
frequency. We visualize the data processing chain in Figure 7.Figure 7
Conceptual illustration of the data acquisition and processing chain.
Data are acquired, processed, and stored on the electrical data
recorder device. Both raw samples and processed data are available
for in-depth analysis on a PC. The EDR also supports remote data
storage, live monitoring, and current measurements, which are not
depicted here. See [18, 45, 46] for details.


#### 2.2.2. Calculation of the Frequency

Signal frequency estimation in the EDR begins by passing the digitized
waveformut through a linear phase low-pass finite impulse response
(FIR) filter. For 25 kHz sampling rate and 50 Hz nominal
frequency, the filter has a length of L=200 and a 3 dB cut-off frequency close to 50 Hz.
Subsequently, the EDR determines periodwise frequency readings fpk by measuring the time between zero-crossings as follows:(15)fpk=1tzck−tzck−1.The series of zero-crossing timestzck are designated by instances where θ=0,2π,…, and k is the index of a specific zero-crossing. Looking at
zero-crossings makes the sampling times tzc dependent on the instantaneous phase of ut. For any practical signal, the elements of tzck are not equidistantly spaced on the time axis.Since data with nonuniform sampling is subsequently harder to analyze and
impossible to compare directly between locations and with other data
sources, the acquired information needs to be mapped to regular time
intervals. We call this desired regular output “frequency
reports”fn, where n points to a specific element of the series. Each report
carries a time stamp trep that is derived from a GPS-synchronized clock. The time
between trepn and trepn+1 is called “report interval.” In the
following, we describe how we calculate fn from fpk.


#### 2.2.3. Frequency Estimation-Increasing Time Resolution

As shown in the previous section, frequency estimates are especially
interesting at time scales around and below one second, as they reveal new
characteristics. Obtaining these estimates requires special attention with
regard to temporal alignment.For report intervals of one second and larger, it seems plausible to directly
aggregatefp and assign the results to tr-ignoring the slight misalignment due to the position of
the zero-crossings. We call this approach the full-period method (fp).
Figure 8 illustrates the aggregation
window used by this approach in red and the report interval in blue. The
frequency estimation ffp for the aggregation interval n is found by averaging over the K periods of the aggregation interval.(16)ffpn=1Kn∑k=1Knfpk,n.Figure 8
Zero-crossing algorithm (ZC algorithm). Illustration of the
zero-crossing positions relative to a 100 ms time
interval.Due to the time misalignment, this kind of estimation causes an error that
depends on both the phaseθt at the interval borders and the derivative of the
instantaneous frequency finst within the interval. The error thereby scales linearly
with the latter. In the worst case, a full period in the average does not
actually belong to the report interval n. For one-second records, this equals 2% of the
involved data and is usually tolerated. For 100 ms reports this
figure increases to 20% and causes intolerable deviations during
transient frequency events.To address this challenge, we introduced the partial period method (pp) in
[21] that performs a local
resampling of the data to the regular report intervals. Contributions of the
border periods (i.e., periods 1 and 6 in Figure 8) are partially included in the average according to
their temporal overlap with the report interval marked in blue as follows:(17)fppn=1Kn−2+w1+wK∑k=1Knwkfpk,n,where K is the number of periods that are in touch with the
aggregation interval n and wk are weights for the individual periods.The first weight is given by(18)w1=tzc1−trepn−1·fp1,and the last weight is given by(19)wK=trepn−tzcK−1·fpK,while all periods in between are equally weighted with wk=2…K−1=1.While the pp-method eliminates the majority of shift-related estimation
errors, subtle deviations remain. The main reason is that the approach
assumes thatfinst is constant during a period. This condition is hardly
satisfied for any real signal. Nevertheless, the estimation error in dynamic
conditions is drastically improved compared to the fp-method [21]. The partial period method,
therefore, practically solves the temporal alignment problem.


#### 2.2.4. Accuracy of Frequency Measurements

Since the individual period, measurements are only as good as the underlying
time reference, we use the GPS receivers’ pulse per second (PPS)
signal as a local frequency standard. The error in frequencyΔf relates to the uncertainty in the time Δt via(20)Δf=−T=1fnom−2Δt.The edge-to-edge error of our receiver (Garmin GPS-18LVC) and PPS acquisition
is less than 50 ns, which translates to a worst-case frequency error
of 125µHz in the 50 Hz case.
Experimentally, we further verify that the frequency estimation error in
steady-state conditions over a 10-minute interval is only
26.7 µHz at 50 Hz within a
95.4% confidence interval [45]. In dynamic conditions, however, the estimation errors are
dominated by the principles of the zero-crossing method. A detailed analysis
of the estimation errors in the time domain can be found in [21]. For the frequency ramps defined in
the standard for Phasor Measurement Units (IEEE 60255-118-1) we found the
error limit of our method to be ±0.6mHz.


#### 2.2.5. Influence of Signal Disturbances

As shown in II A 3. and Figure6, the
observed PSDs can be reproduced using a stochastic process composed of an
Ornstein-Uhlenbeck process (including white noise) and added correlated
fluctuations, where the latter lead to the increasing spectrum on a short
time scale. We thus investigate how our processing chain deals with
disturbances in the signal. By doing so, we also ensure that no artificial
signal is introduced into the system.From reviewing the signal chain, disturbances superimposed onut are the most likely source of the observed fluctuations,
since adding a random signal eut to ut will impact how well the period boundaries can be
determined. The relation between such voltage noise eu and noise in the period duration (time) measurement et is given by the gradient of ut at the zero-crossing as follows:(21)etk=euddtut=tZCk=2πAf0euk.Voltage noise, therefore, translates linearly to noise in period duration
(phase noise), scaled by power-grid frequency and the signal amplitude.
However, ifeu is not random, the influence might vanish. This is
especially the case for harmonics, where eu is periodic in multiples of f.To gain an impression of the influence of noise, we conduct the following
computer experiment. We generate a sinusoidal signal, superimpose it with
white Gaussian noise and process it with the estimation algorithm to obtain
a series of frequency readingsfn. We then determine the Welch PSD estimates Sf of this series as carried out for the measured frequency
time series. Results are shown in Figure 9 from which we make three major observations:(1)
The curves for 30 dB and 60 dB signal-to-noise
ratio (SNR) are similar in shape but vertically shifted by three
orders of magnitude. The shift corresponds exactly to the
30 dB difference in noise excitation levels. Therefore,
the noise levels inf and ut can be directly related.(2)
Sf rises approximately proportional to quadratic
with t−2 which is expected from frequency demodulator
theory [47]. Since the ZC
algorithm averages K period measurements, we would assume that the
total noise amplitude of f scales with 1/K. However, due to the increase in noise
density, the total noise amplitude in f increases linearly with K. Changing the resolution from 1 s to
100 ms translates to a tenfold increase in noise
amplitude (in the 60dB case, and σ=0.48mHzvs.4.9mHz in the 30dB case). Hence, the data with higher temporal
resolution contains more noise in absolute terms.(3)
For smaller values off we observe a slight flattening of S indicating a possible convergence to a
constant value. The shallow part indicates some feed-through of
white phase noise, which is even more pronounced for the
full-period method (shown in Supplementary Material). We find
the source of this effect in the violation of the previously
introduced assumption. The instantaneous frequency is not
constant during each period. Compared to the noise contributions
from the high-frequency end of the spectrum, the feed-through
has a minuscule effect on the results.Figure 9
Power spectral density off for the partial period zero-crossing method. Input
is a sine wave with the additional white noise of different levels.
PSD estimates saturate for a long period duration and rise nearly
quadratic (solid black line) towards shorter periods. For
60 dB SNR the dashed line indicates the maximum level of Sf at
0.01 mHz2·Hz−1.


#### 2.2.6. Quantification of Disturbance Influence

So far, we have characterized the effects of noise on frequency estimation.
We still need to quantify what noise contribution is to be expected in an
actual measurement environment so that we can establish a boundary of trust
for our results. For this, we need to (a) analyze the distortions in the
real-world setting and (b) characterize the susceptibility of the
algorithm.To approach (a), we investigate the properties ofut. In Figure 10, we
plot an excerpt of the waveform and the spectrum of the measured voltage.
The majority of the disturbance is in the nonsinusoidal shape of the
waveform and, therefore, in the harmonic content.Figure 10
Captured voltage waveform in Stockholm. Top: example of the captured
voltage waveformut. Data: Stockholm 9th of May 2019,
20:35. Bottom: power spectrum of one hour of data (Welch estimate,
resolution 1 Hz, Hann window). Signal-to-noise ratio (SNR)
including distortion: 29.6 dB. This is worst case SNR
observed for EDR measurements.From the spectrum ofut we can determine the power contained in the disturbances
to obtain an SNR for the measurements (details in Supplementary Material).
Throughout the campaigns, we find a minimum SNR of ≈30dB. However, if we selectively suppress the contribution of
the first ten harmonics, the minimum SNR rises to ≈60dB.Additionally to the above, short-term transients caused by local switch
operations will occasionally occur inut. From the time-domain analysis of the algorithm in [21], we know that such transients will
cause large (usually 0.1 Hz and above) single-point frequency
deviations. Since a pulse in the Dirac sense translates to a constant in the
frequency domain, those transients might cause slightly elevated Sf levels. However, most of these single-point deviations are
removed from the data via an outlier filter prior to our investigation.For (b), we already know from the processing chain, that components above
50 Hz are rejected by the initial low-pass filtering. It is further
known that static waveform deformations have no impact on the estimation
[21]. We are therefore convinced
that omitting the major harmonic contributions from the SNR estimation is
justified (further discussion in Supplementary Material).Consequently, the expected contribution of the ZC algorithm to the observed
PSD estimates follows the 60 dB curve in Figure9. This means that in all cases except Oldenburg the
algorithm’s contribution explains less than one-tenth of the observed Sf. In the case of Oldenburg, the factor is ≈1/4, while in most cases, it is in the range of ≈1/100. Hence, we are confident that the observed increases of Sf are not a measurement artifact but a signal property.


## 2.1. Emerging Phenomena on Short-Time Scales

The following observations are based on newly available 100 ms frequency
estimates and particularly focus on phenomena emerging at the newly-accessible
time scale below one second. The acquisition of the high-resolution frequency
data is described later in Subsection2.2.As one key quantity, we compute the power spectral density (PSD)Sf, which allows us to analyze the power distribution of
different oscillations in the components of the frequency recorded from the
system in a stochastic process [29]. As
with many processes observed in nature, we also expect the power spectrum of
self-affine signals to decay following a power law Sf∼f−β [30, 31]. PSD is a common tool that allows us to
perform further analysis on data originating from stochastic processes, such as
data obtained from electroencephalogram (EEG) [32], as well as for power-grid frequency signal analysis, which we
are currently considering here [29].
Furthermore, we may also estimate other properties such as the fractality, i.e.,
self-similarity on various scales, or the correlation of the time series, which
are closely connected to the spectrum [30, 33, 34].


### 2.1.1. Power Law Analysis

Computing the PSD for the power-grid frequency recordings reveals new
qualitative dynamics on the short time scale at essentially all measurement
locations, Figure2. We utilize
several independent recordings (see Figure 2(a) for a map) from the Continental European grid (Figures
2(c) and 2(d)) and from three different synchronous areas
Figure 2(b). Let us briefly remark
on two observations we do not discuss in detail.Figure 2
Power-grid dynamics change qualitatively at the 1-second time scale.
We utilize recordings from various locations in Europe (a) to then
plot the spectral densitySf of the power-grid frequency (b)–(d). This
analysis includes independent synchronous areas (b) as well as
Continental Europe data obtained from grid radar (c) and from the
EDR (d). Note that we plot Sf and also report the corresponding time scales t=1/f. While all recordings initially display an
approximately decaying power law β, several recordings display a rising spectrum with
slope γ for time scales below t∼1s.
(a)(b)(c)(d)First, on the longer time scale of∼3…10 seconds, we observe pronounced interarea oscillations,
where different regions within one synchronous area oscillate against each
other. We notice the spatial dependency in the spectrum as different
locations show different magnitudes and positions (in terms of time) of
these peaks. This is explained by the fact that, e.g., north-south and
east-west oscillations have different frequency; see also [15, 27, 35] for a more
detailed discussion of this phenomenon.Second, there are various narrow and distinct peaks on the time scale below
one second in the spectrum; see, e.g., Figure2(d) Karlsruhe and Oldenburg. Curiously, not all
recordings from the Continental European area display such peaks, e.g.,
neither Lisbon nor Istanbul. Other synchronous areas again show mixed
results: recordings from Sweden only show less pronounced peaks, while the
peaks are again present in the Faroe Island recordings. These peaks could
arise from characteristic eigenfrequencies of the power grid when subject to
a disturbance; see Figure 3. While
the short circuit with an active power step of approximately 440MW only leads to a frequency change of the order 10−6 (0.15 mHz) in a realistic grid simulation (see
Figure 3(a)), the dynamics take place
on the second to subsecond scale with eigenfrequencies of about
1.5 Hz. Note the 3 full oscillations between 2 and 4 seconds
in Figures 3(a) and 3(c). Further investigations will be
necessary to fully explore the origin and relevance of these peaks. In the
current study, we instead focus on another key observation: The qualitative
change in the spectrum around 1 second and the emergence of complex
phenomena at time scales <1 second (greyed out).Figure 3
Frequency and power dynamics in a transmission grid display
characteristic subsecond dynamics. We display the local frequencies
at 380/220/110-kV station DAXLA with 3 powerstation-blocks (a), the
real power of one block (b), as well as the absolute voltage
amplitude (c), and the absolute current on the line that displayed a
short circuit to trigger the event (d). The transient dynamics is
triggered via a three-phase short circuit between 0.5 and
0.65 seconds on a 110 kV line between Daxlanden and
Leopoldshafen (connection point of the KIT campus North). Line and
node parameters use realistic values for the South-West German power
grid. Frequency and voltage are given in the per-unit (p.u.) system,
i.e., at 1 p.u., they are at the reference value.On the long time scale oft>1 second, all measurements display a power law-like decay of
the spectral density Sf∼f−β∼tβ, with f the frequency and t=1/f the period. For ease of interpretation, we plot Sf as a function of f but also provide the corresponding 1/t values; see Figure 2. This decay follows approximately f−β, where we determine the exponent −β as the slope of the dashed line in the double-logarithmic
plot for each measurement site separately. As expected, the decay constant β is positive, indicating a decrease in the spectral density
for shorter periods (large frequencies) of the recordings. This decreasing
spectral density, i.e., the power law, has been observed in numerous
systems, ranging from many man-made systems, such as the World Wide Web,
over transport networks to biological systems [5, 6]. Power laws
remain an active field of research, emerging both on microscopic and
macroscopic scales. It often remains open to which extent the individual
units or their collective dynamics contribute to the observed power laws
[36].Meanwhile, the spectral density starts to increase on short time scales oft<1 second. We again characterize the slope of the spectrum
via a linear fit in the double-logarithmic plot and thereby determine the
slope γ, which in this case is smaller than zero for all
recordings from the Continental European power grid (Figures 2(c) and 2(d)). Interestingly, this negative γ is not exclusive to the Continental European region; see
Figure 2(b), where we show the power
spectral density of five other power-grid recordings, namely from the power
grids on the Faroe Islands, Russia (St. Petersburg), and Nordic (Sweden).
The measurements taken on the Faroe Islands and in St. Petersburg used again
the EDR, which was also used to make the recordings in Figure 2(a). In addition, we also include
three recordings provided by the Swedish TSO in three undisclosed distinct
locations in Sweden for the month of January 2020. These locations are
reported simply as north Sweden, middle
Sweden, and south Sweden, as provided by the
Swedish TSO.Let us systematically investigate the different slopesβ and γ of the spectral decays observed in Figure 2 as a scatter plot in Figure 4, see also Supplementary Material for
numerical details. First, we note that all measurements consistently report
a decaying power law with exponent β∼2.5. However, the dynamics at short time scales, given by γ is more varied γ∼−2…2. We identify three regimes: the Swedish recordings display
a continuous decay of the spectrum. Next, we observe several recordings with
small absolute values of γ, implying the spectrum becomes almost flat. Finally, many
recordings, mostly from the Continental European grid, display negative γ values and have an increasing spectrum, pointing towards a
complex process emerging at these time scales.Figure 4
Spectral measurements on short time scales are either positively,
negatively or not correlated. We plot estimates of the power law
exponentsβ (long time scale) and γ (short time scale).Given this surprising increase of the spectrum on short time scales, let us
elaborate on why this is almost certainly not a problem of the measurement
system but a genuine effect of the underlying power system. Subsection2.2. deals with the EDR data
acquisition and analyzes it towards reliability and trust. Additionally, we
make the following observations: First, the EDR measurements within one grid
are consistent: We clearly observe negative γ values in Continental Europe at several independent
measurement points, namely in Istanbul, Oldenburg, Karlsruhe, and Lisbon.
Second, the negative γ values are further confirmed for Continental Europe by
measurements from grid radar (Figure 2(c)), relying on a completely independent hardware and software
solution. Third, the EDR measurements show γ≈0 in some other synchronous areas, such as in Russia (St.
Petersburg) and on the Faroe Islands.Interestingly, the time series provided by the Swedish TSO has a positiveγ, while the EDR-recorded time series in Stockholm (recorded
during a different time period) yields a low negative value of γ. As we do not know the reason for this discrepancy, we
list some potential causes: the measurements by the EDR and the Swedish TSO
were taken at different times and could hence measure different states of
the power system. Furthermore, the EDR was connected to a power plug in a
hotel, i.e., to the low-voltage distribution grid used for residential
areas. Meanwhile, the Swedish TSO has access to measurements at power plants
or directly at the high-voltage transmission grid. While the frequency is
typically a global property, the negative γ values could possibly arise from the complex dynamics of
components at the local distribution layer, such as inverters and nonlinear
loads. Finally, we do not have access to the raw data recorded by the
Swedish TSO and hence cannot infer anything about their measurement setup,
or filtering while we discuss the limits of our measurements in detail.


### 2.1.2. Multifractality and Correlations

In our sample trajectories in Figure1, we already noted a certain multifractal structure to the
power-grid frequency recordings. One natural question arising is: what is
the exact nature of the fractal or multifractal structure that underlies
power-grid frequency?Power-grid systems are complex dynamical systems affected and impacted by
various endogenous and exogenous pressures from both within (control
mechanics, thermal effects, and oscillations) and from outside (changes in
generation and consumption, weather, and failures). Every impactful element
acts at various different time scales, and these time scales range from
years (ever-growing generation), to months (weather cycles), to weeks
(workweek-weekend cycles), to hours (market and dispatch changes), to
minutes (generation adjustments and automated control), to seconds (dynamic
interference, travelling waves, and occasional failures), to milliseconds
(electrodynamic effects) [10, 37].To best describe the power-grid frequency on very short time scales, we model
the time series as a combination of stochastic processes. To that end, we
focus on the incremental structure of our recordingsft, given by(1)Δfr=ft+r−ft,where r is our incremental lag. From a stochastic process
point-of-view it is usually denoted “increments.” These
increments carry information relating directly to the process’
fastest time scales, i.e., a process changes on a time scale of order r, as well as information relating to the existence of
memory and fractality in the recordings.One of the most common characteristics of the incremental time seriesΔfr is the “scaling” of the increments in
relation to the incremental lag r. This is commonly encapsulated in the structure function SqΔfr,r(2)SqΔfr,r=Δfrq∼rqhq,where hq is the generalized Hurst exponent, also known as the
Hölder exponent, which captures the change in the (noncentered and
nonnormalized) statistical moments, where q is the power of the statistical moment. hq is related to the scaling exponent τq, which is often used in turbulence analysis, as follows: τq=qhq−1 [38].We utilize the generalized Hurst exponent for two purposes: to investigate
correlations and to explore multifractality. For an uncorrelated Brownian
motion, we haveh2=0.5, while a correlated process has h2>0.5 and an anticorrelated process h2<0.5. Furthermore, a monofractal process is one where the
function hq is essentially constant, i.e., independent of q, while a multifractal process is characterized by a q-dependent generalized Hurst exponent hq.If the original time series yields aq-dependent hq, i.e., a multifractal process, this could either be due to
different long-range correlations or due to a broad probability
distribution. To ascertain how this multifractality emerges, we shuffle the
power-grid frequency recordings, i.e., randomize the position of each time
series entry. We can then calculate the generalized Hurst exponents hshuffledq for the shuffled data. Since shuffling the data removes
correlations, time series that are not distributional
multifractal will have hq=h2=1/2, i.e., they will be uncorrelated at all scales. If, on the
contrary, we find that hq is nonconstant, this is an indication of distributional
multifractality.Computing the generalized Hurst exponent reveals again a complex structure of
the time series on the short time scale; see Figure5. On a time scale of 0.2 to 1.0 seconds, the
generalized Hurst exponent is clearly dependent on q, indicating a multifractal process coherent with earlier
observations (Figure 1). This
multifractality persists even when shuffling, indicating a distributional
multifractality. For longer time scales, particularly above
4 seconds, the time series becomes monofractal (constant hq). Interestingly, the correlation changes over the time
scales: On short time scales of 0.4 second to 1 second, the
power-grid frequency is anticorrelated as indicated by a small value of h2<0.5. For these time scales, essentially, any disturbance on
the grid away from the reference frequency of 50Hz or 60Hz is followed by an opposite effect. While we do not expect
control actions to play a major role [10] on this time scale, the effect is similar, any disturbance
is counteracted, leading to a frequency time series close to the reference
value. A possible explanation for the emerging anticorrelation could be the
self-regulating effect of frequency-dependent loads [40], i.e., a frequency drop would be met by a reduced
load and thereby counteract the frequency change. Meanwhile, the time series
becomes positively correlated on longer time scales of several seconds and
more, as indicated by a large value of h2>0.5. This persistent behavior aligns with our expectations of
a power system subject to load and generation ramps on time scales of
seconds to minutes and an overall inertial system.Figure 5
Frequency recordings are multifractal and anticorrelated on short and
correlated on longer time scales. We plot the generalized Hurst
exponenthq of power-grid frequency recordings from the Lisbon
recording in Continental Europe at different time scales, from top
to bottom: 0.4∼1 s, 1∼4 s, and 4∼10 s. The smallest time scales show a
nonconstant exponent hq, which points to the existence of a large number
of degrees of freedom in the system. As the time scales become
larger, the system gets progressively more ordered, as fluctuations
are quenched by the strong synchrony in the system. At a time scale
of 4∼10 s the shuffled grid frequency recording
becomes indistinguishable from noise hq=0.5 [39].Let us briefly discuss the role of distributional multifractality before then
utilizing the insights about correlations in the time series to reproduce
the observed spectrum. We do not know for certain why and how the
multifractality arises in the frequency recordings discussed. However, we do
know that one mechanism giving rise to a multifractal time series is the
superposition of several processes at a shorter time scale, formalized under
the concept of superstatistics by Beck and Cohen [41]. Assuming the observed frequency time series arises
from such superstatistics, it would offer further directions of research. In
superstatistics, a system changes from one state to another while we treat
all observations as if they originated from the same state. In terms of the
energy system, these different states of the system could be related to
different generation mixes and demand patterns or due to drivers such as
superstatistical behavior in the wind [42]. Note that superstatistics has been discussed on power-grid
frequency recordings before [11], and
there are further hints that the dynamics of power-grid recordings get more
complex on shorter time scales. In fact, if we consider the temporal scales
wherein we observe the strongest evidence of multifractality, i.e., for t<4 s, this agrees well with previous results [27, 28], wherein power-grid frequency increments are shown to be
well described by a convolution of Gaussian statistics with a parametric
function of varying variance. All of this points to even more research is
necessary to understand the complexity of power-grid frequency time series
on this short time scale.


### 2.1.3. Stochastic Process

Let us return to the observation of an increasing spectrum and try to explain
this via physical processes on a short time scale. We have established from
our generalized Hurst exponent analysis that the frequency is positively
correlated on a longer time scale of several seconds and negatively
correlated on shorter time scales. For now, we consider simple mono-fractal
processes and explore how these can reproduce the observed spectrum. Indeed,
a very simple stochastic process can already reproduce such a spectrum,
namely combining a fractional, positively correlated
Ornstein–Uhlenbeck process with added, negatively correlated
fluctuations; see Figure6. Before
discussing Figure 6 in detail, let us
motivate the usage of an Ornstein–Uhlenbeck process:Figure 6
Reproducing the spectrum using an Ornstein-Uhlenbeck process with
added correlated short-term fluctuations.Following earlier work [11, 43], we may model the frequency
dynamics of the lossless high-voltage grid on the short time scale as follows:(3)ddtθi=ωi,Middtωi=Pi+σiξi−giωi+∑j=1NKijsinθj−θi,where we assume that each node i in the grid behaves approximately like a synchronous
machine whose state is characterized by its voltage phase angle θi and angular velocity ωi=2πfi−fR, which is proportional to the deviation of the frequency
from the reference frequency of fR=50Hz or fR=60Hz. The dynamics of the network are then determined by the
following parameters: inertia Mi, mechanical power Pi (positive for effective generators and negative for
effective consumers), fluctuations ξi with amplitude σi, damping gi, and the coupling matrix Kij.As we typically do not have precise information about the full network
topology, all nodal parameters and states, we want to simplify this to an
equation for the bulk angular velocityω-≔∑i=1NMiωi/∑i=1NMi by assuming a homogeneous damping to inertia ratio, c=gi/Mi [44] and balanced
average power ∑i=1NPi=0. We obtain the aggregated swing equation as follows:(4)ddtω-=−cω-+σ¯ξ¯t,with aggregated fluctuation amplitude σ¯ and aggregated fluctuations ξ¯, which depend on the precise nature of the noise [11]. The bulk frequency dynamics given
in equation (4) describes
the average behavior in one synchronous area, such as the Continental
European grid, neglecting specific local properties. Crucially, we have
simplified the nonlinear network dynamics to a relatively simple
Ornstein–Uhlenbeck process. One main difference to a standard
Ornstein–Uhlenbeck process is that the aggregated fluctuations ξ¯ are typically not well described by ordinary Brownian
motion, but instead, the effective fluctuations acting on the system are
either heavy-tailed itself or arise from a superposition of several
stochastic processes, so-called superstatistics [11, 41],
mentioned earlier.Regardless of the precise nature of the fluctuationsξ¯, for an Ornstein–Uhlenbeck process, we would expect
a spectrum following a simple power law (dotted line in Figure 6). In order to obtain an explicit
expression for the power-spectral density, recall that the energy of a time
series is given by(5)E≔∫−∞∞xt2dt.Parseval’s theorem dictates that(6)∫−∞∞xt2dt=∫−∞∞x∧f2df,with x∧f the Fourier transform of xt, i.e.,(7)x∧f=∫−∞∞xte−i2πftdt,with i the imaginary unit. Using Parseval’s theorem, we
see that the integrand of the energy in equation (5) is simply the spectral
density Sf given as follows:(8)Sf≔x∧f2.In the following, we are interested in a solution to the aggregated swing
equation (4) in Fourier
space. Defining ω-∧ and ξ¯∧ as the Fourier transforms of the bulk angular velocity ω- and the aggregated noise amplitude ξ¯, respectively, and noting that c and σ¯ are constants, yields(9)ifω-∧+cω-∧=σ¯ξ¯∧,which we can express in relation to ω-∧ as follows:(10)ω-∧=σ¯ξ¯∧if+c.Inserting this into our expression of the power spectral density in equation
(8) results in the power
spectral density of ω- as follows:(11)Sf=σ¯2ξ¯∧2f2+c2.Sinceσ¯2 and c2 are constants, the spectral density dependence on the
frequency f is solely determined by the 1/f2-term and the noise contribution ξ¯. If the noise is homogeneous in the frequency, i.e., it is
white noise, then we obtain a power law Sf∼f−2, i.e., Sf decays with an exponent of −2, consistent with our observations.To reproduce the observed spectrum on short time scales, we consider the
addition of negatively correlated Brownian fluctuationsηt so that the full dynamics reads(12)ω-t=−c∫0tω-t′dt′+σ∫0tξ-t′dt′+ηt.Indeed, the experimental spectra can be reproduced by this synthetic process;
see Figure6. We observe that instead
of a continued decaying power law, both the synthetic process but also the
measurements (in this case from Lisbon) do change qualitatively at the
one-second time scale with a positive slope of the spectrum.In this section, we have demonstrated new characteristics of the power-grid
frequency dynamics on a time scale below one second. Our analysis relied on
the availability of high-resolution frequency readings. In the next section,
we discuss how such recordings were obtained and validated.


## 2.1.1. Power Law Analysis

Computing the PSD for the power-grid frequency recordings reveals new
qualitative dynamics on the short time scale at essentially all measurement
locations, Figure2. We utilize
several independent recordings (see Figure 2(a) for a map) from the Continental European grid (Figures
2(c) and 2(d)) and from three different synchronous areas
Figure 2(b). Let us briefly remark
on two observations we do not discuss in detail.Figure 2
Power-grid dynamics change qualitatively at the 1-second time scale.
We utilize recordings from various locations in Europe (a) to then
plot the spectral densitySf of the power-grid frequency (b)–(d). This
analysis includes independent synchronous areas (b) as well as
Continental Europe data obtained from grid radar (c) and from the
EDR (d). Note that we plot Sf and also report the corresponding time scales t=1/f. While all recordings initially display an
approximately decaying power law β, several recordings display a rising spectrum with
slope γ for time scales below t∼1s.
(a)(b)(c)(d)First, on the longer time scale of∼3…10 seconds, we observe pronounced interarea oscillations,
where different regions within one synchronous area oscillate against each
other. We notice the spatial dependency in the spectrum as different
locations show different magnitudes and positions (in terms of time) of
these peaks. This is explained by the fact that, e.g., north-south and
east-west oscillations have different frequency; see also [15, 27, 35] for a more
detailed discussion of this phenomenon.Second, there are various narrow and distinct peaks on the time scale below
one second in the spectrum; see, e.g., Figure2(d) Karlsruhe and Oldenburg. Curiously, not all
recordings from the Continental European area display such peaks, e.g.,
neither Lisbon nor Istanbul. Other synchronous areas again show mixed
results: recordings from Sweden only show less pronounced peaks, while the
peaks are again present in the Faroe Island recordings. These peaks could
arise from characteristic eigenfrequencies of the power grid when subject to
a disturbance; see Figure 3. While
the short circuit with an active power step of approximately 440MW only leads to a frequency change of the order 10−6 (0.15 mHz) in a realistic grid simulation (see
Figure 3(a)), the dynamics take place
on the second to subsecond scale with eigenfrequencies of about
1.5 Hz. Note the 3 full oscillations between 2 and 4 seconds
in Figures 3(a) and 3(c). Further investigations will be
necessary to fully explore the origin and relevance of these peaks. In the
current study, we instead focus on another key observation: The qualitative
change in the spectrum around 1 second and the emergence of complex
phenomena at time scales <1 second (greyed out).Figure 3
Frequency and power dynamics in a transmission grid display
characteristic subsecond dynamics. We display the local frequencies
at 380/220/110-kV station DAXLA with 3 powerstation-blocks (a), the
real power of one block (b), as well as the absolute voltage
amplitude (c), and the absolute current on the line that displayed a
short circuit to trigger the event (d). The transient dynamics is
triggered via a three-phase short circuit between 0.5 and
0.65 seconds on a 110 kV line between Daxlanden and
Leopoldshafen (connection point of the KIT campus North). Line and
node parameters use realistic values for the South-West German power
grid. Frequency and voltage are given in the per-unit (p.u.) system,
i.e., at 1 p.u., they are at the reference value.On the long time scale oft>1 second, all measurements display a power law-like decay of
the spectral density Sf∼f−β∼tβ, with f the frequency and t=1/f the period. For ease of interpretation, we plot Sf as a function of f but also provide the corresponding 1/t values; see Figure 2. This decay follows approximately f−β, where we determine the exponent −β as the slope of the dashed line in the double-logarithmic
plot for each measurement site separately. As expected, the decay constant β is positive, indicating a decrease in the spectral density
for shorter periods (large frequencies) of the recordings. This decreasing
spectral density, i.e., the power law, has been observed in numerous
systems, ranging from many man-made systems, such as the World Wide Web,
over transport networks to biological systems [5, 6]. Power laws
remain an active field of research, emerging both on microscopic and
macroscopic scales. It often remains open to which extent the individual
units or their collective dynamics contribute to the observed power laws
[36].Meanwhile, the spectral density starts to increase on short time scales oft<1 second. We again characterize the slope of the spectrum
via a linear fit in the double-logarithmic plot and thereby determine the
slope γ, which in this case is smaller than zero for all
recordings from the Continental European power grid (Figures 2(c) and 2(d)). Interestingly, this negative γ is not exclusive to the Continental European region; see
Figure 2(b), where we show the power
spectral density of five other power-grid recordings, namely from the power
grids on the Faroe Islands, Russia (St. Petersburg), and Nordic (Sweden).
The measurements taken on the Faroe Islands and in St. Petersburg used again
the EDR, which was also used to make the recordings in Figure 2(a). In addition, we also include
three recordings provided by the Swedish TSO in three undisclosed distinct
locations in Sweden for the month of January 2020. These locations are
reported simply as north Sweden, middle
Sweden, and south Sweden, as provided by the
Swedish TSO.Let us systematically investigate the different slopesβ and γ of the spectral decays observed in Figure 2 as a scatter plot in Figure 4, see also Supplementary Material for
numerical details. First, we note that all measurements consistently report
a decaying power law with exponent β∼2.5. However, the dynamics at short time scales, given by γ is more varied γ∼−2…2. We identify three regimes: the Swedish recordings display
a continuous decay of the spectrum. Next, we observe several recordings with
small absolute values of γ, implying the spectrum becomes almost flat. Finally, many
recordings, mostly from the Continental European grid, display negative γ values and have an increasing spectrum, pointing towards a
complex process emerging at these time scales.Figure 4
Spectral measurements on short time scales are either positively,
negatively or not correlated. We plot estimates of the power law
exponentsβ (long time scale) and γ (short time scale).Given this surprising increase of the spectrum on short time scales, let us
elaborate on why this is almost certainly not a problem of the measurement
system but a genuine effect of the underlying power system. Subsection2.2. deals with the EDR data
acquisition and analyzes it towards reliability and trust. Additionally, we
make the following observations: First, the EDR measurements within one grid
are consistent: We clearly observe negative γ values in Continental Europe at several independent
measurement points, namely in Istanbul, Oldenburg, Karlsruhe, and Lisbon.
Second, the negative γ values are further confirmed for Continental Europe by
measurements from grid radar (Figure 2(c)), relying on a completely independent hardware and software
solution. Third, the EDR measurements show γ≈0 in some other synchronous areas, such as in Russia (St.
Petersburg) and on the Faroe Islands.Interestingly, the time series provided by the Swedish TSO has a positiveγ, while the EDR-recorded time series in Stockholm (recorded
during a different time period) yields a low negative value of γ. As we do not know the reason for this discrepancy, we
list some potential causes: the measurements by the EDR and the Swedish TSO
were taken at different times and could hence measure different states of
the power system. Furthermore, the EDR was connected to a power plug in a
hotel, i.e., to the low-voltage distribution grid used for residential
areas. Meanwhile, the Swedish TSO has access to measurements at power plants
or directly at the high-voltage transmission grid. While the frequency is
typically a global property, the negative γ values could possibly arise from the complex dynamics of
components at the local distribution layer, such as inverters and nonlinear
loads. Finally, we do not have access to the raw data recorded by the
Swedish TSO and hence cannot infer anything about their measurement setup,
or filtering while we discuss the limits of our measurements in detail.


## 2.1.2. Multifractality and Correlations

In our sample trajectories in Figure1, we already noted a certain multifractal structure to the
power-grid frequency recordings. One natural question arising is: what is
the exact nature of the fractal or multifractal structure that underlies
power-grid frequency?Power-grid systems are complex dynamical systems affected and impacted by
various endogenous and exogenous pressures from both within (control
mechanics, thermal effects, and oscillations) and from outside (changes in
generation and consumption, weather, and failures). Every impactful element
acts at various different time scales, and these time scales range from
years (ever-growing generation), to months (weather cycles), to weeks
(workweek-weekend cycles), to hours (market and dispatch changes), to
minutes (generation adjustments and automated control), to seconds (dynamic
interference, travelling waves, and occasional failures), to milliseconds
(electrodynamic effects) [10, 37].To best describe the power-grid frequency on very short time scales, we model
the time series as a combination of stochastic processes. To that end, we
focus on the incremental structure of our recordingsft, given by(1)Δfr=ft+r−ft,where r is our incremental lag. From a stochastic process
point-of-view it is usually denoted “increments.” These
increments carry information relating directly to the process’
fastest time scales, i.e., a process changes on a time scale of order r, as well as information relating to the existence of
memory and fractality in the recordings.One of the most common characteristics of the incremental time seriesΔfr is the “scaling” of the increments in
relation to the incremental lag r. This is commonly encapsulated in the structure function SqΔfr,r(2)SqΔfr,r=Δfrq∼rqhq,where hq is the generalized Hurst exponent, also known as the
Hölder exponent, which captures the change in the (noncentered and
nonnormalized) statistical moments, where q is the power of the statistical moment. hq is related to the scaling exponent τq, which is often used in turbulence analysis, as follows: τq=qhq−1 [38].We utilize the generalized Hurst exponent for two purposes: to investigate
correlations and to explore multifractality. For an uncorrelated Brownian
motion, we haveh2=0.5, while a correlated process has h2>0.5 and an anticorrelated process h2<0.5. Furthermore, a monofractal process is one where the
function hq is essentially constant, i.e., independent of q, while a multifractal process is characterized by a q-dependent generalized Hurst exponent hq.If the original time series yields aq-dependent hq, i.e., a multifractal process, this could either be due to
different long-range correlations or due to a broad probability
distribution. To ascertain how this multifractality emerges, we shuffle the
power-grid frequency recordings, i.e., randomize the position of each time
series entry. We can then calculate the generalized Hurst exponents hshuffledq for the shuffled data. Since shuffling the data removes
correlations, time series that are not distributional
multifractal will have hq=h2=1/2, i.e., they will be uncorrelated at all scales. If, on the
contrary, we find that hq is nonconstant, this is an indication of distributional
multifractality.Computing the generalized Hurst exponent reveals again a complex structure of
the time series on the short time scale; see Figure5. On a time scale of 0.2 to 1.0 seconds, the
generalized Hurst exponent is clearly dependent on q, indicating a multifractal process coherent with earlier
observations (Figure 1). This
multifractality persists even when shuffling, indicating a distributional
multifractality. For longer time scales, particularly above
4 seconds, the time series becomes monofractal (constant hq). Interestingly, the correlation changes over the time
scales: On short time scales of 0.4 second to 1 second, the
power-grid frequency is anticorrelated as indicated by a small value of h2<0.5. For these time scales, essentially, any disturbance on
the grid away from the reference frequency of 50Hz or 60Hz is followed by an opposite effect. While we do not expect
control actions to play a major role [10] on this time scale, the effect is similar, any disturbance
is counteracted, leading to a frequency time series close to the reference
value. A possible explanation for the emerging anticorrelation could be the
self-regulating effect of frequency-dependent loads [40], i.e., a frequency drop would be met by a reduced
load and thereby counteract the frequency change. Meanwhile, the time series
becomes positively correlated on longer time scales of several seconds and
more, as indicated by a large value of h2>0.5. This persistent behavior aligns with our expectations of
a power system subject to load and generation ramps on time scales of
seconds to minutes and an overall inertial system.Figure 5
Frequency recordings are multifractal and anticorrelated on short and
correlated on longer time scales. We plot the generalized Hurst
exponenthq of power-grid frequency recordings from the Lisbon
recording in Continental Europe at different time scales, from top
to bottom: 0.4∼1 s, 1∼4 s, and 4∼10 s. The smallest time scales show a
nonconstant exponent hq, which points to the existence of a large number
of degrees of freedom in the system. As the time scales become
larger, the system gets progressively more ordered, as fluctuations
are quenched by the strong synchrony in the system. At a time scale
of 4∼10 s the shuffled grid frequency recording
becomes indistinguishable from noise hq=0.5 [39].Let us briefly discuss the role of distributional multifractality before then
utilizing the insights about correlations in the time series to reproduce
the observed spectrum. We do not know for certain why and how the
multifractality arises in the frequency recordings discussed. However, we do
know that one mechanism giving rise to a multifractal time series is the
superposition of several processes at a shorter time scale, formalized under
the concept of superstatistics by Beck and Cohen [41]. Assuming the observed frequency time series arises
from such superstatistics, it would offer further directions of research. In
superstatistics, a system changes from one state to another while we treat
all observations as if they originated from the same state. In terms of the
energy system, these different states of the system could be related to
different generation mixes and demand patterns or due to drivers such as
superstatistical behavior in the wind [42]. Note that superstatistics has been discussed on power-grid
frequency recordings before [11], and
there are further hints that the dynamics of power-grid recordings get more
complex on shorter time scales. In fact, if we consider the temporal scales
wherein we observe the strongest evidence of multifractality, i.e., for t<4 s, this agrees well with previous results [27, 28], wherein power-grid frequency increments are shown to be
well described by a convolution of Gaussian statistics with a parametric
function of varying variance. All of this points to even more research is
necessary to understand the complexity of power-grid frequency time series
on this short time scale.


## 2.1.3. Stochastic Process

Let us return to the observation of an increasing spectrum and try to explain
this via physical processes on a short time scale. We have established from
our generalized Hurst exponent analysis that the frequency is positively
correlated on a longer time scale of several seconds and negatively
correlated on shorter time scales. For now, we consider simple mono-fractal
processes and explore how these can reproduce the observed spectrum. Indeed,
a very simple stochastic process can already reproduce such a spectrum,
namely combining a fractional, positively correlated
Ornstein–Uhlenbeck process with added, negatively correlated
fluctuations; see Figure6. Before
discussing Figure 6 in detail, let us
motivate the usage of an Ornstein–Uhlenbeck process:Figure 6
Reproducing the spectrum using an Ornstein-Uhlenbeck process with
added correlated short-term fluctuations.Following earlier work [11, 43], we may model the frequency
dynamics of the lossless high-voltage grid on the short time scale as follows:(3)ddtθi=ωi,Middtωi=Pi+σiξi−giωi+∑j=1NKijsinθj−θi,where we assume that each node i in the grid behaves approximately like a synchronous
machine whose state is characterized by its voltage phase angle θi and angular velocity ωi=2πfi−fR, which is proportional to the deviation of the frequency
from the reference frequency of fR=50Hz or fR=60Hz. The dynamics of the network are then determined by the
following parameters: inertia Mi, mechanical power Pi (positive for effective generators and negative for
effective consumers), fluctuations ξi with amplitude σi, damping gi, and the coupling matrix Kij.As we typically do not have precise information about the full network
topology, all nodal parameters and states, we want to simplify this to an
equation for the bulk angular velocityω-≔∑i=1NMiωi/∑i=1NMi by assuming a homogeneous damping to inertia ratio, c=gi/Mi [44] and balanced
average power ∑i=1NPi=0. We obtain the aggregated swing equation as follows:(4)ddtω-=−cω-+σ¯ξ¯t,with aggregated fluctuation amplitude σ¯ and aggregated fluctuations ξ¯, which depend on the precise nature of the noise [11]. The bulk frequency dynamics given
in equation (4) describes
the average behavior in one synchronous area, such as the Continental
European grid, neglecting specific local properties. Crucially, we have
simplified the nonlinear network dynamics to a relatively simple
Ornstein–Uhlenbeck process. One main difference to a standard
Ornstein–Uhlenbeck process is that the aggregated fluctuations ξ¯ are typically not well described by ordinary Brownian
motion, but instead, the effective fluctuations acting on the system are
either heavy-tailed itself or arise from a superposition of several
stochastic processes, so-called superstatistics [11, 41],
mentioned earlier.Regardless of the precise nature of the fluctuationsξ¯, for an Ornstein–Uhlenbeck process, we would expect
a spectrum following a simple power law (dotted line in Figure 6). In order to obtain an explicit
expression for the power-spectral density, recall that the energy of a time
series is given by(5)E≔∫−∞∞xt2dt.Parseval’s theorem dictates that(6)∫−∞∞xt2dt=∫−∞∞x∧f2df,with x∧f the Fourier transform of xt, i.e.,(7)x∧f=∫−∞∞xte−i2πftdt,with i the imaginary unit. Using Parseval’s theorem, we
see that the integrand of the energy in equation (5) is simply the spectral
density Sf given as follows:(8)Sf≔x∧f2.In the following, we are interested in a solution to the aggregated swing
equation (4) in Fourier
space. Defining ω-∧ and ξ¯∧ as the Fourier transforms of the bulk angular velocity ω- and the aggregated noise amplitude ξ¯, respectively, and noting that c and σ¯ are constants, yields(9)ifω-∧+cω-∧=σ¯ξ¯∧,which we can express in relation to ω-∧ as follows:(10)ω-∧=σ¯ξ¯∧if+c.Inserting this into our expression of the power spectral density in equation
(8) results in the power
spectral density of ω- as follows:(11)Sf=σ¯2ξ¯∧2f2+c2.Sinceσ¯2 and c2 are constants, the spectral density dependence on the
frequency f is solely determined by the 1/f2-term and the noise contribution ξ¯. If the noise is homogeneous in the frequency, i.e., it is
white noise, then we obtain a power law Sf∼f−2, i.e., Sf decays with an exponent of −2, consistent with our observations.To reproduce the observed spectrum on short time scales, we consider the
addition of negatively correlated Brownian fluctuationsηt so that the full dynamics reads(12)ω-t=−c∫0tω-t′dt′+σ∫0tξ-t′dt′+ηt.Indeed, the experimental spectra can be reproduced by this synthetic process;
see Figure6. We observe that instead
of a continued decaying power law, both the synthetic process but also the
measurements (in this case from Lisbon) do change qualitatively at the
one-second time scale with a positive slope of the spectrum.In this section, we have demonstrated new characteristics of the power-grid
frequency dynamics on a time scale below one second. Our analysis relied on
the availability of high-resolution frequency readings. In the next section,
we discuss how such recordings were obtained and validated.


## 2.2. Obtaining High-Quality Recordings from the Power Grid

To explain and investigate the origin of our frequency measurements, we begin by
defining the power-grid frequency. We model the power-grid voltage curveut at a measurement location in the grid by(13)ut=Asinθt,wherein θt is the instantaneous phase at time t and A the amplitude. The instantaneous frequency of this signal is
then defined as follows:(14)finstt=12πdθtdt.Withfinstt we have the definition of the power-grid frequency.
Unfortunately, though, finstt is not directly observable, so we need a measurement process
to determine some repeated approximations of finstt from ut.In the following, we will lay out this process in Subsections2.2.1–2.2.3. Making use of already published developments
and we focus on accuracy in II B 4. We extend the works by addressing the impact
of signal disturbances on Sf in II B 5.


### 2.2.1. Data Acquisition

Our measurement device, the electrical data recorder (EDR), was introduced
some years ago and is being developed continuously; see also [18–21]. In the following, we provide a general overview of
the data acquisition and frequency calculation. We then specify how to
derive frequency data with high temporal resolution.The EDR captures the voltageut at a given connection point with a sampling frequency of
25 kHz. Using the time information provided by a GPS satellite
receiver, captured data are assigned a standardized time and date
information (UTC timestamp). Via the time stamps, the acquisition and the
subsequent data processing inherit the precision of GPS as a local standard
of the second. This allows to directly overlay and compare data from
different locations without the need for manual alignment in postprocessing.
Based on the recording of the ut time series, we extract features, such as the signal
frequency. We visualize the data processing chain in Figure 7.Figure 7
Conceptual illustration of the data acquisition and processing chain.
Data are acquired, processed, and stored on the electrical data
recorder device. Both raw samples and processed data are available
for in-depth analysis on a PC. The EDR also supports remote data
storage, live monitoring, and current measurements, which are not
depicted here. See [18, 45, 46] for details.


### 2.2.2. Calculation of the Frequency

Signal frequency estimation in the EDR begins by passing the digitized
waveformut through a linear phase low-pass finite impulse response
(FIR) filter. For 25 kHz sampling rate and 50 Hz nominal
frequency, the filter has a length of L=200 and a 3 dB cut-off frequency close to 50 Hz.
Subsequently, the EDR determines periodwise frequency readings fpk by measuring the time between zero-crossings as follows:(15)fpk=1tzck−tzck−1.The series of zero-crossing timestzck are designated by instances where θ=0,2π,…, and k is the index of a specific zero-crossing. Looking at
zero-crossings makes the sampling times tzc dependent on the instantaneous phase of ut. For any practical signal, the elements of tzck are not equidistantly spaced on the time axis.Since data with nonuniform sampling is subsequently harder to analyze and
impossible to compare directly between locations and with other data
sources, the acquired information needs to be mapped to regular time
intervals. We call this desired regular output “frequency
reports”fn, where n points to a specific element of the series. Each report
carries a time stamp trep that is derived from a GPS-synchronized clock. The time
between trepn and trepn+1 is called “report interval.” In the
following, we describe how we calculate fn from fpk.


### 2.2.3. Frequency Estimation-Increasing Time Resolution

As shown in the previous section, frequency estimates are especially
interesting at time scales around and below one second, as they reveal new
characteristics. Obtaining these estimates requires special attention with
regard to temporal alignment.For report intervals of one second and larger, it seems plausible to directly
aggregatefp and assign the results to tr-ignoring the slight misalignment due to the position of
the zero-crossings. We call this approach the full-period method (fp).
Figure 8 illustrates the aggregation
window used by this approach in red and the report interval in blue. The
frequency estimation ffp for the aggregation interval n is found by averaging over the K periods of the aggregation interval.(16)ffpn=1Kn∑k=1Knfpk,n.Figure 8
Zero-crossing algorithm (ZC algorithm). Illustration of the
zero-crossing positions relative to a 100 ms time
interval.Due to the time misalignment, this kind of estimation causes an error that
depends on both the phaseθt at the interval borders and the derivative of the
instantaneous frequency finst within the interval. The error thereby scales linearly
with the latter. In the worst case, a full period in the average does not
actually belong to the report interval n. For one-second records, this equals 2% of the
involved data and is usually tolerated. For 100 ms reports this
figure increases to 20% and causes intolerable deviations during
transient frequency events.To address this challenge, we introduced the partial period method (pp) in
[21] that performs a local
resampling of the data to the regular report intervals. Contributions of the
border periods (i.e., periods 1 and 6 in Figure 8) are partially included in the average according to
their temporal overlap with the report interval marked in blue as follows:(17)fppn=1Kn−2+w1+wK∑k=1Knwkfpk,n,where K is the number of periods that are in touch with the
aggregation interval n and wk are weights for the individual periods.The first weight is given by(18)w1=tzc1−trepn−1·fp1,and the last weight is given by(19)wK=trepn−tzcK−1·fpK,while all periods in between are equally weighted with wk=2…K−1=1.While the pp-method eliminates the majority of shift-related estimation
errors, subtle deviations remain. The main reason is that the approach
assumes thatfinst is constant during a period. This condition is hardly
satisfied for any real signal. Nevertheless, the estimation error in dynamic
conditions is drastically improved compared to the fp-method [21]. The partial period method,
therefore, practically solves the temporal alignment problem.


### 2.2.4. Accuracy of Frequency Measurements

Since the individual period, measurements are only as good as the underlying
time reference, we use the GPS receivers’ pulse per second (PPS)
signal as a local frequency standard. The error in frequencyΔf relates to the uncertainty in the time Δt via(20)Δf=−T=1fnom−2Δt.The edge-to-edge error of our receiver (Garmin GPS-18LVC) and PPS acquisition
is less than 50 ns, which translates to a worst-case frequency error
of 125µHz in the 50 Hz case.
Experimentally, we further verify that the frequency estimation error in
steady-state conditions over a 10-minute interval is only
26.7 µHz at 50 Hz within a
95.4% confidence interval [45]. In dynamic conditions, however, the estimation errors are
dominated by the principles of the zero-crossing method. A detailed analysis
of the estimation errors in the time domain can be found in [21]. For the frequency ramps defined in
the standard for Phasor Measurement Units (IEEE 60255-118-1) we found the
error limit of our method to be ±0.6mHz.


### 2.2.5. Influence of Signal Disturbances

As shown in II A 3. and Figure6, the
observed PSDs can be reproduced using a stochastic process composed of an
Ornstein-Uhlenbeck process (including white noise) and added correlated
fluctuations, where the latter lead to the increasing spectrum on a short
time scale. We thus investigate how our processing chain deals with
disturbances in the signal. By doing so, we also ensure that no artificial
signal is introduced into the system.From reviewing the signal chain, disturbances superimposed onut are the most likely source of the observed fluctuations,
since adding a random signal eut to ut will impact how well the period boundaries can be
determined. The relation between such voltage noise eu and noise in the period duration (time) measurement et is given by the gradient of ut at the zero-crossing as follows:(21)etk=euddtut=tZCk=2πAf0euk.Voltage noise, therefore, translates linearly to noise in period duration
(phase noise), scaled by power-grid frequency and the signal amplitude.
However, ifeu is not random, the influence might vanish. This is
especially the case for harmonics, where eu is periodic in multiples of f.To gain an impression of the influence of noise, we conduct the following
computer experiment. We generate a sinusoidal signal, superimpose it with
white Gaussian noise and process it with the estimation algorithm to obtain
a series of frequency readingsfn. We then determine the Welch PSD estimates Sf of this series as carried out for the measured frequency
time series. Results are shown in Figure 9 from which we make three major observations:(1)
The curves for 30 dB and 60 dB signal-to-noise
ratio (SNR) are similar in shape but vertically shifted by three
orders of magnitude. The shift corresponds exactly to the
30 dB difference in noise excitation levels. Therefore,
the noise levels inf and ut can be directly related.(2)
Sf rises approximately proportional to quadratic
with t−2 which is expected from frequency demodulator
theory [47]. Since the ZC
algorithm averages K period measurements, we would assume that the
total noise amplitude of f scales with 1/K. However, due to the increase in noise
density, the total noise amplitude in f increases linearly with K. Changing the resolution from 1 s to
100 ms translates to a tenfold increase in noise
amplitude (in the 60dB case, and σ=0.48mHzvs.4.9mHz in the 30dB case). Hence, the data with higher temporal
resolution contains more noise in absolute terms.(3)
For smaller values off we observe a slight flattening of S indicating a possible convergence to a
constant value. The shallow part indicates some feed-through of
white phase noise, which is even more pronounced for the
full-period method (shown in Supplementary Material). We find
the source of this effect in the violation of the previously
introduced assumption. The instantaneous frequency is not
constant during each period. Compared to the noise contributions
from the high-frequency end of the spectrum, the feed-through
has a minuscule effect on the results.Figure 9
Power spectral density off for the partial period zero-crossing method. Input
is a sine wave with the additional white noise of different levels.
PSD estimates saturate for a long period duration and rise nearly
quadratic (solid black line) towards shorter periods. For
60 dB SNR the dashed line indicates the maximum level of Sf at
0.01 mHz2·Hz−1.


### 2.2.6. Quantification of Disturbance Influence

So far, we have characterized the effects of noise on frequency estimation.
We still need to quantify what noise contribution is to be expected in an
actual measurement environment so that we can establish a boundary of trust
for our results. For this, we need to (a) analyze the distortions in the
real-world setting and (b) characterize the susceptibility of the
algorithm.To approach (a), we investigate the properties ofut. In Figure 10, we
plot an excerpt of the waveform and the spectrum of the measured voltage.
The majority of the disturbance is in the nonsinusoidal shape of the
waveform and, therefore, in the harmonic content.Figure 10
Captured voltage waveform in Stockholm. Top: example of the captured
voltage waveformut. Data: Stockholm 9th of May 2019,
20:35. Bottom: power spectrum of one hour of data (Welch estimate,
resolution 1 Hz, Hann window). Signal-to-noise ratio (SNR)
including distortion: 29.6 dB. This is worst case SNR
observed for EDR measurements.From the spectrum ofut we can determine the power contained in the disturbances
to obtain an SNR for the measurements (details in Supplementary Material).
Throughout the campaigns, we find a minimum SNR of ≈30dB. However, if we selectively suppress the contribution of
the first ten harmonics, the minimum SNR rises to ≈60dB.Additionally to the above, short-term transients caused by local switch
operations will occasionally occur inut. From the time-domain analysis of the algorithm in [21], we know that such transients will
cause large (usually 0.1 Hz and above) single-point frequency
deviations. Since a pulse in the Dirac sense translates to a constant in the
frequency domain, those transients might cause slightly elevated Sf levels. However, most of these single-point deviations are
removed from the data via an outlier filter prior to our investigation.For (b), we already know from the processing chain, that components above
50 Hz are rejected by the initial low-pass filtering. It is further
known that static waveform deformations have no impact on the estimation
[21]. We are therefore convinced
that omitting the major harmonic contributions from the SNR estimation is
justified (further discussion in Supplementary Material).Consequently, the expected contribution of the ZC algorithm to the observed
PSD estimates follows the 60 dB curve in Figure9. This means that in all cases except Oldenburg the
algorithm’s contribution explains less than one-tenth of the observed Sf. In the case of Oldenburg, the factor is ≈1/4, while in most cases, it is in the range of ≈1/100. Hence, we are confident that the observed increases of Sf are not a measurement artifact but a signal property.


## 2.2.1. Data Acquisition

Our measurement device, the electrical data recorder (EDR), was introduced
some years ago and is being developed continuously; see also [18–21]. In the following, we provide a general overview of
the data acquisition and frequency calculation. We then specify how to
derive frequency data with high temporal resolution.The EDR captures the voltageut at a given connection point with a sampling frequency of
25 kHz. Using the time information provided by a GPS satellite
receiver, captured data are assigned a standardized time and date
information (UTC timestamp). Via the time stamps, the acquisition and the
subsequent data processing inherit the precision of GPS as a local standard
of the second. This allows to directly overlay and compare data from
different locations without the need for manual alignment in postprocessing.
Based on the recording of the ut time series, we extract features, such as the signal
frequency. We visualize the data processing chain in Figure 7.Figure 7
Conceptual illustration of the data acquisition and processing chain.
Data are acquired, processed, and stored on the electrical data
recorder device. Both raw samples and processed data are available
for in-depth analysis on a PC. The EDR also supports remote data
storage, live monitoring, and current measurements, which are not
depicted here. See [18, 45, 46] for details.


## 2.2.2. Calculation of the Frequency

Signal frequency estimation in the EDR begins by passing the digitized
waveformut through a linear phase low-pass finite impulse response
(FIR) filter. For 25 kHz sampling rate and 50 Hz nominal
frequency, the filter has a length of L=200 and a 3 dB cut-off frequency close to 50 Hz.
Subsequently, the EDR determines periodwise frequency readings fpk by measuring the time between zero-crossings as follows:(15)fpk=1tzck−tzck−1.The series of zero-crossing timestzck are designated by instances where θ=0,2π,…, and k is the index of a specific zero-crossing. Looking at
zero-crossings makes the sampling times tzc dependent on the instantaneous phase of ut. For any practical signal, the elements of tzck are not equidistantly spaced on the time axis.Since data with nonuniform sampling is subsequently harder to analyze and
impossible to compare directly between locations and with other data
sources, the acquired information needs to be mapped to regular time
intervals. We call this desired regular output “frequency
reports”fn, where n points to a specific element of the series. Each report
carries a time stamp trep that is derived from a GPS-synchronized clock. The time
between trepn and trepn+1 is called “report interval.” In the
following, we describe how we calculate fn from fpk.


## 2.2.3. Frequency Estimation-Increasing Time Resolution

As shown in the previous section, frequency estimates are especially
interesting at time scales around and below one second, as they reveal new
characteristics. Obtaining these estimates requires special attention with
regard to temporal alignment.For report intervals of one second and larger, it seems plausible to directly
aggregatefp and assign the results to tr-ignoring the slight misalignment due to the position of
the zero-crossings. We call this approach the full-period method (fp).
Figure 8 illustrates the aggregation
window used by this approach in red and the report interval in blue. The
frequency estimation ffp for the aggregation interval n is found by averaging over the K periods of the aggregation interval.(16)ffpn=1Kn∑k=1Knfpk,n.Figure 8
Zero-crossing algorithm (ZC algorithm). Illustration of the
zero-crossing positions relative to a 100 ms time
interval.Due to the time misalignment, this kind of estimation causes an error that
depends on both the phaseθt at the interval borders and the derivative of the
instantaneous frequency finst within the interval. The error thereby scales linearly
with the latter. In the worst case, a full period in the average does not
actually belong to the report interval n. For one-second records, this equals 2% of the
involved data and is usually tolerated. For 100 ms reports this
figure increases to 20% and causes intolerable deviations during
transient frequency events.To address this challenge, we introduced the partial period method (pp) in
[21] that performs a local
resampling of the data to the regular report intervals. Contributions of the
border periods (i.e., periods 1 and 6 in Figure 8) are partially included in the average according to
their temporal overlap with the report interval marked in blue as follows:(17)fppn=1Kn−2+w1+wK∑k=1Knwkfpk,n,where K is the number of periods that are in touch with the
aggregation interval n and wk are weights for the individual periods.The first weight is given by(18)w1=tzc1−trepn−1·fp1,and the last weight is given by(19)wK=trepn−tzcK−1·fpK,while all periods in between are equally weighted with wk=2…K−1=1.While the pp-method eliminates the majority of shift-related estimation
errors, subtle deviations remain. The main reason is that the approach
assumes thatfinst is constant during a period. This condition is hardly
satisfied for any real signal. Nevertheless, the estimation error in dynamic
conditions is drastically improved compared to the fp-method [21]. The partial period method,
therefore, practically solves the temporal alignment problem.


## 2.2.4. Accuracy of Frequency Measurements

Since the individual period, measurements are only as good as the underlying
time reference, we use the GPS receivers’ pulse per second (PPS)
signal as a local frequency standard. The error in frequencyΔf relates to the uncertainty in the time Δt via(20)Δf=−T=1fnom−2Δt.The edge-to-edge error of our receiver (Garmin GPS-18LVC) and PPS acquisition
is less than 50 ns, which translates to a worst-case frequency error
of 125µHz in the 50 Hz case.
Experimentally, we further verify that the frequency estimation error in
steady-state conditions over a 10-minute interval is only
26.7 µHz at 50 Hz within a
95.4% confidence interval [45]. In dynamic conditions, however, the estimation errors are
dominated by the principles of the zero-crossing method. A detailed analysis
of the estimation errors in the time domain can be found in [21]. For the frequency ramps defined in
the standard for Phasor Measurement Units (IEEE 60255-118-1) we found the
error limit of our method to be ±0.6mHz.


## 2.2.5. Influence of Signal Disturbances

As shown in II A 3. and Figure6, the
observed PSDs can be reproduced using a stochastic process composed of an
Ornstein-Uhlenbeck process (including white noise) and added correlated
fluctuations, where the latter lead to the increasing spectrum on a short
time scale. We thus investigate how our processing chain deals with
disturbances in the signal. By doing so, we also ensure that no artificial
signal is introduced into the system.From reviewing the signal chain, disturbances superimposed onut are the most likely source of the observed fluctuations,
since adding a random signal eut to ut will impact how well the period boundaries can be
determined. The relation between such voltage noise eu and noise in the period duration (time) measurement et is given by the gradient of ut at the zero-crossing as follows:(21)etk=euddtut=tZCk=2πAf0euk.Voltage noise, therefore, translates linearly to noise in period duration
(phase noise), scaled by power-grid frequency and the signal amplitude.
However, ifeu is not random, the influence might vanish. This is
especially the case for harmonics, where eu is periodic in multiples of f.To gain an impression of the influence of noise, we conduct the following
computer experiment. We generate a sinusoidal signal, superimpose it with
white Gaussian noise and process it with the estimation algorithm to obtain
a series of frequency readingsfn. We then determine the Welch PSD estimates Sf of this series as carried out for the measured frequency
time series. Results are shown in Figure 9 from which we make three major observations:(1)
The curves for 30 dB and 60 dB signal-to-noise
ratio (SNR) are similar in shape but vertically shifted by three
orders of magnitude. The shift corresponds exactly to the
30 dB difference in noise excitation levels. Therefore,
the noise levels inf and ut can be directly related.(2)
Sf rises approximately proportional to quadratic
with t−2 which is expected from frequency demodulator
theory [47]. Since the ZC
algorithm averages K period measurements, we would assume that the
total noise amplitude of f scales with 1/K. However, due to the increase in noise
density, the total noise amplitude in f increases linearly with K. Changing the resolution from 1 s to
100 ms translates to a tenfold increase in noise
amplitude (in the 60dB case, and σ=0.48mHzvs.4.9mHz in the 30dB case). Hence, the data with higher temporal
resolution contains more noise in absolute terms.(3)
For smaller values off we observe a slight flattening of S indicating a possible convergence to a
constant value. The shallow part indicates some feed-through of
white phase noise, which is even more pronounced for the
full-period method (shown in Supplementary Material). We find
the source of this effect in the violation of the previously
introduced assumption. The instantaneous frequency is not
constant during each period. Compared to the noise contributions
from the high-frequency end of the spectrum, the feed-through
has a minuscule effect on the results.Figure 9
Power spectral density off for the partial period zero-crossing method. Input
is a sine wave with the additional white noise of different levels.
PSD estimates saturate for a long period duration and rise nearly
quadratic (solid black line) towards shorter periods. For
60 dB SNR the dashed line indicates the maximum level of Sf at
0.01 mHz2·Hz−1.


## 2.2.6. Quantification of Disturbance Influence

So far, we have characterized the effects of noise on frequency estimation.
We still need to quantify what noise contribution is to be expected in an
actual measurement environment so that we can establish a boundary of trust
for our results. For this, we need to (a) analyze the distortions in the
real-world setting and (b) characterize the susceptibility of the
algorithm.To approach (a), we investigate the properties ofut. In Figure 10, we
plot an excerpt of the waveform and the spectrum of the measured voltage.
The majority of the disturbance is in the nonsinusoidal shape of the
waveform and, therefore, in the harmonic content.Figure 10
Captured voltage waveform in Stockholm. Top: example of the captured
voltage waveformut. Data: Stockholm 9th of May 2019,
20:35. Bottom: power spectrum of one hour of data (Welch estimate,
resolution 1 Hz, Hann window). Signal-to-noise ratio (SNR)
including distortion: 29.6 dB. This is worst case SNR
observed for EDR measurements.From the spectrum ofut we can determine the power contained in the disturbances
to obtain an SNR for the measurements (details in Supplementary Material).
Throughout the campaigns, we find a minimum SNR of ≈30dB. However, if we selectively suppress the contribution of
the first ten harmonics, the minimum SNR rises to ≈60dB.Additionally to the above, short-term transients caused by local switch
operations will occasionally occur inut. From the time-domain analysis of the algorithm in [21], we know that such transients will
cause large (usually 0.1 Hz and above) single-point frequency
deviations. Since a pulse in the Dirac sense translates to a constant in the
frequency domain, those transients might cause slightly elevated Sf levels. However, most of these single-point deviations are
removed from the data via an outlier filter prior to our investigation.For (b), we already know from the processing chain, that components above
50 Hz are rejected by the initial low-pass filtering. It is further
known that static waveform deformations have no impact on the estimation
[21]. We are therefore convinced
that omitting the major harmonic contributions from the SNR estimation is
justified (further discussion in Supplementary Material).Consequently, the expected contribution of the ZC algorithm to the observed
PSD estimates follows the 60 dB curve in Figure9. This means that in all cases except Oldenburg the
algorithm’s contribution explains less than one-tenth of the observed Sf. In the case of Oldenburg, the factor is ≈1/4, while in most cases, it is in the range of ≈1/100. Hence, we are confident that the observed increases of Sf are not a measurement artifact but a signal property.


## 3. Discussion

Overall, we uncovered complex behavior in power-grid frequency recordings from
high-quality measurements. Analyzing the data, we made three key contributions:(i)
First, we have seen that the PSD decays approximately following a power
law until a time scale of 1 second, thereby displaying a typical
complex system property.(ii)
Second, we noted distinct dynamics on the shorter time scale with
increasing spectral amplitude, multifractal behavior and a change from
positive to negative correlations.(iii)
Third, we showed how key properties of the spectrum might be reproduced
by describing the power grid as a complex system following an
Ornstein–Uhlenbeck process with added fluctuations, motivated
mostly by basic principles.Finally, we described the acquisition processing chain used to acquire the frequency
data. The challenges in terms of temporal alignment and measurement noise arising
when moving to the high temporal resolution were outlined, and we discussed how the
EDR-derived data deals with these challenges.What does this mean for our understanding of complex systems and power system
operation and modeling? To enable any of our observations, it was critical to have
access to high-quality, highly resolved data. This stresses the need to monitor
power systems closely and also to make these data, including simultaneous and
independent measurements in multiple synchronous areas, available for broad analysis
and comparison.Furthermore, we have seen that fluctuations, especially on short time scales, are
complex and not simply described as Gaussian white noise but instead seem to follow
an anticorrelated and multifractal process. Possible explanations for the emerging
anticorrelation range from effects in the distribution grid, the self-regulating
effect of frequency-dependent loads [40], and
layered network dynamics [48] to
superstatistics. Regardless of the cause, this complexity should be implemented
explicitly by any modeler when simulating and analyzing power systems, especially as
an increasing share of renewables and reductions of the total inertia will likely
amplify nonstandard statistics only further. Only when we fully model and understand
fluctuations can we design efficient and effective control strategies to ensure
stability.Further research is necessary to fully explain what particular elements in power
systems give rise to multifractal phenomena and investigate their potential
connection to superstatistics. On a similar note, complex spectra are observed in
other systems [49], so we suspect the
existence of more nonstandard spectra in other time series.



---

*Source: 2657039-2023-03-06.xml*