1

Robust Recovery of Missing Data in Electricity
Distribution Systems

arXiv:1708.01583v1 [cs.SY] 4 Aug 2017

Cristian Genes, Iñaki Esnaola, Samir M. Perlaza, Luis F. Ochoa, and Daniel Coca.

Abstract—The advanced operation of future electricity distribution systems is likely to require significant observability of the
different parameters of interest (e.g., demand, voltages, currents,
etc.). Ensuring completeness of data is, therefore, paramount. In
this context, an algorithm for recovering missing state variable
observations in electricity distribution systems is presented. The
proposed method exploits the low rank structure of the state
variables via a matrix completion approach while incorporating
prior knowledge in the form of second order statistics. Specifically, the recovery method combines nuclear norm minimization
with Bayesian estimation. The performance of the new algorithm
is compared to the information-theoretic limits and tested trough
simulations using real data of an urban low voltage distribution
system. The impact of the prior knowledge is analyzed when a
mismatched covariance is used and for a Markovian sampling
that introduces structure in the observation pattern. Numerical
results demonstrate that the proposed algorithm is robust and
outperforms existing state of the art algorithms.
Index Terms—recovery of missing data, distribution systems,
matrix completion, Bayesian estimation

I. I NTRODUCTION

T

HE wide-spread adoption of residential scale low carbon technologies, from PV systems to electric vehicles,
will undoubtedly bring technical challenges to the electricity
distribution systems as they have been designed for passive
loads [1] and [2]. Therefore, and as part of the Smart Grid
vision, electricity distribution systems, including low voltage
circuits, are likely to adopt a more active role so as to costeffectively manage controllable network elements as well as
participants [3]. As a result, monitoring and control procedures
are expected to face increasingly demanding performance requirements posed by the dynamic and unknown scenarios that
the smart grid gives rise to. Advanced control strategies require
timely and accurate data describing the state of the grid. In
this setting, the sensing infrastructure is expected to provide
complete and reliable state information of the distribution
system. However, in practical scenarios, the operator faces
challenges like data injection attacks [4], [5] or missing data
Cristian Genes and Daniel Coca are with the Department of Automatic
Control and Systems Engineering, University of Sheffield, Sheffield S1 3JD,
UK.
Iñaki Esnaola is with the Department of Automatic Control and Systems
Engineering, University of Sheffield, Sheffield S1 3JD, UK, and also with
the Department of Electrical Engineering, Princeton University, Princeton NJ
08540, USA.
Samir M. Perlaza is with the Institut National de Recherche en Informatique
et Automatique (INRIA), Lyon, France, and also with the Department of
Electrical Engineering, Princeton University, Princeton NJ 08540, USA.
Luis F. Ochoa is with the Department of Electrical and Electronic Engineering, The University of Melbourne, Melbourne 3010, Australia, and
also with the School of Electrical and Electronic Engineering, The University of Manchester, Manchester M13 9PL, UK. (cgenes1@sheffield.ac.uk,
esnaola@sheffield.ac.uk, samir.perlaza@inria.fr, luis ochoa@ieee.org, and
d.coca@sheffield.ac.uk).

[6], [7]. Sensor failures, unreliable communication or data
storage issues are some of the causes for incomplete sets of
observations. As a consequence, the state of the grid is not
perfectly known and control mechanisms are difficult to implement. For instance, accurate measurements are necessary to
implement a centralized control scheme for voltage regulation
in distribution systems [8]. In view of this, it is vital to develop
estimation procedures for the missing data using the available
observations.
Missing data recovery can be cast as a minimum mean
square error (MMSE) estimation problem when a probabilistic
description of the underlying process governing the state
variables is available. However, the MMSE estimation relies
on accurate second order statistics which is an unrealistic
assumption in practical scenarios [7], [9]. The increased number of nonlinear loads and the turbulent nature of distributed
generation options in the locally controlled grid affects the
precision of the postulated statistics for the state variables.
For that reason, the efficiency of MMSE estimation is limited
in the smart grid context [7].
Matrix completion (MC) was recently proposed to recover
missing data from partial observations [10]. The main advantage is that the recovery via MC requires mild assumptions
about the setting, e.g. access to second order statistics is not
required. Instead, matrix completion-based recovery exploits
the fact that correlated state variable vectors give rise to
approximately low rank data matrices. That being the case, the
recovery of the missing entries of low rank matrices is feasible
in a convex optimization context provided that a sufficient
fraction of the entries is observed [10], [11], and [12]. The
key theoretical results therein are based on the assumption that
the locations of sampled entries are uniformly distributed. In
practice, however, this assumption is not always satisfied. For
instance, in electricity distribution systems missing data entries
tend to display significant structure across both space and
time [7]. The applicability of MC recovery for non-uniform
sampling is studied in [6], [13]. Not surprisingly, low rank
minimization tools are also used to address the problem of
electricity price forecasting [14] and to develop a framework
for efficient processing of synchrophasor data [15].
This paper compares the performance of different missing data recovery strategies with respect to the informationtheoretic limit and introduces a novel algorithm that addresses
the main shortcomings of existing techniques. The performance of the new algorithm is tested against singular value
thresholding (SVT) recovery [16] and MMSE estimation under
realistic assumptions, i.e., the postulated statistics are not
accurate and the sampling pattern is not uniform. In particular,
a mismatched covariance matrix model is used in [17] and

