# A Model-Agnostic Method for PMU Data Recovery Using Optimal Singular Value
Thresholding

Shuchismita Biswas, and Virgilo A. Centeno

###### Abstract

This paper presents a fast model-agnostic method for recovering noisy Phasor
Measurement Unit (PMU) data streams with missing entries. The measurements are
first transformed into a Page matrix, and the original signals are
reconstructed using low-rank matrix estimation based on optimal singular value
thresholding. Two variations of the recovery algorithm are shown- a) an
offline block-processing method for imputing past measurements, and b) an
online method for predicting future measurements. Information within a PMU
channel (temporal correlation) as well as from different PMU channels in a
network (spatial correlation) are utilized to recover degraded data. The
proposed method is fast, and needs no explicit knowledge of the underlying
system model or measurement noise distribution. Performance of the recovery
algorithms is illustrated using simulated measurements from the IEEE 39-bus
test system as well as real measurements from an anonymized U.S. electric
utility. Extensive numeric tests show that the original signals can be
accurately recovered in the presence of additive noise, consecutive data drop
as well as simultaneous data erasures across multiple PMU channels.

###### Index Terms:

Phasor measurement unit (PMU), synchrophasor data, missing data recovery,
matrix estimation

##  I Introduction

Phasor Measurement Units (PMUs) allow high-resolution insight into power
systems dynamics through precise time-synchronized measurements [1]. In recent
years, electric utilities have made great strides towards deploying PMUs in
their networks and utilizing the reported measurements for wide-area
situational awareness [2]. Phasor measurements are used in both online (state
estimation, remedial action schemes) and offline (model validation,
contingency analysis, post-event diagnosis) applications [3],[4], [5]. As
shown in fig. 1, to reach the point-of-use from the point-of-measurement, PMU
data flows through various communication channels and intermediate data
concentrators. Hence, the data is susceptible to channel congestion or
component malfunction issues which lead to degraded information quality [6].
Reliable measurements are critical to developing robust control and monitoring
algorithms for the power grid, making fast and accurate data recovery critical
as well.

Traditional model-based recovery methods are sensitive to the underlying model
accuracy, and may need the knowledge of real-time network topology.
Oftentimes, such methods assume the knowledge of transmission line parameters,
and are adversely affected by inaccurate information [7]. Another class of
model-based methods use linear state estimation (LSE) for PMU data
conditioning [4, 8]. LSE-based data recovery is contingent on breaker status
information and any error or latency in reporting may deteriorate performance
[9]. On the other hand, measurements from PMUs within a network exhibit
spatial and temporal correlation that can be utilized to estimate missing
samples without explicit knowledge of the power system model. Several recent
papers have addressed model-agnostic PMU data recovery [10, 11, 12, 13, 14,
15, 16]. Some of the proposed methods target imputation of data blocks [13,
10, 11], while others are aimed at step-ahead prediction [12, 14, 15, 16].
Further, these recovery approaches broadly employ strategies based on - a)
filtering [12], b) low rank matrix completion [15, 14, 13, 16], or c) low rank
tensor completion [10, 11].

Figure 1: PMU data flow from point-of-measurement to point-of-use

In [12], the authors propose a kalman filter-based missing data estimation
algorithm that predicts the value of an incoming sample using the last three
measurements. The accuracy of this method deteriorates when more than three
consecutive entries are missing. Moreover, as the algorithm processes one PMU
channel at a time, information from other channels or nearby PMUs cannot be
leveraged to reconstruct segments of consecutive missing data. Matrix (tensor)
estimation methods propose to stack correlated PMU measurement channels
together to construct matrices (tensors), whose low-rank property can be
exploited to recover corrupt data. Tensor estimation methods are more
computationally expensive, limiting their potential for real-time use. Methods
to speed up tensor estimation have been explored in [10, 11].

In the matrix estimation (ME) area, various strategies have been put forth for
stacking measurement channels. In [13], the channels are stacked rowwise, and
an iterative alternating direction method of multipliers is used to fill
missing measurements. In [14], the channels are stacked columnwise, and the
spatial correlation among PMUs is used to recover missing entries using
singular value thresholding. Another variation that uses the temporal
correlation among channels is shown in [15]. Since these methods process
multiple PMU channels together, they are able to accurately recover missing
data sequences on one channel using information from other devices. On the
other hand, data prediction accuracy is severely affected by high noise
content in any one channel. An online recovery method exploiting the low-rank
property of the Hankel matrix constructed by overlapping segments of PMU data
has been proposed in [16]. As columns in a Hankel matrix overlap, this
approach faces the following drawbacks- a) repeated entries greatly increase
the size of the Hankel matrix, thereby increasing computation burden, b) noise
in the matrix elements are highly correlated, c) noisy entries are repeated
multiple times affecting recovery accuracy.

To overcome the limitations discussed above, this paper proposes a novel
technique for PMU data recovery. A sequence of PMU measurements is first
transformed into a Page matrix [17], and then recovered using a variation of
the truncated singular value thresholding algorithm [18]. As the matrix
columns consist of non-overlapping data segments, they are smaller in size
than Hankel matrices and the problem of highly correlated noise in matrix
entries is avoided. An optimal hard threshold is used for singular values;
thus the matrix rank does not need to be explicitly estimated at every step,
resulting in significant computational savings for online prediction. Two
variations of the algorithm are proposed, a) an offline imputation method for
archived PMU data, and b) an online one-step ahead prediction method.
Performance of the Page matrix method is compared to the Hankel matrix method
using the same estimation technique. Extensive numerical tests show that both
methods have similar recovery accuracy, and the Page matrix method is
computationally much faster.

