# Grid-Graph Signal Processing (Grid-GSP):
A Graph Signal Processing Framework
for the Power Grid

Raksha Ramakrishna and Anna Scaglione Work done at the School of Electrical,
Computer and Energy Engineering (ECEE), Arizona State University, Tempe, AZ.
This research was supported in part by the Director, Office of Electricity
Delivery and Energy Reliability, Cybersecurity for Energy Delivery Systems
program, of the U.S. Department of Energy, under contract DOE0000780. Any
opinions, and findings expressed in this material are those of the authors and
do not necessarily reflect those of the sponsors Preliminary work was
presented in [1, 2].

###### Abstract

The underlying theme of this paper is to explore the various facets of power
systems data through the lens of graph signal processing (GSP), laying down
the foundations of the Grid-GSP framework. Grid-GSP provides an interpretation
for the spatio-temporal properties of voltage phasor measurements, by showing
how the well-known power systems modeling supports a generative low-pass graph
filter model for the state variables, namely the voltage phasors. Using the
model we formalize the empirical observation that voltage phasor measurement
data lie in a low-dimensional subspace and tie their spatio-temporal structure
to generator voltage dynamics. The Grid-GSP generative model is then
successfully employed to investigate the problems pertaining to the grid of
data sampling and interpolation, network inference, detection of anomalies and
data compression. Numerical results on a large synthetic grid that mimics the
real-grid of the state of Texas, ACTIVSg2000, and on real-world measurements
from ISO-New England verify the efficacy of applying Grid-GSP methods to
electric grid data.

##  I Introduction

The power grid is one of the foremost examples of a large-scale man-made
network. The nodes of the associated graph are the grid buses and its edges
are its transmission lines. It is therefore natural to see measurements from
the power grid as graph signals [3] and model power grid measurements using
tools from the theory of graph signal processing (GSP) whose goal is to extend
fundamental insights that come from the frequency analysis for time series to
the domain of signals indexed by graphs [3, 4, 5]. One of the factors that
motivate the development of GSP for the power grid is the abundance of high-
quality data that can be acquired using phasor measurement units (PMU)s, the
sensors producing estimates of the voltage and current phasors [6]. With that,
classical signal processing questions pertaining to sampling, interpolation,
denoising and compression and questions that hinge on the underlying structure
of the voltage phasors graph signal arise.

The overarching goal of this paper is to develop GSP based models for power
systems from first principles by building upon the existing system-level
knowledge of power systems to create a solid foundation to analyze power-grid
measurements using tools from GSP. We call this the Grid-GSP framework. By
identifying the correct graph shift operators (GSO), we extend well-known
results in GSP to power system data without losing the associated physical
interpretation.

The core idea is to rewrite the differential algebraic equations (DAE) [7], in
a way often done in transient stability analysis of power systems, to reveal
that the inherent structure in voltage phasors can be explained using a linear
low-pass graph filter as a generative model, whose inputs are the generator
voltages. This input signal is the generators’ response to electric load in
the grid. Through this model the paper shows also that the temporal dynamics
of the input signal, i.e. the generator voltages, can be explained using a
non-linear GSP model defined via another GSO derived from the generator-only
Kron-reduced network. This is done utilizing the well-known classical swing
equations [8, 9].  This spatio-temporal generative model supports the
empirical observation that voltage data obtained using PMUs tend to be
confined to a much smaller dimension compared to the size of the data record
in both space and time [10, 11]. Many papers have leveraged the empirical
observation of the low-rank of phasor data for the interpolation of missing
data [11], correcting bad data [12] and to detect faulty events [13, 14, 10,
15]. Importantly, our framework explicitly puts forth the structure of this
low-dimensional subspace using our GSP-based generative model, directly tying
this subspace to the graph Fourier domain of the GSO.

###  I-A Literature review

We review prior works by dividing the most relevant literature into three
categories: 1) a general survey of works that use concepts from graph theory
and GSP in power systems in the areas of sensor placement, interpolation and
network inference, 2) FDI attack detection and 3) literature pertaining to
data compression for PMU data.

Graph theory for power systems: Several papers have used insights from
spectral and algebraic graph theory. A few applications include optimal
placement [16, 17] and generating statistically accurate topologies [18]. Grid
topology identification is a network inference problem and has been studied by
several works such as in [19, 20, 21, 22, 23]. GSP concepts have been
leveraged in [24, 25] to detect FDI attacks. Prior work in [26] dealt with
performance limits on fault localization with inadequate number of PMUs and
connected it with graph signal sampling theory and optimal placement of PMUs
for best possible resolution of fault localization in this under-sampled
regime.

The Kron-reduced network among the generator buses and the associated
properties are used in [27, 28] to detect low-frequency oscillations as well
as the resulting islanding patterns. In [29], the authors have shed light on
the relationship that exists between graph Laplacian and modes in power
systems. Recently, a comprehensive review of graph-theoretical concepts in
power systems is presented in [7].

Additionally, there have been several papers adopting graphical models for
state estimation [30], topology estimation [31] and optimal power flow [32].
While the modeling approach is valid, the graphical models capture correlation
whereas our method models the underlying cause for that correlation structure
thereby opening the door for statistical and non-statistical approaches.

Note that with the exception of [24, 25], no other papers make the connection
with GSP and even in the aforementioned work, GSP is used in an empirical
manner. On the other hand, in our preliminary work in [1, 2], established a
case for GSP in a more rigorous manner.

Detection of FDI attacks: While anomaly detection can be broadly applied to
identify various events, a large body of prior research has focused on False
Data Injection (FDI) attacks that can bypass classical bad data detection
(BDD) mechanisms [33, 34, 35, 36] and trigger incorrect decisions or hide line
overflows and contingencies [37, 38, 39]). FDI attacks to PMUs can, for
instance, exploit the vulnerability of GPS signals to spoofing attacks [40,
41].
Recent papers on PMU data integrity have proposed leveraging the low rank
spatio-temporal nature of PMU data not only to help with erasures but also to
strengthen conventional BDD mechanisms [42]. In [38] the authors have
suggested an FDI attack strategy that can pass the aforementioned BDD approach
in [42] by generating false samples that approximately preserve the original
subspace structure.

Compression of PMU data: Due to their relatively high sampling rate and their
wide deployment, PMU data have have called for compression. The performance of
several off-the-shelf lossless encoding techniques applied to PMU data were
investigated in [43]. In [44], a lossless compression called slack referenced
encoding (SRE) method of PMU voltage measurements is introduced, by
identifying a slack-bus and differentially encoding the difference between
slack-bus measurements in time and all other buses. Similarly, in [45], phasor
angle data is encoded in a lossless manner by preprocessing data using
techniques from [44] and then using Golomb-Rice entropy encoding.
The idea of slow-variation with time in PMU data such as phase angles is used
in [46] to firstly transform data into frequency domain (in time) and then
using a ‘reverse water-filling’ technique to encode frequency components in
the difference measurements. Many lossy compression techniques utilize the
low-rank structure inherent in voltage phasor data [47, 48, 49]. Several other
wavelet based compression algorithms exist in the literature as well [50].

###  I-B Contributions

The aim of this paper is to establish the framework of Grid-GSP and elucidate
properties of power grid signals using tools from GSP. In particular, we:

  1. 1.

Establish that PMU voltage measurements from the power grid are result of an
excitation to a low-pass graph filter whose graph shift operator (GSO) is
defined using a function of the system admittance matrix. This was partly
explored in our previous work in [1, 2] and used for blind community
detection.

  2. 2.

Study the spatio-temporal structure of the excitation that, at a fast time
scale, is dominated by the generators dynamics. It is shown that this
excitation can be modeled as an auto-regressive graph filter [51] (GF-AR (2))
for the input signals from the generator internal bus. The spatial properties
in quasi-steady state are captured by defining another GSO using the
‘generator-only’ or Kron-reduced [9] network for the generator buses.

These models set the foundations to revisit known GSP based algorithms for
sampling and reconstruction, interpolation and denoising and network inference
in the context of signal-processing PMU data. We harness the GFT for feature
extraction, to detect anomalies (specifically, FDI attacks [2]) and to derive
a lossy PMU voltage data compression algorithm, leveraging the sparsity of the
GFT signal. By elucidating all the steps in modeling power systems data from
graph construction, signal model, identification of the low-pass structure
that is responsible for low-dimensional representation to signal denoising,
network inference and anomaly detection, we illustrate how one can similarly
develop GSP based models in other application domains especially data that can
be modeled as the output of low-pass graph signals [52].

###  I-C Paper organization

Section II reviews concepts from GSP focusing on complex-valued graph signals
and applicable more broadly to bandpass signals whose signal models rely on
complex envelopes or phasors. It also reviews measurements and parameters
pertaining to the grid. Section III lays the foundation for Grid-GSP mapping
the physical laws to a spatio-temporal generative model for voltage signals.
Through these lens, in Section IV, the paper revisits algorithms and tools
from GSP for PMU data pertaining to sampling and reconstruction along with
optimal placement of PMUs, interpolation of missing samples and network
inference. Section V highlights applications of Grid-GSP to detect FDI attacks
and for sequential lossy voltage data compression. The algorithms and methods
are tested numerically in Section VI. Section VII summarizes conclusions and
future research directions.

Notation: Boldfaced lowercase letters are used for vectors, $\bm{x}$ and
uppercase for matrices, $\mathbf{A}$. Transpose is
$\bm{x}^{\top},\mathbf{A}^{\top}$ and conjugate transpose is $\mathbf{A}^{H}$.
$\left[\bm{x}\right]_{\mathcal{M}}$ is the new vector that has elements of the
vector indexed by the set $\mathcal{M}$. The operation $\Re\\{.\\},\Im\\{.\\}$
denote the real and imaginary parts of the argument. Pseudo-inverse of a
matrix is $\mathbf{A}^{\dagger}$. The operation $\texttt{diag}(\bm{x})$
creates a diagonal matrix with elements from a vector $\bm{x}$ and
$\texttt{Diag}(\mathbf{A})$ is the extraction of diagonal values of a matrix
$\mathbf{A}$.

##  II Preliminaries

###  II-A Graph Signal Processing (GSP) in a nutshell

Consider an undirected graph
$\mathcal{G}=\left(\mathcal{N},\mathcal{E}\right)$ with nodes
$i\in\mathcal{N}$ and edges $(i,j)\in\mathcal{E}$. A graph signal
$\bm{x}\in\mathbb{C}^{\lvert\mathcal{N}\rvert}$ is a vector whose $i$th entry
$\left[\bm{x}\right]_{i}$ is associated to node $i\in\mathcal{N}$. The set of
nodes connected to node $i$ is called the neighborhood of $i$ and denoted as
$\mathcal{N}_{i}$. GSP generalizes the notion of discrete time shift for a
time series by introducing the notion of graph shift operator (GSO):

######  Definition 1.

A graph shift operator (GSO) is a linear neighborhood operator, so that each
entry of the shifted graph signal is a linear combination of the graph signal
neighbors’ values [5].

The linear combinations can use complex-valued weights
$s_{ij}\in\mathbb{C},~{}\text{with}~{}s_{ij}=0,(i,j)\notin\mathcal{E}$ and,
clearly, the GSO can be defined as matrix multiplication, with a matrix
$\mathbf{S}\in\mathbb{C}^{\lvert\mathcal{N}\rvert\times\lvert\mathcal{N}\rvert}$.
Although not the only option, one common choice for the GSO, is that of the
graph weighted Laplacian111In truth, the Laplacian should be considered a
graph differential operator as opposed to a shift operator, but we use the
conventional name nonetheless in the rest of the paper to be consistent with
the literature., i.e:

$\left[\mathbf{S}\right]_{i,j}=\begin{cases}\sum_{k\in\mathcal{N}_{i}}s_{i,k},~{}i=j\\\
-s_{i,j},~{}i\neq j\end{cases}$ (1)

In this work, we focus on complex symmetric GSOs,
$\mathbf{S}=\mathbf{S}^{\top}$ as is applicable to the weighted graph
Laplacian for the power grid. Having defined the notion of shift, one can
introduce the notion of shift-invariance:

######  Definition 2.

Given a GSO $\mathbf{S}$ a shift invariant operator ${\cal H}$ acting on a
graph signal is such that:

${\cal H}:\bm{x}\mapsto\bm{v},~{}~{}~{}~{}\Longleftrightarrow~{}~{}~{}~{}{\cal
H}:\mathbf{S}\bm{x}\mapsto\mathbf{S}\bm{v}.$ (2)

Linear shift-invariant operators must be matrix polynomials of the GSO
$\mathbf{S}$ [53].  Therefore, a linear shift-invariant graph filter is a
linear operator and can be defined as:

$\textstyle{\bm{v}}={\cal H}(\mathbf{S})\bm{x},~{}~{}{\cal
H}\left(\mathbf{S}\right)\\!=\\!\textstyle\sum\limits_{k=0}^{K-1}h_{k}\mathbf{S}^{k}$
(3)

Additionally, linear shift-invariant graph filters satisfy the condition:
$\mathbf{S}{\cal H}\left(\mathbf{S}\right)={\cal
H}\left(\mathbf{S}\right)\mathbf{S}$.

Consider the following eigenvalue decomposition of the complex symmetric GSO
$\mathbf{S}$, given by Theorem $4.4.13$ in [54] for diagonalizable complex
symmetric matrices:

$\displaystyle\mathbf{S}$
$\displaystyle=\mathbf{U}\bm{\Lambda}\mathbf{U}^{\top},$
$\displaystyle\mathbf{U}^{\top}\mathbf{U}=\mathbf{U}\mathbf{U}^{\top}=\mathbb{I}.$
(4)

Here $\bm{\Lambda}$ is the diagonal matrix with eigenvalues
$\lambda_{0},\lambda_{1},\dots\lambda_{\lvert\mathcal{N}\rvert-1}$ on the
principal diagonal and $\mathbf{U}$ are complex orthogonal eigenvectors. An
equivalent concept of frequency domain in GSP is defined using eigenvalues and
eigenvectors of the GSO.

