1

Robust and Scalable Power System State
Estimation via Composite Optimization

arXiv:1708.06013v1 [] 20 Aug 2017

Gang Wang, Student Member, IEEE, Georgios B. Giannakis, Fellow, IEEE, and Jie Chen, Senior Member, IEEE

Abstract—In today’s cyber-enabled smart grids, high penetration of uncertain renewables, purposeful manipulation of meter
readings, and the need for wide-area situational awareness, call
for fast, accurate, and robust power system state estimation. The
least-absolute-value (LAV) estimator is known for its robustness
relative to the weighted least-squares (WLS) one. However, due
to nonconvexity and nonsmoothness, existing LAV solvers based
on linear programming are typically slow, hence inadequate
for real-time system monitoring. This paper develops two novel
algorithms for efficient LAV estimation, which draw from recent
advances in composite optimization. The first is a deterministic
linear proximal scheme that handles a sequence of convex
quadratic problems, each efficiently solvable either via off-theshelf algorithms or through the alternating direction method
of multipliers. Leveraging the sparse connectivity inherent to
power networks, the second scheme is stochastic, and updates
only a few entries of the complex voltage state vector per
iteration. In particular, when voltage magnitude and (re)active
power flow measurements are used only, this number reduces
to one or two, regardless of the number of buses in the
network. This computational complexity evidently scales well
to large-size power systems. Furthermore, by carefully minibatching the voltage and power flow measurements, accelerated
implementation of the stochastic iterations becomes possible. The
developed algorithms are numerically evaluated using a variety
of benchmark power networks. Simulated tests corroborate that
improved robustness can be attained at comparable or markedly
reduced computation times for medium- or large-size networks
relative to the “workhorse” WLS-based Gauss-Newton iterations.

Index terms— SCADA measurements, cyberattacks, alternating direction method of multipliers, prox-linear algorithms.
I. I NTRODUCTION
The North American electric grid, the largest machine on
earth, is recognized as the greatest engineering achievement of
the 20th century [1]: thousands of miles of transmission lines
and millions of miles of distribution lines, linking thousands
of power plants to millions of factories and homes. Accurately
monitoring the grid’s operating condition is critical to several
control and optimization tasks, including optimal power flow,
reliability analysis, attack detection, and future network expansion planning [2].
The work was supported in part by NSF grants 1423316, 1442686, 1508993,
and 1509040. G. Wang and G. B. Giannakis are with the Digital Technology
Center and the Electrical and Computer Engineering Department, University
of Minnesota, Minneapolis, MN 55455, USA. G. Wang is also with the
State Key Lab of Intelligent Control and Decision of Complex Systems,
Beijing Institute of Technology, Beijing 100081, P. R. China. J. Chen is
with the State Key Lab of Intelligent Control and Decision of Complex
Systems, Beijing Institute of Technology, Beijing 100081, P. R. China. Emails: {gangwang, georgios}@umn.edu; chenjie@bit.edu.cn.

To enable grid-wide monitoring, power system engineers in
the 1960s pursued voltages at critical buses based on readings
collected from current and potential transformers. But the
power flow equations were never feasible due to timing and
modeling inaccuracies. In a seminal work [3], the statistical
foundation was laid for power system state estimation (PSSE),
whose central task is to obtain the voltage magnitude and
angle information at all buses given network parameters and
measurements acquired from across the grid. Since then, there
have been numerous PSSE contributions; see [2] for a recent
review of PSSE advances, some of which are outlined below.
Power grids are primarily monitored by supervisory control
and data acquisition (SCADA) systems. Parameter uncertainty,
instrument mis-calibration, and unmonitored topology changes
can however, yield grossly corrupted SCADA measurements
(a.k.a. ‘bad data’) [4]. Designed for functionality and efficiency with little attention paid to security, today’s SCADA
systems are vulnerable to cyberattacks [5]. Bad data also
come in the form of purposeful manipulation of smart meter
readings, as asserted by the first hacker-caused power outage:
the 2015 Ukraine blackout [6]. Any of these events can happen
which will cause a given data collection to be much more
inaccurate than is assumed by popular mathematical models.
Efficient robust PSSE approaches against cyber threats are thus
well motivated in the smart grid context [7].
Commonly used PSSE criteria include the weighted leastsquares (WLS) and the least-absolute value (LAV) [8]. The
former coincides with the maximum likelihood criterion when
additive white Gaussian noise is assumed. Unfortunately, WLS
comes with a number of limitations. First, obtaining the
WLS estimate based on SCADA measurements amounts to
minimizing a nonconvex quartic polynomial. As a result, the
quality of the resultant iterative estimators (e.g., the GaussNewton one) relies heavily on the initialization. Further,
the convergence of Gauss-Newton iterations for nonconvex
objectives is hardly guaranteed in general [9]. As least as
important, WLS estimators are sensitive to bad data [4]. They
may yield very poor estimates in the presence of outliers.
These issues were somewhat mitigated by incorporating the
largest normalized residual (LNR) test for bad data removal
[4], or, via reformulating the (possibly regularized) WLS into
a semidefinite program (SDP) via convex relaxation [10],
[11], [12], [13], [14]. The former alternates between the LNR
test and the estimation, while the latter solves SDPs. The
least-median-squares and the least-trimmed-squares estimators
have provably improved performance under certain conditions
[15]. Unfortunately, their computational complexities scale
unfavorably with the number of buses in the network [2].

2