The proposed algorithms can be applied to both univariate and multivariate
time-series (both single-channel and multi-channel cases). Although using
measurements from multiple PMUs translates to higher accuracy in data
recovery, some use cases might warrant processing single channels. For
instance, researchers may have access to limited PMU data due to their
sensitive nature. Moreover, computation may be sped up by parallelly
processing single measurement channels. This may be useful in cases where the
streaming synchrophasor data is of superior quality, and significant data
drops are not expected.

Contributions of this work may be summarized as follows. First, we propose a
model-agnostic method for recovering PMU measurements from noisy signals with
data drops. The proposed methodology is fast, scalable, easy to implement, and
poses minimal memory requirements, making it well-suited for real-time use.
Second, the methodology is extended into two algorithms - a) an offline method
intended for recovering archived data, and b) an online method for predicting
the next measurement, aimed at real-time applications. Third, through
extensive numerical experiments on simulated and real data, effectiveness of
the algorithms is verified. It is shown that the original measurement signals
can be reconstructed with high accuracy even in the presence of additive noise
and simultaneous data erasures across multiple channels.

The remaining paper is organized as follows. Section II describes the missing
data recovery problem. Section III succinctly presents the online and offline
recovery algorithms, and section IV illustrates their performance. Section V
concludes the paper and outlines future research directions.

##  II Problem Set-up

In this section, we establish the mathematical set-up for the PMU data
recovery problem and show how it relates to low-rank ME. Standard mathematical
notations are followed. Calligraphic symbols represent sets, lower case bold
letters represent column vectors, and upper case bold letters denote matrices.
All zero and all one vectors and matrices of appropriate size are denoted by
$\mathbf{0}$ and $\mathbf{1}$ respectively.

###  II-A PMU Data Recovery Using Matrix Estimation

Simply stated, PMUs are sensors deployed at different points of a power
network to measure electrical quantities like voltage and current magnitudes,
angles, frequency and rate of change of frequency. Measurements are time-
synchronized and typically reported at 30 or 60 frames per second (fps). Due
to the physical laws that govern power flow, measurements recorded by a PMU
and its neighbors are correlated. Moreover, data within a channel is
correlated in time. These spatial and temporal correlations may hence be
utilized to recover missing and corrupt measurements.

Formally, the data recovery problem can be posed as follows. Consider a
discrete-time setting with time instants indexed by $t\in\mathds{Z^{+}}$. Let
us say that for each $t\in[1,2,\dots T]$, PMU $i$ records a measurement vector
$\boldsymbol{x}_{i}(t)$ of length $c$, where $c$ is the number of measurement
channels. Measurements may contain observation noise, and it is assumed that
$\mathds{E}[\boldsymbol{x}_{i}(t)]=\boldsymbol{f}_{i}(t)$, where
$\boldsymbol{f}_{i}(t)$ are the true values of system states. Although the
underlying mean signal $\boldsymbol{f}_{i}(t)$ is strongly correlated in time,
it is assumed that the per-step noise are independent mean-zero random
variables with time-varying variance. Given some $\boldsymbol{x}_{i}(t)$, data
recovery algorithms may be designed to address two goals: a) imputation
(estimate $\boldsymbol{f}_{i}(t)$ for $t\in[1,2,\dots T]$), and b) prediction
(estimate $\boldsymbol{f}_{i}(T+1)$).

Time-series data recovery is a well-studied problem that appears in different
domains like econometrics, geosciences and healthcare. Classical methods for
time-series imputation and prediction have employed approaches such as hidden
Markov and state-space models [19]. Different deep neural networks (NN) have
also been used [20, 21]. Recent work has shown that low-rank matrix estimation
methods can provide simple, effective and computationally efficient means for
time-series recovery [22, 23]. This class of methods eliminates the training
data requirement of NN models, and hence provide a generalized framework
suitable for quick deployment.

The objective of ME is to recover a parameter matrix $\mathbf{M}$ from a
partially observed signal matrix $\mathbf{X}$ with corrupt entries, where
$\mathds{E}(\mathbf{X})=\mathbf{M}$. A detailed picture of the state-of-the-
art is available from [24, 25] and references therein. A key observation from
ME literature is that matrix $\mathbf{M}$ can be reconstructed from partial
and noisy observations by considering a low-rank approximation of the observed
matrix. ME algorithms are fairly model-agnostic in terms of the structure of
$\mathbf{M}$ and the distribution of $\mathbf{X}$ given $\mathbf{M}$.
Therefore, PMU measurements can be transformed into matrices and recovered by
applying ME methods to the transformed matrix. Truncated singular value
decomposition (SVD) based matrix estimation methods are popularly used [18].

###  II-B Matrix Transformation

Several methods have been proposed to transform time-series signals into
matrices. A naive method involves simply stacking signals together [14, 13,
15]. Although empirically this approach has shown reasonable effectiveness, it
cannot be used if very few measurement channels are available. Of course, this
is not a pressing concern for a transmission network with many PMUs. But in
some cases, it may be necessary to work with a limited number of measurement
signals. For instance, researchers outside electric utilities may only have
access to limited PMU data. Moreover, it may be desirable to process data from
blocks of few PMUs in a parallel manner to speed up computation during real-
time application.

Figure 2: Page and Hankel matrix construction from PMU measurements

An alternative approach described in [16] uses the Hankel matrix
transformation. Overlapping segments of PMU data are placed side by side to
form a Hankel matrix. That such matrices are approximately low-rank has been
empirically verified in [16]. As the Hankel matrices contain repeated entries,
they are large in size and noisy elements appear multiple times. Repetition of
noisy entries may reduce data recovery accuracy and the large matrix size
increases computation burden.

