# Sneak Preview: PowerDynamics.jl
\--
An Open-Source library for analyzing dynamic stability in power grids with
high shares of renewable energy

Tim Kittel, Sabine Auer and Christina Horn  Potsdam Institute for Climate
Impact Research (PIK)
Member of the Leibniz Association
P.O. Box 60 12 03
D-14412 Potsdam
Germany
Email<EMAIL_ADDRESS>
###### Abstract

PowerDynamics.jl is an Open-Source library for dynamic power grid modeling
built in the latest scientific programming language, Julia. It provides all
the tools necessary to analyze the dynamical stability of power grids with
high share of renewable energy. In contrast to conventional tools, it makes
full use of the simplicity and generality that Julia combines with highly-
optimized just-in-time compiled Code. Additionally, its ecosystem provides
DifferentialEquations.jl, a high-performance library for solving differential
equations with built-in solvers and interfaces to industrial grade solvers
like Sundials. PowerDynamics.jl provides a multitude of dynamics for different
node/bus-types, e.g. rotating masses, droop-control in inverters, and is able
to explicitly model time delays of inverters. Furthermore, it includes
realistic models of fluctuations from renewable energy sources. In this paper,
we demonstrate how to use PowerDynamics.jl for the IEEE 14-bus distribution
grid feeder.

##  I Introduction

Dynamic stability analysis of power grids concerns the investigation of
transient stability and power quality. Especially with more intermittent
renewable generation sudden changes in power followed by the transient
reaction of frequency and voltage dynamics may present a challenge to power
system stability. At the same time, fluctuations also challenge power quality,
which refers to whether voltage and frequency stay within safety bounds and
whether the waveform has undesired harmonic distortions. Furthermore, the wide
deployment of power electronics, witch RES are connected to the grid,
introduces new features such as delayed reaction into future power grids.

These challenges have been the topic of recent research on the stability of
renewable power grids [1, 2, 3]. However, the usage of state-of-the-art Open-
Source software repeatedly posed limitations to the possible modeling scope.
One reason is the lack of solvers for the different types of dynamics which
need to be described a stochastic differential equations, delay differential
equations and differential algebraic equations. At the same time to network
size was constrained by computational performance. This means a realistic
model implementation and hence reliable stability analysis has been
constrained by the available OS software framework.

In this conference paper, we give a sneak preview to PowerDynamics.jl, an new
Open-Source approach to dynamic power grid modeling111This paper is meant as a
sneak preview, as the planned publication date is Oct 15, 2018, i.e. before
the Wind Integration Workshop 2018 will take place.. While there are multiple
modeling suites available they are either proprietary (e.g. PSS® [4]) or need
proprietary tools to be used (e.g. PSAT [5, 6]). With PowerDynamics.jl, we aim
to provide a modeling framework that tackles all the needs a modeler for power
grids with high shares of renewables to execute dynamic analyses and the first
publication will provide tools for time-domain modeling.

At the moment, PowerDynamics.jl is under high development, with the basic
software architecture already being developed and implemented. We are
currently adding new types of dynamics and controls to represent different
buses, e.g. droop control, multiple inverter types, higher-order
representations of synchronous machines, and in particular intermittency
representations of renewable energy sources. A second focus of effort goes
into usability, in order to improve the experience of modelers and decrease
their time they have to concern with software-related problems so they can
focus on their modeling instead. Our vision is to provide a modular, fully
flexible library where modelers can decide whether they (a) want to simply use
repdefined bus models to analyze a power grid quickly, (b) want to into all
the detail of modeling each bus type, or (c) want to find their right position
in-between (a) and (b). All of that is done in the context of an Open-Source
context in order to make it available to and invitre contributions from
research and industry.

PowerDynamics.jl is written in Julia [7], a new, fast developing programming
language targeting the scientific computing community. With its recent release
to Julia 1.0, it has matured enough to become a powerful player competing with
more traditional programming languages like Python and C in the modeling
world. In the wider context, there is a current movement within the Open
Energy Modeling community [8] to switch to Julia and unite modeling efforts.
With PowerDynamics.jl we are aiming to support these efforts by contributing
to dynamic power system analysis.

The rest of the paper is structured as follows: in Section II ‘‘Julia for
Scientific Modeling’’ we give a detailed reasoning why we find Julia to be the
correct choice as the programming language for our modeling effort; in Section
III ‘‘Mathematics of PowerDynamics.jl’’ we introduce the basic mathematical
notions we use in PowerDynamics.jl; in Section IV ‘‘The IEEE 14-bus
distribution grid feeder’’ we provide the details of the example system used
in this paper; in Section V ‘‘Implementation in PowerDynamics.jl’’ we show how
to implement the IEEE 14-bus feeder in PowerDynamics.jl; in Section VI
‘‘Modeling Results’’ we present the results for our example model; and in
Section VII ‘‘Conclusion & Outlook’’ we conclude the paper.

### Convention & Notation

Within this paper, we take up the convention of PowerDynamics.jl where minor
letters are used as variables for dynamic variables and capital letters denote
(constant) parameters and all variables are in defined in the co-rotating
frame.

Further, we take $u_{a}$ to to be the complex voltage, $v_{a}=|u_{a}|$ to be
the voltage magnitude, $\varphi_{a}=\arg(u_{a})$ to be the voltage angle,
$i_{a}=\sum_{b}Y^{L}_{ab}u_{b}$ to be the complex current, $s_{a}=u_{a}\cdot
i_{a}^{*}$ to be the complex power, $p_{a}=\text{Re}(s_{a})$ the active power
and $q_{a}=\text{Im}(s_{a})$ to be the reactive power of the

$a^{\text{th}}bus.Y^{L}istheadmittanceLaplacian.Theimaginaryelementisdenotedasj=\sqrt{-1}.\par$

##  II Julia for Scientific Modeling

The advantages of Julia for scientific modeling can be summarized as (adapted
from [9]):

  * •

Julia is fast! It was designed from the beginning as a high-performance
language. Using the just-in-time (JIT) compilation method, it produces
efficient native code for multiple platforms.

  * •

Julia allows for rapid development! It uses dynamic typing, so it is easy to
use, feels like a scripting language and has good interactive use without
losing any performance. That way, it saves the developer a lot of time.

  * •

Julia is technical! It is developed for numerical computing, the syntax is
focusing on precise mathematics and many datatypes and even parallelism is
available out of the box.

  * •

Julia is general! Using multiple dispatch as the fundamental paradigm, it is
easy to express object-oriented and functional programming patterns at the
same time.

  * •

Julia is composable! Julia has been designed such that independent packages
work well together. Hence, we can use matrices of unit quantities or
differentiation with sparse matrices without having these types ever defined
explicitly before.

  * •

Julia has a growing ecosystem! The aforementioned advantages attract a growing
community222The community growth recently reached the point where the
organizers had to cut off ticket sales for this year’s JuliaCon at some point.
and with that a growing ecosystem of Julia packages. Particularly important
for modeling is the DifferentialEquations.jl library which provides high-
performance solvers and interfaces to industrial grade solvers like SUNDIALS
[10] for solving multiple kinds of differential equations (ordinary,
algebraic, delayed, stochastic).

  * •

Julia allows for metaprogramming! This means, at execution time, the source
code is available as data and one can easily modify it. Furthermore, one can
even add own syntax to Julia. Within PowerDynamics.jl, we make heavy use of
metaprogramming to make it as simple as possible for a user to implement
her/his own kind of bus model (see Section V-A). Also, it allows to write very
simple code that can automatically be modified to highly optimized code.

Of course, this is just a simply overview and we would strongly recommend the
reader to try out Julia and convince
herself/himself333https://julialang.org/learning/.

##  III Mathematics of PowerDynamics.jl

Within this paper, we take the approach to describe the dynamics of a power
grid by a (set of) semi-explicit differential algebraic equation(s). We
explicitly distinguish between the complex voltages $u_{a}$ (of the

$a^{\text{th}}bus),whoseinformationistransmittedbetweenbusesdirectlythroughthecomplexcurrenti,andinternalvariablesx_{ab}(theb^{\text{th}}internalvariableofthea^{\text{th}}bus),thatexistlocallyatanode.Anexampleforaninternalvariablewouldbethefrequencyofasynchronousmachine.Whilethephaseangle\varphi
istransmitteddirectlyviathecurrent,itsderivativeistheanglevelocitywhosedynamicsisdefinedlocallybythetherotatingmass.Weareawarethatamultitudeofapproachesarepossibletowritedownthedynamicsandwehavechosenthisparticularoneasitcanbeeasilytranslatedintosourcecode.Writingthisideainformulasyields\lx@equationgroup@subnumbering@begin\begin{aligned}
i_{a}=\sum_{c}Y^{L}_{ac}\cdot u_{c},\\\
m^{u}_{a}\dot{u}_{a}=f_{a}(u_{a},x_{a},i_{a}),\\\
m^{x}_{ab}\dot{x}_{ab}=g_{ab}(u_{a},x_{a},i_{a}),\end{aligned}\lx@equationgroup@subnumbering@end
wherem^{u}/m^{x}arethemassesforthevoltages/internalvariablesrespectively.Also,f_{a}encodesthevoltagedynamicsofthea^{\text{th}}busandg_{ab}thedynamicsoftheb^{\text{th}}internalvariableofthea^{\text{th}}bus.Notetheconventionofusingminorlettersfordynamicalvariables.\par
Y^{L}istheadmittanceLaplacian.SohavingtheadmittancesencodedasY_{ab}whereY_{ab}=0ifthereisnoline,orY_{ab}withacomplexvalue(with\text{Re}(Y_{ab})\geq
0)forexistinglines,thentheadmittanceLaplacianisdefinedas\begin{aligned}
Y^{L}_{ab}=\delta_{ab}\sum_{c}Y_{ac}-Y_{ab},\end{aligned}where\delta_{ab}isKronecker-\delta.Particularexampleswouldbeanalgebraicslackbusasthea^{\text{th}}buswithafixedcomplexvoltageU_{a}.Therearenointernalvariablesandthevoltagemassism_{a}^{u}=0(asitsanalgebraicequation),soitcansimplybewrittenas\begin{aligned}
0\cdot\dot{u}_{a}=f_{a}(u_{a},x_{a},i_{a})=u_{a}-U_{a},\end{aligned}wherewehaveexplicitlynotremoved\dot{u}inordertokeepastronganalogythetheimplementationshownlater.Similarly,wecantreatanalgebraicloadatthea^{\text{th}}buswithaconstantcomplexpower-
S_{a}flowingout.Again,therearenointernalvariablesandthevoltagemassis0soitreducesto\begin{aligned}
0\cdot\dot{u}_{a}=f_{a}(u_{a},x_{a},i_{a})=s_{a}-S_{a}=u_{a}\cdot
i_{a}^{*}-S_{a}.\end{aligned}Asynchronousmachineasthea^{\text{th}}buscanberepresentedbytheSwingequation(or2^{\text{nd}}-ordersynchronousmachinemodel)\cite[cite]{[\@@bibref{}{sauer1998power,kundur1994power}{}{}]}withaproducedactivepowerP_{a},thedampingconstantD_{a},theratedfrequency\Omega
andtheinertiaH_{a}.Inthiscase,thevoltagemassis1,thereisoneinternalvariable,thefrequency,x_{a1}=\omega_{a}withamassof1becauseitisdynamicandnotanalgebraicconstraint.Hence,theequationsreduceto\lx@equationgroup@subnumbering@begin\begin{aligned}
&\dot{u}_{a}=f_{a}(u_{a},(\omega_{a}),i_{a})=u_{a}j\omega_{a},\\\
&\begin{aligned}
\dot{x}_{a1}=\dot{\omega}_{a}&=g_{a1}(u_{a},(\omega_{a}),i_{a})\\\
&=\frac{2\pi\Omega}{H_{a}}\left(P_{a}-D_{a}\omega_{a}-\text{Re}(u_{a}\cdot
i_{a}^{*})\right)\end{aligned}\end{aligned}\lx@equationgroup@subnumbering@end
wherewewouldliketoremindthereaderofjbeingusedastheimaginaryelementinthispaper.Notethatwithu_{a}=e^{j\varphi_{a}}\
\implies\
du_{a}=ju_{a}d\varphi_{a},\lx@cref{creftypecap~refnum}{eq:swing-u}reducesto\dot{\varphi}_{a}=\omega_{a}andonerecoverstheusualversionintermsofthevoltageangle.\par
Inthefollowingsection,wewillintroducetheexamplesysteminordertodemonstrateafterwards,howthesekindofdynamicsarethendescribedinPowerDynamics.jl{}.\par$