On the other hand, LAV estimators simultaneously identify
and reject bad data while acquiring an accurate estimate of
the state [16]. Recent research efforts have focused on dealing
with the nonconvexity and nonsmoothness in LAV estimation.
Upon linearizing the nonlinear measurement functions at the
most recent iterate, a series of linear programs was solved [16].
Techniques for improving the linear programming by exploiting the system’s structure [17], or via iterative reweighting [18]
have also been reported. LAV estimation based only on PMU
data was studied [19], [20], in which a strategic scaling was
suggested to eliminate the effect of leverage measurements
[15]. Despite these efforts, LAV estimators have not been
widely employed yet in today’s power networks due mostly
to their computational inefficiency [19].
The LAV-based PSSE is revisited in this work from the
viewpoint of composite optimization [21], [22], which considers minimizing functions f (v) = h(c(v)) that are compositions of a convex function h, and a smooth vector function c.
Two novel proximal linear (prox-linear) procedures are developed based upon minimizing a sequence of convex quadratic
subproblems. The first deterministic LAV solver minimizes
functions constructed from a linearized approximation to the
original objective and a quadratic regularization, each efficiently implementable using either off-the-shelf solvers, or,
the alternating direction method of multipliers (ADMM). The
convergence of such deterministic prox-linear algorithms has
been documented [22], [23].
The second LAV solver builds on a stochastic prox-linear
algorithm, and has each iteration minimizing the summation
of a linearized approximation to the LAV loss of a single
measurement and the regularization term [24]. Interestingly,
each iteration of the stochastic LAV solver has a closedform update. Thanks to the sparse connectivity inherent to
power networks, this amounts to updating very few entries
of the state vector. In particular, if the available measurements include voltage magnitudes and power flows only, each
stochastic iteration performs about 10 scalar operations, which
is independent of the network size. If power injections are
included too, the number of operations increases to the order
of the number of neighboring buses. This complexity evidently
scales favorably to large-size power grids. Moreover, even
faster implementation of the stochastic solver is realized by
means of judiciously mini-batching the measurements. Finally,
the novel algorithms are numerically tested using the IEEE 14, 118-bus, and the PEGASE 9, 241-bus benchmark networks.
Simulations corroborate their merits relative to the WLS-based
Gauss-Newton iterations.
Outline. Grid modeling and problem formulation are given
in Section II. Upon reviewing the basics of composite optimization, Section III presents the deterministic LAV solver,
followed by its stochastic alternative in Section IV. Numerical
tests are presented in Section V, while the paper is concluded
in Section VI.
Notation. Matrices (column vectors) are denoted by upper(lower-) case boldface letters. Symbols (·)T and (·)H represent
(Hermitian) transpose, and (·) complex conjugate. Sets are
denoted using calligraphic letters. Symbol <(·) (=(·)) takes the
real (imaginary) part of a complex number. Operator dg(xi )

defines a diagonal matrix holding entries of xi on its diagonal,
while [xi ]1≤i≤N returns a matrix with xH
i being its i-th row.
II. G RID M ODELING AND P ROBLEM F ORMULATION
An electric grid having N buses and L lines is modeled
as a graph G = (N , E), whose nodes N := {1, 2, . . . , N }
correspond to buses and whose edges E := {(n, n0 )} ⊆ N ×E
to lines. The complex voltage per bus n ∈ N is expressed in
rectangular coordinates as vn = <(vn ) + j=(vn ), with all
nodal voltages forming the vector v := [v1 · · · vN ]H ∈ CN .
The voltage magnitude square Vn := |vn |2 = <2 (vn ) +
2
= (vn ) can be compactly expressed as a quadratic function of
v, namely
Vn = vH HVn v,

with HVn := en eTn

(1)

where en denotes the n-th canonical vector in RN . To express
power injections as functions of v, introduce the so-termed
bus admittance matrix Y = G + jB ∈ CN [2]. In rectangular
coordinates, the active and reactive power injections pn and
qn at bus n are given by
pn = <(vn )

N
X

[<(vn0 )Gnn0 − =(vn )Bnn0 ]

n0 =1

+ =(vn )

N
X

[=(vn0 )Gnn0 + <(vn0 )Bnn0 ]

(2)

n0 =1

qn = =(vn )

N
X

[<(vn0 )Gnn0 − =(vn )Bnn0 ]

n0 =1

− <(vn )

N
X

[=(vn0 )Gnn0 + <(vn0 )Bnn0 ]

(3)

n0 =1

which admits a compact representation as
YH en eTn + en eTn Y
(4a)
2
T
H
T
Y en en − en en Y
qn = vH Hqn v, where Hqn :=
. (4b)
2j

pn = vH Hpn v, where Hpn :=

Recognize that the line current from bus n to n0 obeys
Inn0 = eTnn0 if = eTnn0 Yf v, where if ∈ C|E| collects all line
currents, and Yf ∈ C|E|×N relates the bus voltages to all line
currents at the ‘from’ (sending) end. Ohm’s and Kirchhoff’s
laws assert that the ‘from-end’ power flow over line (n, n0 )
f
f
f
∗ f
can be expressed as S nn0 = Pnn
0 − jQnn0 = vn inn0 =
H
T
H
T
(v en )(enn0 if ) = v en enn0 Yf v, yielding
f
H P
Pnn
Hnn0 v, where HP
0 = v
nn0 :=

Q
Qfnn0 = vH HQ
nn0 v, where Hnn0 :=

YfH enn0 en + en eTnn0 Yf
2
(5a)
YfH enn0 en − en eTnn0 Yf
.
2j
(5b)

The active and reactive power flows measured at the ‘to’
t
t
(receiving) ends Pnn
0 and Qnn0 can be written symmetrically
f
f
to Pnn0 and Qnn0 ; and hence, they are omitted here for brevity.
Given line parameters collected in Y and Yf , all SCADA
measurements including squared voltage magnitudes as well

3

as active and reactive power injections and flows can be
expressed as quadratic functions of the voltages v ∈ CN .
This justifies why v is referred to as the system state. If
f
t
SV , Sp , Sq , SPf , SQ
, SPt , and SQ
signify the smart meter locations of the corresponding type, we have available
the following (possibly noisy or even corrupted) measuref
ments: {V̌n }n∈SV , {p̌n }n∈Sp , {q̌n }n∈Sq , {P̌nn
0 }(n,n0 )∈S f ,
P

