Security and Statistics on Power Grids
Mauro Cesar Escobar Santoro

Submitted in partial fulfillment of therequirements for the degreeof Doctor of Philosophy
in the Graduate School of Arts and Sciences

COLUMBIA UNIVERSITY2019c 2019
Mauro Cesar Escobar Santoro
All Rights Reserved

ABSTRACT
Security and Statistics on Power Grids
Mauro Cesar Escobar Santoro
Improving the functioning and the safety of the electrical grids is a topic of great concern,
given its magnitude and importance in today’s world. In this thesis, we focus in these
two subjects.
In the first part, we study undetectable cyber-physical attacks on power grids, which
are attacks that involve physical disruptions, including tripping lines and load modifications, and sensor output alterations. We propose a sophisticated attack model describedunder the full Alternating Current (AC) power flow equations and show its feasibility onlarge grids from a test cases library. As counter-measures, we propose different defensive
strategies that the network’s controller can apply under a suspected cyber attack. These
are random, simple and fast procedures that change the voltages across the network andaim to unmask the current status of the system, assuming that the attacker cannot reactagainst their randomness.
Secondly, with access to data collected through Phasor Measurement Units (PMUs)
by a power utility in the United States, we perform statistical analyses on the frequency
and voltage time series that have been recorded at a rate of 30 Hz. We focus on intervalsof time where the sampled data shows to be in steady-state conditions and, with the useof appropriate signal processing filters, we are able to extract hidden anomalies such as
spatio-temporal correlations between sensors and harmonic distortions.

Table of Contents
List of Figures

iv

List of Tables

vii

Chapter 1 Introduction

1

Chapter 2 Preliminaries

4

2.1

Power Grids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.1.1

Transmission Lines . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.1.2

Generation and Demand . . . . . . . . . . . . . . . . . . . . . . .

9

2.1.3

Operational Considerations . . . . . . . . . . . . . . . . . . . . .

10

Optimal Power Flow Problems . . . . . . . . . . . . . . . . . . . . . . . .

11

2.2.1

AC-OPF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.2.2

DC-OPF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.2.3

Numerical Solutions . . . . . . . . . . . . . . . . . . . . . . . . .

15

2.3

Time-Varying Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.4

State Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

2.2

Chapter 3 Cyber-Physical Attacks

20

3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

3.2

Previous Work. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

3.3

Attack Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

3.3.1

30

Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

i

3.3.2Computational Viability . . . . . . . . . . . . . . . . . . . . . . .

33

3.3.3

Computational Implementation . . . . . . . . . . . . . . . . . . .

34

3.3.4

Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

3.4

The Follow-Up Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

3.5

Defense . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

3.5.1

Controlling Voltages Through Generation Changes . . . . . . . .

49

3.5.2

Overcoming Sensor Error, and the Current-Voltage Defense . . . .

54

3.5.3

Covariance Defense . . . . . . . . . . . . . . . . . . . . . . . . . .

59

Chapter 4 Learning from PMU data

64

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

4.2

Logic and Main Steps . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

4.3

Description and Averaging of the Time Series . . . . . . . . . . . . . . .

69

4.3.1

Moving Average and Covariance . . . . . . . . . . . . . . . . . . .

70

4.3.2

Averaging over Sliding Time Horizon . . . . . . . . . . . . . . . .

71

4.3.3

Quiet Periods . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

Fourier Filtering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

4.4.1

Cutting off high frequency components . . . . . . . . . . . . . . .

83

4.4.2

High Frequency Filter . . . . . . . . . . . . . . . . . . . . . . . .

96

4.4.3

Band-pass and Band-stop Filters . . . . . . . . . . . . . . . . . .

96

4.4

4.5

Covariance Matrix and Principal Component Analysis . . . . . . . . . . . 100
4.5.1

4.6

4.7

Singular Value Decomposition . . . . . . . . . . . . . . . . . . . . 102

Accounting for Temporal Correlations . . . . . . . . . . . . . . . . . . . . 106
4.6.1

Auto-Correlation Functions . . . . . . . . . . . . . . . . . . . . . 106

4.6.2

Cross-Correlation Residue (CCR) . . . . . . . . . . . . . . . . . . 111

Outcomes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

Chapter 5 Conclusion

116

Bibliography

131iiAppendices
Appendix

132
Correlation Plots

133

A.1 Auto-Correlation: Voltage Phase Angle . . . . . . . . . . . . . . . . . . . 133
A.2 Auto-Correlation: Voltage Magnitude . . . . . . . . . . . . . . . . . . . . 137
A.3 Cross-Correlation: Voltage Phase Angle . . . . . . . . . . . . . . . . . . . 141
A.4 Cross-Correlation: Voltage Magnitudeiii. . . . . . . . . . . . . . . . . . . 142

List of Figures
Figure 2.1

Branch Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

Figure 2.2

Turbine-Generator Model. . . . . . . . . . . . . . . . . . . . . .

17Figure 2.3

Frequency Contingency . . . . . . . . . . . . . . . . . . . . . . . .

18

Figure 3.1

Attacked zone A1 and its neighborhood . . . . . . . . . . . . . . .

36

Figure 3.2

Attacked zone A2 and its neighborhood . . . . . . . . . . . . . . .

38

Figure 3.3

Attacked zone A3 and its neighborhood . . . . . . . . . . . . . . .

41

Figure 4.1

Geographical location of PMU’s (anonymized coordinates), eachmark represents the position of a sensor. . . . . . . . . . . . . . . . . . .

68

Figure 4.2

Moving average of frequency for different α . . . . . . . . . . . . .

73

Figure 4.3

Moving average of voltage angle for different α . . . . . . . . . . .

74

Figure 4.4

Moving average of voltage magnitude for different α . . . . . . . .

75

Figure 4.5

Average over sliding time horizon of frequency for different S . . .

76

Figure 4.6

Average over sliding time horizon of voltage angle for different S .

77

Figure 4.7

Average over sliding time horizon of voltage magnitude for different S 78

Figure 4.8

Fourier transform of PMU signals . . . . . . . . . . . . . . . . . .

81

Figure 4.9

Fourier transform of complex voltage . . . . . . . . . . . . . . . .

82

Figure 4.10 Cut-off filter vector, λ = 4 Hz . . . . . . . . . . . . . . . . . . . .

84

Figure 4.11 Filtered time series for cut-off filter vector with λ = 4 Hz . . . . .

84

Figure 4.12 Cut-off filter vector, λ = 10 Hz . . . . . . . . . . . . . . . . . . . .

85

Figure 4.13 Filtered time series for cut-off filter vector with λ = 6 Hz . . . . .

85

Figure 4.14 Desired frequency response and its IIR . . . . . . . . . . . . . . .

87iv

Figure 4.15 Rectangular and Hann windows . . . . . . . . . . . . . . . . . . .

88

Figure 4.16 Low-pass filter vector, λ = 4 Hz . . . . . . . . . . . . . . . . . . .

89

Figure 4.17 Filtered time series for low-pass filter vector with λ = 4 Hz . . . .

89

Figure 4.18 Low-pass filter vector, λ = 1 Hz . . . . . . . . . . . . . . . . . . .

90

Figure 4.19 Filtered time series for low-pass filter vector with λ = 1 Hz . . . .

90

Figure 4.20 Low-pass filter vector, λ = 2 Hz . . . . . . . . . . . . . . . . . . .

91

Figure 4.21 Filtered time series for low-pass filter vector with λ = 2 Hz . . . .

91

Figure 4.22 Low-pass filter vector, λ = 3 Hz . . . . . . . . . . . . . . . . . . .

92

Figure 4.23 Filtered time series for low-pass filter vector with λ = 3 Hz . . . .

92

Figure 4.24 Low-pass filter vector, λ = 4 Hz . . . . . . . . . . . . . . . . . . .

93

Figure 4.25 Filtered time series for low-pass filter vector with λ = 4 Hz . . . .

93

Figure 4.26 Low-pass filter vector, λ = 5 Hz . . . . . . . . . . . . . . . . . . .

94

Figure 4.27 Filtered time series for low-pass filter vector with λ = 5 Hz . . . .

94

Figure 4.28 Low-pass filter vector, λ = 6 Hz . . . . . . . . . . . . . . . . . . .

95

Figure 4.29 Filtered time series for low-pass filter vector with λ = 6 Hz . . . .

95

Figure 4.30 High-pass filter vector, λ = 5 Hz . . . . . . . . . . . . . . . . . . .

97

Figure 4.31 Filtered time series for high-pass filter vector with λ = 5 Hz . . .

97

Figure 4.32 Band-pass filter vector, λ = 5 Hz . . . . . . . . . . . . . . . . . .

98

Figure 4.33 Filtered time series for band-pass filter vector with λ = 5 Hz . . .

98

Figure 4.34 Band-stop filter vector, λ = 5 Hz . . . . . . . . . . . . . . . . . .

99

Figure 4.35 Filtered time series for band-stop filter vector with λ = 5 Hz . . .

99

Figure 4.36 Eigen-values of frequency covariance matrix . . . . . . . . . . . . 103
Figure 4.37 Eigen-vector of frequency covariance matrix . . . . . . . . . . . . 103
Figure 4.38 Eigen-values of voltage angle covariance matrix . . . . . . . . . . 104
Figure 4.39 Eigen-vector of voltage angle covariance matrix . . . . . . . . . . 104
Figure 4.40 Eigen-values of voltage magnitude covariance matrix . . . . . . . 105
Figure 4.41 Eigen-vector of voltage magnitude covariance matrix . . . . . . . 105
Figure 4.42 Auto-correlation functions for frequency . . . . . . . . . . . . . . 107v

Figure 4.43 Auto-correlation functions for frequency using band-stop filter . . 108
Figure 4.44 Auto-correlation functions for frequency using band-pass filter . . 109
Figure 4.45 Auto-correlation residual maps for frequency . . . . . . . . . . . . 111
Figure 4.46 Cross-correlation residual matrices for frequency . . . . . . . . . . 112
Figure 4.47 Cross-correlation with sensor k = 134 for frequency . . . . . . . . 113
Figure A.1 Auto-correlation functions for voltage phase angle . . . . . . . . . 133
Figure A.2 Auto-correlation functions for voltage phase angle using band-stop
filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
Figure A.3 Auto-correlation functions for voltage phase angle using band-passfilter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
Figure A.4 Auto-correlation residual maps for voltage phase angle . . . . . . 136
Figure A.5 Auto-correlation functions for voltage magnitude . . . . . . . . . . 137
Figure A.6 Auto-correlation functions for voltage magnitude using band-stop
filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
Figure A.7 Auto-correlation functions for voltage magnitude using band-passfilter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
Figure A.8 Auto-correlation residual maps for voltage magnitude . . . . . . . 140
Figure A.9 Cross-correlation for voltage magnitude . . . . . . . . . . . . . . . 141
Figure A.10 Cross-correlation for voltage magnitude . . . . . . . . . . . . . . . 142

vi

List of Tables
Table 3.1

True and reported flow for attacked lines. (Overflow in bold) . . .

36

Table 3.2

Load and generation before and after the attack. . . . . . . . . . .

37

Table 3.3

True and reported flow for attacked lines. (Overflow in bold) . . .

39

Table 3.4

Load and generation before and after the attack. . . . . . . . . . .

40

Table 3.5

True and reported flow for attacked lines. (Overflow in bold) . . .

42

Table 3.6

Load and generation before and after the attack. . . . . . . . . . .

43

Table 3.7

AC Voltage Changes . . . . . . . . . . . . . . . . . . . . . . . . . .

54

Table 3.8

Current-voltage defense. . . . . . . . . . . . . . . . . . . . . . . . .

59

Table 4.1

Selected quiet periods. . . . . . . . . . . . . . . . . . . . . . . . . .

72

vii

Acknowledgments
I would like to thank my advisor, Prof. Daniel Bienstock, for being a great and dedicatedmentor and from whom I have learned a lot during my time at Columbia University.
Thanks to him I started to explore the world of optimization applied to the power grid,without previous knowledge of electricity. I feel very fortunate of receiving his devotion
to instructing me through the research that we have done.
I would also like to thank the wisdom, expertise and contribution that I have received
from faculty members of the department, specially from Profs. Garud Iyengar, Ali Hirsa,
Yuri Faenza, and Donald Goldfarb. I have truly enjoyed their classes and I appreciate
every comment and feedback that they have given me. I believe that, together with
Prof. Daniel Bienstock, they are a central piece in the optimization community of thedepartment.
I want to thank Prof. Michael Chertkov that hosted me for three weeks in Los Alamos
National Lab. Thanks to him, my time there was very productive and I could get a
new perspective of the analyses that we were doing. In addition, I am grateful forthe collaboration held with Dr. Jonatan Ostrometzky and Prof. Gil Zussman, from theElectrical Engineering department. The Fourier analysis that we did are based on theirexpertise and feedback.
IEOR has been a wonderful place to be part of, I will never forget the people that
I met in the department. Thanks to them, which many have become truly friends, theyears that I spent in New York City has been pleasant and enjoyable. Special mention
to Gonzalo, Camilo, Liz, Enrique, Julien, and Nouri.
I want to thank Eduard for the company through most of this process, I receiveviii

constant motivation and energy from admiring what he does, I am very grateful for what
we have. Lastly, to my family for being always supportive and encouraging, I constantly
feel them close and present even though they are far away. They are a fundamental pillarof who I am and the path that I have followed.

ix

CHAPTER 1. INTRODUCTION

1

Chapter 1
Introduction
The electrical grid is one the largest and most complex infrastructures built by mankind.
It has a great importance in our lives, moreover, we have become daily dependent on
it: most jobs make use of electricity, as well as modern systems of telecommunication,means of transportation, and important institutions —such as governments or hospitals.
In summary, the world’s economies rely on a solid system to supply power. Because ofthis reason, it is crucial to maintain a stable and safe system that provides energy to its
customers in a reliable manner.
The existing power grid in the United States nowadays presents several challenges. It
has exponentially grown across the last century without planning its future development
and, consequently, faces weaknesses when it has to satisfy the increasing demand seen
nowadays. As detailed in Gretchen Bakke’s book [10] from 2016, more than 70 percentof the transmission lines and transformers in the United States are twenty-five years
old and the average age of the power plants is thirty-four. As a result, the system haslarge inefficiencies that have caused an increasing number of power outages across thelast decades: 15 in 2001, 78 in 2007, 307 in 2011. Additionally, the unpredictability of
renewable energy sources has contributed to the instability of the grid —where production
has to meet demand at all times. The current grid is not prepared to deal with thisuncertainty.

CHAPTER 1. INTRODUCTION

2

These current challenges together with the complexity of the physics that supports
the generation and transmission of electricity throughout a system has motivated decadesof research on different topics related with the power network. The equations behind themodel of the power flows are not simple —they are, in its full and most realistic form,
non-linear and non-convex— and, therefore, difficult to understand and compute with.
Nevertheless, decades of investigation have lead us to have a trustworthy system thatworks well most of the time.

The main goal of this thesis is to contribute towards a better understanding of thepower grid’s functioning and improving its safety. For that purpose we first study cyberphysical attacks on the electrical grid, meaning that there is an adversary that has the potential ability to cause physical damage or disruption in the system —such as, a tripping aline, modifying demand in certain points, or changing the functioning of a transformer—
coordinated with a cyber attack. A cyber attack corresponds to a hack that blocks oralters the data that the grid’s controller obtains from sensors distributed at several locations in the network. These sensors provide information used to estimate the status
and the state variables of the system and, this system is used by the controller to makecontrol decisions. The purpose of the hack perpetrated by the adversary, is to hide fromthe controller his or her physical actions (the physical attack) in order that the controller
remains unaware of them and, potentially, executes inopportune decisions.
Diverse models of cyber and cyber-physical attacks have been proposed and studied
in the last years, most of them using a linear approximation of the power flow equations.
In this thesis, we analyze a new model where the adversary can modify demand of asubset of nodes of the grid and at the same time alter the measurements sampled in thatzone. We use the full description of the power flow equations and show that the attack
can cause important overloads that, if kept hidden from the controller for long enoughtime, would have unfavorable consequences.
We complement this section by also describing different strategies that the grid’s con-

CHAPTER 1. INTRODUCTION

3

troller (under the suspicion of an attack) can perform in order to unmask the real statusof the variables of the system. These strategies are simple generic procedures that canbe used in different kind of cyber attacks.

As a second theme of research, we center our attention on the measurements sampledover the grid. With new technologies and the digital age, modernization of equipment andcontrol of the network has improved, which is giving an opportunity to enhance processes
and have a better understanding of the network behavior. One particular example is the
upgrade of the sensors that measure electro-physical quantities across the grid (voltage,
frequency, current): from the old RTUs (remote terminal units) that produce samples
every few seconds to the new PMUs (phasor measurement units) that sample as fast as
30 times per second. With this new volume of information being generated, we might beable to get a new insight about the nature of the power flowing through the network.
We have at our disposal historical PMU measurements that cover a period of 15
months obtained from a US Independent System Operator; these data correspond tosampling collected in approximately 200 locations at a rate of 30 measurements per
second. With this information we performed statistical analyzes using covariance (correlation) matrices in order to understand several questions regarding the structure of thedata: Can we obtain useful information from the spetral decomposition of these matrices?
Or is the data temporarily correlated? We aim to answer these question by observing
graphical results of the statistical tools over selected periods of measurements.

This document is organized as follows: we begin by describing the basic properties
and models of the electrical grid and its power flow equations in Chapter 2, we continue
with the cyber-attack study in Chapter 3, and finalize with the statistical analysis of the
PMU sampling in Chapter 4.

CHAPTER 2. PRELIMINARIES

4

Chapter 2
Preliminaries
In this section we provide a description of a power system that we will use in the next
chapters. The mathematical model that is presented is widely used in the literature.

2.1

Power Grids

A power transmission system can be characterized by a graph, where the set of nodes
or buses represent the physical locations where the power is generated, consumed or
redistributed (see [13, 41]). We denote by N this set of buses and distinguish the set
G ⊂ N of generators; buses that do not generate power are called load buses. The set
E of edges of the graph, also called branches, represent the transmission lines of the
network. For bus k ∈ N , we denote by δ(k) ⊂ E the set of branches that are incident
with k. For a set A ⊂ N of buses, AC denotes the complement of A with respect to N ,that is N \A, and ∂A denotes the set of buses in A that are connected with some bus in
AC , i.e. the boundary of A. Notation: when x is a complex number, x will denote its
complex conjugate; and if X is a complex matrix, X H denotes its Hermitian transpose(i.e. conjugate transpose, X > ).
While complex power is generated and consumed throughout the system, transmission
lines are in charge of carrying the power across the different locations. In order toCHAPTER 2. PRELIMINARIES

5

understand the equations that describe a power flow system, we need to introduce threephysical concepts present in the model: voltage, current and power.
The amount of power flowing from bus k to bus m connected by the branch km
depends on the physical characteristics of the branch and the voltages Vk and Vm , respectively associated with both buses. Voltages at buses are the state variables of thesystem, they are the result of the whole network configuration (physical characteristics,
generation, demand). The voltage at bus k is represented as a complex number, in anyof the three different forms:
Vk = |Vk |ejθk = |Vk | θk = |Vk | cos(θk ) + j|Vk | sin(θk ) ,| {z }| {z }|{z}exponential

polar

(2.1)

rectangularwhere |Vk | is the magnitude of the voltage and θk is its (phase) angle. Voltage, also
called electric potential difference, makes sense when it is compared between two points
or with respect to a reference point, this difference is measured in volts or [V]. One volt
(as defined in [48]) “is the potential difference between two points of a conducting wire
carrying a constant current of one ampere, when the power dissipated between thesepoints is equal to one watt.”
The electric current —or simply current— is the flow of electric charge, in our caseelectrons, moving in a wire. The unit for measurement is the ampere or [A] definedas “constant current which, if maintained in two straight parallel conductors of infinite
length, of negligible circular cross-section, and placed one metre apart in vacuum, would
produce between these conductors a force equal to 2 · 10−7 [mks] unit of force (newton)
per metre of length” [48]. We denote by Ikm the complex current that flows from bus k
to bus m through branch km.
Finally, the electric power is the rate, per unit time, at which the electric energy istransferred through a transmission line. The unit of measurement is the watt or [W] that“is the power which in one second gives rise to energy of one joule” [48]. We denote by
Skm = Pkm + jQkm the complex power that flows from bus k to bus m, where the realpart Pkm is called active power and the imaginary part Qkm is known as reactive power .

CHAPTER 2. PRELIMINARIES

6

It is usually convenient to work with normalized units for the previously mentioned
quantities. For this purpose, we define a per-unit normalization as
quantity in per unit =actual quantity,
base value of quantity(2.2)

where the base value for quantities are picked to satisfy the same of relationship as theactual variables [13].
Before making explicit the power flow equations we need to model the branches ofthe system.

2.1.1

Transmission Lines

We now describe a model that characterizes the physical attributes of the transmission
lines (see [18, 101]), including the presence of transformers and shunt admittance. Thisis also known as π-model.
Consider a branch between buses k and m, see Figure 2.1. Safety operational func.
tioning indicates that the difference of voltage phase angles θkm = θk − θm between busesminmax
k and m should not larger that a pre-established bounds θkm
and θkm. In addition,

branch km has a series impedance
zkm = rkm + jxkm

(2.3)—that measures the opposition to the flow of current, where rkm is the resistance andxkm is the reactance— and series admittanceykm = (zkm )−1 = gkm + jbkm ,

(2.4)

2
2
where gkm = rkm /(rkm+ x2km ) is the conductance and bkm = −xkm /(rkm+ x2km ) is thesusceptance. There is also a shunt —such as a capacitor or inductor, that creates lowresistance— with shunt admittance
sh
sh
ykm
= gkm
+ jbshkm(2.5)

CHAPTER 2. PRELIMINARIES

7

and a transformer with ratio
Nkm = τkm e jσkm ,(2.6)

where τkm > 0 is its magnitude and σkm is the phase shift angle. The transformer islocated on the branch next to bus k, it scales the voltages by a factor of 1/Nkm andcurrents by a factor of Nkm .
∗
Nkm
Ikm

Ikm

Imk
s
ykm
=

Vk
Nkm

Vk

1
rkm + jxkm

y sh
2

y sh
2

Vm

Nkm = τkm ejσkm

Figure 2.1: Branch Model.

In the rest of this section we will omit the subscript km so as to simplify the equations.
With these elements, the branch admittance matrix is defined as
 



y sh  11y+−y −jσ ykk ykm2 τ2
τ e sh  .
 = 
Ykm = 1yymk ymm
y+
−y jσ
τe
2

(2.7)

Therefore, the complex current injections Ikm and Imk at buses k and m of the branch,respectively, can be expressed in terms of the voltages Vk and Vm and the admittancematrix Ykm :y sh  11y+Vk − y −jσ Vm IV
2 τ2 
τ esh  km  = Ykm  k  = 

.
y
1
Imk
Vm
−y jσ Vk + y +
Vm
τe
2








 

(2.8)

According to Omh’s law, the power injected at a branch is defined as the voltage ofCHAPTER 2. PRELIMINARIES
the bus times the complex conjugate of the injected current, in other words:
y sh  1
1
Skm = Vk Ikm = y +|Vk |2 − y s jσ Vk Vm ,
2
2 τ
τe

y sh 
1
Smk = Vm Imk = −y −jσ Vk Vm + y +|Vm |2 ,τe
28(2.9a)(2.9b)

from where we can deduce expressions for active and reactive power injected at bothends:
g sh  |Vk |2 |Vk |
g+
|V
|
g
cos(θ
−
σ)
+
b
sin(θ
−
σ)
, (2.10a)
−
m
km
km
2
τ2
τ


bsh  |Vk |2 |Vk |= − b+
+
|V
|
b
cos(θ
−
σ)
−
g
sin(θ
−
σ)
, (2.10b)
m
km
km
2
τ2
τ


|Vk |
g sh 
|Vm |2 −|Vm | g cos(θkm − σ) − b sin(θkm − σ) , (2.10c)
= g+
2
τ


bsh 
|Vk |
= − b+|Vm |2 +|Vm | b cos(θkm − σ) + g sin(θkm − σ) . (2.10d)
2
τ

Pkm =
Qkm
Pmk
Qmk



As we have seen, power can be specified as a function of the voltages at both ends ofa branch. Therefore, when convenient, we will denote it as Skm (|Vk |, |Vm |, θk , θm ) andSmk (|Vm |, |Vk |, θm , θk ) the complex power injected at each bus connected by branch km,
respectively.
Power losses can also be computed taking into account the active power. The active
power loss at branch km is defined as
Losskm = Pkm + Pmkg sh  |Vk |2
|Vk |
2
= g+
+
|V
|
−2
|Vm |g cos(θkm − σ).
m
2
2
τ
τ

(2.11)

Finally, we define the apparent power to be the magnitude (as complex number) of
the power injected by k and m to branch km:q2Pkm+ Q2km = ||(Pkm , Qkm )|||Skm | =q2Pmk+ Q2mk = ||(Pmk , Qmk )||,
and|Smk | =(2.12a)(2.12b)respectively. The apparent power has the same units as active or reactive power, however
it is common to use the unit voltampere or [VA]. Moreover, the branch has a capacitymax
Skmthat limits the amount of apparent power.

CHAPTER 2. PRELIMINARIES