##  IV The IEEE 14-bus distribution grid feeder

Figure 1: Technical representation of the IEEE 14-bus test system, extracted
from [13]. It contains five synchronous machines of which two are generators.
The one at bus 2 is taken as the slack bus for the modeling in this paper.
Furthermore, there are 8 loads (buses 4, 5, 8, 9, 10, 11, 12, 13, and 14). Bus
7 is for the representation of the transformer.

The IEEE 14-bus system is a representation of a medium-voltage distribution
grid. A graphical representation is in Figure 1, which is extracted from [13].
There are two generators, where the one at bus 2 is the slack bus, three
synchronous compensators, and eight (complex) loads. Bus 7 is for the
representation of the transformer only.

The generator at bus 1 and the synchronous compensators will be modeled with
the swing equation, the slack bus 2 with the algebraic constraint introduced
in Section III and the loads with the algebraic constraints for complex loads
(see Section III also). The bus parameters haven been extracted from [13] and
are presented in Table I. The resistance $R$ and reactance $X$ for the lines
(taken from [13] as well) are summarized in Table II.

Having a description of the power grid in question enables us to have a look
at the implementation in PowerDynamics.jl now.

TABLE I: Bus parameters of the IEEE 14-bus test system. INERTIA UNIT IS? Bus | Type | $U\,[p.u.]$ | $P\,[p.u.]$ | $Q\,[p.u.]$ | $D\,[p.u.]$ | $H[p.u.\cdot s]$
---|---|---|---|---|---|---
1 | Generator |  | 2.32 |  | 2 | 5.148
2 | Slack | 1 |  |  |  |
3 | Syn. Comp. |  | -0.942 |  | 2 | 6.54
4 | Load |  | -0.478 | 0 |  |
5 | Load |  | -0.076 | -0.016 |  |
6 | Syn. Comp. |  | -0.112 |  | 2 | 5.06
7 | Load |  | 0 | 0 |  |
8 | Syn. Comp. |  | 0 |  | 2 | 5.06
9 | Load |  | -0.295 | -0.166 |  |
10 | Load |  | -0.09 | -0.058 |  |
11 | Load |  | -0.035 | -0.018 |  |
12 | Load |  | -0.061 | -0.016 |  |
13 | Load |  | -0.135 | -0.058 |  |
14 | Load |  | -0.149 | -0.05 |  |
TABLE II: Line parameters of the IEEE 14-bus test system. line number | from bus | to bus | $R\,[p.u.]$ | $X\,[p.u.]$
---|---|---|---|---
1 | 1 | 2 | 0.01938 | 0.05917
2 | 1 | 5 | 0.05403 | 0.22304
3 | 2 | 3 | 0.04699 | 0.19797
4 | 2 | 4 | 0.05811 | 0.17632
5 | 2 | 5 | 0.05695 | 0.17388
6 | 3 | 4 | 0.06701 | 0.17103
7 | 4 | 5 | 0.01335 | 0.04211
8 | 4 | 7 | 0.0 | 0.20912
9 | 4 | 9 | 0.0 | 0.55618
10 | 5 | 6 | 0.0 | 0.25202
11 | 6 | 11 | 0.09498 | 0.1989
12 | 6 | 12 | 0.12291 | 0.25581
13 | 6 | 13 | 0.06615 | 0.13027
14 | 7 | 8 | 0.0 | 0.17615
15 | 7 | 9 | 0.0 | 0.11001
16 | 9 | 10 | 0.03181 | 0.0845
17 | 9 | 14 | 0.12711 | 0.27038
18 | 10 | 11 | 0.08205 | 0.19207
19 | 12 | 13 | 0.22092 | 0.19988
20 | 13 | 14 | 0.17093 | 0.34802

##  V Implementation in PowerDynamics.jl

In this section, we will show first how to implement the dynamics for the
different buses presented in Section III in PowerDynamics.jl and then how
create a power grid model from these bus dynamics definitions, in this case
the IEEE 14-bus system.

###  V-A Defining the dynamics

The implementation of such a power grid in PowerDynamics.jl is strongly based
on the $\mathtt{@DynamicNode}$ macro. A macro is a metaprogramming function
that takes part of the source codes and modifies it at parsing time before
giving it to the compiler, that creates the actual bitcode. In case of
$\mathtt{@DynamicNode}$, this is used to hide all the complicated internals of
PowerDynamics.jl and make the definition of a type of dynamics as easy as
possible. The following examples are already provided by default in
PowerDynamics.jl but using them as an example, any kind of dynamics can be
implemented.

Using $\mathtt{@DynamicNode}$, the slack bus constraint from Section III is
implemented as

⬇

1@DynamicNode SlackAlgebraic(U) <: OrdinaryNodeDynamicsWithMass(m_u=false)
begin

2end [] begin

3 du = u - U

4end

where line breaks and line numbers have been added for easier orientation. The
left part on the first line states the name of the new type
$\mathtt{SlackAlgebraic}$ and the parameter name $\mathtt{U}$. Separated by
the subtyping operator ($\mathtt{<:}$) comes the statement that it can be
represented by an ordinary differential equation with a mass term, where the
mass for the voltage `m_u` is set to $\mathtt{false}$, meaning $m^{u}=0$ from
Section III as it is a constraint. The `[]` in line 2 states that there are no
internal variables. And line 3 simply provides the formula for the dynamics
and assigns it to $\mathtt{du}$. As `m_u = false` is given, this simply
translates it to a constraint.

In a similar manner, the load is defined as an algebraic PQ-constraint

⬇

1@DynamicNode PQAlgebraic(S) <: OrdinaryNodeDynamicsWithMass(m_u=false) begin

2end [] begin

3 s = u*conj(i)

4 du = S - s

5end

In this case, an additional line was added to calculate the outflowing complex
power $\mathtt{s}$ from the voltage $\mathtt{u}$ and the current $\mathtt{i}$
and then the definition of the constraint is written in line 4.

Concerning the generator and the synchronous compensators, we implemented the
swing equation as

⬇

1@DynamicNode SwingEq(H, P, D, SOmega) <: OrdinaryNodeDynamics() begin

2 @assert D > 0 "damping (D) should be >0"

3 @assert H > 0 "inertia (H) should be >0"

4 SOmega_H = SOmega * 2pi / H

5end [[Somega, dSomega]] begin

6 p = real(u * conj(i))

7 du = u * im * Somega

8 dSomega = (P - D*Somega - p)*SOmega_H

9end

Again, the left part of the first line provides the name of the new type and
the name of the parameters. Lines 2 to 4 are code that should be run only
once. In this case these are consistency checks, that the damping and inertia
are positive, and the reduction of the rated frequency $\mathtt{\Omega}$ and
the inertia $\mathtt{H}$ to a single variable $\mathtt{\Omega\\_H}$. In line
5, the variable name of the internal variable $\mathtt{\omega}$ and the name
for its derivative $\mathtt{d\omega}$ are given. Finally, lines 6 to 8
implement Sections III and III by simply writing down the mathematical terms.

###  V-B Instantiating the power grid model

In order to create the grid model, we first have to instantiate the bus models
simply by calling them with the corresponding parameter values from Table I,
e.g.:

⬇

1SwingEq(H=5.148, P=2.32, D=2, SOmega=50) # for bus 1

2PQAlgebraic(S=-0.295-0.166im) # for bus 9

Within the actual code we simply loaded the data from a $\mathtt{.csv}$ file
into and automatized this instantiation777The full source code is not part of
this paper as it would be too long, but it will be published along with
PowerDynamics.jl.. The instantiated bus models should then be saved in an
array called e.g. $\mathtt{nodes}$. Similarly, the admittance Laplacian should
be generated from the line data in Table II and saved in a matrix called e.g.
$\mathtt{LY}$. The actual grid model instatiation is then simply one line
where the model is saved in the variable $\mathtt{g}$:

⬇

1g = GridDynamics(node_list, LY)

Now, $\mathtt{g}$ contains all the information of the power grid and in the
following section we will show how to solve it.

##  VI Modeling Results

(a)

(b)

Figure 2: The two subfigures show the outflowing power (top) at each node and
the angular frequency (bottom) for the buses modeled by the swing equation (1,
3, 6, 8) in for the two perturbation scenarios: (left) frequency perturbation
at bus 1 and (right) line tripping of line 2 (between bus 1 and 5). The letter
in bracets refers to the modeling type, i.e. $G$ for generator / synchronous
compensator as swing equation, $S$ for the slack bus, and $L$ for a load as
algebraic PQ-constraint.

Within this section, we will analyze two simple cases: (a) a frequency
perturbation at the largest generator at bus 1 and (b) a line tripping event
at of line 2 (between bus 1 and 5) (cmp. Figure 1).

Before we find the normal point of operation for the power grid, i.e. the
fixed point of Sections III, III and III or synchronous state for the IEEE
14-bus system. For that, we use the grid model $\mathtt{g}$ generated in the
previous Section V and use the function provided by PowerDynamics.jl:

⬇

1fp = operationpoint(g, ones(SystemSize(g)))

where $\mathtt{ones(SystemSize(g))}$ is a vector of the correct length for the
initial condition of the fixed point search. $\mathtt{fp}$ is now a
$\mathtt{State}$ object that we can use as initial condition for the solving
the differential equations corresponding to the power grid.

###  VI-A Frequency perturbation

In order to model a frequency perturbation, one can simply take a copy of the
fixed point as found before and adjust the initial frequency value:

⬇

1x0 = copy(fp)

2x0[1, :int, 1] += 0.2

The second line can be read as ‘‘add 0.2 to the 1${}^{\text{st}}$ internal
variable of the 1${}^{\text{st}}$ node’’ which is the frequency $\omega$ as
there is only one internal variable. Note that the first $\mathtt{1}$ refers
to the node and then second $\mathtt{1}$ to the internal variable counter. The
power grid with the initial condition $\mathtt{x0}$ can then be solved for a
time span of 0.5 seconds by calling:

⬇

1sol = solve(g, x0, (0.0,.5))