The limitations of Hankel matrices may be overcome by using Page matrices [22,
23]. A Page matrix is constructed from observation vector $\boldsymbol{x}(t)$
by placing contiguous segments of size $L>1$ side by side as non-overlapping
columns of the resultant matrix [17]. A schematic description of the matrix
transformation process is shown in fig. 2. It can be seen that when a time-
series of the same length is transformed, the Hankel matrix is much larger
than the corresponding Page matrix. The low-rank property of Page matrices has
also been examined in [22]. It is established that, in expectation, for a
large class of processes, the generated Page matrix is either exactly or
approximately low-rank. These proceeses include linear recurrent functions
(LRF) described by $f(t)=\sum_{g=1}^{G}\alpha_{g}f(t-g)$, for some $G\geq 1$.
That power systems quantities such as bus voltages, line currents and
frequency follow LRF has been previously concluded in literature. For example,
it is posited that steady state power systems measurements at any instant are
a linear combination of the last three measurements [26]. This result has also
been used to deisgn the kalman filter-based data prediction and smoothing
algorithm in [12].

The LRF nature of PMU measurements also helps in formulating a prediction
algorithm. Once the low-rank approximation of the transformed Page matrix is
obtained via some ME technique, the last row of this matrix can be expressed
as a linear combination of the other rows. Therefore, future values can be
forecast by applying linear regression to the approximate Page matrix.

###  II-C Optimal Singular Value Thresholding (OSVT)

Various techniques exist for low-rank ME [25, 24]. In the proposed model-
agnostic recovery framework, no information about the rank of the Page
matrices is available beforehand. Hence, a variation of the truncated SVD
method that does not need matrix rank information has been used in this paper
[18]. The optimal singular value threshold is chosen based on findings
reported in [27]. The main steps in the estimation algorithm are detailed next
in algorithm 1. It is assumed that missing data points in the observation
matrix are preliminarily filled in by the last available observation.

Algorithm 1 Optimal Singular Value Thresholding (OSVT)

1:Scaling:  Entries of the observation matrix $\mathbf{X}$ are scaled to lie
in the interval -1 to 1. Let the scaled observation matrix be called
${\mathbf{Y}}$ with individual entries $y_{i,j}$. Mathematically,
${y}_{i,j}=({x}_{i,j}-0.5(a+b))/0.5(b-a)$, where $a$ and $b$ are the minimum
and maximum entries of $\mathbf{X}$ respectively.

2:Singular value decomposition:  Let
${\mathbf{Y}}=\sum_{i=1}^{m}\sigma_{i}\boldsymbol{u}_{i}\boldsymbol{v}_{i}^{T}$
be the singular value decomposition of ${\mathbf{Y}}$. The singular values are
given by $\sigma_{i}$; and $\boldsymbol{u}_{i}$ and $\boldsymbol{v}_{i}$ are
the corresponding left and right singular vectors respectively.

3:Singular value thresholding:  Choose a set $S$ of thresholded singular
values such that : $S:=\\{\sigma_{i}>\sigma_{th}\\}$, where the optimality
threshold $\sigma_{th}$ is given by:
$\sigma_{th}=\sqrt{2(\zeta+1)+\frac{8\zeta}{(\zeta+1)+\sqrt{\zeta^{2}+14\zeta+1}}}$
Here, $\zeta=m/n$, where ${\mathbf{Y}}$ is a $m\times n$ matrix. Moreover,
$m\leq n$. In the case that $m>n$, the estimation algorithm must be applied to
$\mathbf{X}^{T}$ to obtain an estimate of $\mathbf{M}^{T}$.

4:Low rank approximation:  The low rank approximation of matrix $\mathbf{Y}$
is given by $\mathbf{\hat{Y}}=\sum_{\sigma_{i}\in
S}\sigma_{i}\boldsymbol{u}_{i}\boldsymbol{v}_{i}^{T}$. The final estimation
$\hat{\mathbf{M}}$ of the parameter matrix is obtained by scaling back the
values of $\hat{\mathbf{Y}}$ to the interval $[a,b]$.

The hard singular value threshold proposed in [27] is optimal in an asymptotic
sense. It is postulated that for large low-rank matrices, when a data singular
value $\sigma_{i}$ is too small, the corresponding singular vectors
$\boldsymbol{u}_{i}$ and $\boldsymbol{v}_{i}$ are very noisy and the component
$\sigma_{i}\boldsymbol{u}_{i}\boldsymbol{v}_{i}^{T}$ should not be used in
approximating matrix $\hat{\mathbf{Y}}$ from $\mathbf{Y}$. The cutoff for
singular values is determined to be $\sigma_{th}$ as described in algorithm 1.
The alternative method used in literature for estimating matrix rank involves-
a) selecting a threshold for rank approximation error, and b) choosing the
lowest rank for which the approximate matrix does not violate the
predetermined error threshold. The choice of the approximation error threshold
is somewhat arbitrary and has been empirically decided in works like [16].
Using a hard threshold eliminates the need for repeated calculations of rank
approximation error at every step, thereby significantly improving computation
speed. It is further shown in [27] that the optimally tuned thresholding
method outperforms (in terms of mean squared error) classical truncated SVD
when signal noise content is low to moderate; and the methods perform roughly
similarly when noise content is high.

###  II-D Multivariate Time-Series Recovery

Figure 3: Stacked Page matrix for multivariate data recovery

As mentioned before, PMU measurements in a network are correlated, and hence
information from different PMUs can be utilized to recover degraded data.
Readings from multiple PMU channels can be transformed into a ‘stacked’ Page
matrix by concatenating individual Page matrices columnwise, as shown in fig.
3. The low-rank property of stacked Page matrices for a large class of
processes including LRF has been verified in [23]. Let $N$ be the total number
of PMU measurement channels available. The data recovery process can now be
summarized with the following steps: a) transform observation vectors
$\boldsymbol{x}_{i}(t),i\in[1,2,\dots,N],t\in[1,2,\dots,T]$ to Page matrices
$\mathbf{P}_{i},i\in[1,2,\dots,N]$, each with dimension $L\times T/L$, by
placing non-overlapping contiguous data segments of length $L>1$ (an
algorithmic hyperparameter) side-by-side as columns; b) concatenate the
individual Page matrices columnwise to form a stacked Page matrix $\mathbf{X}$
of dimension $L\times(NT/L)$; c) perform OSVT (algorithm 1) on matrix
$\mathbf{X}$ to obtain a denoised estimate $\hat{\mathbf{X}}$; d) learn a
forecasting model that expresses the last row of $\hat{\mathbf{X}}$ as a
linear combination of its remaining rows. The in-sample imputation and
prediction error for this method scales as $1/\sqrt{NT}$ [23].