Graph frequencies are the eigenvalues of the GSO and the order of frequencies
is based on the total variation (TV) criterion [5, 55] defined using the
discrete $p$ Dirichlet form $S_{p}\\!\\!\left(\bm{x}\right)$ with $p=1$ as in
[56] as:

$\displaystyle
S_{1}\\!\\!\left(\bm{x}\right)\\!=\\!\lVert\mathbf{S}\bm{x}\rVert_{1}\implies
S_{1}\\!\left(\bm{u}_{i}\right)=\lVert\mathbf{S}\bm{u}_{i}\rVert_{1}=\lvert\lambda_{i}\rvert\lVert\bm{u}_{i}\rVert_{1}$
(5)

After normalizing the eigenvectors such that
$\lVert\bm{u}_{i}\rVert_{1}=1~{}\forall i$, it is clear that
$S_{1}\left(\bm{u}_{i}\right)>S_{1}\left(\bm{u}_{j}\right)\implies\lvert\lambda_{i}\rvert>\lvert\lambda_{j}\rvert.$
Hence, the ascending order of eigenvalues corresponds to increase in
frequency,
$\lvert\lambda_{0}\rvert=0\leq\lvert\lambda_{1}\rvert\leq\lvert\lambda_{2}\rvert\dots\lvert\lambda_{\lvert\mathcal{N}\rvert}\rvert$.
This ordering is not unique since two distinct complex eigenvalues can have
the same magnitude.

The Graph Fourier Transform (GFT) basis is the complex orthogonal basis
$\mathbf{U}$ in (4). Hence, the GFT of a graph signal $\bm{x}$,
$\tilde{\bm{x}}$ and the inverse GFT are given by
$\tilde{\bm{x}}=\mathbf{U}^{\top}\bm{x}$ and $\bm{x}=\mathbf{U}\tilde{\bm{x}}$
respectively where $\left[\tilde{\bm{x}}\right]_{m}$ is the frequency
component that corresponds to the $m$-th eigenvalue $\lambda_{m}$222It is
worth noting that the graph shift operator and Fourier transforms do not have
in general important properties that are found in their conventional
counterparts for time series. One notable fact is that the spectrum of
${\bm{S}}\bm{x}$ does not have the same amplitude as the spectrum of $\bm{x}$.
In fact, the GSO effect is closer to that of a derivative, since each of the
GFT coefficients is rescaled by the corresponding frequency. For complex
symmetric, rather than Hermitian operators, unfortunately also Parseval
theorem is not valid.. Also, we can define the graph-frequency response of the
graph filter, $\tilde{\bm{h}}$, by writing

$\displaystyle\textstyle\mathcal{H}\left(\mathbf{S}\right)$
$\displaystyle\\!=\\!\\!\textstyle\sum\limits_{k=0}^{K}h_{k}\mathbf{S}^{k}\\!=\\!\mathbf{U}\texttt{diag}(\tilde{\bm{h}})\mathbf{U}^{\top},$
(6) $\displaystyle[\tilde{\bm{h}}]_{i}$
$\displaystyle\\!\\!=\\!\\!\mathrm{H}(\lambda_{i}),~{}~{}~{}~{}\mathrm{H}(\lambda):=\sum\limits_{k=0}^{K}h_{k}\lambda^{k}$
(7)

$h_{k}\in\mathbb{C},~{}i=1,2,\dots\lvert\mathcal{N}\rvert$. The frequency
response of the filter is given by elements in $\tilde{\bm{h}}$. Subsequently,
the input and output of a graph filter in graph-frequency domain are related
as

$\displaystyle{\bm{v}}={\cal
H}(\mathbf{S})\bm{x}~{}~{}\rightarrow~{}~{}\tilde{\bm{v}}=\texttt{diag}(\tilde{\bm{h}})\tilde{\bm{x}},$
(8)

which is analogous to convolution theorem for time-domain signals. Naturally,
this leads to the extension of notions such as low-pass, high-pass and band-
pass filters and signals that are at the heart of sampling and interpolation
schemes.

###  II-B GSP for time series of graph signals

So far, only the nodal index for the graph signal $\bm{v}$ was considered.
However, one can also encounter graph signal processes i.e. temporal
variations in a graph signal $\left\\{\bm{v}_{t}\right\\}_{t\geq 0}$. Since we
are interested in the temporal characterization of voltage graph signals, we
revise GSP concepts that are applied to time series of graph signals [57, 58]
in this subsection. Then, we utilize these concepts while modeling the
temporal dynamics at generator buses in Section III-B.

In order to characterize graph signal process $\\{{\bm{v}}_{t}\\}_{t\geq 0}$,
a joint time-vertex domain is considered in the literature by defining filters
whose response is shift invariant with respect to the time series shift
operator $z^{-1}$ and an appropriately chosen GSO [59]. To study the same, map
the time series of graph signal $\bm{v}_{t}$ in both the graph frequency (GF)
and $z-$domain by the application of $z$-transform to the GFT of the graph
signal process:

$\displaystyle\textstyle\bm{V}(z)=\sum_{t=0}^{+\infty}\bm{v}_{t}z^{-t},\widetilde{\bm{V}}(z)=\mathbf{U}^{T}\bm{V}(z),$
(9)

We refer to $\widetilde{\bm{V}}(z)$ as the $z$-GFT.  A graph temporal filter’s
[60] impulse response ${\cal H}_{t}(\mathbf{S})$ and output $\bm{v}_{t}$ are

$\displaystyle{\cal
H}_{t}(\mathbf{S})=\sum_{k=0}^{K}h_{k,t}\mathbf{S}^{k},~{}\bm{v}_{t}=\sum_{\tau=0}^{t}{\cal
H}_{t-\tau}(\mathbf{S})\bm{x}_{\tau},$ (10)

respectively. Graph filter output $\bm{v}_{t}$ in the $z$-domain is:

$\displaystyle\\!\\!\\!\\!\bm{V}(z)\\!=\\!{\cal
H}(\mathbf{S}\\!\otimes\\!z)\bm{X}(z),\text{where}~{}{\cal
H}(\mathbf{S}\\!\otimes\\!z)\\!:=\\!\sum_{t=0}^{+\infty}\\!{\cal
H}_{t}(\mathbf{S})z^{-t}$ (11)

when the input is $\bm{x}_{t}$ with $z$-transform ${\bm{X}}(z)$ and ${\cal
H}_{t}(\mathbf{S})$ are matrix polynomials of the GSO operator:

$\textstyle{\cal
H}_{t}(\mathbf{S})=\sum_{k=0}^{K}h_{k,t}\mathbf{S}^{k}~{}~{}\leftrightarrow~{}~{}{\cal
H}(\mathbf{S}\otimes z)=\sum\limits_{k=0}^{K}H_{k}(z)\mathbf{S}^{k}.$ (12)

Here $H_{k}(z)$ is the $z$-transform of the filter $h_{k,t}$. We can define
also the following impulse response in the GF domain:

$\textstyle[\tilde{\bm{h}}_{t}]_{i}=\mathrm{H}_{t}(\lambda_{i}),~{}~{}~{}~{}\mathrm{H}_{t}(\lambda):=\sum_{k=0}^{K}h_{k,t}\lambda^{k}$
(13)

and the graph-temporal joint transfer function in the $z$ and GF domain as:

$\textstyle[\tilde{\bm{h}}(z)]_{i}=\mathbb{H}(\lambda_{i},z),~{}~{}\mathbb{H}(\lambda,z)=\sum_{t=0}^{+\infty}\sum_{k=0}^{K}h_{k,t}\lambda^{k}z^{-t}$
(14)

With that, we obtain following input-output relationship:

$\displaystyle\tilde{\bm{V}}(z)=\texttt{diag}\left(\tilde{\bm{h}}(z)\right)\tilde{\bm{X}}(z),$
(15)

by applying GFT to $z$-domain in (11).

In this work, we focus on a class of graph-temporal filters called GF-ARMA
$(q,r)$ filter [60, 51] . The input-output relation in both time and $z$-GFT
domain are described below, respectively:

$\begin{split}&{\bm{v}}_{t}\\!-\\!{\cal
A}_{1}(\mathbf{S}){\bm{v}}_{t-1}\\!\cdots\\!-\\!{\cal
A}_{q}(\mathbf{S}){\bm{v}}_{t-q}\\!\\!=\\!\\!{\cal
B}_{0}(\mathbf{S}){\bm{x}}_{t}\\!\\!+\cdots\\!+\\!{\cal
B}_{r}(\mathbf{S}){\bm{x}}_{t\\!-\\!r},\\\
&\texttt{diag}\left(\tilde{\bm{a}}(z)\right)\tilde{\bm{V}}(z)=\texttt{diag}\left(\tilde{\bm{b}}(z)\right)\tilde{\bm{X}}(z),\end{split}$

where $\widetilde{\bm{a}}(z)=1-\sum_{t=1}^{q}\widetilde{\bm{a}}_{t}z^{-t}$ and
$\widetilde{\bm{b}}(z)=\sum_{t=0}^{r}\widetilde{\bm{b}}_{t}z^{-t}$ are the
$z$-transform of the graph frequency responses of the graph filter taps
$\\{{\cal A}_{t}(\mathbf{S})\\}_{t=1}^{q}$, $\\{{\cal
B}_{t}(\mathbf{S})\\}_{t=0}^{r}$ for the GF-ARMA $(q,r)$ filter. Particularly,
the GF-AR (2) filter is used in Section III-B to describe generator temporal
dynamics.

###  II-C Measurements and parameters of the electric grid

The electric grid network can be represented by an undirected graph
$\mathcal{G}=(\mathcal{N},\mathcal{E})$ where nodes are buses and its edges
are its transmission lines. The vertex set is a union between set of
generator, $\mathcal{N}_{G}$ and non-generator/ load buses, $\mathcal{N}_{L}$,
$i\in\left\\{\mathcal{N}_{G}\cup\mathcal{N}_{L}\right\\}=\mathcal{N}$ and the
edge set $(i,j)\in\mathcal{E}$ depicts electrical connections. To obtain Ohm’s
law for a network of transmission lines, one starts from the telegrapher
equations for a single line to obtain the so-called ABCD parameters that
relate input-output currents and voltages in the Fourier domain. The equations
are then rearranged and the so-called $\pi$-model is attained, which is an
equivalent circuit containing a series impedance element and parallel
susceptance elements. The $\pi$-model leads to the $2\times 2$ branch
admittance matrix that relates current and voltage injections at the from and
to ends of a transmission line [61]. From the branch admittance matrix of the
network, applying Kirchhoff’s law, one can relate the current and voltage
phasors for the entire network, introducing a system admittance matrix,
$\bm{Y}\in\mathbb{C}^{\lvert\mathcal{N}\rvert}$[61] thus obtaining the network
version of Ohm’s law (see (18)). The matrix $\bm{Y}$ is defined as:

$\displaystyle\left[\bm{Y}\right]_{i,j}=\begin{cases}\sum_{k\in\mathcal{N}_{i}}y_{i,k},~{}i=j\\\
-y_{i,j},~{}i\neq j\end{cases}$ (16)

where $y_{i,j}$ is the admittance of the branch between buses $i$ and $j$ if
$(i,j)\in\mathcal{E}$. The system admittance matrix $\bm{Y}$ is a complex
symmetric matrix and it is equivalent to the complex-valued graph Laplacian
matrix associated with the power grid.  Next, we will partition the nodes or
buses into generator and non-generators, so that:

$\displaystyle\bm{Y}=\begin{bmatrix}\bm{Y}_{gg}&\bm{Y}_{g\ell}\\\
\bm{Y}^{\top}_{g\ell}&\bm{Y}_{\ell\ell}\end{bmatrix},$ (17)

where $\bm{Y}_{gg}$ is the generator buses-only network, $\bm{Y}_{g\ell}$
includes the portion connecting generators and loads and $\bm{Y}_{\ell\ell}$
corresponds to the section of the grid connecting the loads buses among
themselves. The shunt (fixed admittance to ground at a bus) elements at all
generator buses are denoted by
$\bm{y}_{sh}^{g}\in\mathbb{C}^{\lvert\mathcal{N}_{G}\rvert}$ and at all load
buses by $\bm{y}_{sh}^{\ell}\in\mathbb{C}^{\lvert\mathcal{N}_{L}\rvert}$.

The state of the system, from which all other physical quantities of interest
can be derived, are the voltage phasors at each bus. In the following we
assume that a PMU installed on node/bus $i\in\mathcal{N}$ provides a noisy
measurement of voltage and current phasors at time $t$ where $v(t,i)=\lvert
v(t,i)\rvert e^{j\theta(t,i)}$. With some abuse of notation, we will refer to
the PMU data as $v(t,i)$ as well. Let the vector of voltage phasors collected
at time $t$ be $\bm{v}_{t}\in\mathbb{C}^{\lvert\mathcal{N}\rvert}$. After
$\bm{v}_{t}$ is partitioned into voltages at generator and non-generator
buses, let $\bm{i}^{g}_{t}\in\mathbb{C}^{\lvert\mathcal{N}_{G}\rvert}$ be the
generator current and
$\bm{i}^{\ell}_{t}\in\mathbb{C}^{\lvert\mathcal{N}_{L}\rvert}$ the load
current.  Ohm’s law for a network is333Note that the admittances values are
frequency responses evaluated at $60$Hz (for the US) the voltage and current
signals are the corresponding envelopes at the same frequency; hence the
assumption is the voltage and currents are narrowband and the convolution can
be approximated by gain and phase rotation equal to the Fourier response at
$60$ Hz.:

$\displaystyle\left(\\!\bm{Y}\\!+\\!\texttt{diag}\\!\left(\begin{bmatrix}\bm{y}_{sh}^{g}\\\
\bm{y}_{sh}^{\ell}\end{bmatrix}\right)\\!\right)\bm{v}_{t}\\!=\\!\bm{i}_{t},\text{where}~{}\bm{v}_{t}\\!=\\!\begin{bmatrix}\bm{v}^{g}_{t}\\\
\bm{v}^{\ell}_{t}\end{bmatrix},\\!\bm{i}_{t}\\!=\\!\begin{bmatrix}\bm{i}^{g}_{t}\\\
\bm{i}^{\ell}_{t}\end{bmatrix}$ (18)

To describe the operating conditions of the system we introduce a few more
quantities. In power systems transient dynamic analysis the impact of
generating units is modeled as an internal bus characterized by a generator
impedance (or admittance)
$\bm{y}_{g}\in\mathbb{C}^{\lvert\mathcal{N}_{G}\rvert}$ for
$g\in\mathcal{N}_{G}$ connected to an ideal voltage source called internal
voltage[61]; we denote its value at time $t$ by $E(t,i)=\lvert E(t,i)\rvert
e^{\delta_{(}t,i)},i\in\mathcal{N}_{G}$ and the corresponding vector as
$\bm{e}_{t}\in\mathbb{C}^{\lvert\mathcal{N}_{G}\rvert}$ so that
$\left[\bm{e}_{t}\right]_{i}=E(t,i)$. The current at generator bus in (18),
$\bm{i}^{g}_{t}$, is obtained as the multiplication of generator admittance
and the difference in voltage at the internal bus and the generator bus [7] :

$\bm{i}^{g}_{t}=\texttt{diag}(\bm{y}_{g})\left(\bm{e}_{t}-\bm{v}^{g}_{t}\right)$
(19)

As mentioned in Section. I, the generators respond to electric load in the
grid. In order to model the generators response, a commonly used approximation
is that at the load buses $\ell\in\mathcal{N}_{L}$ are slowly varying
admittances [7]. We denote them as $\bm{y}_{\ell}(t)\in\mathbb{C}^{N}$.

##  III Graph Signal Processing for the grid

Having described the relevant GSP concepts and introduced grid quantities and
parameters of interest, we are ready to introduce the Grid-GSP framework444Our
preliminary GSP modeling effort can be found in [1]. Firstly, we define the
GSO for the grid, then support the definition by introducing the graph-filter
model that justifies it, and finally characterize its temporal dynamics. All
of the above yields a GSP generative model for the voltage phasor measurements
as a low-pass GSP model, as detailed next.

###  III-A Grid graph generative model

Grid-GSP for voltage phasors data relies on the following definitions:

######  Definition 3.

The graph shift operator (GSO) is a complex symmetric matrix equal to a
diagonal perturbation of the system admittance matrix with generator
admittance values,

$\displaystyle\mathbf{S}\triangleq\bm{Y}+\texttt{diag}\left(\begin{bmatrix}\bm{y}_{g}+\bm{y}_{sh}^{g}\\\
\bm{y}_{sh}^{\ell}\end{bmatrix}\right)$ (20)

From the definition of the GSO it follows that:

######  Definition 4.

The grid Graph Fourier Transform (GFT) basis for voltage phasors is the
orthogonal matrix $\mathbf{U}$ given by the eigenvalue decomposition of the
GSO in Definition 3:

$\displaystyle\mathbf{S}=\mathbf{U}\bm{\Lambda}\mathbf{U}^{\top},~{}~{}\lvert\lambda_{\text{min}}\rvert>0$
(21)

Here, the GSO $\mathbf{S}$ is a complex-symmetric matrix that has the same
support as the electric-grid graph Laplacian as $\bm{Y}$ with the diagonal
addition of generator admittances. Note that unlike the graph Laplacian, this
GSO is invertible, $\lvert\lambda_{\text{min}}\rvert>0$. Even when shunt
elements $\bm{y}_{sh}^{g},\bm{y}_{sh}^{\ell}$ are ignored as conventionally
done to solve power-flow problems in power systems, a diagonal term with the
generator admittances $\bm{y}_{g}$ that is added to the principal diagonal of
$\bm{Y}$, makes the GSO $\mathbf{S}$ invertible (see (20)).

With the GSO $\mathbf{S}$ is defined as in (20), one can rewrite (18) and
substitute for $\bm{i}^{g}_{t}$ from (19) :

$\displaystyle\left(\\!\bm{Y}\\!+\\!\texttt{diag}\\!\left(\begin{bmatrix}\bm{y}_{sh}^{g}\\\
\bm{y}_{sh}^{\ell}\end{bmatrix}\right)\\!\right)\bm{v}_{t}\\!=\\!\begin{bmatrix}\texttt{diag}(\bm{y}_{g})\left(\bm{e}_{t}-\bm{v}^{g}_{t}\right)\\\
\bm{i}^{\ell}_{t}\end{bmatrix}$ (22)

From now on with slight abuse of notation we denote $\bm{v}_{t}$ as voltage
phasor measurements that are noisy therefore we add measurement noise
$\eta_{t}$ which yields the following equation,

$\displaystyle\boxed{\bm{v}_{t}=\mathcal{H}(\mathbf{S})\begin{bmatrix}\texttt{diag}(\bm{y}_{g})\bm{e}_{t}\\\
\bm{i}^{\ell}_{t}\end{bmatrix}+\bm{\eta}_{t}}$ (23)

The Grid-GSP generative model for voltage phasor measurements $\bm{v}_{t}$ is
given by (23). The linear shift-invariant graph filter is
$\mathcal{H}(\mathbf{S})=\mathbf{S}^{-1}$.

######  Remark 1.

Shift-invariance of $\mathcal{H}(\mathbf{S})=\mathbf{S}^{-1}$ can be directly
verified from (2) i.e.
$\mathcal{H}(\mathbf{S})\mathbf{S}=\mathbf{S}\mathcal{H}(\mathbf{S})=\mathbb{I}$.
Since $\mathcal{H}(\mathbf{S})$ is shift invariant, it can be expressed as a
matrix polynomial in $\mathbf{S}$ (Theorem 1 in [53]). Also,
$\mathcal{H}(\mathbf{S})=\mathbf{S}^{-1}$ can be written as in (3) where
coefficients $h_{k}$ can be determined by the application of Cayley-Hamilton
theorem for inverse matrices [62].

$\mathcal{H}(\mathbf{S})$ is approximately a low-pass graph filter [52] due to
the inversion of GSO since the graph frequency response of the filter can be
written from (7) as
$\texttt{diag}\left(\tilde{\bm{h}}\right)=\bm{\Lambda}^{-1}$. This implies
that as the graph frequency decreases, the magnitude of the filter response
declines. More importantly, since generic power grids tend to be organized as
communities system admittance matrix $\bm{Y}$ tends to be sparse [63].
Therefore the GSO $\mathbf{S}$ has a high condition number and the graph
frequency response of $\mathcal{H}(\mathbf{S})$ is such that it tapers off
after a certain $\lambda_{k}$.

To visualize this more explicitly, consider $\bm{\Lambda}_{\mathcal{K}}$ to be
the diagonal matrix with entries
$\lambda_{i},~{}i\in\mathcal{K}=\\{1,\ldots,k\\}$. Define a low-pass filter
with $k$ frequency components and consequently the voltage phasor measurements
as

$\displaystyle\mathcal{H}_{k}(\mathbf{S})$
$\displaystyle\triangleq\\!{\mathbf{U}}~{}\texttt{diag}\left(\tilde{\bm{h}}_{k}\right){\mathbf{U}}^{\top},\left[\tilde{\bm{h}}_{k}\right]_{i}\\!\\!=\\!\\!\begin{cases}\lambda_{i}^{-1},~{}i\in\mathcal{K}\\\
0,~{}~{}\text{else}\end{cases}$ $\displaystyle{\bm{v}}_{t}$
$\displaystyle\approx\mathcal{H}_{k}(\mathbf{S})\begin{bmatrix}\texttt{diag}(\bm{y}_{g})\bm{e}_{t}\\\
\bm{i}^{\ell}_{t}\end{bmatrix}+{\bm{\eta}}_{t},$ (24)

where $\mathcal{H}_{k}(\mathbf{S})$ will represent the principal subspace of
the voltage phasors whose dimensionality is the number of graph-frequencies
$|\mathcal{K}|$. Therefore (24) defines the low-dimensional generative model
for quasi-steady state voltage phasor measurements. The error term
$\bm{\eta}_{t}$ now also captures modeling approximation.

To provide insights on the temporal dynamics of the voltage phasors, we need
to capture the structure of the excitation term. As a matter of fact,
$\bm{e}_{t}$ and $\bm{i}^{\ell}_{t}$, have different dynamics, as discussed in
the subsequent subsections.

###  III-B A GSP model for generator dynamics: $\bm{e}_{t}$

The excitation term corresponding to generator currents has elements as
$\left[\bm{e}_{t}\right]_{i}$ coming from each generator $i\in G$. We
illustrate a non-linear dynamical model for the generators internal voltages,
namely $\bm{e}_{t}\in\mathbb{C}^{\mathcal{N}_{G}}$ utilizing a GF-AR(2) graph
temporal filter from Section II-B. The model is inspired by the classical
swing equations[8, 9] that describes the coupled dynamics of the generators
phase angles, ${\delta}_{i}(t),~{}i\in G$ and the resulting variation in
frequency, ${\omega}_{i}(t)\triangleq\dot{\delta_{i}}-\omega_{0}$ where
${\omega}_{0}=2\pi f_{0}$ with $f_{0}$ being the grid frequency ($50$ or $60$
Hz).

Our model, relies on two steps. First, we model the dynamics of a signal
$\bm{x}_{t}$ obtained through the following non-linear transformation of the
internal generator voltages:

$\displaystyle\bm{x}_{t}$
$\displaystyle\triangleq(\texttt{diag}(\bm{m}))^{\frac{1}{2}}\ln(\bm{e}_{t})$
(25) $\displaystyle\bm{\delta}_{t}$
$\displaystyle=(\texttt{diag}(\bm{m}))^{-\frac{1}{2}}\Im\\{\bm{x}_{t}\\},\lvert\bm{e}\rvert_{t}\\!=\\!(\texttt{diag}(\bm{m}))^{-\frac{1}{2}}\Re\\{\bm{x}_{t}\\}$

where the vector $\bm{m}$ entries are the so-called generators masses,
$\bm{\delta}_{t}$ are the generators angles that appear in the swing equations
and $\lvert\bm{e}\rvert_{t}$ are internal generator voltage magnitudes.
Second, like in the swing equations, to describe the generators interactions,
we resort to a Kron-reduction [8, 9] of the network, in which generators are
all adjacent. To define this generator-only network and the corresponding GSO,
consider the following admittance matrix, $\bm{Y}_{\text{all}}$ that describes
the network topology consisting of the generator internal buses, generator
buses and non-generator buses like done in [9]:

$\displaystyle\bm{Y}_{\text{all}}\\!=\\!\begin{bmatrix}\texttt{diag}\left(\bm{y}_{g}\\!+\\!\bm{y}_{sh}^{g}\right)&-\begin{bmatrix}\\!\texttt{diag}\left(\bm{y}_{g}\\!+\\!\bm{y}_{sh}^{g}\right)\\!&\\!\bm{0}\end{bmatrix}\\\
-\begin{bmatrix}\\!\texttt{diag}\left(\bm{y}_{g}\\!+\\!\bm{y}_{sh}^{g}\right)\\!&\\!\bm{0}\end{bmatrix}^{\top}&\mathbf{S}_{\text{ci}}\\\
\end{bmatrix}$
$\displaystyle\mathbf{S}_{\text{ci}}\triangleq\mathbf{S}+\\!\texttt{diag}\left(\begin{bmatrix}\bm{0}&\bm{y}_{\ell}^{\top}\end{bmatrix}\right)$

In order to model $\bm{Y}_{\text{all}}$, it is assumed that the loads are
varying very slowly in time i.e.
$\bm{y}_{\ell}(t)\approx\bm{y}_{\ell}~{}~{}\forall t$. Then, let us denote by
$\mathfrak{Sh}(\bm{A},\bm{B})$ the Schur complement of block $\bm{B}$ of
matrix $\bm{A}$. We compute the Schur complement of block
$\mathbf{S}_{\text{ci}}$ of the matrix $\bm{Y}_{\text{all}}$ which is nothing
but Kron reduction. The Schur complement of the $\mathbf{S}_{\text{ci}}$ in
$\bm{Y}_{\text{all}}$ has two contributions:

$\mathfrak{Sh}(\bm{Y}_{\text{all}},\mathbf{S}_{\text{ci}})=j\bm{Y}_{\text{red}}+{\bm{E}}_{\text{red}}$
(26)

where ${\bm{E}}_{\text{red}}$ is a real diagonal dominated matrix, and the
imaginary part $\bm{Y}_{\text{red}}$ has the structure of a graph Laplacian.
The proposed dynamical model for the graph signal $\bm{x}_{t}$ relies on the
following definition for the GSO of the Kron-reduced generator-only graph:

######  Definition 5.

A GSO is defined for the Kron-reduced generator only network as

$\mathbf{S}_{\text{red}}=(\texttt{diag}(\bm{m}))^{-\frac{1}{2}}\bm{Y}_{\text{red}}(\texttt{diag}(\bm{m}))^{-\frac{1}{2}}\in\mathbb{R}^{\lvert\mathcal{N}_{G}\rvert}$
(27)

with the following eigenvalue decomposition,

$\displaystyle\mathbf{S}_{\text{red}}=\mathbf{U}_{\text{red}}\bm{\Lambda}_{\text{red}}\mathbf{U}_{\text{red}}^{\top}$
(28)

and the orthonormal GFT basis being $\mathbf{U}_{\text{red}}$.

We introduce the GSP based dynamical model for the complex-valued generator
internal voltages $\bm{e}_{t}$ via graph temporal filter GF-AR (2) as follows,

GSP-based dynamics for generator internal voltages

$\bm{e}_{t}=\exp\left((\texttt{diag}(\bm{m}))^{-\frac{1}{2}}\bm{x}_{t}\right)$
(29)