t
{Q̌fnn0 }(n,n0 )∈S f , {P̌nn
and {Q̌tnn0 }(n,n0 )∈SQt ,
0 }(n,n0 )∈S t ,
P
Q

henceforth concatenated in the vector z ∈ RM , where M
denotes the total number of measurements.
In this paper, the following corruption model is considered
[25]: If {ξi } ⊆ R models an arbitrary attack (or outlier)
sequence, given the measurement matrices {Hm }M
m=1 in (1)(5a), we observe for 1 ≤ m ≤ M the samples
 H
v Hm v if m ∈ I nom
zm ≈
(6)
ξm
if m ∈ I out
where additive measurement noise can be included if ≈
is replaced by equality, and I nom , I out ⊆ {1, 2, . . . , M }
collect the indices of nominal data and outliers, respectively.
In other words, I out is the set of meter indices that can
be compromised. The indices in I out are assumed chosen
randomly from {1, 2, . . . , M }. Instrument failures occur at
random, although the attack sequence {ξm } may rely on
{Hm } (even adversarially). Specifically, two models will be
considered for the attacks.
M
M1 Matrices {Hm }M
m=1 are independent of {ξm }m=1 .
M2 Nominal measurement matrices {Hm }m∈I nom are independent of {ξm }m∈I out .
It is worth highlighting that M1 requires full independence
between the corruption and measurements. That is, the attacker
may only corrupt ξm without knowing Hm . On the contrary,
M2 allows completely arbitrary dependence between ξm and
Hm for m ∈ I out , which is natural as the type of corruption
may also rely on the individual measurement Hm being taken.
Having elaborated on the system and corruption models,
the PSSE problem can be stated as follows: Given matrices
Y, Yf , and the available measurements z ∈ RM , with
entries as in (6) obeying M1 or M2, recover the voltage
vector v ∈ CN . The first attempt may be seeking the WLS
estimate, or the ML one when assuming independent Gaussian
noise [3]. It is known however that the WLS criterion is
sensitive to outliers, and may yield very bad estimates even if
there are few grossly corrupted measurements [4]. As is well
documented in statistics and optimization, the `1 -based losses
yield median-based estimators [26], and handle gross errors
in the measurements z in a relatively benign way. Prompted
by this, we will consider here minimizing the `1 loss of the
residuals, which leads to the so-called LAV estimate [16]
M
1 X H
v Hm v − zm .
minimize f (v) :=
M m=1
v∈CN

(7)

Because of {vH Hm v}M
m=1 and the absolute-value operation,
the LAV objective in (7) is nonsmooth, nonconvex, and not
even locally convex near the optima ±v∗ . This is clear from
the scalar case f (v) = |v ∗ v−1|, where v ∈ R. A local analysis
based on convexity and smoothness is thus impossible, and

f (v) is in general difficult to minimize. For this reason, GaussNewton is not applicable to minimize (7). Nevertheless, the
criterion f (v) possesses several unique structural properties,
which we exploit next to develop efficient algorithms.
III. D ETERMINISTIC P ROX -L INEAR LAV S OLVER
In this section, we will develop a deterministic solver of
(7). To that end, let us start rewriting the objective in (7) as
minimize f (v) := h(c(v))
v∈CN

(8)

the composition of a convex function h : RM → R, and
a smooth vector function c : CN → RM , a structure that
is known to be amenable to efficient algorithms [21], [22].
It is clear that this general form subsumes (7) as a special
case,
 H for which we can take h(u) = (1/M )kuk1 and c(v) =
v Hm v − zm 1≤m≤M . The compositional structure lends
itself well to the proximal linear (prox-linear) algorithm, which
is a variant of the Gauss-Newton iterations [21]. Specifically,
define close to a given v the local “linearization” of f as

fv (w) := h c(v) + ∇H c(v)(w − v)
(9)
where ∇c(v) ∈ CN ×M denotes the Jacobian matrix of c at
v based on Wirtinger derivatives for functions of complexvalued variables [27, Appendix]. In contrast to the nonconvex
f (v), function fv (w) in (9) is convex in w, which is the key
behind the prox-linear method. Starting with some point v0 ,
which can be the flat-voltage profile point, namely the all-one
vector, construct the iteration


1
kv − vt k22
(10)
vt+1 := arg min fvt (v) +
2µt
v∈CN
where µt > 0 is a stepsize that can be fixed in advance, or be
determined by a line search [22].
Evidently, the subproblem (10) to be solved at every iteration of the prox-linear algorithm is convex, and can be
handled by off-the-shelf solvers such as CVX [28]. However,
these interior-point based solvers may not scale well when
the matrices {Hm } are large. For this reason, we derive next
a more efficient iterative procedure using ADMM iterations
[29].
When specifying f to be the LAV objective of (7), the
minimization in (10) becomes
vt+1 = arg min kAt (v − vt ) − ct k1 +
v∈CN

1
2
kv − vt k2 (11)
2

with coefficients given by


At := (2µt /M )vtH Hm 1≤m≤M


ct := (µt /M )(zm − vtH Hm vt ) 1≤m≤M .

(12a)
(12b)

For brevity, let w := v − vt , and rewrite (11) equivalently as
a constrained optimization problem
1
ku − ct k1 + kwk22
2
subject to At w = u.
minimize

u∈CM ,

w∈CN

(13a)
(13b)

4

To decouple constraints and also facilitate the implementation
of ADMM, introduce an auxiliary copy ũ and w̃ for u and w
accordingly, and rewrite (13) into
1
2
kw̃k2
2
subject to ũ = u, w̃ = w, At w = u.
minimize kũ − ct k1 +

(14a)

ũ, w̃, u, w

(14b)