##  III Data Recovery Algorithms

This section shows how the individual pieces described in the previous section
are put together for recovering degraded PMU data. The algorithms can be
applied to both single-channel and multi-channel data using suitable Page
matrices.

###  III-A Offline Data Imputation

Let us consider the problem of denoising already recorded PMU measurements and
imputing missing readings. Assume $n$ number of PMUs, each with $c$
measurement channels. Then the total number of data channels available is
$N=n\times c$. The observations are partitioned into say $k$ windows of length
$T$ each. Note here that for the imputation task, hours of data can be
processed at once. Thus, $T$ can be quite large. It is assumed here that $T$
is perfectly divisible by the value of $L$ chosen. The sequential steps to be
performed for offline data imputation are listed in algorithm 2.

Algorithm 2 Offline PMU data imputation

1:Initialization: Set $T\leftarrow$ window length, $N\leftarrow$ number of PMU
channels, $k\leftarrow$ number of data windows.

2:For $j=1:k$, do

3: Matrix transformation: Construct a stacked Page matrix using data from $N$
PMU channels in the $j^{th}$ measurement window with $T$ observations, as
described in section II-D.

4: Low-rank matrix estimation: Obtain a low-rank approximation of the stacked
Page matrix $\hat{\mathbf{X}}$ using the OSVT method outlined in algorithm 1.

5: Recover estimated measurements: Reshape matrix $\hat{\mathbf{X}}$ to
recover the estimated measurements.

6:End for

###  III-B Online Data Prediction

The online forecast problem pertains to predicting the signal value
$\boldsymbol{f}(T+1)$ given past observations
$\boldsymbol{x}(t),t\in[1,2,\dots,T]$. This is analaogous to performing
regression with noisy data. The online forecasting algorithm proposed in this
paper- first, denoises and imputes past observations (algorithm 2), and
second, uses linear regression to learn the relationship between the last row
and remaining rows of the imputed observation matrix. Next, the learned
regression parameters are used to predict the next sample from a Page matrix
shifted by one sample. Of course, the first data window of length $T$ needs to
be filled before next-step prediction can proceed. The forecast procedure is
described with better clarity in algorithm 3.

Figure 4: Visual description of matrices
$\mathbf{G},\mathbf{H},\mathbf{G^{\prime}}$ and $\mathbf{H^{\prime}}$ for a
Page matrix using data from a single channel Algorithm 3 Online PMU data
prediction

1:Initialization: Set $T\leftarrow$ window length, $N\leftarrow$ total number
of PMU channels, $j\leftarrow 0$.

2:While PMU data streams are available, do:

3: Matrix formation: Convert measurement vectors $\boldsymbol{x_{i}}(t+j)$,
$i\in[1,2,\dots N]$, $t\in[1,2,\dots T]$ to a stacked Page matrix
$\mathbf{X_{T+j}}$, say.

4: Matrix imputation: Using algorithm 2, denoise and impute $\mathbf{X_{T+j}}$
to obtain ${\mathbf{\hat{X}_{T+j}}}$.

5: Learn the linear forecast model: Partition matrix $\mathbf{\hat{X}_{T+j}}$
into two parts $\mathbf{G}$ and $\mathbf{H}$ such that $\mathbf{G}$ comprises
of the first $L-1$ rows in $\mathbf{\hat{X}_{T+j}}$ and $\mathbf{H}$ contains
the last. Linear regression here pertains to estimating the parameter vector
$\boldsymbol{\beta_{j}}$ for which
$\mathbf{H}=\mathbf{G}{\boldsymbol{\beta_{j}}}+\boldsymbol{\epsilon}$ in the
least squares sense.

6: Forecast the one-step-ahead data: Construct matrix $\mathbf{G^{\prime}}$
with the last $L-1$ rows of $\mathbf{\hat{X}_{T+j}}$. Estimate
$\mathbf{H^{\prime}}$ from $\mathbf{G^{\prime}}$ as
$\mathbf{H^{\prime}}=\mathbf{G^{\prime}}\boldsymbol{\beta_{j}}$ using the
value of $\boldsymbol{\beta_{j}}$ learnt in step 5. For the univariate
prediction case, the last entry of $\mathbf{H^{\prime}}$ is the prediction of
the next measurement $x(T+j+1)$. For the multivariate case, pertinent entries
from $\mathbf{H^{\prime}}$ need to be extracted. For clarity, a visual
description of the matrices $\mathbf{G},\mathbf{H},\mathbf{G^{\prime}}$ and
$\mathbf{H^{\prime}}$ is provided in fig. 4.

7: Update: Set $j\leftarrow j+1$

8:End while

Figure 5: PMU measurements: Positive sequence voltage magnitude

Figure 6: PMU measurements: Positive sequence voltage angle (unwrapped)

Figure 7: PMU measurements: Frequency

###  III-C Choice of Hyperparameters

Selecting good hyperparameters is essential for achieving high-accuracy data
recovery while minimizing computation time. For the data recovery algorithms
proposed in this work, two main hyperparameters need to be chosen- a) $L$ or
the number of rows in the Page matrix, and b) $T$ or length of data window.