The solution is shown in Figure 2a. It shows that the system is very stable
against frequency perturbations. The actual dynamics is not so exciting as the
system is very stable. Please note that this system was taken as an example to
present on how easily one can model a power grid using PowerDynamics.jl, not
to find new exciting dynamics.

###  VI-B Line tripping

To show some more dynamic behavior we simulated a line tripping as well. We
model this effect by taking the operation point of the full power grid
($\mathtt{fp}$) as initial condition but defining a new admittance Laplacian
where line 2 (between bus 1 and 5, see Tables II and 1) has been taken out
(i.e. the admittance is set to $0$). Running the model with this new Laplacian
yields Figure 2b.

In the frequency plot we identify how the frequency of bus 1, where the line
tripping happened, compensates for the momentarily excess power at the bus.
The lacking power in the rest of the grid is matched by the synchronous
compensators whose frequency decreases in turn. Note that the angular
frequency $\omega$ is shown, so with a division by $2\pi$ the maximal
frequency deviation $f$ is $\approx 0.05\,$Hz. After about one second, the
system recovers to the normal state of operation.

##  VII Conclusion & Outlook

Within this paper, we have seen how one one can use the Open-Source library
PowerDynamics.jl in order to model the dynamics of a power grid with just a
few lines of code. We have seen how the fundamental mathematical equations
(given in Section III) translate to source code that reads exactly the same.
We employed PowerDynamics.jl for the IEEE 14-bus distribution grid feeder in
order to demonstrate how one can easily simulate faults and analyze the
transient reaction of the power grid dynamics. As example scenarios we used a
frequency perturbation and a line tripping.

Finally, this paper is really just a sneak preview with a simple example the
publication of PowerDynamics.jl is planned for October 15, 2018. By then, we
will have added more inverter control schemes and stochastic descriptions of
intermittency due to renewable energy sources.

## Acknowledgment

This paper was presented at the 19th Wind Integration Workshop and published
in the workshop’s proceedings.

We would like to thank the German Academic Exchange Service for the
opportunity to participate at the Wind Integration Workshop 2018 in Stockholm
via the funding program ‘‘Kongressreisen 2018’’. This paper is based on work
developed within the Climate-KIC Pathfinder project ‘‘elena -- electricity
network analysis’’ funded by the European Institute of Innovation &
Technology. This work was conducted in the framework of the Complex Energy
Networks research group at the Potsdam Institute for Climate Impact Research.

We would like to thank Frank Hellmann and Paul Schultz for the discussions on
structuring an Open-Source library for dynamic power grid modeling.

## References

  * [1] S. Auer, F. Hellmann, M. Krause, and J. Kurths, ‘‘Stability of synchrony against local intermittent fluctuations in tree-like power grids,’’ _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 27, no. 12, p. 127003, 2017.
  * [2] S. Auer, ‘‘The stability and control of power grids with high renewable energy share,’’ Ph.D. dissertation, 2018.
  * [3] B. Schäfer, C. Grabow, S. Auer, J. Kurths, D. Witthaut, and M. Timme, ‘‘Taming instabilities in power grid networks by decentralized control,’’ _The European Physical Journal Special Topics_ , vol. 225, no. 3, pp. 569--582, 2016.
  * [4] Siemens AG, ‘‘PSS® power system simulation and modeling software,’’ https://www.digsilent.de/en/powerfactory.html, accessed: 2018-08-28.
  * [5] F. Milano, ‘‘PSAT,’’ http://faraday1.ucd.ie/psat.html, accessed: 2018-08-28.
  * [6] F. Milano, L. Vanfretti, and J. C. Morataya, ‘‘An open source power system virtual laboratory: The psat case and experience,’’ _IEEE Transactions on Education_ , vol. 51, no. 1, pp. 17--23, 2008.
  * [7] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, ‘‘Julia: A fresh approach to numerical computing,’’ _SIAM review_ , vol. 59, no. 1, pp. 65--98, 2017.
  * [8] J. C. Richstein, ‘‘open energy modeling initiative,’’ http://www.openmod-initiative.org/, accessed: 2018-08-28.
  * [9] JuliaLang, ‘‘The Julia Language,’’ https://julialang.org/, accessed: 2018-08-28.
  * [10] A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker, and C. S. Woodward, ‘‘SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers,’’ _ACM Transactions on Mathematical Software (TOMS)_ , vol. 31, no. 3, pp. 363--396, 2005.
  * [11] P. W. Sauer and M. Pai, ‘‘Power system dynamics and stability,’’ _Urbana_ , 1998\.
  * [12] P. Kundur, N. J. Balu, and M. G. Lauby, _Power system stability and control_. McGraw-hill New York, 1994, vol. 7.
  * [13] S. K. M. Kodsi and C. A. Canizares, ‘‘Modeling and simulation of ieee 14 bus system with facts controllers,’’ _Technical Report 2003-3_ , 2003.


# Uncovering Load-Altering Attacks Against $N-1$ Secure Power Grids: A Rare-
Event Sampling Approach

Maldon Patrice Goodridge, Subhash Lakshminarayana , and Alessandro Zocca M.P.
Goodridge is with Global Development Initiatives, UK, (email:
gdi.uklimited@gmail.com). S. Lakshminarayana (Corresponding author) is with
the School of Engineering, University of Warwick, Coventry, UK, CV47AL (email:
subhash.lakshminarayana@warwick.ac.uk). A. Zocca is with the Department of
Mathematics, Vrije Universiteit Amsterdam, NL (email: a.zocca@vu.nl). A
preliminary version of the paper was published at PMAPS 2022 [1].

###### Abstract

Load-altering attacks targetting a large number of IoT-based high-wattage
devices (e.g., smart electric vehicle charging stations) can lead to serious
disruptions of power grid operations. In this work, we aim to uncover
spatiotemporal characteristics of LAAs that can lead to serious impact. The
problem is challenging since existing protection measures such as $N-1$
security ensures that the power grid is naturally resilient to load changes.
Thus, strategically injected load perturbations that lead to network failure
can be regarded as _rare events_. To this end, we adopt a rare-event sampling
approach to uncover LAAs distributed temporally and spatially across the power
network. The key advantage of this sampling method is the ability of sampling
efficiently from multi-modal conditional distributions with disconnected
support. Furthermore, we systematically compare the impacts of static (one-
time manipulation of demand) and dynamic (attack over multiple time periods)
LAAs. We perform extensive simulations using benchmark IEEE test simulations.
The results show (i) the superiority and the need for rare-event sampling in
the context of uncovering LAAs as compared to other sampling methodologies,
(ii) statistical analysis of attack characteristics and impacts of static and
dynamic LAAs, and (iii) cascade sizes (due to LAA) for different network sizes
and load conditions.

##  I Introduction

The growing penetration of internet-of-things (IoT) enabled high-wattage
appliances (e.g., electric vehicle charging stations) may create new malicious
attack opportunities against power grids. _Load-altering attacks_ (LAAs) are
large-scale load fluctuations that can be generated by changing the
operational settings of IoT-based load devices under a botnet-type attack.
Such LAAs can pose a serious threat to power grid operations [2, 3, 4, 5, 6,
7]. LAAs can be classified into two categories: (i) _static_ LAAs (S-LAAs),
which refers to a sudden one-time manipulation of power grid demand, and (ii)
_dynamic_ LAAs (D-LAAs), which refers to a series of load alterations over
time.

Analyzing the impact of LAAs on power grid operations is an important step
toward developing measures to enhance system resilience. It was shown [3, 4]
that the imbalance caused due to the mismatch between the supply and demand
under S-LAAs can lead to unsafe frequency excursions causing line outages,
generator trips, and/or increase the grid’s operational costs. D-LAAs, on the
other hand, do not have an instantaneous impact, but however, can potentially
destabilize the power grid’s frequency control loop over a period of time [6].
In [5], the authors analyzed the LAA impact under more realistic operational
regimes including protection measures such as under-frequency load shedding,
over-frequency protections, etc., and showed that large-scale LAAs can still
cause islanding and generator disconnections. Subsequent work has investigated
the impact of spatial (i.e., location of the attack within a power grid
network structure) and temporal (i.e., the time of the attack at different
load conditions) aspects of the attack [7, 8, 9]. A framework to determine the
locations of the grid (i.e., victim buses) from which an attacker can launch
the most impactful attacks was formulated in [7] using the theory of second-
order dynamical systems. Furthermore, reference [9] showed that low inertia
conditions due to high penetration of renewable energy resources can
exacerbate the consequences of LAAs.

Recent work has also focused on defending against LAAs. These measures can be
classified into (i) offline strategies, such as the optimal deployment of
protection-enabled loads [6] or finding generator operating points that can
prevent the destabilizing effects of LAAs [10], and (ii) online measures such
as real-time detection and localization of LAAs by monitoring the load data
[11] or frequency/phase angle fluctuations from phasor measurement units
following an LAA [12, 13]. Reference [14] proposes an online strategy to
mitigate the destabilizing effects of LAAs by redispatching fast-acting
inverter-based resources.

Despite the growing work on LAAs, existing literature lacks a framework to
understand the full extent of the consequences of LAAs in realistic power grid
conditions and a formal comparison of different attack strategies. In this
work, we address this critical research gap. Firstly, we identify specific
LAAs that can lead to power grid failures, even when the system is scheduled
for $N-1$ contingencies. Secondly, we extensively compare the impact of two
types of LAAs (S-LAAs and D-LAAs) to better understand their relative
importance in power grid stability. We elaborate on the importance and
challenge of these two aspects in the following.

First, $N-1$ scheduling policies provide resilience to power grids against
various contingencies, including LAAs. To the best of our knowledge, the
majority of the work in this area, with the exception of [5] and our precursor
work [1], has ignored $N-1$ scheduling. However, [5] considers only a few
specific attack scenarios (e.g., “IoT demand attack even distributed across
all load points”) and only a few power grid operational conditions (e.g.,
system load). Again to the best of our knowledge, there is no work that has
considered different distributions of the attack across the victim nodes and
grid conditions. However, discovering risk-posing LAAs is challenging since
$N-1$ scheduling implies that component disconnections due to LAAs are
extremely rare. To overcome this challenge, we employ a specialized sampling
approach known as the _skipping sampler_ (see [15]), which has been used in
the literature to obtain samples of high-impact events of low likelihood in
low-inertial power systems [16]. The application in our case is significantly
more challenging, however, since simulating D-LAAs requires sampling multiple
load changes at each node in the network, which can result in a challenging
high-dimension sampling problem. To address this, we extend the skipping
sampler algorithm to utilize multiple skipping trajectories per proposal
instead of the single trajectory used in the works cited. The points sampled
on each trajectory can be understood to represent the nodal load changes
applied to the network. Thus, D-LAAs comprised of $N$ load changes over time
can be generated from $N$ skipping trajectories.

Second, as noted before, S-LAAs and D-LAAs are qualitatively different attack
strategies. However, it remains uncertain which of the two strategies has a
greater impact. To address this crucial question, we leverage the skipping
sampler framework introduced earlier. To this end, we introduce a realistic
and generalized model of D-LAAs. In particular, the original D-LAA proposed in
[6] involves (i) continuous change in the system load with infinitesimally
small intervals between successive attacks, and (ii) follows a pre-defined
pattern, namely, the attack increases the load when the system frequency drops
below the setpoint and decreases the load when the system frequency is above
the set-point (thus worsening the generation-load imbalance at all times).
Both these assumptions have limitations. Continuous load changes are
unrealistic considering the attacker’s limited access and load response times.
Moreover, the predefined attack pattern fails to capture other sophisticated
strategies that attackers might employ to destabilize the system. We
generalize both these aspects by (i) considering a discretized attack model
with a periodic and non-zero interval between successive attacks, and (ii)
independently sampling load attacks across time and space. By varying the
interval between the attack, we can transition from dynamic (short intervals
between attacks) and static attacks (long interval between attacks).