2.1.2

9

Generation and DemandEvery generator in the system produces power that is injected into the network. We
denote by
Skg = Pkg + jQgk

(2.13)

the amount of complex power produced by generator k ∈ G, which must lie inside its capacities of generation. Depending on the physical capabilities of the generator, generated
active power has as lower and upper bounds Pkg,min and Pkg,max , respectively. Analogously,
generated reactive power has as lower and upper bounds Qg,min
and Qg,max, respectively.
kkElectrical power is also consumed across buses of the network. It is common to saythat the load of a bus is the amount of power that the bus demands, which is denotedby
Skd = Pkd + jQdk

(2.14)

for bus k ∈ N .
According to [7], power systems are split into balancing authority areas that are in
charge of maintaining the balance between load, generation, and interchange of power between the different balancing authority areas. This goal is accomplished through several
processes, including• Unit Commitment (UC): determines which generators will be on and off in anhourly based schedule, these generators are always ready to satisfy the forecasted
demand;• Economic Dispatch (ED): assigns the production level of the generators that wereselected by the UC procedure, minimizing the cost, subject to meeting the demandand physical constraints. The ED algorithm is run every hour using the day-aheadload prediction, but also every 5-10 minutes, using the minute-ahead forecast;• Primary Generation Control (PGC): performs on-line adjustments caused by generation outages, line tripping, demand fluctuations, or any spontaneous change inthe system;

CHAPTER 2. PRELIMINARIES

10• Automatic Generation Control (AGC): returns the frequency —the angular speed ofgenerator’s rotors (see Section 2.3 for further details)— to its nominal value of 60 Hz
after PGC actions; this is obtained by reassigning the generation output acrossdifferent balancing authority areas; AGC updates commands every 2-4 seconds.

2.1.3

Operational Considerations

As described in [13, p. 327], operational consideration must be taken into account withthe objective of understanding state and control variables of the system:• For generators, the active power Pkg and the voltage magnitude |Vk | can be specified—by varying turbine power and generator field current.
• Complex load Skd at every bus is known (by previous estimation), thus, is considered
as a parameter of the system.
We see then that active power can be set for all buses, however, this cannot be done
independently, the sum of all the Pkg ’s must equal the sum of active power loads andlosses. One generator is designated to have unspecified active power generation, and thisamount comes as the result of the system’s active power balance. This specific bus iscalled a slack bus or swing bus (for which we will reserve the index k = 0), that instead of
having specified P0g , both magnitude |V0 | and angle θ0 voltage values are indicated. Since
power flow equations (2.10) always depend on the difference between the angle phases of
two buses, by specifying the slack bus phase angle, we just set an angle reference for the
rest of the buses.
In summary, we have three types of buses, depending on the specified variables:
1. A voltage V source. At the slack bus.
2. P, |V | sources. At the other generators, also called voltage control buses.
3. P, Q sources. At the load buses.

CHAPTER 2. PRELIMINARIES

2.2

Optimal Power Flow Problems

2.2.1

AC-OPF

11

The Alternative Current Optimal Power Flow (AC-OPF) problem is formulated in order
to find an optimal solution of a power flow that minimizes the generated power subject
to physical laws of energy conservation. AC-OPF was introduced in 1962 as the problemof economic dispatch [26].
Summarizing the previous sections, we have as parameters of the system:• Vkmin and Vkmax , voltage magnitude bounds for every bus k ∈ N ;• Skd = Pkd + jQdk , the complex load at every bus k ∈ N ;• Pkg,min and Pkg,max , active power generation bounds for every generator k ∈ G;
• Qg,minand Qg,max, reactive power generation bounds for every generator k ∈ G;kk
kk ykm
• Ykm = [ yymk
ymm ], the branch admittance matrix, for every transmission line km ∈ E;

max
• Skm
, the apparent power capacity for every branch km ∈ E;maxmin
, bounds of the voltage angle difference for every branch km ∈ E;
and θkm• θkm• slack bus phase angle, usually specified as 0.
The control variables of the systems are:• |Vk |, voltage magnitude at every generator k ∈ G;
• θ0 , voltage phase angle of the slack bus;
• Pkg , active power generation at every non-slack generator k ∈ G\{0}.
And the state variables are:• |Vk |, voltage magnitude at every load bus k ∈ N \G;CHAPTER 2. PRELIMINARIES

12• θk , voltage phase angle at every non-slack bus k ∈ N \{0};• P0g , active power generation at the slack bus;
• Qgk , reactive power generation at every generator k ∈ G;
• Pkm and Pmk , active power injection at both ends of every branch km ∈ G;
• Qkm and Qmk , reactive power injection at both ends of every branch km ∈ G.
The AC-OPF problem is then formulated as follows:
minimizesubject tof ({Pkg }k∈G )(2.15a)X∀k ∈ G, (Pkg − Pkd ) + j(Qgk − Qdk ) =(Pkm + jQkm )(2.15b)km∈δ(k)Pkg,min ≤ Pkg ≤ Pkg,max ,Qg,min≤ Qgk ≤ Qg,max
kkX∀k ∈ N \G, −(Pkd + jQdk ) =(Pkm + jQkm )(2.15c)(2.15d)km∈δ(k)∀km ∈ E, Pkm + jQkm = ykk |Vk |2 + ykm |Vk ||Vm |ej(θk −θm )(2.15e)Pmk + jQmk = ymm |Vm |2 + ymk |Vk ||Vm |e−j(θk −θm )(2.15f)

2
max 2
Pkm
+ Q2km ≤ (Skm),

(2.15g)

2
max 2
Pmk
+ Q2mk ≤ (Skm
)

min
max
θkm
≤ θk − θm ≤ θkm(2.15h)∀k ∈ N , Vkmin ≤ |Vk | ≤ Vkmax(2.15i)θ0 = 0.

(2.15j)The objective function f depends on the active power generation at all generators, convex
functions —typically quadratic— are commonly used so as to minimize the amount of
power generated.
Constraints (2.15b) and (2.15d) describe the complex power balance at generatorsand load buses, respectively; generated power minus demand must equal the sum of
power injected into the branches incident to the bus. Whereas constrain (2.15c) states
the bounds for the power that a generator can produce.

CHAPTER 2. PRELIMINARIES

13

Constrains (2.15e)-(2.15f) summarize the power equations that buses inject at theirincident branches. These equations are explicitly stated in (2.10). Constraint (2.15g)
limits the amount of apparent power that can be injected at each branch, while constraint (2.15h) states the bounds for the difference of voltage phase angles between two
adjacent buses.
Finally, constraint (2.15i) established the bounds for the voltage magnitude and constraint (2.15j) sets the reference phase angle for the slack bus.
The AC-OPF problem formulated as in (2.15) is non-linear and non-convex; thesource of non-linearity is mainly provided by the power flow equations (2.15e)-(2.15f),
since they depends though trigonometric functions and products of the problem variables{|Vk |, θk }k∈N . It is proved in [24] that finding feasible flows for the AC-OPF problem
in general graphs is strongly NP-hard. However, numerical solvers can find approximate
solutions for test networks with few thousand buses.
Approximations and relaxations of the AC-OPF problem are widely analyzed. The
most simple and used one is the linear approximation, the so-called DC approximation,that will be described in the next section. [55] and [54] provide several relaxations using
Second Order Cone Programming (SOCP) and Semidefinite Programming (SDP).

2.2.2

DC-OPFIn this section we describe a linear approximation of the AC-OPF problem (2.15). Given
its simplicity and fast implementation, the DC-OPF is the problem that the economic
dispatch procedure solves every 5-10 minutes.
Under normal conditions, the following observations are made in order to simplify
constraints (2.15e)-(2.15f):• for every bus k ∈ N , |Vk | ≈ 1 (using per-unit system),• for every branch km ∈ E,
– θkm = θk − θm is small so that sin(θkm ) ≈ θkm and cos(θkm ) ≈ 1,

CHAPTER 2. PRELIMINARIES

14– resistance is much smaller that reactance, meaning that rkm  xkm . In consequence, gkm ≈ 0 and bkm ≈ −1/xkm ,sh– shunt admittance and transformer effects are ignored, that is, ykm= 0 andNkm = 1.
With these observations, expressions (2.10) for active and reactive power become:θk − θm,xkm
θk − θm
Pmk = bkm (θk − θm ) = −,xkm

Pkm = −bkm (θk − θm ) =Qkm = 0,(2.16a)Qmk = 0.

(2.16b)In other words, in the DC approximation the transmission lines are lossless (Pkm +Pmk = 0) and carry only active power and, therefore, the reactive power —and imaginary
part of variables and equations— is completely ignored. The balance equations (2.15b)and (2.15d) are simplified to∀k ∈ N ,Pkg − Pkd =X θk − θm,
xkm(2.17)

km∈δ(k)where we define as Pkg = 0 when k ∈/ G. For each choice of nonnegative values Pkg and Pkd
PP
such that k∈N Pkg = k∈N Pkg , the system (2.16)-(2.17) has a unique solution in thePkm (Lemma 1.1 in [23]).
An equivalent formulation is given by defining the vectors θ = (θk : k ∈ N ), P g =(Pkg : k ∈ N ) and P d = (Pkd : k ∈ N ) and the |N | × |N | bus susceptance matrix B as∀k ∈ N , Bkk =X
km∈δ(k)1
xkm,∀km ∈ E, Bkm = Bmk = −

1
xkm

,

(2.18)

and Bkm = 0, otherwise, we can restate the power balance system (2.17) as the |N | × |N |linear system
P g − P d = B θ.

(2.19)

Because of the construction of B, the set of equations of the linear system (2.19) is linearlyPdependent and is feasible only if k∈N (Pkg − Pkd ) = 0. By adding the equation θ0 = 0that sets the phase angle for the slack bus, the system becomes linearly independent

CHAPTER 2. PRELIMINARIES

15

when the underlying network is connected. In this case, the solution to system (2.19)
has a unique solution, of the form
θ = B̆0 (P g − P d )(2.20)

where B̆0 is an appropriate pseudo-inverse of B which depends on the choice of thereference (slack) bus.
As in the case of the AC formulation, Pkg for k ∈ G\{0} are control variables, and thePPgeneration at the slack bus is set to be P0g = k∈N Pkd − k∈G\{0} Pkg .
The DC-OPF problem is formulated as the following linearly constrained system:
minimizef (P g )subject to

P g − P d = B θ,

(2.21a)θ0 = 0(2.21b)∀k ∈ G, Pkg,min ≤ Pkg ≤ Pkg,max(2.21c)max
max∀km ∈ E, −Skm
≤ (θk − θm )/xkm ≤ Skm(2.21d)

min
max
θkm
≤ θk − θm ≤ θkm.

(2.21e)If f is a convex function, then program (2.21) is a convex optimization problem.

2.2.3

Numerical Solutions

There exist solvers that find the solution (or an approximate optimum when the problem
is difficult) to the optimization problems stated above; either in their OPF form —that
finds the optimum flow— or in their PF version, that just finds an (approximate) feasible
solution to a power flow.
Matpower [101] is a Matlab toolbox widely used in the literature that solves
optimization power flow instances and has a library of cases —usually known as IEEE test
cases— with different number of buses, from 10 to 10,000, and serves to test algorithms,
theory and their implementation.

CHAPTER 2. PRELIMINARIES

2.3

16

Time-Varying Elements

In this section we describe how the elements detailed in the previous sections evolve as
time varies. The optimization problems stated in Section 2.2 that find optimal power
flows are solved considering one particular instant of time, however, as generation mustmatch demand —plus losses— at every moment, variables and parameters evolve dynamically with time.
We introduce a time dependent argument to the following variables and parameters:• voltage Vk (t) = |Vk (t)|ejθk (t) ,• demand Skd (t) = Pkd (t) + jQdk (t) and• generation Skg (t) = Pkg (t) + jQgk (t).
In consequence, since current and power depend on the voltages, they are also timedependent, denoted as Ikm (t) and Skm (t) = Pkm (t) + jQkm (t), respectively.
In steady state, voltages and currents are sinusoidal functions of time. That is, voltageof bus k can be written as
Vk (t) = |Vk (t)|ejθk (t)with θk (t) = ωk (t) · t + δk (t),(2.22)

where ωk (t) is the voltage angular frequency —or, simply, frequency— and δk (t) the voltage phase angle referenced to ωk (t)·t. For steady state analyses it is usual to consider thatωk (t) and δk (t) slowly change through time or do not change at all and, in consequence,
they are assumed to be constant functions. Moreover, the frequency is considered equal
for all buses at ωk (t) = ω0 = 60 Hz. This is explained by the fact the power injected into
the system is produced by generators that spin at approximately 60 Hz, see Figure 2.2(detailed description of generators can be found in [13]).
The grid has been built in such a way that the generators’ rotor must synchronously
spin at a regulated frequency —60 Hz— to avoid mechanical resonances. However, inpractice, the frequency deviates from 60 Hz by small amounts, that is, ωk (t) = 60 + k (t)CHAPTER 2. PRELIMINARIES

17
statorconductor
rotor

turbinesteam

ω, Voltage, Current

energysourceFigure 2.2: Turbine-Generator Model. A turbine converts different type of energy —such as
steam or burnt natural gas— into rotation. The turbine shaft is connected to a generator shaft.
A generator typically consists in a rotor (rotating part) that spins at a frequency ω(t) inside a
stator (static part). The rotation of the rotor creates a magnetic field that induces current into
the conductor at frequency ω(t).

with |k (t)| small. In normal operation, the frequency is controlled such that |k (t)| < 0.05
Hz. Nevertheless, frequency can deviate as much as 5% for short periods without causing
load shedding, generator tripping, damaged equipment, or threatening the stability of thesystem [53].
Fluctuations of frequency are caused by generation outages, tripping lines, intermittent generation or fluctuations of demand. When demand suddenly increases, the inertia
maintained by spinning rotors —mechanical energy— is transferred into electrical energy
to overcome the power imbalance, therefore, rotors slow down and frequency decreases.
Figure 2.3 shows a contingency where 2,600 [MW] of generation is lost, as a consequence,
the frequency drops almost instantaneously but it is afterwards recovered by automatedmeasures (Automatic Generation Control).

CHAPTER 2. PRELIMINARIES

18

60.04

Frequency (Hz)

60.02
60.00

2,600 [MW]59.98

Generation Lost

59.96

Turbine−Generator inertia
instantly resists the drop
in frequency while the

59.94
59.92
59.90
5:50

AGC response

turbine speed governors
start to increase turbine
output power

6:00

6:10

6:20

6:30

Time (PM)Figure 2.3: Contingency reserves rebalance the system after the sudden loss of generation [53].

2.4State Estimation

State estimation (SE) has served as a tool to supervise and control the functioning ofa power system. Thanks to SE, the grid’s controller can take on-line decisions with theobjective of optimizing processes, adjusting generation, and detecting bad data [46].
As we have seen, the state of the system can be described by their state variables(voltages at all buses) together with the physical characteristics of the grid. A set of
sensors located through the network make measurements of some of these variables, andsend this information to the control center. A computational tool uses this data to perform SE in the following way: first, it undertakes an observability analysis that concludes
if the received data is enough to get a good state estimation of the system; next, badmeasurements are identified and eliminated of the data set; finally, an optimization problem is solved to extrapolate the measured and trusted sampling to the rest of variables
that cannot be observed.
The most advanced sensors are the PMUs (phasor measurement units) that are located on branches next to a bus. They are synchronized by a GPS clock and report,
30 times a second, data that includes time of the measurement, frequency ωk (t), voltage
phase angle θk (t) and voltage magnitude |Vk (t)| of the adjacent bus, status of the branch

CHAPTER 2. PRELIMINARIES

19(active or inactive), complex current Ikm (t) of the branch, and other data.
Traditional methods for determining the state variables are solved by using theweighted least-squares (WLS) method. For that purpose, if z is a vector of observable
quantities, then the solution x of the systemz = h(x) + e,

(2.23)

gives us the state variables (voltage phase angles and voltage magnitudes) of the network,where h is a set of nonlinear functions of the state variables and e is a zero-mean Gaussian
measurement noise vector with covariance matrix C.
Then, the SE problem findsx̂ = arg min [z − h(x)]> W −1 [z − h(x)],xwhere the weighting matrix W is usually taken diagonal and it is related with the noise
covariance matrix C.

CHAPTER 3. CYBER-PHYSICAL ATTACKS

20

Chapter 3
Cyber-Physical Attacks
In this section we introduce the notion of a cyber-physical attack on power systems.
We also model a complex attack that we demonstrate to be feasible on large networks.
Finally, we describe defense mechanisms against these and more general attacks.

3.1

Introduction

Cyber-physical attacks have been widely analyzed in the last decade. The main purposeof this study is to investigate complex weaknesses of the grid and elaborate mechanisms
to strengthen the security of the system.
The cyber adjective of the attacks correspond to the notion that sensor measurements
used to estimate the state of the system are perturbed or deleted. Voltage and current
throughout the network are periodically measured with sensors (PMUs, RTUs) and sent
to the control center of the grid. Using this data, the control center performs state
estimation, that is, using incomplete knowledge of the all the variables of the system,
unknown parameters and variables are estimated (typically by minimizing square errors)with the objective of supervising normal operation and plan-ahead future control actions(see Section 2.4 for details). If the information received by the control center is modified,actions that can be damaging to the network stability might be taken.

CHAPTER 3. CYBER-PHYSICAL ATTACKS

21

A physical attack consists on a physical alteration of the system or a portion of it, suchas changes in generation or demand (e.g. by shutting down a generator or coordinatedactivation of many air conditioners), line tripping, alteration of physical properties oflines or renewable energy generation, among others.
Thus, a cyber-physical attack has both events happening together and it is typically
performed by an adversarial agent. One of the iconic examples of a cyber-physical attack
that has recently taken place is the attack on the Ukrainian electrical system on 23
December 2015 where attackers gained access to the data acquisition system, were able
to delete and modify data, and also changed the grid topology by activating circuit
breakers [34]. In consequence, 225,000 customers lost power for a period of 1 to 6 hours.
In this work, we propose defensive techniques to be deployed when a high-fidelity
attack on a power grid is suspected. The attack is assumed to be partial in the sensethat only a subset of buses and lines are attacked, but this subset is unknown by thesystem controller. These techniques involve two ideas:(a) using network resources to randomly change power flow quantities, especially voltages and, in particular(b) changing the covariance structure of e.g. voltages in a manner unpredictable by theattacker. The specific version of this idea that we analyze introduces a low-rank
adjustment to the covariance of phase angles.
A precise definition of the attack is given in Section 3.3. These defensive techniques
focus on the phase following the initial attack, and aim to expose inconsistencies in themodified sensor data stream which is output by the attacker. We describe conditions
under which the defense succeeds in discovering the boundary of the attacked zone.
We justify such defenses by pointing out that the possibility of dangerous “cyberphysical” attacks of high-fidelity and with sparse signatures has already been indicatedin the literature (see Section 3.2). The data component of the attacks is designed to passa stringent test, namely that the falsified data satisfies the full AC power flow equations

CHAPTER 3. CYBER-PHYSICAL ATTACKS

22

at every bus and line. The data attack is coordinated with a physical attack encompassing various types (in particular, line tripping, or load modifications as considered inthis paper) that results in a dangerous system condition, e.g. a line overload. The data
modification hides this overload, with the result that sensor data received by operators isboth unimpeachable and portrays safe system operation. We term these attacks “ideal”
because, while sparse, they do assume technical sophistication and the ability to coordinate physical action and computation. Sparsity is a goal for the attacker because it
increases the likelihood of undetectability long enough for the overload to lead to line
tripping (typically several minutes). Putting aside the actual feasibility of such attacks,
the computational challenge is significant.
We will start by giving a summary of related work. Then, we develop a new optimization procedure that successfully computes cyber-physical attacks on large transmission
systems with thousands of buses. Finally, we will focus on the defensive mechanisms
against ideal cyber-physical attacks.

3.2

Previous Work

Early work on full observability of grid sensor data [66] and bad data detection has been
done in [11,30,61,68,69]. These works have developed techniques to detect bad measurements and remove them from the data set. The principal tool that they use is based on
the fact that the objective function minimized (e.g. squared errors) during state estimation becomes significantly high when arbitrary bad measurements are present. As shownin [64], these techniques do not work well when an adversary —having enough knowledgeof the system— introduces false data, since it could make the injected information lookcoherent with some state of the grid. [64] was the first to introduce the false data injection(FDI) attacks, where they use the DC power flow model for attack analysis, arguing that
it is less accurate, but more simpler and robust than the AC model [90].

CHAPTER 3. CYBER-PHYSICAL ATTACKS

23

False Data Injection Attacks
As explained before, this kind of attacks rely on the modification of the measurements
received by the control center affecting a subset of observable variables. [64, 87] study
cyber attacks where an attack vector a is added by the adversary to the observed vector
z in the system of equations (2.23) to create a malicious measurement vector z a = z + a(the vector a has non-zero entries only in the components related with the sensors thatthe attacker has access to). Ideally, the attacker expects that a vector xa minimizes theerrors of the system, with low weighted least-squares error, so that xa appears to be areasonable state variables vector. They describe results that apply when the attacker
has limited access to certain sensors or when he has a budget of how large the changes
can be. [51] also proves the existence of attack vectors that belong to a subspace of theobservable vector domain, which becomes relevant when the attacker does not have fullobservability of the sensor’s readings. In [47], this class of attacks are extended to AC
state estimation, creating unobservable attacks on specific sub-graphs of the network.
In [35] and [73] consider an attacker with no information about the grid topology.
Under certain assumptions, the attacker can learn the topology from the power flow
observations an still launch undetectable cyber attacks.
In [50] the authors study an FDI attack including topology modification, in particular,
it is consistent with a solution where branches are either not working but in reality
they are, and viceversa. The paper considers cases when either the attacker has full
information (observability) of the network or just some local knowledge. Meanwhile,
the attacks proposed in [59] based in the AC SE can lead to generation re-dispatch andline overflows. The work presented in [63, 95] analyzes FDI attacks that cause a load
redistribution in the SE using the DC power flow model.
Authors in [32] and [52] propose different algorithms to establish a minimum-cost
set of sensors should be secured in order to prevent the kind of cyber attacks describedabove. Whereas [91] defines two metrics to quantify the importance of substations (howmany different attack vectors can involve a particular sensor) and the cost of the attacks

CHAPTER 3. CYBER-PHYSICAL ATTACKS

24

on specific measurements (how many sensors are needed to be perturbed in order to alter
the measurement); and proposes strategies improve the system security with respect tothese metrics.

Cyber-Physical Attacks
The study in [80] considers a cyber-physical attack under the DC model where an adversary attacks a zone H of nodes by disconnecting some lines (without disconnecting thegraph) and obstructing the measurements within H. (They assume that after the attack
there is no edge {i, j} ∈ EH such that θi0 = θj0 .) As a first step, they show that if there isa matching between nodes outside H and inside H that covers the nodes in H, the phaseangles can be recovered almost surely (in the sense that for very specific combination ofreactance values of the edges in the matching the phase angles could not be recovered).
Secondly, the authors prove that the support of a linear system solution coincides withthe set of lines that failed (disconnected lines), and this solution is unique if and only if
the subgraph induced by nodes in the attacked zone H is acyclic. Also, if H is a cycle
and strictly less than half of its edges are disconnected, the support of a linear program
gives the set of failed lines. Other sufficient conditions are developed for planar graphs.
They extend these results in the presence of noisy data.
In [82] the authors propose an algorithm to estimate the phase angle after attacks(that trip “any number” of lines and obstruct information from the attacked zone) using
a convex relaxation of the estimation for DC power flows (relaxing equalities to second
order cone inequalities), no specific/theoretical use of AC power flow equations.
The work presented in [83] also study cyber-physical attacks affecting the physicalinfrastructure and the SCADA (Supervisory Control And Data Acquisition) system of the
network. This work is an extension of the results showed in [80] but using the AC power
flow model: an attacker disconnects lines from an area and obstructs the information
measured in that area, here generation and load remain the same after the attack for all
nodes. The paper introduces the EXPOSE Algorithm to detect line failures using the AC

CHAPTER 3. CYBER-PHYSICAL ATTACKS

25

power flow equations, the algorithm is independent of the size of the grid and the numberof line failures. Together with [82], the authors claim to have the only methods that scalewith the size of the grid to detect any number of line failures with the AC power flow
equations. Using a similar idea than in [80], they show that the complex voltages of
the attacked zone can be recovered almost surely if there is a matching that covers theattacked zone nodes by solving a linear system. Also, the authors state the same result
than in [80] when the attacked zone is a cycle. In more general cases, they propose the
EXPOSE Algorithm that solves a convex minimization problem that assuming that thevoltage magnitude of the attacked buses remain constant after the attack.
The authors in [96] study cyber-physical attacks with the following characteristics:
an intelligent attacker has access and gets knowledge of the network (topology, operation
costs, historical data, observable measurements), the attacker performs an attack that
trips a specific transmission line with the objective of overload another chosen line, thisphysical attack is coordinated with a cyber attack that masks the tripped branch andavoids its detection. The proposed strategy has two steps: first, the detection of a set ofbuses able to attack and a target line to trip; and second, the computation of an attack
vector to inject as false data. As in [60], the authors compute first a DC attack vector bysolving a linear optimization problem, and then, using AC state estimation corresponding
to the DC attack vector, they compute a more accurate AC attack vector. They test
their results in the IEEE 24-bus test case. An attack of similar characteristics using
the DC power flow model is proposed in [33] together with countermeasures to detect
the intrusion: the availability of trusted PMUs and tracking of power system equivalent
impedance. Tests in the IEEE 9-bus, 14-bus, 30-bus, 118-bus, and 300-bus test powersystems are shown.
[58] also study cyber-physical attacks that, in this case, the physical attack preserves
the grid topology but redistributes loads (see also [63, 95]) within certain budget acrossthe network and hides it through manipulated injection of data, under the DC model.
The paper bases the analysis on the computation of the generation shift factors that