$\bullet$ Choice of $L$: Empirically, it was observed that choosing a $L$
value in the 5-10 range yielded good results for both the online prediction
and offline imputation tasks. Keeping the parameter $L$ small enables
capturing the short-term temporal patterns in PMU data. However, for very
noisy data, increasing $L$ can help in obtaining smoother estimates.

$\bullet$ Choice of $T$: The in-sample imputation and prediction error for the
multivariate data recovery method scales as $1/\sqrt{NT}$. Therefore, the
longer the data window, the better the prediction accuracy. On the other hand,
choosing a long data window increases computation time. As computation time is
not a prime concern for offline data imputation, a fairly long window can be
chosen. In the numerical results section of this paper, data imputation with
$T=54000$ (30 minutes of data) has been demonstrated .

Computation time is of greater concern when it comes to online prediction, as
every prediction step requires a matrix estimation and linear regression
operation. Hence, choosing a shorter time window is beneficial. However, the
window length $T$ should be selected carefully. If $T$ is too small,
prediction accuracy will be impacted by measurement noise. If $T$ is too
large, the data window may contain obsolete modes thereby degrading prediction
accuracy. It was empirically observed that a data window of about $\sim$ 30-45
samples yielded good performance without unduly increasing prediction time.

###  III-D Scalability

The size of the Page matrix will also be determined by the number of PMU
channels available. In a real system with hundreds of field PMUs, computation
may be sped up by dividing the PMUs into groups with similar modal signatures
and processing the groups in parallel. Grouping together PMUs geographically
close to each other will also enhance visibility into ‘local’ dynamics which
might have been obscured by aggregating signals over a wide-area grid. Some
strategies for grouping PMU signals to ensure low-rankness of the measurement
window have been presented in literature [28].

During a disturbance, the system deviates from its predicted behavior, and the
low-rank property of the stacked Page matrix may not hold true. Therefore, at
the onset of a disturbance, online predictions may vary greatly from observed
measurements. The difference between the actual observations and algorithm
predictions decreases gradually.

##  IV Numerical Results

This section describes the numerical tests performed to validate the
performance of the proposed imputation and prediction methodologies. The first
set of tests are performed on simulated data, artificially distorted by the
injection of additive noise and random data drop. The next set of tests
considers real noisy PMU measurements recorded by a U.S. utility. in this
case, the true values of the measurements are unknown, but visual inspection
reveals that reasonable values are predicted for swathes of missing data. All
computations are performed on a PC with 16 GB RAM and 2.6 GHz Intel core
i7-9750HF processor.

###  IV-A Simulated Measurements Dataset

Numerical tests were performed on 86.6 seconds of measurements from PMUs
installed at ten generator bus terminals of the IEEE 39-bus transmission
model[29]. Quasi-steady state operations and three-phase faults were recorded
using RTDS power systems simulator and GTNETx2 based PMUs at 60 fps reporting
rate. During the length of the simulation, quasi steady-state conditions were
simulated by modulating the mechanical torque of generator $G_{1}$ every 200
ms by a random perturbation within $\pm 1\%$ of the nominal value. The data
also shows three disturbances. At 18.33 seconds, a self-clearing three-phase
fault is followed by tripping of the faulted line, leading to a topology
change. At 55.67 seconds, the tripped line is reconnected, restoring the
initial network topology. At 78.13 seconds, another three-phase self-clearing
three-phase fault takes place.

In this study, three measurement channels from each PMU were used - positive
sequence voltage magnitude, positive sequence voltage angle and frequency.
Data from all ten PMUs are shown in figures 7, 7 and 7. PMUs are referred to
by the generator terminal they are installed at. For example, the PMU at
generator $G_{1}$ terminal is called PMU $G_{1}$. Voltage angle at PMU $G_{1}$
is considered as the reference angle.

Simulated measurements have been used for evaluating the data recovery
algorithms as the ‘ground truth’ data is available for comparison. On the
other hand, in real PMU measurements, some readings may already be missing or
corrupt, and there is no way of exactly knowing what those measurements should
have been. Further, as this simulation records network topology changes, we
can investigate if varying topologies affect recovery accuracy.

###  IV-B Data transformation

The PMU measurements are scaled before being transformed into the stacked Page
matrix described in section II-D. The scaling process used in this paper is
described below:
$\bullet$ Voltage magnitude: Measurements were transformed into the per unit
(p.u.) system.
$\bullet$ Voltage angle: The reference voltage angle is subtracted from
individual channel data. In the dataset used, angle readings were already
unwrapped. For unwrapping angles in real-time, the strategy outlined in [30]
may be followed.
$\bullet$ Frequency: Frequency measurements were scaled as follows:
$f_{scaled}=(f_{measured}-60)\times 10$.

The scaling method described above empirically showed good results; however
other approaches may also be used.

###  IV-C Error Metric

Mean absolute percentage error (MAPE) has been used as the error metric to
evaluate the accuracy of the proposed algorithms. Mathematically, MAPE for a
time-series of length $n$ maybe expressed as:

$\displaystyle
MAPE=\frac{1}{n}\sum_{t=1}^{n}\bigg{|}\frac{A_{t}-\hat{A_{t}}}{A_{t}}\bigg{|}$

Here, $A_{t}$ is the actual time-series value and $\hat{A_{t}}$ is the
corresponding prediction. As it intuitively conveys relative error, MAPE is
widely used in regression problems and model evaluation tasks.

Figure 8: Error in $G_{2}$ PMU channels with simultaneous data drop in all PMU
channels (median over 20 runs)

Figure 9: Error in $G_{2}$ PMU channels with zero-mean additive noise on all
$G_{2}$ PMU channels (median over 20 runs)

###  IV-D Offline Data Imputation

(a) Positive sequence voltage magnitude

(b) Positive Sequence voltage angle (referenced)

(c) Frequency