We address these two research questions by performing extensive simulations
using the skipping sampler. We simulate the power grid’s transient dynamics
using a third-order model that accounts for both the frequency as well as the
voltage dynamics [17]. We treat the cumulative magnitude of LAA load changes
that perturb the dynamical system as random variables whose realisations are
sampled independently at each target node according to a log-normal
distribution111Note that the proposed framework works independently of the
chosen underlying cyber-attack distribution., which was assumed to reflect a
realistic magnitude occurrence distribution based on the empirical data from
data breaches [18]. We perform simulations including appropriate emergency
responses (ERs), such as generation/load shedding, etc., for different system
load conditions (evening, night, etc.), and various degrees of network
vulnerability to load-altering attacks.

Our results provide several novel insights into different strategies that an
adversary can potentially adopt to cause a grave impact on the system by
manipulating large-scale demand. The results show that adversaries employing
D-LAA strategies can induce the disconnection of critical network components
when manipulating a smaller magnitude of loads in the network as compared to
S-LAAs. Furthermore, while S-LAA strategies cause the most damage to the power
grid during low-demand periods, the system is particularly vulnerable to
D-LAAs during the peak demand period. We also identify the dominant failure
modes in these different load conditions (in terms of emergency response that
is activated most frequently).

This paper significantly extends our preliminary work [1], providing an in-
depth extended analysis. Specifically, the contributions of this paper can be
summarized as follows:

  * •

We propose a rare-event sampling approach to uncover LAAs that lead to power
system emergency actions in $N-1$ secure power grids. We compare its
performance with several other sampling methodologies and show that the
skipping sampler identifies significantly more attack strategies (in terms of
load attack distributions over the victim nodes and attack time window) that
can lead to network component disconnections.

  * •

Our previous analysis [1] only considered S-LAAs, but in this paper, we expand
our investigation to include D-LAAs. Unlike S-LAAs, these dynamic attacks
enable adversaries to exploit both spatial and temporal distributions of load
changes to elicit network disruptions. However, analyzing the risks posed by
D-LAAs comes with increased sampling complexity. To address this challenge, we
extend the functionality of the skipping sampler algorithm to include
multiple, simultaneous skipping trajectories.

  * •

We exhaustively compare the impacts of S-LAAs and D-LAAs in terms of (i) the
proportion of vulnerable loads on the network, (ii) the average magnitude of
each load change prior to the first disconnection event, (iii) the interval
between attacks, and (iv) the diurnal load balance.

The rest of the paper is organised as follows. Section II introduces the
system model; Section III presents the statistical model and details of the
rare-event sampling method. Section IV discusses the simulation results and
Section V concludes.

##  II System model

The power system model for rare-event sampling is based on the 3rd-order model
for generators [19]. In addition, our mathematical model accounts for
generator governor action, automatic voltage regulation, and protection system
emergency responses which disconnect loads, lines, and generators when network
parameters exceed pre-defined thresholds.

###  II-A Power System Model

We consider a third-order model for generator dynamics, which includes both
frequency and voltage dynamics [17, 16], enabling reliable estimation of
frequency dynamics for longer transient periods. We model subsequent power
grid emergency actions occurring after the initial LAA, allowing us to capture
critical actions during both primary and secondary frequency control over the
entire simulation duration.

Consider a graph theoretic formulation of a power system
$\mathcal{G}=(\mathcal{N},\mathcal{T})$, where $\mathcal{N}$ is the set of
network buses and $\mathcal{T}$ is the set of transmission lines connecting
the buses. The set $\mathcal{N}$ has cardinality $|\mathcal{N}|=N+L$, and
contains $N$ generation buses and $L$ load buses. At each generation bus
$i=1,\dots N$, the dynamics for the voltage phase angle $\delta_{i}$, voltage
magnitude $E_{i}$, and governor action $\rho_{i}$ are described by:

$\displaystyle
M(\psi)\ddot{\delta}_{i}+D\dot{\delta}_{i}=\psi_{i}P^{G}_{i}-P^{L}_{i}(\iota_{i})$
$\displaystyle\hskip 66.86414pt-
E_{i}\sum_{j=1}^{N+L}B_{ij}\Omega_{ij}E_{j}\sin(\delta_{ij})$ (1a)
$\displaystyle
S_{i}\dot{E}_{i}=\psi_{i}(E_{f,i}-\text{v}_{i})-E_{i}+X_{i}\sum_{j=1}^{N+L}B_{ij}\Omega_{ij}E_{j}\cos(\delta_{ij})$
(1b)
$\displaystyle\dot{\rho}_{i}=-A_{i}\dot{\delta}_{i}(1-1_{\mathcal{\mathcal{\mathcal{W}}}}[\dot{\delta}_{i}]).$
(1c)

In a similar manner, the dynamics for $\delta_{i}$ and $E_{i}$ at each load
bus $i=N+1,\dots,N+L$ are given by:

$\displaystyle
M(\psi)\ddot{\delta}_{i}+D\dot{\delta}_{i}=-P^{L}_{i}(\iota_{i})-E_{i}\sum_{j=1}^{N+L}B_{ij}\Omega_{ij}E_{j}\sin(\delta_{ij})$
(2a) $\displaystyle
S_{i}\dot{E}_{i}=\psi_{i}E_{f,i}-E_{i}+X_{i}\sum_{j=1}^{N+L}B_{ij}\Omega_{ij}E_{j}\cos(\delta_{ij}).$
(2b)

In equations (1) and (2), the indicator variables $\psi_{i}$, $\Omega_{ij}$
and $\iota_{i}$ represent the state of relays which control generator, line,
and load disconnections respectively (discussed in Section II-C). Table I
provides descriptions of model variables.

TABLE I: Variables used in (1) and (2).  Symbol | Meaning | Units
---|---|---
$A_{i}$ | Governor’s droop response | MW/rad
$B_{ij}$ | Susceptance matrix | p.u.
$P_{i}^{G}$ | Net generation at node $i$ | p.u.
$P_{i}^{L}(\iota)$ | Net loads at node $i$ | p.u.
$D$ | System damping | %
$\delta_{i}$ | Phase angle at node $i$ | p.u
$\delta_{ij}$ | $\delta_{i}-\delta_{j}$ | p.u.
$\dot{\delta}_{i}$ | Frequency at node $i$ | p.u
$\ddot{\delta}_{i}$ | Rate of change of frequency (RoCoF) | p.u.
$E_{i}$ | Voltage at node $i$ | p.u.
$E_{f,i}$ | Machine $i$ rotor field voltage | p.u.
$M(\psi)$ | System angular momentum | Ws2
$\Omega_{ij}$ | Line $ij$ disconnection indicator | -
$\psi_{i}$ | Generator shed indicator | -
$\iota_{i}$ | UFLS counter | -
$S_{i}$ | Machine $i$ transient time constant | s
$X_{i}$ | Machine $i$ equivalent reactance | ohms
$\mathcal{W}$ | Governor’s deadband frequency range | Hz

At nodes $i=[1,\dots N]$, net generation is given by
$P^{G}_{i}\coloneqq\min\\{P_{i}^{\text{max}},P^{e}_{i}+\rho_{i}\\}$, where
$P_{i}^{\text{max}}$ is the nominal maximum power output of generator $i$,
$P_{i}^{e}$ is the equilibrium power of the generator at $t=0^{-}$ before any
disturbance is applied to the network, and $\rho_{i}$ models the action of the
generator governor, with dynamics given in (1c). The variable vi accounts for
automatic voltage regulation (cf. [16] for further details).

###  II-B Modelling network loads

For each network, we model a fixed, maximum magnitude of loads $P_{i}^{TL}$ at
each node. We model nominal equilibrium loads at $t=0^{-}$ as a fraction of
total loads are modelled to be online and active $P_{i}^{eq}<P_{i}^{TL}$,
reflecting typical power balances observed during critical points of the
diurnal cycle. We introduce the network variable $\nu\in[0,1]$ to represent
the proportion of total loads susceptible to LAAs (i.e., connectivity-enabled
loads that are vulnerable to cyber attacks). Thus, we decompose the total
loads at each node into a vulnerable component $\nu P_{i}^{TL}$, comprised of
unsecured, IoT devices, and a secure component $(1-\nu)P_{i}^{TL}$, comprised
of non-IoT or protected devices.

Load-altering attack model: LAA refers to the manipulation of the vulnerable
component of network loads by an attacker to disrupt the power balance of the
network with the intention of leading the system to an unsafe state. Existing
literature models two different types of LAAs – S-LAAs consisting of a one-
time load manipulation) and D-LAAs consisting of multiple changes to loads in
the network over a period of time [4, 5, 6, 7].

In this work, we present a unified model for S- and D-LAAs consisting of
_discrete-time_ load changes. In particular, we model load manipulations
occurring at regular discrete times $t_{j}$ spaced $\mathcal{I}$ seconds
apart, i.e., $t_{j}=(j-1)\mathcal{I}$ for $j=1,2,\dots,n$, with $n$ load
changes during the simulation duration denoted by $T_{\max}$ seconds. We
assume that $\mathcal{I}\geq\mathcal{I}_{\min},$ where $\mathcal{I}_{\min}>0$
represents the physical limits on how quickly load magnitudes can change.
Under this model, an S-LAA can be considered a limiting case of a dynamic
attack with $n=1$. Without a loss of generality, the first load change is
modelled to occur at $t=t_{1}=0$ for both strategies.

We model the $j^{th}$ commanded nodal load change at time $t_{j}$ to be a
time-dependent fraction of the equilibrium loads, namely
$\eta_{i}(t_{j})P_{i}^{eq}$, for some factor
$\eta_{i}(t_{j})\in[\eta_{\min},\eta_{\max}]$, where $\eta_{i}(t_{j})<0$
($\eta_{i}(t_{j})>0$) can be understood as a command to reduce (respectively,
increase) the load of node $i$ at time $t_{j}$222Note that this model also
presents a generalized version of D-LAAs as compared to prior works [6, 7]
which modelled them as _continous_ load changes and the load change pattern
was restricted to being proportional to the instantaneous frequency
deviations.. The variables $\eta_{\min}$ and $\eta_{\max}$ represent the
minimum and maximum allowed values of $\eta_{i}(t_{j})$ in order to model the
physical limits of the load change. For a fixed number $n$ of load changes, a
load-altering attack is thus fully described by the $L\times n$ matrix $\eta$
whose entries describe all the commanded load changes, more specifically
$\eta_{i,j}:=\eta_{i}(t_{j})$. The commanded load change at time $t$ at node
$i$ can be modelled via the step-function
$\hat{u}_{i}:\mathbb{R}^{+}\mapsto\mathbb{R}^{+}$:

$\hat{u}_{i}(t)=P_{i}^{eq}\sum_{j=1}^{n}\eta_{i}(t_{j})\mathbbm{1}_{[t_{j},t_{j+1})}(t),$
(3)

