The wide-area measurement system (WAMS) using phasor measurement units (PMUs) has been regarded as one of the key enabling technologies in monitoring, control, and protections of the next-generation power grids 
. With continuous increase in PMU deployment and the resulting explosion in data volume, the design and deployment of an efficient wide area communication and computing infrastructure, especially from the point of view of resilience against a large number of missing data, is evolving as one of the greatest challenges to the power system and IT communities. With thousands of networked PMUs being scheduled to be installed in the United States by 2020, exchange of synchrophasor data between balancing authorities for any type of wide-area control will involve an enormous number of data flow in real-time per event, thereby opening up a wide spectrum of probabilities of data losses and data quality degradations in an unpredictable way. Data missing makes the system unobservable, degrades the performance of the state estimates, and weakens the security and stability of the system. Therefore, recovering missing PMU measurements has become a significant and inevitable problem in power systems.
PMU data can be structured as a matrix with each column and row representing the measurements of one channel and sample instant, respectively. Since large amounts of PMU data exhibit heavily correlated property [3, 4, 5], the matrix is approximately low-rank, and the problem of recovering the missing PMU data can be formulated as a low-rank matrix-completion problem. Studies on matrix completion algorithms are extensive, including atomic decomposition of minimum rank approximation (ADMiRa) 
, singular value projection (SVP), information cascading matrix completion (ICMC) , among which nuclear-norm-regularized matrix approximation [9, 10, 11, 12] and maximum-margin matrix factorization (MMMF)  are widely adapted. Using nuclear-norm-regularized matrix approximation, a singular value threshold has to be designed which influences the estimate accuracy. Developing the nuclear-norm-regularized matrix approximation, an alternating direction method (ADM) is provided for solving the matrix completion problem [14, 15]
. However, the calculation of the singular value decomposition (SVD) in ADM approach increases the computational time and complexity. Based on MMMF, Jain et al. and Hardt  proposed alternating least squares (ALS) schemes for solving the matrix-completion problem. Further, softImpute-ALS is provided for reducing the computational complexity . Gao et al. applied the MMMF approach on recovering the missing PMU data  firstly. Most of the existing approaches rely on an estimation of the rank of the data matrix, which is typically unavailable and time variant in practice. Inaccurate estimation of introduces modelling errors in the matrix completion problem. The computational complexity is lower with a smaller . On the other hand cannot be too small for estimation accuracy. Therefore, design of an adaptive and scalable online algorithm of PMU data recovery is an open challenge.
Motivated by these insights, we develop an algorithm that can recover the missing PMU measurement with low computational complexity and less operating time. The fundamental set-up for this optimization was based on MMMF and alternating direction method of multipliers (ADMM) [20, 21, 22]. Firstly, the observed PMU data is structured as a matrix whose columns and rows represent the measurements from one channel and the same sampling instant, respectively. Then we formulate the data recovery as an optimization problem in which we minimize the rank of the estimated matrix while keeping elements in the same as the corresponding ones in if they are present. An ADMM algorithm is proposed to solve the optimization problem in an iterative way. In the update equations there is no matrix inverse computation, which immensely reduces the computational complexity. In addition, it is not necessary to estimate the rank of the original data matrix without missing elements, which significantly cuts down the influence of the uncertain factor into the performance. Furthermore, we consider the case of missing data from all PMU channels. In this case, all elements in one row of the observed matrix are missing. One efficient algorithm is presented to reshape the observed matrix, and the lost data from all the channels can be recovered using ADMM approach. We illustrate the results using simulations of the IEEE 68-bus system model.
Ii Problem Formulation
Persistent model is one simple and traditional method to recover the missing PMU data. It utilizes the temporal correlation of the PMU measurements to recover the lost data in one channel. However, if in the disturbance scenario the measurements in the same channel are missing during a long time, the recovery with persistent model is not an advisable choise. In this section, we process a spatial-temporal blocks of PMU data, present the low-rank property of PMU data, and formulate the data recovery as a matrix completion problem.
Ii-a Low-rank property of PMU measurements
Denote as the PMU measurement matrix without data missing. Each column and row correspond to a sequence of measurements of one PMU channel, and the PMU measurments at the same sampling instant, respectively. Due to the noise, all the singular values of are larger than zero. An approximating rank approach, referred to Frobenius norm proportion , is stated as follows.
where are the singular values of the matrix and , , is the proportion factor. in (1) denotes the approximate rank of the matrix. Since the PMU measurements of voltage or current phasors or magnitudes from different lines or buses are strongly correlated, the approximate rank of is much smaller than [3, 4, 5]. Due to the low-rank property of PMU data, missing PMU measurement estimation can be converted into a low-rank matrix completion problem.
Ii-B An ADMM based approach for PMU data estimation
Let and denote the observed PMU measurements with missing data and the recovered matrix, respectively. Since should be a low-rank matrix, the matrix completion problem is formulated as follows:
where denotes the Hadamard product, i.e., . is the structural identity with its entry defined as
where the nuclear norm is the sum of the singular values of .
Using MMMF to further change the optimization problem (4), let , in which and . Without loss of generality, we assume . Since is equivalent to with Frobenius norm , the optimization function is equivalent to
In the previous work [16, 17], people estimated the rank of , set and , and applied ALS to solve (5). The computational complexity is . If , the computational complexity is , which is a biquadrate function of . With smaller the computational complexity is reduced. However, the value of cannot be too small to guarantee the estimation accuracy. For reducing the influence of the uncertain factor into the performance, we set the sizes of matrices and only depend on the size of observed matrix . In addition, we apply the ADMM method to solve (5) in an iterative way using the Lagrangian multiplier approach.
The augmented Lagrangian for (5) can be formulated as
where and are the matrices of the primal variables, is the matrix of the dual variables or the Lagrange multipliers associated with (5), and denotes a penalty weight.
After some algebraic, the augmented Lagrangian can be rewritten as
The gradients of the augmented Largrangian in (7) with respect to and are respectively given by
The updates in Algorithm 1 requires no matrix inverse, and the computational complexity is , which is a quadratic function of . In addition, it is not necessary to estimate the rank of matrix , which reduces the influence of uncertain factor into the performance. Penalty weight denotes the step size of the dual variable update. In general, large results in fast convergent rate.
Compared to approaches like interpolations and persistent models, ADMM algorithm utilizes the spatial and temporal correlations of PMU data to improve accuracy. In the persistent model, it replaces the missing data by the previous available data point. The persistent method recovers the lost data only based on temporal correlation. If the data from one channel are missing during a long time, and if there exists a dynamic in the time, then the estimation using persistent method doubtlessly is a nightmare. On the other hand, based on the spatial correlation, the missing data can be recovered using ADMM. We will compare the estimates using ADMM and persistent model with IEEE 68-bus power system simulation in SubsectionIII-C.
Ii-C Special case: missing data from all the channels
The power system often suffers natural and artificial disturbances during operation. It is possible that the data from all the channels are missing simultaneously under communication failure. In this case, no existing algorithms can recover the missing data. For solving this problem, the observed matrix has to be reshaped to avoid some of its rows missing. Our goal is that the proportion of missing elements in one row of the reshaped observed matrix is as small as possible. Meanwhile the corresponding reshaped recovery matrix is still low-rank.
We provide an alternative method, called cut-column reshaping method (CCRM), for reshaping the observed matrix. Using CCRM each column with length is separated into shorter columns with a length of . Thus, the -by- matrix is reshaped to a -by- matrix, and the original column correlation is held. The length of the new column should be larger than the row length of the original matrix, i.e., . also satisfies that , where denotes the smallest integer number which is larger than . Thus the numbers of rows and columns of reshaped matrix are both larger than . Due to the size, the rank of is no more than . Using CCRM the rank of reshaped matrix will not be reduced by the new size. In addition, with holding the column correlation, CCRM minimizes the proportion of zero elements in one row of reshaped matrix.
Consider a simple example to illustrate the reshaping method. A -by- matrix can be expressed as:
whose fifth row is missing. Using CCRM with and matrix is reshaped into a -by- matrix:
Now for each column and row, not all measurements are missing. If and are strongly correlated, and , and , and and are strongly correlated in pairs. The ranks of matrice and are both no more than . The proportion of missing elements to the first row is ; while it is to the fifth row of . CCRM is illustrated in Algorithm 2.
The missing PMU measurements from all the channels can be recovered using ADMM in Algorithm 1 after reshaped matrix using CCRM in Algorithm 2. Notice that if all the elements in one column of the reshaped observed matrix are missing, they cannot be recovered using ADMM. The recovery accuracy using ADMM will be declined sharply, if the measurements in one channel are missing more than successive sampling instants. With less lost data, the recovery accuracy will be enhanced.
Iii Simulation results
The IEEE 68-bus system is used to carry out the simulation to verify the proposals. We build up a PMU measurement matrix whose column and row corresponding to a sequence voltage phasors on lines and the sampling instants, respectively. The simulated measurements are obtained using the power systems toolbox (PST) nonlinear dynamics simulation routine and the data file data16m.m . A three-phase fault is imposed at the line connecting buses and . The fault starts at s, and clears on bus at s and on bus at s. For approaching to the true measurements, we add white Gaussian noise () into the PMU data. The measurements are observed during s and there are samples in one second. The -by- matrix is with no missing measurements and its approximate rank is with in (1). To test the recovery accuracy of the presented ADMM algorithm, some observed data in is set to be lost. Since the PMU data are missing arbitrary and unpredictable, in this paper we consider two cases of missing data: (1) Missing data randomly. The delivery of PMU measurements from multiple remote locations of power grids to monitoring centers can result in the random unavailability of PMU measurements; (2) Missing data in all channels simultaneously. The transform link malfunctions may result in data missing in all channels. We choose the penalty weight using ADMM, and the dual parameter and the estimated rank of filled completion matrix using ALS for comparison. In the paper, the computational time is obtained by operating Matlab programming.
Iii-a Case 1: Missing data randomly
In this case, we assume an independent and identical distribution (i.i.d) of the missing rate. For each data point, with a probability the measurement is missing and set to zero in artificially. Notice that it is different from the data which is equal to zero. If the actual data is zero, the corresponding element in is equal to . While if the data is missing, the corresponding element in is equal to . Table I compares some properties of ALS and ADMM in Case 1.
|# iterations||time||Sensitivity of parameters|
Though the number of convergence iterations using ADMM is larger than the one using ALS, the computational time using ADMM is less than , which is much smaller than using ALS.
Fig. 1 shows the statistic, maximum, and minimum values of mean absolute errors MAEs using ADMM and ALS with different observed data probabilities, respectively. The observed data probability denotes the likelihood of the observed data occurrence, i.e., .The statistic, maximum, and minimum values of MAEs are obtained by Monte Carlo method with independent times. With larger probability of observed measurements, MAE becomes smaller. The statistic values of MAEs using ADMM and ALS are close with each observed data probability. The difference between the maximum and minimum values of MAEs using ADMM is larger than the one using ALS.
Iii-B Case 2: Missing data in all channels
In this case, one row of data in matrix is lost. The -by- matrix which contains voltage phasor measurements can be treated as sub-matrices with a size of -by-. The observed data probability denotes the proportion of the observed sub-matrices to the total ones. For recovering the missing data in one row, firstly we reshape the observed matrix using CCRM. Since the rank of the orginal matrix is no more than , the number of the rows and columns of the reshaped matrix should be more than for avoiding reducing the rank artificially. In addition, with holding the column correlation, one purpose of reshaping is minimizing the proportion of zero elements in one row of the reshaped matrix. Thus using CCRM, the original -by- matrix is reshaped to a -by- matrix with . With in (1), the approximate rank of the reshaped observed matrix is . Since the size of the transposed reshaped observed matrix is similar to the observed matrix, the computational time using ADMM and ALS is similar to the results in Table I, respectively.
Fig. 2 shows the statistic, maximum, and minimum values of MAEs using ADMM and ALS with different observed data probabilities, respectively. The statistic values of MAEs using ADMM and ALS are still close. Compared with Case 1, the MAEs using both ADMM and ALS are larger. Though approximate rank of reshaped matrix is still 1, the minimum singular value becomes larger, whose influence into the recovery accuracy cannot be ignored. If the observed matrix is not reshaped, the missing row cannot be recovered using neither ADMM nor ALS, and the MAEs with different observed data probabilities are all around .
Iii-C Comparison among ADMM, ALS, and persistent model approaches
In the persistent model, it replaces the missing data at the sampling instant with the data at the if it is available. Only based on the temporal correlation of the PMU data, in a disturbed scenario the data which are lost during several successive sampling instants cannot be recovered successfully using persistent method. In this subsection, we let the data be lost from the sampling instant to the sampling instant on lines. Fig. 3 shows the estimated measurements using ADMM, ALS, and persistent methods from sampling instant to on Line 1. The blue line shows the values of actual measurements. Using both ADMM and ALS approaches, estimated measurements are close to the actual one. While using the persistent model, the estimate deviates from the actual data due to the dynamics in the measurements.
In this paper, we presented ADMM algorithm for missing PMU measurement recovery. We illustrated our results with noisy measurements from the IEEE 68-bus power system model. Compared with the ALS algorithm, the computational complexity and operating time are much smaller using the ADMM algorithm. In addition, the ADMM algorithm avoids to estimate the rank of filled completion matrix, which reduces the influence of the uncertain factor into the performance. We also consider the case of missing data in all the channels simultaneously and provide one approach to reshape the observed matrix for the recovery. Our future work in this area will include recovering continuous several rows of the observed matrix with all missing elements and testing the proposal using actual PMU data.
-  A. G. Phadke and J. S. Thorp, “Synchronized Phasor Measurements and Their Applications,” New York, Springer, 2008.
-  Y. Chen, L. Xie, and P. Kumar, “Dimensionality Reduction and Early Event Detection Using Online Synchrophasor Data,” in Proc. IEEE Power and Energy Society General Meeting, pp. 1-5, 2013.
-  N. Dahal, R. L. King, and V. Madani, “Online Dimension Reducion of Synchrophasor Data,” in Proc. IEEE PES Transmission and Distribution Conf. Expo. (T& D), pp. 1-7, 2012.
-  P. Gao, M. Wang, S. Ghiocel, and J. H. Chow, “Modeless Reconstruction of Missing Synchrophasor Measurements,” in Proc. IEEE PES General Meeting, pp. 1-5, 2014.
-  K. Lee and Y. Bresler, “ADMiRA: Atomic Decomposition for Minimum Rank Approximation,” IEEE Trans. Inf. Theory, vol. 56, no. 9, pp. 4402-4416, 2010.
-  P. Jain, R. Meka, and I. S. Dhillon, “Guaranteed Rank Minimization via Singular Value Projection,” in Adv. Neural Inf. Process. Syst., pp. 937-945, 2010.
-  R. Meka, P. Jain, and I. S. Dhillon, “Matrix Completion from Power-Law Distributed Samples,” in Adv. Neural Inf. Process. Syst., pp. 1258-1266, 2009.
-  E. J. Candes and T. Tao, “The Power of Convex Relaxation: Near-Optimal Matrix Completion,” IEEE Trans. on Information Theory, vol. 56, no. 5, pp. 2053-2080, 2010.
R. Mazumder, H. Trevor, and T. Robert, “Spectral Regularization Algorithms for Learning Large Incomplete Matrices,”
Journal of machine learning research, vol. 11, pp. 2287-2322, 2010.
-  J. Cai, E. J. Candes, and Z. Shen, “A Singular Value Thresholding Algorithm for Matrix Completion,” SIAM Journal on Optimization, vol. 20, no. 4, pp. 1956-1982, 2010.
-  S. Ma, D. Goldfarb, and L. Chen, “Fixed Point and Bregman Iterative Methods for Matrix Rank Minimization,” Mathematical Programming, vol. 128., no. 1, pp. 321-353, 2011.
-  J. D. Rennie and S. Nathan, “Fast Maximum Margin Matrix Factorization for Collaborative Prediction,” in Proc. the 22nd international conference on Machine learning(ACM), pp. 713-719, 2005.
-  C. Chen, B. He, and X. Yuan, “Matrix Completion via an Alternating Direction Method,” IMA Journal of Numerical Analysis, vol. 32, no. 1, pp. 227-245, 2012.
-  F. Xu, and P. Pan, “A New Algorithm for Positive Semidefinite Matrix Completion,” Journal of Applied Mathematics, vol. 3, pp. 1-5, 2016.
P. Jain, N. Praneeth, and S. Sujay, “Low-Rank Matrix Completion Using Alternating Minimization.” in
Proc. the forty-fifth annual ACM symposium on Theory of computing (ACM), 2013.
-  M. Hardt, “Understanding Alternating Minimization for Matrix Completion,” Foundations of Computer Science (FOCS), 2014 IEEE 55th Annual Symposium, 2014.
-  T. Hastie, R. Mazumder, J. D. Lee, and R. Zadeh, “Matrix Completion and Low-Rank SVD via Fast Alternating Least Squares,” J. Mach. Learn. Res, vol. 16, no. 1, pp. 3367-3402, 2015.
-  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, 2016.
-  R. Glowinski and A. Marroco, “Sur l’approximation, par èlèments finis d’ordre un, et la rèsolution, par pènalisation-dualitè d’une classe de problèmes de Dirichlet non linèaires.” ESAIM: Mathematical Modelling and Numerical Analysis - Mod lisation Math matique et Analyse Num rique 9.R2 pp. 41-76, 1975.
-  D. Gabay and B. Mercier, “A Dual Algorithm for the Solution of Nonlinear Variational Problems via Finite Element Approximation”. Computers & Mathematics with Applications, vol. 2, no. 1, pp. 17-40, 1976.
-  R. Glowinski and L. T. Patrick, “Augmented Lagrangian and operator-splitting methods in nonlinear mechanics,” Society for Industrial and Applied Mathematics, 1989.
-  X. Zhang, “Matrix Analysis and Applications,” Tsinghua and Springer Publishing house, Beijing, pp. 71-100, 2004.
-  M. Fazel, “Matrix Rank Minimization with Applications,” Ph. D. dissertation, Stanford Univ., Stanford, CA, 2002.
-  J. H. Chow and K. W. Cheung, “A Toolbox for Power System Dynamics and Control Engineering Education and Research,” IEEE Transactions on Power Systems, vol. 7, no. 4, pp. 1559-1564, 1992.