Figure 10: Imputed measurements with 50% simultaneous missing data and 2%
additive noise on all $G_{2}$ PMU channels

(a) Positive sequence voltage magnitude

(b) Positive Sequence voltage angle (referenced)

(c) Frequency

Figure 11: Predicted measurements with 50% simultaneous missing data and 2%
additive noise on all $G_{2}$ PMU channels

As communication channel congestion may impact all PMU channels, data erasures
may be correlated. To capture the efficacy of the proposed imputation method
under realistic data degradation conditions, the following modes were checked:

$\bullet$ Data drop: Simultaneous data drops on all PMU channels in the
network were considered. A fraction of timestamps (determined by the chosen
data drop rate) were randomly selected and corresponding measurements were
dropped from all PMU channels. As the data drop rate is increased, the chance
of missing consecutive data segments also increases. Error on all channels of
PMU $G_{2}$ as the data drop rate is varied is shown in fig. 9. These are
median values of observations over 20 runs. It is evident that the
measurements can be reconstructed with acceptable accuracy. Time taken for
imputation did not vary significantly with data drop rate and median time
taken over 100 runs was 0.0184 seconds. In the interest of brevity, only
results for PMU $G_{2}$ has been included in this paper, but similar results
were obtained for other PMU channels as well.

Results for the extreme scenario where readings from all PMUs are missing are
shown here. For less extreme cases, i.e. when data from only some of the PMU
channels are missing, higher accuracy in signal reconstruction may be
expected.

$\bullet$ Additive noise: Noise in PMU measurements may arise due to errors in
calibration, instrumentation and quantization. Existing studies suggest that
the zero-mean gaussian distribution is a suitable model to characterize this
noise. The signal-to-noise ratio (SNR) for real transmission-level PMUs is
estimated to be around $\sim$45 dB, while for distribution-level PMUs the SNR
is estimated to be lower [31]. Similar noise was injected into PMU $G_{2}$
channels for the numeric tests in this work.

Zero-mean gaussian noise was added to all PMU $G_{2}$ channels. The standard
deviation of the noise distribution on each channel was given by a percentage
of the median of true steady-state data (let us call this percentage the noise
rate). Data recovery error with varying noise rates is shown in fig. 9. Median
computation time over 100 runs was 0.0192 seconds.

For both the cases above, the number of rows in the stacked Page matrix was
10. Fig. 11 shows both the corrupt and imputed data when 2% noise was added to
the PMU $G_{2}$ channels and 50% of the readings were missing. It can be seen
that the PMU signals are reconstructed with reasonable accuracy even when
consecutive data segments are missing. Further, the data recovery accuracy is
robust to topology changes in the power network.

Figure 12: Prediction error in $G_{2}$ PMU channels with simultaneous missing
data on all $G_{2}$ PMU channels (median over 20 runs)

Figure 13: Prediction error in $G_{2}$ PMU channels with additive noise on all
$G_{2}$ PMU channels (median over 20 runs)

###  IV-E Online Data Prediction

Similar numerical tests were conducted with the same PMU dataset for assessing
online data prediction accuracy. Performance of the stacked Page-matrix based
prediction algorithm has been compared with the Hankel-matrix based prediction
method put forth in [16]. For better comparison, the same ME technique is used
for both algorithms.

$\bullet$ Prediction error: Figures 13 and 13 show the prediction error for
different data drop rates (simultaneous) and additive noise on all $G_{2}$ PMU
channels. The results are median observations over 20 runs. Data drops and
noise were introduced in the data in the same manner as discussed in section
IV-D. It is observed that the prediction error of the Page matrix method is
similar to/slightly better than the Hankel matrix based method. The corrupt
and reconstructed signals when 50% data drop and 2% noise is added to all
$G_{2}$ PMU channels is shown in fig. 11. In the experiments, number of rows
used was 5, and window length considered was 25.

$\bullet$ Prediction time: One-step ahead predictions have multiple
applications. Missing samples can be filled in with predicted values.
Similarly, irregularities in data may be detected looking at how far
measurements stray from their predictions. Now, for any real-time algorithm
implementation, computation time is of prime concern. U.S. electric utilities
typically use PMUs with reporting rates of 30 or 60 fps. For these PMUs, the
time intervals between two consecutive samples is 0.0333 or 0.0167 seconds.
Therefore, in order to predict a sample before it arrives at the control
center, the prediction time must be much lower.

Our numeric tests showed that the prediction time for the proposed Page matrix
algorithm was $\sim$0.001 seconds. In comparison, time taken by the Hankel
matrix based method was $\sim$0.003 seconds. All times are median values
recorded over 100 runs. The prediction times did not vary significantly with
varying rates of data drop or noise. Thus, we see that the method proposed in
this paper is much faster than the Hankel-matrix based method, and provides
similar/slightly better accuracy.

$\bullet$ Window length of Hankel matrix: The preceding set of experiments
showed that the proposed Page matrix-based data recovery method provided
accuracy levels similar to the Hankel matrix-based method, while speeding up
computations. The computational savings are largely due to the smaller size of
the Page matrix for the same measurement window. A natural question arises
here: how would using a shorter time window for the Hankel matrix method
affect data recovery accuracy? It is expected that predictions obtained using
a smaller number of observations will be noisy and less accurate. This notion
is experimentally verified in this work.

Figure 14: Prediction error in $G_{2}$ PMU channels with simultaneous missing
data on all $G_{2}$ PMU channels (median over 20 runs)