where $\mathbbm{1}_{[t_{j},t_{j+1})}(t)=1$ if $t_{j}\leq t<t_{j+1}$ and $0$
otherwise. However, despite the commanded load change, the fixed quantity of
loads on the network limits the realized load change experienced. Therefore,
at time $t$ during the simulation, the actual net load at node $i$, including
the realized LAA but not the emergency responses, is given by:

$\hat{P}_{i}^{L}(t)\coloneqq(1-\nu)P_{i}^{eq}+\nu\min\big{[}0,\max\big{(}P_{i}^{TL},P_{i}^{eq}+\hat{u}_{i}(t)\big{)}\big{]}.$
(4)

The second term on the RHS represents the realized load change, which
restricts the authority of the attacker to manipulate the vulnerable loads
present at node $i$ between a minimum of $0$ MW and an upper limit of $\nu
P_{i}^{TL}$. This upper limit also highlights that the attacker can bring
inactive vulnerable loads online during the LAA. To record the realised load
change applied at $t=t_{j}$, at each node, we define the set of time-dependent
realized load changes
$\\{\lambda_{i}(t_{1}),\lambda_{i}(t_{2}),\dots,\lambda_{i}(t_{n})\\}$, where
each $\lambda_{i}(t_{j})\in\mathbb{R}$ at node $i$ can be calculated as:

$\lambda_{i}(t_{j})\coloneqq\hat{P}_{i}^{L}(t_{j})-\hat{P}_{i}^{L}(t_{j}-\delta
t),\quad j=1,\dots,n,$ (5)

where $\delta t$ represents a small time increment. At each node $i$, the
realized load changes $\lambda_{i}(t_{j})$’s can be aggregated over time to
describe the full extent of the nodal attack, yielding the nodal cumulative
realised attack $\Sigma_{i}$ defined as

$\Sigma_{i}\coloneqq\sum_{j=1}^{n}\left|\lambda_{i}(t_{j})\right|.$ (6)

This is a measure valid for both S- and D-LAAs of the total magnitude of
realised load changes deployed against each node (in units of power) over the
entire simulation.

###  II-C Emergency Responses

We model independent network protection systems intended to emulate the
systems designed to protect sensitive power system components from transients
in network frequency and voltage following a change in the active power
balance. We provide a brief description of these protection systems.

####  II-C1 Generation shedding

Two schemes are modelled to disconnect generators from the network: (1)
excessive RoCoF induced generation shedding (RIGS), where generators are
disconnected if the locally measured RoCoF $|\ddot{\delta}_{i}|$ exceeds a
pre-specified upper threshold, and (2) over frequency generation shedding
(OFGS), which disconnects generators when nodal frequency $\dot{\delta_{i}}$
deviations surpass a pre-specified upper threshold.

####  II-C2 Load shedding

Loads in the network can be shed via two independent schemes: (1) Under
Voltage Load Shedding (UVLS) and (2) Under Frequency Load Shedding (UFLS). We
model a simplified UVLS scheme that disconnects 5% of nodal loads if nodal
voltage magnitude $E_{i}$ falls below the pre-specified threshold of 0.9 p.u
for more than 5 seconds. This mechanism activates only once and is designed to
prevent voltage collapse. We introduce an indicator variable
$\mathcal{V}_{i}$, which equals $1$ when the nodal voltage magnitude at node
$i$ meets the activation criterion for the UVLS relay and $0$ otherwise.

The UFLS scheme is designed to prevent frequency collapse by progressively
disconnecting loads when nodal frequency $\dot{\delta}_{i}$ falls below a
strictly decreasing sequence of thresholds
$F^{U}\coloneqq\\{F^{U}_{1},\dots,F^{U}_{4}\\}$ where $F^{U}_{j-1}>F^{U}_{j}$
[16]. At each threshold, 10% of equilibrium loads $P^{eq}_{i}$ is disconnected
to arrest the decline in nodal frequency. Letting
$\mathcal{F}_{i}\in\\{0,1,2,3,4\\}$ count the cumulative number of UFLS relay
activations at node $i$ at each time step $t$ in the simulation, the net load
inclusive of any UFLS and UVLS activations up to time $t$, is a dynamic
variable in the power system model:

$P_{i}^{L}(t,\iota)=\iota_{i}\cdot\hat{P}_{i}^{L}(t),$ (7)

where
$\iota_{i}\coloneqq\big{(}1-0.1\mathcal{F}_{i}-0.05\mathcal{V}_{i}\big{)}.$
Correspondingly, load-shedding events reduce the effectiveness of LAAs by
disconnecting compromised loads, degrading the attacker’s ability to further
manipulate the network’s power balance.

####  II-C3 Line disconnection

We model the disconnection of inter-area transmission lines, as their loss can
exacerbate frequency instability and lead to islanding. When the power flow
$\phi_{ij}\coloneqq B_{ij}E_{i}E_{j}\sin(\delta_{ij})$ through any
interconnector line $(i,j)$ exceeds a pre-set power threshold $P^{\phi}$, the
line $(i,j)$ is tripped, and power is no longer allowed to flow through it for
the rest of the simulation. The indicator $\Omega_{ij}$ represents the state
of the interconnector between nodes $i$ and $j$, with the value $1$ indicating
nominal operation. Once the conditions for line disconnection are triggered,
$\Omega_{ij}$ switches to $0$ for the rest of the simulation and the power
flow on this line stops [16].

As the simulation evolves, the ER model inspects $\dot{\delta}_{i}(t)$,
$\ddot{\delta}_{i}(t)$, $E_{i}(t)$ and $\phi_{ij}(t)$ from the power system
model (1). Once the activation criterion for an ER is observed, the
corresponding ER is activated. This results in a discontinuity in (1), where
changes to the relevant line, generator, or load are applied. The simulation
is then resumed with the new network parameters.

##  III Statistical model and rare-event sampling approach for LAAs

Direct prediction of _network-threatening_ LAA characteristics is particularly
challenging, as they can exploit unforeseen vulnerabilities in network
architecture or operational practices, inducing emergency responses in
unforeseen ways. Instead, this work employs a sampling approach to identify
LAAs that lead to at least one subsequent emergency response, by applying
randomly generated load changes and network conditions (referred to as
‘proposals’) to the power system described in Section II-A. We then construct
the sample of successful LAAs by only retaining proposals that (1) trigger at
least one emergency response and (2) follow a desired distribution.

###  III-A The statistical model

To enable the aforementioned sampling approach, we account for the uncertainty
in key LAA characteristics (magnitude and frequencies) and network parameters
(load vulnerability and equilibrium conditions) as follows:

  * •

Distribution of network variables: To explore the interplay between LAAs and
different network operational configurations, we sample two network
parameters: (i) the network vulnerability ratio $\nu$ from a uniform
distribution $\nu\sim\mathcal{U}([0,1])$ and (ii) the diurnal power balance
scenarios $\tau$ uniformly from the set $\\{1,2,3,4\\}$. These scenarios
correspond to expected load magnitudes observed during night, morning,
afternoon, and evening hours, respectively. The nodal equilibrium loads
$P_{i}^{eq}$ under these load scenarios are proportional adjustments to the
published power balance of the networks investigated [20, 19]. The relative
factors $\\{0.4,1,0.85,1.3\\}$ for $\tau=1,\dots,4$ reflect the peaks in the
UK’s daily load cycle [21].

  * •

Distribution of LAA frequency and magnitudes: For dynamic LAAs, we assume the
time between attack $\mathcal{I}$ is a uniform random variable with
$\mathcal{I}\sim\mathcal{U}([1,T_{\max}])$, excluding intervals of less than
$1$ second to account for the time taken to process and execute programmed
load changes (i.e., $\mathcal{I}_{\min}=1$ second), and $T_{\max}=60$ seconds
is the total duration of each simulation. Once $\mathcal{I}$ has been sampled,
we determine the corresponding number $n$ of attacks with $n=\lfloor
T_{\max}/\mathcal{I}\rfloor$ and assume they are equally interspersed at times
$t_{j}=(j-1)\mathcal{I}$ for $j=1,\dots,n$.

  * •

Distribution of LAA magnitudes: We generate a set of commanded proportional
load changes by sampling uniformly at random from
$[\eta_{\min},\eta_{\max}]^{L\times n}$ and storing them in a $L\times n$
matrix $\eta$. Each entry $\eta_{i,j}$ represents the load change factor at
time $t_{j}$ with respect to the equilibrium load $i$. In particular, we set
$\eta_{\min}=-1$ and $\eta_{\max}=5$. The minimum value of $-1$ signifies
reducing the entire vulnerable part of the load at node $i$ to zero, while
positive values indicate load increases. The maximum load change factor is set
to $5$ to account for the fixed total load magnitude in the network.

Therefore, in our setting, a load-altering attack is fully characterized by
the quadruple $(\eta,\nu,\tau,\mathcal{I})$, where $\nu$ is the proportion of
vulnerable loads, $\tau$ describes the network power equilibrium state,
$\mathcal{I}$ is the interval between attacks and the $L\times n$ matrix
$\eta$ whose entries describe all the commanded load changes (recall that $n$
is automatically determined by $\mathcal{I}$). Using the power system model
(PSM) (1)–(2), for each attack $(\eta,\nu,\tau,\mathcal{I})$, we calculate the
nodal cumulative realized attack $\mathbf{\Sigma}$ and determine if it was
successful $S$ by

$(S,\mathbf{\Sigma})=\textbf{PSM}(\eta,\nu,\tau,\mathcal{I}),$

where $S=1$ if the attack was successful and $0$ otherwise.

Ultimately, we desire the nodal LAAs in the final sample reflect a realistic
distribution of LAA sizes expected in real-world attacks, where LAAs with
larger cumulative magnitudes $\mathbf{\Sigma}$ are rare. To achieve this bias
for each nodal LAA, we set a target probability density $\rho$ on
$\mathbb{R}^{L}_{+}$ to be
$\rho\sim\prod_{i=1}^{L}\text{LogNormal}(\mu,\sigma^{2})$, with parameters
$\mu=1$ and $\sigma=5$. This joint distribution is supported on
$[0,+\infty)^{L}$, and renders large magnitude LAAs at each node rare.

Thus, letting $A$ denote the set (or ‘event’) of LAA magnitudes that results
in the activation of at least one emergency response, our goal is to sample
from the conditional distribution
$\pi=\frac{\rho(\mathbf{\Sigma})\mathbbm{1}_{A}}{\rho(A)},$ where
$\rho(\mathbf{\Sigma})$ is the probability density of an LAA with nodal
magnitudes $\mathbf{\Sigma},$ (where
$\mathbf{\Sigma}=(\Sigma_{i})_{i=1}^{L}$), and $\rho(A)$ is the likelihood
that an LAA triggers an emergency response. The indicator
$\mathbbm{1}_{A}=S\in\\{0,1\\}$ takes a value of 1 if the LAA with nodal
magnitudes $\Sigma$ is associated with the activation of at least one
emergency response, and 0 otherwise. Together, $\pi$ describes a probability
distribution conditioned on the activation of at least one emergency response
occurring (event $A$).

###  III-B The Need for Rare-Event Sampling