where $\bm{x}_{t}$ is a GF-AR (2) process, i.e. the $z$-GFT
$\tilde{\bm{X}}(z)$ satisfies the following:

$\displaystyle\texttt{diag}\left(\tilde{\bm{a}}(z)\right)\tilde{\bm{X}}(z)$
$\displaystyle=\tilde{\bm{W}}(z)$ (30) $\displaystyle\tilde{\bm{a}}(z)$
$\displaystyle=1-\tilde{\bm{a}}_{1}z^{-1}-\tilde{\bm{a}}_{2}z^{-2}$ (31)

The GSP based dynamical model takes inspiration from swing equation for
generator angles and we empirically choose to model generator internal voltage
magnitude $\lvert\bm{e}\rvert_{t}$ also using a GF-AR (2) model although in
most power system models, the dynamics of the amplitudes of the generators
internal voltages are typically ignored. The swing equation for the generators
angles [64] are a key tool for power systems dynamical analysis:

$\displaystyle\texttt{diag}\left(\bm{m}\right)\ddot{\bm{\delta}}+\texttt{diag}\left(\bm{d}\right)\dot{\bm{\delta}}=\overline{\bm{w}}-\bm{Y}_{\text{red}}\bm{\delta},$
(32)

where $\bm{m}$ are the generators masses, introduced previously, $\bm{d}$ are
the damping coefficients of generators (often neglected) and
$\overline{\bm{w}}-\bm{Y}_{\text{red}}\bm{\delta}$ is the imbalance between
the electrical and mechanical power that triggers the change in generator
angular velocity and acceleration. Note that
$\Im\\{\bm{x}\\}=\left(\texttt{diag}\left(\bm{m}\right)\right)^{\frac{1}{2}}\bm{\delta}$.
We can manipulate (32) to prove the following:

######  Proposition 1.

Let $\bm{w}$ be such that
$\Im\\{\bm{w}\\}=(\texttt{diag}(\bm{m}))^{-\frac{1}{2}}\overline{\bm{w}}$.
Using the approximation:

$\texttt{diag}\left(\bm{d}\right)\texttt{diag}\left(\bm{m}\right)^{-1}\approx\chi\mathbb{I},$

i.e. the homogeneous simplification as in [64, 65], the dynamics of
$\Im\\{\bm{x}_{t}\\}$ can be justified with an GF-AR (2) model with GSO
$\mathbf{S}_{\text{red}}$:

$\displaystyle\Im\\{\bm{x}_{t}\\}-{\cal
A}_{1}(\mathbf{S}_{\text{red}})\Im\\{\bm{x}_{t-1}\\}-{\cal
A}_{2}(\mathbf{S}_{\text{red}})\Im\\{\bm{x}_{t-2}\\}\approx\Im\\{\bm{w}\\},$
$\displaystyle{\cal
A}_{1}(\mathbf{S}_{\text{red}}):=(2\\!-\\!\chi)\mathbb{I}-\mathbf{S}_{\text{red}},~{}~{}{\cal
A}_{2}(\mathbf{S}_{\text{red}}):=(\chi\\!-\\!1)\mathbb{I}$ (33)

###### Proof.

Simple algebra on (32) allows us to recast the equations in the following
form:

$\displaystyle\Im\\{\ddot{\bm{x}}\\}+\chi\Im\\{\dot{\bm{x}}\\}=\texttt{diag}\left(\bm{m}\right)^{-\frac{1}{2}}\overline{\bm{w}}-\mathbf{S}_{\text{red}}\Im\\{\bm{x}\\}.$
(34)

Assuming that the sampling rate is fast enough, and normalizing it to 1, the
finite difference approximations for the derivatives are
$\dot{\bm{x}}\approx\bm{x}_{t}-\bm{x}_{t-1},~{}~{}\ddot{\bm{x}}\approx\bm{x}_{t+1}-2\bm{x}_{t}+\bm{x}_{t-1}$
and can be used to obtain AR (2) GF equations for the samples
$\Im\\{\bm{x}_{t}\\}$ in (33). ∎

The model that we introduce is simply extending the GF-AR (2) model to capture
both the real and imaginary part of $\bm{x}_{t}$ i.e. internal generator
voltage magnitudes$\lvert\bm{e}\rvert_{t}$ and angles $\bm{\delta}_{t}$
respectively and suggesting to search the $2|{\cal N}_{\cal G}|$ parameters to
fit the model with $\tilde{\bm{a}}_{1},\tilde{\bm{a}}_{2}$ rather than
exploring a general MIMO filter response. For simplicity of representation, we
write the dynamical equation for $\bm{x}_{t}$ in the GF domain,

$\displaystyle\tilde{\bm{x}}_{t}=\texttt{diag}\left(\tilde{\bm{a}}_{1}\right)\tilde{\bm{x}}_{t-1}+\texttt{diag}\left(\tilde{\bm{a}}_{2}\right)\tilde{\bm{x}}_{t-2}+\tilde{\bm{w}}_{t}$
(35)

such that the impulse response of the filter at graph-frequency
$\lambda_{\text{red},i}$ is defined by
$\left[\tilde{\bm{a}}_{1}\right]_{i},\left[\tilde{\bm{a}}_{2}\right]_{i}$.

Figure 1: Block diagram showing generative model for voltage phasor
measurements.

Note from (13) that
$\left[\tilde{\bm{a}}_{1}\right]_{i},\left[\tilde{\bm{a}}_{2}\right]_{i}$ can
be written as polynomials in graph frequency $\lambda_{\text{red},i}$,

$\displaystyle\textstyle\left[\tilde{\bm{a}}_{1}\right]_{i}=\sum_{k=0}^{K_{1}-1}a_{k,1}\lambda_{\text{red},i}^{k},~{}\left[\tilde{\bm{a}}_{2}\right]_{i}=\sum_{k=0}^{K_{2}-1}a_{k,2}\lambda_{\text{red},i}^{k},$
(36)

which characterizes the poles of the AR system using graph frequencies. From
(33), a first order polynomial may suffice to characterize
$\left[\tilde{\bm{a}}_{1}\right]_{i},\left[\tilde{\bm{a}}_{2}\right]_{i}$ in
most cases.

###  III-C Load dynamics: $\bm{i}^{\ell}_{t}$

There are several papers in the literature that deal with load forecasting and
modeling [66]. We adopt a simple AR-2 model per node or load bus to describe
the dynamics of the load,

$\displaystyle\bm{i}^{\ell}_{t}=\texttt{diag}\left(\bm{b}_{1}\right)\bm{i}^{\ell}_{t-1}+\texttt{diag}\left(\bm{b}_{2}\right)\bm{i}^{\ell}_{t-2}+\bm{\epsilon}_{t}$
(37)

where parameters $\bm{b}_{1},\bm{b}_{2}$ are estimated load data time series.
The block diagram in Fig.1 summarize our modeling efforts.

The unique nature of voltage phasor measurements allows us to describe a
similar model for any subset of measurements on a graph. This is discussed
next.

###  III-D Low-pass property of down-sampled voltage graph signal

Let $\bm{v}_{\mathcal{M}}$ (time index $t$ is ignored for simplicity) be the
down-sampled voltage graph signal where $\mathcal{M}\subset\mathcal{N}$ is the
set of node indices at which measurements are available.  It can be shown that
any down-sampled graph signal with arbitrary graph frequency response is low-
pass in the reduced-graph frequency domain. It suggests that one can utilize
all the methods for low-pass graph signals onto down-sampled versions of the
graph signal as well. This is summarized in lemma 1 below.

######  Lemma 1.

Let $\bm{v}_{\mathcal{M}}$ be any graph signal down-sampled in the vertex-
domain with $\lvert\mathcal{M}\rvert$ samples. Let the GSO defined with
respect to the full graph $\mathbf{S}$ be invertible. Then, with the GSO
defined with respect to the reduced-graph of $\mathcal{M}$ vertices as
$\mathbf{S}_{\text{red},\mathcal{M}}$, graph signal $\bm{v}_{\mathcal{M}}$ is
the output of low-pass graph filter
$\mathcal{H}\left(\mathbf{S}_{\text{red},\mathcal{M}}\right)\triangleq\mathbf{S}_{\text{red},\mathcal{M}}^{-1}$

$\displaystyle\bm{v}_{\mathcal{M}}=\mathcal{H}\left(\mathbf{S}_{\text{red},\mathcal{M}}\right)\bm{\varphi}$
(38)

where the GSO for the reduced-graph is given by Kron-reduction of
$\mathbf{S}$,
$\mathbf{S}_{\text{red},\mathcal{M}}=\mathfrak{Sh}\left(\mathbf{S},\mathbf{S}_{\mathcal{M}^{c}\mathcal{M}^{c}}\right)$.

###### Proof.

Consider a graph signal $\bm{v}$ with arbitrary graph frequency response with
respect to GSO $\mathbf{S}$,

$\displaystyle\bm{v}=\mathcal{H}\left(\mathbf{S}\right)\bm{x}=\mathbf{S}^{-1}\left(\mathbf{S}\mathcal{H}\left(\mathbf{S}\right)\bm{x}\right)$
(39)

The GSO $\mathbf{S}$ is rewritten in a $2\times 2$ block form

$\displaystyle\mathbf{S}=\begin{bmatrix}\mathbf{S}_{\mathcal{M}\mathcal{M}}&\mathbf{S}_{\mathcal{M}\mathcal{M}^{c}}\\\
\mathbf{S}_{\mathcal{M}\mathcal{M}^{c}}^{\top}&\mathbf{S}_{\mathcal{M}^{c}\mathcal{M}^{c}}\end{bmatrix},$
(40)

and $\mathbf{S}^{-1}$ can be written using inverse formula for block matrices.
When graph signal $\bm{v}$ is down-sampled, only $\mathcal{M}$ rows are
considered on both sides of (39). Thus we have,

$\displaystyle\bm{v}_{\mathcal{M}}$
$\displaystyle=\mathbf{S}_{\text{red},\mathcal{M}}^{-1}\overbrace{\begin{bmatrix}\mathbb{I}_{\lvert\mathcal{M}\rvert}&-\mathbf{S}_{\mathcal{M}\mathcal{M}^{c}}\mathbf{S}_{\mathcal{M}^{c}\mathcal{M}^{c}}^{-1}\end{bmatrix}\left(\mathbf{S}\mathcal{H}\left(\mathbf{S}\right)\bm{x}\right)}^{\bm{\varphi}}$
(41)

where $\mathbf{S}_{\text{red},\mathcal{M}}$ is the Schur complement of the
block $\mathbf{S}_{\mathcal{M}^{c}\mathcal{M}^{c}}$ in the GSO $\mathbf{S}$
i.e.,

$\displaystyle\\!\\!\\!\mathbf{S}_{\text{red},\mathcal{M}}\\!=\\!\mathfrak{Sh}\\!\left(\mathbf{S},\\!\mathbf{S}_{\mathcal{M}^{c}\mathcal{M}^{c}}\\!\right)\\!=\\!\mathbf{S}_{\mathcal{M}\mathcal{M}}\\!-\\!\mathbf{S}_{\mathcal{M}\mathcal{M}^{c}}\mathbf{S}_{\mathcal{M}^{c}\mathcal{M}^{c}}^{-1}\mathbf{S}_{\mathcal{M}\mathcal{M}^{c}}^{\top}~{}~{}~{}~{}\qed$
(42)

Lemma 1 translates to an interesting self-similarity/fractional property for
voltage graph signals in that the down-sampled version $\bm{v}_{\mathcal{M}}$
is still a low-pass graph signal. The self-similarity is due to
$\mathcal{H}\left(\mathbf{S}\right)=\mathbb{I}$. In summary, for voltage graph
signals,

$\displaystyle\bm{v}=\mathbf{S}^{-1}\bm{i},~{}~{}\bm{v}_{\mathcal{M}}=\mathbf{S}_{\text{red},\mathcal{M}}^{-1}\left(\bm{i}_{\mathcal{M}}\\!-\\!\mathbf{S}_{\mathcal{M}\mathcal{M}^{c}}\mathbf{S}_{\mathcal{M}^{c}\mathcal{M}^{c}}^{-1}\bm{i}_{\mathcal{M}^{c}}\right)$
(43)

In the power grid, this property has been illustrated empirically in several
papers [68, 12] that highlight low-dimensionality of measurements from a
subset of buses. Although the reduced-graph is denser compared to the original
graph, it still helps to infer faults or events that occurred in a subset of
nodes where sensors are not installed as long as correct placement strategies
are devised i.e. that of choosing the subset $\mathcal{M}$. Work in [26]
explored the optimal placement for fault localization in the under-sampled
regime and also made connections with GSP theory.

##  IV Revisiting algorithms from GSP for PMU data

In this section we study some of the implications Grid-GSP has while
understanding sampling, optimal placement of measurement devices in power
systems, interpolation of missing samples and network inference. The
underlying generative model responsible for low-rank nature of data that has
been established in the previous section helps explaining the success that
many past works, such as [12, 69, 70], have attained in recovering missing PMU
data using matrix completion methods. The low-pass nature of the voltage graph
signals discussed in Section III provides the theoretical underpinning that
support the arguments made in the literature.

###  IV-A Sampling and recovery of grid-graph signals

From the approximation in (24) we see that voltage graph signals have graph
frequency content that drops as $\lambda_{k}$ grows. This characteristic
renders the signal approximately band-limited in the GFT domain [71] which
means that there is a cut-off frequency $\lambda_{k}$ such that frequency
content corresponding to $\lambda_{k+1}$ and higher is negligible. Let the GFT
basis corresponding to the first dominant $k$ graph frequecies be
$\mathbf{U}_{\mathcal{K}}$. The bandlimiting operator is,
$\bm{\mathcal{B}}_{\mathcal{K}}=\mathbf{U}_{\mathcal{K}}\mathbf{U}^{\top}_{\mathcal{K}}\in\mathbb{C}^{\lvert\mathcal{N}\rvert\times\mathcal{K}}$
and the low frequency component of $\bm{v}_{t}$ is:

$\displaystyle\bm{\mathcal{B}}_{\mathcal{K}}\bm{v}_{t}=\mathbf{U}_{\mathcal{K}}\mathbf{U}^{\top}_{\mathcal{K}}\bm{v}_{t}$
(44)

Similarly, a vertex limiting operator (with $\lvert\mathcal{M}\rvert$)
vertices is
$\bm{\mathcal{D}}_{\mathcal{M}}=\bm{P}_{{\mathcal{M}}}\bm{P}^{\top}_{{\mathcal{M}}}$
where $\bm{P}_{{\mathcal{M}}}$ has columns that are coordinate vectors such
that each column chooses a vertex/node. When the voltage measurements on the
electrical network are from a few nodes, $i\in\mathcal{M}$ at time $t$, it can
be written as
$\left[\bm{v}_{t}\right]_{\mathcal{M}}=\bm{P}^{\top}_{{\mathcal{M}}}\bm{v}_{t}$.
For reconstruction, results in [71] dictate the necessary condition be that
$\lvert\mathcal{M}\rvert\geq\lvert\mathcal{K}\rvert$. In the presence of
modeling error relative to the perfect band-limited definition, optimal
sampling pattern i.e. the best placement for PMUs on the grid to minimize the
worst-case reconstruction error is closely tied to the grid topology and the
model mismatch relative to a strictly band-limited graph signal [72]. An
optimal placement strategy of PMUs that minimizes the worst-case
reconstruction error in the presence of model mismatch due to imperfect band-
limited nature of the voltage graph signal, also known as the E-optimal
design[72], is sought by maximizing the smallest singular value,
$\sigma_{\min}\left(\bm{\mathcal{D}}_{\mathcal{M}}\mathbf{U}_{\mathcal{K}}\right)$,
i.e. choose rows of $\mathbf{U}_{\mathcal{K}}$ such that they are as
uncorrelated as possible and the resulting matrix has the highest condition
number [71, 72]. Consider then the spatial sampling mask
$\bm{\mathcal{D}}_{\mathcal{M}}=\texttt{diag}(\bm{1}_{\mathcal{M}})$ that
selects $M$ locations.

####  IV-A1 Sampling

The optimal placement of $M$ PMUs maximizes
$\sigma_{\min}(\bm{\mathcal{D}}_{\mathcal{M}}{\mathbf{U}}_{\mathcal{K}})$
which amounts to choosing the rows of $\mathbf{U}_{\mathcal{K}}$ with the
smallest possible coherence (as close as possible to being orthogonal). In
[71] and references therein, a greedy method is employed to find $M$ rows from
${\mathbf{U}}_{\mathcal{K}}$ so that the least singular value is maximized.

Power systems topologies exhibit naturally a community structure that is
reflected in the system admittance matrix $\bm{Y}$ [1] due to population
density or clusters of loads. It is known that a method to determine $k$
communities in a graph is to minimize the Ratio Cut [73] and spectral
clustering performs a relaxed Ratio Cut minimization via $k-$means algorithm
on rows of the eigenvectors $\mathbf{U}_{\mathcal{K}}$ [74]. Thus, choosing
rows of $\mathbf{U}_{\mathcal{K}}$ to be uncorrelated is intuitively putting
PMUs in different graph-clusters or communities. This fact was also discussed
in [26] in the context of sensor placements for fault localization. The PMUs
sampling rate in time exceeds the needs for reconstructions in a quasi-steady
state conditions by a significant margin and it is designed to help detect
sharp transients in the system.

####  IV-A2 Reconstruction

Voltage data samples are obtained down-sampling in space after the optimal
placement of PMUs and also uniformly down-sampling in time. At time $t$ when
$\lvert\mathcal{M}\rvert$ samples, $\left[\bm{v}_{t}\right]_{\mathcal{M}}$ are
available, the following model applies

$\displaystyle\left[\bm{v}_{t}\right]_{\mathcal{M}}\approx\bm{P}_{{\mathcal{M}}}^{\top}\mathbf{U}_{\mathcal{K}}\tilde{\bm{v}}_{t}$
(45)

where $\tilde{\bm{v}}_{t}$ is the GFT of graph signal $\bm{v}_{t}$ Therefore,
reconstruction in spatial domain is done via GFT basis as

$\displaystyle\hat{\bm{v}}_{t}=\mathbf{U}_{\mathcal{K}}\left(\bm{P}_{{\mathcal{M}}}^{\top}\mathbf{U}_{\mathcal{K}}\right)^{\dagger}\left[\bm{v}_{t}\right]_{\mathcal{M}}$
(46)

Reconstruction in temporal domain can be done independently by up-sampling,
i.e. via the windowed inverse Fourier transform of the up-sampled signal
created from uniformly time-decimated data.

###  IV-B Interpolation of missing samples

When voltage measurements are missing or corrupted, denoising and
interpolation of such data can be cast as a graph signal recovery problem by
regularizing the total variation, (TV). Overall, the problem resembles time-
vertex graph signal recovery [59]. Let
$\mathbf{V}=\begin{bmatrix}\bm{v}_{1}&\bm{v}_{2}&\dots&\bm{v}_{T}\end{bmatrix}$
represent the voltage phasor measurements matrix collected over $T$ time
instants. Let $\mathbb{P}_{\Omega}\left(\hat{\mathbf{V}}\right)$ be the set of
available measurements that have samples in entries of set $\Omega$ and are
noisy,

$\displaystyle\min_{\mathbf{V}}~{}~{}$
$\displaystyle\textstyle\lVert\mathbb{P}_{\Omega}\left(\hat{\mathbf{V}}\\!-\\!\mathbf{V}\right)\rVert_{F}^{2}\\!+\\!c_{g}\sum_{t=1}^{T}\lVert\mathbf{S}\bm{v}_{t}\rVert_{1}$
$\displaystyle+\textstyle
c_{t}\sum_{t=2}^{T}\lVert\bm{v}_{t}-\bm{v}_{t-1}\rVert_{2}^{2}$ (47)

where the two regularizing terms measure the variation in the graph and time
domain and $c_{g},c_{t}$ are the corresponding regularization constants.
Importantly, one can use the GSO of the reduced graph,
$\mathbf{S}_{\text{red},\mathcal{M}}$ if we only have access to a subset of
measurements on the grid, $\mathcal{M}$ and employ the same formulation as in
(47) for interpolation of missing samples.

###  IV-C Network inference as graph Laplacian learning

The problem of estimation of GSO $\mathbf{S}$ from voltage phasor measurements
can be cast as a solving a problem similar to graph Laplacian learning [75]
which seeks the GSO that minimizes the total variation of the observed voltage
phasors. If current measurements $\bm{i}_{t}$ are available, then another
regularization term $\lVert\mathbf{S}\bm{v}_{t}-\bm{i}_{t}\rVert_{2}^{2}$ can
be added such that Ohm’s law is satisfied. Therefore, estimation of GSO can be
accomplished by solving the following problem:

$\displaystyle\min_{\mathbf{S}}~{}~{}~{}~{}\textstyle\sum_{t=1}^{T}\lVert\mathbf{S}\bm{v}_{t}\rVert_{1}+\gamma\lVert\mathbf{S}-\texttt{diag}\left(\texttt{Diag}\left(\mathbf{S}\right)\right)\rVert_{F}^{2}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}+\sum_{t=1}^{T}\lVert\mathbf{S}\bm{v}_{t}-\bm{i}_{t}\rVert_{2}^{2}$
(48) $\displaystyle\text{subject
to}~{}~{}\Re{\left[\mathbf{S}\right]_{i,j}}=\Re{\left[\mathbf{S}\right]_{j,i}},\Im{\left[\mathbf{S}\right]_{i,j}}=\Im{\left[\mathbf{S}\right]_{j,i}},i\neq
j,$ (49)
$\displaystyle\Re{\text{Tr}\left(\mathbf{S}\right)}=\alpha\lvert\mathcal{N}\rvert,\Im{\text{Tr}\left(\mathbf{S}\right)}=\beta\lvert\mathcal{N}\rvert$
(50)

Additional constraints on the GSO can be imposed based on the properties of
complex-symmetry (see (49)), sparse off-diagonal entries via the term
$\lVert\mathbf{S}-\texttt{diag}\left(\texttt{Diag}\left(\mathbf{S}\right)\right)\rVert_{F}^{2}$
and dominant diagonal values (see (50)). Also, $\mathbf{S}$ tends to have
larger imaginary values than real especially on the diagonal. and
$\alpha,\beta>1$ control the amplitude of real and imaginary values on the
diagonal. As before, the problem above can be recast with down-sampled voltage
graph signals to infer the Kron-reduced GSO
$\mathbf{S}_{\text{red},\mathcal{M}}$ with the approximation that the term
$\mathbf{S}_{\mathcal{M}\mathcal{M}^{c}}\mathbf{S}_{\mathcal{M}^{c}\mathcal{M}^{c}}^{-1}\bm{i}_{\mathcal{M}^{c}}$
in (41) is treated as additive Gaussian noise. Simulation results for network
inference can be found in Section VI.

##  V Applications of Grid-GSP

The goal of this section is to showcase the benefits of casting problems in
the Grid-GSP framework through two exemplary applications, namely anomaly
detection and data compression. The common thread between them is the use of
the Grid-GFT as a tool to extract informative features from PMU data.

###  V-A Detection of FDI attacks on PMU measurements

This application is based on our preliminary work in [2]. Note that, even
though we cast the problem as that of FDI attacks detection, the idea can be
easily extended to unveil sudden changes due to physical events (like fault-
currents, or topology changes) that similarly excite high GF content. We
assume that we have access to PMU measurements of voltage and current from the
buses they are installed on. Let $\mathcal{A}$ be the set of available
measurements where PMUs are installed and $\mathcal{U}$ be set of unavailable
ones. A measurement model can be written using ‘state’ to be the voltage as

$\displaystyle\underbrace{\begin{bmatrix}\hat{\bm{i}}_{\mathcal{A}}\\\
\hat{\bm{v}}_{\mathcal{A}}\end{bmatrix}}_{\bm{z}_{t}}=\underbrace{\begin{bmatrix}\bm{Y}_{\mathcal{A}\mathcal{A}}&\bm{Y}_{\mathcal{A}\mathcal{U}}\\\
\mathbb{I}_{\lvert\mathcal{A}\rvert}&\bm{0}\end{bmatrix}}_{\bm{H}}\underbrace{\begin{bmatrix}\bm{v}_{\mathcal{A}}\\\
\bm{v}_{\mathcal{U}}\end{bmatrix}}_{\bm{v}}+\bm{\varepsilon}$ (51)

The attacker follows the strategy of FDI attack to manipulate both current and
voltage on the set of malicious buses, $i\in\mathcal{C}\subset\mathcal{A}$ by
introducing a perturbation

$\displaystyle\delta\bm{v}^{T}_{t}=\begin{bmatrix}\delta\bm{v}_{\mathcal{C}}^{T}&\bm{0}_{\lvert\mathcal{P}\rvert+\lvert\mathcal{U}\rvert}^{T}\end{bmatrix},~{}\text{such
that}~{}\bm{Y}_{\mathcal{P}\mathcal{C}}\delta\bm{v}_{\mathcal{C}}=\bm{0}$ (52)

where $\mathcal{P}$ is the set of honest nodes. This requires special
conditions and placement, since $\bm{Y}_{\mathcal{P}\mathcal{C}}$ is tall.
Nonetheless, since the system admittance matrix $\bm{Y}$ is generally sparse
[63], $\bm{Y}_{\mathcal{P}\mathcal{C}}$ does not have full column-rank for a
sufficient number of attackers $\mathcal{C}$ even when all the measurements
are available with $\mathcal{A}=\mathcal{N}$. Our detection problem entails
deciding between the hypotheses of attack $\textit{H}_{1}$ and no attack
$\textit{H}_{0}$. To this end, we can leverage the low-dimensional generative
model for the voltage graph signal that comes from (24), which imposes
additional constraint on the perturbation along with that in (52). In short,
for the attacker to be successful and undetected, she needs to have knowledge
of system parameters and the graph filter with $k$ frequency components
$\mathcal{H}_{k}(\mathbf{S})$. However, since the attacker does not have all
this knowledge, a typical FDI attack as studied in literature is launched
using (52). Using the generative model in (24), we know that under normal
operating conditions in quasi-steady state, the received data $\bm{z}_{t}$
under the no-attack and attack hypotheses $\textit{H}_{0},\textit{H}_{1}$
respectively have the structure:

$\displaystyle\bm{z}_{t}\\!=\\!\\!\begin{cases}\\!\textit{H}_{0}:\\!\\!&\bm{H}\mathcal{H}_{k}(\mathbf{S})\\!\\!\begin{bmatrix}\left(\texttt{diag}(\bm{y}_{g})\bm{e}_{t}\right)^{\top}&\left(\bm{i}^{\ell}_{t}\right)^{\top}\end{bmatrix}^{\top}\\!\\!+\\!\bm{\varepsilon}_{t}\\\
\\!\textit{H}_{1}:\\!\\!&\bm{H}\mathcal{H}_{k}(\mathbf{S})\\!\\!\begin{bmatrix}\left(\texttt{diag}(\bm{y}_{g})\bm{e}_{t}\right)^{\top}&\left(\bm{i}^{\ell}_{t}\right)^{\top}\end{bmatrix}^{\top}\\!\\!\\!\\!+\\!\\!\bm{H}\delta\bm{v}_{t}\\!+\\!\bm{\varepsilon}_{t}\end{cases}$
(53)

Therefore, we project $\bm{z}_{t}$ onto the subspace orthogonal to columnspace
of $\bm{H}\mathcal{H}_{k}(\mathbf{S})$ to get a test statistic, $d({\bm{z}})$.
The projector is:

$\displaystyle\bm{\Pi}_{\perp\bm{H}\mathcal{H}_{k}(\mathbf{S})}\triangleq\mathbb{I}-\left[\bm{H}\mathcal{H}_{k}(\mathbf{S})\right]\left[\bm{H}\mathcal{H}_{k}(\mathbf{S})\right]^{\dagger}$
(54)

and under the no attack hypothesis $\textit{H}_{0}$, energy in the orthogonal
subspace is less than when there is an attack, $\textit{H}_{1}$. This can be
converted to the following test,

$\displaystyle
d(\bm{z}_{t})\triangleq\lVert\bm{\Pi}_{\perp\bm{H}\mathcal{H}_{k}(\mathbf{S})}\bm{z}_{t}\rVert_{2}^{2}~{}~{}\mathop{\gtreqless}_{\textit{H}_{0}}^{\textit{H}_{1}}~{}~{}\tau$
(55)

where $\tau$ is a threshold that can be chosen based on an empirical receiver
operator characteristics (ROC) curve. Note that, since
$\mathcal{H}_{k}(\mathbf{S})$ is a low pass filter, the projector
$\bm{\Pi}_{\perp\bm{H}\mathcal{H}_{k}(\mathbf{S})}$ in (54) is filtering high
graph frequencies and the detection measures the energy on such frequencies as
a signature for anomalies.

Isolation of compromised buses or estimate of $\delta{\bm{v}}$ can also
undertaken with a similar logic. Firstly, using the assumptions in the
previous section we can solve the following regression problem to recover
$\delta\bm{v}_{t}$, formulating a LASSO relaxation of the sparse support
recovery problem:

$\displaystyle\min_{\delta{\bm{v}_{t}}}~{}\lVert\bm{\Pi}_{\perp\bm{H}\mathcal{H}_{k}(\mathbf{S})}({\bm{z}}\\!-\\!\bm{H}\delta{\bm{v}})\rVert^{2}_{2}~{}~{}\text{subject
to}~{}~{}\lVert\delta{\bm{v}_{t}}\rVert_{1}\\!\leq\\!\mu$ (56)

Constraint on the $\ell_{1}$ norm is used to incorporate the prior knowledge
that the attacker has access to a few measurement buses,
$\mathcal{C}\ll\mathcal{N}$. Note that the performance of the algorithm is
also dependent on the number of graph-frequency components i.e. $k$
considered.

###  V-B Compression of PMU measurements

The proposed compression algorithm leverages both (23) and (35). The measure
of distortion we use is the mean-squared error (MSE):

$\displaystyle\textstyle
d(\bm{v},\hat{\bm{v}})\triangleq(\lvert\mathcal{N}\rvert
T)^{-1}\sum_{t=0}^{T}\lVert\bm{v}_{t}-\hat{\bm{v}}_{t}\rVert_{2}^{2}$ (57)

where $T$ denotes the time instant at which samples are stopped collecting.

Algorithm 1 Encoding algorithm for compression

1:${\tilde{\bm{x}}}_{0},{\tilde{\bm{x}}}_{1},\bm{i}_{0}^{\ell},\bm{i}_{1}^{\ell},\left\\{\bm{v}_{t}\right\\}_{t=2}^{T}$

2:for $t=2:T$ do

3: States $\hat{\tilde{\bm{x}}}_{t},\hat{\bm{i}}_{t}^{\ell}$ from (60) and
(61) respectively.

4: Voltage estimate:
$\displaystyle\hat{\bm{v}}^{0}_{t}\\!=\\!\mathcal{H}(\mathbf{S})\\!\begin{bmatrix}\texttt{diag}\\!\left(\bm{y}_{g}\right)\exp\left\\{\\!\texttt{diag}(\bm{m})^{-\frac{1}{2}}\mathbf{U}_{\text{red}}\hat{\tilde{\bm{x}}}_{t}\\!\right\\}\\\
\hat{\bm{i}}^{\ell}_{t}\end{bmatrix}$ (58)

5: Compute GFT of modeling error:
$\tilde{\bm{\xi}}_{t}=\mathbf{U}^{\top}\left(\bm{v}_{t}-\hat{\bm{v}}^{0}_{t}\right)$

6: Quantize:
$\hat{\tilde{\bm{\xi}}}_{t}=\mathcal{Q}\left\\{\tilde{\bm{\xi}}_{t}\right\\},~{}\hat{\bm{v}}_{t}=\hat{\bm{v}}^{0}_{t}+\mathbf{U}\hat{\tilde{\bm{\xi}}}_{t}$

7: Update states,
$\displaystyle\hat{\tilde{\bm{x}}}_{t}\leftarrow\mathbf{U}_{\text{red}}^{\top}\texttt{diag}(\bm{m})^{\frac{1}{2}}\ln(\hat{\bm{e}}_{t}),\hat{\bm{i}}^{\ell}_{t}\leftarrow\left[\mathbf{S}\hat{\bm{v}}_{t}\right]_{\mathcal{N}_{L}}$
(59)
$\displaystyle\text{where}~{}\hat{\bm{e}}_{t}=\left(\texttt{diag}(\bm{y}_{g})\right)^{-1}\left[\mathbf{S}\hat{\bm{v}}_{t}\right]_{\mathcal{N}_{G}}$

8:end for

9:$\left\\{\hat{\tilde{\bm{\xi}}}_{t}\right\\}_{t=2}^{T}$

Since we have a temporal dynamical model for the evolution of voltage in time,
we use differential encoding [76] to quantize the residuals in both generator
and load dynamics, $\tilde{\bm{w}}_{t}$ and $\bm{\epsilon}_{t}$ respectively.
The voltage at time $t$ is:

$\displaystyle{\bm{v}}_{t}=\mathcal{H}(\mathbf{S})\begin{bmatrix}\texttt{diag}\left(\bm{y}_{g}\right)\exp\left\\{\texttt{diag}(\bm{m})^{-\frac{1}{2}}\mathbf{U}_{\text{red}}\left({\tilde{\bm{x}}}^{0}_{t}+\tilde{\bm{w}}_{t}\right)\right\\}\\\
{\bm{i}}^{\ell,0}_{t}+\bm{\epsilon}_{t}\end{bmatrix}$
$\displaystyle\text{where}~{}~{}~{}~{}\tilde{\bm{x}}_{t}^{0}={\texttt{diag}\left(\bm{a}_{1}\right){\tilde{\bm{x}}}_{t-1}\\!+\texttt{diag}\left(\bm{a}_{2}\right)}{\tilde{\bm{x}}}_{t-2}$
(60)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\bm{i}}^{\ell,0}_{t}=\texttt{diag}\left(\bm{b}_{1}\right)\hat{\bm{i}}_{t-1}^{\ell}\\!+\texttt{diag}\left(\bm{b}_{2}\right)\hat{\bm{i}}_{t-2}^{\ell}$
(61)

Thus, $\bm{v}_{t}$ can be approximated as:

$\displaystyle{\bm{v}}_{t}\\!\\!\approx\\!\mathcal{H}(\mathbf{S})\\!\left(\\!\begin{bmatrix}\\!\texttt{diag}\left(\\!\bm{y}_{g}\\!\right)\exp\left\\{\\!\texttt{diag}(\bm{m})^{-\frac{1}{2}}\mathbf{U}_{\text{red}}{\tilde{\bm{x}}}^{0}_{t}\right\\}\\\
{\bm{i}}^{\ell,0}_{t}\end{bmatrix}\\!\\!\\!+\\!\bm{\xi}_{t}\\!\\!\right)$ (62)

Note that, the vector $\bm{\xi}_{t}$ GFT, $\mathbf{U}^{\top}\bm{\xi}_{t}$ has
energy mostly in lower frequency components and is therefore an appropriate
term to quantize using an optimal rate allocation. Specifically, we allocate
bits to each component by setting a desired level of total distortion,
applying the reverse water-filling result [77] which is optimum for a random
vector $\mathbf{U}^{\top}\bm{\xi}_{t}$ whose entries are circularly symmetric
complex independent Gaussian random variables and then quantize the components
accordingly.555The covariance matrix is not diagonal and ideally one would
first whiten the vector $\mathbf{U}^{\top}\bm{\xi}_{t}$ and then quantize the
individual components with bit-allocation akin to reverse water-filling. Since
the statistics of $\mathbf{U}^{\top}\bm{\xi}_{t}$ are time-varying, one has to
perform the whitening transform at each time instant which is a cumbersome
operation. Therefore we make the assumption of a diagonal covariance matrix
while sacrificing the benefit of modeling the underlying correlations among
the random variables.  Then, we use the quantized vector
$\mathbf{U}^{\top}\bm{\xi}_{t}$ to update the state i.e. to estimate
$\tilde{\bm{x}}_{t}$ and $\bm{i}^{\ell}_{t}$. Algorithms 1 and 2 describe the
encoding and decoding algorithms respectively.

Algorithm 2 Decoding algorithm for reconstruction

1:$\tilde{\bm{x}}_{0},\tilde{\bm{x}}_{1},\bm{i}_{0}^{\ell},\bm{i}_{1}^{\ell},$
$\left\\{\hat{\tilde{\bm{\xi}}}_{t}\right\\}_{t=2}^{T}$

2:for $t=2:T$ do

3: States $\hat{\tilde{\bm{x}}}_{t},\hat{\bm{i}}_{t}^{\ell}$ from (60) and
(61) respectively.

4: Reconstruction of voltage from (58)
$\hat{\bm{v}}_{t}=\hat{\bm{v}}_{t}^{0}+\mathbf{U}\hat{\tilde{\bm{\xi}}}_{t}$

5: Update states $\hat{\tilde{\bm{x}}}_{t},\hat{\bm{i}}_{t}^{\ell}$ from (59)

6:end for

7:$\left\\{\hat{\bm{v}}_{t}\right\\}_{t=2}^{T}$

Note that the proposed scheme of compression is sequential unlike others in
literature. Several corrections can be made as data is collected in time such
as the update of parameters
$\tilde{\bm{a}}_{1},\tilde{\bm{a}}_{2},\bm{b}_{1},\bm{b}_{2}$.

Figure 2: Support of the GSO $\mathbf{S}$ of the network.

##  VI Numerical Results

The numerical results in this section are mostly obtained using data from the
synthetic ACTIVSg2000 case [78], a realistic model emulating the ERCOT system,
which includes $2,000$ buses-with $432$ generators and the rest non-generator
buses. The ACTIVSg2000 case data include a realistic PMU data time series, in
which $392$ generators are dispatched to meet variable load demand. The
sampling rate, as for real PMUs, is $30$ samples per second. As all the system
related parameters are known, it is easier to verify the proposed modeling
strategy through the ACTIVSg2000 PMU data set. Fig. 2 shows the support of the
graph Laplacian or the $\bm{Y}$ matrix when ordered into generator and non-
generator buses. The block-diagonal structure is notable, and is the result of
the population distribution in the state of Texas, which is concentrated in
$8$ metropolitan areas.

Grid-GSP model: In Fig. 3, magnitude of GFT of voltage graph signal
$\bm{v}_{t}$ and the input $\bm{x}_{t}$ are plotted for a single time instant
with respect to their corresponding normalized graph frequencies
$\lvert\lambda_{i}\rvert/\max_{i}\lvert\lambda_{i}\rvert$ and shown in log-
scale. From the linear decay, it is evident that the magnitude of GFT
coefficients $\lvert\tilde{\bm{v}}_{t}\rvert$ corresponding to lower
frequencies are more significant as compared to higher frequencies. Similarly,
the GFT of the exponent in the input,
$\tilde{\bm{x}}_{t}=\mathbf{U}_{\text{red}}^{\top}\bm{x}_{t}$ with the
generator GSO $\mathbf{S}_{\text{red}}$, is plotted with respect to the graph
frequencies in Fig. 3. The decay in GFT coefficients with respect to frequency
is less pronounced confirming that graph signal $\bm{x}_{t}$ is not
necessarily low-pass and in general depends on the topology of the generator
only network.

In Fig. 4, magnitude of GFT of the downsampled voltage graph signal,
$\tilde{\bm{v}}_{\mathcal{M}}=\mathbf{U}^{\top}_{\text{red},\mathcal{M}}{\bm{v}}_{\mathcal{M}}$
for $\lvert\mathcal{M}\rvert=867$ and $\lvert\mathcal{M}\rvert=261$ with two
different down-sampling strategies: with PMUs placed at buses in few
communities within the GSO $\mathbf{S}$ and the other being optimal placement
for graph signal reconstruction.

Figure 3: Magnitude of Graph Fourier Transform (GFT) for voltage graph signal,
$\lvert\tilde{\bm{v}_{t}}\rvert$ (left) and input to generator-only network
$\lvert\tilde{\bm{x}_{t}}\rvert$ (right) plotted with respect to normalized
graph frequency.  Figure 4: Magnitude of GFT for spatially downsampled voltage
graph signal, $\lvert\tilde{\bm{v}}_{\mathcal{M}}\rvert$ plotted with respect
to normalized graph frequency.

The placement strategy has an effect on the low-pass nature of the downsampled
signal. The steeper attenuation of GFT magnitude with placement strategy being
community-wise is a result of loss in spatial-resolution.

To highlight the temporal variation in the GFT domain of input exponent,
$\tilde{\bm{x}}_{t}$, a short time-series of the real and imaginary parts
along with the fit of the AR model are shown in Fig. 5. As expected, the AR
model fits well. Fig. 6 shows the similar AR-2 model fit to the load current
at a bus that had the highest absolute value of load.

Figure 5: AR model fit to the GFT of $\bm{x}_{t}$. Component corresponding to
smallest graph frequency, $\lambda_{\text{red},1}$,
$\left[\tilde{\bm{x}}_{t}\right]_{1}$ Figure 6: AR model for load bus,
$j=1312$, i.e. bus with highest absolute value of load.