Letting λ ∈ CN and ν ∈ CM be the Lagrange multipliers
corresponding to the w- and u-consensus constraints, respectively, the augmented Lagrangian after leaving out the last
equality in (14b) can be expressed as
1
2
L(w̃, ũ,w, u; λ, ν) := kũ − ct k1 + kw̃k2 + λH (w̃ − w)
2
ρ
ρ
2
2
+ ν H (ũ − u) + kw̃ − wk2 + kũ − uk2 (15)
2
2
where ρ > 0 is a predefined learning rate. With k ∈ N
denoting the iteration index for solving (13), or equivalently
(10), ADMM cycles through the following recursions


1
ρ
k
k 2
k+1
2
w̃ − (w − λ ) 2 (16a)
w̃
:= arg min kw̃k2 +
2
2
w̃


1
ρ
k
k 2
k+1
ũ −(u −ν ) 2 (16b)
ũ
:= arg min kũ−ct k1 +
2
2
ũ
 k+1 k+1
w
,u
:=
arg min

w − (w̃

k+1

w, u



k

+λ )

2
2

k+1

+ u − (ũ

k

+ν )

s. to At w = u
  k

λ
λ + (w̃k+1 − wk+1 )
=
ν k+1
ν k + (ũk+1 − uk+1 )

2
2

(16c)

k+1

(16d)

where all the dual variables have been scaled by the factor
ρ > 0 [29]. Interestingly enough, the solutions of (16a)-(16c)
can be provided in closed form, as we elaborate next.
Proposition 1. The solutions of problems (16a) and (16b) are

ρ
w k − λk
1+ρ

:= ct + S1/ρ uk − ν k − ct

w̃k+1 :=

(17a)

ũk+1

(17b)

where the shrinkage operator Sτ (x) : CN × R+ → RN is
Sτ (x) := sign(x) max(|x| − τ 1, 0) [30], with
and | · |
denoting the entrywise multiplication and modulus operators,
nx
/|x|, x 6= 0
respectively, and sign(x) :=
provides an
0,
x=0
entrywise definition of sign(x).
The constrained minimization of (16c) essentially projects
the pair (w̃k+1 +λk , ũk+1 +ν k ) onto the convex set specified
by the linear equality constraint, namely {(w, u) : At w = u}.
Its solution is derived in a simple closed form next (see the
Appendix for a proof).
Proposition 2. Given b ∈ CN and d ∈ CM , the solution of
minimize

w∈CN , u∈CM

subject to

1
1
2
2
kw − bk2 + ku − dk2
2
2
Aw = u

TABLE I: Deterministic LAV Solver Using ADMM
1: Input data {(zm , Hm )}M
m=1 , stepsizes µ, ρ > 0, and initialization v0 = 1.
2: For t = 0, 1, . . . , do
3: Initialize w0 , u0 , λ0 , and ν 0 to zero.
4: Evaluate At and ct as in (12).
5: For k = 0, 1, . . . , do
6:
Update (w̃k+1 , ũk+1 ), (wk+1 , uk+1 ), and (µk+1 , ν k+1 )
using (17), (19), and (16d), respectively.
7: End for
8: Update vt+1 via (20).
9: End for

is given as
w∗ := I + AH A

−1


b + AH d

u∗ := Aw∗ .

(18a)
(18b)

Using Proposition 2, the minimizer of (16c) is found as
 k+1

−1
wk+1 := (I +AH
w̃
+ λk + AH (ũk+1 + ν k )
t At )
(19a)
uk+1 := Awk+1 .

(19b)

The four updates in (16) are computationally simple except
for the matrix inversion of (19a), which nevertheless can be
cached once computed during the first iteration. In addition,
variables ũ0 , λ0 , and ν 0 of ADMM can be initialized to zero.
Finally, the solution of (11) can be obtained as
vt+1 := vt + w∗

(20)

where w∗ is the converged w-iterate of the ADMM iterations
in (17), (19), and (16d).
The novel deterministic LAV solver based on ADMM is
summarized in Table I, in which the inner loop consisting of
Steps 3-8 can be replaced with off-the-shelf solvers to solve
(11) for vt+1 .
As far as performance is concerned, if h is L-Lipschitz and
∇c is β-Lipschitz, then taking a constant stepsize µ ≤ 1/(Lβ)
in (10) guarantees that [23, Sec. 5]:
i) the proposed solver in Table I is a descent method; and,
ii) the iterate sequence {vt } converges to a stationary point
of the LAV objective in (7).
The computational burden of the ADMM based deterministic solver is dominated by the projection step of (19), which
incurs per-iteration computational complexity on the order
of O(M N 2 ). This complexity can be afforded in small- or
medium-size PSSE tasks, but may not be efficient enough for
nowadays increasingly interconnected power networks. This
motivates our stochastic alternative of the ensuing section that
relies on very inexpensive iterations.
IV. S TOCHASTIC P ROX - LINEAR LAV S OLVER
Finding the minimizer of (11) exactly per iteration of the
deterministic LAV solver may be computationally expensive,
and can be intractable when the network size grows very large.
Considering the wide applicability of LAV estimation as well

5

as the increasing interconnection of microgrids, scalable online
and stochastic approaches become of substantial interest. In
this section, a stochastic linear proximal algorithm of [24] is
adapted to our PSSE task, and enables the prox-linear method
to efficiently solve the LAV estimation problem at scale. Advantages of the stochastic approaches over their deterministic
counterparts include oftentimes simple closed-form updates as
well as faster convergence to yield an (approximately) optimal
solution. Leveraging the sparsity structure of the measurement
matrices and judiciously grouping measurements into small
mini-batches can considerably speedup implementation of the
stochastic solver.
A. Stochastic LAV solver
Instead of dealing with the quadratic subproblems in (11),
each iteration of the stochastic LAV solver samples a datum
mt ∈ {1, 2, . . . , M } uniformly at random from the total M of
measurements, and substitutes functions (h, c) by (hmt , cmt )
associated with the sampled datum in the local linearization
of (9), hence also in (11), yielding