TABLE II: Comparison of sampling methodologies when drawing realizations of network threatening nodal static-LAAs, against two versions of the IEEE 39 test network – one with protection systems and network characteristics tuned to satisfy the $N-1$ contingency criterion and a second one that is $N-0$ secure that does not consider contingencies in its design. Additionally, we investigate two uncertainty models for LAAs, the first, a uniform distribution (all LAAs are equally likely), and a LogNormal distribution (large magnitude LAAs are rare). All simulations were conducted for 100 hours, at which point the resulting sample size was recorded. Network Security | Sampling Methodology | Underlying Distribution | Accepted samples | Acceptance Rate, % | range $\nu$
---|---|---|---|---|---
$N-0$ Secure | Monte Carlo | $\mathrm{LogNormal}(1,5)$ | 0 | 0 | n/a
Random walk Metropolis | $\mathrm{LogNormal}(1,5)$ | 5000 | 5 | $[0.15,0.99]$
Skipping Sampler | $\mathrm{LogNormal}(1,5)$ | 57000 | 50 | $[0.12,0.99]$
$N-1$ Secure | Monte Carlo | $\mathcal{U}[0,5]$ | 0 | 0 | n/a
Monte Carlo | $\mathrm{LogNormal}(1,5)$ | 0 | 0 | n/a
Random walk Metropolis | $\mathcal{U}[0,5]$ | 200 | $0.0019$ | $[0.95,0.99]$
Random walk Metropolis | $\mathrm{LogNormal}(1,5)$ | 0 | 0 | n/a
Skipping Sampler | $\mathcal{U}[0,5]$ | 37000 | 37 | $[0.64,0.99]$
Skipping Sampler | $\mathrm{LogNormal}(1,5)$ | 9000 | 16 | $[0.66,0.99]$

Directly sampling $\pi$ using popular algorithms like Monte Carlo samplers can
be challenging due to several reasons: (1) Since networks are designed and
operated to maximize resilience against contingencies, the event $A$ of
interest is a rare event with a low probability density in realistic LAA
distributions. Consequently, evaluating a large number of proposed attacks is
necessary to identify a single relevant attack. (2) Understanding the
characteristics of network-threatening LAAs requires sampling individual load
changes at multiple nodes. For large networks, this leads to a high-
dimensional sampling problem that can hinder algorithm performance. This
problem is exacerbated for D-LAAs, which require sampling multiple load
changes at each node. (3) Emergency responses depend on complex interactions
between network conditions and the magnitude, spatial distribution, and timing
of load changes. As a result, the set $A$ may exhibit complex geometry and
possibly disconnected regions in the high-dimensional sample space of all
possible LAAs (cf. Figure 1).

The results summarised in Table II highlight the sampling challenge the
statistical model presents. Consider the case of sampling network-threatening,
static LAAs (S-LAAs) against the IEEE 39 test network.A Monte Carlo sampler
was unable to draw any samples of network-threatening LAAs in the stipulated
time, regardless of network resilience or LAA uncertainty model. This is a
consequence of the extreme rareness of such events in the high dimension
sample space of the problem, leading the algorithm to spend the entire
computational effort evaluating the more common, but irrelevant LAAs which do
not threaten the network. This renders these approaches inefficient or
ineffective when attempting to sample LAAs from $\pi$.

Markov Chain Monte Carlo (MCMC) sampling techniques can overcome many of these
challenges and enable efficient sampling of rare events. In general, MCMC
algorithms can be summarised as a two-step procedure: (1) a Proposal step,
where a new, proposed state $Z$ for consideration is generated from a
distribution of the user’s choice (with few restrictions); (2) an
Acceptance/rejection step\- using the target distribution $\pi$, the density
of the proposed state $\pi(Z)$ is compared to the density of the most recently
accepted state, denoted $Y$ (thus, $\pi(Y)$), always accepting $Z$ into the
sample if its density is higher (i.e. it is more likely to occur than an
already accepted state), and sometimes accepting $Z$ if its density is lower.

The generality of the MCMC methodology has led to multiple variants, including
the well-known Random-Walk Metropolis (RWM) algorithm. However, standard
formulations of RWM may also be ineffective when the event of interest (in our
case, event A) is both rare with respect to some underlying distribution and
disconnected in the sample space. In these cases, RWM algorithms are liable to
become ‘stuck’ in one component of the event, unable to escape to draw samples
from others [15]. We observe evidence of this phenomenon when an RWM is
applied to draw samples of network-threatening S-LAAs against different case
studies involving the IEEE 39 network, cf. Table II. The RWM algorithm
generates valid S-LAA samples only in simplified sampling cases, where the set
of threatening LAAs is much less rare in the state space, namely in the cases
of: (1) an insecure $N-0$ IEEE 39 network with Log-Normal distribution for
attack factors or (2) a $N-1$ secure IEEE 39 network with uniform distribution
for attack factors. However, the limited range of $\nu$ in the sample suggests
the algorithm became trapped in a local subset of $A$, and did not fully
explore the state space. When applied to the desired, more challenging case of
sampling threatening S-LAAs with a Log-Normal distribution against a ${N-1}$
secure network, the RWM also failed to generate any samples. Consequentially,
these results justify the need for specialized algorithms to efficiently
sample realistically distributed LAAs which may induce disconnections in
large, ${N-1}$ secure power systems.

###  III-C Sampling Procedure

The flexibility in the specification of the proposal step of MCMC algorithms
allows users to devise variant methodologies to tackle unique sampling
challenges. In this work, we apply one such variant, i.e., the skipping
sampler that ‘skips’ over $A^{c}$ until a point from $A$ is sampled or a
specific halting index $K$ is reached (in our simulations, we use $K=5$ for
each skipping path). We further amend this MCMC variant to provide a unified
framework to generate both S-LAAs ($n=1$) and D-LAAs ($n>1$) load changes
leading to the methodology detailed in Algorithm 1. In the special case of
S-LAA, it is enough to consider the interval between attacks constant and
equal to $T_{\max}$, which results in $n=1$ and a much smaller dimensional
sampling problem. The main ideas underlying the procedure presented in
Algorithm 1 can be summarised as follows:

1 Input: halting index $K$, number of steps $m$;

2 Sample a network state triplet as
$(\nu,\tau,\mathcal{I})\sim\mathcal{U}[0,1]\times\mathcal{U}[0,4]\times\mathcal{U}[1,T_{\max}]$;

3 Calculate the number of load attacks $n=\lfloor T_{max}/\mathcal{I}\rfloor$;

4 Generate a $L\times n$ matrix $\eta^{(0)}$ with i.i.d. entries sampled from
$\mathcal{U}[-1,5]$ describing an initial set of commanded load changes;

5 $(S_{0},\mathbf{\Sigma}_{0})=\textbf{PSM}(\eta^{(0)},\nu,\tau,\mathcal{I})$;

6 $Y_{0}\coloneqq[S_{0},\mathbf{\Sigma}_{0},\nu,\tau,\mathcal{I}]$;

7 for  _$w=0,1,\dots,m-1$_ do

8 Set $k=1$ and $\hat{\eta}^{(1)}\coloneqq\eta^{(t)}$;

9 Generate $n$ skipping directions $\mathbf{\Phi}_{j}\in\mathbb{R}^{L}$ for
$j=1,\dots,n$;

10while  _$S_{k}\neq 1\textbf{{ and }}k <K$_ do

11 Generate $n$ random distance increments $R_{j}^{(k)}$ conditionally on
$\mathbf{\Phi}_{j}$ for $j=1,\dots,n$;

12 Calculate the new matrix of commanded load changes by setting
$\hat{\eta}^{(k+1)}_{i,j}\coloneqq\hat{\eta}^{(k)}_{i,j}+R_{j}^{(k)}\mathbf{\Phi}_{j}$
for $j=1,\dots,n$;

13
$(S_{k+1},\mathbf{\Sigma}_{k+1})=\textbf{PSM}(\hat{\eta}^{(k+1)},\nu,\tau,\mathcal{I});$

14 Increment $k$ by $1$;

15

16if  _$S_{k}=1$_ then

17 Set $Z\coloneqq[S_{k},\mathbf{\Sigma}_{k},\nu,\tau,\mathcal{I}]$;

18 Evaluate the acceptance probability:
$\alpha(Y_{t},Z)=\begin{cases}\min\left(1,\frac{\pi(Z)}{\pi(Y_{t})}\right)&\text{
if }\pi(Y_{t})\neq 0,\\\ 1,&\text{ otherwise, }\end{cases}$ Generate a uniform
r.v. $V$ on $(0,1)$;

19 if  _$V\leq\alpha(Y_{i},Z)$_ then

20 $Y_{w+1}=Z$;

21 else

22 $Y_{w+1}=Y_{w}$;

23 else

24 $Y_{w+1}=Y_{w}$;

25Output: $[Y_{1},Y_{2},\dots,Y_{m}]$

Algorithm 1 Adjusted skipping sampler algorithm

The proposal step: We first sample (i) the parameters governing the network,
i.e., vulnerability ratio ($\nu$) and the load scenario ($\tau$), and (ii) the
interval between attacks ($\mathcal{I}$), which to compute the number of load
changes as $n=\lfloor T_{max}/\mathcal{I}\rfloor$. We then generate $n$
distinct $L$-dimensional commanded load changes
$\eta^{(0)}_{1},\dots\eta^{(0)}_{n}$ uniformly at random from $[-1,5]^{L}$
(recall that $\eta^{(0)}_{j}$ encodes the attack at time $t_{j}$). These,
along with the network parameters, are evaluated in the power system model. If
no emergency response is activated, we return to and update each commanded
load change using the key feature of the skipping sampler- we linearly
displace and thus update each
$\eta^{(0)}_{j}\to\eta^{(1)}_{j}=\eta^{(0)}_{j}+\Phi_{j}r_{j}$ by generating
$n$ trajectories $\Phi_{1},\dots,\Phi_{n}$ from a unit Gaussian density as
well as $n$ random distance increments $r_{1},\dots,r_{j}$ from a suitable
distribution with strictly positive support. The updated load changes
$\eta^{(1)}_{1},\dots\eta^{(1)}_{n}$, along with the previous network
conditions, are applied to the power system as before. This linear update for
$\eta^{(k)}_{j}$ is repeated by sampling further random distance increments
until either an emergency response is reported by the power system model, or
the halting index $K$ is reached. Taking the final outcome of this subroutine,
the resulting cumulative LAA magnitudes $\Sigma_{k}$, emergency response
indicator $S_{k}$ as well as the network parameters ($\nu,\tau,\mathcal{I}$)
are reported as the proposal $Z$ to be accepted or rejected in the next step.

The acceptance/rejection step: The proposal
$Z\coloneqq(\Sigma_{k},\nu,\tau,\mathcal{I})$ is evaluated using the target
density $\pi$ to determine whether it should be admitted to the final sample.
As described earlier for all MCMC algorithms, if the density of $Z$ ($\pi(Z)$)
is greater than the density of the most recently accepted state $\pi(Y_{w})$,
then $Z$ is added to the final sample and can be understood as the algorithm
moving to a more likely state in the density $\pi$. Conversely, if its density
is lower, it is sometimes accepted (as controlled by the uniform r.v. $V$).
This move to a less desirable state allows the algorithm to escape high-
density regions of the state space and promotes greater exploration.