CHAPTER 3. CYBER-PHYSICAL ATTACKS

26

capture sensitivity of power flow changes on lines when bus injections are altered. Twostage mixed-integer linearized optimization problems are proposed to obtain the set oflines and buses that will be affected by the attack. Given the dimensionality of the
problem, the authors show results of coordinated cyber-physical attacks on small IEEE
test systems with 14 and 118 buses.

Other Related Literature
The work in [28] proposes a game theory model in where an attacker and a defender ofthe network have limited budget that allow them to attack targeted elements of the grid
(trip a line, break a transformer) o to secure them, respectively. They implement analgorithm and show that when defenders deploy an strategy before an attack is initiated,
the loss can be predictable and limited to a minimum level.
Other works have focused on the system vulnerability to line outage, see e.g. [14,18, 29, 39, 72, 85, 86, 97, 98, 100]. The objective of these studies is to identify from thedata when line tripping has happened, but also analyze the consequence that they could
caused, for example, cascade line failures. In particular, [72] uses a bilevel mixed integer
nonlinear programming to formulate the problem of finding a small set of lines whose
removal from the network would cause severe blackouts.

3.3

Attack Model

We will now describe and formulate the cyber-physical attacks that we considerate, combining physical disruption and data intrusion. We will show that the numerical solution
to this problem scales well to systems with thousands of buses, with running times inthe tens of seconds or less on a standard computer. These experiments do not prove thatthe attacks, though sparse, could easily be executed. Rather they show that the attack
computation is tractable, thus providing added justification for studying sophisticatedand scalable defense mechanisms.

CHAPTER 3. CYBER-PHYSICAL ATTACKS

27

Attack models considered in prior work allow the attacker different capabilities. Regardless of the model, Template 1 given below (similar to one in [96]) broadly outlines the
structure of an AC-undetectable attack. We use the term “initial” to indicate that theattack comprises actions taken at one point in time. Later we will discuss a “follow-up”
phase that follows the initial attack.

Template 1. Initial Attack(a) It is assumed that at each bus k there is a sensor measuring voltage at k andcurrent at each line km ∈ δ(k).
(b) The attacker has selected a (sparse) subset A of buses, as well as a target line
uv within A that will be overloaded.
(c) For any bus k ∈ A, the attacker can modify data provided by a sensor locatedat k.
(d) The attacker’s physical actions are of two types. First, the attacker can modify
loads at buses in A. Additionally the attacker can disconnect lines with bothends in A.
(e) Actions (c)-(d) are performed in a single step.
(f) The data received by the control center satisfies complete fidelity as per AC
power flow equations (2.7)-(2.10) and shows all system limits being satisfied,
while in actuality line uv is overloaded.
(g) When the attack includes load changes, secondary response (i.e. AGC response) is taken into account by the attacker.

Conditions (a), (f) and (g) amount to a strong form of undetectability. Nevertheless,
we provide examples of large scale systems that are susceptible to attacks of the form

CHAPTER 3. CYBER-PHYSICAL ATTACKS

28(a)-(g). Note that we allow loads to be modified, but not generation. In our numericalexamples we enforce that G ∩ A = ∅, out of a perception that generator sites are more
carefully protected.
Point (e) requires some discussion. Immediately following any physical modification
to a system, we can expect a change in voltages (magnitudes and phase angles) and evento system frequency, the latter especially when net loads are changed. More properly,system dynamics will undergo a change. Understanding the precise nature of that change
is a substantial computational task. The current state-of-the-art involves a numericalsimulation that alternates between simulation of true dynamical behavior at generators(the so-called swing equation) with AC power flow updates. This combined computation
will typically run much slower than the actual dynamics, and assumes correct knowledgeof the underlying transmission system. Under adversarial attack that e.g. modifies the
topology, the rapid success of such a computational approach to identifying the current
grid state seems uncertain. And once action (d) (i.e. modification of sensor data) is taken
a completely falsified, and consistent view of the system is being presented.
We next present conditions that we will impose so as to guarantee undetectability.
True data will be the true physical data. In contrast, reported data is that which isactually received by the control center and includes the attacker’s modifications. The
true data will be given by the (voltage, current) pair of vectors (V T , I T ) whereas thereported data will be given by (V R , I R ).
An important requirement for the reported data is




R
R
I
V
 km  = Ykm  k  ,current-voltage consistency:
R
Imk
VmR

(3.1)

i.e. equation (2.8), should be satisfied at all branches. This condition will be enforced in
the computation given below in an indirect fashion (also see [67] for a different use of thisrequirement). In general, of course, an attacker might only seek approximate consistency,
using ambient noise to hide errors. Additionally to (3.1):(s.1) On a bus k ∈/ A the true and reported data agree (no data modification outside A,

CHAPTER 3. CYBER-PHYSICAL ATTACKS

29

by definition).
(s.2) At a boundary bus1 k ∈ ∂A the attacker is constrained by the condition VkR =VkT . This condition is applied to avoid attack detection, given (a) and the second
equation in (3.1) applied to a line km where m ∈/ A.
(s.3) On buses k ∈ A\∂A we may have VkR 6= VkT and on lines with at least one end inA\∂A the true and reported currents may also differ.
(s.4) The reported voltages and currents must be consistent with meaningful (complex)PR∗
equals thepower injections. Specifically, consider a bus k. Then km∈δ(k) VkR Ikm
power injected into the system at bus k, according to the reported data. If k ∈/APT∗which isby definition (of reported and true data) this sum equals km∈δ(k) VkT Ikm
the true power injected by bus k. On the other hand if k ∈ A the sum may differ
from the true injection at k.
(s.5) If the attack causes a net change in the sum of loads, the resulting AGC-mandatedchange in generator output must be taken into account.
We call condition (s.4) power-injection consistency, that is:
XR∗
VkR Ikm= net injection at k

∀k,

(3.2)

km∈δ(k)where “net injection” is the reported net injection on buses in A and the true net injection
for buses not in A.
Subject to these requirements, the attacker seeks to create a (true) line overload on
uv with both ends in A, while the reported data shows safe system operation (voltage,
angle, and power flow limits are satisfied). In other words, for η > 1 (e.g. η = 1.2 or
1Recall that k ∈ A is a boundary bus if it has a neighbor outside of A.

CHAPTER 3. CYBER-PHYSICAL ATTACKS30

η = 1.5), the attacker would try to satisfy
T 2
2
T 2
max 2
(Puv
) + (QTuv ) = (Suv ) > (ηSuv )∀km ∈ E,∀k ∈ N ,