Figure 14 shows how the recovery accuracy varies when data is dropped
simultaneously from all channels of PMU $G_{2}$ using the same methodology as
section IV-D. Three matrix transformations were checked, a) Page matrix
(dimensions $5\times 6$), b) Hankel matrix with longer window (dimensions
$5\times 26$), and c) Hankel matrix with shorter window (dimensions $5\times
6$). For the first two cases, measurement window length used for predicting
the next sample is $T=30$, while for the third case it becomes $T=10$. As the
Page matrix and Hankel matrix with shorter window have the same dimensions,
time taken to predict the next sample using these matrices is almost the same.
However, as evident from figure 14, the Page matrix method yields higher
recovery accuracy. When using a Hankel matrix with longer window, the
prediction accuracy improves, but computation time increases as well.

Figure 15: Matrix rank for different data windows

$\bullet$ Verifying the low rank-property: The low rank properties of both the
stacked Page and Hankel matrices were checked for all the data windows, as
shown in fig. 15. It can be seen that both the Page and Hankel matrices are
generally low-rank, and the matrix ranks increase at the beginning of events.
During the first event, the Hankel matrix becomes full-rank, and the low-rank
property does not hold.

###  IV-F Real PMU measurements dataset

Next, we check how well the proposed algorithms perform with real PMU data
from an anonymized U.S. electric utility. The data corresponds to 30 minutes
of measurements (54000 samples) from four PMUs of the utility, each of which
reports three channels- positive sequence voltage magnitude, positive sequence
voltage angle, and frequency. PMU reporting rate is 30 fps. The voltage angle
of PMU4 is considered as reference, since it has the least amount of missing
entries. Voltage angles have been unwrapped using the algorithm from [30].

The recorded measurements were quite noisy, and had large patches of missing
data. The percentage of missing entries, and maximum length of consecutive
missing data segments in each channel is shown in fig. 16. It can be seen that
the PMU channels have data missing in the $\sim$20-40% range. Maximum length
of missing data segments is $\sim$100\. The noisy and intermittent PMU
measurements are shown in fig. 19, and the imputed measurements are shown in
fig. 19. Since there is no way to know what the actual measurements should
have been, objectively evaluating the imputation algorithm is not possible.
However, visual inspection shows that the proposed algorithm is able to impute
the measurements very well. Median computation time to impute 54000
measurement samples over 20 runs was 0.109 seconds.

Figure 19 shows the imputed measurements when using the Hankel matrix method
(in algorithm 2, measurements are transformed into a stacked Hankel matrix in
place of a Page matrix). Visually, it appears that the Hankel method does not
provide better estimates than the Page method. Moreover, time taken for
imputation using the Hankel method was $\sim 1.74$ seconds, much higher than
the Page method.

Further, additional $2\%$ noise was injected into the channels of PMU1 and
recovered using the proposed offline imputation method using Page matrix. The
noisy and imputed data is shown in figure 20. It can be inferred that the
proposed imputation method is able to recover data from real noisy PMU
archives.

Figure 16: Data quality in the real PMU dataset

Figure 17: Thirty minutes of measurements from four PMUs in an anonymized U.S.
electric utility. Data is reported at 30 fps. The vertical lines show data
drops.

Figure 18: Imputed PMU data using Page matrix

Figure 19: Imputed PMU data using Hankel matrix Figure 20: Imputed readings
for PMU1 after additional noise is injected into its channels

##  V Conclusion

In this work, a simple model-agnostic data recovery method based on low-rank
matrix approximation has been proposed for improving the quality of phasor
measurements with additive noise and data drop. The method is applicable for
both single and multiple measurement channels; and can deal with simultaneous
and consecutive data drop on all channels. Two variations of the recovery
algorithm are shown- a) an offline block-processing method for imputing past
measurements, and b) an online method for predicting future measurements. The
performance of the proposed algorithms have been illustrated through extensive
numeric experiments on simulated measurements on the IEEE 39-bus test system.
It is seen that the proposed methodology has high accuracy, has low memory
requirement and is computationally faster than other methods in literature.
The performance of the algorithm is independent of the underlying system
model, topology changes, and measurement noise distribution. Tests on real PMU
data validate the performance of the data recovery strategy put forth in this
work. The fast computation speed and ease of implementation make the
algorithms developed in this work suitable for quick deployment.

## References

  * [1] J. De La Ree, V. Centeno, J. S. Thorp, and A. G. Phadke, “Synchronized phasor measurement applications in power systems,” _IEEE Trans. on Smart Grid_ , vol. 1, no. 1, pp. 20–27, Jun. 2010.
  * [2] “Advancement of synchrophasor technology in projects funded by the American Recovery and Reinvestment Act of 2009,” Oak Ridge National Laboratory, Tech. Rep., Mar. 2016.
  * [3] S. Brahma, R. Kavasseri, H. Cao, N. R. Chaudhuri, T. Alexopoulos, and Y. Cui, “Real-time identification of dynamic events in power systems using PMU data, and potential applications—models, promises, and challenges,” _IEEE Trans. on Power Delivery_ , vol. 32, no. 1, pp. 294–301, 2017.
  * [4] K. D. Jones, “Synchrophasor-only dynamic state estimation & data conditioning,” Ph.D. dissertation, Virginia Tech, 2013.
  * [5] X. Fan, M. A. Elizondo, P. V. Etingov, M. R. Vallem, S. Biswas, S. Pal, C. Erroa, C. Munoz, D. Polanco, and V. Villeda, “Model validation study for Central American regional electrical interconnected system,” in _IEEE Power & Energy Soc. Innov. Smart Grid Technol. Conf. (ISGT)_, Washington D.C, USA, Feb. 2021, pp. 01–05.
  * [6] NASPI PMU Applications requirements Task Force, “A Framework for the Attributes of PMU Data Quality and a Methodology for Examining Data Quality Impacts to Synchrophasor Applications,” NASPI, Tech. Rep., 2017.
  * [7] D. Shi, D. J. Tylavsky, and N. Logic, “An adaptive method for detection and correction of errors in PMU measurements,” _IEEE Trans. on Smart Grid_ , vol. 3, no. 4, pp. 1575–1583, 2012.
  * [8] T. Faris and L. Zhang, “LSE for synchrophasor data quality – implementation and performance at BPA,” Mar. 2016. [Online]. Available: https://www.naspi.org/sites/default/files/2016-10/bpa_faris_lse_synchrophasor_data_quality_20160323.pdf
  * [9] J. Thorp, “Linear state estimator and synchrophasor data conditioning and validation,” in _WECC Joint Synchronized Information Subcommittee Meeting_ , Tempe, AZ, Jan. 2013. [Online]. Available: https://www.wecc.org/Reliability/WECC%20JSIS%20LSE%20and%20