This particular formulation of an MCMC algorithm allows more efficient
sampling of rare events which may be disconnected in state space, as the
skipping proposal mechanism prioritizes exploration of the state space. This
emphasis allows the algorithm to avoid being trapped in local modes of $\pi$
depending on the initial starting point of the algorithm, a potential weakness
of other MCMC algorithms (e.g. a random walk Metropolis algorithm). Further,
the use of multiple skipping paths for each load change over time offers both
a general and practical mathematical framework to represent both S-LAAs and
D-LAAs. This technique eliminates the need to change the dimensions of the
state space between proposals when more (or fewer) load changes are required.

As per Table II, in contrast to MC and RMW algorithms, the skipping approach
generated sufficiently large samples in all cases and was the only algorithm
to generate S-LAAs of interest in the most realistic scenario where a
resilient, $N-1$ compliant network was subjected to LAAs with rare large
attacks (LogNormal distribution), thus showing the effectiveness of the
methodology in sampling relevant LAAs.

Figure 1: Illustration of a skipping proposal for a D-LAA in the commanded
load change $\eta$ state space (where the ‘skipping’ update occurs) [left] and
the realised load change ($\Sigma$) state space (where proposals are
accepted/rejected) [right]. The set $A$ of realised load changes associated
with emergency responses is denoted in orange. Each sampled D-LAA is composed
of $n=3$ load changes in $L=2$ nodes and is thus described by three vectors
$\hat{\eta}_{1},\hat{\eta}_{2},\hat{\eta}_{3}\in\mathbb{R}^{2}$, which are
thus represented as distinct points [left]. Instead of simply discarding the
initial sampled unsuccessful LAA $\eta^{(1)}$ ($S_{1}=0$ and the corresponding
point $\mathbf{\Sigma}_{1}$ is not in the target region $A$ in orange), each
commanded load change (at $t_{1}$, $t_{2}$, and $t_{3}$) is updated in a
linear fashion by randomly generating new points along the random directions
$\Phi_{1}$, $\Phi_{2}$, and $\Phi_{3}$ respectively, determining
$\hat{\eta}^{(2)},\hat{\eta}^{(3)},\dots$ and thus of
$\mathbf{\Sigma}_{2},\mathbf{\Sigma}_{3},\dots$ (via the power system model,
denoted PSM). This skipping process continues until the resulting realised
load changes trigger at least one emergency response, corresponding to
$\mathbf{\Sigma}_{4}$ (which is the first in the target region $A$ in orange,
i.e., the first with $S_{4}=1$).

##  IV Simulations

###  IV-A Simulation Settings

We consider two well-known test networks in the literature: the Kundur 4-bus,
two-area system (KTAS) [19] and the IEEE 39 network [20]. The Kron-reduced
versions of these networks are used, with the reduced KTAS comprising $N=4$
nodes with generation units and $L=2$ pure loads and the reduced IEEE 39
network consisting of $N=10$ generator buses (two of which are modelled to
include loads as well), and $L=17$ pure load buses. In the KTAS network, the
interconnectors correspond to the transmission lines connecting Areas 1 and 2
(see [22]). The larger IEEE 39 network was divided into 3 areas, connected by
a total of 4 interconnector lines (see [20]). The emergency response
parameters are adjusted to ensure $N-1$ security for each system.

Simulations are conducted independently for each network and for each LAA type
(dynamic and static). The system is initialized in one of 4 power equilibrium
states specified by the sampled load scenario at $t=0^{-}$. The initial states
of each system333Note: Parameters for each network can be found at
https://github.com/maldongoodridge/power_system_model_MG.git variable
($\delta_{i}\left(0\right)$, $E_{i}\left(0\right)$, $\rho_{i}(0)$) are
determined numerically such that $\ddot{\delta}_{i}\approx 0$. The skipping
sampler is employed to generate a sample of LAA magnitudes and network
parameters, conditioned on at least one activation of an emergency response.
In each proposal step of the algorithm, we generate a matrix of commanded load
changes $\hat{\eta}^{(k)}\in\mathbb{R}^{L\times n}$. Along with the network
variables, these inputs are fed into the power system model (1), which we use
to simulate frequency and voltage dynamics for $T_{\max}=60$ seconds using
MATLAB. For each experiment, we conducted $50,000$ proposals.

The proposal densities of the sampling procedures are tuned separately for
IEEE 39 and KTAS simulations to achieve an acceptance rate (the proportion of
accepted proposals with respect to the total number of proposals) between 10%
and 40%, as recommended in the literature [23] for adequate exploration of the
sample space. This results in sample sizes of $\sim$14000 for D-LAAs, and
$\sim$9000 for S-LAAs, providing sufficient data for meaningful inferences.

###  IV-B Comparison of S- and D-LAAs via Statistical Analysis

We perform statistical analysis (of the impactful LAAs uncovered by the
skipping sampler in Section III-B) to compare S-LAAs and D-LAAs in terms of
the average magnitude of load change before the first disconnection event,
denoted
$\mu_{\lambda^{-}}:=\sum_{i=1}^{L}\int_{0}^{r}\lambda_{i}(t_{j})\mathbbm{1}_{\\{t=t_{j}\\}}dt$
where $t=r$ is the response time – the time of first disconnection for S-LAAs
and D-LAAs. Simulation results for the IEEE 39 network (Figure 2) reveal that
both the average and minimum of $\mu_{\lambda^{-}}$ are lower for D-LAAs
(Figure 2(a)) as compared to S-LAAs (Figure 2(c)), with a hypothesis test
comparing the averages revealing the difference, is statistically significant
($p$-value = $2\times 10^{-25}$). The reason behind this is that static
attacks must deliver the entire load change needed to destabilize the network
in a single event at time $t=0$, leading to a bias towards larger magnitude
load changes. Conversely, the ability to induce multiple load changes provides
additional mechanisms to trigger disconnections by varying the temporal
characteristics of the LAA. This reduces the reliance on magnitude-dominant
approaches in LAA design and enables the use of lower-magnitude D-LAA
strategies.

(a) Distribution of $\mu_{\lambda^{-}}$ for high-frequency D-LAAs
($\mathcal{I}<30$s).

(b) Distribution of $\mu_{\lambda^{-}}$ for low-frequency D-LAAs
($\mathcal{I}\geq 30$s).

(c) Distribution of $\mu_{\lambda^{-}}$ for S-LAAs .

Figure 2: Distributions of the average load change magnitude during response
time $\mu_{\lambda^{-}}$ for D-LAAs and S-LAAs applied against the IEEE 39
network, conditioned on the occurrence of at least one disconnection event
$X$. Illustrations are provided for high-frequency D-LAAs, low-frequency
D-LAAs, and static LAAs.

In the IEEE 39 network, we observe bi-modality of the conditional distribution
of $\mu_{\lambda^{-}}$ for both low-frequency D-LAA and S-LAA strategies.
These modes correspond to different attack strategies during varying load
balance scenarios ($\tau$). The lower mode is associated with LAAs leading to
emergency responses during the peak demand period (Evening), while the higher
mode is linked to LAAs during the nadir demand period (Night). Further details
explaining these observations are provided in Section IV-C. In contrast, the
distribution of $\mu_{\lambda^{-}}$ for high-frequency LAAs lacks bi-modality,
indicating the inherent threat posed by this attack strategy to network
operations. Such attacks induce disconnections in all scenarios with load
changes of much smaller average magnitude.

###  IV-C Impact of Intra-Day Power Equilibrium

Figure 3: Average cascade size for the IEEE 39-bus network.

We examined the vulnerability of both networks to LAAs under different load
balance scenarios, reflecting diurnal demand cycles. In the KTAS network
(Figure 3), the risk of disconnections is proportional to the magnitude of
equilibrium loads. This risk is maximum during periods of peak demand when
adversaries can induce a disconnection event by manipulating only $\sim 6\%$
of equilibrium loads (see Table III). This low threshold is due to the larger
absolute magnitude of loads online during peak demand and the network’s high-
stress state in this period. Thus, even marginal increases in loads can
overwhelm network generation and trigger a localized load-shedding event.
Conversely, the risk of LAA-induced disconnections is minimized during nadir
demand periods (referred to as Night). Although disconnections are rare during
this scenario, our analysis provides insights into how such events could
occur. In this scenario, successful attacks strategies are power-dominated
instead of period-dominated: observing $\mu_{\lambda^{-}}$ in Table III, both
S-LAA and D-LAA strategies require attackers to manipulate large proportions
of active loads to trigger disconnections. While an attacker can induce
disconnections using D-LAAs with smaller load change magnitudes, this
dramatically extends the response $r$, underscoring the dominance of low-
frequency, large-magnitude strategies in this scenario.

TABLE III: Summary of successful attacks in terms of the average magnitude of each load change ($\mu_{\lambda}$), average response time $r$, and average attack period $\mathcal{I}$ for each network and LAA strategy. Network | LAA | Metric | Night | Morning | Afternoon | Evening
---|---|---|---|---|---|---
IEEE 39 | D-LAA | Avg. $\mu_{\lambda}$ | 35.70 | 24.35 | 31.37 | 16.07
Avg. $r$, s | 34.75 | 48.82 | 39.99 | 50.79
Avg. $\mathcal{I}$, s | 32.25 | 31.45 | 30.46 | 33.96
S-LAA | Avg. $\mu_{\lambda}$ | 44.1 | 43.91 | 48.96 | 29.69
Avg. $r$, s | 29.60 | 48.62 | 40.73 | 50.17
KTAS | D -LAA | Avg. $\mu_{\lambda}$ | 55.29 | 18.86 | 32.62 | 6.31
$r$, s | 18.37 | 8.28 | 10.93 | 8.39
Avg. $\mathcal{I}$,s | 30.33 | 31.89 | 31.11 | 31.63
S-LAA | Avg.$\mu_{\lambda}$ | 72.36 | 20.17 | 36.78 | 6.56
$r$, s | 5.88 | 7.23 | 6.64 | 8.09

LAA-induced emergency responses exhibit fundamental differences in the larger
and more connected IEEE 39 network. Firstly, successful attack strategies in
this network require multiple load changes, with the average response time $r$
exceeding the average attack period $\mathcal{I}$ as shown in Table III. This
contrasts with the KTAS network, where an initial load change at $t=0$ was
typically sufficient to induce the first disconnection event. Secondly, the
results highlight how successful LAA strategies in the IEEE 39 network depend
on the load scenario. During low-demand scenarios, larger magnitude load
changes are required on average to trigger disconnections, resulting in a
significantly reduced response time. Conversely, during peak demand periods,
the preferred strategy is to manipulate a smaller proportion of loads to
trigger a disconnection, providing more response time.

These distinct characteristics of successful attack strategies based on the
load scenario are also evident in the induced emergency responses, see Figure
4. During low demand periods, RoCoF generation shedding dominated power
disconnected over the entire simulation, leading to the loss of a significant
proportion of generation in the network. This highlights the network’s
vulnerability to power-focused LAA strategies during the nadir demand period
(Night). In higher demand scenarios, the damping effect of loads, coupled with
the relatively smaller magnitude of load changes, limits RoCoF effects.
Instead, successful LAAs in this regime exploit the cumulative effect of
multiple load changes to gradually overwhelm already stressed generation
capabilities, resulting in UFLS dominating emergency responses. Generation
losses, a common failure mechanism profile of network overload, are rare in
this regime.