max 2
2
R 2
R 2
) + (QR(Pkm
km ) = (Skm ) ≤ (Skm )

(3.3)
(3.4)

min
R
max
θkm
≤ θkR − θm
≤ θkm

(3.5)

Vkmin ≤ |VkR | ≤ Vkmax

(3.6)

R
R
T
, QT
where (Pkm
km ) and (Pkm , Qkm ) are the true and reported flows given by the expres-sions in (2.10) using {VkT }k and {VkR }k , respectively. In the next section we present amathematical formulation for this problem.

3.3.1

Formulation

As input to the problem we have a set A ⊂ N \G of buses, a set of lines L to bedisconnected, all with both ends in A, and a line uv ∈/ L with both ends in A. Write
AC = N \ A. Let (Ŝkg = P̂kg + j Q̂gk )k∈N and (Ŝkd = P̂kd + j Q̂dk )k∈N be the complex power
generation and loads, respectively, at the time of the attack. We assume that the attacker
observes all these quantities.
The initial attack problem uses the following real-valued variables, where “T” indicates true and “R”, reported:• ∀ bus k ∈ N : |VkT |, θkT denote the true voltage magnitude and voltage angle;• ∀ bus k ∈ N : |VkR |, θkR denote the reported voltage magnitude and voltage angle;• ∀ bus k ∈ A: Pkd,T , Qd,T
denote the active and reactive true loads in A;k• ∀ bus k ∈ A: Pkd,R , Qd,R
denote the active and reactive reported loads in A;k• ∀ bus k ∈ R: Pkg , Qgk denote active and reactive generation at participating buses;T• ∀ line km ∈ E\L : Pkm, QT
km denotes the active and reactive true power flows;R• ∀ line km ∈ E: Pkm, QR
km denotes the active and reactive reported power flows;

CHAPTER 3. CYBER-PHYSICAL ATTACKS

31• ∆ denotes the net change in active power generation.
Therefore, the problem is formulated as follows:
T 22
maximize (Puv
) + (QT
uv )(3.7a)

subject to ∀k ∈ AC ∪ ∂A, |VkT | = |VkR |, θkT = θkR
X
R
)
=
(Pkm
+ jQR
∀k ∈ A, −(Pkd,R + jQd,Rkm ),k(3.7b)(3.7c)

km∈δ(k)− (Pkd,T + jQd,T
k ) =

X

T
+ jQT
(Pkm
km ),

(3.7d)

km∈δ(k)\L

Pkd,R ≥ 0, Pkd,T ≥ 0(3.7e)X∀k ∈ AC \R, P̂kg − P̂kd + j(Q̂gk − Q̂dk ) =

T
(Pkm
+ jQT
km )

(3.7f)

km∈δ(k)X∀k ∈ R, Pkg − P̂kd + j(Qgk − Q̂dk ) =

T
(Pkm
+ jQT
km )

(3.7g)

km∈δ(k)Pkg − P̂kg = αk ∆(3.7h)∀k ∈ G, Pkg,min ≤ Pkg ≤ Pkg,max , Qg,min≤ Qgk ≤ Qg,maxkk(3.7i)∀k ∈ N , Vkmin ≤ |VkT |, |VkR | ≤ Vkmax(3.7j)

Rmax
Tmax∀km ∈ E, |θkR − θm| ≤ θkm; |θkT − θm| ≤ θkm
if km ∈
/ L,

(3.7k)

R
R
R
max
max{ k(Pkm, QRkm )k , k(Pmk , Qmk )k } ≤ Skm ,(3.7l)

T
T
T
T T
Pkm
+ jQT
/L
km = Skm (|Vk |, |Vm |, θk , θm ) km ∈(3.7m)

T
T
T
T
T
Pmk
+ jQT
/L
mk = Smk (|Vm |, |Vk |, θm , θk ) km ∈(3.7n)

R
R
R R
R
Pkm
+ jQR
km = Skm (|Vk |, |Vm |, θk , θm )(3.7o)

R
R
R
R
R
Pmk
+ jQR
mk = Smk (|Vm |, |Vk |, θm , θk )(3.7p)In this formulation, power flows are represented through the quadratic functions
Skm , Smk (see equations (2.9)-(2.10)) which appear in the formulation as (3.7m)-(3.7p).
Note that we include voltage variables but no current variables. However, having solved
the above optimization problem, the attacker reports, for each line km with both ends

CHAPTER 3. CYBER-PHYSICAL ATTACKS

32

R
R) computed using the formula, Imk
in A, a current pair (Ikm




R
R
Ikm
|VkR |ejθk = Ykm 
,
R
R
R jθm
Imk
|Vm |ethereby attaining current-voltage consistency (3.1). Note that if either k ∈ ∂A or m ∈ ∂A
the true and reported voltage values are identical —see Lemma 3 below.
Lemma 1. Consider a feasible solution to problem (3.7). Let H denote either T or R
H

(i.e. true or reported). Then the voltages |VkH |ejθk for all k ∈ N yield a solution to thepower flow problem where(1) Bus k has load Pkd,H + jQd,H
for k ∈ A and P̂kd + j Q̂dk if k ∈ AC .
k(2) Bus k ∈ G has generation Pkg + jQgk if k ∈ R and P̂kg + j Q̂gk if k ∈ G \ R.
H(3) Line km has power flow Pkm+ jQH
km when H = R and also when H = T and

km ∈
/ L.
(4) When H = R (reported data) the solution is fully feasible, i.e. it satisfies voltage,
generation, phase angle and power flow limits.
(5) When H = T (true data) the solution satisfies voltage, generator and phase angle
limits, but only satisfies power flow limits on lines km with both k, m ∈ AC ∪ ∂A.
The solution is also consistent with lines in L being cut.
Proof. Property (3) follows from constraints (3.7m)-(3.7p). Hence, (1) and (2) follow
from constraints (3.7c)-(3.7f). Properties (4)-(5) follow from constraints (3.7i)-(3.7l).
As a corollary to (1)-(2) of Lemma 1, a feasible solution to problem (3.7) satisfies,
exactly, power-injection consistency, i.e. condition (s.4) above.
Lemma 2. Consider a feasible solution to problem (3.7). The solution is consistent witha secondary-response adjustment of active power generator amounting to ∆ units.
Proof. Follows from constraint (3.7h).

CHAPTER 3. CYBER-PHYSICAL ATTACKS

33

Lemma 3. Consider a feasible solution to problem (3.7). Then (a) the true and reported
voltages agree on AC ∪ ∂A. Further, (b) the true and reported currents on a line km are
identical if k, m ∈ AC ∪ ∂A.
Proof. (a) Follows from constraint (3.7b), and (b) is a consequence of (a).
Corollary 4. Suppose we compute a feasible solution to problem (3.7) whose objective
max 2) . Then the reported solution amounts to an undevalue is strictly greater than (Suv

tectable attack that hides an overload on line uv.

3.3.2

Computational Viability

Above we have presented a mathematically correct version of the initial attack problem
that would lead to an (initially) undetectable attack, via problem (3.7) which is a nonlinear, nonconvex optimization problem, and thus, in principle, a challenging computationaltask. Nevertheless this problem is similar to the standard ACOPF or PF problem and(at least) a local optimum should be efficiently computable; this expectation is borneout by our experiments. Strict maximization in (3.7) is not required for an attack to
be successful (all that is needed is an overload of the line uv), it is easy to see that theobjective function (3.7a) can be replaced by a feasibility constraint
T 2
2
max 2
(Puv
) + (QT> (ηSuv),
uv )for some η > 1, and solve a feasibility problem instead of an optimization problem.
A broader issue concerns the selection of the sets A and L. This is a combinatorial
problem which is bound to be intractable. In fact [81] describes a number of strong
NP-hardness results in the DC setting, e.g. given vectors of phase angles θ and θ0 it is
NP-hard to compute a set L such that B 0 θ0 = Bθ where B 0 is the bus susceptance matrixof the network with L removed.
Nevertheless, as discussed in the literature, an attacker may be willing to incur significant computational costs in order to compute a successful attack. While it is reasonable

CHAPTER 3. CYBER-PHYSICAL ATTACKS

34

to assume that an attacker’s ability to take physical action or to modify data is limited(see the discussion in [52,64]), not assuming computational intelligence on the part of an
attacker amounts to a limitation on the part of the defender.
We separate two distinct issues here: first the identification of the set A, which
is done in advance and may be computationally intensive, and second, the solution to
problem (3.7) which only requires a few seconds. Let us assume that the attacker has had(undetected) access to system and sensor data long enough to identify a weak sector of thetransmission system, i.e. the set A. In this task the attacker would rely on the fact thattypical (time- and day-dependent) load and generation profiles for transmission systems
are statistically predictable with some accuracy. This fact would help the attacker inthe computation of a target set A, perhaps using enumeration, using load estimates inproblem (3.7).
Having identified a particular set A, problem (3.7) would be run once again just prior
to the attack, now using close estimates of the loads obtained from ambient conditions.
Assuming that the attack is perpetrated during a period of slowly changing loads, and not
close in time to a generator redispatch2 , the attack will likely be sufficiently numericallyaccurate so as to become difficult to detect.

3.3.3

Computational Implementation

Our experiments and implementations are based on the Matpower test cases. Matpower [101] is an open-source package for Matlab that solves optimization problems andpower system simulations, such as the DC- and AC-OPF. Matpower has a built-in optimization solver called MIPS (Matpower Interior Point Solver), a primal/dual interior
point method, but other solvers and toolboxes are also supported, such as fmincon (fromthe Matlab Optimization Toolbox), IPOPT [92], GUROBI [43], CPLEX, MOSEK [6],and others. These solver can be used with Matpower if they are properly installed for2

A generator redispatch is a request made by the control center to a generator to adjust its real power

injection to the network.

CHAPTER 3. CYBER-PHYSICAL ATTACKS

35

running in Matlab.
In other to implement the type of attacks that have previously been described, we
have coded our own optimization/feasibility problems that find solutions to power flow
instances. Our code reads test case files that are structured in the same shape as Matpower test cases, and creates a Matlab file where the state variables are the ones
described in the previous sections and the objective function and constraints, with theircorresponding Jacobian and Hessian, are explicit (in contrast with Matpower functions that are difficult to read and modify). The code can be found in the website
github.com/me2533/acopf and its performance is in some cases better than Matpower.

3.3.4

ExamplesIn the following instances we consider the case2746wp and the case1354pegase (thathave 2746 and 1354 buses, respectively) from the Matpower case library [101].
For the first case, the adversary attacks the set of buses
A = {1137, 1138, 1139, 1141, 1361, 1491}
with

A1 − ∂A1 = {1137, 1138, 1141, 1491}.

See Figure 3.1. In this attack the quantity ∆ in (3.7h) equals 135.09. We also have L = ∅(no lines are cut). The set of generators participating in secondary response is R1 = {17,18, 55, 57, 150, 383, 803, 804, 1996} with participating factors αk = 1/9 for all k ∈ R1 .
Table 3.1 shows the true and reported flow for lines where the solutions differ, witha strong overload on line (1361, 1141) and (1138, 1141).
Table 3.2 displays the load and generation of the buses involved in the attack.

CHAPTER 3. CYBER-PHYSICAL ATTACKS

36

1491
1139

1137

1361

1141
1138

1110

85

1625

1844

1512

1287

1252

1295

104

1651

1358

Figure 3.1: Attacked zone A1 (in red) and its neighborhood. Generators are shown in blue.

Table 3.1: True and reported flow for attacked lines. (Overflow in bold)

bus k

bus m

1139

1137

1361

1141

1141

1491

1141

1138

pT
km
pR
km

T
qkm
R
qkm

3.36 2.66
3.36 2.66
229.01 10.49
108.51 10.49
13.46 2.41
6.20 2.39
209.25 4.44
98.06 5.24Tk(pTkm , qkm )kRk(pR
km , qkm )k

4.29
4.28
229.25109.02
13.68
6.64
209.29
98.20

max
Skm

114.00114.00114.00114.00CHAPTER 3. CYBER-PHYSICAL ATTACKS

37

Table 3.2: Load and generation before and after the attack.

bus k
1137
1138
1139
1141
1361
1491

Before Attack
P̂kd
Q̂dk
0
0
103.29
29.84
0
0
0
0
0
0
4.76
1.12

3.36
208.91
3.36
5.97
1.58
13.45

After Attack
Qd,T
Pkd,R Qd,Rkk
2.68 14.74 1.37
3.32 20.82 1.46
2.68 15.80 1.37
2.64 24.20 1.45
2.39 91.90 1.45
2.65 20.58 1.48

gen k
17
18
55
57
150
383
803
804
1996

P̂kg
Q̂gk
140.00 120.00
140.0041.01
130.00 −20.00
130.00 −20.00
90.00
0
21.27
10.54
0
0
0
10.00
90.00
62.86

Pkg
155.01155.01145.01145.01105.01
36.28
15.01
15.01
105.01

Qgk
116.90
61.33
19.29
38.33
26.65
13.78
6.47
8.27
77.28

Pkd,T

αk ∆
15.01
15.01
15.01
15.01
15.01
15.01
15.01
15.01
15.01

CHAPTER 3. CYBER-PHYSICAL ATTACKS38

For the second instance we consider the case1354pegase and two different attacks.
First
A2 = {44, 367, 1027, 1172, 1833, 1923, 1973, 2372, 2458, 2644, 2919, 2928, 3243,
3543, 3610, 3657, 3855, 4885, 5308, 5477, 5648, 7148, 7865, 8104, 8722, 9191}
with responding generator set R2 = {352, 757, 1794, 2421, 2816, 4918, 7267, 7808}. See
Figure 3.2. The participating factors are the same for all generators.
4432
8361
196

3654

8507
5131

7865

2928

1394
1973

3855

3657

4885

4157

5477

2273

3962
2644

9191

1172

8722

8104

3243

7373

3610

2372
5308

8628

1027

8879
8347

2458
8560

1545

2230
1838

3543

367
2919

1833

2421

44
7148

3541
5648

1923

7273
3830

6313
5837

6989

6427

8030

5610
1758
7267

7537
1465
306

Figure 3.2: Attacked zone A2 (in red) and its neighborhood. Generators are shown in blue.

Table 3.3 shows the true and reported flow for lines where the solutions differ —with
overloads of 53%—, while Table 3.4 displays the load and generation of the buses involvedin the attack.

CHAPTER 3. CYBER-PHYSICAL ATTACKS

39

Table 3.3: True and reported flow for attacked lines. (Overflow in bold)

bus k

bus m

5477

3657

5477

2928

367

1172

367

9191

367

8722

3657

3855

3657

3855

3657

2928

3657

2928

3657

1172

3657

9191

3657

8722

8104

8722

2644

9191

4885

2928

3855

1973

3855

1973

2928

7865

pT
km
pR
km

T
qkm
R
qkm

−75.2415.42−119.01
18.14
75.24
−6.50
74.35
−8.66
384.37 −106.48
401.00 −104.03
341.06 −89.26
404.27 −80.36
809.96 −66.76
418.41 −81.99
146.49 −48.88
192.83 −55.15
162.87 −48.21
214.05 −53.35
108.03 −27.27
127.77 −32.80
112.07 −31.88
132.54 −38.27−382.70
118.57−298.52
124.10−339.8099.37−271.11103.01117.8099.39−277.51105.57−0.29
0.24−67.70
−6.08
0.00
0.23
−66.39
−6.01
0.00
−0.00
−13.49−0.63
195.73
6.21
215.31−5.72
113.54−6.52
124.26 −14.44
−0.00
0.00
14.67
−2.73

T
k(pT
km , qkm )kRk(pR
km , qkm )k76.80
120.3975.52
74.85
398.85414.28352.55412.18812.70426.37154.43200.56169.85220.59111.42131.92116.52137.95400.65323.29354.04290.02
154.12296.91
0.37
67.97
0.23
66.67
0.00
13.50
195.82215.39113.73125.09
0.00
14.92

max
Skm

491.00433.00529.00529.00529.00529.00491.00395.00376.00491.00453.00414.00491.00∞281.00529.00491.00395.00CHAPTER 3. CYBER-PHYSICAL ATTACKS

40

Table 3.4: Load and generation before and after the attack.

bus k
367
1027
1172
1923
1973
2458
2644
2928
3610
3657
3855
4885
5308
5477
7148
7865
8104
8722
9191

Before Attack
P̂kd
Q̂dk
65.17 −51.50−5.91 −2.33
00192.1145.50
210.1066.20
132.00
26.00
0
0
0
0
185.8031.80
131.001.30−38.73 −15.32
166.40
50.40
0
0
0
0
0
0
84.99
19.60
127.20
40.80
0
0
0
0

0
0
0
166.170
128.88
0
0
0
0
0
0
0
0
0
0
0.29
924.58
0

gen k
352
757
1794
2421
2816
4918
7267
7808

P̂kg
870.56120.00777.83100.00746.1980.00
100.00797.28Pkg
Qgk
867.39 −297.68
116.8251.30
774.66 −288.66
96.82
37.39
743.02 −297.74
76.82 −12.57
96.82
77.34
794.10926.35

Q̂gk
22.55
51.30
−3.27
−12.45
7.62
37.27
23.18
−1.84

Pkd,T

After Attack
Qd,T
Pkd,R Qd,Rkk
0 311.70 3.89
0000 100.80 8.15
0 166.17049.00 30.27 29.06
0 128.88
00 66.39 6.24
14.22 11.10 2.03
0000 60.91 10.61
20.56 67.17 28.79
1.49 13.49 2.12
11.08
0 11.08
12.88 44.67 12.32
63.08
0 63.08
5.23 14.67 2.49
0 67.70 6.32
0 71.99 6.16
0 65.55 6.07
αk ∆
−3.18
−3.18
−3.18
−3.18
−3.18
−3.18
−3.18
−3.18

CHAPTER 3. CYBER-PHYSICAL ATTACKS41

For the last attack example, also in the case1354pegase, we set
A3 = {174, 305, 953, 1035, 1311, 1817, 1965, 2365, 2526, 3579, 3613, 3649, 3697,
3794, 4504, 4874, 5106, 5469, 6555, 6901, 7903, 7905, 8180, 8373, 8748, 8931}
with responding generator set R3 = {564, 1001, 7466}. See Figure 3.3. The participating
factors are the same for all generators.
7138
6337

1237

8128

3649

5379
3579

5469

7466

1483

8931
8494

4504
6901
1730

953

7903
6199

5837

4874

2526

564

8180

6738

2252

1817
5106

2365

1035
305

8316
3613

1311

3794
8475

174
4141

8486

5589

8748
6555

6430

5213

891

1502

8373

7905
3608
3697

7541

7530

5926
1965

3022
1552

6952

4426
1265

Figure 3.3: Attacked zone A3 (in red) and its neighborhood. Generators are shown in blue.

Table 3.5 shows the true and reported flow for lines where the solutions differ —withline overloads of 66%—, while Table 3.6 displays the load and generation of the buses
involved in the attack.

CHAPTER 3. CYBER-PHYSICAL ATTACKS

42

Table 3.5: True and reported flow for attacked lines. (Overflow in bold)

bus k

bus m

4874

4504

4874

3579

4874

3579

4874

7903

4504

7903

3579

5469

6901

4874

pT
km
pR
km

T
qkm
R
qkm

−54.84127.62−241.102.65−254.86
3.09−202.23304.73−56.8918.03−498.23−120.65
953.11569.95−17.66
−6.05
21.74
16.17
27.21
17.07
−56.32−37.51−13.90−16.3934.98
21.26
236.71
88.70

T
k(pT
km , qkm )kRk(pR
km , qkm )k57.61
127.76242.0816.39
256.3117.35
209.92307.0358.56
24.37
499.46122.51982.06576.81max
Skm

453.00376.00∞338.00∞491.00591.00CHAPTER 3. CYBER-PHYSICAL ATTACKS

43

Table 3.6: Load and generation before and after the attack.

bus k
174
305
953
1035
1311
1817
1965
2365
2526
3579
3613
3649
3697
3794
4504
4874
5106
5469
6555
6901
7903
7905
8180
8373
8748
8931

Before Attack
P̂kd
Q̂dk
370.50 −17.00
103.0012.00
108.70 −12.15
382.50
57.90
37.94
15.20
108.9133.70
248.70
8.10
37.94
15.40
238.4546.86
277.727.35
152.5014.40
00−32.71 −13.76
165.307.40
170.10 −9.30−0.58 −53.03
39.94
13.00
106.60 −18.80−37.81 −21.64−37.84 −19.23
00347.1087.20
291.6047.00
216.30 −6.90
279.70
4.60
−0.56
0

0
0
0
0
0
0
0
0
0
0
0
155.13110.410
1.831702.23
0
0
0
0
59.13
676.25000712.85After Attack
Qd,T
Pkd,R
k
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0 126.38
000 155.13
40.10 110.41
000 108.55
0 130.51
000 385.57
000 383.16
0 639.63
0 676.25
0
0
0
0
0
0
0 712.85

gen k
564
1001
7466

P̂kg
100.00120.00100.00Q̂gk−19.3555.02−59.77

Pkg
60.04
80.04
60.04

Qgk
αk ∆
18.47 −39.96
55.02 −39.96
75.85 −39.96

Pkd,T

Qd,R
k
0
0
0
0
0
0
0
0
0
15.21
0
0
40.10
0
7.80
12.96
0
36.48
0
148.01
6.27
0
0
0
0
0

CHAPTER 3. CYBER-PHYSICAL ATTACKS

3.4

44

The Follow-Up Phase

Following the initial attack, the attacker needs to dynamically perturb the attack dataso as to produce a realistic data stream that is both consistent and continues to hide theoverload.
In analogy to our notation for the initial problem, at time t > 0 following the attack,
we denote by VkR (t) and VkT (t) be reported and true voltages at t and similarly with currents. Reported data for A will be manufactured by the attacker aiming to approximately
satisfy current-voltage consistency (3.1) and power-injection consistency (3.2).
In addition, in this work we assume that the attack is perpetrated when ambientconditions (in particular loads) are, approximately, constant. Let us denote by VkR (0)the voltage at a bus k computed by the initial attack (3.7), i.e.
R.
VkR (0) = |VkR |ejθkRand likewise define the current Ikm(0) on line km. The statement that ambient conditions

are approximately constant, post-attack, can be informally rephrased asRR
VkR (t) ≈ VkR (0) ∀ k, and Ikm(t) ≈ Ikm(0) ∀ km.

(3.8)

If ambient conditions are approximately constant, (3.8) will hold (statistically) for any
bus k and line km not in the attacked zone A but are otherwise a requirement for theattacker.
Two types of attack have been used in the literature. First, the “noisy data” attackin our setting works as follows:

CHAPTER 3. CYBER-PHYSICAL ATTACKS

45

Template 2. Noisy Data Attack
At time t > 0 the attacker reports at each bus k ∈ A a voltage
VkR (t) = VkR (0) + νk (t).
Here νk (t) is a random phasor drawn from small variance, zero mean distributiona .
Likewise the attacker reports for each line km with both ends in A, currents



 

R
R
Ikm
(t)
Ikm
(0)
µkm (t)

 = 
+

(3.9)
R
R
Imk (t)Imk (0)µmk (t)where µkm (t), µmk (t) are drawn from zero mean distributions with small variance.
a

We use boldface to indicate random variables.

Note that these definitions satisfy requirement (3.8), and approximately satisfy currentvoltage consistency. As a functionally equivalent alternative to (3.9) the attacker could
simply set



R
Ikm
(t)
R
Imk
(t)





 = Ykm VkR (t)VmR (t),(3.10)

our analyses below apply to either form.
A second form of attack that has been considered is the data replay attack. Herethe attacker supplies a previously observed (or computed) pair of time series V R (t) andI R (t) for buses and lines within the set A.

Discussion
The reader may recall that in the initial attack computation we enforced that reported
voltages in ∂A are exact, i.e. equal to the true voltages (constraint (3.7b)). In the timevarying phase this condition is necessarily relaxed by the attacker, though this action carries the risk (to the attacker) that current-voltage consistency will not hold, statistically,

CHAPTER 3. CYBER-PHYSICAL ATTACKS

46

for some line km with k ∈ ∂A and m ∈/ A. Thus e.g. in the noisy data attack template
given above the distributions for the νk (t), µmk (t) and µkm (t) should have sufficiently
small variance relative to the variance of ambient conditions. Further requirements onsuch variances will be discussed in Section 3.5.3. In any case, when ambient conditions(e.g. loads) are nearly constant, the noisy-data attack may continue to approximately
satisfy current-voltage and power-injection consistency and thus remain numerically undetectable. The same holds for the data replay version provided the replayed voltages in∂A closely approximate ambient conditions.
In the next section we present defensive mechanisms that dynamically change voltages
in a way that is unpredictable by the attacker. The key observation is that a substantialchange to voltages in ∂A will cause the noisy-data attack, applied verbatim as in Template
2, to fail, because of large current-voltage inconsistencies on lines km with k ∈ ∂A andm ∈/ A. Of course, the template need not be applied verbatim, and in particular theattacker may seek to leverage the possibility of sensor error. We will consider this pointin the next section.
In [83] current-voltage consistency is used in a different setting: (i) the attacked zone
A is known by the defender, (ii) the attacker only disconnects lines. Under a number ofassumptions, in particular that there is a matching between AC and A that covers all
buses in A it is shown that the attack can be accurately recovered.

3.5

DefenseIn the above sections we showed that, conceptually at least, it is possible to computehigh-fidelity attacks that disguise dangerous network conditions. Other attacks are also
potentially conceivable, e.g. impedance changes, transformer tap changes, etc. In thissection we describe a generic randomized defense strategy that can deployed when acomplex attack is suspected. We will assume that the attack impacts a proper subset A
of the system that is unknown to the control center, as was the case above, though theCHAPTER 3. CYBER-PHYSICAL ATTACKS

47

generic defense strategy applies under more general attacks as well. The strategy can besummarized by the following template:
Procedure 3. Random Injection Defense.
Iterate:D1: Choose, for each k ∈ G a (random) value δk such thatPg∈G δk ≈ 0. Commandeach generator k ∈ G to change its output to Pkg + δk .
D2: Following the generation change in step D1 identify inconsistencies in theobserved sensor readings.

Here, an “inconsistency” is a difference between a voltage value that is reported by theattacker, and a corresponding value that predicted by the defender. We will describe
several concrete versions of this idea below. See Procedures 4, 5 and 6.
Each iteration would last several seconds, and statistically significant inconsistencies
identified by this scheme would be flagged as potential evidence of an attack. Below wewill describe several specific implementations of the random ingredient; randomness isused because the attacker cannot anticipate the random injections and thus will not beable to instantaneously update the sensor readings within A. The above strategy could
be AGC-like if only generators k ∈ R (the responding generators) are allowed to have
δk 6= 0 and in general it amounts to a generator redispatch. The strategy in Procedure 3
is likely to succeed, in particular against the noisy data or data replay attacks, if thegeneration changes result in significant voltage changes across the system. Lemma 5
given below explains why a particular implementation of Procedure 3 attains this goal.
An additional point is that an implementation of step D1 should guarantee safe system
operation; this consideration leads to computation of the δk in step D1 by means of an
OPF-like problem.
We note that there is an existing literature on using network resources so as to change
power flow physics in order to detect structure or faults. See [15–17,89,99]. Indeed, even

CHAPTER 3. CYBER-PHYSICAL ATTACKS

48

though the description of our random defense focuses on power injections, one could
also consider other random probing strategies that change power flows, such as adjustingtransformer settings, controlled line tripping, and the use of DER (distributed energy
resources)3 , storage and FACTS (flexible AC transmission system) devices.
There are several implementations of the generic strategy. Generally the defenderwants to make the |δk | large because to first order changes in voltage angles are proportional to kδk2 , and a large change in phase angles is likely to give rise to a significant
current-voltage or power-injection inconsistencies in sensor readings in ∂A, as discussed
above. This idea forms the basis for a simple, current-consistency based version of Procedure 3 given in Section 3.5.2.
An attacker aware that the random defense strategy is applied may try to replacee.g. the noisy data attack with a more careful manipulation of reported data. For example, the attacker could react to a significant change to voltages in ∂A by solving anonlinear, nonconvex system of inequalities designed to guarantee approximate currentvoltage and power-injection consistency. In addition, any implied load change within A
must be very small (or it would contradict observed frequencies). Finally the attacker
would need to perform this computation very quickly, and repeatedly (because the defense will be applied repeatedly).
This online complex computation could in principle be bypassed by the attacker byconsidering changes to readings of voltages at buses in ∂A only; with the remaining
voltages in A computed as in Template 2. We will term this the enhanced noisy dataattack. We remark that the adversary would still have to maintain AC consistency forlines within A, which is nontrivial. Nevertheless, the ability to adjust readings in ∂A
beyond what is prescribed by the noisy data attack may provide some flexibility for theattacker. However, in Section 3.5.2 we will show that when the random injection defense
causes large-enough voltage changes in ∂A, the enhanced noisy-data attack fails. See
Lemma 10.
3

Small local artifacts that generate or store energy, located close to the load point that they serve.

CHAPTER 3. CYBER-PHYSICAL ATTACKS

49

A more practicable alternative (for the attacker) would be to consider arbitrary
changes to voltages in buses in ∂A, with the remaining voltages in A obtained as inTemplate 2. We will term this the enhanced noisy data attack. In Section 3.5.2 wewill show that when the random defense causes large-enough voltage changes in ∂A, theenhanced noisy-data attack fails. See Lemma 10.
A more sophisticated defensive idea, given in Section 3.5.3, changes the stochastics ofpower flow data, in particular voltage covariance, and probes the corresponding propertiesof the reported data.
Our defensive strategies can be easily adjusted if sensors are not available throughout
the system. Of course, the fewer the sensors the more limited the impact of the defense. Indeed, some interesting work (using the standard, DC-equation state estimation)
precisely seeks to perform system identification post-attack when only limited sensor
information is available [80–83]. Note that the attack problem becomes easier (for theattacker) if sensors are not widespread. The attacks computed in Section 3.3.4 do assume
sensors at every bus, and yet succeed even in large-scale cases.

3.5.1

Controlling Voltages Through Generation Changes

As discussed above, a goal of the defense is to produce large voltage angle changes in busesin A, with the intention of revealing inconsistencies in reported data on lines between ∂A
and AC . The defender, of course, does not know the set A and thus it is of interest tounderstand when the voltage at any given bus can be changed by appropriately choosing
the injections δ.
In this section we address the task of changing voltage angles through injections.
First, we will argue by using the DC power flow approximation (2.19), that Procedure 3
does succeed in changing phase angles (see Lemmas 5 and 6). In Section 3.5.1.1 we will
present experiments under the AC power model that verify the DC-based results. And
in Section 3.5.2 we further argue that the voltage changes are large enough to overcome
sensor error.

CHAPTER 3. CYBER-PHYSICAL ATTACKS50

The defensive strategy that we develop, as a specific implementation of the random
injection defense Procedure 3, assumes that there is a known set T of generator buses
that are known to be “trusted”, that is to say, we can assume that data from buses in
T is known to be unmodified. This concept is not new; see [25, 33, 40, 52] for related
discussions. Without such an assumption the entire suite of signals received by thecontrol centers could be falsified and it is questionable whether any meaningful attack
reconstruction can be performed. The following template describes the strategy:
Procedure 4. Pairs-Driven version of Procedure 3.
At each execution of step D1, select a random pair of generator buses s and t,
both in T , as well as random Γ > 0, and use the injections
δs = Γ, δt = −Γ, and δk = 0 ∀k 6= s, t.

(3.11)In an application of this defense, defineP̂ g = P g + δ
We denote by B̂ the bus susceptance matrix of the network, after the attack. Thismatrix will be different from the original bus susceptance matrix B in case of a topology
or susceptance attack; thus the control center does not know B̂. Recall that as statedabove we are assuming that the network remains connected after the attack.
Lemma 5. Suppose B̂θ = P g − P d , and B̂ θ̂ = P̂ g − P d . Let k 6= t be a bus such thatthe post-attack network contains a path between s and k that does not include t. Thenθ̂k − θ̂t > θk − θt .

(3.12)

Proof. Equation (3.12) does not change if we subtract from every θ̂h any constant, andlikewise with the θh . Thus, without loss of generality θ̂t = θt = 0. Under this assumption(3.12) reads:θ̂k − θk > 0.
Let M be the set of buses p 6= t such that

(3.13)

CHAPTER 3. CYBER-PHYSICAL ATTACKS

51(1) The network contains a path from s to p that avoids t, and(2) Subject to (1), θ̂p − θp is minimum.
Aiming for a contradiction, we will assume thatθ̂p − θp ≤ 0 for p ∈ M.

(3.14)Showing that (3.14) is false yields (3.13). For any line km define the flow value fkm =(θ̂k − θ̂m − θk + θm )/xkm . Since B̂(θ̂ − θ) = P̂ g − P g , the flow vector f corresponds(under the DC power flow model) to a power flow with Γ units of generation at s, Γunits of load at t, and zero generation and load elsewhere. Note that for any line km,
fkm > 0 if and only ifθ̂k − θk > θ̂m − θm .

(3.15)

θ̂s − θs > 0.

(3.16)

This observation implies[To obtain this fact, decompose the flow vector f into a set of path flows from s to t
and telescope (3.15) along any such path.] Pick any p ∈ M and let P be a path froms to p that avoids t. Say P = v0 , v1 , . . . , vi where v0 = s and vi = p, and let h be
smallest such that vh ∈ M . By (3.16) s ∈/ M , i.e., h > 0. Then by definition of h,θ̂vh−1 − θvh−1 > θ̂vh − θvh , i.e. fvh−1 ,vh > 0. But by assumption vh 6= t. So there exists
some line vh , m such that fvh ,m > 0. Therefore using the assumption θ̂k − θk ≤ 0 for allk ∈ M , vh ∈ M , and (3.15),
0 ≥ θ̂vh − θvh > θ̂m − θm .

(3.17)

So m 6= t, and as a result by construction there is a path from s to m that avoids t. But
then (3.17) contradicts the fact that vh ∈ M .
Lemma 6. Suppose k is any bus and that there are at least two generators available toimplement Procedure 3. Then a pair s, t satisfying the assumptions of Lemma 5 exists.

CHAPTER 3. CYBER-PHYSICAL ATTACKS

52

Proof. Choose s ∈ R such that s is closest to k.
Note: if Γ is chosen negative in (3.11), then instead of (3.12) we obtain θ̂k − θ̂t < θk − θtthrough essentially the same proof.
For future reference, we state the following analogue of Lemma 5, with similar proof(omitted).
Lemma 7. Let suppose θ and θ̂ be as in Lemma 5. Let k 6= t be a bus such that in the
post-attack network every path between s and k must include t. Thenθ̂k − θ̂t = θk − θt .

(3.18)To analyze the pairs-driven defense we will express phase angles using t as the reference bus; see equation (2.20). Thus B̂ θ̂ = P̂ g − P d has a unique solution with θ̂t = 0 of
the formθ̂ = B̆t (P̂ g − P d ) = θ + B̆t δ(3.19)

where B̆t is an appropriate pseudo-inverse of B̂. (In this expression θ is also written withrespect to t). As a result of Lemma 5 we have:
Lemma 8. Let k be any bus. Then, with high probability we will have θ̂k = θk + βk Γ forsome value βk > 0.
Proof. Let k be any bus. Then with high probability, at multiple iterations of the pairsdriven defense the buses k, s and t will satisfy the conditions of Lemma 5. Let us assume
that bus t is the reference bus. Then θ̂t = θt = 0, and thus θ̂k > θk . Thus, (by (3.19))
θ̂k is as desired.
This results suggests the following detection paradigm:

CHAPTER 3. CYBER-PHYSICAL ATTACKS

53

Procedure 5. Pairs-driven detection criterion
As Procedure 4 iterates, for each bus k, estimate the correlation between θ̂kand Γ. The defensive procedure terminates when all these estimates are stable.
At that point, any bus k whose correlation coefficient is nonpositive is flagged
as suspicious.

Lemma 9. With high probability the pairs-driven defense will defeat the noisy data anddata replay attacks in the sense that each bus whose data is modified will be flagged, andany bus that is not attacked will not be flagged.
Proof. Let k be any bus. Then with high probability by Lemma 8 the control center
expects that θ̂k = θk + βk Γ for some value βk > 0 unknown to the control center. This
fact yields the desired result since in either attack case, if k ∈ A the signals produced by
the attacker at bus k will not have the stated form.
Lemmas 5 through 9 assume the DC power flows model, which is only a first-orderapproximation to the AC model we consider here. In the next section we perform numerical experiments, under the AC power flows model, of the random injection defense 3.
A separate issue concerns ambient noise; we need voltage changes to overcome currently
found noise levels in measurements. This issue is taken up in Section 3.5.2.
3.5.1.1

Numerical experiments using AC power flows

The above discussion concerns DC power flows. In order to investigate how voltages
change under injection changes, under AC power flows, we perform a experiments using
examples from the Matpower library [101]. For each system we perform ten experiments.
In each experiment we compute an AC power flow which is constrained to satisfying
the given voltage bounds at all generator buses, but not at load buses, as well as power
injection constraints and generator limits, while allowing large injection changes in arandom subset of generators. For a non-generator bus k, let Vkb be its voltage in theCHAPTER 3. CYBER-PHYSICAL ATTACKS

54

base case (i.e. the Matpower case), and let Vi,k be its voltage in experiment i = 1, . . . , 10.
Finally, define|Vi,k − Vkb |.
score(k) = max.
1≤i≤10|Vkb |In Table 3.7, “Min Score” is the minimum score across all non-generator buses. Thus
the table provides experimental verification for substantial AC voltage changes under
random generator injections.
Table 3.7: AC Voltage Changes

Case

3.5.2

Min Score

Average Score

case118

11.61%

32.77%

case1354pegase

7.62%

51.00%

case2746wp

5.00%

10.09%

Overcoming Sensor Error, andthe Current-Voltage Defense

If sensor misestimation (i.e., error ) is present, a strategy based on Procedure 3 may failto detect data inconsistencies if the random power injections cause voltage changes thatare too small as compared to the error. In order to derive a version of Procedure 3 that
deals with this issue, we next describe a particular implementation of step D2 which
relies on the current-voltage consistency condition (3.1) which takes into account thepossibility of sensor error. This implementation will take into account the possibility of
sensor error. Whereas above a phasor (voltage or current quantity) φ had a true value φT(the physical value) and a reported value φR (the value received by the control center),
now we will have the sensed value φS which is the value actually produced by the sensor.
.
Due to sensor error, sensed and true data may differ. For a phasor φ define err(φ) =φS − φT . In the PMU setting, the TVE (total vector error) criterion [70, 84] guarantees

CHAPTER 3. CYBER-PHYSICAL ATTACKS

55

that
|err(φ)| < τ |φT |,(3.20)

where 0 < τ < 1 is a tolerance. Standards enforce τ = 1%, though experimental testingof PMUs shows far smaller errors [36]. From (3.20) we have that
|φS | − |φT | ≤ |φS − φT | < τ |φT |and|φT | − |φS | ≤ |φS − φT | < τ |φT |,therefore, we obtain(1 − τ )|φT | < |φS | < (1 + τ )|φT |τ
|err(φ)| <|φS |.
1−τ(3.21a)(3.21b)

We will describe three sensor-error-aware voltage-current consistency criteria. An important point is that the current-voltage consistency condition (3.1), combined withestimations of possible sensor error, yields a nonlinear relationship, and appropriately
reformulation of this relationship can render useful benefits. To simplify notation we will
drop the “(t)” from phasors though it should be understood throughout. For a line kmwriteYkm = ykk

ykm

ymk ymm.

.
T
Criterion 1. We have that Imk= ymk VkT + ymm VmT . Write zmk = [ymk ]−1 . HenceS− ymm VmS ) = err(Vk ) − zmk (err(Imk ) − ymm err(Vm ))VkS − zmk (Imk

(3.22)

which yields, using (3.20), (3.21), and the triangle inequality
S|VkS − zmk (Imk− ymm VmS )| ≤ |err(Vk )| + |zmk |(|err(Imk )| + |ymm ||err(Vm )|)τ |zmk | S(|I | + |ymm | |VmS |)1 − τ mkτ |zmk | ST= τ |zmk (Imk− ymm VmT )| +(|Imk | + |ymm | |VmS |)1−τ
2τ |zmk | S≤(|Imk | + |ymm | |VmS |).
1−τ< τ |VkT | +CHAPTER 3. CYBER-PHYSICAL ATTACKS

56

S
ove lines km must
In summary, Criterion 1 states that the sensed values VkS , VmS and Imksatisfy the following inequality:S− ymm VmS )| <|VkS − zmk (Imk2τ |zmk | S(|Imk | + |ymm | |VmS |).
1−τ(3.23)

R
Under Criterion 1, if, statistically, the reported phasors VkR , VmR , Imk
fail to satisfy (3.23)
Rline km is flagged as suspicious. A similar analysis concerns VkR , VmR , Ikm. Remark: By

construction, if k, m ∈/ A then line km will not be flagged.
Criterion 2. Proceeding as above we haveS
|Ikm
− ykk VkS − ykm VmS | = |err(Ikm ) − ykk err(Vk ) − ykm err(Vm )|τ<(|I S | + |ykk ||VkS | + |ykm ||VmS |).
1 − τ km(3.24)(and similarly with Imk ), andCriterion 3. When line km is a pure impedance line (no transformer), from (2.7) wehave that
Ykm = 

ykk
ymk

sh
ykm
y
+
ykm
km
2
 = 

ymm
−ykm






−ykm


sh  ,
ykm
ykm +
2

sh
with branch admittance ykm and shunt admittance ykm
S
S
S
T
S
T
|Ikm
+ Imk
| = |(Ikm− Ikm) + (Imk− Imk) + (ykk + ymk )VkT + (ykm + ymm )VmT |= |err(Ikm ) + err(Imk ) + (ykk + ymk )VkT + (ykm + ymm )VmT |shykm(VkT + VmT )|2shτ |ykm
|
τ
S
S
|) +
<
(|Ikm
| + |Imk(|V S | + |VmS |).
1−τ
2(1 − τ ) k= |err(Ikm ) + err(Imk ) +(3.25)

If the reported phasors do not satisfy (3.24) or (3.25) then the line is flagged.
Discussion
Note that a line not attacked will not be flagged, as per the TVE condition. Additional
criteria can be developed to handle power-injection consistency. To analyze the effectiveness of these criteria, we turn to the enhanced noisy data attack discussed in Section 3.4.

CHAPTER 3. CYBER-PHYSICAL ATTACKS57

To remind the reader, in this type of attack the voltage readings in ∂A can be arbitrarily
adjusted. While this action may create inconsistencies on lines with just one end in ∂A
the attacker may be able to “hide” such inconsistencies if they are small enough relative
to sensor error.
We next show that Criterion 1 alone can suffice to defeat the enhanced noisy dataattack (i.e. uncover inconsistencies) when voltage angles are sufficiently changed under
our random injection defense.
To understand this point, consider a bus k ∈ ∂A such that there is a line km withm∈
/ A and also a line ka where a ∈ A−∂A. We study an iteration of the random defense
which (to simplify notation) we assume begins at time t = 0. Consider line ka first. To
avoid having line ak flagged, the attacker t will need to manufacture a time series VkR (t),R(t) that (statistically) satisfy (3.23). But under the noisy data attack, on
VaR (t) and IakRRaverage VaR (t) = VaR (0) and Iak(t) = Iak(0). Hence the attacker needs (on average) that2τ |zak | RR(|Iak (0)| + |yaa | |VaR (0)|) > |VkR (t) − zak (Iak(0) − yaa VaR (0))|
1−τ= |VkR (t) − VkR (0)|.
RS
Now consider line km. Since m ∈/ A, VmR (t) = VmS (t) and Imk(t) = Imk(t). Also, denote:• VkT (∗) = the true voltage at k at the start of the current iteration of the random
defense, i.e. the voltage resulting from the injection changes in step D1. Then,assuming unbiased sensor errors and zero-mean ambient noise, VmT (∗) will equal
the expectation of VmT (t) during the iteration.
T• Likewise define the current Imk(∗).

Hence the attacker needs (on average) that2τ |zmk | TT(∗) − ymm VmT (∗))|(|Imk (∗)| + |ymm | |VmT (∗)|) > |VkR (t) − zmk (Imk1−τ
= |VkR (t) − VkT (∗)|.
As a result of these observations we have:CHAPTER 3. CYBER-PHYSICAL ATTACKS

58

Lemma 10. Consider buses k, a, m as described above. Suppose that
|VkT (∗) − VkR (0)|>2τ |zak | R
2τ |zmk | T(|Iak (0)| + |yaa | |VaR (0)|) +(|Imk (∗)| + |ymm | |VmT (∗)|)
1−τ
1−τ

(3.26)

Then it is impossible for the enhanced noisy data attacker to statistically satisfy Criterion
1 on both lines ka and km.
Comment: This lemma highlights how large changes in voltages caused by the random
defense challenge the attacker.
Experiment
Next we describe a set of experiments involving the current-voltage defense applied tothe attack given in Section 3.3.4. The current defense was implemented as follows:• For any generator bus k ∈/ R, |δk | ≤ Pkg . We used values  = 0.01, 0.05.
• The set of responding generators, R, was of cardinality 200. For k ∈ R |δk | can bearbitrarily large. We chose δk > 0 with probability 1/2.
• No generator may exceed its limits (voltage or generation), but subject to all theseP
conditions we maximize k∈G |δk |.
In Table 3.8, we perform the above analysis on the lines (k = 1139, a = 1137) and(k = 1139, m = 1110) with τ = 0.01. “Ratio” is the ratio of the left-hand side to theright-hand side of expression (3.26). We see that the condition for Lemma 10 is amply
satisfied. A similar analysis pertains to line (1141, 1361), the other line connecting A to
its complement.

CHAPTER 3. CYBER-PHYSICAL ATTACKS

59

Table 3.8: Current-voltage defense.

Experiment 1
0.01
289.01
174.47P
δ+Pk∈G k−k∈G δkExperiment 2
0.05
964.77256.04

Line (k = 1139, a = 1137)|VaR (0)|∠θaR (0)
R
Iak
(0)

1.0919∠ − 6.993◦−0.0275 + 0.0281j

1.0919∠ − 6.993◦−0.0275 + 0.0281j

T|VmT (∗)|∠θm
(∗)
T
Imk
(∗)

Line (k = 1139, m = 1110)
1.0309∠ − 7.822◦
1.0391∠ − 7.848◦
0.0905 − 0.4976j
0.1289 − 0.4901j|VkR (0)|∠θkR (0)|VkT (∗)|∠θkT (∗)

Voltages at k = 1139
1.0919∠ − 6.991◦
1.0919∠ − 6.991◦
1.0104∠ − 7.822◦
1.0187∠ − 7.936◦

Lemma 10 applied to bus k = 1139
Ratio
1.913

3.5.3

1.732

Covariance Defense

In this section we describe an elaboration of the random-pairs defense Procedure 4; theelaboration is motivated by the fact that real-world PMU data streams exhibit nongeneric stochastic structure in (for example) voltage angles [19, 21, 93]. In particular,
covariance matrices across several time scales have very low rank (typically smaller than
10). Our defense will defeat both the noisy-data and data-replay attacks, under appropriate assumptions.
As before, we assume that the buses in a certain set T are trusted. The emphasisof the methods in this section is that we aim to modify the covariance matrix of phaseangles, whereas the random injection defense in Procedures 3 or 4 change the averagevoltage values. Such a change should prove more difficult for the attacker to correctly
counteract since such a correction involves an estimation that requires time, during which

CHAPTER 3. CYBER-PHYSICAL ATTACKS

60

the attacker will be producing incorrect data.
We additionally assume that the attacker’s data stochastics are stationary (i.e. theparameters of the stochastic process do not change as a function of time). This implies
in particular that the attacker does not react to the covariance defense by changing thestochastics. Below we will discuss, however, why reacting to the defense would provevery difficult. We also assume that ambient conditions are also stationary.
In order to describe the defense we need some definitions. For a pair of buses s, t ∈ T ,define the vectors us,t and v s,t by



1,
if k = s,



.
us,t
= −1, if k = t,
k





0,
otherwise(3.27a).
v s,t = B̆t us,t .

(3.27b)Formally, the covariance defense works as follows. Let t1 , t2 be two fixed trusted
buses, and let P be a real-valued probability distribution, with zero-mean and varianceσΓ2 > 0. The defense has two phases.
(I) During an initial phase, after suspecting the attack, for i = 1, 2, we compute thematrix.
σθ2R,i = covariance matrix of observed phase angles, with respect to reference bus ti .
(II) After the initial phase, we perform iterations as in Procedure 4, as follows. We
randomly (uniformly) choose one pair of buses of the form (s, t) where s ∈ T and t =t1 or t2 , and we draw a random value Γ from the distribution P, independently fromthe stochastics of the attacker and ambient stochastics. We then apply step D1 ofProcedure 4 using this triple (s, t, Γ).
Throughout the second phase, we compute the following matrix, for i = 1, 2,.
• σθ̂2R,i = covariance matrix of observed phase angles, with respect to reference bus ti .

CHAPTER 3. CYBER-PHYSICAL ATTACKS61

The defense concludes when the estimates for these two matrices become stable.
Procedure 6. Covariance-driven detection criterion
At termination of the defense, we flag a bus k as suspicious if, for both i = 1
and 2, the difference between the (k, k) entry of σθ̂2R,i and the corresponding
entry of σθ2R,i is smaller than λ, defined byσΓ2ω, where|T | − 1(3.28a).
ω = min{(vjs,t )2 : vjs,t 6= 0}.

(3.28b)

.
λ =

s,t,j

This concludes the description of the defense, with analysis given in Lemmas 11, 13 and 9.
In preparation for those results, suppose that at some point in phase (II) the pair(s, ti ) has been selected. Let(a.1) θ̂ T,i be the vector of true voltage phase angles, using ti as the reference bus.
(a.2) θ T,i describe the true vector of voltage phase angles, had the power injections in
the defense not been applied at that point of time. It is also given using ti as thereference bus.
Lemma 11 given next concerns the relationship between these last two vectors.
Lemma 11. Suppose that at some point in (II) the pair (s, ti ) is being used. Then,θ̂ T,i = θ T,i + Γv s,ti .
Proof. Given that the pair (s, ti ) is being used, the injections in the random defense, perequation (3.11) are given by δ = Γus,ti . The result follows from equation (3.19).
The next result presents a key feature of the covariance of phase angles. Recall that
Γ is drawn independent of all other stochastics.
Lemma 12. For i = 1, 2,σθ̂2T,i = σθ2T,i +XσΓ2v s,ti (v s,ti )> .
|T | − 1 s∈T −ti(3.29)

CHAPTER 3. CYBER-PHYSICAL ATTACKS

62

Proof. We proceed by conditioning on the pair (s, ti ) being selected by the defense.
Subject to this conditioning, by Lemma 11 the covariance of θ̂ T,i equalsσθ2T,i + σΓ2 v s,ti (v s,ti )> + covar(θ T,i , Γv s,ti ).
The last term in this expression is zero, by the independence assumption on Γ. The
result follows since each pair is chosen with probability (|T | − 1)−1 .
Lemma 13. Let k be any bus. Then, for at least one of i = 1 or 2, the (k, k) entry ofσθ̂2T,i is at least as large as the corresponding entry of σθ2T,i , plus λ —defined as in (3.28a).
Proof. Without loss of generality, there is a path between k and t1 that avoids t2 . Notethat the pair (s, t2 ) with s = t1 is one of the pairs available for the defense. By Lemma 5,
we have that vkt1 ,t2 > 0 and so vkt1 ,t2 ≥ ω 1/2 . Considering (3.29) for i = 2 we see that oneof the terms in the sum corresponds to s = t1 . As just argued, the (k, k) entry of thisterm is at least ω. The (k, k) entries in the remaining terms of the sum are nonnegative(since each term is a positive-semidefinite matrix). Thus the result follows.
Lemma 14. The suspicious labels computed by the covariance defense are correct.
Proof. Consider first a bus k that is not attacked. For such a bus, by definition, θ̂kT,i =θ̂kR,i , for both i = 1, 2. Thus, by Lemma 13, bus k is not flagged as suspicious. On
the other hand, suppose k is attacked. Then under either the noisy-data or data-replay
attacks the (k, k) entry of σθ̂2R,i will be equal to the corresponding entry of σθ2R,i , by thestationarity assumption (the attacker does not change stochastics when the defense is
implemented). Hence bus k is flagged.
Remarks:(1) Recall (3.28a) and (3.28b). The quantity ω depends on the bus susceptance matrix
B, only. Hence by choosing σΓ2 large enough we can make λ large. Also note that in theabove proofs, we can restrict the set T to a subset of size 2, again helping λ attain large
values.

CHAPTER 3. CYBER-PHYSICAL ATTACKS

63(2) Given a pair (s, ti ) used in the defense, by Lemma 5 any entry vks,ti is positive if there
is a path from bus k to s that avoids ti . Let A denote the set of such buses. By Lemma 7,
for k ∈/ A, vks,ti = 0.
Thus, in the term v s,ti (v s,ti )> in (3.29) the entire submatrix with rows and columns
in A is positive, and the remaining entries in v s,ti (v s,ti )> are zero. By adjusting the proof
of Lemma 13 we conclude that for k and m in A, the entry (k, m) of σθ̂2T,i is at least
as large as the corresponding entry of σθ2T,i , plus λ. Thus a submatrix of the covariancematrix will change via the defense (and not just the diagonal entries). If the network
is guaranteed to be 2-connected [5] one can in fact prove that the entire matrix mustchange.
The covariance defense has an additional important feature, namely that the sum on
the right-hand side of (3.29) has rank |T | − 1 (as shown next) whereas we expect theleft-hand side of (3.29) to have low rank.
Lemma 15. For i = 1, 2 the vectors v s,ti = B̆ti us,ti as in (3.27b) are linearly independent.
Hence the second term in (3.29) has rank at least |T | − 1.
Proof. The |T | − 1 vectors us,ti arising from all pairs (s, ti ) under consideration arelinearly independent, by construction in (3.27a). Hence, the corresponding vectors v s,ti
are also linearly independent.
Lemma 15 highlights the challenges faced by the attacker, even if the attacker is awarethat the covariance defense is being deployed. The attacker will have to alter reported
data in a way consistent with an appropriate rank change, but the attacker does notknow the pairs (si , t) being used (or the distribution P). Such “learning” would require
data observations, i.e. time, during which the attacker is still expected to produce datareadings, producing an error trail.

CHAPTER 4. LEARNING FROM PMU DATA

64

Chapter 4
Learning from PMU data
This chapter describes a statistical study that we have performed using a dataset of PMU
measurements that we have obtained from an industrial partner. We detail the natureof the data together with the statistical tools that we have used. Plots and videos help
visualize properties of the sampling.
This study was mainly performed under the supervision of Prof. Michael Chertkov
while I was visiting Los Alamos National Lab, NM, in September 2018. Dr. Jonatan
Ostrometzky has helped to improve the Fourier filters that we have applied.

4.1

Introduction

Data-driven techniques in power systems have at least fifty years of history, starting withstatic state estimations developed by Schweppe and co-authors [77–79], then transitioning
to dynamic state estimation analysis and applications, see e.g. [45, 62, 88] and references
therein, and most recently discussed under the umbrella of “big data” as the most significant enabler of power system operations, security and resiliency in the future [9,49]. (Seealso related discussion in the description of the US Department of Energy new funding
opportunity to “explore the use of big data, artificial intelligence, and machine learning technology to leverage the power of grid sensors” [4].) Many specific questions andCHAPTER 4. LEARNING FROM PMU DATA

65

approaches, including but not limited to detection of modes of oscillation and related
analysis of stability [1–3, 88], dimension reduction for faster processing and analysis [31],early event detection [94], missing data recovery [38], identification of cyber attacks [37]and real time (online) event detection [57] are among the most recent research thrusts.
On the methodology side data-driven methods developed in other engineering disciplines have been adopted, modified and used for many (e.g. aforementioned) powersystem applications. Principal component analysis [31, 37, 38, 57, 94], auto-correlation
analysis of memory effects [74], and linear model driven spectral analysis of the dynamic
state matrix [1–3, 27, 62, 65, 75, 88] are arguably the most popular data-driven techniquescurrently in use in the power system research.
Even though the sophistication level of the methods already used in power systemapplications is impressive, coherence and understanding of the potential of new generationof big data methods, driven during the last decade largely through heavy investment of
IT industry, is still lacking [4]. We anticipate that many of the most modern methods,
especially Deep Learning (DL) and related techniques linked to Machine Learning (ML)and Artificial Intelligence (AI), revolutionary advances in data science and more generally
theoretical engineering [8, 12, 42, 56, 76], will impact the power-system operation-room
reality in an even more significant ways. However, one problem with applications of thenovel methods of DL and alike in sciences and engineering is that they are application
agnostic/generic – very effective for many business cases, but lacking “explainability”,
i.e. intuitive physical/engineering explanations. This significant handicap of the mostadvanced and recent ML & AI methods slows down development of related applicationsin power systems. Indeed, power system practitioners would generally not consider aspractical any new methods lacking “power systems informed” explanations.
This manuscript takes a step towards closing the gap between the rich variety of methods already developed and utilized in power systems and the yet to be unleashed power ofthe upcoming Big Data revolution. Specifically, we start walking towards exciting sophistication of DL slowly, from the well-established and intuitive trenches of practical system

CHAPTER 4. LEARNING FROM PMU DATA

66

engineering. We develop in this manuscript a pragmatic “phasor-detective” approach toanalysis of the streaming Phasor Measurement Unit (PMU) data which allows to extractand interpret spatial and temporal correlations in a computationally light fashion andwithout making any constraining assumptions about origin of the correlations.

Our Contribution
We analyze synchronized historical PMU data recorded at ≈ 200 most significant locations of a US Independent System Operator (ISO) over the course of two years. At each
PMU location the data includes complex current and voltage recorded with a millisecond resolution. Given geo-spatial locations of the PMU, but no information about thegrid characteristics and layout, we pose the following principal questions: What can we
possibly reconstruct from the data stream about the system current ambient behavior?
To answer the question we utilize available statistical tools. In relation to preliminary
data processing we apply to the raw signal three filtering techniques: moving average,
sliding time horizon and Fourier analysis pre-processing. This allows to provide robust
identification of the “quiet” periods and also prepare data for subsequent statistical analysis by means of Principal Component Analysis (PCA) and Auto-Correlation Analysis(ACA). We show that the two complementary tools, applied to the raw as well as topre-processed signal, allow to separate scales and also provide compressed, thus easy tovisualize, descriptors for online tracking of current state of the grid in a much broader way
than what the current Energy Management System (EMS) actually uses. PCA provides
a robust set of indicators which record slow/adiabatic changes on the scale of seconds toten of seconds and slower. We have observed that only a very few principal modes are
significant at any moment of time, even though these modes may be different for voltage
amplitude, phase difference and frequency (the three main characteristics) we track. The
results do not change when PCA is applied to the filtered signal, consistently with thefact that PCA averages over time but does not catch different-time correlations. ACA isthe tool used to analyze the latter, in particular identifying significant, persistent corre-

CHAPTER 4. LEARNING FROM PMU DATA

67

lations, missed by PCA, at shorter time scales - subseconds-to-seconds. Following ACA
curves at different spatial locations we were able to identify nodes where correlations do
not decay with time showing significant memory-effects. Remarkably, these nodes withsignificant memory cluster geographically. We observe two areas in the grid which show
especially strong sustainable temporal correlations. We then proceed with ACA analysisof the Fourier-filtered signal. This helps us to identify and localize different harmonics. In particular, we observe (for a particular quiet period) emergence of significant
oscillations in the 4-6 Hz range at a small number of nodes. Interestingly, nodes withsignificant sustainable oscillations are either wind farms, big aggregated loads or mid-size
generators. We conjecture that the sustainable oscillations are indicators of malfunction
at these critical elements of the grid. We also observe that sustainable oscillations, seenclearly through emergence of a residue in the ACA analysis of the raw signal, disappear
when applied to the Fourier-filtered signal (cutting off the 4-6 Hz oscillations). Finally,
analyzing spatial cross-correlations (of the residue) we were able to identify a group of
nodes with significant inter-dependency.

4.2

Logic and Main StepsIn this part of the work we study the data streams reported by over 200 PMUs operated
by an ISO and spanning a period over one year long. Figure 4.1 displays a rendition of thelocations of the PMUs using anonymized coordinates. As mentioned in Section 2.4, thedata stream from each PMU includes frequency, (complex) current and voltage, reported
30 times per second. Using this data one can obtain real-time estimates of complex powerat each location. Working with a data set this large (on the order of 28 TB) presents
some obvious challenges; additionally there are specific artifacts that can arise in thedata. For example, not all PMUs are always reporting, and occasionally some PMUs
exhibit what appears to be errant behavior.
Our work has centered on performing statistical analysis aimed at inferring “struc-

CHAPTER 4. LEARNING FROM PMU DATA

68

Figure 4.1: Geographical location of PMU’s (anonymized coordinates), each mark represents
the position of a sensor.

ture” in the underlying transmission system as well as identifying complex behaviors,such as resonance and oscillations. In this manuscript we focus specifically on identification and characterization of “quiet” periods, also known as ambient conditions periods.
Oversimplifying (see related discussion below) we consider a period quiet if fluctuations
around the mean (e.g. characterized in terms of the standard deviations) are smaller thana reasonable pre-defined threshold. This focus on the quiet/ambient periods is motivatedby the following considerations:• The development of a strong understanding of quiet periods and (in particular)
efficient online algorithms for recognition of such periods is a necessary step prior
to studying less-quiet or even anomalous regimes, for otherwise we risk significant
misinterpretation, i.e. errors in online detection of anomalies.
• As will be seen in this paper, the quiet regimes display informative patterns andcorrelations, all (slowly) time-evolving. Identifying such features is important withregards to:– Developing fast and reliable identification techniques.
– Uncovering hidden malfunction of assets thus providing significant contribu-

CHAPTER 4. LEARNING FROM PMU DATA

69

tion towards forecasting most probable (and destined to occur) failures.
• The richness of correlations observed in the quiet regime, in fact, suggests that
separation of what is normal/quiet from what is anomalous/atypical will be challenging. Even though we observe that quiet regimes dominate, relatively abrupt
jumps of moderate size (i.e. jumps exceeding tracked standard deviation by factorof two or three) are rather frequent although overall they account for a relatively
short fraction of the stream. As a result, it is rather difficult to find sufficiently
long entirely quiet periods in the available data.
• Clearly, understanding quiet vs. volatile behavior will be helpful toward building
predictive models for better optimization, control and planning.
• “Cyber-physical” attacks on power systems —like the ones described in Section 3—
are a venue where fast and effective learning of (changing) stochastics may prove
useful in identifying attacks.
The methodology adopted in this manuscript to identify the quiet periods is explainedin Section 4.3.3.

4.3

Description and Averaging of the Time Series

The available data encompasses the period from January 1st 2013 to March 21st 2014.
Each of the N = 240 PMUs records the following measurements 30 times per second:
time of the measurement (GPS tagged), bus ID, voltage amplitude, voltage phase angle,
current magnitude, current angle, and frequency. Additionally, the 2-dimensional coordinates of the PMU locations is also available, together with their corresponding nominal
voltages. We note that PMUs report 30 times/second, but they sample at a far higher
rate and perform filtering (e.g. anti-aliasing1 ) before reporting.
1

Aliasing is an effect that causes the measured/sampled signal to be different than the real one. Thismight be caused for example by precision inaccuracy of the measurement instrument.

CHAPTER 4. LEARNING FROM PMU DATA70

We denote a generic scalar or complex measurement (e.g. complex voltage) at PMU
location k, at time t by mk (t). The parameter t will be used to refer to the discrete
time sequence with each temporal data point separated from preceding one by the same
duration ∆ (1/30th of a second in our case).
Typical pre-processing steps in the statistical analysis of data (especially with thegoal of analyzing correlations) involve modifications through de-trending, offsetting (subtraction) of moving average, and normalizations. We will apply such techniques below.
Specific details are provided next.

4.3.1

Moving Average and Covariance

The moving average µ(m) and moving variance σ (m) vectors of a sampled series m iscomputed in the following way,∀k ∈ {1, . . . , N }, ∀t ≥ 1 :
(m)

(m)

µk (t; α) = α · mk (t) + (1 − α) · µk (t − 1; α),
(m)

(4.1)

2

(m)

(m)

σk (t; α) = α · mk (t) − µk (t − 1; α) + (1 − α) · σk (t − 1; α), (4.2)
(m)

(m)

with some initial values, say µk (0; α) = 0 and σk (0; α) = 1. The parameter α ∈ (0, 1)
represents the degree of weighting decrease. Note that in (4.2) we are using the moving
average defined in (4.1). One can likewise define moving covariance parameters.
Figures 4.2, 4.3 and 4.4 show the sampling of the frequency, voltage angle and voltage
magnitude (respectively) for a particular PMU during one minute; together with thecorresponding moving average and moving standard deviation for different values of α.
We introduce the zero mean and normalized zero mean data streams
(m)

(m)

(m)

(m)
mk (t) − µk (t − 1; α)m̄k (t; α) = mk (t) − µk (t − 1; α),m̂k (t; α) =

q

,

(4.3)
(4.4)

(m)
σk (t − 1; α)obtained from the input stream by making use of the moving average and variance. The
zero mean parameters will help us identify quiet periods, as discussed next.

CHAPTER 4. LEARNING FROM PMU DATA

4.3.2

71

Averaging over Sliding Time Horizon

Consider a number S of measurements, this value corresponds to the memory budget.
We define N dimensional vectors of means and variances (averaged over the (last) sliding
time horizon of duration S) as follows:∀k ∈ {1, . . . , N }, ∀t,(s)µk (t; S)t
1 X=
mk (τ ),
S τ =t−S+1

(4.5)

(s)
σk (t; S)

t
2
1 X
(s)
mk (τ ) − µk (t; S) .
=S τ =t−S+1(4.6)

Figures 4.5, 4.6 and 4.7 show the sampling of the frequency, voltage phase angle andvoltage magnitude (respectively) for a particular PMU during one minute; together withthe corresponding average and standard deviation over sliding time horizon for different
values of S.
As before, we also define the re-scaled zero-mean data vectors:
(s)

(s)

m̄k (t; S) = mk (t) − µk (t − 1; S),

(4.7)

(s)

(s)

m̂k (t; S) =mk (t) − µk (t − 1; S)q.
(s)σk (t − 1; S)(4.8)

Numerically, we found that a reasonable choice for S, that allows to separate powerelectronics (milliseconds) and electro-mechanical (seconds) time scales, is the number ofreadings in 1 second, that is S = 30. Similar normalizations are obtained when α = 0.05.

4.3.3

Quiet Periods

Given a reference time t and a length parameter Q consider the N × Q matrix of normalized measurements(m)

M (t; α; Q) = [ m̂k (τ ; α) | ∀k ∈ {1, . . . , N }, τ ∈ {t − Q + 1, . . . , t} ],corresponding to the last Q measurements before time t for all buses. We define theperiod (t − Q, t] as quiet if the absolute value of all entries of the matrix M (t; α; Q) isCHAPTER 4. LEARNING FROM PMU DATA

72

below some preset threshold. The reason behind this definition is that sudden jumps inthe data appear as large values in the normalized time series, whereas normalized values
close to zero mean that the data is behaving in a steady way. Effectively, a quiet period is
an interval of time where all sensor-reported data behave in a stationary way. Moreover
it is relatively cheap to compute M (t; α; Q) from the stream data in an online fashion,as we just have to complete a column at every time that samples are received and keep
track and update of the moving variance of each sensor.
Again, in our analysis, we used α = 0.05 as the length of the sliding time horizon,
and we consider quiet periods spanning 15 minutes. Over a selection of five differentdays across the database, we compute we compute the matrix M (t; α; Q) and record its
maximum absolute value when t spans over the complete day. Just in few cases the
maximum was below 10 units, we selected 2-3 intervals between these recorded cases.
Table 4.1 displays the specific time intervals selected as quiet periods that we haveanalyzed.
Table 4.1: Selected quiet periods.

#
1
2
3
4
5
6
7
8
9
10

Date
January 15, 2013
March 10, 2013
March 10, 2013
March 10, 2013
April 3, 2013
April 3, 2013
April 3, 2013
July 30, 2013
July 30, 2013
July 30, 2013

Time Window
12:13 – 12:28 AM
12:09 – 12:24 AM
4:29 – 4:44 AM
2:46 – 3:01 PM
1:52 – 2:07 AM
7:19 – 7:34 AM
7:26 – 7:41 PM
1:26 – 1:41 AM
4:39 – 4:54 PM
9:24 – 9:39 PM

CHAPTER 4. LEARNING FROM PMU DATA

73

Figure 4.2: Frequency measurements with respect to a reference bus (in blue) of PMU k = 156
during one minute, between 9:37 and 9:38 PM on July 30, 2013. Moving average (in red) anda band of one standard deviation (in purple) are also shown for (a) α = 0.01, (b) α = 0.05, (c)α = 0.10, and (d) α = 0.25.

CHAPTER 4. LEARNING FROM PMU DATA

74

Figure 4.3: Voltage phase angle measurements with respect to a reference bus (in blue) of PMU
k = 156 during one minute, between 9:37 and 9:38 PM on July 30, 2013. Moving average (inred) and a band of one standard deviation (in purple) are also shown for (a) α = 0.01, (b)α = 0.05, (c) α = 0.10, and (d) α = 0.25.

CHAPTER 4. LEARNING FROM PMU DATA

75

Figure 4.4: Voltage magnitude measurements with respect to bus nominal value (in blue) of
PMU k = 156 during one minute, between 9:37 and 9:38 PM on July 30, 2013. Moving average(in red) and a band of one standard deviation (in purple) are also shown for (a) α = 0.01, (b)α = 0.05, (c) α = 0.10, and (d) α = 0.25.

CHAPTER 4. LEARNING FROM PMU DATA

76

Figure 4.5: Frequency measurements with respect to a reference bus (in blue) of PMU k = 156
during one minute, between 9:37 and 9:38 PM on July 30, 2013. Average over sliding time
horizon (in red) and a band of one standard deviation (in purple) are also shown for (a) S = 10,(b) S = 30, (c) S = 300, and (d) S = 900; corresponding to the number of samples made over
1/3, 1, 10, and 30 seconds, respectively.

CHAPTER 4. LEARNING FROM PMU DATA

77

Figure 4.6: Voltage phase angle measurements with respect to a reference bus (in blue) of PMU
k = 156 during one minute, between 9:37 and 9:38 PM on July 30, 2013. Average over sliding
time horizon (in red) and a band of one standard deviation (in purple) are also shown for (a)S = 10, (b) S = 30, (c) S = 300, and (d) S = 900; corresponding to the number of samples
made over 1/3, 1, 10, and 30 seconds, respectively.

CHAPTER 4. LEARNING FROM PMU DATA

78

Figure 4.7: Voltage magnitude measurements with respect to bus nominal value (in blue) of
PMU k = 156 during one minute, between 9:37 and 9:38 PM on July 30, 2013. Average oversliding time horizon (in red) and a band of one standard deviation (in purple) are also shownfor (a) S = 10, (b) S = 30, (c) S = 300, and (d) S = 900; corresponding to the number ofsamples made over 1/3, 1, 10, and 30 seconds, respectively.

CHAPTER 4. LEARNING FROM PMU DATA

4.4

79

Fourier Filtering

In the following analysis we assume that Q is an even integer. Let F[mk (·)] ∈ CQ be thediscrete Fourier transform of mk (·), where F[mk (·)]q is the amplitude corresponding tofrequency ωq , q ∈ {0, . . . , Q − 1}. Since the frequency of the readings is 30 Hz, we canrepresent the frequency domain as Q equidistant points between 0 and 30 Hz. In otherwords,
Q−1
n
qωt o
. X
mk (t) exp −2πjF[mk (·)]q =
,
30
t=0

∀q ∈ {0, . . . , Q − 1},t.
where ωt = 30 · .
QGiven a vector x = (xq )q∈{0,...,Q−1} , the inverse Fourier transform is defined as∀t ∈ {0, . . . , Q − 1},

F

−1



Q−1
tωq
. 1 X
[x]t =
xq exp 2πj
Q q=0
30

As it is well known, we have the following
Observation 16 ( [71]). The inverse Fourier transform of the Fourier transform recovers
the original time series, that is∀t ∈ {0, . . . , Q − 1},

mk (t) = F −1 [F[mk (·)]] .

These relations give a one-to-one correspondence between the original time series(characterized in the time domain t) and its Fourier transform series (characterized in
the frequency domain ωq ). Moreover, given the periodicity of the sinusoidal functions,
we can represent the time series as


Q/2−1
tωq
1 X
shift
F
[mk (·)]q exp 2πj,
mk (t) =Q
30
q=−Q/2where the frequencies ωq spans from −15 and 15 Hz, andF[·]Q+q , if q < 0,shiftF[·]q =

F[·]q ,if q ≥ 0.

(4.9)

CHAPTER 4. LEARNING FROM PMU DATA

80

F shift represents the shifted Fourier transform, placing the zero-frequency component
to the center of the spectrum. This representation of the Fourier transform is widelyused in signal processing since the symmetry with respect to zero ease its visualizationand the application of different filters.
Observation 17. (x(t))t∈{0,...,Q−1} is a real-valued time series, if and only if, F[x]0 and
F[x]Q/2 are real and for any q ∈ {1, . . . , Q2 − 1}, F[x]q = F[x]Q−q .
Proof. (⇒) By definition of the Fourier transform, it is easy to see that F[x]0 and F[x]Q/2
are real when x is real. Let q > 0,
Q−1

n
1 X
qωt o
F[x]q =
mk (t) exp −2πj
Q t=0
30


Q−1
1 X
(Q − q)ωt=
mk (t) exp +2πj
= F[x]Q−q ,
Q t=0
30
t
t
The last equality holds, since by definition Qω
= t, and therefore exp {2πj Qω
} = 1.
30
30(⇐) Let t ∈ {0, . . . , Q − 1}, using Observation 16 and the fact that exp {2πj

tωQ−q
}=
30

q
exp {2πj t(Q−q)} = exp {−2πj tq} = exp {−2πj tω
},
Q
Q
30


tωQ/2
Q · x(t) = F[x]0 + F[x]Q/2 exp 2πj
30




Q/2−1 
X
tωQ−q
tωq
+
F[x]q exp 2πj+ F[x]Q−q exp 2πj
30
30
q=1

= F[x]0 + F[x]Q/2 exp {πtj}
Q/2−1 Xtωqtωq+
+ F[x]Q−q exp −2πj.
F[x]Q−q exp 2πj
30
30
q=1
Since the two first terms in the last expression are real, and the summation goes overterms that add a quantity and its conjugate, we conclude that Q · x(t) is real.
The property that a time series satisfies xq = xQ−q will be called conjugate-symmetry.
Corollary 18. If (x(t))t∈{0,...,Q−1} is a real-valued time series, then F shift [x]q = |F[x]q | =|F[x]Q−q | = F shift [x]−q for any q ∈ {1, . . . , Q2 − 1}.

CHAPTER 4. LEARNING FROM PMU DATA

81

Remark. An analogous result to Observation 17 can be obtain for the inverse Fouriertransform. Therefore, we conclude that a symmetric real-valued —and, therefore, conjugatesymmetric— time series has a Fourier transform that is also real-valued and (conjugate-)
symmetric.
Figure 4.8 shows the absolute value of the Fourier transform for a particular PMU
during a 15-minute quiet period. We observe predominant peaks near 5 Hz in the frequency domain of the signal, this particularity repeats over different PMUs and differentdays and hours. Figure 4.9 shows the absolute value of the Fourier transform of thecomplex voltage of the same time series as Figure 4.8.

Figure 4.8: Absolute value of the Fourier transform over the positive frequency domain (0 to
15 Hz) of (a) frequency, (b) voltage phase angle, and (c) voltage magnitude of PMU k = 156
over quiet period #10 (see Table 4.1). Mean have been subtracted from the raw signal —toprevent a high peaks at frequency 0.

Consider a vector ϕ̃ = ϕ̃q : q ∈ {− Q2 , . . . , Q2 − 1} ∈ CQ .

CHAPTER 4. LEARNING FROM PMU DATA

82

Figure 4.9: Absolute value of the Fourier transform of the complex voltage of PMU k = 156
over quiet period #10. Again, mean have been subtracted from the raw signal.

Definition 19. We denote by (fk [ϕ̃](t))t∈{0,...,Q−1} the ϕ̃-filtered time series of mk (·) thatis defined as the vector whose Fourier transform is the result of multiplying componentwise each of the terms of ϕ̃ and the shifted Fourier transform of mk (·).
In mathematical terms, the previous definition can be written as∀q ∈ {0, . . . , Q − 1},
We will use the operator

to denote the component-wise product of two vectors, there-

fore we can write F[fk [ϕ̃](·)] = ϕ̃∀t ∈ {0, . . . , Q − 1},F[fk [ϕ̃](·)]q = ϕ̃q · F shift [mk (·)]q .

F shift [mk (·)]. Equivalently, we have:



Q/2−1
tωq
. 1 X
shift
fk [ϕ̃](t) =ϕ̃q · F[mk (·)]q exp 2πj.
Q30
q=−Q/2(4.10)
If we extend the definition of the vector ϕ̃, in such a way that it repeats its values
periodically, meaning that for any integer n, ϕ̃nQ+q = ϕ̃q for all q ∈ {− Q2 , . . . , Q2 − 1}; we
can equivalently write (4.10) as


Q−1
1 X
tωq
fk [ϕ̃](t) =ϕ̃q · F[mk (·)]q exp 2πj
.
Q q=0
30
.
Let us denote ϕ = F −1 [ϕ̃] the inverse Fourier transform of ϕ̃. Then,
Lemma 20 ( [71]). ∀t ∈ {0, . . . , Q − 1}, fk [ϕ̃](t) =

Q−1
P
s=0

mk (s) · ϕ(t−s) .

(4.11)

CHAPTER 4. LEARNING FROM PMU DATA

83

Therefore, we can express the filtered series as in Definition 19 or as a convolution
between the original series mk (·) and ϕ.
Analogously to equations (4.5)–(4.8), we can compute (averaging over a sliding time
horizon) the zero-mean data vector and the re-scaled zero-mean data vector for fk [ϕ̃](·)(s)(s)
instead of mk (·), we will call them f¯k [ϕ̃](·; S) and fˆk [ϕ̃](·; S), respectively, when the

average is made over the last S measurements.
We will define next several filter vectors and compare the time series that result after
applying the corresponding Fourier filters.

4.4.1

Cutting off high frequency components

We want to understand and visualize how the filtered data looks like when high frequencies are suppressed. For that purpose, we will consider the following family of cut-off(CO) filter vectors:∀q ∈ {− Q2 , . . . , Q2 − 1},
1, if |ωq | ≤ λ,
.
COλ
ϕ̃q
= (χ[−λ,λ] )q =0, otherwise,(4.12)

where λ ∈ (0, 15) is the cut-off parameter.
Figure 4.10 shows the filter vector ϕ̃COλ and its inverse Fourier transform ϕCOλ ,for λ = 4 Hz. Figure 4.11 shows (on plots on the left) the original time series forfrequency, voltage phase angle and voltage magnitude of PMU k = 156 across 10 secondsof measurements (300 samples). The plots on the right show the resulting filtered timeseries when the filter vector ϕ̃CO4 is applied. As a remainder, the filter has been applied
to the frequency domain of the Fourier filter obtained on an interval of 15 minutes.
We also show that when we cut off frequencies higher than the peaks around 5 Hzthat we see in Figures 4.8 and 4.9, for example, cutting off frequencies larger that 6 Hz,
we obtain almost no change between the original and the filtered time series. Figure 4.12
show the cut-off filter vector ϕ̃CO6 and the filtering results are shown in Figure 4.13.

CHAPTER 4. LEARNING FROM PMU DATA

84

Figure 4.10: (a) Absolute value and (c) real and imaginary parts of filter vector ϕ̃CO4 ; (b) and(d) show the component-wise absolute value of the inverse Fourier transform of ϕ̃CO4 .

Figure 4.11: Original sampling (in blue) for (a) frequency, (c) voltage phase angle, and (e)voltage magnitude across 10 seconds on period #10 (see Table 4.1). Corresponding filteredseries are shown (in purple) in plots (b), (d), and (f) using filtering vector ϕ̃CO4 .

CHAPTER 4. LEARNING FROM PMU DATA

85

Figure 4.12: (a) Absolute value and (c) real and imaginary parts of filter vector ϕ̃CO6 ; (b) and(d) show the component-wise absolute value of the inverse Fourier transform of ϕ̃CO6 .

Figure 4.13: Original sampling (in blue) for (a) frequency, (c) voltage phase angle, and (e)voltage magnitude across 10 seconds on period #10 (see Table 4.1). Corresponding filteredseries are shown (in purple) in plots (b), (d), and (f) using filtering vector ϕ̃CO6 .

CHAPTER 4. LEARNING FROM PMU DATA

86

More sophisticated filters are used in signal processing, so that the cut-off process
can be made in a online and smooth fashion. The cut-off vector presented above is an
example of a low-pass filter (LPF), meaning that the filter passes lower frequencies thana selected cut-off frequency.
In the following examples we present LPF vectors that are first described by its inverse
Fourier transform ϕLPFλ —i.e. the coefficients that will participate in the convolution to
get the filtered series— for cut-off frequency λ. In this filtering process, the convolution
is made with a small number (compared with Q) of measurements, the rest of coefficientsare set to be zero; this property allows to have an online computation of the filtered timeseries. These filters are called finite impulse response (FIR) systems [71].
On the other hand, the infinite impulse response (IIR) systems consist of infinite
vectors of non-zero coefficients that are convoluted with an infinite time series (in our case,we might think about extending our original time series by repeating it periodically out of
the [0, Q − 1] interval). For illustrative purposes, consider the cut-off filter vectors (4.12)
in the continuum set-up:∀ω ∈ [−15, 15],

Hd (ω) =
1, if |ω| ≤ λ,(4.13)

0, otherwise,
this is the ideal desired frequency response and can be represented as


∞
X
ωt
,
Hd (ω) =hd (t) exp −2πj30t=−∞(4.14)

where hd is the corresponding IIR sequence, that expressed in terms of Hd becomes


Z 15
tω
1
∀t ∈ Z,
hd (t) =Hd (ω) exp 2πjdω.
(4.15)
30 −1530
Figure 4.14 shows the desired frequency response Hd for λ = 5 Hz and its partial IIR
sequence. As expected, the IIR vector assimilates the ones computed above for ϕCOλ (seeFigures 4.12d and 4.10d).
Given the impractibality of convoluting a times series by an infinite vector, the computation is usually truncated or, equivalently, the IIR sequence is truncated to a chosenbounded window and setting equal to zero the terms that fall out of the selected window.

CHAPTER 4. LEARNING FROM PMU DATA

87

Figure 4.14: (a) Desired frequency response Hd for cut-off frequency λ = 5 Hz and (b) its
infinite impulse response (showing for t between −100 and 100).

As explained in [71], one technique of designing FIR filters is by windowing. That is,given a finite-duration window sequence w, wherew(t) = 0 for t < −M1 or t > M2 ,(with M1 , M2 ∈ N),

a FIR filter is obtained by the componen-wise product between hd and w:∀t ∈ Z,(or, equivalently, h = hd.
h(t) = hd (t)w(t),(4.16)

w). Since filtering the original time series mk (·) by the filter

vector h(·) corresponds to the convolution in the time domain of this two sequences, inorder to obtain the filtered element at some time t we need M1 future measurements.
Therefore, we can obtain an online filtered series with a delay of M1 time steps.
Different type of windows have been proposed in the literature and give appropriate
results depending on the nature and the pyshics behind of the time series. The simplest
one if the rectangular window
1, if − M1 ≤ t ≤ M2 ,.
w1 (t) =
0, otherwise,(4.17)

CHAPTER 4. LEARNING FROM PMU DATA

88

that simply truncates the IIR to the [−M1 , M2 ] window. The average over a sliding
time horizon defined in (4.5) can be computed as a filtered time series with FIR filter
h = hdw, where hd (t) = 1 for all t and w is a rectangular window with M1 = 0 and

M2 = S − 1. In the following examples, we will use the Hann window [44] defined as
 





1
2π(t + M )π(t + M )2
 1 − cos
= sin
, if − M ≤ t ≤ M,
.
2
2M
2M
w2 (t) =(4.18)
0,
otherwise,
for M ∈ N, since it has given us satisfactory results. Figure 4.15 shows the positivecoefficients of a rectangular window and a Hann window.

Figure 4.15: (a) Rectangular window and (b) Hann window.

By increasing M we obtaing a filter that is closer to the desired frequency response,
but at the same time, we increase the delay of the online filtering process. Therefore, wemust look for a small parameter M that gives useful filtered time series.
Figure 4.16 shows the low-pass filter for λ = 4 Hz described by 101 coefficients (bysetting M = 50 in (4.18)), shown as black dots in plot (d). The Fourier transform of the
FIR filter h = ϕLPF4 results as a smoother version of the one obtained by ϕCO4 (compare
Figures 4.10a and 4.16a).
Finally, we show that the amount of coefficients used to describe the LPF vector isrelated with the slope of the vector from the region of frequencies that are kept and theones that are cut; the more coefficients the steepest the slope is. Figures 4.18 to 4.29
show the LPF vector ϕLPFλ with 31 coefficients (by setting M = 15 in (4.18)) and theresulting filtered time series for λ = 1, 2, 3, 4, 5, and 6 Hz.

CHAPTER 4. LEARNING FROM PMU DATA

89

Figure 4.16: (a) Absolute value, (c) real and imaginary parts of filter vector ϕ̃LPF4 (M = 50);(b) and (d) show the component-wise absolute value of the inverse Fourier transform of ϕ̃LPF4 .

Figure 4.17: Original sampling (in blue) for (a) frequency, (c) voltage phase angle, and (e)voltage magnitude across 10 seconds on period #10 (see Table 4.1). Corresponding filteredseries are shown (in purple) in plots (b), (d), and (f) using filtering vector ϕ̃LPF4 .

CHAPTER 4. LEARNING FROM PMU DATA

90

Figure 4.18: (a) Absolute value and (c) real and imaginary parts of filter vector ϕ̃LPF1 (M = 15);(b) and (d) show the component-wise absolute value of the inverse Fourier transform of ϕ̃LPF1 .

Figure 4.19: Original sampling (in blue) for (a) frequency, (c) voltage phase angle, and (e)voltage magnitude across 10 seconds on period #10 (see Table 4.1). Corresponding filteredseries are shown (in purple) in plots (b), (d), and (f) using filtering vector ϕ̃LPF1 .

CHAPTER 4. LEARNING FROM PMU DATA

91

Figure 4.20: (a) Absolute value and (c) real and imaginary parts of filter vector ϕ̃LPF2 (M = 15);(b) and (d) show the component-wise absolute value of the inverse Fourier transform of ϕ̃LPF2 .

Figure 4.21: Original sampling (in blue) for (a) frequency, (c) voltage phase angle, and (e)voltage magnitude across 10 seconds on period #10 (see Table 4.1). Corresponding filteredseries are shown (in purple) in plots (b), (d), and (f) using filtering vector ϕ̃LPF2 .

CHAPTER 4. LEARNING FROM PMU DATA

92

Figure 4.22: (a) Absolute value and (c) real and imaginary parts of filter vector ϕ̃LPF3 (M = 15);(b) and (d) show the component-wise absolute value of the inverse Fourier transform of ϕ̃LPF3 .

Figure 4.23: Original sampling (in blue) for (a) frequency, (c) voltage phase angle, and (e)voltage magnitude across 10 seconds on period #10 (see Table 4.1). Corresponding filteredseries are shown (in purple) in plots (b), (d), and (f) using filtering vector ϕ̃LPF3 .

CHAPTER 4. LEARNING FROM PMU DATA

93

Figure 4.24: (a) Absolute value and (c) real and imaginary parts of filter vector ϕ̃LPF4 (M = 15);(b) and (d) show the component-wise absolute value of the inverse Fourier transform of ϕ̃LPF4 .

Figure 4.25: Original sampling (in blue) for (a) frequency, (c) voltage phase angle, and (e)voltage magnitude across 10 seconds on period #10 (see Table 4.1). Corresponding filteredseries are shown (in purple) in plots (b), (d), and (f) using filtering vector ϕ̃LPF4 .

CHAPTER 4. LEARNING FROM PMU DATA

94

Figure 4.26: (a) Absolute value and (c) real and imaginary parts of filter vector ϕ̃LPF5 (M = 15);(b) and (d) show the component-wise absolute value of the inverse Fourier transform of ϕ̃LPF5 .

Figure 4.27: Original sampling (in blue) for (a) frequency, (c) voltage phase angle, and (e)voltage magnitude across 10 seconds on period #10 (see Table 4.1). Corresponding filteredseries are shown (in purple) in plots (b), (d), and (f) using filtering vector ϕ̃LPF5 .

CHAPTER 4. LEARNING FROM PMU DATA

95

Figure 4.28: (a) Absolute value and (c) real and imaginary parts of filter vector ϕ̃LPF6 (M = 15);(b) and (d) show the component-wise absolute value of the inverse Fourier transform of ϕ̃LPF6 .

Figure 4.29: Original sampling (in blue) for (a) frequency, (c) voltage phase angle, and (e)voltage magnitude across 10 seconds on period #10 (see Table 4.1). Corresponding filteredseries are shown (in purple) in plots (b), (d), and (f) using filtering vector ϕ̃LPF6 .

CHAPTER 4. LEARNING FROM PMU DATA96

We observe that there is prominent difference is the filtered series when frequencies
larger than 5 Hz are erased from the original sampling. We will consider other two filter
vectors to study and visualize the effect of the peaks near 5 Hz in the frequency domainof the original signal.

4.4.2

High Frequency Filter

Contrary to LPF we can also describe filter that keep the frequencies in the high range,these filters are called high-pass filters (HPF). Figure 4.30 describe the filter ϕ̃HPF5 thatcuts frequencies lower that 5 Hz, the filter is built given the M = 15 convolution coefficients depicted in plot (d). Figure 4.31 shows the resulting filtered series, note that now
the filtered series have mean close to 0 (since the mean of the time series is represented
by the zero-frequency coefficient of the Fourier transform, which is erased). Also notethat the plots of the filtered series (on the right) have the same range length as theoriginal time series plots (on the left). The oscillation showed on the filtered series might
be obtained by the residue from the peaks around 5 Hz in the frequency domain of theoriginal sampling, since the filter vector has mild slope.

4.4.3

Band-pass and Band-stop Filters

The two last filter that we will describe are called band-pass filter (BPF) and bandstop filter (BSF). The first one keeps the frequencies of the original time series (in thefrequency domain) belonging to a particular band around a selected frequency; while thelatter does the opposite, it suppresses the frequencies in that band. We will use thesefilters to specifically observe the contribution of the frequencies around the observed
peaks at 5 Hz using ϕ̃BPF5 (see Figures 4.32 and 4.33), and we remove them with ϕ̃BSF5(see Figures 4.34 and 4.35). In both cases, the filter vector have been described using
M = 30 coefficients.

CHAPTER 4. LEARNING FROM PMU DATA

97

Figure 4.30: (a) Absolute value and (c) real and imaginary parts of filter vector ϕ̃HPF5 (M = 15);(b) and (d) show the component-wise absolute value of the inverse Fourier transform of ϕ̃HPF5 .

Figure 4.31: Original sampling (in blue) for (a) frequency, (c) voltage phase angle, and (e)voltage magnitude across 10 seconds on period #10 (see Table 4.1). Corresponding filteredseries are shown (in purple) in plots (b), (d), and (f) using filtering vector ϕ̃HPF5 .

CHAPTER 4. LEARNING FROM PMU DATA

98

Figure 4.32: (a) Absolute value and (c) real and imaginary parts of filter vector ϕ̃BPF5 (M = 30);(b) and (d) show the component-wise absolute value of the inverse Fourier transform of ϕ̃BPF5 .

Figure 4.33: Original sampling (in blue) for (a) frequency, (c) voltage phase angle, and (e)voltage magnitude across 10 seconds on period #10 (see Table 4.1). Corresponding filteredseries are shown (in purple) in plots (b), (d), and (f) using filtering vector ϕ̃BPF5 .

CHAPTER 4. LEARNING FROM PMU DATA

99

Figure 4.34: (a) Absolute value and (c) real and imaginary parts of filter vector ϕ̃BSF5 (M = 30);(b) and (d) show the component-wise absolute value of the inverse Fourier transform of ϕ̃BSF5 .

Figure 4.35: Original sampling (in blue) for (a) frequency, (c) voltage phase angle, and (e)voltage magnitude across 10 seconds on period #10 (see Table 4.1). Corresponding filteredseries are shown (in purple) in plots (b), (d), and (f) using filtering vector ϕ̃BSF5 .

CHAPTER 4. LEARNING FROM PMU DATA

4.5

100

Covariance Matrix andPrincipal Component Analysis

In this section we will describe a method to study the correlation of the measurements.
We start by defining a covariance matrix of the sampling, which then is analyzed in termsof its spectrum (eigen-values and eigen-vectors), specifically we will observe how theseparameters change over time.
Let T ≤ Q be the number of measurements of a time period. (Recall that Q is thenumber of measurements in a quiet period, that in our analyses has been set to be of
15-minute long.) The N × N covariance (correlation) matrix of the signal, which is basedon the last T measurements, is defined by:∀ time t ∈ {T, . . . , Q},.
Σ0 (t; T ; m (·)) =
∗

h 1

t
X

T τ =t−T +1m∗k (τ )m∗` (τ )∀k, ` ∈ {1, . . . , N }t
1 Xm∗ (τ )m∗ (τ )H ,=T τ =t−T +1i(4.19)

where m∗ (·) could be replaced by any of the normalized series defined above, that is,m̄(m) (·; α), m̂(m) (·; α), m̄(s) (·; S), m̂(s) (·; S), f¯(s) [ϕ̃](·; S), or fˆ(s) [ϕ̃](·; S). It is easy tocheck that Σ0 (t; T ; m∗ (·)) is a positive semidefinite matrix: let z ∈ CN , thenz H Σ0 (t; T ; m∗ (·)) z =

t
t
1 X
1 X
2
z H m∗ (τ )m∗ (τ )H z =z H m∗ (τ ) ≥ 0.
T τ =t−T +1
T τ =t−T +1

We aim to choose a (fixed) value T that is sufficiently large so that the covariance isweakly dependent on T , but also not too large to keep memory as small as possible (withan eye on streaming applications). Empirical experiments show that with T corresponding to the number of measurements in three minutes we can obtain stable covariancematrices when t varies.
In what follows we will drop the inputs T and m∗ (·) from Σ0 (t; T ; m∗ (·)) since they
become fixed in the following analysis. We perform an eigen-decomposition on the cor-

CHAPTER 4. LEARNING FROM PMU DATA

101

relation matrix Σ0 (t) and track the results as the function of t:
Σ0 (t) =NX

λp (t)ξp (t)ξp (t)H ,(4.20)

p=1

where ξp (t) and λp (t) are the orthonormal eigenvectors and corresponding eigenvalues (indecreasing order) of Σ0 (t) respectively, i.e.
∀p ∈ {1, . . . , N } : Σ0 (t)ξp (t) = λp (t)ξp (t),∀p, q ∈ {1, . . . , N } : ξp (t)H ξq (t) = δpq ,λ1 (t) ≥ λ2 (t) ≥ · · · ≥ λN (t) ≥ 0,

(4.21)(4.22)(4.23)

where δ is the Kronecker delta function (δpq = 1 if and only if p = q). Since Σ0 (t) is
a positive semidefinite matrix, (4.23) is justified. In our tests, Σ0 (t) is (numerically)
rank-deficient, thus justifying the Principal Component Analysis (PCA) approach. PCA
may be considered exact, if P < N principal components are tracked or approximated.
The results of the PCA analysis are presented in Figures 4.36–4.41 which are screenshots of the movies available in [20]. The screenshots and the movies show at time tfor each of the three variables (frequency, voltage angle, and voltage magnitude) four
indicators:• normalized vector m∗ (t) of measurements: for each sensor, we plot at its geographical position the value of the normalized data using different colors for different
values;• first 40 eigenvalues of Σ0 (t), in decreasing order;• largest 4 eigenvalues λ1 (·), λ2 (·), λ3 (·), λ4 (·) for the last minute before t;• corresponding eigenvectors ξ1 (t), ξ2 (t), ξ3 (t), ξ4 (t) at time t, component values ofthe vector are plotted geographically, in blue values close to −1 and in red values
close to 1.

CHAPTER 4. LEARNING FROM PMU DATA

102

For the numerical experiments we used traditional methods to compute the eigendecomposition of Σ0 (t), typically taking O(N 3 ) time and using O(N T ) space. However,for streaming data that falls into a subspace of the original space (like in our case), we
can consider lighter and faster algorithms to compute PCA, see [22].
From the simulations we can conclude that with few eigen-vector (say, less than 10)
we can characterize more than 80% of the spectrum of the correlation matrix. Moreover,the largest eigen-values and the corresponding eigen-vectors are stable over time during
the quiet periods. All these characteristics (correlation matrix, eigen-values and eigenvectors) show spatial relationship between the measurements (that is, the relation thatexists between different PMU sensors for a period of time). In what follows, we argue
that we can also obtain temporal correlation across sensors.

4.5.1

Singular Value Decomposition

Consider the N × T measurement matrix
M (t) = [ Mk,τ = m∗k (τ ) | ∀k ∈ {1, . . . , N }, ∀τ ∈ {t − T + 1, . . . , t} ] .
The singular value decomposition (SVD) of the measurement matrix M (t) is
M (t) = U (t)D(t)W (t)H ,(4.24)

where the “spatial” matrix U (t) is an N × N matrix the columns of which are orthogonalunit vectors of length N which are called left singular vectors of M (t), the “temporal”
matrix, W (t), is S × S whose columns are right singular vectors of M (t) and D(t) is the
N × S rectangular diagonal matrix of positive numbers. The covariance matrix, Σ0 (t),is related to the measurement matrix, M (t), according to
Σ0 (t) = M (t)M (t)H = U (t)D(t)D(t)H U (t)H ,(4.25)

where we took into account that W (t)H W (t) = I. One observes that the covariancematrix does not depend on the temporal matrix W (t). To investigate time-related correlation effects, in the next section we also study auto-correlations.

CHAPTER 4. LEARNING FROM PMU DATA

103

Figure 4.36: Left: normalized frequency m̂(s) (t; S) representing each PMU in its geographicallocation. Right: First 40 eigen-values of Σ0 (t; T ; m̂(s) (·; S)) in blue, and its spectral contribution
in red. [t=21:29:00 on July 30, 2013; S = 30; T = 5400]

Figure 4.37: Largest four eigen-values from Figure 4.36 across the last minute before t, together
with their corresponding normalized eigen-vectors at time t (positive values are shown in red,
negative values are shown in blue). [t=21:29:00 on July 30, 2013]

CHAPTER 4. LEARNING FROM PMU DATA

104

Figure 4.38: Left: normalized voltage phase angle m̂(s) (t; S) representing each PMU in its
geographical location. Right: First 40 eigen-values of Σ0 (t; T ; m̂(s) (·; S)) in blue, and its spectral
contribution in red. [t=21:29:00 on July 30, 2013; S = 30; T = 5400]

Figure 4.39: Largest four eigen-values from Figure 4.38 across the last minute before t, together
with their corresponding normalized eigen-vectors at time t (positive values are shown in red,
negative values are shown in blue). [t=21:29:00 on July 30, 2013]

CHAPTER 4. LEARNING FROM PMU DATA

105

Figure 4.40: Left: normalized voltage magnitude m̂(s) (t; S) representing each PMU in its geographical location. Right: First 40 eigen-values of Σ0 (t; T ; m̂(s) (·; S)) in blue, and its spectral
contribution in red. [t=21:29:00 on July 30, 2013; S = 30; T = 5400]

Figure 4.41: Largest four eigen-values from Figure 4.40 across the last minute before t, together
with their corresponding normalized eigen-vectors at time t (positive values are shown in red,
negative values are shown in blue). [t=21:29:00 on July 30, 2013]

CHAPTER 4. LEARNING FROM PMU DATA

4.6

106

Accounting for Temporal Correlations

In this section we describe a method to obtain temporal correlation between the PMU
sensors.
Consider the delayed covariance matrix generalizing Eq. (4.19):∀∆ ≥ 0, ∀t ≥ T + ∆,.
Σ∆ (t; T ; m (·)) =
∗

=

h 1

t
X

T τ =t−T +1m∗k (τ )m∗` (τ − ∆) ∀k, ` ∈ {1, . . . , N }t
1 Xm∗ (τ )m∗ (τ − ∆)H .
T τ =t−T +1i(4.26)We are interested to study how Σ∆ (t; T ; m∗ (·)) changes when ∆ increases and then track
evolution with t. Note that, in comparison with Σ0 , Σ∆ is not necessarily positive
semidefinite when ∆ > 0.
However, evolution of the spectrum (in the two-dimensional (∆, t) space) is challenging. Instead, we study two surrogate objects, introduced in the following two subsections,which are easier to visualize. Again, since T and m∗ (·) are fixed, we will omit them as
input of the Σ∆ function.

4.6.1

Auto-Correlation Functions

To study the correlation between measurements on a sensor and its own past we introduce
the auto-correlation functions.
Definition 21. The normalized PMU’s auto-correlation functions at different nodes aredefined as follows:∀k ∈ {1, . . . , N },Ak (∆; t) =[Σ∆ (t)]kk.
[Σ0 (t)]kk(4.27)We choose to normalize the function with respect to the term [Σ0 (t)]kk so we can havea relative comparison tool for different sensor. It might be useful to think about theseCHAPTER 4. LEARNING FROM PMU DATA

107

auto-correlation functions for a fixed t and see it as a function of ∆, that is Ak (·; t). Thisobject is of interest because of the following two reasons• Dependence of the auto-correlation function on ∆ indicates whether fluctuations
around the mean at a particular node decay or not with time. Stated differently
this analysis tests if there are significant memory effects or if memory is lost.
• It accounts for the part of the measurement matrix which is ignored in the PCA
analysis, as discussed above –that is, it accounts for the temporal matrix, W .
We show the auto-correlation function for frequency at a fixed t for four selected
PMU sensors in Figure 4.42. We obtain different patters and shapes of the resulting
auto-correlation functions: some of them have large amplitude, other ones have sinusoidal
shape.

(a)

(b)

(c)

(d)

Figure 4.42: Auto-correlation functions for frequency for PMU’s (a) k = 2, (b) k = 6, (c)k = 112, and (d) k = 139; at t=21:29:00 on July 30, 2013. The correlation matrices areconstructed with the normalized time series m̂(s) (·; S), S = 30, T = 5400.

CHAPTER 4. LEARNING FROM PMU DATA

108

We will compare these functions with the one that we obtain if we first apply a bandstop and band-pass Fourier filters at 5 Hz to the original sampling, with the objectiveof suppressing the oscillations that we observe at that frequency, and then observe theircontribution on the these auto-correlation functions.

(a)

(b)

(c)

(d)

Figure 4.43: Auto-correlation functions for frequency for PMU’s (a) k = 2, (b) k = 6, (c)k = 112, and (d) k = 139; at t=21:29:00 on July 30, 2013. The correlation matrices areconstructed with the normalized time series fˆ(s) [ϕ̃BSF5 ](·; S), S = 30, T = 5400.

Figure 4.43 shows the auto-correlation functions when a band-stop Fourier filter isapplied (see details of the filter in Figure 4.34). We note that the oscillations and theamplitude of three out of the four cases that are depicted drop down, therefore, thereason of the high auto-correlation of the measurements might be explained by the 5 Hz
oscillation of the sampling —as we have seen, shown as peaks in the frequency domainof the Fourier transform.
Figure 4.44 illustrates the auto-correlation functions when a band-pass Fourier filter
is applied, the frequencies that are kept in the spectrum are shown in Figure 4.32. We

CHAPTER 4. LEARNING FROM PMU DATA

109

observe that for the three sensors where the band-stop filter made a difference, the bandstop filter keeps the same shape of the functions as the original series (without filtering)
amplifying their amplitude; however, for sensor k = 112 on Figure 4.44c, the autocorrelation functions changes completely, implying that the cause of the auto-correlation
is out of the 4-6 Hz range in the frequency domain of the signal.

(a)

(b)

(c)

(d)

Figure 4.44: Auto-correlation functions for frequency for PMU’s (a) k = 2, (b) k = 6, (c)k = 112, and (d) k = 139; at t=21:29:00 on July 30, 2013. The correlation matrices areconstructed with the normalized time series fˆ(s) [ϕ̃BPF5 ](·; S), S = 30, T = 5400.

We show study dependence of the auto-correlations on time, i.e. dependence of Ak (·; t)on t. The movies are available in [20], showing examples of buses that have large autocorrelation amplitude. In general, the auto-correlation functions are stable on time (whent varies), indicating that the correlation of the sampling are consistent on quiet periods.
We are interested in understanding the fact that some buses have larger auto-correlations
than others. This can be measured by the amplitude of the auto-correlation functions,
we formalize this concept in the following definition:

CHAPTER 4. LEARNING FROM PMU DATA

110

Definition 22. The residue of the auto-correlation functions is defined as∀k ∈ {1, . . . , N },.
ρk ([∆min , ∆max ]; t) =

∆max
X−1
1
|Ak (∆; t)|,∆max − ∆min ∆=∆(4.28)

minwhere ∆min and ∆max are the initial and ending points of contribution to the residue.
We plot the residue for each sensor in a geographical map, see Figure 4.45. In thesecomputations, we have ignored the first second of the auto-correlation function —that
represent the correlation of the sampling with its immediate past—, by setting ∆min = 30.
The maximum point of contribution has been set to be ∆max = 1800 (1 minute), since
we have empirically verified that with this value of ∆max the results that are obtained
are robust with O(1) changes). In the figure, we show the residual map obtained withthe original normalized time series, using a band-stop filter, and using a band-pass filter.
As seen in the movies of the residue maps shown in Figure 4.45, residue values at anumber of special nodes do not decay with time. Moreover we observe that nodes with asignificant residue cluster, specifically, there are two sets of nodes that repeat their highvalues across different days, months and time of the day. Those groups are also shownin Figure 4.45a: in the bottom-left part of the map and the upper-right zone (closer tothe center); they might show up together or by separate. We also note that the residue
drastically drops when a band-stop filter is applied, as can be seen in Figure 4.45b; andit is accentuated when a band-pass filter is used, see Figure 4.45c.
Observation that sustainable correlations stay for sufficiently long period of timesuggest to analyze spatio-temporal features of the sustainable correlations via the crosscorrelation residue described in next.
Plots for the same auto-correlation analysis for the voltage phase angle and voltage
magnitude measurements are shown in Figures A.1–A.8 in Appendices A.1 and A.2.