To emphasize the temporal nature of the input, the $2$-dimensional frequency
response (in both graph and time domains) is plotted for the input
$\tilde{\bm{x}}_{t}$ in Fig. 7. The figure provides evidence of the coupling
between the graph frequencies and time series Fourier power spectrum, and the
variability of the temporal response depending on what GFT frequency mode is
excited with Fourier spectra that are more or less concentrated towards low
frequencies depending on the GFT mode. Hence, temporal dynamics can inform
about what is happening in space (i.e. the trends are coupled).

Figure 7: $2-$D frequency response for $\tilde{\bm{x}}_{t}$ Figure 8: Optimal
placement of $\lvert\mathcal{M}\rvert=100$ PMUs.
$\lvert\mathcal{K}\rvert=100$.  Figure 9: Reconstruction performance after
optimal placement [71] of $\lvert\mathcal{M}\rvert$ PMUs. Number of
frequencies used: $\lvert\mathcal{K}\rvert=\lvert\mathcal{M}\rvert$. For
random placement, $\lvert\mathcal{K}\rvert=100$ used.  Figure 10:  The map of
PMUs placed in ISO-NE test case $3$ [79] (left) and the support of estimated
GSO via (48) (right) shown. Note that the community structure corresponds to
groups of PMUs in the actual system as highlighted in the figure for a few
clusters.  Figure 11: Interpolation of missing measurements for an ISO-NE case
using GSO based regularization. Note the contiguous missing of samples and our
ability to interpolate. The relative noise level used is,
$(\lvert\mathcal{M}\rvert T)\sigma^{2}/\lVert\mathbf{V}\rVert_{F}^{2}=10^{-4}$
Normalized MSE for this run is $6.22\times 10^{-4}$.

Revising GSP tools: sampling and optimal placement Fig. 8 shows the placement
of $\lvert\mathcal{M}\rvert=100$ PMUs super-imposed on the support of the
ordered Grid-GSO, $\mathbf{S}$ when $\lvert\mathcal{K}\rvert=100$ graph
frequency components are considered. Note the distribution of PMUs to
different communities as well as on the generator buses as they belong to
different graph communities. Fig. 9 exhibits the performance of the GSP based
reconstruction method on optimally placed $\lvert\mathcal{M}\rvert$ PMUs that
provide down-sampled measurements, $\bm{v}_{\mathcal{M}}$. The number of graph
frequencies considered for reconstruction are
$\lvert\mathcal{K}\rvert=\lvert\mathcal{M}\rvert$. Even with just $5\%$ of
measurements ($100$ PMUs), the reconstruction error is extremely low. For
random placement, $\lvert\mathcal{M}\rvert$ PMUs are chosen at random and
$\lvert\mathcal{K}\rvert=100$ graph frequencies are used for reconstruction.
The trial of random placement is repeated $1,000$ times and the most
frequently occuring error (estimate of mode of the error distribution) is
plotted. As expected, the reconstruction error for random placement is orders
of magnitude higher than optimal placement.

To illustrate that the proposed modeling holds and algorithms work well also
for real PMU data, in the next numerical experiments we used a real-world
dataset of measurements from $35$ PMUs placed in ISO New-England grid (ISO-NE)
[79]. The data corresponds to a period of $180$ seconds when a large generator
near Ln:2 and Ln:4 introduces oscillations in the system. We decimated in time
the PMU signals down to sampling frequency $1$ sample/s.

Network inference: As the underlying GSO is unknown, it is estimated via (48)
with the goal of recovering the underlying reduced-GSO. Since admittance
values are not given, we only compare the support of the estimated GSO with
the community of PMUs in the network. Fig. 10 shows the support of the
estimated GSO and compares it with the map of PMUs highlighting a few clusters
of correspondence. From Fig. 10 we see that the block-diagonal nature of the
estimated GSO captures the community structure in the map.

Interpolation of missing measurements: Once the GSO is estimated, we consider
the interpolation problem in (47) for the same ISO-NE dataset. We delete data
at random and add noise. We solve the problem in (47) to recover missing
measurements. In Fig.11 we compare the original, corrupted and recovered
measurements. Corrupted measurements have missing samples not just at random
but also contiguous in time. The normalized MSE,
$\lVert\mathbf{V}-\hat{\mathbf{V}}\rVert_{F}^{2}/\lVert\mathbf{V}\rVert_{F}^{2}$
is the metric used to gauge the reconstruction performance. As a comparison,
we tested on the same data the AM-FIHT algorithm proposed in [80], which
regularizes the reconstruction task assuming that the Hankel matrix formed
with the columns of $\mathbf{V}$, i.e. ${\mathbb{H}(\mathbf{V})}$, has low
rank $r$,

$\displaystyle\min_{\mathbf{V}}\lVert\mathbb{P}_{\Omega}\left(\hat{\mathbf{V}}-\mathbf{V}\right)\rVert_{F}^{2}~{}~{}\text{subject
to}~{}~{}\text{rank}\left({\mathbb{H}(\mathbf{V})}\right)=r$ (63)

Figure 12: Comparison of AM-FIHT algorithm in [80] ($r=10,n1=3,\beta=0$) and
the proposed GSO based interpolation with $50\%$ of missing measurements in
the ISO-NE dataset.

The plot comparing the two methods is shown in Fig. 12. As seen, the GSO based
method outperforms the AM-FIHT for this dataset, indicating that the
regularization using the GSO is more effective at capturing the low-rank
nature of the data, compared to seeking an arbitrary low rank structure in the
the Hankel matrix of the data.

Figure 13: Components of projection of received measurement $\bm{z}$ onto the
orthogonal subspace,
$\bm{\Pi}_{\perp\bm{H}\mathcal{H}_{k}(\mathbf{S})}{\bm{z}}$,
$\lvert{\mathcal{A}}\rvert=500,\lvert\mathcal{K}\rvert=200,\lvert{\mathcal{C}}\rvert=250$.
Figure 14:  Empirical ROC curve for different $\lvert{\mathcal{C}}\rvert$ when
a percentage of them are malicious,
$\lvert{\mathcal{C}}\rvert/\lvert{\mathcal{A}}\rvert\times 100$ with
$\lvert{\mathcal{A}}\rvert=500$ (out of $2,000$) are available.  Figure 15:
Reconstruction of attack vector.  Figure 16: Empirical ROC curve for methods
proposed here and by Drayer & Routtenberg [25] when all voltage measurements
are available, $\lvert\mathcal{A}\rvert=2,000$. A percent of the measurements,
$\lvert{\mathcal{C}}\rvert/\lvert{\mathcal{A}}\rvert\\!\times\\!100$ are
malicious. The relative noise level is $10^{-2}$.

Detection of FDI attacks: Fig.13 shows the magnitude of the projection of the
received measurement $\bm{z}$ on the orthogonal subspace
$\bm{\Pi}_{\perp\bm{H}\mathcal{H}_{k}(\mathbf{S})}$. From Fig. 13 it is
evident that when there is no attack, the magnitude of the projected component
is orders of magnitude lower than when the measurements are under the FDI
attack. This validates the idea of using high GFT frequency activity as an
indicator of anomalies. Fig. 14 shows the empirical receiver operator
characteristics (ROC) curve highlighting the detection performance of the
proposed FDI attack detection scheme. The detection performance remains good,
even when very few buses are attacked. We compare the performance of the
proposed FDI attack detection with that of the method in [25] when the full
state i.e. when all voltage measurements are available,
$\lvert\mathcal{A}\rvert=2,000$. The underlying principle to detect the attack
in [25] is to look at the magnitude of graph frequency components at higher
graph frequencies which is similar in principle to the detection test we
undertake. They use the real and imaginary parts of the system admittance
matrix as $2$ GSOs, $\Re\\{\bm{Y}\\}$ and $\Im\\{\bm{Y}\\}$ respectively.
Their test statistic is comprised of four components that are the frequency
response of high-pass filtered real and imaginary voltage measurements (see
Algorithm.2 in [25]). Fig. 16 shows the empirical ROC curves that compare the
performance of the proposed method and the one in [25] when all voltage
measurements are available and are noisy. The relative noise level used is
$10^{-2}$. As evident from the curves, the proposed method performs better
than the method in [25]. This is because our test statistic is more robust to
noise and also significantly more sensitive in detecting the attack vectors,
even when only few buses are attacked.

Fig. 15 shows the reconstruction of magnitude of the attack vector
$\delta\bm{v}$ when $500$ measurements are available and number of attacked
buses $\lvert\mathcal{C}\rvert=50$.

Compression based results:

Figure 17: Empirical rate distortion (RD) curve for the proposed compression
method compared with singular value thresholding and quantization.

For voltage data compression, we compared with two schemes: scalar
quantization and singular value thresholding (SVT) from [49]. Fig. 17 plots
the empirical rate-distortion (RD) curve and shows the comparison between all
$3$ schemes. As expected, scalar quantization does poorly compared to the
other schemes. The SVT scheme simply uses few of the largest singular vectors
for data reconstruction. Considering that it is indicative of voltage graph
signal lying in a low-dimensional subspace, it is not surprising that the SVT
scheme does well. However, the SVT curve rate-distortion curve eventually
saturates. Note that the performance of the proposed method are comparable to
those of the SVT. However, the latter is a batch method, while the proposed
method is sequential, which has important implication for the online
communications of PMU data.

##  VII Conclusions

In this paper, we proposed the framework of Grid-GSP for the power grid that
highlights the inherent spatio-temporal structure in the voltage phasors by
employing concepts from GSP. Grid-GSP revisits the concepts of sampling and
reconstruction, interpolation, network inference and applications, to
detection of FDI attacks and a lossy sequential data compression, were
introduced using the lens of GSP. The resulting algorithms were tested on data
from both synthetic and real-world datasets. The paper opens the door to
leverage the GSP foundations for all types of grid data analytical tasks.