Overall, these findings suggest the IEEE 39 network is particularly vulnerable
to LAA-induced disconnections in periods of extreme demand. This is further
supported by the bi-modal distribution of disconnections conditioned on load
scenarios in Figure 3, with increased vulnerability in nadir (Night) and peak
(Evening) demand scenarios.

(a) Average cascade sizes during nadir demand Night scenario.

(b) Average cascade sizes during peak-demand Evening scenario.

Figure 4: Average cascade size for the IEEE 39-bus network.

###  IV-D Analysis of D-LAA Attack Patterns

In the literature, a common descriptor of D-LAAs which potentially threaten
network security is the so-called ‘reverse governor’ model, characterized by a
negative relationship between network frequency and the effective network load
change. Such load changes act in opposition to and therefore exacerbate
frequency deviations, eventually exceeding preset protection systems
thresholds for activation [6].

We use a simple linear model,
$\hat{\Sigma}=\hat{\beta}_{0}-\hat{\beta}_{1}\dot{\underline{\delta}}$, where
$\hat{\Sigma}\coloneqq\sum_{i=1}^{L}\lambda_{i}^{k}$ is the net load change
over the network, $\dot{\underline{\delta}}$ is the system frequency at the
time of each load change event, and $\beta_{1}$ is the slope coefficient
measuring the impact of load changes on frequency. The reverse governor model
posits that $\beta_{1}$ should be negative. We investigate this claim using
the samples generated from our simulations.

Figure 5: Summary of linear regression analyses between network frequency and
network load changes in each sample. The chart illustrates the proportion of
results where the relationship between load changes and frequency conforms to
the reverse governor model (blue), the proportion which does not (yellow), and
the proportion where there is no statistically significant relationship
between loads and frequency (red). Samples with fewer than 3 load changes were
excluded to remove trivial regression cases. We take a $p$-value less than 5%
as an indicator of statistical significance.

The pie-chart in Figure 5 suggests a reverse governor model was statistically
valid in only 40% of cases sampled where LAAs led to subsequent emergency
responses. In these cases, negative slope estimators ($\hat{\beta}<0$) with
p-values less than $0.05$ indicate there is less than a 5% chance this
relationship is spurious. However, in 52% of the sample, there was no
statistically significant relationship between network frequency and LAA load
changes. Surprisingly, in 8% of samples, there was even a positive
relationship between load changes and frequency deviations, directly
contradicting the assumption of a reverse governor model for D-LAAs.

These results reflect the literature in that when a meaningful relationship
between load changes and frequency can be extrapolated, the reverse governor
model dominates, describing $80\%$ of such cases. However, in general, the
reverse governor model is applicable only to a minority of all potential
failure scenarios. The statistical analysis highlights variations in network
frequency alone explained 22–30% of the variation in network loads, as
measured by the average coefficient of determination $R^{2}$. This suggests
the reverse governor model is actually only one of many potential dynamic
attack strategies.

##  V Conclusions and Future Research

In this work, we present a rare-event sampling framework to uncover impactful
LAAs in power grids designed with $N-1$ security criteria and provide a
rigorous comparison between S-LAAs and D-LAAs. Our results highlight the
importance of adopting specialized sampling methodologies such as the skipping
sampler in discovering impactful LAAs. Furthermore, we draw the following
conclusions: (i) D-LAAs can induce power grid failures by manipulating a
smaller fraction of IoT loads (cumulatively) as compared to S-LAAs. (ii) In
low-demand scenarios, larger magnitude load changes are required on average to
trigger disconnections (i.e., the dominant attack strategy being S-LAAs).
Conversely, during peak-demand periods, manipulation of a smaller and high-
frequency load manipulation triggers a disconnection (i.e., the dominant
attack strategy being D-LAAs). These findings are corroborated by the time
taken and the nature of emergency responses. (iii) We discovered a significant
number of attack samples that do not conform to the ‘reverse governor’ attack
model proposed in [6], highlighting the need for advanced sampling
methodologies to identify impactful LAAs. In the future, we will investigate
defence measures to enhance the resilience of power systems to LAAs using the
insight derived from this work.

## References

  * [1] M. P. Goodridge, S. Lakshminarayana, and C. Few, “Analysis of load-altering attacks against power grids: A rare-event sampling approach,” in _Proc. PMAPS_ , 2022, pp. 1–6.
  * [2] A.-H. Mohsenian-Rad and A. Leon-Garcia, “Distributed internet-based load altering attacks against smart power grids,” _IEEE Transactions on Smart Grid_ , vol. 2, no. 4, pp. 667–674, 2011.
  * [3] A. Dabrowski, J. Ullrich, and E. R. Weippl, “Grid shock: Coordinated load-changing attacks on power grids: The non-smart power grid is vulnerable to cyber attacks as well,” in _Proc. ACSAC_ , 2017.
  * [4] S. Soltan, P. Mittal, and H. V. Poor, “BlackIoT: IoT botnet of high wattage devices can disrupt the power grid,” in _Proc. USENIX Security Symposium_ , Aug. 2018, pp. 15–32.
  * [5] B. Huang, A. A. Cardenas, and R. Baldick, “Not everything is dark and gloomy: Power grid protections against IoT demand attacks,” in _Proc. USENIX Security Symposium_ , Aug. 2019, pp. 1115–1132.
  * [6] S. Amini, F. Pasqualetti, and H. Mohsenian-Rad, “Dynamic load altering attacks against power system stability: Attack models and protection schemes,” _IEEE Trans Smart Grid_ , vol. 9, no. 4, pp. 2862–2872, 2018.
  * [7] S. Lakshminarayana, S. Adhikari, and C. Maple, “Analysis of IoT-based load altering attacks against power grids using the theory of second-order dynamical systems,” _IEEE Trans Smart Grid_ , vol. 12, no. 5, pp. 4415–4425, 2021.
  * [8] J. Ospina, X. Liu, C. Konstantinou, and Y. Dvorkin, “On the feasibility of load-changing attacks in power systems during the covid-19 pandemic,” _IEEE Access_ , vol. 9, pp. 2545–2563, 2021.
  * [9] S. Lakshminarayana, J. Ospina, and C. Konstantinou, “Load-altering attacks against power grids under covid-19 low-inertia conditions,” _IEEE Open Access Journal of Power and Energy_ , vol. 9, pp. 226–240, 2022.
  * [10] S. Soltan, P. Mittal, and H. V. Poor, “Protecting the grid against mad attacks,” _IEEE Transactions on Network Science and Engineering_ , vol. 7, no. 3, pp. 1310–1326, 2020.
  * [11] E.-N. S. Youssef, F. Labeau, and M. Kassouf, “Detection of load-altering cyberattacks targeting peak shaving using residential electric water heaters,” _Energies_ , vol. 15, no. 20, 2022.
  * [12] S. Amini, F. Pasqualetti, M. Abbaszadeh, and H. Mohsenian-Rad, “Hierarchical location identification of destabilizing faults and attacks in power systems: A frequency-domain approach,” _IEEE Transactions on Smart Grid_ , vol. 10, no. 2, pp. 2036–2045, 2019.
  * [13] S. Lakshminarayana, S. Sthapit, H. Jahangir, C. Maple, and H. V. Poor, “Data-driven detection and identification of IoT-enabled load-altering attacks in power grids,” _IET Smart Grid_ , vol. 5, no. 3, pp. 203–218, 2022\.
  * [14] Z. Chu, S. Lakshminarayana, B. Chaudhuri, and F. Teng, “Mitigating load-altering attacks against power grids using cyber-resilient economic dispatch,” _IEEE Transactions on Smart Grid_ , pp. 1–1, 2022.
  * [15] J. Moriarty, J. Vogrinc, and A. Zocca, “A Metropolis-class sampler for targets with non-convex support,” _Statist. Comput._ , vol. 31, no. 6, 2021\.
  * [16] M. P. Goodridge, J. Moriarty, and A. Pizzoferrato, “A rare-event study of frequency regulation and contingency services from grid-scale batteries,” _Philos. Trans. Royal Soc. A_ , vol. 379, no. 2202, p. 20190433, 2021.
  * [17] K. Schmietendorf, J. Peinke, R. Friedrich, and O. Kamps, “Self-organized synchronization and voltage stability in networks of synchronous machines,” _The European Physical Journal Special Topics_ , vol. 223, no. 12, pp. 2577–2592, 2014.
  * [18] B. Edwards, S. Hofmeyr, and S. Forrest, “Hype and heavy tails: A closer look at data breaches,” _Journal of Cybersecurity_ , vol. 2, no. 1, pp. 3–14, 2016.
  * [19] P. Kundur and N. Balu, _Power System Stability and Control_. McGraw-Hill, 1994.
  * [20] T. M. Athay, R. Podmore, and S. Virmani, “A practical method for the direct analysis of transient stability,” _IEEE Transactions on Power Apparatus and Systems_ , vol. PAS-98, pp. 573–584, 1979.
  * [21] A. J. Pimm, T. T. Cockerill, and P. G. Taylor, “The potential for peak shaving on low voltage distribution networks using electricity storage,” _Journal of Energy Storage_ , vol. 16, pp. 231–242, 2018.
  * [22] M. P. Goodridge, S. Lakshminarayana, and C. Few, “Analysis of load-altering attacks against power grids: A rare-event sampling approach,” in _Proc. PMAPS_. IEEE, 2022, pp. 1–6.
  * [23] D. Gamerman and H. Lopes, _Markov Chain Monte Carlo, Stochastic Simulation for Bayesian Inference_ , 2nd ed. Chapman & Hall, 2006.

## Appendix: Case Studies

We present two case studies below, one which exemplifies the reverse governor
LAA model, and one which does not.

(a) Example of a high-frequency D-LAA which exhibits a ‘reverse governor’
relationship between frequency and load changes.

(b) Example of a low-frequency D-LAA where there is no significant
relationship between loads and frequency.

Figure 6: Network frequency profiles (black, with nodal frequencies super-
imposed in colour), network generation profiles (green), and network load
profiles (blue) during simulated D-LAA in the IEEE 39 network. Load changes
are depicted with a blue marker and emergency responses as vertical red dashed
lines at the time of their occurrence.

Figures 6(a) and 6(b) displays frequency, power, and load dynamics during
simulated D-LAAs against the IEEE 39 network. At $t=0^{-}$, the system is in
equilibrium, with loads equal to generation and network frequency of 60Hz.
Figure 6(a) illustrates a high-frequency D-LAA where load changes oscillate in
opposition to generation to exacerbate both frequency deviations and RoCoF,
eventually triggering a loss of load event and the disconnection of a
generator at 51 and 59 seconds respectively. This attack was modelled to occur
in the low-demand night scenario, thus generators had excess capacity to
attempt to stabilize frequency excursions. However, the adversary exploited
this response mechanism by commanding alternating load changes which
exacerbate the frequency dynamics instead of dampening them.

In Figure 6(b), the D-LAA is instead deployed against the IEEE 39 network
during a peak loading period when generation is near capacity. The attacker
deploys a low-frequency D-LAA strategy, where at $t=0$, the initial attack
increases loads beyond the capacity of the online generation. Thus, even after
ramping to full capacity, there is net demand in the network, leading to a
collapse in frequency. With no recourse, under frequency load shed relays
independently activate across the network to arrest the decline in frequency.
As the failure sequence is primarily induced by the initial LAA, there is
little relationship between load changes and frequency deviations, violating
the tenets of the reverse governor model.