2

[18] for the case in which imperfect second order statistics
are available. Non-uniform sampling is considered to account
for structure on the observation pattern [15]. Numerical results
show a significant gain in performance for both cases when
compared to SVT recovery. Remarkably, the new algorithm is
robust to mismatched statistics and to non-uniform sampling
patterns.
II. S YSTEM MODEL
Consider a low voltage distribution system with L feeders. Each feeder includes a sensing unit that measures the
electrical magnitudes of operational interest at predetermined
time instants. These measurements that include phase active
power, phase reactive power and phase voltage support the
operator in controlling, monitoring, and managing the network.
In practice, the acquisition process provides the operator with a
noisy and incomplete set of state variables. For that reason, the
operator needs to recover the missing data using the available
observations.
A. Source Model for State Variables
(s)

For a given electrical magnitude s, let mi,j be the value of
that particular electrical magnitude at feeder i at time j. The
matrix of state variables for s, denoted by M(s) ∈ RN ×L ,
contains the aggregated state variable vectors from all feeders,
∆
(s)
(s)
(s)
i.e. M(s) = [m1 , m2 , ..., mL ]. The state variable vector
for s, in feeder i, contains the N state variables generated in
(s)
(s)
(s)
(s)
the feeder and is given by mi = [mi,1 , mi,2 , ..., mi,N ]T ∈
RN . Without loss of generality the analysis is carried out for
a particular electrical magnitude, and therefore, the index s
is dropped. The resulting data matrix M contains the state
variable of interest at time instants 1, 2, ..., N for all L feeders.
Real data is used to model the statistical structure of the data
generated in an electricity distribution system. The real data
set under consideration contains values from 200 residential
secondary substations across the North West of England collected from June 2013 to January 2014. The data collection is
part of the “Low Voltage Network Solutions” project run by
Electricity North West Limited [19]. Each substation creates a
daily file containing values of voltage, current and power levels
for all three phases. An analysis of the distribution and sample
covariance matrix of the voltage data set under consideration
is presented in [7]. Therein it is shown that state variables can
be modelled as a multivariate Gaussian random process
mi ∼N (µ, Σ),

(1)

and {mi }L
i=1 is an independent and identically distributed sequence of random variables. Consequently, M is a realization
of the random process describing the value of the state variable
of interest across the grid.
Fig. 1 describes the distribution system monitoring model.
In this setting, the electrical magnitudes describing the state
of the system are modelled as a random process that outputs
a realization M ∈ RN ×L every N time instants. The state
of the grid is fully described by the entries of the matrix
M. However, the operator observes a subset of the complete
set of state variables, i.e. measurements are lost during the
acquisition process. The aim of the estimation process is to
recover the missing entries.

Low Voltage
Distribution
System

M

R

Observation

P⌦ (R)

P⌦

Missing
Data
Recovery

c
M

N
Figure 1. Block diagram describing the system model.

B. Acquisition
The sensing infrastructure introduces additive white Gaussian noise (AWGN) as a result of the thermal noise present at
each sensor. The resulting measurements are given by
R = M + N,

(2)

(N)i,j ∼ N (0, σ 2 ),

(3)

where
for i ∈ {1, 2, ..., N } and j ∈ {1, 2, ..., L}. Moreover, it
is also assumed that only a fraction of the complete set
of measurements (entries in R) are communicated to the
operator. Denote by Ω the subset of observed entries, i.e.,
Ω ⊆ {1, 2, ...N } × {1, 2., ..., L}. By definition it follows that
Ω is given by
∆

Ω = {(i, j) : (R)i,j is observed}.

(4)

Formally, the acquisition process is modelled by the function
f : RN ×L → R|Ω| with f (M) = PΩ (R) where
PΩ (R) = (R)Ω ,

(5)

and |Ω| denotes the cardinality of Ω. The observations given
by (5) describe all the data that is available to the operator for
estimation purposes and therefore, the recovery of the missing
data is performed from the observations PΩ (R).
C. Estimation
The estimation process of the complete matrix of state
variables based on the available observations is modelled by
c = g(f (M))
the function g : R|Ω| → RN ×L . The estimate M
is obtained by solving an optimization problem based on an
optimality criterion. In this paper, the optimality criterion is
the mean square error (MSE) given by


E kM − g(f (M))k2F
,
(6)
MSE (M; g) =
NL
where k·kF denotes the Frobenius norm. The normalized mean
square error (NMSE) is defined as
NMSE (M; g) = MSE (M; g)

NL
.
kMk2F

(7)

For this optimality criterion, the optimal estimate of the
missing data is given by the MMSE estimate
c MMSE = E[M|f (M), Σ],
M

(8)

where Σ ∈ RN ×N is the covariance matrix defined in (1).
c MMSE
Note that, in general, obtaining the optimal estimate M
requires knowledge of the probability distribution describing
the state variables. If the state variables follow a Gaussian
distribution it boils down to knowledge of the second order
moments, i.e. the covariance matrix Σ which needs to be
known prior to the estimation process. In practice, the operator

3

relies on postulated statistics that typically do not match the
actual statistics. Consequently, the accuracy of the estimate is
a function of the difference between the real and the postulated
statistics.