1
2
H
kv − vt k2
vt+1 := arg min amt (v − vt ) − cmt +
2µt
v∈CN
(21)
where the coefficients are given by
amt := 2Hmt vt
cmt := zmt −

(22a)

vtH Hmt vt .

(22b)

Different from iteratively seeking solutions of (11) based on
ADMM iterations, the minimization of (21) admits a simple
closed-form minimizer presented in the next result, which is
proved in the Appendix.
Proposition 3. Given a ∈ CN and c ∈ R, the solution of

TABLE II: Stochastic LAV Solver
1: Input data {(zm , Hm )}M
m=1 , stepsize {µt > 0}, and initialization v0 = 1.
2: For t = 0, 1, . . . do
3: Draw mt uniformly at random from {1, 2, . . . , M }, or, by
cycling through {1, 2, . . . , M }.
4: Evaluate amt and cmt as in (22).
5: Update vt+1 using (24).
6: End for

For instance, one can choose µt = αt−β with appropriately
selected constants α > 0 and β ∈ (0.5, 1]. Then the sequence
{vt } converges to a stationary point of the LAV objective in
(7) almost surely [24, Thm. 1].
In terms of computational complexity, it can be verified that
each Hm matrix corresponding to a square voltage magnitude
or (re)active power flow measurement [cf. (1) and (5)] has
exactly one or three nonzero entries, respectively. As such, if
the available measurements include only these two types of
measurements, evaluating (amt , cmt ) as well as updating vt
per stochastic iteration requires just a small number (≤ 10) of
scalar multiplications and additions, therefore incurring periteration complexity of O(1), which is independent of the
network size N . This complexity evidently scales favorably
to very large interconnected power networks. It is also worth
highlighting that only one or two entries of vt are updated
depending on whether a voltage or power flow measurement
is processed at each iteration. On the other hand, even if
the power injections are measured too, the number of scalar
operations per iteration increases to the order of the number
of neighboring buses, which still remains much smaller than
N in most real-world networks.

(23)

Remark 1. PMU measurements can be easily accounted for
in (7). The developed prox-linear LAV schemes apply without
any algorithmic modification.

is given by ŵ := projτ (c/kak22 ) · a, where the projection
operator projτ (x) : R × R+ → R returns the real number in
interval [−τ, τ ] closest to any given x ∈ R.

Remark 2. An auxiliary sparse variable o ∈ RM has been
introduced to model explicitly the outliers in [31]. These robust
schemes with WLS estimation regularized by kok1 can also
benefit from our efficient and scalable algorithms.

minimize
w∈CN

aH w − c +

1
kwk22
2τ

Based on Proposition 3, the solution of (21) is given by
vt+1 := vt + projµt(cmt /kamt k22 ) · amt .

(24)

Intuitively, measurements with a relatively small absolute
residual, namely |cmt | ≤ kamt k22 , are deemed ‘nominal,’ and
vt is updated with a step of cmt /kamt k22 along the current
direction of amt . On the other hand, the measurements of
larger absolute residuals are likely to be outliers, so vt is
updated along its direction amt by only a step of τt as opposed
to cmt /kamt k22 .
The proposed stochastic prox-linear LAV solver is listed
in Table II. For convergence, a diminishing stepsize sequence
{µt } is required. Specifically, we consider stepsizes that are
square summable but not summable; that is,
∞
X
t=0

µt = ∞,

and

∞
X
t=0

µ2t < ∞.

(25)