## References

  * [1] R. Ramakrishna and A. Scaglione, “On Modeling Voltage Phasor Measurements as Graph Signals,” in _IEEE Data Science Workshop (DSW) 2019_ , June 2019, pp. 275–279.
  * [2] R. Ramakrishna and A. Scaglione, “Detection of False Data Injection Attack using Graph Signal Processing for the Power Grid,” in _2019 IEEE Global Conference on Signal and Information Processing (GlobalSIP). IEEE_ , 2019\.
  * [3] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,” _IEEE Signal Processing Magazine_ , 2013.
  * [4] N. Tremblay, P. Gonçalves, and P. Borgnat, “Design of graph filters and filterbanks,” in _Cooperative and Graph Signal Processing_. Elsevier, 2018, pp. 299–324.
  * [5] A. Sandryhaila and J. M. F. Moura, “Discrete Signal Processing on Graphs: Frequency Analysis,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 12, June 2014.
  * [6] A. Phadke and J. Thorp, “History and Applications of Phasor Measurements,” in _2006 IEEE PES Power Systems Conference and Exposition_. IEEE, 2006, pp. 331–335.
  * [7] T. Ishizaki, A. Chakrabortty, and J.-I. Imura, “Graph-Theoretic Analysis of Power Systems,” _Proceedings of the IEEE_ , vol. 106, no. 5, pp. 931–952, May 2018.
  * [8] F. Dörfler and F. Bullo, “Kron Reduction of Graphs With Applications to Electrical Networks,” _IEEE Transactions on Circuits and Systems-I: Regular Papers_ , vol. 60, no. 1, 2013.
  * [9] F. Dorfler and F. Bullo, “Spectral Analysis of Synchronization in a Lossless Structure-Preserving Power Network Model,” in _2010 First IEEE International Conference on Smart Grid Communications_ , 2010.
  * [10] L. Xie, Y. Chen, and P. R. Kumar, “Dimensionality reduction of synchrophasor data for early event detection: Linearized analysis,” _IEEE Trans. on Power Systems_ , vol. 29, no. 6, pp. 2784–2794, Nov 2014.
  * [11] 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_ , vol. 31, no. 2, pp. 1006–1013, 2016.
  * [12] M. Wang, “Data quality management of synchrophasor data in power systems by exploiting low-dimensional models,” in _2017 51st Annual Conference on Information Sciences and Systems (CISS)_ , 2017.
  * [13] M. Jamei, A. Scaglione, C. Roberts, E. Stewart, S. Peisert, C. McParland, and A. McEachern, “Anomaly Detection Using Optimally-Placed PMU Sensors in Distribution Grids,” _IEEE Transactions on Power Systems_ , vol. 33, no. 4, pp. 3611–3622, July 2018.
  * [14] W. Li, M. Wang, and J. H. Chow, “Real-Time Event Identification Through Low-Dimensional Subspace Characterization of High-Dimensional Synchrophasor Data,” _IEEE Transactions on Power Systems_ , vol. 33, no. 5, Sept 2018\.
  * [15] 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_ , vol. 63, no. 5, pp. 1102–1114, March 2015.
  * [16] P. Du, Z. Huang, R. Diao, B. Lee, and K. K. Anderson, “PMU Placement for enhancing dynamic observability of a power grid,” in _2010 IEEE Conference on Innovative Technologies for an Efficient and Reliable Electricity Supply_ , Sept 2010, pp. 15–21.
  * [17] A. Pal, G. A. Sanchez-Ayala, V. A. Centeno, and J. S. Thorp, “A PMU Placement Scheme Ensuring Real-Time Monitoring of Critical Buses of the Network,” _IEEE Transactions on Power Delivery_ , vol. 29, no. 2, pp. 510–517, April 2014.
  * [18] Z. Wang, A. Scaglione, and R. Thomas, “Generating Statistically Correct Random Topologies for Testing Smart Grid Communication and Control Networks,” _IEEE Transactions on Power Systems_ , vol. 1, no. 1, pp. 28–39, 2010.
  * [19] Li, Xiao and Poor, H Vincent and Scaglione, Anna, “Blind topology identification for power systems,” in _2013 IEEE International Conference on Smart Grid Communications (SmartGridComm)_. IEEE, 2013, pp. 91–96.
  * [20] S. Grotas, Y. Yakoby, I. Gera, and T. Routtenberg, “Power Systems Topology and State Estimation by Graph Blind Source Separation,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 8, pp. 2036–2051, 2019.
  * [21] Z. Xiang, K. Huang, W. Deng, and C. Yang, “Blind Topology Identification for Smart Grid Based on Dictionary Learning,” in _2019 IEEE Symposium Series on Computational Intelligence (SSCI)_ , 2019, pp. 1319–1326.
  * [22] D. Deka, M. Chertkov, and S. Backhaus, “Joint Estimation of Topology and Injection Statistics in Distribution Grids With Missing Nodes,” _IEEE Transactions on Control of Network Systems_ , vol. 7, no. 3, pp. 1391–1403, 2020.
  * [23] S. Talukdar, D. Deka, H. Doddi, D. Materassi, M. Chertkov, and M. V. Salapaka, “Physics informed topology learning in networks of linear dynamical systems,” _Automatica_ , vol. 112, 2020.
  * [24] E. Drayer and T. Routtenberg, “Detection of False Data Injection Attacks in Power Systems with Graph Fourier Transform,” in _2018 IEEE Global Conference on Signal and Information Processing (GlobalSIP). IEEE, 2018_ , 2018, pp. 890–894.
  * [25] E. Drayer and T. Routtenberg, “Detection of False Data Injection Attacks in Smart Grids Based on Graph Signal Processing,” _IEEE Systems Journal_ , August 2019.
  * [26] M. Jamei, R. Ramakrishna, T. Tesfay, R. Gentz, C. Roberts, A. Scaglione, and S. Peisert, “Phasor Measurement Units Optimal Placement and Performance Limits for Fault Localization,” _IEEE Journal on Selected Areas in Communications_ , 2019.
  * [27] L. Fan, “Interarea Oscillations Revisited,” _IEEE Transactions on Power Systems_ , vol. 32, no. 2, pp. 1585–1586, 2017.
  * [28] T. Huang, N. M. Freris, P. R. Kumar, and L. Xie, “Localization of Forced Oscillations in the Power Grid Under Resonance Conditions,” in _2018 52nd Annual Conference on Information Sciences and Systems (CISS)_ , 2018.
  * [29] L. Guo, C. Zhao, and S. H. Low, “Graph Laplacian Spectrum and Primary Frequency Regulation,” in _2018 IEEE Conference on Decision and Control (CDC)_ , 2018, pp. 158–165.
  * [30] Y. Weng, R. Negi, and M. D. Ilić, “Graphical model for state estimation in electric power systems,” in _2013 IEEE International Conference on Smart Grid Communications (SmartGridComm)_. IEEE, 2013, pp. 103–108.
  * [31] D. Deka, M. Chertkov, and S. Backhaus, “Topology Estimation Using Graphical Models in Multi-Phase Power Distribution Grids,” _IEEE Transactions on Power Systems_ , vol. 35, no. 3, pp. 1663–1673, May 2020.
  * [32] K. Dvijotham, M. Chertkov, P. Van Hentenryck, M. Vuffray, and S. Misra, “Graphical models for optimal power flow,” _Constraints_ , vol. 22, no. 1, pp. 24–49, 2017.
  * [33] R. B. Bobba, K. M. Rogers, Q. Wang, H. Khurana, K. Nahrstedt, and T. J. Overbye, “Detecting false data injection attacks on DC state estimation,” in _Preprints of the First Workshop on Secure Control Systems, CPSWEEK_ , vol. 2010, 2010.
  * [34] G. Dan and H. Sandberg, “Stealth attacks and protection schemes for state estimators in power systems,” in _Smart Grid Communications (SmartGridComm), 2010 First IEEE International Conference on_ , 2010, pp. 214–219.
  * [35] O. Kosut, L. Jia, R. J. Thomas, and L. Tong, “Malicious data attacks on the smart grid,” _IEEE Transactions on Smart Grid_ , vol. 2, no. 4, pp. 645–658, 2011.
  * [36] G. Liang, J. Zhao, F. Luo, S. R. Weller, and Z. Y. Dong, “A review of false data injection attacks against modern power systems,” _IEEE Transactions on Smart Grid_ , vol. 8, no. 4, pp. 1630–1638, July 2017.
  * [37] I. Esnaola, S. M. Perlaza, H. V. Poor, and O. Kosut, “Maximum Distortion Attacks in Electricity Grids,” _IEEE Transactions on Smart Grid_ , vol. 7, no. 4, pp. 2007–2015, Jul. 2016.
  * [38] J. Zhang, Z. Chu, L. Sankar, and O. Kosut, “False data injection attacks on phasor measurements that bypass low-rank decomposition,” in _IEEE International Conference on Smart Grid Communications (SmartGridComm)_ , Oct. 2017\.
  * [39] Y. He, G. J. Mendis, and J. Wei, “Real-time detection of false data injection attacks in smart grid: A deep learning-based intelligent mechanism,” _IEEE Transactions on Smart Grid_ , vol. 8, no. 5, pp. 2505–2516, Sept 2017\.
  * [40] D. P. Shepard, T. E. Humphreys, and A. A. Fansler, “Evaluation of the vulnerability of phasor measurement units to gps spoofing attacks,” _International Journal of Critical Infrastructure Protection_ , vol. 5, no. 3-4, pp. 146–153, 2012.
  * [41] L. Heng, J. J. Makela, A. D. Dominguez-Garcia, R. B. Bobba, W. H. Sanders, and G. X. Gao, “Reliable GPS-based timing for power systems: A multi-layered multi-receiver architecture,” in _Power and Energy Conference at Illinois (PECI), 2014_ , 2014, pp. 1–7.
  * [42] 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_ , vol. 64, no. 21, pp. 5557–5570, Nov 2016\.
  * [43] P. Top and J. Breneman, “Compressing phasor measurement data,” in _2013 IEEE Power Energy Society General Meeting_ , July 2013, pp. 1–4.
  * [44] R. Klump, P. Agarwal, J. E. Tate, and H. Khurana, “Lossless compression of synchronized phasor measurements,” in _IEEE PES General Meeting_ , July 2010, pp. 1–7.
  * [45] J. E. Tate, “Preprocessing and Golomb-Rice encoding for lossless compression of phasor angle data,” _IEEE Trans. on Smart Grid_ , vol. 7, no. 2, pp. 718–729, March 2016.
  * [46] S. Kirti, Z. Wang, A. Scaglione, and R. Thomas, “On the Communication Architecture for Wide-Area Real-Time Monitoring in Power Networks,” in _2007 40th Annual Hawaii International Conference on System Sciences (HICSS’07)_ , 2007.
  * [47] P. H. Gadde, M. Biswal, S. Brahma, and H. Cao, “Efficient compression of PMU data in WAMS,” _IEEE Trans. on Smart Grid_ , vol. 7, no. 5, pp. 2406–2413, Sept 2016.
  * [48] Y. Ge, A. J. Flueck, D. K. Kim, J. B. Ahn, J. D. Lee, and D. Y. Kwon, “Power system real-time event detection and associated data archival reduction based on synchrophasors,” _IEEE Trans. on Smart Grid_ , vol. 6, no. 4, pp. 2088–2097, July 2015.
  * [49] J. C. S. de Souza, T. M. L. Assis, and B. C. Pal, “Data Compression in Smart Distribution Systems via Singular Value Decomposition,” _IEEE Transactions on Smart Grid_ , vol. 8, no. 1, January 2017.
  * [50] R. Mehra, V. Patel, F. Kazi, N. M. Singh, and S. R. Wagh, “Modes preserving wavelet based multi-scale PCA algorithm for compression of smart grid data,” in _2013 International Conference on Advances in Computing, Communications and Informatics (ICACCI)_ , Aug 2013, pp. 817–821.
  * [51] E. Isufi, A. Loukas, A. Simonetto, and G. Leus, “Autoregressive Moving Average Graph Filtering,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 2, pp. 274–288, 2017.
  * [52] R. Ramakrishna, H. T. Wai, and A. Scaglione, “A User Guide to Low-Pass Graph Signal Processing and Its Applications: Tools and Applications,” _IEEE Signal Processing Magazine_ , vol. 37, no. 6, pp. 74–85, November 2020.
  * [53] A. Sandryhaila and J. M. Moura, “Discrete Signal Processing on Graphs,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 7, pp. 1644–1656, 2013.
  * [54] R. A. Horn and C. R. Johnson, _Matrix Analysis_. Cambridge University Press, 1990.
  * [55] S. Mallat, _A Wavelet Tour of Signal Processing: The Sparse Way_. Academic press, 2009.
  * [56] R. Singh, A. Chakraborty, and B. Manoj, “Graph Fourier transform based on direcetd Laplacian,” in _2016 International Conference on Signal Processing and Communications (SPCOM)_ , 2016.
  * [57] E. Isufi, G. Leus, and P. Banelli, “2-Dimensional finite impulse response graph-temporal filters,” in _2016 IEEE Global Conference on Signal and Information Processing (GlobalSIP)_ , 2016.
  * [58] E. Isufi, A. Loukas, A. Simonetto, and G. Leus, “Filtering Random Graph Processes Over Random Time-Varying Graphs,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 16, August 2017.
  * [59] F. Grassi, A. Loukas, N. Perraudin, and B. Ricaud, “A Time-Vertex Signal Processing Framework: Scalable Processing and Meaningful Representations for Time-Series on Graphs,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 3, pp. 817–829, February 2018.
  * [60] E. Isufi, A. Loukas, A. Simonetto, and G. Leus, “Separable autoregressive moving average graph-temporal filters,” in _2016 24th European Signal Processing Conference (EUSIPCO)_. IEEE, 2016, pp. 200–204.
  * [61] J. D. Glover, M. S. Sarma, and T. J. Overbye, _Power system analysis and design_. Cengage Learning, 2008.
  * [62] H. P. Decell, Jr, “An application of the cayley-hamilton theorem to generalized matrix inversion,” _SIAM Review_ , vol. 7, no. 4, pp. 526–528, 1965.
  * [63] N. Sato and W. F. Tinney, “Techniques for Exploiting the Sparsity or the Network Admittance Matrix,” _IEEE Transactions on Power Apparatus and Systems_ , vol. 82, no. 69, pp. 944–950, 1963.
  * [64] P. Huynh, H. Zhu, Q. Chen, and A. E.Elbanna, “Data-Driven Estimation of Frequency Response From Ambient Synchrophasor Measurements,” _IEEE Transactions on Power Systems_ , vol. 33, no. 6, 2018.
  * [65] F. Paganini and E. Mallada, “Global Analysis of Synchronization Performance for Power Systems: Bridging the Theory-Practice Gap,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 7, pp. 3007–3022, 2020.
  * [66] F. Gao, J. S. Thorp, A. Pal, and S. Gao, “Dynamic state prediction based on Auto-Regressive (AR) Model using PMU data,” in _2012 IEEE Power and Energy Conference at Illinois_. IEEE, 2012, pp. 1–5.
  * [67] G. Golub and C. V. Loan, _Matrix Computations_. The Johns Hopkins University Press, 1990.
  * [68] M. Wang, J. H. Chow, P. Gao, X. T. Jiang, Y. Xia, S. G. Ghiocel, B. Fardanesh, G. Stefopolous, Y. Kokai, N. Saito, and M. Razanousky, “A low-rank matrix approach for the analysis of large amounts of power system synchrophasor data,” in _2015 48th Hawaii International Conference on System Sciences_ , Jan 2015, pp. 2637–2644.
  * [69] N. Dahal, R. L. King, and V. Madani, “Online dimension reduction of synchrophasor data,” in _IEEE PES Transmission and Distribution Conference and Exposition (T & D)_, 2012.
  * [70] M. Liao, D. Shi, Z. Yu, Z. Yi, Z. Wang, and Y. Xiang, “An Alternating Direction Method of Multipliers Based Approach for PMU Data Recovery,” _IEEE Transactions on Smart Grid_ , 2018.
  * [71] M. Tsitsvero, S. Barbarossa, and P. Di Lorenzo, “Signals on Graphs: Uncertainty Principle and Sampling,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 18, pp. 4845–4860, 2016.
  * [72] A. Anis, A. Gadde, and A. Ortega, “Efficient Sampling Set Selection for Bandlimited Graph Signals Using Graph Spectral Proxies,” _IEEE Transactions on Signal Processing_.
  * [73] L. Hagen and A. B. Kahng, “New Spectral Methods for Ratio Cut Partitioning and Clustering,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , vol. 11, no. 9, pp. 1074–1085, 1992.
  * [74] U. Von Luxburg, “A tutorial on spectral clustering,” _Statistics and computing, Springer_ , vol. 17, no. 4, pp. 395–416, 2007.
  * [75] X. Dong, D. Thanou, M. Rabbat, and P. Frossard, “Learning graphs from data: A signal representation perspective,” _IEEE Signal Processing Magazine_ , vol. 36, no. 3, pp. 44–63, 2019.
  * [76] W. Weber, “Differential Encoding for Multiple Amplitude and Phase Shift Keying Systems,” _IEEE Transactions on Communication_ , vol. 26, no. 3, 1978.
  * [77] T. M. Cover and J. Thomas, _Elements of Information Theory_. John Wiley, 1991.
  * [78] A. B. Birchfield, T. Xu, K. M. Gegner, K. S. Shetye, and T. J. Overbye, “Grid Structural Characteristics as Validation Criteria for Synthetic Networks,” _IEEE Transactions on Power Systems_ , vol. 32, no. 4, pp. 3258–3265, July 2017.
  * [79] S. Maslennikov, B. Wang, Q. Zhang, F. Ma, X. Luo, K. Sun, and E. Litvinov, “A Test Cases Library for Methods Locating the Sources of Sustained Oscillations,” in _IEEE PES General Meeting, Boston, MA_ , 2016.
  * [80] S. Zhang, Y. Hao, M. Wang, and J. H. Chow, “Multi-Channel Hankel Matrix Completion through Nonconvex Optimization,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 12, no. 4, pp. 617–632, 2018.