Data%20ValidationVirginaTech.pdf

  * [10] D. Osipov and J. H. Chow, “PMU missing data recovery using tensor decomposition,” _IEEE Trans. on Power Systems_ , vol. 35, no. 6, pp. 4554–4563, Nov. 2020.
  * [11] A. Ghasemkhani, I. Niazazari, Y. Liu, H. Livani, V. A. Centeno, and L. Yang, “A regularized tensor completion approach for PMU data recovery,” _IEEE Trans. on Smart Grid_ , pp. 1–1, 2020.
  * [12] K. D. Jones, A. Pal, and J. S. Thorp, “Methodology for performing synchrophasor data conditioning and validation,” _IEEE Trans. on Power Systems_ , vol. 30, no. 3, pp. 1121–1130, May 2015.
  * [13] 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 Trans. on Smart Grid_ , vol. 10, no. 4, pp. 4554–4565, Jul. 2019.
  * [14] 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 Trans. on Power Systems_ , vol. 31, no. 2, pp. 1006–1013, Mar. 2016.
  * [15] G. M. De Mijolla, S. Konstantinopoulos, P. Gao, J. H. Chow, and M. Wang, “An evaluation of algorithms for synchrophasor missing data recovery,” in _Power Systems Computation Conference (PSCC)_ , Dublin, Ireland, 2018, pp. 1–6.
  * [16] Y. Hao, M. Wang, J. H. Chow, E. Farantatos, and M. Patel, “Modelless data quality improvement of streaming synchrophasor measurements by exploiting the low-rank hankel structure,” _IEEE Trans. on Power Systems_ , vol. 33, no. 6, pp. 6966–6977, Jun. 2018.
  * [17] A. Damen, P. Van den Hof, and A. Hajdasinski, “Approximate realization based upon an alternative to the hankel matrix: the page matrix,” _Systems & Control Letters_, vol. 2, no. 4, pp. 202 – 208, Dec. 1982.
  * [18] S. Chatterjee, “Matrix estimation by universal singular value thresholding,” _Ann. Statist._ , vol. 43, no. 1, pp. 177–214, Feb. 2015.
  * [19] W. Palma, _Time Series Analysis_. Hoboken, NJ: Wiley, 2016.
  * [20] Y. Luo, X. Cai, Y. Zhang, J. Xu, and X. Yuan, “Multivariate time series imputation with generative adversarial networks,” in _Proceedings of the 32nd International Conference on Neural Information Processing Systems_. Red Hook, NY, USA: Curran Associates Inc., 2018, p. 1603–1614.
  * [21] J. T. Connor, R. D. Martin, and L. E. Atlas, “Recurrent neural networks and robust time series prediction,” _IEEE Transactions on Neural Networks_ , vol. 5, no. 2, pp. 240–254, 1994.
  * [22] A. Agarwal, M. J. Amjad, D. Shah, and D. Shen, “Model agnostic time series analysis via matrix estimation,” _Proc. ACM Meas. Anal. Comput. Syst._ , vol. 2, no. 3, Dec. 2018. [Online]. Available: https://doi.org/10.1145/3287319
  * [23] A. Agarwal, A. Alomar, and D. Shah, “On multivariate singular spectrum analysis,” 2020.
  * [24] Y. Chen and Y. Chi, “Harnessing structures in big data via guaranteed low-rank matrix estimation: Recent theory and fast algorithms via convex and nonconvex optimization,” _IEEE Signal Processing Magazine_ , vol. 35, no. 4, pp. 14–31, 2018.
  * [25] L. T. Nguyen, J. Kim, and B. Shim, “Low-rank matrix completion: A contemporary survey,” _IEEE Access_ , vol. 7, pp. 94 215–94 237, Jul. 2019\.
  * [26] Fenghua Gao, J. S. Thorp, A. Pal, and Shibin Gao, “Dynamic state prediction based on auto-regressive (AR) model using PMU data,” in _IEEE Power and Energy Conference at Illinois_ , Champaign, IL, Feb. 2012, pp. 1–5.
  * [27] M. Gavish and D. L. Donoho, “The optimal hard threshold for singular values is $4/\sqrt{3}$,” _IEEE Trans. on Information Theory_ , vol. 60, no. 8, pp. 5040–5053, Aug. 2014.
  * [28] K. Chatterjee, N. R. Chaudhuri, and G. Stefopoulos, “Grouping PMU signals for guaranteed recovery under corruption: Insights and recommendations,” in _IEEE PES General Meeting (PESGM)_ , 2020, pp. 1–5.
  * [29] M. Naglic, “PMU measurements of IEEE 39-bus power system model,” 2019. [Online]. Available: https://dx.doi.org/10.21227/vkz3-2e96
  * [30] V. Venkatasubramanian, “Real-time strategies for unwrapping of synchrophasor phase angles,” _IEEE Trans. on Power Systems_ , vol. 31, no. 6, pp. 5033–5041, 2016.
  * [31] M. Brown, M. Biswal, S. Brahma, S. J. Ranade, and H. Cao, “Characterizing and quantifying noise in PMU data,” in _IEEE PES General Meeting (PESGM)_ , 2016, pp. 1–5.