III. I NFORMATION -T HEORETIC L IMIT
In order to assess the performance of the missing data
recovery techniques in absolute terms, this section introduces
the optimal performance theoretically attainable (OPTA) by
an estimator g when the state variables follow a multivariate
Gaussian distribution. For a given number of observations, the
minimum distortion achievable by any estimation method is
determined by the rate-distortion function [20]. In the electricity distribution setting described above, the observations
entries are corrupted by additive white Gaussian noise which
determines the finite rate at which information about the state
variables is obtained from the observations. Consequently, the
optimal performance is bounded by the capacity of the AWGN
channel
R(D) < C,

(9)

where R is the rate at which the source needs to be observed
to achieve a distortion D, and C is the capacity of the parallel
AWGN channels modelling the observation process. In view
of this, the OPTA for a multivariate Gaussian source is given
by
|Ω|
log (1 + snr),
(10)
R(D) ≤
2N L 10
where the signal to noise ratio, denoted by snr, is defined as
1
∆ N Tr(Σ)
,
σ2

snr =

(11)

Linear MMSE (LMMSE) estimation achieves the optimal
performance in the recovery of missing data for a given set of
observations Ω when the data is generated by a multivariate
Gaussian source and the optimality criteria is the MSE.
However, this estimation procedure relies on access to second
order statistics of the state variables. In particular, each column
of the matrix PΩ (R) is given by
PΩ (ri ) = Ai (mi + ni ),

where R is the source rate in nats/symbol, D is the mean
square error distortion per entry, λi is the i th largest eigenvalue of Σ, and θ is a parameter. The NMSE theoretically
attainable, NMSE(M; OPTA), follows from combining (7) and
(12) and is determined by
NL
.
kMk2F

(13)

IV. R ECOVERY OF MISSING DATA
In this section, the information-theoretic limit for missing
data recovery presented in Section III, is compared with
MMSE estimation and the singular value thresholding (SVT)
recovery.

(14)

where Ai is defined such that Ai mi = PΩ (mi ) and i ∈
{1, 2, ..., L}. Consequently, the LMMSE estimate for each
state variable vector is given by
b i = µ + Γi (PΩ (ri ) − Ai µ),
m

(15)

where µ is defined in (1) and

Γi = ΣATi (Ai ΣATi + σ 2 I)−1 ,

(16)

i ∈ {1, 2, ..., L}. The normalized error achieved by the
LMMSE estimator is given by
NMSE(M; LMMSE) =

c LMMSE k2
kM − M
F
,
kMk2F

(17)

c LMMSE = [m
b 1, m
b 2 , ..., m
b L ], with m
b i defined in (15).
where M
B. Singular Value Thresholding

Low rank matrices are recovered from a subset of the
entries via rank minimization techniques under mild coherence
conditions on the set of observations [10]. Specifically, the
missing entries are recovered by solving the following rank
minimization problem:
minimize
X