CHAPTER 4. LEARNING FROM PMU DATA

111

(a)

(b)

(c)

Figure 4.45: Residue of auto-correlation functions for frequency using (a) m̂(s) (·; S), (b)fˆ(s) [ϕ̃BSF5 ](·; S), and (c) fˆ(s) [ϕ̃BPF5 ](·; S); at t=21:29:00 on July 30, 2013, with S = 30,
T = 5400, ∆min = 30, and ∆max = 1800. Geometrical figures show the position of the sensors
depicted in Figures 4.42–4.44.

4.6.2

Cross-Correlation Residue (CCR)

The cross-correlation version of equations (4.27)–(4.28) is described in the following:
Definition 23. The cross-correlation functions and their residue are defined, respectively, as∀k, ` ∈ {1, . . . , N },.
Bk` (∆; t) = [Σ∆ (t)]k` ,.
Rk` ([∆min , ∆max ]; t) =(4.29)∆max
X−1
1|Bk` (∆; t)|.
∆max − ∆min ∆=∆(4.30)

min

Here, we drop normalization in (4.29) to avoid singularities associated with signals
at different nodes which are not correlated. To visualize the CCR (4.30), we plot theCHAPTER 4. LEARNING FROM PMU DATA

112

cross-correlation residual matrixhi.
R([∆min , ∆max ]; t) = Rk` ([∆min , ∆max ]; t) | ∀k, ` ∈ {1, . . . , N } .

(4.31)In Figure 4.46 we plot the cross-correlation matrix (4.31) component-wise, using darker
colors for higher values.
We observe that for the normalized original time series (see Figure 4.46a), there aresome sensors that have high correlation with others, specifically, given the block structureof the matrix, there is a cluster of PMU sensors that have high correlation between eachothers. Figure 4.46b shows the cross-correlation matrix when a band-stop filter at 5 Hz isapplied, dropping drastically the cross residue, which is consistent with what we observedwith the auto-correlation. In Figure 4.46c we have applied a band-pass filter, and in thiscase, we see that the amount of sensors that are cross-correlated increases.

(a)

(b)

Figure 4.46:(c)

Cross-correlation residual matrix for frequency using (a) m̂(s) (·; S), (b)fˆ(s) [ϕ̃BSF5 ](·; S), and (c) fˆ(s) [ϕ̃BPF5 ](·; S); at t=21:29:00 on July 30, 2013, with S = 30,
T = 5400, ∆min = 30, and ∆max = 1800.

CHAPTER 4. LEARNING FROM PMU DATA

113

For each of the three matrices plotted in Figure 4.46 we plot the components of therow corresponding to sensor k = 134 in order to see the where are the locations of thesensors that are cross-correlated with it. See Figure 4.47.
In the normalized original sampling, Figure 4.47a, we observe a high cross-correlation
between the selected sensor and two clusters located at the bottom-left and upper-rightof the map. Again, this structures repeats from the one observed in Figure 4.45a. As
expected, no high cross-correlation is shown in the measurements filtered with a bandstop filter (in Figure 4.47b); however, we observe high cross-correlation between sensor
k = 134 and most of the other sensor in Figure 4.47c. This might be related with thefact the the oscillations observed at 5 Hz are coordinated between most of the PMU’s.

(a)

(b)

(c)

Figure 4.47: Cross-correlation residue for frequency, geographically located, between sensor
k = 134 (indicated with an extra circle around it) and the remaining sensors using (a) m̂(s) (·; S),(b) fˆ(s) [ϕ̃BSF5 ](·; S), and (c) fˆ(s) [ϕ̃BPF5 ](·; S); at t=21:29:00 on July 30, 2013, with S = 30,
T = 5400, ∆min = 30, and ∆max = 1800.

Plots for the same cross-correlation analysis for the voltage phase angle and voltage

CHAPTER 4. LEARNING FROM PMU DATA

114

magnitude measurements are shown in Figures A.9 and A.10 in Appendices A.3 and A.4.

4.7

Outcomes

We finalize this chapter by summarizing observations that we have made from the statistical analyzis that has been performed.
Averaged over time PMU signal shows interesting spatial correlations. Correlations
are different for different objects of interest (frequency, phase and voltage) and also
different for different quiet periods, altough general structures of the set of sensors thatare correlated and their correlation functions are consistent over time. The matrix ofcorrelations is sparse, also revealing that number of the high-intensity contributions to thecorrelations is small. (Note that the statement of sparseness is consistent with previous
studies of the measurement matrix, see e.g. [57].) Each of the contributions characterizes
a mode localized on a relatively few nodes (PMU positions) within the system. Principal
modes, computed over the quiet periods, are almost frozen in time, however responding
fast to any significant perturbation, thus suggesting them as efficient features/indicatorsfor changes.
Fourier-analysis of the signal reveals interesting spatial patterns. Extracting modesin the 4-6 Hz range one observes significant contributions from only few PMU nodes.
Whereas observing the contribution that the modes at the same range explains the hightemporal auto- and cross-correlation of a significant portion of the sensors. Nodes showing
large 4-6 Hz contributions were identified as aggregated loads, mid-size generators andlarge wind-farms.
In general the signals are long-correlated in time. However the memory effects becomes significantly less prominent when higher frequencies are removed. Nodes withsignificant auto-correlations (memory) have a spatial pattern is adiabatic/frozen (changing slowly during the quiet periods). Like other adiabatic patterns mentioned above, thepattern changes from one quiet period to another and we also observed slight different

CHAPTER 4. LEARNING FROM PMU DATA

115

patterns for different characteristics (frequency, phase, voltage).
Analyzing time-delayed cross-correlations between different nodes we observe that,like in the case of the auto-correlations, correlations between some nodes have long memory. Nodes which mutual inference shows a long memory form a sparse pattern. These
patterns, like others described above are adiabatic and evolving from characteristic tocharacteristic and from one quiet period to another quiet period.

CHAPTER 5. CONCLUSION

116

Chapter 5Conclusion
During the first part of this thesis, we study cyber-physical attacks, attacks that arecomposed by a physical perturbation or disturbance in electrical grid simultaneously
accompanied by a hack on the data produced by the sensors (PMUs) that measure
the status of physical quantities of the network. Cyber-physical attacks have gained
increased attention during the last decade since new technologies have become available
and, therefore, their security systems are not strong enough yet. Hackers have discovered
vulnerabilities and used them to harm the stability of the networks.
We introduce a sophisticated attack model, that includes load modifications and tripping of transmission lines in a small zone of the network, and avoids standard detection
methods by injecting false measurements reflecting realistic safe flows. We are able tocompute such attacks on large systems in seconds of CPU time, with hidden overloads ofmore than 50% of the lines capacity. If these attacks remain undetected for long periodsof time, they could cause catastrophic consequences such as cascade line failures andextended blackouts.
We propose different stochastic defense mechanisms to augment standard detection
tools for generic cyber attacks. Our defenses change the stochastics of system data using
intelligent procedures that modify the power injection at specific generators in a waythat is recognizable by the defender but difficult to anticipate by the attacker. The firstCHAPTER 5. CONCLUSION

117

mechanism reveals the lines that connect the buses that are attacked with the rest of the
network by using the change in generation to induce a change in the voltage phase angleof most of the nodes. Branches that show inconsistencies in terms of the voltage-current
relationship are candidates for being at the boundary of the affected buses.
The second mechanism causes low-rank changes in the covariance (correlation) matrix of phase angles, affecting all the components of the matrix. The randomness during
the iterations of the strategy and the complexity of the power flow equations makes anattacker’s counteraction very difficult and, therefore, the nodes that were hacked will be
exposed by not showing the proper low-rank covariance matrix correction.

The second half of this thesis is dedicated to the analysis of PMU measurements
performed by an Independent System Operator in the United States. The available data
ranges for a period of 15 months, where readings were made at a rate of 30 times per second, across approximately 240 sensors. We searched for intervals of 15 minutes where thesampling (specifically frequency, voltage phase angle and voltage magnitude) reflected anambient conditions behavior and applied signal processing tools and statistical analyses
on these periods.
By focusing on the frequency domain of the time series —obtained by its discretetime Fourier transform— we observe anomalous peaks at 5 Hz, consistently present in alarge portion of the buses across different days and moments of the day. We use Fourier
filtering to isolate the contribution of these frequency modes or to suppress them.
The statistical tools used are principal component analysis (PCA) of the covariance(correlation) matrices, temporal auto-correlation and temporal cross-correlation of thenormalized or filtered time series. From the first one, we conclude that the correlation matrix has low-rank, the largest 10 eigen-values account for more than 80% of the spectrum,with the first one representing more than 25%. We also notice that the eigen-values andeigen-vectors do not change dramatically over time, that is, they have a steady or a slow
changing behavior during the analyzed quiet periods. Leading eigen-vectors are struc-

CHAPTER 5. CONCLUSION

118

tured in such a way that sensors (associated with the components of the eigen-vectors)
having the same value are geographically closed.
The correlation functions that we define account for the temporal correlation thatmight exist between a sensor’s readings and its own past (auto-correlation), and also
between a sensor’s readings and the past of a different sensor (cross-correlation). The
correlation functions that are obtained usually show a periodic behavior and, interestingly, it is mainly due to the modes around 5 Hz observed in the frequency domain of theseries —when the modes near 5 Hz are suppressed, the correlation functions drop theirvalue to zero.
The sensors with high auto-correlation are located in two specific areas of the map.
This behavior is extended through the different observed periods. Moreover, the sensors
located in these areas show high cross-correlation as well, which might indicated thatthey all function in a coordinated fashion.

As an ongoing work, we continue improving the filters that isolate different modes
from the frequency domain of the times series. We will use the filtered times series to infer
the branch admittance parameters across different transmission lines of the network.

BIBLIOGRAPHY

119

Bibliography[1]

Report of WECC Joint Synchronized Information Subcommittee on Modes ofInter-Area Power Oscillations in Western Interconnection. https://www.wecc.
biz/Reliability/WECCJSISModesofInter-AreaOscillations-2013-12-REV1.
1.pdf, 2013. Accessed: 2018-09-11.

[2]

Power System Oscillatory Behaviors:

Sources, Characteristics, & Analyses.

https://www.pnnl.gov/main/publications/external/technical_reports/
PNNL-26375.pdf, 2017. Accessed: 2018-09-11.
[3]

Reliability

Guideline

ForcedOscillation

Monitoring&

Mitigation.

https://www.nerc.com/pa/RAPA/rg/ReliabilityGuidelines/Reliability_
Guideline_-_Forced_Oscillations_-_2017.pdf, 2017. Accessed: 2018-09-11.
[4]

Department
sor

Data.

of

Energy:Big

Data

AnalysisofSynchropha-

https://grantbulletin.research.uiowa.edu/

netl-big-data-analysis-synchrophasor-data, 2018.
[5]

R. Ahuja, T. Magnanti, and J. Orlin. Network Flows: Theory, Algorithms, andApplications. Prentice Hall, 1993.

[6]

E. D. Andersen and K. D. Andersen. The Mosek Interior Point Optimizer for
Linear Programming: An Implementation of the Homogeneous Algorithm, pages
197–232. Springer US, Boston, MA, 2000.

BIBLIOGRAPHY
[7]

120

D. Apostolopoulou, P. W. Sauer, and A. D. Domı́nguez-Garcı́a. Automatic Generation Control and Its Implementation in Real Time. In 2014 47th Hawaii International Conference on System Sciences, pages 2444–2452, January 2014.

[8]

I. Arel, D. Rose, and T. Karnowski. Deep Machine Learning - A New Frontier in Artificial Intelligence Research [Research Frontier]. IEEE Computational Intelligence
Mag., 5(4):13–18, November 2010.

[9]

R. Arghandeh and Y. Zhou. Big Data Application in Power Systems. Elsevier,
2018.

[10] G. Bakke. The Grid: The Fraying Wires between Americans and our Energy Future.
Bloomsbury USA, 2016.
[11] R. Baldick, K. A. Clements, Z. Pinjo-Dzigal, and P. W. Davis. Implementing
Nonquadratic Objective Functions for State Estimation and Bad Data Rejection.
IEEE Transactions on Power Systems, 12(1):376–382, February 1997.
[12] Y. Bengio. Learning Deep Architectures for AI. Foundations and Trends in Machine
Learning, 2(1):1–127, January 2009.
[13] A. R. Bergen and V. Vittal. Power Systems Analysis. Prentice Hall, Upper Saddle
River, NJ, Second edition, 2000.
[14] A. Bernstein, D. Bienstock, D. Hay, M. Uzunoglu, and G. Zussman. Power Grid
Vulnerability to Geographically Correlated Failures Analysis and Control Implications. In IEEE INFOCOM 2014 - IEEE Conference on Computer Communications,pages 2634–2642, April 2014.
[15] S. Bhela, V. Kekatos, and S. Veeramachaneni. Enhancing Observability in Distribution Grids Using Smart Meter Data. IEEE Transactions on Smart Grid,9(6):5953–5961, November 2018.

BIBLIOGRAPHY

121[16] S. Bhela, V. Kekatos, and S. Veeramachaneni. Smart inverter grid probing for learning loads: Part i - identifiability analysis. IEEE Transactions on Power Systems,
to appear, 2019.
[17] S. Bhela, V. Kekatos, and S. Veeramachaneni. Smart inverter grid probing forlearning loads: Part ii - probing injection design. IEEE Transactions on Power
Systems, to appear, 2019.
[18] D. Bienstock. Electrical Transmission System Cascades and Vulnerability: An
Operations Research Viewpoint. Society for Industrial and Applied Mathematics,
Philadelphia, PA, USA, 2015.
[19] D. Bienstock. Machine learning with PMU data. In 2017 NASPI Work Group
Meeting, Gaithersburg, MD, March 2017.
[20] D. Bienstock, M. Chertkov, and M. Escobar. Learning from ISO-scale PMU data
stream, 2018.
[21] D. Bienstock, M. Chertkov, and M. Escobar. Learning from power system data
stream: phasor-detective approach. arXiv:1902.03223, 2018.
[22] D. Bienstock, A. Shukla, and S. Yun.

Non-Stationary Streaming PCA.

arXiv:1902.03223, 2019.
[23] D. Bienstock and A. Verma. The N − k Problem in Power Grids: New Models, Formulations, and Numerical Experiments. SIAM Journal on Optimization,
20(5):2352–2380, 2010.
[24] D. Bienstock and A. Verma. Strong NP-hardness of AC power flows feasibility.
arXiv:1512.07315, December 2015.
[25] R. Bobba, K. Rogers, Q. Wang, H. Khurana, K. Nahrstedt, and T. Overbye. De-BIBLIOGRAPHY

122

tecting False Data Injection Attacks on DC State Estimation. In Proceedings ofthe First Workshop on Secure Control Systems, 2010.
[26] J. Carpentier. Contribution á l’étude du dispatching économique. Bulletin de laSociété Française des Électriciens, 8(3):431–447, 1962.
[27] C. Chang and Z. Li. Recursive stochastic subspace identification for structuralparameter estimation. In Proceedings of SPIE, volume 7292, March 2009.
[28] G. Chen, Z. Y. Dong, D. J. Hill, and Y. S. Xue. Exploring Reliable Strategies for
Defending Power Systems Against Targeted Attacks. IEEE Transactions on Power
Systems, 26(3):1000–1009, August 2011.
[29] L. Cheng, C. You, and L. Chen. Identification of Power Line Outages Based on
PMU Measurements and Sparse Overcomplete Representation. In 2016 IEEE 17th
International Conference on Information Reuse and Integration, pages 343–349,
July 2016.
[30] T. V. Cutsem, M. Ribbens-Pavella, and L. Mili. Hypothesis Testing Identification:
A New Method For Bad Data Analysis In Power System State Estimation. IEEE
Transactions on Power Apparatus and Systems, PAS-103(11):3239–3252, November
1984.
[31] N. Dahal, R. L. King, and V. Madani. Online dimension reduction of synchrophasor
data. In PES T D 2012, pages 1–7, May 2012.
[32] G. Dán and H. Sandberg. Stealth Attacks and Protection Schemes for State Estimators in Power Systems. In 2010 First IEEE International Conference on Smart
Grid Communications, pages 214–219, October 2010.
[33] R. Deng, P. Zhuang, and H. Liang. CCPA: Coordinated Cyber-Physical Attacks
and Countermeasures in Smart Grid. IEEE Transactions on Smart Grid, 8(5):2420–

BIBLIOGRAPHY

1232430, September 2017.
[34] Electricity-Information Sharing and Analysis Center (E-ISAC). Analysis of the
Cyber Attack on the Ukrainian Power Grid. Technical report, SANS Industrial
Control Systems, March 2016.
[35] M. Esmalifalak, H. Nguyen, and R. Zheng. Stealth False Data Injection using
Independent Component Analysis in Smart Grid. In 2011 IEEE International
Conference on Smart Grid Communications, pages 244–248, October 2011.
[36] G. Frigo, C. Narduzzi, D. Colangelo, M. Pignati, and M. Paolone. Definition
and assessment of reference values for PMU calibration in static and transientconditions. In 2016 IEEE International Workshop on Applied Measurements forPower Systems (AMPS), pages 1–6, September 2016.
[37] P. Gao, M. Wang, J. H. Chow, S. G. Ghiocel, B. Fardanesh, G. Stefopoulos, and
M. P. Razanousky. Identification of Successive “Unobservable” Cyber Data Attacks
in Power Systems Through Matrix Decomposition. IEEE Transactions on Signal
Processing, 64(21):5557–5570, November 2016.
[38] P. Gao, M. Wang, S. G. Ghiocel, J. H. Chow, B. Fardanesh, and G. Stefopoulos. Missing Data Recovery by Exploiting Low-Dimensionality in Power System
Synchrophasor Measurements. IEEE Transactions on Power Systems, 31(2):1006–
1013, March 2016.
[39] M. Garcia, T. Catanach, S. V. Wiel, R. Bent, and E. Lawrence. Line Outage Localization Using Phasor Measurement Data in Transient State. IEEE Transactions
on Power Systems, 31(4):3019–3027, July 2016.
[40] A. Giani, E. Bitar, M. Garcia, M. McQueen, P. Khargonekar, and K. Poolla. Smart
Grid Data Integrity Attacks. IEEE Transactions on Smart Grid, 4(3):1244–1253,

BIBLIOGRAPHY

124

September 2013.
[41] J. D. Glover, M. S. Sarma, and T. J. Overbye. Power Systems Analysis and Design.
Cengage Learning, Stamford, CT, Fourth edition, 2008.
[42] I. Goodfellow, Y. Bengio, and A. Courville. Deep Learning. The MIT Press, 2016.
[43] Gurobi Optimization, LLC.

Gurobi Optimizer Reference Manual, 2018.

http://www.gurobi.com.
[44] F. J. Harris. On the Use of Windows for Harmonic Analysis with the Discrete
Fourier Transform. Proceedings of the IEEE, 66(1):51–83, January 1978.
[45] J. F. Hauer, C. J. Demeure, and L. L. Scharf. Initial results in Prony analysis ofpower system response signals. IEEE Transactions on Power Systems, 5(1):80–89,
February 1990.
[46] Y. Huang, S. Werner, J. Huang, N. Kashyap, and V. Gupta. State Estimation inElectric Power Grids: Meeting New Challenges Presented by the Requirements of
the Future Grid. IEEE Signal Processing Magazine, 29(5):33–43, September 2012.
[47] G. Hug and J. A. Giampapa. Vulnerability Assessment of AC State Estimation
With Respect to False Data Injection Cyber-Attacks. IEEE Transactions on Smart
Grid, 3(3):1362–1370, September 2012.
[48] International Bureau of Weights and Measures (BIPM). The International System
of Units (SI). Brochure, 2006.
[49] M. Kezunovic, L. Xie, and S. Grijalva. The role of big data in improving powersystem operation and protection. In 2013 IREP Symposium Bulk Power System
Dynamics and Control - IX Optimization, Security and Control of the Emerging
Power Grid, pages 1–9, August 2013.

BIBLIOGRAPHY

125[50] J. Kim and L. Tong. On Topology Attack of a Smart Grid: Undetectable Attacks and Countermeasures. IEEE Journal on Selected Areas in Communications,31(7):1294–1305, July 2013.
[51] J. Kim, L. Tong, and R. J. Thomas. Subspace Methods for Data Attack on State
Estimation: A Data Driven Approach. IEEE Transactions on Signal Processing,
63(5):1102–1114, March 2015.
[52] T. T. Kim and H. V. Poor. Strategic Protection Against Data Injection Attacks
on Power Grids. IEEE Transactions on Smart Grid, 2(2):326–333, June 2011.
[53] B. J. Kirby, J. Dyer, C. Martinez, R. A. Shoureshi, R. Guttromson, and J. Dale.
Frequency Control Concerns in the North American Electric Power System. Technical Report ORNL/TM-2003/41, U.S. Department of Energy, December 2002.
[54] B. Kocuk, S. S. Dey, and X. A. Sun. Inexactness of SDP Relaxation and Valid
Inequalities for Optimal Power Flow.

IEEE Transactions on Power Systems,

31(1):642–651, January 2016.
[55] B. Kocuk, S. S. Dey, and X. A. Sun. Strong SOCP Relaxations for the Optimal
Power Flow Problem. Operations Research, 64(6):1177–1196, 2016.
[56] Y. LeCun, Y. Bengio, and G. Hinton. Deep Learning. Nature, 521(7553):436–444,
May 2015.
[57] W. Li, M. Wang, and J. Chow. Real-Time Event Identification Through LowDimensional Subspace Characterization of High-Dimensional Synchrophasor Data.
IEEE Trans. on Power Systems, 33(5):4937–4947, September 2018.
[58] Z. Li, M. Shahidehpour, A. Alabdulwahab, and A. Abusorrah. Bilevel Model for
Analyzing Coordinated Cyber-Physical Attacks on Power Systems. IEEE Transactions on Smart Grid, 7(5):2260–2272, September 2016.

BIBLIOGRAPHY

126[59] J. Liang, O. Kosut, and L. Sankar. Cyber Attacks on AC State Estimation: Unobservability and Physical Consequences. In 2014 IEEE PES General Meeting —Conference Exposition, pages 1–5, July 2014.
[60] J. Liang, L. Sankar, and O. Kosut. Vulnerability Analysis and Consequences of
False Data Injection Attack on Power System State Estimation. IEEE Transactions
on Power Systems, 31(5):3864–3872, September 2016.
[61] J. Lin and H. Pan. A Static State Estimation Approach Including Bad Data
Detection and Identification in Power Systems. In 2007 IEEE Power Engineering
Society General Meeting, pages 1–7, June 2007.
[62] G. Liu, J. Quintero, and V. M. Venkatasubramanian. Oscillation monitoring system
based on wide area synchrophasors in power systems. In 2007 iREP Symposium- Bulk Power System Dynamics and Control - VII. Revitalizing Operational Reliability, pages 1–13, August 2007.
[63] X. Liu and Z. Li. Local Load Redistribution Attacks in Power Systems With
Incomplete Network Information. IEEE Transactions on Smart Grid, 5(4):1665–
1676, July 2014.
[64] Y. Liu, P. Ning, and M. K. Reiter. False Data Injection Attacks Against State
Estimation in Electric Power Grids. In Proceedings of the 16th ACM Conference
on Computer and Communications Security, CCS ’09, pages 21–32, New York, NY,
USA, 2009. ACM.
[65] A. Lokhov, M. Vuffray, D. Shemetov, D. Deka, and M. Chertkov. Online Learningof Power Transmission Dynamics. PSCC 2018, arXiv:1710.10021, 2017.
[66] N. M. Manousakis, G. N. Korres, and P. S. Georgilakis. Taxonomy of PMU Placement Methodologies. IEEE Transactions on Power Systems, 27(2):1070–1077, May

BIBLIOGRAPHY

127

2012.
[67] D. K. Molzahn and J. Wang. Detection and Characterization of Intrusions to
Network Parameter Data in Electric Power Systems. IEEE Transactions on Smart
Grid, to appear, 2018.
[68] A. Monticelli and A. Garcia. Reliable Bad Data Processing for Real-Time State Estimation. IEEE Transactions on Power Apparatus and Systems, PAS-102(5):1126–
1139, May 1983.
[69] A. Monticelli, F. F. Wu, and M. Yen. Mutiple Bad Data Identification for State
Estimation by Combinatorial Optimization. IEEE Transactions on Power Delivery,
1(3):361–369, July 1986.
[70] K. Narendra, D. Rangana, and A. Rajapakse. Dynamic Performance Evaluation
and Testing of Phasor Measurement Unit (PMU) as per IEEE C37.118.1 Standard. In Doble Client Committee Meetings & International Protection Testing
Users Group (PTUG), October 2010.
[71] A. V. Oppenheim, R. W. Schafer, and J. R. Buck. Discrete-Time Signal Processing.
Prentice-Hall, Inc., Upper Saddle River, NJ, USA, Second edition, 1999.
[72] A. Pinar, J. Meza, V. Donde, and B. Lesieutre. Optimization Strategies for the
Vulnerability Analysis of the Electric Power Grid. SIAM Journal on Optimization,
20(4):1786–1810, February 2010.
[73] M. A. Rahman and H. Mohsenian-Rad. False Data Injection Attacks with Incomplete Information Against Smart Power Grids. In 2012 IEEE Global Communications Conference, pages 3153–3158, December 2012.
[74] S. Roy and B. Lesieutre. Frequency Band Decomposition of a Dynamic Persistence
Measure Using Ambient Synchrophasor Data. In 2018 Power Systems Computation

BIBLIOGRAPHY

128

Conference (PSCC), pages 1–7, June 2018.
[75] S. Sarmadi and V. Venkatasubramanian. Inter-Area Resonance in Power Systems
From Forced Oscillations. IEEE Trans. on Power Systems, 31(1):378–386, January
2016.
[76] J. Schmidhuber. Deep Learning in Neural Networks. Neural Networks, 61(C):85–
117, January 2015.
[77] F. C. Schweppe. Power System Static-State Estimation, Part III: Implementation.
IEEE Transactions on Power Apparatus and Systems, PAS-89(1):130–135, January
1970.
[78] F. C. Schweppe and D. B. Rom. Power System Static-State Estimation, Part II:
Approximate Model. IEEE Transactions on Power Apparatus and Systems, PAS89(1):125–130, January 1970.
[79] F. C. Schweppe and J. Wildes. Power System Static-State Estimation, Part I: Exact
Model. IEEE Transactions on Power Apparatus and Systems, PAS-89(1):120–125,
January 1970.
[80] S. Soltan, M. Yannakakis, and G. Zussman. Power Grid State Estimation Followinga Joint Cyber and Physical Attack. IEEE Transactions on Control of Network
Systems, 5(1):499–512, March 2018.
[81] S. Soltan, M. Yannakakis, and G. Zussman. REACT to Cyber Attacks on Power
Grids. IEEE Transactions on Network Science and Engineering, to appear, 2018.
[82] S. Soltan and G. Zussman. Power Grid State Estimation after a Cyber-PhysicalAttack under the AC Power Flow Model. In 2017 IEEE Power Energy Society
General Meeting, pages 1–5, July 2017.

BIBLIOGRAPHY

129[83] S. Soltan and G. Zussman. EXPOSE the Line Failures following a Cyber-PhysicalAttack on the Power Grid. IEEE Transactions on Control of Network Systems, toappear, 2018.
[84] Y. Tang and G. N. Stenbakken. Calibration of phasor measurement unit at NIST.
In 2012 Conference on Precision electromagnetic Measurements, pages 414–415,
July 2012.
[85] J. E. Tate and T. J. Overbye. Line Outage Detection Using Phasor Angle Measurements. IEEE Transactions on Power Systems, 23(4):1644–1652, November 2008.
[86] J. E. Tate and T. J. Overbye. Double Line Outage Detection Using Phasor Angle
Measurements. In 2009 IEEE Power Energy Society General Meeting, pages 1–5,
July 2009.
[87] A. Teixeira, S. Amin, H. Sandberg, K. H. Johansson, and S. S. Sastry. Cyber
Security Analysis of State Estimators in Electric Power Systems. In 49th IEEE
Conference on Decision and Control (CDC), pages 5991–5998, December 2010.
[88] D.
theTrudnowski.
WECC

Properties
Interconnect.

of

the

Dominant

Inter-Area

Modes

in

https://www.wecc.biz/Reliability/

WECCmodesPaper130113Trudnowski.pdf, 2012. Accessed: 2018-09-11.
[89] D. Trudnowski and J. Pierre. Signal Processing Methods for Estimating SmallSignal Dynamic Properties from Measured Responses, pages 1–36. Springer US,
Boston, MA, 2009.
[90] D. Van Hertem, J. Verboomen, K. Purchala, R. Belmans, and W. L. Kling. Usefulness of DC Power Flow for Active Power Flow Analysis with Flow Controlling
Devices. In The 8th IEE International Conference on AC and DC Power Transmission, pages 58–62, March 2006.

BIBLIOGRAPHY

130[91] O. Vuković, K. C. Sou, G. Dán, and H. Sandberg. Network-layer Protection
Schemes against Stealth Attacks on State Estimators in Power Systems. In 2011
IEEE International Conference on Smart Grid Communications, pages 184–189,
October 2011.
[92] A. Wächter and L. T. Biegler. On the implementation of an interior-point filter
line-search algorithm for large-scale nonlinear programming. Mathematical Programming, 106(1):25–57, March 2006.
[93] L. Xie, Y. Chen, and P. R. Kumar. Dimensionality Reduction of Synchrophasor
Data for Early Event Detection: Linearized Analysis. IEEE Transactions on Power
Systems, 29(6):2784–2794, November 2014.
[94] L. Xie, Y. Chen, and P. R. Kumar. Dimensionality Reduction of Synchrophasor
Data for Early Event Detection: Linearized Analysis. IEEE Transactions on Power
Systems, 29(6):2784–2794, November 2014.
[95] Y. Yuan, Z. Li, and K. Ren. Modeling Load Redistribution Attacks in Power
Systems. IEEE Transactions on Smart Grid, 2(2):382–390, June 2011.
[96] J. Zhang and L. Sankar. Physical System Consequences of Unobservable State-andTopology Cyber-Physical Attacks. IEEE Transactions on Smart Grid, 7(4):2016–
2025, July 2016.
[97] Y. Zhao, J. Chen, A. Goldsmith, and H. V. Poor. Identification of Outages in Power
Systems With Uncertain States and Optimal Sensor Locations. IEEE Journal of
Selected Topics in Signal Processing, 8(6):1140–1153, December 2014.
[98] Y. Zhao, A. Goldsmith, and H. V. Poor. On PMU Location Selection for Line
Outage Detection in Wide-area Transmission Networks. In 2012 IEEE Power andEnergy Society General Meeting, pages 1–8, July 2012.

BIBLIOGRAPHY

131[99] N. Zhou, D. J. Trudnowski, J. W. Pierre, and W. A. Mittelstadt. Electromechanical Mode Online Estimation Using Regularized Robust RLS Methods. IEEE
Transactions on Power Systems, 23(4):1670–1680, November 2008.
[100] H. Zhu and G. B. Giannakis. Sparse Overcomplete Representations for Efficient
Identification of Power Line Outages.

IEEE Transactions on Power Systems,

27(4):2215–2224, November 2012.
[101] R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas. MATPOWER:Steady-State Operations, Planning, and Analysis Tools for Power Systems Researchand Education. IEEE Transactions on Power Systems, 26(1):12–19, February 2011.

132

Appendices

APPENDIX. CORRELATION PLOTS

133

Appendix
Correlation Plots
A.1

Auto-Correlation: Voltage Phase Angle

(a)

(b)

(c)

(d)

Figure A.1: Auto-correlation functions for voltage phase angle for PMU’s (a) k = 2, (b) k = 6,(c) k = 112, and (d) k = 139; at t=21:29:00 on July 30, 2013. The correlation matrices areconstructed with the normalized time series m̂(s) (·; S), S = 30, T = 5400.

APPENDIX. CORRELATION PLOTS

134

(a)

(b)

(c)

(d)

Figure A.2: Auto-correlation functions for voltage phase angle for PMU’s (a) k = 2, (b) k = 6,(c) k = 112, and (d) k = 139; at t=21:29:00 on July 30, 2013. The correlation matrices
are constructed with the normalized time series fˆ(s) [ϕ̃BSF5 ](·; S) that has been filtered by aband-stop Fourier filter, S = 30, T = 5400.

APPENDIX. CORRELATION PLOTS

135

(a)

(b)

(c)

(d)

Figure A.3: Auto-correlation functions for voltage phase angle for PMU’s (a) k = 2, (b) k = 6,(c) k = 112, and (d) k = 139; at t=21:29:00 on July 30, 2013. The correlation matrices
are constructed with the normalized time series fˆ(s) [ϕ̃BPF5 ](·; S) that has been filtered by aband-pass Fourier filter, S = 30, T = 5400.

APPENDIX. CORRELATION PLOTS

136

(a)

(b)

(c)

Figure A.4: Residue of auto-correlation functions for voltage phase angle using (a) m̂(s) (·; S),(b) fˆ(s) [ϕ̃BSF5 ](·; S), and (c) fˆ(s) [ϕ̃BPF5 ](·; S); at t=21:29:00 on July 30, 2013, with S = 30,
T = 5400, ∆min = 30, and ∆max = 1800. Geometrical figures show the position of the sensors
depicted in Figures A.1–A.3.

APPENDIX. CORRELATION PLOTS

A.2

137

Auto-Correlation: Voltage Magnitude

(a)

(b)

(c)

(d)

Figure A.5: Auto-correlation functions for voltage magnitude for PMU’s (a) k = 2, (b) k = 6,(c) k = 112, and (d) k = 139; at t=21:29:00 on July 30, 2013. The correlation matrices areconstructed with the normalized time series m̂(s) (·; S), S = 30, T = 5400.

APPENDIX. CORRELATION PLOTS

138

(a)

(b)

(c)

(d)

Figure A.6: Auto-correlation functions for voltage magnitude for PMU’s (a) k = 2, (b) k = 6,(c) k = 112, and (d) k = 139; at t=21:29:00 on July 30, 2013. The correlation matrices
are constructed with the normalized time series fˆ(s) [ϕ̃BSF5 ](·; S) that has been filtered by aband-stop Fourier filter, S = 30, T = 5400.

APPENDIX. CORRELATION PLOTS

139

(a)

(b)

(c)

(d)

Figure A.7: Auto-correlation functions for voltage magnitude for PMU’s (a) k = 2, (b) k = 6,(c) k = 112, and (d) k = 139; at t=21:29:00 on July 30, 2013. The correlation matrices
are constructed with the normalized time series fˆ(s) [ϕ̃BPF5 ](·; S) that has been filtered by aband-pass Fourier filter, S = 30, T = 5400.

APPENDIX. CORRELATION PLOTS

140

(a)

(b)

(c)

Figure A.8: Residue of auto-correlation functions for voltage magnitude using (a) m̂(s) (·; S),(b) fˆ(s) [ϕ̃BSF5 ](·; S), and (c) fˆ(s) [ϕ̃BPF5 ](·; S); at t=21:29:00 on July 30, 2013, with S = 30,
T = 5400, ∆min = 30, and ∆max = 1800. Geometrical figures show the position of the sensors
depicted in Figures A.5–A.7.

APPENDIX. CORRELATION PLOTS

A.3

141

Cross-Correlation: Voltage Phase Angle

(a)

(b)

(c)

(d)

(e)

(f)

Figure A.9: Cross-correlation residual matrix (on the left plots) and cross-correlation residue
for voltage phase angle, geographically located, between sensor k = 134 —indicated with an
extra circle around it— and the remaining sensors (on the right plots) using m̂(s) (·; S) in (a)and (b), fˆ(s) [ϕ̃BSF5 ](·; S) in (c) and (d), and fˆ(s) [ϕ̃BPF5 ](·; S) in (e) and (f); at t=21:29:00 on
July 30, 2013, with S = 30, T = 5400, ∆min = 30, and ∆max = 1800.

APPENDIX . CORRELATION PLOTS

A.4

142

Cross-Correlation: Voltage Magnitude

(a)

(b)

(c)

(d)

(e)

(f)

Figure A.10: Cross-correlation residual matrix (on the left plots) and cross-correlation residue
for voltage magnitude, geographically located, between sensor k = 134 —indicated with an
extra circle around it— and the remaining sensors (on the right plots) using m̂(s) (·; S) in (a)and (b), fˆ(s) [ϕ̃BSF5 ](·; S) in (c) and (d), and fˆ(s) [ϕ̃BPF5 ](·; S) in (e) and (f); at t=21:29:00 on
July 30, 2013, with S = 30, T = 5400, ∆min = 30, and ∆max = 1800.