1. Introduction

This paper summarizes our recent progress on missing data recovery by exploiting low-dimensional structures in high-dimen- sional Phasor Measurement Unit (PMU) datasets. In previous study ⁽¹⁾ a missing data recovery method has been developed, referred to as OLAP(Online Analytical Processing), to reconstruct missing PMU data and obtain promising results on actual PMU datasets. This study extends our previous efforts in this line of research. Specially, data recovery under extreme conditions is considered, e.g., almost all the measurements in all PMU channels are lost simultaneously. Existing low-rank-based methods cannot handle the pattern of simultaneous data loss. It is demon- strated that not only the matrix of spatial-temporal PMU data is low rank, but also the Hankel matrix of the PMU data is low rank. This property of low-rank Hankel matrix results from the power system dynamics and does not hold for general low-rank matrices. The central idea of this study is to further exploit the low-rank property of the Hankel matrix to recover missing points under extreme conditions and correspondingly modify and enhance the existing OLAP method, resulting in a new method named as OLAP-H. OLAP-H has been tested numerically using a KEPCO PMU dataset. In order to simulate extreme conditions, four models of missing data patterns are proposed and the performance of OLAP-H with OLAP was compared. Numerical experiments demonstrate that the recovery accuracy is improved when OLAP-H is employed. This paper is structured as follows. Section 2 describes the low-rank property of the spatial-temporal blocks of PMU data. The low-rank property is the key to missing data recovery methods. Section 3 presents the research motivation. Section 4 describes the Hankel matrix and its low-rank property. Section 5 presents the developed method OLAP-H to conduct online missing PMU data recovery. Section 6 records the simulation results with different missing data modes. Conclusions are provided in Section 7.

2. Low-rank property of PMU data

The PMU data used in this paper are from five 345 KV substations recorded by KEPCO(Fig. 1). The recording rate of PMUs is 60 samples per second. With the provided measurements during 10 minutes, we conduct pre-processing and obtain the measurements during the common period, which starts from 2016-02-27 05:15:30:967 and ends at 05:25:29:983. The data of SMITH is excluded in the analysis because it is misaligned with other PMUs in time synchronization. Fig. 1shows the geographical location of six substations at which the selected PMUs are installed in KEPCO system and demonstrates the recorded frequencies of all six substations. The data of SMITH leads the data of the remaining five substations about 24s, while they share the same time tags in the provided data files. One can also see that from the voltage magnitudes in Fig. 2. Clearly, SMITH does not record the disturbance because of the time differences. Therefore, the data of substation SMITH are excluded in the remaining part of the paper.

The data of voltage angles of the remaining five substations are shown in Fig. 3. Since the event occurs at about 25s from the beginning of 05:15:30:967, in this paper, we focus on the data from the beginning to 40s that contains the event.

Let a 5 by 2400 complex matrix M contain the PMU data in the first 40 seconds. Each row corresponds to a sequence of voltage phasor measurements of one substation. Each column corresponds to the PMU measurements at the same sampling instant. Fig. 4shows the singular values of . The five singular values are 38241.1, 50.0, 31.7, 10.6, 4.1.

If one wants to approximate matrix with a rank-r matrix with a minimum approximation error, the problem can be formulated as:

$$ \begin{array}{c} \min _{\hat{M}}\|M-\hat{M}\|_{F} \\ \text {s.t. } \operatorname{rank}(\hat{M})=r \end{array} $$ Let denote singular value decomposition of matrix M, where denotes the conjugate transpose of complex matrix V. Let denote the r dominant left singular vectors. Let denote the diagonal matrix with only the r dominant singular values in the diagonal. contains the r dominant right singular vectors. Then, the solution to the above optimization problem is . The approximation error is defined as follows:

$$ e_{a}=\frac{\|M-\hat{M}\|_{F}^{2}}{\| M_{F}^{2}} $$ Based on the analysis, it was found that the KEPCO data matrix M can be approximated by a rank-1 matrix with an approximation error less than 0.01%. Thus M can be approximated well with a rank-1 matrix. The low-dimensionality of PMU measurements can be employed to conduct PMU data compression ⁽²⁾, system event detection ^(3,6), missing data recovery ^(1,5), cyber data attack detection ⁽⁴⁾, etc. This study focuses on online missing data recovery. Examples of existing methods on online subspace tracking and missing data recovery include Parallel Estimation and Tracking by Recursive Least Squares (PETRELS)⁽²⁾ and Grassmannian Rank-One Update Subspace Estimation (GROUSE)⁽³⁾. These two methods assume the dimension of the subspace is known and fixed. The dimension of the subspace, however, usually changes when a disturbance happens in power systems. In past study an online algorithm for PMU data processing (OLAP) ⁽¹⁾ was proposed to track the underlying subspace with varying dimensionalities and fill in the missing observations.