where σ 2 is defined in (3). The rate-distortion function of a
multivariate Gaussian process is computed using the following
parametric equations [21]
(
PN −1
R(θ) = N1 i=0 max(0, 12 log λθi )
(12)
PN −1
D(θ) = N1 i=0 min(θ, λi ),

NMSE(M; OPTA) = D

A. Minimum Mean Squared Error Estimation

rank(X)

subject to PΩ (X) = PΩ (M).

(18)

Unfortunately, this rank minimization problem is NP-hard.
Favorably, in [10] it is shown that when the entries on Ω
are sampled uniformly at random, the solution of the rank
minimization problem in (18) is obtained with high probability
by solving the nuclear norm minimization problem in (20).
SVT is a matrix completion algorithm [16] which produces
a sequence of matrices Xk that converges to the unique
solution of the following optimization problem:
1
τ kXk∗ + kXk2F
2
subject to PΩ (X) = PΩ (M),
minimize
X

(19)

where kXk∗ is the nuclear norm of the matrix X. Note that
when τ → ∞, the optimization problem in (19) converges to
the nuclear norm minimization problem proposed in [10]
minimize
X

kXk∗

subject to PΩ (X) = PΩ (M).

(20)

For large values of τ , SVT provides the solution to the
nuclear norm minimization problem. Compared to alternatives
like SeDuMi [22] or SDPT3 [23], SVT features a lower

4

computational cost per iteration. This is achieved by exploiting
the sparsity of Yk and the low-rank property of Xk to reduce
storage requirements. The low computational cost results in
the possibility of using larger matrices. Simulation results in
[16] show that SVT recovers matrices with nearly a billion
entries. In comparison, SeDuMi and SDPT3 produce accurate
solutions for squared matrices with dimension close to fifty.
In [24] the structure of the problem is exploited to reduce
the memory requirements and increase the matrix size up to
350. Because of the dimension of the data sets produced by
low voltage distribution systems, the remaining of the paper
focuses on the SVT as a benchmark MC-based recovery. The
main idea in SVT consists in the following iteration steps:
(
Xk = Dτ (Yk−1 ),
(21)
Yk = Yk−1 + δs PΩ (M − Xk ),
where Y0 = 0, δs is the step size that obeys 0 < δs < 2, and
the soft-thresholding operator, Dτ , applies a soft-thresholding
rule to the singular values of Yk−1 , shrinking these towards
zero. For a matrix Y ∈ RN ×L of rank r with singular value
decomposition given by
Y = USVT ,

S = diag({σi (Y)}1≤i≤r ),

(22)

where U and V are unitary matrices of size N × r and L × r,
respectively, and σi (Y) are the singular values of the matrix
Y, the soft-thresholding operator is defined as

10-2

10-4

10-6

10-8
0

0.2

0.4

C. Performance Evaluation with Real Data
This subsection presents a comparison between LMMSE
and SVT, and the theoretical limit, OPTA, using real electricity
distribution system data. The test matrix, M, is a square
matrix of size 500, i.e. N = L = 500, and contains voltage
measurements covering the state of the grid for a period of 2
hours. Each column represents a different state variable vector
that describes the grid on a different day and for a different

0.8

1

Figure 2. Real data recovery performance using SVT, LMMSE estimation,
for different levels of mismatch, and the OPTA, when SNR = 20 dB.

feeder. The entries in Ω are sampled uniformly at random with
probability
1
P[(i, j) ∈ Ω] =
E[|Ω|],
(24)
NL
and the performance of the SVT-based recovery is defined in
terms of the NMSE given by
NMSE(M; SVT) =

∆

Dτ (Y) = UDτ (S)VT , with Dτ (S) = diag({(σi (Y)−τ )+ }),
(23)
where t+ = max(0, t). Interestingly, the choice of τ is
important to guarantee a successful recovery, since large values
guarantee a low-rank matrix estimate but for values larger
than max (σi (Y)) all the singular values vanish. In [16],
i
the proposed threshold is τ = 5N . However, simulation
results presented in [7] show that τ = 5N gives suboptimal
performance when the number of missing entries is large.
Unfortunately, finding the optimal threshold when the matrix
is sparse is still an open problem. In general, the value of
the threshold for soft-thresholding based recovery algorithms
is obtained via numerical optimization in [7] and [12]. The
same soft-thresholding operator, Dτ , is used in a different
framework for denoising [12], [25], and [26]. In this context,
the performance of the denoiser, measured in MSE, is estimated using Stein's unbiased risk estimate (SURE) [27]. In
[28] a closed-form expression for the unbiased risk estimate
is presented for the operator Dτ .

0.6

c SVT k2
kM − M
F
,
kMk2F

(25)

c SVT is the SVT estimate of M based on PΩ (R). Let
where M
γ be the expected value of the ratio of missing entries for the
matrix M, that is:
∆

γ =1−

1
E[|Ω|].
NL

(26)

Since the performance of the LMMSE estimator depends on
the covariance matrix Σ, a mismatched covariance matrix
model is introduced to account for the difference between
the postulated and actual statistics. Specifically, the postulated
covariance matrix is given by
Σ∗ = Σ +

1 kΣk2F
∆,
SMR k∆k2F

(27)

where Σ is the actual covariance matrix, ∆ = HHT with
H ∈ RN ×N , the entries of H are distributed as N (0, 1).
The strength of the mismatch is determined by the signal
to mismatch ratio (SMR), which is defined such that for
SMR = 1 the norm of the mismatch is equal to the norm
of the real covariance matrix, i.e., kΣk2F = kα∆k2F .
Fig. 2 shows the performance, measured in NMSE, for
the SVT-based recovery compared to the performance of the
LMMSE estimator when different levels of mismatch are
introduced and to the theoretical limit given by the OPTA.
Numerical results in this section are obtained for a signal to
∆
noise ratio value of SNR= 20 dB, where SNR = 10 log10 snr.
It can be seen that the performance of the SVT algorithm
is closer to the theoretical limit when the number of missing

5

entries is large. Interestingly, the LMMSE estimator gives better performance when SMR ≥ 100. However, when SMR=10
and γ ≤ 0.55 the SVT algorithm outperforms the LMMSE
estimator. Moreover, the SVT provides a better recovery for
SMR=1 for almost all values of γ. In view of this, the LMMSE
estimation requires accurate second order statistics to perform
competitively in this setting which is an unrealistic assumption
in a practical scenario. Moreover, the performance of the SVT
algorithm depends of the threshold τ [7] which is difficult to
optimize for this case.
V. M AIN R ESULT
This section introduces a novel algorithm for missing data
recovery that incorporates imperfect second order information
statistics. The new approach is based on the SVT algorithm
but it exploits the information about the second order statistics
to optimize the threshold τ at each iteration k.

for i ∈ {1, 2, ..., N } and j ∈ {1, 2, ..., L}. In this setting, the
SURE [27] is given by
min(N,L)

SURE(Dτ )(Z) = −

2
N LσZ

+

X

min(τ 2 , σi2 (Z))

i=1

2
+ 2σZ
div(Dτ (Z)),

where σi (Z) are is the i-th singular value of Z for i ∈
{1, 2, . . . , N }. A closed-form expression for the divergence
of this estimator is obtained in [28]. For the case in which
Z ∈ RN ×L the divergence is given by

min(N,L) 
X
(σi (Z) − τ )+
div(Dτ (Z)) =
I(σi (Z) > τ ) + |N − L|
σi (Z)
i=1
min(N,L)

+2

X

i6=j,i,j=1

σi (Z)(σi (Z) − τ )+
,
σi2 (Z) − σj2 (Z)

A. Soft-thresholding parameter
The main shortcoming of the SVT algorithm is the lack
of guidelines for tuning the threshold τ . Numerical results in
[7] show that the value 5N proposed in [16] is not optimal
for every scenario. In order to provide better recovery it is
essential to tune the value of τ for each iteration of the
algorithm. In SVT the soft-thresholding operator is applied
on a sparse matrix which increases the difficulty of the tuning
process.
B. Exploiting second order statistics
In order to overcome the limitation imposed by the sparse
structure of the matrix Yk , the proposed algorithm estimates
the missing entries prior to the soft-thresholding step. Thus, the
available prior knowledge is exploited to produce an estimate
of the entries not contained in Ω. In this case, at each iteration
k of the proposed algorithm the matrix Zk is computed as
Zk = Yk + Lk ,

(28)

where Yk is defined as in the SVT algorithm and Lk is the
LMMSE estimate given by
k
Lk = PΩc (µ) + ΣΩc Ω Σ−1
ΩΩ (PΩ (Y ) − PΩ (µ)),

(29)

c

where Ω is the set of observed entries, Ω is the set of missing
entries, ΣΩc Ω is the covariance matrix between the entries in
Ωc and the entries in Ω and ΣΩΩ is the covariance matrix
of the entries in Ω. In a nutshell, the unknown entries are
estimated using the LMMSE-based recovery at each iteration
k. The result is a complete matrix Zk for which the tuning of
the threshold is feasible.
C. Optimization of thresholding parameter
Using the main result in [28], the performance of the softthresholding operator can be estimated when the input matrix
accepts the following model
Z = M + W,

(30)

where the entries of W are
iid

2
(W)i,j ∼ N (0, σZ
),

(31)

(32)

(33)

when Z has no repeated singular values and is zero otherwise.
Therefore, combining (32) and (33) gives a closed-form expression for the performance of the soft-thresholding operator
2
for different values of τ and different noise levels σZ
.
2
The proposed algorithm approximates σZ with the weighted
2
sum of the noise in Ω and in Ωc . Consequently, σZ
k is
calculated as
kYk − PΩ (M)k2F + |Ωc |DLMMSE
,
(34)
NL
where DLMMSE represents the average noise per entry in Ωc .
The optimal threshold for the matrix Zk is denoted by τ∗k and
it is calculated using
2
σZ
k =

τ∗k = arg min SURE(Dτ )(Zk ),

(35)

τ

2
where σZ
k is given by (34). Therefore, the iterations of the
proposed algorithm are

k
k−1

),
X = Dτ (Z
k
k−1
(36)
Y =Y
+ δb PΩ (M − Xk ),

 k
k
k
Z =Y +L ,

where the Dτ is defined by (23) and the step size δb is similar
to the step size δs in the SVT algorithm. The initial conditions
are Z0 = 0, Y0 = 0 and τ = 0. The stopping criteria is
similar to the SVT algorithm, namely
kPΩ (Xk − M)kF
≤ .
kPΩ (M)kF

(37)

The main advantage of the proposed algorithm is that the
threshold is optimized at each iteration facilitated by the prior
knowledge incorporated into the structure of the algorithm.
First, an initial guess of the unavailable entries is formed,
at each iteration k, based on Yk and the covariance matrix
Σ. The results are aggregated in the matrix Zk which is
approximated by the model in (30). In this case, an estimate
2
of the noise level, σZ
k , is needed to compute the SURE.
The optimal value of τ for Zk is obtained by minimizing
SURE(Dτ )(Zk ) in (32). Admittedly, the optimization of the
threshold is only possible as long as second order statistics

6

Algorithm 1 Bayesian Singular Value Thresholding
Input: observations set Ω, and observed entries PΩ (R), mean
µ, covariance matrix Σ, step size δb , tolerance , and
maximum iteration count kmax
c BSVT
Output: M
1: Set Y 0 = 0
2: Set Z0 = 0
3: Set τ = 0
4: Set Ωc = {1, 2, ..., N } × {1, 2, ..., L} \ Ω
5: for k = 1 to kmax do
6:
Compute [U, S, V] = svd(Zk−1 )
Pmin(N,L)
7:
Set Xk = j=1
(max(0, σj (Zk−1 ) − τ )uj vj
k
8:
if kPΩ (X − M)kF /kPΩ (M)kF ≤  then break
9:
end if
10:
Set Yk = Yk−1 + δb PΩ (M − Xk )
k
11:
Set Lk = PΩc (µ) + ΣΩc Ω Σ−1
ΩΩ (Y − PΩ (µ))
k
k
k
12:
Set Z = Y + L
2
k
2
c
13:
Set σZ
k = (kY − PΩ (M)kF + |Ω |DLMMSE )/N L
k
14:
Set τ = arg min SURE(Dτ )(Z )
τ

15:
16:

end for
c BSVT = Xk
Set M

are available. Therefore, the new approach requires additional
knowledge that is not necessary when using the SVT algorithm. That being said, the SVT algorithm requires setting the
value for the threshold which in general is difficult to tune. The
same amount of prior knowledge, i.e., covariance matrix, is
required by the LMMSE estimator. Still, when the postulated
statistics are not accurate, the performance of the LMMSEbased recovery reduces by up to an order of magnitude in
NMSE (See Fig. 2). For the proposed algorithm, the tradeoff between the performance and the accuracy of the prior
knowledge is studied in Section VI.
VI. N UMERICAL A NALYSIS
This section analyzes the performance of the BSVT algorithm using the real data set presented in Section II-A. The data
matrix M, utilized to assess the performance of the proposed
algorithm, is the same used in Section IV-C and contains the
voltage measurements from the electricity distribution system.
Similarly, the simulations in this section assume a signal to
noise ratio value of SNR= 20 dB. Moreover, the performance
of the BSVT algorithm is also measured in terms of NMSE
given by
NMSE(M; BSVT) =

c BSVT k2
kM − M
F
,
kMk2F

(38)

c BSVT is the output of the BSVT recovery. The
where M
performance of each recovery technique is averaged over 20
realisations of Ω for each ratio of missing entries. Numerically,
the proposed algorithm is evaluated on three aspects. First, the
gain in performance for the optimized threshold is assessed.
The Section VI-A compares the performance of the SVT-based
recovery with the BSVT algorithm when accurate second
order statistics are available. Secondly, the robustness of the
BSVT recovery when perfect prior knowledge is not available

10-1

10-2

10-3
0

0.2

0.4

0.6

0.8

1

Figure 3. Real data recovery performance using SVT, LMMSE estimation
and BSVT for different levels of mismatch, when SNR = 20 dB.

is evaluated. A comparison between the SVT algorithm, the
LMMSE estimator and the BSVT recovery is presented for
different SMR values. The case in which perfect second-order
statistics are available is also included. Finally, the robustness of the BSVT recovery to different sampling patterns is
evaluated using Markov-chain-based sampling. The numerical
performance of the new algorithm is compared to the SVT
algorithm for the case in which the positions of the missing
entries are not uniformly distributed.
A. Performance of the optimized threshold
In this section, the performance of the new algorithm is
compared to the SVT-based recovery using same data matrix
M and the same sets of available entries, Ω, for a particular
ratio of missing entries γ as defined in (26). The positions of
the missing entries are sampled uniformly at random from the
set of all entries.
Fig. 3 depicts the performance of both algorithms when
applied in identical scenarios. Clearly, the optimized threshold
and the Bayesian estimation step increase the performance
of the algorithm when accurate second order statistics are
available. When the postulated statistics, i.e., those available
to the operator are identical to the real statistics, the BSVT
algorithm provides a better performance for all values of
γ. The gain in performance is larger when the ratio of
missing entries is smaller than 0.4. Interestingly, the boost in
performance is substantial in the region in which SVT is least
efficient when compared to the fundamental limit (See Fig.
2). However, in practical scenarios the postulated and actual
statistics are different. The impact of mismatched statistics is
considered in the following section.
B. Robustness with respect to mismatched statistics
In order to address the problem of missing data recovery
in a realistic scenario, a level of mismatch between the real
covariance matrix and the one available to the operator is
considered. The mismatch covariance matrix model presented
in (27) is used in this section to assess the sensibility of the

7

10-1

0

20

40
10-2
60

80

100

10-3
0

20

40

60

80

100

Figure 4. Positions of the observed entries, Ω, generated by the Markovian
model for a 100 × 100 matrix, when E[L0 ] = N and E[γ] = 0.8.

proposed algorithm to inaccurate prior knowledge. Hence, the
LMMSE estimator and the BSVT algorithm are compared in
the no-mismatch regime and for a SMR value of 100 and
10. The performance of the SVT-based recovery is included
as a benchmark for comparing rank minimization based approaches.
Fig. 3 depicts the performance of the different estimation methods when mismatched second order statistics are
available. Remarkably, the proposed algorithm is robust to
mismatch in the second order statistics. In contrast with the
LMMSE estimator, the performance of the BSVT algorithm
does not change significantly when mismatch occurs. Moreover, the BSVT algorithm gives better recovery than the SVTbased recovery in all mismatch regimes throughout the range
of γ. In comparison with the LMMSE estimation, the BSVT
algorithm performs better for SMR = 100 when γ ≤ 0.65.
Furthermore, for SMR = 10 the proposed approach is the
best performing recovery method for almost all values of γ.
In practical scenarios, when the mismatch regime is difficult to
establish, the choice between LMMSE and SVT is difficult to
make. BSVT is a robust alternative and gives better recovery
in a wide range of missing date regimes.

0

0.2

0.4

0.6

0.8

1

Figure 5. Real data recovery performance for the Markov-chain-based
sampling model, using SVT and BSVT for different levels of mismatch, when
E[L0 ] = N and SNR = 20 dB.

Figure 6. State diagram for the Markovian sampling model.

to the operator, or in state S2 in which case the entry is not
available. As before, the set Ω contains all the entries from
the matrix M that are available to the operator. In Fig. 6,
p1 is the transition probability from state S1 to S2 and p2 is
the transition probability from S2 to S1 . Hence, the expected
value of the ratio of missing entries is given by the steady state
probability of being in S2 . Consequently, the expected value
of the ratio of missing entries for the Markovian sampling
model is
p1
E[γ] =
.
(39)
p1 + p2
The expected number of consecutive missing entries, E[L0 ],
is:
n
X
p1
E[L0 ] =
l
(1 − p2 )l .
(40)
p1 + p2
l=0

C. Robustness with respect to different sampling patterns
The problem of recovering missing data when the subset of
missing entries is not uniformly sampled is addressed in this
section. In practical scenarios, a sensor failure or a downtime
in the communication line provides the operator with a number
of consecutive unavailable measurements in the state variable
vectors. Let L0 be the number of consecutive missing entries.
The expected value of L0 varies depending on the reliability of
the sensing infrastructure. In the uniform sampling model this
scenario is not possible. In contrast, a more general sampling
procedure is introduced.
The proposed sampling model is based on a two-state
Markov Chain. In this setting, for each entry (M)i,j of the
matrix M, the finite state machine depicted in Fig. 6, is
either in state S1 in which case the entry (i, j) is available

Solving (40) for n → ∞ and combining with (39) leads to
E[L0 ] =

1 − E[γ]
(1 − p2 ).
p22

(41)

Therefore, for any given γ and L0 , using (39) and (41), p1
and p2 are identified such that on average the sampling model
in Fig. 6 has a ratio of missing entries γ and the length of
the vectors with consecutive missing entries L0 . Note that
the case E[L0 ] = 1 reduces to the uniform sampling model
with probability P[(i, j) ∈ Ω] = 1 − γ. In this framework, a
comparison between the SVT and the BSVT-based recoveries
is presented for the case in which the sampling pattern is not
uniform. In order to consider the case in which a particular
feeder does not provide any measurements, the expected length
of the vectors with missing data is selected to be equal to

8

the length of the state variable vectors, i.e., E[L0 ] = N . Fig.
4 shows an example of a sampling pattern generated by the
Markov-chain-based model, when E[L0 ] = N and E[γ] = 0.8.
Fig. 5 compares the performance of the SVT-based recovery
with the BSVT-based recovery for the case in which the matrix
M is sampled using the Markov-chain-based sampling model
with E[L0 ] = N . Different levels of mismatch are introduced
to assess the robustness of the new algorithm to mismatched
prior knowledge when the sampling pattern is not uniform.
Remarkably, the performance of the proposed approach is not
significantly affected by the amount of prior knowledge in
any of the missing data regimes. Moreover, BSVT performs
better than SVT when the sampling pattern is not uniform. A
significant gain in performance is observed for small values of
γ. Consider the following example for the sake of discussion,
for a fixed tolerance of 10−2 in NMSE, the SVT algorithm
recovers up to 4% of the entries of the matrix M while
BSVT recovers 40% (See Fig. 5). The improvement in the
data recovering performance for the same level of tolerance is
significant. Numerical results in this section show that BSVT
is not only providing better performance than SVT when the
entries are not uniformly sampled but it is also robust to
mismatched statistics. The robustness of the new algorithm
extends to different sampling patterns. In view of this, BSVT
represents a better alternative for recovering missing data in
practical scenarios than SVT and LMMSE estimation.
VII. C ONCLUSION
A novel algorithm for recovering missing data in data sets
that admit a low rank description has been presented. The
proposed approach, BSVT, combines the low computational
cost of SVT with the optimality of the LMMSE estimator
when the data source is modelled as a multivariate Gaussian
random process and second order statistics are available. The
combined new approach addresses the issues of individual recovery methods. The robustness of the new algorithm on both
mismatched statistics and sampling patterns was demonstrated
through simulations. In respect to the SVT algorithm the new
approach addresses the issue of choosing the value of τ by
calculating the optimal threshold at each iteration. Compared
with the standard LMMSE estimator the new algorithm is
robust to inaccurate second order statistics. In order to assess
practical scenarios, a sampling model that incorporated missing state variable vectors, is illustrated. The performance gain
compared to SVT was significant for both uniform and nonuniform sampling models. Ultimately, the proposed algorithm
is shown to provide a robust and low complexity method to
recover missing data in low voltage distribution systems.
R EFERENCES
[1] A. Navarro-Espinosa and L. F Ochoa, “Probabilistic impact assessment
of low carbon technologies in LV distribution systems,” IEEE Trans.
Power Syst., vol. 31, no. 3, pp. 2192–2203, May 2016.
[2] R.A. Walling, R. Saint, R. C. Dugan, J. Burke, and L.A. Kojovic,
“Summary of distributed resources impact on power delivery systems,”
IEEE Trans. Power Del., vol. 23, no. 3, pp. 1636–1644, Jul. 2008.
[3] C. Long and L. F. Ochoa, “Voltage control of PV-rich LV networks:
OLTC-fitted transformer and capacitor banks,” IEEE Trans. Power Syst.,
vol. 31, no. 5, pp. 4016–4025, Sep. 2016.

[4] M. Ozay, I. Esnaola, F. T. Y. Vural, S. R. Kulkarni, and H. V. Poor,
“Sparse attack construction and state estimation in the smart grid:
Centralized and distributed models,” IEEE Journal on Selected Areas
in Communications, vol. 31, no. 7, pp. 1306–1318, Jul. 2013.
[5] O. Kosut, L. Jia, R. J. Thomas, and L. Tong, “Malicious data attacks
on smart grid state estimation: Attack strategies and countermeasures,”
in Proc. of the First IEEE International Conference on Smart Grid
Communications, Oct. 2010, pp. 220–225.
[6] 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. Power Syst.,
vol. 31, no. 2, pp. 1006–1013, Mar. 2016.
[7] C. Genes, I. Esnaola, S. M. Perlaza, L. F. Ochoa, and D. Coca, “Recovering missing data via matrix completion in electricity distribution
systems,” in Proc. of the IEEE International Workshop on Signal
Processing Advances in Wireless Communications, Jul. 2016, pp. 1–6.
[8] Y. Isozaki, S. Yoshizawa, Y. Fujimoto, H.i Ishii, I. Ono, T. Onoda,
and Y. Hayashi, “Detection of cyber-attacks against voltage control
in distribution power grids with PVs,” IEEE Trans. Smart Grid, vol. 7,
no. 4, pp. 1824–1835, 2016.
[9] I. Esnaola, A. Tulino, and J. Garcia-Frias, “Linear analog coding of
correlated multivariate Gaussian sources,” IEEE Trans. Commun., vol.
61, no. 8, pp. 3438–3447, Aug. 2013.
[10] E. J. Candès and B. Recht, “Exact matrix completion via convex
optimization,” Communications of the ACM, vol. 55, no. 6, pp. 111–119,
2012.
[11] E. J. Candès and T. Tao, “The power of convex relaxation: Near-optimal
matrix completion,” IEEE Trans. Inf. Theory, vol. 56, no. 5, pp. 2053–
2080, May 2010.
[12] E. J. Candès and Y. Plan, “Matrix completion with noise,” Proceedings
of the IEEE, vol. 98, no. 6, pp. 925–936, Jun. 2010.
[13] R. Meka, P. Jain, and Inderjit S. D., “Matrix completion from power-law
distributed samples,” in Proc. of the Advances in Neural Information
Processing Systems 22, Dec. 2009, pp. 1258–1266.
[14] V. Kekatos, Y. Zhang, and G. B. Giannakis, “Electricity market
forecasting via low-rank multi-kernel learning,” IEEE J. Sel. Topics
Signal Process., vol. 8, no. 6, pp. 1182–1193, Dec. 2014.
[15] P. Gao, M. Wang, S. Ghiocel, J. Chow, B. Fardanesh, and G. Stefopoulos,
“Missing data recovery by exploiting low-dimensionality in power
system synchrophasor measurements,” IEEE Trans. Power Syst., vol.
31, no. 2, pp. 1006–1013, 2016.
[16] J.-F. Cai, E. J. Candès, and Z. Shen, “A singular value thresholding
algorithm for matrix completion,” SIAM Journal on Optimization, vol.
20, no. 4, pp. 1956–1982, Mar. 2010.
[17] I. Esnaola, A. M. Tulino, and H. V. Poor, “Mismatched MMSE
estimation of multivariate Gaussian sources,” in Proc. of the IEEE
International Symposium on Information Theory, Jul. 2012, pp. 716–
720.
[18] S. Verdú, “Mismatched estimation and relative entropy,” IEEE Trans.
Inf. Theory, vol. 56, no. 8, pp. 3712–3720, Aug. 2010.
[19] Electricity North West Limited, “Low voltage network solutions,”
[Online]. Available: http://www.enwl.co.uk/lvns.
[20] T. Cover and J. Thomas, Elements of information theory, John Wiley
& Sons, 2012.
[21] A. Kolmogorov, “On the Shannon theory of information transmission
in the case of continuous signals,” IRE Transactions on Information
Theory, vol. 2, no. 4, pp. 102–108, Dec. 1956.
[22] J.F. Sturm, “Using SeDuMi 1.02, a MATLAB toolbox for optimization
over symmetric cones,” Optimization Methods & Software, vol. 11, no.
1-4, pp. 625–653, 1999.
[23] K.C. Toh, M.J. Todd, and R.H. Tutuncu, “SDPT3 - A MATLAB software
package for semidefinite programming, version 1.3,” Optimization
Methods & Software, vol. 11, no. 1-4, pp. 545–581, 1999.
[24] Z. Liu and L. Vandenberghe, “Interior-point method for nuclear norm
approximation with application to system identification,” SIAM Journal
on Matrix Analysis and Applications, vol. 31, no. 3, pp. 1235–1256,
Nov. 2009.
[25] D. Donoho and M. Gavish, “Minimax risk of matrix denoising by
singular value thresholding,” Ann. Statist., vol. 42, no. 6, pp. 2413–
2440, Dec. 2014.
[26] X. Jia, X. Feng, and W. Wang, “Rank constrained nuclear norm
minimization with application to image denoising,” Signal Processing,
vol. 129, pp. 1–11, Dec. 2016.
[27] C. M. Stein, “Estimation of the mean of a multivariate normal
distribution,” The Annals of Statistics, pp. 1135–1151, Nov. 1981.
[28] E. J. Candès, C. A. Sing-Long, and J. D. Trzasko, “Unbiased risk
estimates for singular value thresholding and spectral estimators,” IEEE
Trans. Signal Process., vol. 61, no. 19, pp. 4643–4657, Oct. 2013.