B. Accelerated implementation using mini-batches
Although it involves only simple closed-form updates, the
stochastic solver in Table II may require a large number of
iterations to converge for high-dimensional power networks.
Stochastic approaches based on mini-batches of measurements
have been recently popular in large-scale machine learning
tasks, because they offer a means of accelerating the stochastic
algorithms. However, the naive way of designing mini-batches
by grouping measurements randomly would yield a sequence
of quadratic programs as in (11) of the deterministic LAV
solver, which does not have closed-form solutions due to the
`1 term. The novelty here is to fully exploit the sparsity of Hm
matrices to group measurements into mini-batches in a way
that closed-form solutions of the resulting quadratic programs
are possible.

6

12

13

6

11

1

5

2

14

10

9

4

7

3

8

Groups

Power flow measurements (Bus index)

1

(1, 2), (3, 4), (5, 6), (7, 9)

2

(1, 5), (2, 4), (6, 11), (12, 13)

3

(2, 3), (4, 5), (6, 12), (9, 14)

4

(2, 5), (4, 7), (6, 13), (10, 11)

5

(4, 9), (7, 8), (9, 10), (13, 14)

Fig. 1: The IEEE 14-bus benchmark configuration.

TABLE III: Mini-batching power flow measurements.

Suppose that active and reactive power flows over all lines
and the square voltage magnitude at all buses are measured.
Since HVn in (1) has exactly one nonzero entry at the (n, n)th position, the corresponding updating vector an in (22a)
has (at most) one nonzero entry at the n-th position. Updating vt using (24) modifies the n-th position only. Hence,
refining v using the N voltage measurements sequentially
in N stochastic iterations is equivalent to updating v using
all N measurements simultaneously in a single iteration. For
power flows, every Hm has three nonzero entries in two rows.
0
For instance, HP
nn0 has nonzero entries indexed by (n, n ),
Q
0
(n, n), and (n , n); and so does Hnn0 . Processing the active
or reactive power flow measurement over line (n, n0 ) amounts
to updating the n-th and n0 -th entries of vt . Hence, so long
as each of a mini-batch of measurements does not share
common indices with the remaining ones, processing a minibatch of such measurements one after the other boils down to
processing all measurements simultaneously in one iteration.
For illustration, consider the IEEE 14-bus test system depicted in Fig. 1 [32]. Consider a total of 54 measurements,
which include 14 square voltage magnitudes, as well as 20
‘from-end’ active and reactive power flows each. All voltage
magnitudes can be grouped as a single mini-batch, or into
several mini-batches by any means. One way of mini-batching
each type of power flow measurements is suggested in Table
III, where 20 active (reactive) power flows yield 5 minibatches of equal size. It can be easily verified that any two
measurements within a group (mini-batch) are measured over
two lines of non-overlapping indices.
Let the entire measurements be divided into B mini-batches
denoted by {Bb }B
b=1 ⊆ {1, 2, . . . , M }. If a mini-batch Bbt of
measurements is drawn uniformly at random from {Bb }B
b=1 at
iteration t, the accelerated implementation by means of minibatching, updates the state estimate according to [cf. (24)]
X
vt+1 := vt +
projµt (cm /kam k22 ) · am
(26)

using MATLAB R2016a. Three power network benchmarks
including the IEEE 14-, 118-bus, and the PEGASE 9, 241bus systems were simulated, following the MATLAB-based
toolbox MATPOWER [32], [33]. The ‘workhorse’ GaussNewton iterations for the WLS-based PSSE were adopted
as a baseline, and were implemented using the SE function
‘doSE.m’ in MATPOWER. Moreover, when all squared voltage magnitudes are measured, the initial point is taken to be
the voltage magnitude vector, unless otherwise specified. In
order to fix the phase ambiguity, the phase generated at the
reference bus was set to 0 in all tests. For numerical stability,
and to eliminate the effect of leverage measurements [19],
the developed solvers were implemented using the normalized
Hm
M
m
data, namely {( kHzm
k2 , kHm k2 )}m=1 .

m∈Bbt

which is in sharp contrast to that of using ADMM iterations to
deal with the quadratic subproblems (10) in the deterministic
LAV solver.
V. N UMERICAL T ESTS
The proposed linear proximal LAV solvers were numerically
tested on an Intel CPU @ 3.4 GHz (32 GB RAM) computer

A. Noiseless case
The first experiment simulates the noiseless data to evaluate
the convergence and runtime of the novel algorithms relative
to the WLS-based Gauss-Newton iterations on the IEEE 14bus test system. The default voltage profile was employed.
Measurements including all active and reactive power flows,
as well as all squared voltage magnitudes were obtained
from MATPOWER [33]. The deterministic solver in Table I
was implemented with 200 ADMM iterations and stepsizes
µ = ρ = 0.35. The stochastic algorithm in Table II used
the diminishing stepsize 1/t0.8 . The accelerated scheme in
(26) was implemented with stepsize 0.8 using a total of 11
mini-batches: 1 for all voltage magnitudes, and 5 of equal
size for active and reactive power flows,
PM each grouped as in
Table III. The WLS objective (1/M ) m=1 (zm − vH Hm v)2
was evaluated at every Gauss-Newton iteration, while the
LAV objective in (7) was computed per iteration of the
deterministic, and every M iterations of the stochastic and
accelerated schemes. Figure 2 plots the objective values for
the Gauss-Newton, deterministic, stochastic, and accelerated
schemes, whose corresponding runtime is 0.034s, 0.148s,
0.053s, and 0.042s. The Gauss-Newton stopped after a few
iterations, but at a near-optimal point of relative root meansquare error (RMSE) kv̂ − vk2 /kvk2 = 1.6 × 10−3 , where
v denotes the true voltage profile, and v̂ the estimate. The
accelerated implementation was faster than the deterministic
and stochastic LAV solvers, and yielded an accurate solution

7

10-1

101

Gauss-Newton
Deterministic
Stochastic
Accelerated

100

RMSE (p.u.)

Objective value

10-2

10-3

10-4

10-5

10-1

10-2

10-3

10-6

Gauss-Newton
Deterministic
Stochastic

10-4
1

10

20

30

40

50

60

Iterations

3

4

5

6

7

Measurement size

Fig. 2: Convergence performance for the IEEE 14-bus system.

Fig. 3: Robustness to outliers for the IEEE 118-bus system.

with RMSE = 4.88 × 10−5 in time comparable to the GaussNewton iterations.

data were measured with additive noise described in the last
experiment, 5% of which were compromised under model M2.
The corrupted data ξm := ṽH Hm ṽ relying on the individual
Hm were generated using ṽ ∈ RN from the standardized multivariate Gaussian distribution. In words, a number of 18, 481
variables were estimated, a total of 91, 919 measurements
were obtained, 4, 595 of which were purposefully manipulated
by adversaries. Initialized with the flat voltage profile point,
the WLS-based Gauss-Newton iterations yielded an estimate
of RMSE 0.9656, while the stochastic scheme with stepsize
µt = 100/t0.8 attained an RMSE of 0.0163. Moreover, the
novel stochastic LAV scheme is at least three times faster than
the Gauss-Newton iterations.

B. Presence of outliers
The second experiment was set to assess the robustness of
the novel deterministic and stochastic solvers to measurements
with outliers using the IEEE 118-bus benchmark network
[32]. The WLS-based Gauss-Newton iterations were simulated
as a baseline. The actual voltage magnitude (in p.u.) and
angle (in radians) of each bus were uniformly distributed over
[0.95, 1.05], and over [−0.05π, 0.05π]. To assess the PSSE
performance versus the measurement size, an additional type
of measurements was included in a deterministic manner, as
described next. All seven types of SCADA measurements were
f
t
first enumerated as {|Vn |2 , Pnk
, Qfnk , Pn , Qn , Pnk
, Qtnk }.
Each x-axis value in Fig. 3 signifies the number of ordered
types of measurements used in the experiment to yield the
corresponding RMSEs, which are obtained by averaging over
100 independent realizations. For example, 5 implies that
f
the first 5 types of data (i.e., |Vn |2 , Pnk
, Qfnk , Pn , Qn over
all buses and lines) were measured. Additive noise was
independently generated from Gaussian distributions having
zero-mean and standard deviation 0.004, 0.008, and 0.01 p.u.
for the voltage magnitude, line flow, and power injection
measurements, respectively [17]. Five percent of the measured
data were corrupted according to model M1, which were
chosen uniformly at random from line flows and bus injections.
The outlying data {ξm } were drawn independently from a
Gaussian distribution with zero-mean and variance 5. The
subproblems (10) with µ = 0.1 of the deterministic scheme
were solved with CVX [28], while the stochastic one was run
with a diminishing stepsize µt = 10/t0.8 . It is evident from
Fig. 3 that our novel prox-linear LAV schemes are resilient to
outlying measurements under corruption model M1.
The last experiment tests the scalability and efficacy of the
stochastic iterations on a larger power network of 9, 241 buses
available in MATPOWER [33]. The true voltage magnitude
of each bus was uniformly distributed over [0.9, 1.1], and
its angle over [−0.1π, 0.1π]. All seven types of SCADA

VI. C ONCLUSIONS
Robust power system state estimation was pursued using
contemporary tools of composite optimization. Building upon
recent algorithmic advances, two novel solvers were put forward to efficiently handle the LAV-based PSSE. Specifically,
a deterministic scheme was developed based on a linear
proximal method, yielding a sequence of convex quadratic
subproblems that can be efficiently solved using off-the-shelf
solvers, or, through ADMM iterations. Inspired by the sparse
connectivity inherent to power networks, a highly scalable
stochastic scheme that can afford simple closed-form updates
was also devised. When only line flows and voltage magnitudes are measured, each stochastic iteration performs merely
a few (less than 10) scalar operations, incurring per-iteration
computational runtime of O(1), regardless of the number of
buses in the network. If the power injections are included as
well, this complexity goes down to the order of the number of
neighboring buses, which still remains much smaller than the
network size in general. A mini-batching technique was also
suggested to further accelerate the stochastic iterations by judiciously leveraging the sparsity structure of the measurement
matrices. Numerical tests on a variety of benchmark networks
showcase the robustness and computational efficiency of the
novel approaches.
Devising decentralized and parallel implementations of the
novel approaches constitutes interesting future research di-

8

rections. Coping with the Y - and ∆-connection, as well as
extending the approaches to multiphase unbalanced distribution systems are practically relevant future research topics
too. Another possibility includes leveraging recent advances
in solving systems of quadratic equations to tackle PSSE [34].
A PPENDIX
Proof of Prop. 2: Letting χ denote the dual variable
associated with the constraint u = Aw, the KKT optimality
conditions are given by [29]
w∗ − b + AH χ∗ = 0
u∗ − d − χ∗ = 0
Aw∗ − u∗ = 0
or in the following compact form

 ∗  

w
b
IN
0
AH
 0
IM
−IM   u∗  =  d  .
χ∗
0
A −IM
0
Eliminating the dual variable via χ∗ = u∗ − d from the
KKT system, yields

 ∗  

w
IN A H
b + AH d
.
(27)
=
u∗
A −IM
0
By further eliminating d∗ and solving for b∗ , the solution
to (27) and also to the minimization in (18) can be found in
two steps as
−1

w∗ := I + AH A
b + AH d
u∗ := Aw∗
which completes the proof of the claim.
Proof of Prop. 3: For any constant c ∈ R, it is easy
to verify that aH w at the optimum w∗ takes real values by
contradiction. The optimality condition for (23) is
 1
1
0 ∈ ∂ aH w − c + w ⇐⇒ 0 ∈ ∂ aH w − c · a + w
τ
τ
or equivalently, 0 ∈ ∂ aH w − c · (τ a) + w, where ∂ denotes
the subdifferential. Let us first examine the case where aH w−
c 6= 0. We thus have ∂|aH w − c| = sign(aH w − c), which
yields the optimum
w∗ = −τ sign(aH w − c) · a.
≥
τ , or aH w∗ − c
=
Note that if c/kak22
−τ kak22 sign(aH w∗ − c) − c < 0, then w∗ = τ a.
Equivalently, if c/kak22 ≤ −τ , or aH w∗ − c =
−τ kak22 sign(aH w∗ − c) − c > 0, then w∗ = −τ a.
If aH w − c = 0, the subdifferential of the absolute operator
is the interval [−1, 1]; hence, the optimality condition becomes
0 ∈ −[−1, 1] · (τ a) + w ⇐⇒ w ∈ [−τ, τ ]a.
Upon letting projτ (x) denote the projection of a real
number x onto the interval [−τ, τ ], one can combine the aforementioned three cases, and write collectively the optimum as
w∗ := projτ (c/kak22 ) · a
which concludes the proof of the proposition.

R EFERENCES
[1] W. A. Wulf, “Great achievements and grand challenges,” The
Bridge, vol. 30, no. 3/4, pp. 5–10, Fall 2010. [Online]. Available:
http://www.greatachievements.org/.
[2] V. Kekatos, G. Wang, H. Zhu, and G. B. Giannakis, “PSSE redux:
Convex relaxation, decentralized, robust, and dynamic approaches,”
arXiv:1708.03981, 2017.
[3] F. C. Schweppe, J. Wildes, and D. Rom, “Power system static state
estimation: Parts I, II, and III,” IEEE Trans. Power App. Syst., vol. 89,
pp. 120–135, Jan. 1970.
[4] H. M. Merrill and F. C. Schweppe, “Bad data suppression in power
system static state estimation,” IEEE Trans. Power App. Syst., no. 6, pp.
2718–2725, Nov. 1971.
[5] P. Fairley, “Cybersecurity at US utilities due for an upgrade: Tech to
detect intrusions into industrial control systems will be mandatory,”
IEEE Spectr., vol. 53, no. 5, pp. 11–13, May 2016.
[6] D. U. Case, “Analysis of the cyber attack on the Ukrainian power grid,”
2016.
[7] S. Zonouz, K. M. Rogers, R. Berthier, R. B. Bobba, W. H. Sanders, and
T. J. Overbye, “SCPSE: Security-oriented cyber-physical state estimation
for power grid critical infrastructures,” IEEE Trans. Smart Grid, vol. 3,
no. 4, pp. 1790–1799, Dec. 2012.
[8] E. Caro and A. Conejo, “State estimation via mathematical programming: A comparison of different estimation algorithms,” IEE Proc.Gener. Transm. Distrib., vol. 6, no. 6, pp. 545–553, Jun. 2012.
[9] D. P. Bertsekas, Nonlinear Programming, 2nd ed.
Belmont, MA:
Athena Scientific, 1999.
[10] H. Zhu and G. B. Giannakis, “Power system nonlinear state estimation
using distributed semidefinite programming,” IEEE J. Sel. Topics Signal
Process., vol. 8, no. 6, pp. 1039–1050, Dec. 2014.
[11] S.-J. Kim, G. Wang, and G. B. Giannakis, “Online semidefinite programming for power system state estimation,” in Proc. IEEE Conf. on
Acoustics, Speech and Signal Process., Florence, Italy, May 2014, pp.
6024–6027.
[12] G. Wang, S.-J. Kim, and G. B. Giannakis, “Moving-horizon dynamic
power system state estimation using semidefinite relaxation,” in Proc.
IEEE PES General Meeting, Washington, DC, July 2014, pp. 1–5.
[13] R. Madani, M. Ashraphijuo, J. Lavaei, and R. Baldick, “Power system
state estimation with a limited number of measurements,” Preprint, 2016.
[14] Y. Zhang, R. Madani, and J. Lavaei, “Conic relaxations for power system
state estimation with line measurements,” IEEE Trans. Control Netw.
Syst., 2017 (to appear).
[15] L. Mili, M. G. Cheniae, and P. J. Rousseeuw, “Robust state estimation
of electric power systems,” IEEE Trans. Circuits Syst. I, vol. 41, no. 5,
pp. 349–358, May 1994.
[16] W. W. Kotiuga and M. Vidyasagar, “Bad data rejection properties of
weighted least-absolute-value techniques applied to static state estimation,” IEEE Trans. Power App. Syst., no. 4, pp. 844–853, Apr. 1982.
[17] A. Abur and M. K. Celik, “A fast algorithm for the weighted leastabsolute-value state estimation (for power systems),” IEEE Trans. Power
Syst., vol. 6, no. 1, pp. 1–8, Feb. 1991.
[18] R. Jabr and B. Pal, “Iteratively re-weighted least-absolute-value method
for state estimation,” IEE Proc.-Gener. Transm. Distrib., vol. 150, no. 4,
pp. 385–391, July 2003.
[19] M. Gol and A. Abur, “LAV based robust state estimation for systems
measured by PMUs,” IEEE Trans. Smart Grid, vol. 5, no. 4, pp. 1808–
1814, Jul. 2014.
[20] C. Xu and A. Abur, “A fast and robust linear state estimator for very large
scale interconnected power grids,” IEEE Trans. Smart Grid, vol. PP,
no. 99, pp. 1–1, 2017.
[21] J. V. Burke, “Descent methods for composite nondifferentiable optimization problems,” Math. Program., vol. 33, no. 3, pp. 260–279, Dec. 1985.
[22] A. S. Lewis and S. J. Wright, “A proximal method for composite
minimization,” Math. Program., vol. 158, no. 1-2, pp. 501–546, July
2016.
[23] D. Drusvyatskiy and A. S. Lewis, “Error bounds, quadratic growth, and
linear convergence of proximal methods,” arXiv:1602.06661, 2016.
[24] J. C. Duchi and F. Ruan, “Stochastic methods for composite optimization
problems,” arXiv:1703.08570, 2017.
[25] ——, “Solving (most) of a set of quadratic equalities: Composite
optimization for robust phase retrieval,” arXiv:1705.02356, 2017.
[26] P. J. Huber, “Robust Statistics,” in International Encyclopedia of Statistical Science. Springer, 2011, pp. 1248–1251.

9

[27] G. Wang, A. S. Zamzam, G. B. Giannakis, and N. D. Sidiropoulos,
“Power system state estimation via feasible point pursuit: Algorithms
and Crmér-Rao bound,” arXiv:1705.04031, 2017.
[28] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex
programming,” 2008. [Online]. Available: http://cvxr.com/cvx.
[29] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed
optimization and statistical learning via the alternating direction method
of multipliers,” Found. Trends Mach. Learn., vol. 3, pp. 1–122, 2010.
[30] A. Maleki, L. Anitori, Z. Yang, and R. G. Baraniuk, “Asymptotic
analysis of complex LASSO via complex approximate message passing
(CAMP),” IEEE Trans. Inf. Theory, vol. 59, no. 7, pp. 4290–4308, July
2013.

[31] V. Kekatos and G. B. Giannakis, “Distributed robust power system state
estimation,” IEEE Trans. Power Syst., vol. 28, no. 2, pp. 1617–1626,
May 2013.
[32] Power systems test case archive. Univ. of Washington. [Online].
Available: http://www.ee.washington.edu/research/pstca.
[33] R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas, “MATPOWER: Steady-state operations, planning and analysis tools for power
systems research and education,” IEEE Trans. Power Syst., vol. 26, no. 1,
pp. 12–19, Feb. 2011.
[34] G. Wang, G. B. Giannakis, and Y. C. Eldar, “Solving systems of random
quadratic equations via truncated amplitude flow,” IEEE Trans. Inf.
Theory, 2017 (to appear); see also arXiv:1605.08285, 2016.