3. Limitations of existing low-rank methods in PMU data recovery

One major advantage of low-rank-based methods is that they are data oriented and do not require the modeling of power systems. One limitation of these methods is that they degrade significantly under extreme conditions. For example, if almost all the PMU measurements are lost simultaneously at a given time instant, existing methods would either fail to recover or estimate the missing points by previous measurements. Moreover, because the measurements are usually noisy, the recovery error generally increases when the noise increases. In addition, due to a lack of regularization terms and the influence of the measurement noise, the recovered data always contains spikes, which can be shown in Fig. 5, where the data are erased randomly and the data loss percentage is 10%. Building on the above analysis, we propose to exploit the Hankel structure to conduct online missing data recovery so that missing points can be recovered accurately under extreme conditions. The proposed low-rank-based methods have been also enhanced by adding a low-pass filter to reduce the impact of noise on the data recovery.

4. Hankel matrix and the low-rank property

Let a vector in denote the measurements at instant, where m is the number of PMU measurement channels. The measurement matrix from time 1 to is represented by

$$ Y(1, t)=\left[\begin{array}{llll} y_{1} & y_{2} & \cdots & y_{t} \end{array}\right] $$ The rows and the columns of the measurement matrix are related the circuit laws and the dynamics of power systems, as is shown in Fig. 6.

Hankel matrices are usually employed to analyze the order and parameters of linear systems based on the obtained inputs and outputs. A Hankel matrix is defined as

$$ Y_{H}^{j}(1, t)=\left[\begin{array}{ccc} y_{1} & y_{2} & \cdots & y_{t-j+1} \\ y_{2} & y_{3} & \cdots & y_{t-j+2} \\ \vdots & \vdots & \ddots & \vdots \\ y_{j} & y_{j+1} & \cdots & y_{t} \end{array}\right] $$ where j denotes the number of observation vectors in each column of Hankel matrix. The size of the constructed Hankel matrix is . Using the KEPCO measurements during the first 40 seconds, the corresponding Hankel matrices were constructed and the approximation errors were computed with low-rank matrices, as is shown in Fig. 7. It can be observed that the constructed Hankel matrices are still low-rank with different choices of j. Note that when , the size of the constructed Hankel matrix almost increases linearly with j. Fig. 8shows the approximation errors if is approximated with a matrix of rank jr, where r is the rank of the approximated matrix to . A decreasing trend in approximation error with the increase of j can be observed in Fig. 8. That means the approximate rank does not increase proportionally with the size of Hankel matrix of the PMU data, which in turn demonstrates that the Hankel structure captures the additional dynamic information in the datasets compared with the original data matrix.

It should be noted that the low-rank property of the Hankel matrix results from the dynamics of power systems. This property does not hold for general low-rank matrices. In order to further illustrate this concept, the columns of are randomly permutated, and let denote the matrix after permutation. One can easily check that has the same rank as and thus, is still low-rank. However, since the columns in do not correspond to consecutive measure- ments, the corresponding Hankel matrix may no longer be low-rank, as is shown in Fig. 9. From the results shown in Fig. 7and Fig. 9, it can be observed that a rank-k approximation matrix to that always achieves a smaller error. For instance, if k=3 and j=5, the approximation error to is while the approximation error to is 0.233. Existing recovery methods are based on the low-rank property of the original observation matrix. The novel idea of this study is to further exploit the low-rank property of the Hankel matrix to improve the accuracy of data recovery. Numerical experiments using the KEPCO data show that the new method performs well under extreme conditions.

5. OLAP-H: an online missing data recovery method

The procedures of the original OLAP are summarized as follows:

1)With the past L output vectors at current instant t, conduct SVD, determine the approximate rank r and obtain the r dominant left singular vectors .

2)Receive new data with missing points, compute , where denotes the support set of the observed entries in .

3)Fill in the missing entries of with , where denotes the complement set of . Next, OLAP-H is described to conduct online missing data recovery based on the low-rank property of the Hankel matrix. The developed method has advantages over OLAP in the following aspects:

a)It can fill in the missing entries with a higher accuracy under extreme conditions by exploiting the low-rank property of the Hankel matrix;

b)The algorithm reduces the influence of measurement noise on missing data recovery performance by (1) adding one regularization term and (2) assigning different weights to PMU channels with different noise levels when determining the coefficients of the method.

The procedures of the developed algorithm are as follows:

Input: Approximation error threshold e_a, coefficient $\lambda \in[0,1]$, number of observation vectors in each column of Hankel matrix $j$, window length $L$.

For $t=1,2,3, \cdots$ do

1.Construct Hankel matrix $Y_{H}^{j}(t-L, t-1)$ from the past L observation vectors, conduct SVD on $Y_{H}^{j}(t-L, t-1)$, $$ Y_{H}^{j}(t-L, t-1)=U \Sigma V^{*}, $$ where $x^{*}$ denotes the conjugate transpose of complex vector or matrix $x$.

2. Find the smallest r satisfying $$ \frac{\|\Sigma-\Sigma\|_{F}^{2}}{\|\Sigma\|_{F}^{2}} \leq e_{a}, $$ then determine the corresponding r dominant left singular vectors $U^{r}$.

3. Receive new data $y_{t} \in C^{m \times 1}$ with erasures, construct a new column of Hankel matrix $\beta_{H}^{t}=\left[\begin{array}{lll}y_{t-j+1}^{*} & y_{t-j+2}^{*} & \cdots y_{t}^{*}\end{array}\right]^{*}$, the support set of observed entries in $\beta_{H}^{t}$ is denoted as $\Psi_{t}$.

4. Compute

where represents the assigned weight for each PMU, $w$ is smaller for a PMU with a higher noise level; $P_{\Psi_{t}}$ is a diagonal matrix, if, $P_{\Psi_{t}}(i, i)=1$ otherwise $P_{\Psi_{t}}(i, i)=0$.

5. Fill in the missing entries of $y_{t}$ with the corresponding entries in $U^{r} \tilde{v}$.

End for

Note that if the rank of underlying subspace is assumed to be fixed as r and known, then step 1 and step 2 can be simplified to conduct SVD (Singular Value Decomposition) on $Y_{H}^{j}(t-L, t-1)$ to obtain the $r$ dominant left singular vectors $U^{r}$.

In practice, $w$ can be estimated as $$ w=\left[\begin{array}{cccc} \sigma_{1}^{-1} & 0 & \cdots & 0 \\ 0 & \sigma_{2}^{-1} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \sigma_{m}^{-1} \end{array}\right], $$ where $\sigma_{i}$ denotes the standard variation of the recorded voltage magnitude of PMU i when no event occurs. In this case, the measurements with higher levels of noise are assigned to smaller weights. That allows a larger gap between the observed data and the recovered data, so as to reduce the impact of noise on determining the coefficient $\tilde{v}$ and the recovered data $U^{r} \tilde{v}$.

It is noticed that due to the high-rate sampling of PMUs, the measurement at time $t$ is not far away from that at $t$-1. In addition, when applying the original OLAP algorithm to recover the missing data, there may exist spikes in the recovered data. In step 4, the first term $\left|P_{\bar{N}_{\mathbf{e}}} w\left(\beta_{H}^{t}-U^{r} v\right)\right|_{F}$ and the second term $\left|P_{\Psi_{t-1}} w\left(\beta_{H}^{t-1}-U^{\tau} v\right)\right|_{F}$ are used to evaluate the distance of the recovered data $U^{r} v$ to the observed data in $\beta_{H}^{t}$ and $\beta_{H}^{t-1}$, respectively. Taking both terms into account, the recovered data $U^{r} v$ are required not only to fit well with the observed data in $\beta_{H}^{t}$ but also to stay close to $\beta_{H}^{t-1}$. In this case, the second term in (1) serves as a low-pass filter with the aim to smooth the recovered data.

6. Simulation Results

7. Conclusions

In this study, the limitations of existing modeless data recovery methods were analyzed, including their degraded performance in missing data recovery under extreme conditions, and the influence of measurement noise on the recovered data was discussed. In order to handle these limitations and reduce the impact of measurement noise, the low-rank property of the Hankel matrix was studied and an improved algorithm, OLAP-H, was developed to recover missing data in real-time. To summarize, the following accomplishments have been achieved:

1.The low-rank properties of original measurement matrix and the constructed Hankel matrix are studied with the KEPCO PMU dataset.

2.A new method OLAP-H is developed based on the original OLAP method to exploit the low-rank property of Hankel matrix to recover the missing data. In order to reduce the influence of noise on the recovered data, PMU channels with higher noise levels were penalized in the data recovery. In addition, the recovered data were smoothened by adding a regularization term of the change in the measurements.

3.In order to simulate the extreme conditions where the performance of existing methods degrades significantly, four modes of data loss were considered. The developed algorithm is applied to fill in unobserved data in all four modes. Numerical experiments show that the recovery performance on the Hankel matrix is much better than the performance on the original measurement matrix.