arXiv:1708.03981v1 [cs.SY] 14 Aug 2017

PSSE Redux: Convex Relaxation, Decentralized,
Robust, and Dynamic Approaches
Vassilis Kekatos, Gang Wang, Hao Zhu, and Georgios B. Giannakis

∗

August 15, 2017

1

Introduction

With the advent of digital computers, power system engineers in the 1960s
tried computing the voltages at critical buses based on readings from current and potential transformers. Local personnel manually collected these
readings and forwarded them by phone to a control center. Nevertheless,
due to timing, modeling, and instrumentation inaccuracies, the power flow
equations were always infeasible. In a seminal contribution [61], the statistical foundations were laid for a multitude of grid monitoring tasks, including
topology detection, static state estimation, exact and linearized models, bad
data analysis, centralized and decentralized implementations, as well as dynamic state tracking. Since then, different chapters, books, and review articles have nicely outlined the progress in the area; see for example [74, 55, 1].
The revolutionary monitoring capabilities enabled by synchrophasor units
have been put forth in [58].
This chapter aspires to glean some of the recent advances in power system state estimation (PSSE), though our collection is not exhaustive by
any means. The Cramér-Rao bound, a lower bound on the (co)variance
∗

V. Kekatos is with the Bradley Department of Electrical and Computer Engineering,
Virginia Tech, Blacksburg, VA 24061, USA. G. Wang is with the Department of Electrical and Computer Engineering, University of Minnesota, Minneapolis, MN 55455, USA,
and also with the State Key Lab of Intelligent Control and Decision of Complex Systems, Beijing Institute of Technology, Beijing 100081, P. R. China. H. Zhu is with the
Department of Electrical and Computer Engineering, University of Illinois, Urbana, IL
61801. G. B. Giannakis is with the Department of Electrical and Computer Engineering, University of Minnesota, Minneapolis, MN 55455, USA. E-mails: kekatos@vt.edu;
{gangwang, georgios}@umn.edu; haozhu@illinois.edu.

1

of any unbiased estimator, is first derived for the PSSE setup. After reviewing the classical Gauss-Newton iterations, contemporary PSSE solvers
leveraging relaxations to convex programs and successive convex approximations are explored. A disciplined paradigm for distributed and decentralized
schemes is subsequently exemplified under linear(ized) and exact grid models. Novel bad data processing models and fresh perspectives linking critical
measurements to cyber-attacks on the state estimator are presented. Finally, spurred by advances in online convex optimization, model-free and
model-based state trackers are reviewed.
Notation: Lower- (upper-) case boldface letters denote column vectors (matrices), and calligraphic letters stand for sets. Vectors 0, 1, and
en denote respectively the all-zero, all-one, and the n-th canonical vectors
of suitable dimensions. The conjugate of a complex-valued object (scalar,
vector or matrix) x is √
denoted by x∗ ; <{x} and ={x} are its real and imaginary parts, and j := −1. Superscripts T and H stand for transpose and
conjugate-transpose, respectively, while Tr(X) is the trace of matrix X. A
diagonal matrix having vector x on its main diagonal is denoted by dg(x);
whereas, the vector of diagonal entries of X is dg(X). The range space of X
is denoted by range(X); and its null space (kernel) by null(X). The notation
N (µ, Σ) represents the Gaussian distribution with mean µ and covariance
matrix Σ.

2

Power Grid Modeling

This section introduces notation and briefly reviews the power flow equations; for detailed exposition see e.g., [1], [19] and references therein. A
power system can be represented by the graph G = (B, L), where the node
set B comprises its Nb buses, and the edge set L its Nl transmission lines.
Given the focus on alternating current (AC) power systems, steady-state
voltages and currents are represented by their single-phase equivalent phasors per unit.
A transmission line (n, k) ∈ L running across buses n, k ∈ B is modeled
by its total series admittance ynk = gnk + jbnk , and total shunt susceptance
jbsnk . If Vn is the complex voltage at bus n, the current Ink flowing from
bus m to bus n over line (m, n) is
Ink = (ynk + jbsnk /2)Vn − ynk Vn .

(1)

The current Inm coming from the other end of the line can be expressed
symmetrically. That is not the case if the two buses are connected via a
2

transformer with complex ratio ρnk followed by a line, where
ynk + jbsnk /2
ynk
Vn − ∗ Vn
2
|ρnk |
ρnk
ynk
= (ynk + jbsnk /2)Vn −
Vn .
ρnk

Ink =

(2a)

Inm

(2b)

Kirchhoff’s
current law dictates that the current injected into bus n is
P
In = n∈Bn Ink , where Bn denotes the set of buses directly connected to
bus n. If vector i ∈ CNb collects all nodal currents, and v ∈ CNb all nodal
voltages, the two vectors are linearly related through the bus admittance
matrix Y = G + jB as
i = Yv.
(3)
Similar to (3), line currents can be stacked in the 2Nl -dimensional vector if ,
and expressed as a linear function of nodal voltages
if = Yf v

(4)

for some properly defined 2Nl × Nb complex matrix Yf [cf. (1)–(2)].
The complex power injected into bus n will be denoted by Sn := Pn +
jQn . Since by definition Sn = Vn In∗ , the vector of complex power injections
s = p + jq can be expressed as
s = dg(v)i∗ = dg(v)Y∗ v∗ .

(5)

∗ .
The power flowing from bus n to bus n over line (m, n) is Snk = Vn Ink
jθ
If voltages are expressed in polar form Vn = Vn e n , the power flow
equations in (5) per real and imaginary entry can be written as
Nb
X

Pn =

n=1
Nb
X

Qn =

Vn Vn [Gnk cos(θn − θn ) + Bnk sin(θn − θn )]

(6a)

Vn Vn [Gnk sin(θn − θn ) − Bnk cos(θn − θn )] .

(6b)

n=1

Since power injections are invariant if voltages are shifted by a common
angle, the voltage phase is arbitrarily set to zero at a particular bus called
the reference bus.
Alternatively to (6), if voltages are expressed in rectangular coordinates
Vn = Vr,m + jVi,m , power injections are quadratically related to voltages
Pn = Vr,m

Nb
X

(Vr,n Gnk − Vi,n Bnk ) + Vi,m

n=1

Nb
X

(Vi,n Gnk + Vr,n Bnk )

n=1

3

(7a)

Qn = Vi,m

Nb
X

(Vr,n Gnk − Vi,n Bnk ) − Vr,m

Nb
X

(Vi,n Gnk + Vr,n Bnk ).

(7b)

n=1

n=1

To compactly express (7), observe that Sn∗ = Vn∗ In = (vH en )(eTn i) =
vH en eTn Yv from which it readily follows that
Pn = vH HPn v

(8a)

Qn = vH HQn v

(8b)

where the involved matrices are defined as

1
en eTn Y + YH en eTn
2

1
en eTn Y − YH en eTn .
:=
2j

HPn :=

(9a)

HQn

(9b)

Similar expressions hold for the squared voltage magnitude at bus n:
Vn2 = vH HVn v, where HVn := en eTn .

(10)

Realizing that a line current can also be provided as Ink = eTnk if , the
∗ = V ∗I
power flow on line (n, k) as seen from bus n is expressed as Snk
n nk =
H
T
H
T
(v en )(enk if ) = v en enk Yf v, from which it follows that
Pnk = vH HPnk v
H

Qnk = v HQnk v

(11a)
(11b)

where HPnk and HQnk are defined by substituting eTn Y and YH en by eTnk Yf
and YfH enk in (9), accordingly.
Equations (8), (10) and (11) explain how power injections, flows, and
squared voltage magnitudes are quadratic functions of voltage phasors as
described by vH Hm v for certain complex Nb × Nb matrices Hm . Regardless
if Y and/or Yf are symmetric or Hermitian, Hm are Hermitian by definition.
T = <{H }, and
This means that Hm = HH
m
m , or equivalently, <{Hm }
={Hm }T = −={Hm }. It can be easily verified that the quadratic functions
can be expressed in terms of real-valued quantities as
vH Hm v = v̄T H̄m v̄

(12)

for the expanded real-valued voltage vector v̄ := [<{v}T ={v}T ]T , and the
real-valued counterpart of Hm , namely


<{Hm } −={Hm }
H̄m :=
.
(13)
={Hm } <{Hm }
4

3

Problem Statement

It was seen in Section 2 that given grid parameters collected in Y and Yf , all
power system quantities can be expressed in terms of the voltage vector v,
which justifies its term as the system state. Meters installed across the grid
measure electric quantities, and forward their readings via remote terminal
units to a control center for grid monitoring. Due to lack of synchronization,
conventional meters cannot utilize the angle information of phasorial quantities. For this reason, legacy measurements involve phaseless power injections and flows along with voltage and current magnitudes at specific buses.
The advent of the global positioning system (GPS) facilitated a precise timing signal across large geographical areas, thus enabling the revolutionary
technology of synchrophasors or phasor measurement units (PMUs) [58].
Recovering bus voltages given network parameters and the available measurements constitutes the critical task of power system state estimation.
This section formally states the problem, provides the Cramér-Rao bound
on the variance of any unbiased estimator, and reviews the Gauss-Newton
iterations. Solvers based on semidefinite relaxation and successive convex
approximations are subsequently explicated, and the section is wrapped up
with issues germane to PMUs.

3.1

Weighted Least-Squares Formulation

Consider M real-valued measurements {zm }M
m=1 related to the complex
power system state v through the model
zm = hm (v) + m

(14)

where hm (v) : CNb → R is a (non)-linear function of v, and m captures
the measurement noise and modeling inaccuracies. Collecting measurements
and noise terms in vectors z and  accordingly, the vector form of (14) reads
z = h(v) + 

(15)

for the mapping h : CNb → RM . Model (15) is instantiated for different
types of measurements next.
Traditionally, the system state v is expressed in polar coordinates, namely
nodal voltage magnitudes and angles. Then h(v) maps the 2Nb -dimensional
state vector to SCADA measurements through the nonlinear equations (6).
Expressing the states in polar form has been employed primarily due to two
reasons. First, the Jacobian matrix of h(v) is amenable to approximations.
5

Secondly, voltage magnitude measurements are directly related to states.
Nevertheless, due to recent computational reformulations, most of our exposition models voltages in the rectangular form. Then, as detailed in (12),
the m-th SCADA measurement zm involves the quadratic function of the
state hm (v) = vH Hm v for a Hermitian matrix Hm .
Expressing voltages in rectangular coordinates is computationally advantageous when it comes to synchrophasors too. As evidenced by (3)–(4),
PMU measurements feature linear mappings hm (v). If PMU measurements
are expressed in rectangular coordinates, the model in (14) simplifies to
z = Hv + 

(16)

for an M × Nb complex matrix H, and complex-valued noise . Following
the notation of (12)–(13), the linear measurement model of (16) can be
expressed in terms of real-valued quantities as
z̄ = (H∗ )v̄ + ¯.

(17)

The random noise vector  in (15) is usually assumed independent of
h(v), zero-mean and circularly symmetric, that is E[H ] = Σ and E[T ] =
0. The last assumption holds if for example the real and imaginary components of  are independent and have identical covariance matrices. This is
true for a PMU measurement, where the actual state lies at the center of a
spherically-shaped noise cloud on the complex plane.
Moreover, the entries of  are oftentimes assumed uncorrelated yielding
2 }) with σ 2 being the variance of the
a diagonal covariance Σ = dg({σm
m
m-th entry m . However, that may not always be the case. For example,
active and reactive powers at the same grid location are derived as products
between the readings of a current transformer and a potential transformer.
Further, noise terms may be correlated between the real and imaginary parts
of the same phasor in a PMU.
Adopting the weighted least-squares (WLS) criterion, power system state
estimation can be formulated as
minimize kΣ−1/2
(z − h(v))k22

v∈CNb

(18)

−1/2

where Σ
is the matrix square root of the inverse noise covariance matrix.
If the noise is independent across measurements, then (18) simplifies to
M
X
(zm − hm (v))2
minimize
.
2
v
σm
m=1

6

(19)

Either way, the PSSE task boils down to a (non)-linear least-squares (LS) fit.
When the mapping h(v) is linear or when the entries of  are uncorrelated,
the measurement model in (15) can be prewhitened. For example, the linear
measurement model z = Hv +  can be equivalently transformed to
Σ−1/2
z = (Σ−1/2 H)v + Σ−1/2 


(20)

−1/2

so that the associated noise Σ
 is now uncorrelated. To ease the presentation, the noise covariance will be henceforth assumed Σ = IM , yielding
v̂ := arg min
v

M
X

(zm − hm (v))2 .

(21)

m=1

For Gaussian measurement noise  ∼ N (0, IM ), the minimizer of (21) coincides with the maximum likelihood estimate (MLE) of v [29].

3.2

Cramér-Rao Lower Bound Analysis

According to standard results in estimation theory [29], the variance of
any unbiased estimator is lower bounded by the Cramér-Rao lower bound
(CRLB). Appreciating its importance as a performance benchmark across
different estimators, the ensuing result shown in the Appendix derives the
CRLB for any unbiased power system state estimator based on the so-termed
Wirtinger’s calculus for complex analysis [41].
Proposition 1. Consider estimating the unknown state vector v ∈ CNb
from the noisy SCADA data {zm }M
m=1 of (15), where the Gaussian measurement error m is independent across meters with mean zero and variance
2 . The covariance matrix of any unbiased estimator v̂ satisfies
σm
Cov(v̂)  [F† (v, v∗ )]1:Nb ,1:Nb

(22)

where the Fisher information matrix (FIM) is given as
" PM
#
PM
1
1
H
∗ ∗ H
2 (Hm v)(Hm v)
2 (Hm v)(Hm v )
m=1 σm
m=1 σm
∗
PM
F(v, v )= PM
.
1
1
∗ ∗
H
∗ ∗
∗ ∗ H
2 (Hm v )(Hm v)
2 (Hm v )(Hm v )
m=1 σm
m=1 σm
(23)
∗
In addition, matrix F(v, v ) has at least rank-one deficiency even when all
possible SCADA measurements are available.
Although rank-deficient, the pseudo-inverse of F(v, v∗ ) qualifies as a
valid lower bound on the mean-square error (MSE) of any unbiased estimator [63]. Rank deficiency of the FIM originates from the inherent voltage
7

angle ambiguity: SCADA measurements remain invariant if nodal voltages
are shifted globally by a unimodular phase constant. Fixing the angle of
the reference bus waives this issue. It is also worth stressing that the CRLB
in Prop. 1 is oftentimes attainable and benchmarks the optimal estimator
performance [63]. Having derived the CRLB for the PSSE task, our next
subsection deals with PSSE solvers.

3.3

Gauss-Newton Iterations

Consider for specificity model (16), though the real-valued model in (17)
or the model involving polar coordinates could be employed as well. When
the noise covariance matrix Σ = IM , the PSSE task in (18) reduces to the
nonlinear LS problem
minimize kz − h(v)k22
(24)
v∈CNb

for which the Gauss-Newton iterations are known to offer the “workhorse”
solution [3, Ch. 1], [1, Ch. 2]. According to the Gauss-Newton method,
the function h(v) is linearized at a given point vi ∈ CNb using Taylor’s
expansion as
h̃(v, vi ) := h(vi ) + Ji (v − vi )
where Ji := ∇h(vi ) is the M × Nb Jacobian matrix of h evaluated at vi ,
whose (m, n)-th entry is given by the Wirtinger derivative ∂hm /∂Vn ; see
e.g., [41] for Wirtinger’s calculus. The Gauss-Newton method subsequently
approximates the nonlinear LS fit in (24) with a linear one of h̃, and relies
on its minimizer to obtain the next iterate as
vi+1 ∈ arg min
v

= arg min
v

z − h̃(v, vi )
z − h(vi )

2

2

− 2(v − vi )H (Ji )H (z − h(vi ))

+ (v − vi )H (Ji )H Ji (v − vi ).

(25)

When matrix (Ji )H Ji is invertible, vi+1 can be found in closed form as

−1
vi+1 = vi + (Ji )H Ji (Ji )H (z − h(vi )).
(26)
The state estimate is iteratively updated using (26) until some stopping
criterion is satisfied.
If, on the other hand, the WLS cost (18) is minimized, the Gauss-Newton
−1/2
−1/2
iterations can be similarly obtained by treating Σ z as z and Σ h(v)
as h(v) in (24), yielding


i −1 i H −1
(J ) Σ (z − h(vi )).
(27)
vi+1 = vi + (Ji )H Σ−1
 J
8

It is well known that the pure Gauss-Newton iterations in (26) or (27)
may not guarantee convergence, which in fact largely depends on the starting
point v0 [3, Ch. 1.5]. A common way to improve convergence and ensure
descent of the cost in (18) consists of including a backtracking line search in
(27) to end up with


i −1 i H −1
vi+1 = vi + µi (Ji )H Σ−1
(J ) Σ (z − h(vi ))
 J

(28)

where the step size µi > 0 is found through the backtracking line search
rule [3, Ch. 1.2]. Due to its intimate relationship with ordinary gradient descent alternatives for nonconvex optimization however, this Gauss-Newton
iterative procedure can be trapped by local solutions [3, Ch. 1.5]. In a
nutshell, the grand challenge remains to develop PSSE solvers capable of attaining or approximating the global optimum at manageable computational
complexity. A few recent proposals in this direction are presented next.

3.4

Semidefinite Relaxation

A method to tackle the nonlinear measurement model that can convert the
PSSE problem of (21) to a convex semidefinite program (SDP) has been
introduced in [83], [84]. Consider first expressing each measurement in z
linearly in terms of the outer-product matrix V := vvH . In this way, the
quadratic models in (8), (10), and (11) can be transformed to linear ones in
terms of the matrix variable V. Thus, each noisy measurement in (14) can
be written as zm = vH Hm v + m = Tr(Hm V) + m . Rewriting the PSSE
task in (21) accordingly in terms of V reduces to
V̂1 := arg

min

V∈CNb ×Nb

M
X


2
zm − Tr(Hm V)

(29a)

m=1

s. to V  0, and rank(V) = 1

(29b)

where the positive semi-definite (PSD) and the rank-1 constraints jointly
ensure that for any V obeying (29b), there always exists a vector v such
that V = vvH .
Although zm and V are linearly related as in (29), nonconvexity is still
present in two aspects: (i) the cost function in (29a) has degree 4 in the
entries of V; and (ii) the rank constraint in (29b) is nonconvex. Aiming
for an SDP reformulation of (29), Schur’s complement lemma, see e.g., [5,
Appx. 5.5], can be leveraged to tightly bound each summand in (29a) using

9

an auxiliary variable χm > 0. Collecting all χm ’s in χ ∈ Rm , the problem
in (29) can be expressed as
{V̂2 , χ̂2 } := arg min 1T χ

(30a)

V, χ

s. to V  0, and rank(V) = 1
(30b)


χm
zm − Tr(Hm V)
 0, ∀m. (30c)
zm − Tr(Hm V)
1
The equivalence among all three SE problems (21), (29), and (30) has been
shown in [83], where their optimal solutions satisfy:

2
V̂1 = V̂2 = v̂v̂H , and χ̂2,m = zm − Tr(Hm V̂2 ) , ∀m.
(31)
The only source of nonconvexity in the equivalent SE problem of (30)
comes from the rank-1 constraint. Motivated by the technique of semidefinite relaxation (SDR); see e.g., the seminal work of [21], one can obtain the
following convex SDP upon dropping the rank constraint
{V̂, χ̂} := arg min 1T χ
V, χ

s. to V  0, and (30c).

(32a)
(32b)

For the SDR-PSSE formulation in (32), a few assumptions have been
made in [84] to establish its global optimality in a specific setup.
(as1) The graph G = (B, L) has a tree topology;
(as2) Every bus is equipped with a voltage magnitude meter; and
(as3) All measurements in z are noise-free, that is  = 0.
Proposition 2. Under (as1)-(as3), solving the relaxed problem (32) attains the global optimum of the original PSSE problem (30) or (21); that is,
rank(V̂) = 1.
Assumptions (as1)-(as3) may offer a close approximation of the realistic
PSSE scenario, thanks to characteristics of transmission systems such as
sparse connectivity, almost flat voltage profile, and high metering accuracy.
Although they do not hold precisely in realistic transmission systems, nearoptimality of the relaxed problem (32) has been numerically supported by
extensive tests [84]. A more crucial issue is to recover a feasible SE solution
from the relaxed problem (32), as V̂ is very likely to have rank greater than
1. This is possible either by finding the best rank-1 approximation to V̂ via
eigenvalue decomposition, or via randomization [46].
10

SDR endows SE with a convex SDP formulation for which efficient
schemes are available to obtain the global optimum using for example the
interior-point solver. The computational complexity for eigen-decomposition
is in the order of matrix multiplication, and thus negligible compared to
solving the SDP; see [46] and references therein. However, the polynomial
complexity of solving the SDP could be a burden for real-time power system monitoring, which motivates well the distributed implementation of
Section 4.2.

3.5

Penalized Semidefinite Relaxation

Building on an alternative formulation of the power flow problem, a penalized
version of the aforementioned SDP-based state estimator has been devised
in [48, 47, 78]. Commencing with the power flow task, it can be interpreted
as a particular instance of PSSE, where:
• Measurements (henceforth termed specifications) are noiseless;
• Excluding the reference bus, buses are partitioned into the subset BPV
for which active injections and voltage magnitudes are specified, and
the subset BPQ , for which active and reactive injections are specified.
The power flow task can be posed as the feasibility problem; that is,
find v ∈ CNb
s. to Pn = vH HPn v,
H

Qn = v HQn v,
Vn2 = vH HVn v,

(33)
∀n ∈ BPV ∪ BPQ
∀n ∈ BPQ
2
= V0 .
∀n ∈ BPV , and Vref

Using the SDP reformulation presented earlier, the power flow task can be
equivalently expressed as
find V ∈ CNb ×Nb

(34)

s. to zm = Tr(Hm V),

m = 1, . . . , 2Nb − 1

V  0, and rank(V) = 1
where the specifications (constraints) of (33) have been generically captured
b −1
by the pairs {(zm ; Hm )}2N
m=1 .
Although the optimization in (34) is non-convex, a convex relaxation
can be obtained by dropping the rank constraint. To promote rank-one
solutions, the feasibility problem is further reduced to [48]
V̂ := arg min Tr(H0 V)

(35)

V0

11

s. to zm = Tr(Hm V), m = 1, . . . , 2Nb − 1.

(36)

If the Hermitian matrix H0 is selected such that H0  0, rank(H0 ) = Nb −1,
and H0 1 = 0; then any state v close to the flat voltage profile 1 + j0 can
be recovered from the rank-one minimizer V̂ of (35); see [48, Thm. 2]. The
stated conditions exclude the case H0 = INb that would have led to the
nuclear-norm heuristic commonly used in low-rank matrix completion.
Spurred by this observation, the PSSE task can be posed as [47]
V̂µ := arg min Tr(H0 V) + µ
V0

M
X

fm zm − Tr(Hm V)



(37)

m=1

for some regularization parameter µ ≥ 0, where M now can be larger than
2Nb − 1. The second term in the cost of (37) is a data-fitting term ensuring
that the recovered state is consistent with the collected measurements based
on selected criteria. Two cases for fm of special interest are the LS fit
LS () := 2 , and the least-absolute value (LAV) one f LAV () := ||, ∀m =
fm
m
1, . . . , M . On the other hand, the first term in (37) can be understood as a
regularizer to promote rank-one solutions V̂µ ; see more details in [47].
In the noiseless setup, where all measurements comply with the model
zm = v0H Hm v0 , the minimization in (37) has been shown to possess a rankLS and f LAV ; see
one minimizer V̂µ = v̂µ v̂µH for all µ ≥ 0 under both fm
m
details in [47]. Interestingly, the solution v̂µ obtained under the LS fit does
not coincide with v0 for any µ ≥ 0, whereas the LAV solution v̂µ provides
the actual state v0 for a sufficiently large µ. Error bounds between V̂µ
and v0 v0H under the regularized LAV solution for noisy measurements are
established in [78].

3.6

Feasible Point Pursuit

The feasible point pursuit (FPP) method studied in [51] offers another
computationally manageable solver for approximating the globally optimal
PSSE. As a special case of the convex-concave procedure [44], FPP is an iterative algorithm for handling general nonconvex quadratically constrained
quadratic programs (QCQPs) [77]. It approximates the feasible solutions of
a nonconvex QCQP by means of a sequence of convexified QCQPs obtained
with successive convex inner-restrictions of the original nonconvex feasibility
set [77].
The first step in applying FPP to PSSE is a reformulation of (18) into

12

a standard QCQP [71, 70]
minimize

v∈CNb ,

χ∈RM

χT Σ−1 χ

(38a)

s. to vH Hm v ≤ zm + χm ,
H

v (−Hm ) v ≤ −zm + χm ,

1≤m≤M

(38b)

1≤m≤M

(38c)

where vector χ ∈ RM collects the auxiliary variables {χm ≥ 0}M
m=1 . For
power flow and power injection measurements, the corresponding Hermitian
measurement matrices {Hm } are indefinite in general; thus, both constraints
(38b) and (38c) are nonconvex. On the contrary, squared voltage magnitude
measurements relate to positive semidefinite matrices {Hm }, so that only
the related constraint (38c) is nonconvex. Either way, problem (38) is NPhard, and hence computationally intractable [57].
Using eigen-decomposition, every measurement matrix Hm can be expressed as the sum of a positive and a negative semidefinite matrix as
−
Hm = H+
m + Hm , so that the constraints in (38) are rewritten as
H −
vH H+
m v + v Hm v ≤ zm + χm

v

H

H+
mv

+v

H

H−
mv

≥ z m − χm

(39a)
(39b)

for m = 1, . . . , M . Observe now that since vH H−
m v is a concave function of
v, it is upper bounded by its first-order (linear) approximation at any point
vi ; that is,
i H −
i H − i
vH H−
m v ≤ 2<{(v ) Hm v} − (v ) Hm v .
The concave function vH (−H+
m )v can be upper bounded similarly.
Based on this observation, the FPP technique replaces the concave functions in the constraints of (39) with their linear upper bounds. The point of
linearization at every iteration is chosen to be the previous state estimate.
Specifically, initializing with v0 , the FPP produces the iterates
{vi+1 , χi+1 } := arg min χT Σ−1
 χ
v, χ≥0

(40)

i H −
i H − i
s. to vH H+
m v + 2<{(v ) Hm v} ≤ zm + (v ) Hm v + χm , ∀m
i H +
i H + i
vH H−
m v + 2<{(v ) Hm v} ≥ zm + (v ) Hm v − χm , ∀m.

At every iteration, the FPP technique solves the now convex QCQP in
(40). The procedure has been shown to globally converge to a stationary
point of the WLS formulation (18) of the PSSE task [71]. Extensions of the
developed FPP solver to cope with linear(ized) measurements and bad data
are straightforward; see [71] for details.
13

100

Mean-square error (p.u.)

Angle error (degrees)

102

102

10-2

Magnitude error (p.u.)

10-4
100

10-5

GN-based SE
SDR-based SE
FPP-based SE

5

10

15

20

25

10

100

10-1

10-2

10

30

GN-based SE
SDR-based SE
FPP-based SE
Cramer-Rao bound

1

-3

1

2

3

4

5

Types of measurements

bus index

Figure 1: Left: Voltage magnitude and angle estimation errors per bus for
the IEEE 30-bus system. Right: MSEs and CRLB versus types of measurements used for the IEEE 14-bus system using: (i) Gauss-Newton iterations;
ii) the SDR-based PSSE; and iii) the FPP-based PSSE.
Figure 1 compares Gauss-Newton iterations, the SDR-based solver, and
the FPP-based solver on the IEEE 14- and 30-bus systems [64]. The actual nodal voltage magnitudes and angles were generated uniformly at random over [0.9, 1.1] and [−0.4π, 0.4π], respectively. Independent zero-mean
Gaussian noise with standard deviation 0.05 for power and 0.02 for voltage
measurements was assumed, and all reported results were averaged over 100
independent Monte Carlo realizations. The measurements for the IEEE 30bus system include all nodal voltage magnitudes and the active power flows
at both sending and receiving ends. The left panel of Fig. 1 depicts that
the magnitude and angle estimation errors attained by the FPP solver are
consistently below its competing alternatives [71].
The second experiment examines the MSE performance of the three approaches relative to the CRLB of (22) for the IEEE 14-bus test system.
Initially, all voltage magnitudes as well as all sending- and receiving-end active power flows were measured, which corresponds to the base case 1 in the
x-axis of the right panel of Fig. 1. To show the MSE performance relative for
an increasing number of measurements, additional types of measurements
were included in a deterministic manner. All types of SCADA measurements
were ordered as {Vn2 , Pkn , Pnk , Qkn , Qnk , Pn , Qn }. Each x-axis value in the
right panel of Fig. 1 implies that the number of ordered types of measurements was used in the experiment to obtain the corresponding MSEs.

14

3.7

Synchrophasors

To incorporate synchrophasors into the PSSE formulation, let ζn = Φn v+εn
collect the noisy PMU data at bus n [cf. (16)]. The related measurement
matrix is Φn , and the measurement noise εn is assumed to be complex zeromean Gaussian, independent from the noise  in legacy meters and across
buses. Following the normalization convention in (21), the noise vector εn
is assumed prewhitened, such that all PMU measurements exhibit the same
accuracy. The PSSE task now amounts to estimating v given both z and
{ζn }n∈P , where P ⊆ B denotes the subset of the PMU-instrumented buses.
Hence, the MLE cost in (21) needs to be augmented by the log-likelihood
induced by PMU data as
v̂ := arg min
v

M
X

(zm − hm (v))2 +

m=1

X

kζn − Φn vk22 .

(41)

n∈P

The SDR methodology is again well motivated to convexify the augmented
PSSE problem (41) into
X

H
minimize 1T χ +
Tr(ΦH
(42a)
n Φn V) − 2<{ζn Φn v}
V, v, χ

n∈P


s. to

V v
vH 1


 0, and (30c).

(42b)

By Schur’s complement, the left SDP constraint in (42b) can be expressed
equivalently as V  vvH . If the latter constraint is enforced with equality,
the matrix V becomes rank-one. Imposing a rank-one constraint in (42)
renders it equivalent to the augmented PSSE task of (41). The SDP here
also offers the advantages of (32), in terms of the near-optimality and the
distributed implementation deferred to Section 4.2. To recover a feasible
solution, one can again use the best rank-1 approximation or adopt the
randomization technique as elucidated in [84].
Alternatively, the two types of measurements can be jointly utilized upon
interpreting the SCADA-based estimate as a prior for PMU-based estimation [58, 32]. Specifically, if v̂s is the SCADA-based estimate, the prior
probability density function of the actual state can be postulated to be a
circularly symmetric complex Gaussian with mean v̂s and covariance Σ̂s .
Given PMU data and the SCADA-based prior, the state can be estimated
following a maximum a-posteriori probability (MAP) approach as
X
v̂ := arg min (v − v̂s )H Σ̂−1
kζn − Φn vk22
(43)
s (v − v̂s ) +
v

n∈P

15

where the first summand is the negative logarithm of the prior distribution,
and the second one is the negative log-likelihood from the PMU data. In
essence, the approach in (43) treats the SCADA-based estimate as pseudomeasurements relying on the model v̂s = v +  with circularly symmetric
zero-mean noise having E[H ] = Σ̂s .

4

Distributed Solvers

Upcoming power system requirements call for decentralized solvers. Measurements are now collected at much finer spatio-temporal scales and the
number of states increases exponentially as monitoring schemes extend to
low-voltage distribution grids [13]. Tightly interconnected power systems
call for the close coordination of regional control centers [22], while operators and utilities perform their computational operations on the cloud.
This section reviews advances in distributed PSSE solvers. As the name
suggests, distributed PSSE solutions spread the computational load across
different processors or control centers to speed up time, implement memoryintensive tasks, and/or guarantee privacy. A network of processors may be
coordinated by one or more supervising control centers in a hierarchical
fashion, or completely autonomously, by exchanging information between
processors. To clarify terminology, the latter architecture will be henceforth
identified as decentralized.
Distributed solvers with a hierarchical structure have been proposed
since the statistical formulation of PSSE [61, Part III]. Different versions
of this original scheme were later developed in [10], [27], [80], [22], [38].
Decentralized schemes include block Jacobi iterations [43], [9]; an approximate algorithm building on the related optimality conditions [18]; or matrixsplitting techniques for facilitating matrix inversion across areas running
Gauss-Newton iterations [53]. Most of the aforementioned approaches presume local identifiability (i.e., each area is identifiable even when shared
measurements are excluded) or their convergence is not guaranteed. Assuming a ring topology, every second agent updates its state iteratively through
the auxiliary problem principle in [16]. Local observability is waived in the
consensus-type solver of [75], where each control center maintains a copy of
the entire high-dimensional state vector resulting in slow convergence. For
a relatively recent review on distributed PSSE solves, see also [23].

16

Figure 2: Left: The IEEE 14-bus system partitioned into four areas [64, 38].
Dotted lassos show the buses belonging to area state vectors vk ’s. PMU bus
voltage (line current) measurements depicted by green circles (blue squares).
Right: The matrix structure for the left system with the distributed SDR
solver: green square denotes the overall V, while dark and light blue ones
correspond to the four area submatrices {Vk } and their overlaps.

4.1

Distributed Linear Estimators

Consider an interconnected system partitioned in K areas supervised by
separate control centers. Without loss of generality, an area may be thought
of as an independent system operator region, a balancing authority, a power
distribution center, or a substation [58]. Area k collects Mk measurements
obeying the linear model
zk = Hk vk + k
(44)
where vector vk ∈ CNk collects the system states related to zk through
the complex matrix Hk . The random noise vector k is zero-mean with
identity covariance upon prewhitening, if measurements are uncorrelated
across areas. The model in (44) is exact for PMU measurements, but it
may also correspond to a single Gauss-Newton iteration as explained in
Section 3.3.
Performing PSSE locally at area k amounts to solving
minimize fk (vk )
vk ∈Xk

(45)

where the convex set Xk captures possible prior information, such as zero17

injection buses or short circuits [55], [1]. If fk (vk ) = kzk − Hk vk k22 /2, the
minimizer of (45) is the least-squares estimate (LSE) of vk , which is also
the MLE of vk for Gaussian k .
As illustrated in Fig. 2, the per-area state vectors {vk }K
k=1 overlap partially. Although area 4 supervises buses {9, 10, 14}, it also collects the current reading on lines (10, 11). Thus, its state vector v4 extends to bus
{11} that is nominally supervised by area 3. To setup notation, define Skl
as the shared states for a pair of neighboring areas (k, l). Let also vk [l]
(vl [k]) denote the sub-vector of vk (vl ) consisting of their overlapping variables ordered as they appear in v. For example, v3 [4] = v4 [3] contain the
bus voltages of {11}. Solving the K problems in (45) separately is apparently suboptimal since the estimates of shared states will disagree, tie-line
measurements have to be ignored, and boundary states may thus become
unobservable.
Coupling the per-area PSSE tasks can be posed as
minimize
{vk ∈Xk }

K
X

fk (vk )

(46)

k=1

s. to vk [l] = vl [k],

∀l ∈ Bk , ∀k

where Nk is the set of areas sharing states with area k. The equality constraints of (46) guarantee consensus over the shared variables. Although
(46) is amenable to decentralized implementations (cf. [16]), areas need a
coordination protocol for their updates. To enable a truly decentralized solution, we follow the seminal approach of [60, 82]. An auxiliary variable vkl
is introduced per pair of connected areas (k, l); the symbols vkl and vlk are
used interchangeably. The optimization in (46) can then be written as
minimize

{vk ∈Xk }, {vkl }

K
X

fk (vk )

(47)

k=1

s. to vk [l] = vkl ,

∀l ∈ Bk , k = 1, . . . , K.

Problem (47) can be solved using the alternating direction method of
multipliers (ADMM) [4], [20]. In its general form, ADMM tackles convex
optimization problems of the form
minimize f (x) + g(z)
x∈X , z∈Z

s. to Ax + Bz = c

18

(48a)
(48b)

for given matrices and vectors (A, B, c) of proper dimensions. Upon assigning a Lagrange multiplier λ for the coupling constraint in (48a), the (x, z)
minimizing (48a) are found through the next iterations for some µ > 0
µ
kAx + Bzi − c + λi k22
2
µ
:= arg min g(z) + kAxi+1 + Bz − c + λi k22
z∈Z
2
:= λi + Axi+1 + Bzi+1 − c.

xi+1 := arg min f (x) +

(49a)

zi+1

(49b)

x∈X

λi+1

(49c)

Towards applying the ADMM iterations to (46), identify variables {vk }
as x in (48a) and {vkl } as z with g(z) = 0. Moreover, introduce Lagrange
multipliers λk,l for each constraint in (47). Observe that λk,l and λl,k correspond to the distinct constraints vk [l] = vkl and vl [k] = vkl , respectively.
According to (49a), the per-area state vectors {vk } can be updated separately as
µX
i
vki+1 := arg min fk (vk ) +
kvk [l] − vkl
+ λik,l k22 .
(50)
vk ∈Xk
2
l∈Bk

From (49b) and assuming every state is shared by at most two areas, the
auxiliary variables vkl can be readily found in closed form given by
i+1
vkl
=


1 i+1
vk [l] + vli+1 [k] + λik,l + λil,k
2

(51)

while the two related multipliers are updated as
i+1
i+1
i
λi+1
k,l := λk,l + (vk [l] − vkl )

λi+1
l,k

:=

λil,k

+

(vli+1 [k]

−

i+1
).
vkl

(52a)
(52b)

i+1
Adding (52) by parts and combining it with (51) yields λi+1
k,l = −λl,k at
all iterations i if the multipliers are initialized at zero. Hence, the auxiliary
variable vkl ends up being the average of the shared states; that is,
i+1
vkl
=


1 i+1
vk [l] + vli+1 [k] .
2

(53)

To summarize, at every iteration i:
(i) Each control area solves (50). If fk (vk ) is the LS fit and for unconstrained problems, the per-area states are updated as the LSEs using
legacy software. The second summand in (50) can be interpreted as
pseudo-measurements on the shared states forcing them to consent
across areas.
19

−2

10

Error curves averaged over areas

LSE
D−RPSSE
−3

10

−4

10

−5

10

−6

10

0

10

20

30

40

50

60

70

80

Iteration

P300 t
P
t
Figure 3: Average error curves 300
k=1 ek,o /300
k=1 ek,c /300 (bottom) and
(top) for the LSE, and its robust counterpart D-RPSSE (see Section 5.1) on
a 4,200-bus grid.
(ii) Neighboring areas exchange their updated shared states. This step involves minimal communication, and no grid models need to be shared.
Every area updates its copies of the auxiliary variables vkl using (53).
(iii) Every area updates the Lagrange multipliers λkl based on the deviation of the local from the auxiliary variable as in (52).
For convex pairs {fk (vk ), Xk }K
k=1 , the aforementioned iterates reach the
optimal cost in (46), under mild conditions. If the overall power system is
observable, the ADMM iterates converge to the unique LSE. The approach
has been extended in [33] for joint PSSE and breaker status verification.
The decentralized algorithms were tested on a 4,200-bus power grid synthetically built from the IEEE 14- and 300-bus systems. Each of the 300
buses of the latter was assumed to be a different area, and was replaced by
a copy of the IEEE 14-bus grid. Additionally, every branch of the IEEE
300-bus grid was an inter-area line whose terminal buses are randomly selected from the two incident to this line areas. Two performance metrics
were adopted: the per area error to the centralized solution of (46), denoted
(c)
by etk,c :=kvk − vkt k2 /Nk , and the per-area error to the true underlying
state defined as etk,o :=kvk − vkt k2 /Nk . Figure 3 shows the corresponding
error curves averaged over 300 areas. The decentralized LSE approached
the underlying state at an accuracy of 10−3 in approximately 10 iterations
20

or 6.2 msec on an Intel Duo Core @ 2.2 GHz (4GB RAM) computer using
MATLAB; while the centralized LSE finished in 93.4 msec.

4.2

Distributed SDR-based Estimators

Although the SDR-PSSE approach incurs polynomial complexity √
when im4
plemented as a convex SDP, its worst-case complexity is still O(M Nb log(1/))
for a given solution accuracy  > 0 [46]. For typical power networks, the
number of measurements M is on the order of the number of buses Nb , and
thus the worst-case complexity becomes O(Nb4.5 log(1/)). This complexity
could be prohibitive for large-scale power systems, which motivates accelerating the SDR-PSSE method using distributed parallel implementations.
Following the area partition in Fig. 2, the m-th measurement per area k
can be written as
zk,m = hk,m (vk ) + k,m = Tr(Hk,m Vk ) + k,m , ∀k, m
where Vk denotes a submatrix of V formed by extracting the rows and
columns corresponding to buses in area k; and likewise for each Hk,m . Due
to the overlap among the subsets of buses, the outer-product Vk of area k
overlaps also with Vl for each neighboring area l ∈ Nk , as shown in Fig. 2.
By reducing the measurements atP
area k to submatrix Vk , one can de2
k
fine the PSSE error cost fk (Vk ) := M
m=1 [zk,m − Tr(Hk,m Vk )] per area
k, which only involves the local matrix Vk . Hence, the centralized PSSE
problem in (32) becomes equivalent to
V̂ = arg min

K
X

V0

fk (Vk ).

(54)

k=1

This equivalent formulation effectively expresses the overall PSSE cost as the
superposition of each local cost fk . Nonetheless, even with such a decomposable cost, the main challenge to implement (54) in a distributed manner
lies in the PSD constraint that couples the overlapping local matrices {Vk }
(cf. Fig. 2). If all submatrices {Vk } were non-overlapping, the cost would
be decomposable as in (54), and the PSD of V would boil down to a PSD
constraint per area k, as in
V̂ = arg min

{Vk 0}

K
X

fk (Vk ).

(55)

k=1

Similar to PSSE for linearized measurements in (46), the formulation in
(55) can be decomposed into sub-problems, thanks to the separable PSD
21

constraints. It is not always equivalent to the centralized (54) though, because the PSD property of all submatrices does not necessarily lead to a
PSD overall matrix. Nonetheless, the decomposable problem (55) is still
a valid SDR-PSSE reformulation, since with the additional per-area constraints rank(Vk ) = 1, it is actually equivalent to (30). While it is totally
legitimate to use (55) as the relaxed SDP formulation for (30), the two
relaxed problems are actually equivalent under mild conditions.
The fresh idea here is to explore valid network topologies to facilitate
such PSD constraint decomposition. To this end, it will be instrumental to leverage results on completing partial Hermitian matrices to obtain
PSD ones [24]. Upon obtaining the underlying graph formed by the specified entries in the partial Hermitian matrices, these results rely on the
so-termed graph chordal property to establish the equivalence between the
positive semidefiniteness of the overall matrix and that of all submatrices
corresponding to the graph’s maximal cliques. Interestingly, this technique
was recently used for developing distributed SDP-based optimal power flow
(OPF) solvers in [28, 12, 42].
Construct first a new graph B 0 over B, with all its edges corresponding
to the entries in {Vk }. The graph G 0 amounts to having all buses within
each subset Nk to form a clique. Furthermore, the following are assumed:
(as4) The graph with all the control areas as nodes, and their edges defined
by the neighborhood subset {Nk }K
k=1 forms a tree.
(as5) Each control area has at least one bus that does not overlap with any
neighboring area.
Proposition 3. Under (as4)-(as5), the two relaxed problems (54) and (55)
are equivalent.
Proposition 3 can be proved by following the arguments in [84] to show
that the entire PSD matrix V can be “completed” using only the PSD
submatrices Vk . The key point is that in most power networks even those
not obeying (as4) and (as5), (55) can achieve the same accuracy as the
centralized one. At the same time, decomposing the PSD constraint in (55)
is of paramount importance for developing distributed solvers. One can
adopt the consensus reformulation to design the distributed solver for (55) as
in (46) of Section 4.1. Accordingly, the ADMM iterations can be employed
to solve (55) through iterative information exchanges among neighboring
areas, and this is the basis of the distributed SDR-PSSE method.
This distributed SDR-PSSE method was tested on the IEEE 118-bus
system using the three-area partition in [31]. All three areas measure their
22

1

10

Area 1
Area 2
Area 3

Local matrix error

10

0

Area 1
Area 2
Area 3

10

Local estimation error

0

−1

10

−2

−1

10

10

−3

10

0

20

40

60

0

80

20

40

60

80

Iteration index

Iteration index

Figure 4: (Left) Per area state matrix error and (Right) state vector estimation error, versus the number of ADMM iterations for the distributed
SDR-PSSE solver using the IEEE 118-bus system.
local bus voltage magnitudes, as well as real and reactive power flow levels
at all lines. The overlaps among the three areas form a tree communication
graph used to construct the area-coupling constraints. To demonstrate convergence of the ADMM iterations to the centralized SE solution V̂ of (55),
i − V̂ k is plotted versus the
the local matrix Frobenius error norm kV̂(k)
(k) F
iteration index i in the left panel of Fig. 4 for every control area k. Clearly,
all local iterates converge to (approximately with a linear rate) their counterparts in the centralized solution. As the task of interest is to estimate the
i − v k is also
voltages, the local estimation error for the state vector kv̂(k)
(k) 2
depicted in the right panel of Fig. 4, where v̂(k) is the estimate of bus volti using the eigen-decomposition
ages at area k obtained from the iterate V̂(k)
method. Interestingly, the estimation error costs converge within the estimation accuracy of around 10−2 after about 20 iterations (less than 10
iterations for area 1), even though the local matrix has not yet converged.
In addition, these error costs decrease much faster in the first couple of iterations. This demonstrates that even with only a limited number of iterations,
the PSSE accuracy can be greatly boosted in practice, which in turn makes
inter-area communication overhead more affordable.

5

Robust Estimators and Cyber Attacks

Bad data, also known as outliers in the statistics parlance, can challenge
PSSE due to communication delays, instrument mis-calibration, and/or line

23

parameter uncertainty. In today’s cyber-enabled power systems, smart meter and synchrophasor data could be also purposefully manipulated to mislead system operators. This section reviews conventional and contemporary
approaches to coping with outliers.

5.1

Bad Data Detection and Identification

Bad data processing in PSSE relies mainly on the linear measurement model
z = Hv + , where H ∈ RM ×N . Recall that this model is exact for PMU
measurements [cf. (16)–(17)], but approximate per Gauss-Newton iteration
or under the linearized grid model. In addition, the aforementioned model
assumes real-valued states and measurements, slightly abusing the symbols
introduced in (17). This is to keep the notation uncluttered and cover both
cases of exact and inexact grid models. Albeit the nominal measurement
noise vector is henceforth assumed zero-mean with identity covariance, results extend to colored noise as per (20).
To capture bad data, the measurement model is now augmented as
z = Hv + o + 

(56)

where o ∈ RM is an unknown vector whose m-th entry om is deterministically non-zero only if zm is a bad datum [30, 39, 15]. Therefore, vector o is
sparse, i.e., many of its entries are zero. Under this outlier-cognizant model
in (56), the unconstrained LSE given as v̂LSE = (HT H)−1 HT z, yields the
residual error
r := z − Hv̂LSE = Pz = P(o + )
(57)
with P := IM − H(HT H)−1 HT being the so-called projection matrix onto
the orthogonal subspace of range(H). The last equality in (57) stems from
the fact that PH = 0. As a projection matrix, P is idempotent, that is,
P = P2 ; Hermitian PSD with (M − N ) eigenvalues equal to one and N
zero eigenvalues; while its diagonal entries satisfy Pm,m ∈ [0, 1] for m =
1, . . . , M ; see e.g., [6].
For  ∼ N (0, IM ), it apparently holds that P ∼ N (0, P). The meansquared residual error is (see also [39] for its Bayesian counterpart)
E[krk22 ] = E[kPk22 ] + kPok22 = (M − N ) + kPok22 .

(58)

In the absence of bad data, or if o ∈ range(H), the squared residual error
follows a χ2 distribution with mean (M − N ). The χ2 -test compares krk22
against a threshold to detect the presence of bad data [1, 55].

24

Finding both v and o from measurements in (56) may seem impossible, given that the number of unknowns exceeds the number of equations.
Leveraging the sparsity of o though, interesting results can be obtained [30].
If τ0 bad data are expected, one would ideally wish to solve
o
n1
{v̂, ô} ∈ arg min
(59)
kz − Hv − ok22 : kok0 ≤ τ0 .
v, o
2
But the `0 -(pseudo) norm kok0 counting the number of non-zero entries of
o, renders (59) NP-hard in general; see also Definition 2 later in Section 5.2.
For the special case of τ0 = 1, problem (59) can be efficiently handled.
Consider the scenario where the only non-zero entry of ô is the m-th one,
and denote the related v̂ minimizer by v̂(m) . Apparently, the m-th entry of
the ô minimizer is ôm := zm − hTm v̂(m) . This choice nulls the m-th residual
(zm − hTm v̂(m) − ôm = 0). With the m-th residual zeroed, the cost in (59)
becomes kr(m) k22 := kz(m) − H(m) v̂(m) k22 , where z(m) is obtained from z
upon dropping its m-th entry and H(m) by removing the m-th row of H.
The problem in (59) is then equivalent to
minimize
m

1
kr k2 .
2 (m) 2

(60)

Problem (60) can be solved by exhaustively finding all M LSEs excluding
one measurement at a time. Fortunately, a classical result from the adaptive
filtering literature relates the error kr(m) k22 to the error attained using all
outlier-free measurements krk22 := kPzk22 ; see e.g., [25, Ch. 9]
krk22 = kr(m) k22 + rm ôm .

(61)

The same result links the a-posteriori error rm to the a-priori error ôm as
rm = Pm,m ôm . Through these links, solving (60) is equivalent to:
|rm |
rmax := maximize p
.
m
Pm,m

(62)

In words, a single bad datum can be identified by properly normalizing the
entries of the original residual vector r = Pz.
Interestingly, the task in (62) coincides with the largest normalized residual (LNR) test that compares rmax to a prescribed threshold to identify a
single bad datum [1, Sec. 5.7]. The p
threshold is derived after recognizing
that in the absence of bad data, rm / Pm,m is standard normal for all m.
The LNR test does not generalize for multiple bad data and problem
(59) becomes computationally intractable for larger τ0 ’s. Heuristically, if
25

a measurement is deemed as outlying, PSSE is repeated after discarding
this bad datum, the LNR test is re-applied, and the process iterates till no
corrupted data are identified. Alternatively, the least-median squares and
the least-trimmed squares estimators have provable breakdown points and
superior efficiency under Gaussian data; see e.g., [52] and references therein.
Nevertheless, their complexity scales unfavorably with the network size.
Leveraging compressed sensing [7], a practical robust estimator can be
found if the `0 -pseudonorm is surrogated by the convex `1 -norm as [30, 31]
n1
o
minimize
kz − Hv − ok22 : kok1 ≤ τ1
(63)
v, o
2
for a preselected constant τ1 > 0, or in its Lagrangian form
{v̂, ô} ∈ arg min
v, o

1
kz − Hv − ok22 + λkok1
2

(64)

for some tradeoff parameter λ > 0. The estimates of (64) offer joint state
estimation and bad data identification. Even when some measurements are
deemed as corrupted, their effect has been already suppressed. The optimization task in (64) can be handled by off-the-shelf software or solvers
customized to the compressed sensing setup. When λ → ∞, the minimizer
ô becomes zero, and thus v̂ reduces to the LSE. On the contrary, by letting λ → 0+ , the solution v̂ coincides with the least-absolute value (LAV)
estimator [40, 6, 17, 68]; presented earlier in (37), namely
v̂LAV := arg min kz − Hvk1 .
v

(65)

For finite λ > 0, the v̂ minimizer of (64) is equivalent to Huber’s Mestimator; see [30] and references therein. Based on this connection and
for Gaussian , parameter λ can be set to 1.34, which makes the estimator 95% asymptotically efficient for outlier-free measurements [50, p. 26].
Huber’s estimate can be alternatively expressed as the v-minimizer of [49],
minimizev, ω 21 kωk22 + λkv − Hv − ωk1 . The bad data identification performance of this minimization has been analyzed in [76].
Table 1 compares several bad data analysis methods on the IEEE 14bus grid of Fig. 2 under the next four scenarios: (S0) no bad data; (S1)
bad data on line (4, 7); (S2) bad data on line current (4, 7) and bus voltage 5; and (S3) bad data on bus voltage 5 and line currents (4, 7) and
(10, 11). In all scenarios, bad data are simulated by multiplying the real
and imaginary parts of the actual measurement by 1.2. The performance
metric here is the `2 -norm between the true state and the PSSE, which is
26

Table 1: Mean-Square Estimation Error in the Presence of Bad Data
Method

GA-LSE

LSE

LNRT

Huber’s

0.0278
0.0313
0.0336
0.0367

0.0278
0.0318
0.1431
0.1434

0.0286
0.0331
0.0404
0.0407

0.0281
0.0322
0.0390
0.0390

(S0)
(S1)
(S2)
(S3)

averaged over 1,000 Monte Carlo runs. Four algorithms were tested: (a)
an ideal but practically infeasible genie-aided LSE (GA-LSE), which ignores
the corrupted measurements; (b) the regular LSE; (c) the LNR test-based
(LNRT) estimator with the test threshold set to 3.0 [1]; and (d) Huber’s
estimator of (64) with λ = 1.34. For (S0)-(S1), the estimators perform
comparably. The few corrupted measurements in (S2)-(S3) can deteriorate
LSE’s performance, while Huber’s estimator performs slightly better than
LNRT. Computationally, Huber’s estimator was run within 1.3 msec, while
the LNRT required 1.5 msec. The computing times were also measured for
the IEEE 118-bus grid without corrupted data. Interestingly, the average
time on the IEEE 118-bus grid without corrupted data are 3.2 msec and
81 msec, respectively.
Towards a robust decentralized state estimator, the ADMM-based framework of Section 4.1 can be engaged here too. If the measurement model for
the k-th area is zk = Hk vk + ok + k , the centralized problem boils down to
minimize

{vk ∈Xk , ok }

K
X
1
k=1

2

kzk − Hk vk − ok k22 + λkok k1 .

(66)

To allow for decentralized implementation, the optimization in (66) can be
reformulated as
minimize

K
X
1
k=1

2

kzk − Hk vk − ok k22 + λkωk k1

(67a)

over {vk ∈ Xk , ok , ωk }, {vkl }

(67b)

s. to vk [l] = vkl , for all l ∈ Bk , k = 1, . . . , K.

(67c)

ok = ωk , for all k = 1, . . . , K.

(67d)

As in Section 4.1, the constraints in (67c) and the auxiliary variables {vkl }
enforce consensus of shared states. On the other hand, the variables {ok }
27

are duplicated as {ωk } in (67d). Then, variables {vk , ok } are put together
in the x-update of ADMM in (49a), whereas {vkl , ωk } fall into the z-update
in (49b). In this fashion, costs are separable over variable groups, and the
minimization involving the `1 -norm enjoys a closed-form solution expressed
in terms of the soft thresholding operator [31].

5.2

Observability and Cyber Attacks

In the cyber-physical smart grid context, bad data are not simply unintentional errors, but can also take the form of malicious data injections [54].
Amid these challenges, the intertwined issues of critical measurements and
stealth cyber-attacks on PSSE are discussed next.
It has been tacitly assumed so far that the power system is observable. A
power system is observable if distinct states v 6= v0 are mapped to distinct
measurements h(v) 6= h(v0 ) under a noiseless setup. Equivalently, if the
so-called measurement distance function is defined as [81]
D(h) := minimize
kh(v) − h(v0 )k0
0
v6=v

(68)

the power system is observable if and only if D(h) ≥ 1. Given the network
topology and the mapping h(v), the well-studied topic of observability analysis aims at determining whether the system state is uniquely identifiable,
at least locally in a neighborhood of the current estimate [1, Ch. 4]. If not,
mapping observable islands, meaning maximally connected sub-grids with
observable internal flows, is important as well.
Observability analysis relies on the decoupled linearized grid model, and
is accomplished through topological or numerical tests [8, 56]. Apparently,
under the linear or linearized model h(v) = Hv, the state v is uniquely
identifiable if and only if H is full column-rank. Phase shift ambiguities can
be waived by fixing the angle at a reference bus.
In the presence of bad data and/or cyber attacks, observability analysis
may not suffice. Consider the noiseless measurement model z = h(v) + o,
where the non-zero entries of vector o correspond to bad data or compromised meters; and let us proceed with the following definitions.
Definition 1 (Observable attack [81]). The attack vector o is deemed as
observable if for every state v there is no v0 6= v, such that h(v)+o = h(v0 ).
Definition 2 (Identifiable attack [81]). The attack vector o is identifiable
if for every v there is no (v0 , o0 ) with v0 6= v and ko0 k0 ≤ kok0 , such that
h(v) + o = h(v0 ) + o0 .
28

If the outlier vector o is observable, the operator can tell that the collected measurements do not correspond to a system state, and can hence
decide that an attack has been launched. Nevertheless, the attacked meters
can be pinpointed only under the stronger conditions of Definition 2.
The resilience of the measurement mapping h(v) against attacks can be
characterized through D(h) in (68): The maximum number of counterfeited
meters for an attack to be observable is Ko = D(h)−1 and to be identifiable,
it is Ki = b D(h)−1
c; see [81, 76]. Here, the floor function bxc returns the
2
greatest integer less than or equal to x.
Consider the linear mapping h(v) = Hv. Measurement m is termed
critical if once removed from the measurement set, it renders the power
system non-identifiable. In other words, although H is full column-rank, its
submatrix H(m) is not. It trivially follows that D(h) = 1, and the system
operator can be arbitrarily misled even if only measurement m is attacked.
Due to the typically sparse structure of H, critical measurements or multiple
simultaneously corrupted data do exist [1]. It was pointed out in [45] that if
an attack o can be constructed to lie in the range(H), it comprises a ‘stealth
attack.’ Although finding D(h) is not trivial in general, a polynomial-time
algorithm leveraging a graph-theoretic approach is devised in [39].

6

Power System State Tracking

The PSSE methods reviewed so far ignore system dynamics and do not
exploit historical information. Dynamic PSSE is well motivated thanks to its
improved robustness, observability, and predictive ability when additional
temporal information is available [26]. Recently proposed model-free and
model-based state tracking schemes are outlined next.

6.1

Model-free State Tracking via Online Learning

In complex future power systems, one may not choose to explicitly commit
to a model for the underlying system dynamics. The framework of online
convex optimization (OCO), particularly popular in machine learning, can
account for unmodeled dynamics and is thus briefly presented next [62].
The OCO model considers a multi-stage game between a player and an
adversary. In the PSSE context, the utility or the system operator assumes
the role of the player, while the loads and renewable generations can be
viewed as the adversary. At time t, the player first selects an action Vt from
a given action set V, and the adversary subsequently reveals a convex loss

29

function ft : V → R. In this round, the player suffers a loss ft (Vt ). The
ultimate goal for the player is to minimize the regret Rf (T ) over T rounds:
Rf (T ) :=

T
X

ft (Vt ) − minimize
V∈V

t=1

T
X

ft (V).

(69)

t=1

The regret is basically the accumulated cost incurred
the player relative
Pby
T
0
to that by a single fixed action V := arg minV∈V t=1 ft (V). This fixed
action is selected with the advantage of knowing the loss functions {ft }Tt=1
in hindsight. Under appropriate conditions, judiciously designed online optimization algorithms can achieve sublinear regret; that is, Rf (T )/T → 0 as
T → +∞.
Building on the SDR-PSSE formulation of Section 3.4, the ensuing method
considers streaming data for real-time PSSE. The data referring to and colt
lected over the control period t are {(zmt ; Hmt )}M
mt =1 with t = 1, . . . , T .
The number and type of measurements can change over time, while the matrix corresponding to measurement m may change over time as indicated by
t
{Hmt }M
mt =1 due to topology reconfigurations. The online PSSE task can be
now formulated as
T
X
minimize
ft (V)
(70)
V0

t=1

PMt

where ft (V) := mt =1 [zmt − Tr(Hmt V)]2 . Online PSSE aims at improving
the static estimates by capitalizing on previous measurements as well as
tracking slow time-varying variations in generation and demand.
Minimizing the cost in (70) may be computationally cumbersome for
real-time implementation. An efficient alternative based on online gradient
descent amounts to iteratively minimizing a regularized first-order approximation of the instantaneous cost instead [37]
Vt+1 := arg min Tr(VH ∇ft (Vt )) +
V0

1
kV − Vt k2F
2µt

(71)

for t = 1, . . ., and suitably selected step sizes µt > 0. Interestingly, the
optimization in (71) admits a closed-form solution given by
Vt+1 = ProjS+ [Vt − µt ∇ft (Vt )]

(72)

with ProjS+ denoting the projection onto the positive semidefinite cone,
which can be performed using eigen-decomposition followed by setting negative eigenvalues to zero. It is worth mentioning that the online PSSE in
30

(72) enjoys sublinear regret [37]. Upon finding Vt , a state estimate vt can be
obtained by eigen-decomposition or randomization as in Section 3.4. With
an additional nuclear-norm regularization term promoting low-rank solutions in (71), online ADMM alternatives were devised in [36]. Interestingly,
online learning tools has recently been advocated for numerous real-time
energy management tasks in [35], [34], [69].

6.2

Model-based State Tracking

Although the previous model-free solver can recover slow time-varying states,
model-based approaches facilitate tracking of fast time-varying system states.
A typical state-space model for power system dynamics is [65]
vt+1 = Ft vt + gt + ωt

(73a)

zt = h(vt ) + t

(73b)

where Ft denotes the state-transition matrix, gt captures the process mismatch, and wt is the additive noise. The nonlinear mapping h(·) comes from
conventional SCADA measurements. Values {(Ft , gt )} can be obtained in
real-time using for example Holt’s system identification method [11]. Two
common dynamic tracking approaches to cope with the nonlinearity in the
measurement model of (73b) include the (extended or unscented) Kalman
filters and moving horizon estimators [65, 14, 26, 73, 72], and they are outlined in order next.
The extended Kalman filter (EKF) handles the nonlinearity by linearizing h(v) around the state predictor. To start, let v̂t+1|t stand for the predicted estimate at time t + 1 given measurements {zτ }tτ =1 up to time t.
Let also v̂t+1|t+1 be the filtered estimate given measurements {zτ }t+1
τ =1 . If
the noise terms ωt and t in (73) are assumed zero-mean Gaussian with
known covariance matrices Qt  0 and Rt  0, respectively, the EKF can
be implemented with the following recursions


v̂t+1|t+1 = v̂t+1|t + Kt+1 zt+1 − h(v̂t+1|t )
(74)
where the state predictor v̂t+1|t and the Kalman gain Kt+1 are given by
v̂t+1|t = Ft v̂t|t + gt

(75a)

H
Kt+1 = Pt+1|t JH
t+1 Jt+1 Pt+1|t Jt+1 + Rt+1

Pt+1|t+1 = Pt+1|t − Kt+1 Jt+1 Pt+1|t
Pt+1|t = Ft Pt|t FH
t + Qt
31

−1

(75b)
(75c)
(75d)

with Jt+1 being the measurement Jacobian matrix of h evaluated at v̂t+1|t ,
and Pt+1|t+1  0 (Pt+1|t  0) denoting the corrected (predicted) state
estimation error covariance matrix at time t + 1. To improve on the approximation accuracy of the EKF, extended Kalman filters (UKF) have been
reported in [65]; see also [79] for their robust versions. Particle filtering
may also be useful if its computational complexity can be supported during
real-time power systems operations.
Because the EKF and UKF are known to diverge for highly nonlinear
dynamics, moving horizon estimation (MHE) has been suggested as an accurate yet tractable alternative with proven robustness to bounded model
errors [59]. Different from Kalman filtering, the initial state v0 , and noises
ωt and t in MHE are viewed as deterministic unknowns taking values from
given bounded sets S, W, and E, respectively. The sets W and E model
disturbances with truncated densities [59].
The idea behind MHE is to perform PSSE by exploiting useful information present in a sliding window of the most recent observations. Consider
here a sliding window of length L + 1. Let v̂t−L|t denote the smoothed
estimate at time t − L given L past measurements, as well as the current one, namely {zτ }tτ =t−L . MHE aims at obtaining the most recent L
t
state estimates {v̂t−L+s|t }L
s=0 based on {zτ }t−L and the available estimate
v̌t−L := v̂t−L|t−1 from time t − 1 and for t ≥ L. A key simplification is that
once v̂t−L|t becomes available, the other L recent estimates at time t can be
recursively obtained through ‘noise-free’ propagation based on the dynamic
model (73a); that is,
v̂t−L+s|t = Ft−L+s−1 v̂t−L+s−1|t

(76)

for s = 1, . . . , L. By relating all recent estimates to v̂t−L|t via successive
multiplications of transition matrices, the update in (76) simplifies to
v̂t−L+s|t = Tt−L+s v̂t−L|t

(77)

where Tt−L+s := Ft−L+s−1 Tt−L+s−1 for s = 1, . . . , L, with Tt−L = I. The
MHE-based state estimate v̂t−L|t is then given by
v̂t−L|t := arg min
v

L
X

zt−L+s − h(Tt−L+s v)

2
2

+ λkv − v̌t−L k22

(78)

s=0

where λ > 0 can be tuned relying on our confidence in the state predictor
v̌t−L , and the measurements {zτ }tt−L . Given the quadratic dependence of
32

the SCADA measurements {h(vt )} and the state v, the optimization problem in (78) is non-convex.
Finding the MHE-based state estimates in real time entails online solutions of dynamic optimization problems. The MHE formulation can be
convexified by exploiting the semidefinite relaxation: vector v is lifted to the
matrix V := vvH  0, and the m-th entry of h(vt−L+s ) for s = 0, . . . , L, is
expressed as
H
hm (Tt−L+s v) = vH TH
t−L+s Hm Tt−L+s v = Tr(Tt−L+s Hm Tt−L+s V).

Upon dropping the nonconvex rank constraint rank(V) = 1, the SDP-based
MHE yields
V̂t−L|t := arg min

V0

L
X

zt−L+s −Tr TH
t−L+s Hm Tt−L+s V



2
+λkv− v̌t−L k22
2

s=0

which can be solved in polynomial time using off-the-shelf toolboxes. Rankone state estimates can be obtained again through eigen-decomposition or
randomization. The complexity of solving the last problem is rather high
in its present form on the order of Nb4.5 [46]. Therefore, developing faster
solvers for the SDP-based MHE by exploiting the rich sparsity structure in
{Hm } matrices is worth investigating. Decentralized and localized MHE
implementations are also timely and pertinent. Devising FPP-based solvers
for the MHE in (78) constitutes another research direction.

7

Discussion

This chapter has reviewed some of the recent advances in PSSE. After developing the CRLB, an SDP-based solver, and its regularized counterpart
were discussed. To overcome the high complexity involved, a scheme named
feasible point pursuit relying on successive convex approximations was also
advocated. A decentralized PSSE paradigm put forth provides the means
for coping with the computationally-intensive SDP formulations, it is tailored for the interconnected nature of modern grids, while it can also afford
processing PMU data in a timely fashion. A better understanding of cyber
attacks and disciplined ways for decentralized bad data processing were also
provided. Finally, this chapter gave a fresh perspective to state tracking
under model-free and model-based scenarios.
Nonetheless, there are still many technically challenging and practically
pertinent grid monitoring issues to be addressed. Solving power grid data
33

processing tasks on the cloud has been a major trend to alleviate data storage, communication, and interoperability costs for system operators and
utilities. Moreover, with the current focus on low- and medium-voltage distribution grids, solvers for unbalanced and multi-phase operating conditions
[2] are desirable. Smart meters and synchrophasor data from distribution
grids (also known as micro-PMUs [67]) call for new data processing solutions. Advances in machine learning and statistical signal processing, such
as sparse and low-rank models, missing and incomplete data, tensor decompositions, deep learning, nonconvex and stochastic optimization tools, and
(multi)kernel-based learning to name a few, are currently providing novel
paths to grid monitoring tasks while realizing the vision of smarter energy
systems.

Acknowledgements
G. Wang and G. B. Giannakis were supported in part by NSF grants 1423316,
1442686, 1508993, and 1509040; and by the Laboratory Directed Research
and Development Program at the NREL. H. Zhu was supported in part by
NSF grants 1610732 and 1653706.

Appendix
2 })).
Proof of Prop. 1. Consider the AGWN model (15) with  ∼ N (0, dg({σm
The data likelihood function is


M
Y
1
(zm − vH Hm v)2
p
p(z; v) =
(79)
exp −
2
2
2σm
2πσm
m=1

and the negative log-likelihood function denoted by f (v) = − ln p(z; x) is
f (v) =


M 
X
2 1

1
H
2
z
−
v
H
v
+
ln
2πσ
.
m
m
m
2
2σm
2

(80)

m=1

The Fisher information matrix (FIM) is defined as the Hessian of the
real-valued function f (v) with respect to the complex vector v ∈ CNb . Deriving the CRLB amounts to finding the Hessian of a real-valued function
with respect to a complex-valued vector. Wirtinger’s calculus confirms that
f (v) can be equivalently rewritten as f (v, v∗ ); see e.g., [41]. Upon introducing the conjugate coordinates [vT (v∗ )T ]T ∈ C2Nb , the Wirtinger derivatives, namely the first-order partial differential operators of functions over
34

complex domains, are given by [41]
∂f
∂f (v, v∗ )
:=
∂v
∂vT
∂f
∂f (v, v∗ )
:=
∂v∗
∂(v∗ )T


∂f
∂f
···
=
∂V1
∂VN constant v∗
v∗


∂f
∂f
.
=
···
∗
∂V1∗
∂VN
v
constant v


constant

constant

These definitions follow the convention in multivariate calculus that derivatives are denoted by row vectors, and gradients by column vectors. Define
for notational brevity φm (v, v∗ ) := zm − (v∗ )T Hm v for m = 1, . . . , M .
Accordingly, the Wirtinger derivatives of f (v, v∗ ) in (80) are obtained as
M
L
X
X
∂f
1
∂φm
∂f
1
∂φm
=
φ
and
=
φ
2 m ∂vT
2 m ∂(v∗ )T
∂v
σm
∂v∗
σm
m=1

(82)

m=1

and the Wirtinger derivatives of φm (v, v∗ ) can be found likewise
∂φm
∂φm
= −(Hm v)H and
= −(H∗m v∗ )H .
∂vT
∂(v∗ )T

(83)

In the conjugate coordinate system, the complex Hessian of f (v, v∗ ) with
respect to the conjugate coordinates [vT (v∗ )T ]T is defined as


Hvv Hv∗ v
∗
2
∗
H(v, v ) := ∇ f (v, v ) =
(84)
Hvv∗ Hv∗ v∗
whose blocks are given as


Hvv

∂
:=
∂vT



Hvv∗

∂
:=
∂vT

∂f
∂v

H

∂f
∂v∗

,



Hv∗ v

∂f
∂v



Hv ∗ v ∗

∂
:=
∂(v∗ )T

∂f
∂v∗

H
,

H

∂
:=
∂(v∗ )T

H
.

Substituting (82) and (83) into the last equations and after algebraic manipulations yields
Hvv =
Hv ∗ v =

M
X
m=1
M
X



−2
σm
Hm v(Hm v)H − φm Hm

(85a)

−2
σm
Hm v(H∗m v∗ )H

(85b)

m=1

35

Hvv∗ =
Hv ∗ v ∗ =

M
X
m=1
M
X

−2 ∗ ∗
σm
Hm v (Hm v)H

(85c)



−2
σm
H∗m v∗ (H∗m v∗ )H − φm H∗m .

(85d)

m=1

Evaluating the Hessian blocks of (85) at the true value of v, and taking
the expectation with respect to , yields E[φm ] = 0. Hence, the φm -related
terms in (85) disappear, and the FIM F(v, v∗ ) := E[H(v, v∗ )] simplifies to
the expression in (23); see also [66].
To show that the FIM is rank-deficient,
define gm := [(Hm v)H (H∗m v∗ )H ]H ,
PM
−2 g gH . Observe now that the nonso that the FIM becomes F = m=1 σm
m m
T
∗
T
T
zero vector d(v) := [v − (v ) ] is orthogonal to gm for m = 1, . . . , M ;
that is,
H
gm
d = vH Hm v − (vH Hm v)∗ = 0.
Based on the latter, it is not hard to verify that Fd = 0, which proves that
the null space of F is non-empty.

References
[1] A. Abur and A. Gómez-Expósito, Power System State Estimation: Theory and Implementation. New York, NY: Marcel Dekker, 2004.
[2] M. Bazrafshan and N. Gatsis, “Comprehensive modeling of three-phase
distribution systems via the bus admittance matrix,” arXiv:1705.06782,
2017.
[3] D. P. Bertsekas, Nonlinear Programming, 2nd ed.
Athena Scientific, 1999.

Belmont, MA:

[4] 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.
[5] S. Boyd and L. Vandenberghe, Convex Optimization. New York, NY:
Cambridge University Press, 2004.
[6] M. K. Celik and A. Abur, “A robust WLAV state estimator using transformations,” IEEE Trans. Power Syst., vol. 7, no. 1, pp. 106–113, Feb.
1992.
36

[7] S. S. Chen, D. L. Donoho, Michael, and A. Saunders, “Atomic decomposition by basis pursuit,” SIAM J. Sci. Comput., vol. 20, pp. 33–61,
Jul. 1998.
[8] K. A. Clements, G. R. Krumpholz, and P. W. Davis, “Power system state estimation with measurement deficiency: An observability/measurement placement algorithm,” IEEE Trans. Power App.
Syst., vol. 102, no. 7, pp. 2012–2020, Jul. 1983.
[9] A. J. Conejo, S. de la Torre, and M. Canas, “An optimization approach
to multiarea state estimation,” IEEE Trans. Power Syst., vol. 22, no. 1,
pp. 213–221, Feb. 2007.
[10] T. V. Cutsem and M. Ribbens-Pavella, “Critical survey of hierarchical
methods for state estimation of electric power systems,” IEEE Trans.
Power App. Syst., vol. 102, no. 10, pp. 3415–3424, Oct. 1983.
[11] A. L. Da Silva, M. Do Coutto Filho, and J. De Queiroz, “State forecasting in electric power systems,” IEE Proc. Generat. Transm. Distrib.,
vol. 130, no. 5, pp. 237–244, Sep. 1983.
[12] E. Dall’Anese, H. Zhu, and G. B. Giannakis, “Distributed optimal
power flow for smart microgrids,” IEEE Trans. Smart Grid, vol. 4,
no. 3, pp. 1464–1475, Sep. 2013.
[13] J. De La Ree, V. A. Centeno, J. Thorp, and A. Phadke, “Synchronized phasor measurement applications in power systems,” IEEE Trans.
Smart Grid, vol. 1, no. 1, pp. 20–27, Jun. 2010.
[14] A. S. Debs and R. E. Larson, “A dynamic estimator for tracking the
state of a power system,” IEEE Trans. Power App. Syst., no. 7, pp.
1670–1678, Sep. 1970.
[15] D. Duan, L. Yang, and L. L. Scharf, “Phasor state estimation from
PMU measurements with bad data,” in Proc. IEEE Workshop on
Comp. Adv. in Multi-Sensor Adaptive Proc., San Juan, Puerto Rico,
Dec. 2011.
[16] R. Ebrahimian and R. Baldick, “State estimation distributed processing,” IEEE Trans. Power Syst., vol. 15, no. 4, pp. 1240–1246, Nov.
2000.

37

[17] A. A. El-Keib, J. Nieplocha, H. Singh, and D. Maratukulam, “A decomposed state estimation technique suitable for parallel processor implementation,” IEEE Trans. Power Syst., vol. 7, no. 3, pp. 1088–1097,
Aug. 1992.
[18] D. M. Falcao, F. F. Wu, and L. Murphy, “Parallel and distributed state
estimation,” IEEE Trans. Power Syst., vol. 10, no. 2, pp. 724–730, May
1995.
[19] G. B. Giannakis, V. Kekatos, N. Gatsis, S.-J. Kim, H. Zhu, and B. Wollenberg, “Monitoring and optimization for power grids: A signal processing perspective,” IEEE Signal Process. Mag., vol. 30, no. 5, pp.
107–128, Sept. 2013.
[20] G. B. Giannakis, Q. Ling, G. Mateos, I. D. Schizas, and H. Zhu, “Decentralized learning for wireless communications and networking,” Splitting Methods in Communication, Imaging, Science, and Engineering,
pp. 461–497, Springer, 2016.
[21] M. X. Goemans and D. P. Williamson, “Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite
programming,” J. ACM, vol. 42, no. 6, pp. 1115–1145, Nov. 1995.
[22] A. Gómez-Expósito, A. Abur, A. de la Villa Jaén, and C. Gómez-Quiles,
“A multilevel state estimation paradigm for smart grids,” Proc. IEEE,
vol. 99, no. 6, pp. 952–976, Jun. 2011.
[23] A. Gómez-Expósito, A. de la Villa Jaén, C. Gómez-Quiles,
P. Rousseaux, and T. V. Cutsem, “A taxonomy of multi-area state
estimation methods,” Electr. Pow. Syst. Res., vol. 81, pp. 1060–1069,
Apr. 2011.
[24] R. Grone, C. R. Johnson, E. M. Sá, and H. Wolkowicz, “Positive definite
completions of partial hermitian matrices,” Linear Algebra & its Appl.,
vol. 58, pp. 109–124, Apr. 1984.
[25] S. Haykin, Adaptive Filter Theory. New York, NY: Prentice Hall, 2002.
[26] S.-J. Huang and K.-R. Shih, “Dynamic-state-estimation scheme including nonlinear measurement function considerations,” IEE Proc. Generat. Transm. Distrib., vol. 149, no. 6, pp. 673–678, Nov. 2002.

38

[27] S. Iwamoto, M. Kusano, and V. H. Quintana, “Hierarchical state estimation using a fast rectangular-coordinate method,” IEEE Trans.
Power Syst., vol. 4, no. 3, pp. 870–880, Aug. 1989.
[28] R. A. Jabr, “Exploiting sparsity in SDP relaxations of the OPF problem,” IEEE Trans. Power Syst., vol. 27, no. 2, pp. 1138–1139, May
2012.
[29] S. M. Kay, Fundamentals of Statistical Signal Processing, Vol. I: Estimation Theory. Prentice Hall, 1993.
[30] V. Kekatos and G. B. Giannakis, “From sparse signals to sparse residuals for robust sensing,” IEEE Trans. Signal Process., vol. 59, no. 7,
pp. 3355–3368, Jul. 2011.
[31] ——, “Distributed robust power system state estimation,” IEEE Trans.
Power Syst., vol. 28, no. 2, pp. 1617–1626, May 2013.
[32] V. Kekatos, G. B. Giannakis, and B. F. Wollenberg, “Optimal placement of phasor measurement units via convex relaxation,” IEEE Trans.
Power Syst., vol. 27, no. 3, pp. 1521–1530, Aug. 2012.
[33] V. Kekatos, E. Vlachos, D. Ampeliotis, G. B. Giannakis, and
K. Berberidis, “A decentralized approach to generalized power system state estimation,” in Proc. IEEE Wkshp. on Comp. Adv. in MultiSensor Adaptive Proc., Saint Martin, Dec. 2013.
[34] V. Kekatos, G. Wang, A. J. Conejo, and G. B. Giannakis, “Stochastic reactive power management in microgrids with renewables,” IEEE
Trans. Power Syst., vol. 30, no. 6, pp. 3386–3395, Nov. 2015.
[35] S. J. Kim, G. B. Giannakis, and K. Y. Lee, “Online optimal power flow
with renewables,” in Asilomar Conf. on Signals, Syst. and Comput.,
Pacific Grove, CA, Nov 2014, pp. 355–360.
[36] S.-J. Kim, “Online power system state estimation using alternating
direction method of multipliers,” in Proc. IEEE PES Summer Meeting,
Denver, CO, 2015.
[37] S.-J. Kim, G. Wang, and G. B. Giannakis, “Online semidefinite programming for power system state estimation,” in Proc. IEEE Intl. Conf.
on Acoustics, Speech, and Signal Process., Florence, Italy, May 2014,
pp. 6024–6027.
39

[38] G. N. Korres, “A distributed multiarea state estimation,” IEEE Trans.
Power Syst., vol. 26, no. 1, pp. 73–84, Feb. 2011.
[39] O. Kosut, L. Jia, J. Thomas, and L. Tong, “Malicious data attacks on
the smart grid,” IEEE Trans. Smart Grid, vol. 2, no. 4, pp. 645–658,
Dec. 2011.
[40] 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.
[41] K. Kreutz-Delgado, “The complex gradient operator and the CRcalculus,” arXiv:0906.4835, 2009.
[42] A. Y. Lam, B. Zhang, and T. David, “Distributed algorithms for optimal power flow problem,” in Proc. IEEE Conf. on Decision and Control,
Maui, HI, Dec. 2012, pp. 430–437.
[43] S.-Y. Lin and C.-H. Lin, “An implementable distributed state estimator and distributed bad data processing schemes for electric power systems,” IEEE Trans. Power Syst., vol. 9, no. 3, pp. 1277–1284, Aug.
1994.
[44] T. Lipp and S. Boyd, “Variations and extension of the convex–concave
procedure,” Optim. Eng., vol. 17, no. 2, pp. 263–287, Jun. 2016.
[45] Y. Liu, P. Ning, and M. K. Reiter, “False data injection attacks against
state estimation in electric power grids,” ACM Trans. Info. and System
Security, vol. 14, no. 1, pp. 13:1–13:33, May 2011.
[46] Z.-Q. Luo, W.-K. Ma, A. M.-C. So, Y. Ye, and S. Zhang, “Semidefinite
relaxation of quadratic optimization problems,” IEEE Signal Process.
Mag., vol. 27, no. 3, pp. 20–34, May 2010.
[47] R. Madani, A. Ashraphijuo, J. Lavaei, and R. Baldick, “Power system
state estimation with a limited number of measurements,” in Proc.
IEEE Conf. on Decision and Control, Las Vegas, NV, Dec. 2016.
[48] R. Madani, J. Lavaei, and R. Baldick, “Convexification of power flow
problem over arbitrary networks,” in Proc. IEEE Conf. on Decision
and Control, Osaka, Japan, Dec. 2015.
[49] O. L. Mangasarian and D. R. Musicant, “Robust linear and support
vector regression,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 22,
no. 9, pp. 950–955, Sep. 2000.
40

[50] R. A. Maronna, R. D. Martin, and V. J. Yohai, Robust Statistics: Theory and Methods. Wiley, 2006.
[51] O. Mehanna, K. Huang, B. Gopalakrishnan, A. Konar, and
N. Sidiropoulos, “Feasible point pursuit and successive approximation
of non-convex QCQPs,” vol. 22, no. 7, pp. 804–808, Nov. 2015.
[52] 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.
[53] A. Minot, Y. M. Lu, and N. Li, “A distributed Gauss-Newton method
for power system state estimation,” IEEE Trans. Power Syst., vol. 31,
no. 5, pp. 3804–3815, Sep. 2016.
[54] B. Y. Mo, T. H.-J. Kim, K. Brancik, D. Dickinson, H. Lee, A. Perrig,
and B. Sinopoli, “Cyber-physical security of a smart grid infrastructure,” Proc. IEEE, vol. 100, no. 1, pp. 195–209, Jan. 2012.
[55] A. Monticelli, “Electric power system state estimation,” Proc. IEEE,
vol. 88, no. 2, pp. 262–282, Feb. 2000.
[56] A. Monticelli and F. F. Wu, “Network observability: Theory,” IEEE
Trans. Power App. Syst., vol. 104, no. 5, pp. 1042–1048, May 1985.
[57] P. M. Pardalos and S. A. Vavasis, “Quadratic programming with one
negative eigenvalue is NP-hard,” J. Global Optim., vol. 1, no. 1, pp.
15–22, Mar. 1991.
[58] A. G. Phadke and J. S. Thorp, Synchronized Phasor Measurements and
Their Applications. New York, NY: Springer, 2008.
[59] C. V. Rao, J. B. Rawlings, and D. Q. Mayne, “Constrained state estimation for nonlinear discrete-time systems: Stability and moving horizon approximations,” IEEE Trans. Autom. Contr., vol. 48, no. 2, pp.
246–258, Feb. 2003.
[60] I. D. Schizas, A. Ribeiro, and G. B. Giannakis, “Consensus in ad hoc
WSNs with noisy links: Part I: Distributed estimation of deterministic
signals,” IEEE Trans. Signal Process., vol. 56, no. 1, pp. 350–364, Jan.
2008.
[61] 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.
41

[62] S. Shalev-Shwartz, “Online learning and online convex optimization,”
Found. and Trends Mach. Learning, vol. 4, no. 2, pp. 107–194, Mar.
2011.
[63] P. Stoica and T. L. Marzetta, “Parameter estimation problems with
singular information matrices,” IEEE Trans. Signal Process., vol. 49,
no. 1, pp. 87–90, Jan. 2001.
[64] Power systems test case archive. Univ. of Washington. [Online].
Available: http://www.ee.washington.edu/research/pstca/
[65] G. Valverde and V. Terzija, “Unscented Kalman filter for power system
dynamic state estimation,” IEE Proc. Generat. Transm. Distrib., vol. 5,
no. 1, pp. 29–37, Jan. 2011.
[66] A. Van den Bos, “A Cramér-Rao lower bound for complex parameters,”
IEEE Trans. Signal Process., vol. 40, no. 10, Oct. 1994.
[67] A. Von Meier, D. Culler, A. McEachern, and R. Arghandeh, “Microsynchrophasors for distribution systems,” in Proc. IEEE Conf. on Innovative Smart Grid Technologies, 2014, pp. 1–5.
[68] G. Wang, G. B. Giannakis, and J. Chen, “Robust and scalable power
system state estimation via composite optimization,” preprint, 2017.
[69] G. Wang, V. Kekatos, A. J. Conejo, and G. B. Giannakis, “Ergodic energy management leveraging resource variability in distribution grids,”
IEEE Trans. Power Syst., vol. 31, no. 6, pp. 4765–4775, Nov. 2016.
[70] G. Wang, A. S. Zamzam, G. B. Giannakis, and N. D. Sidiropoulos,
“Power system state estimation via feasible point pursuit,” in Proc.
IEEE Global Conf. on Signal and Inf. Process., Washington, D.C., Nov.
2016.
[71] ——, “Power system state estimation via feasible point pursuit: Algorithms and Crmér-Rao bound,” arXiv:1203.0853, 2017.
[72] 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, Jul. 2014, pp. 1–5.
[73] S. Wang, W. Gao, and A. S. Meliopoulos, “An alternative method for
power system dynamic state estimation based on unscented transform,”
IEEE Trans. Power Syst., vol. 27, no. 2, pp. 942–950, May 2012.
42

[74] A. J. Wood and B. F. Wollenberg, Power Generation, Operation, and
Control, 2nd ed. New York, NY: Wiley & Sons, 1996.
[75] L. Xie, D.-H. Choi, and S. Kar, “Cooperative distributed state estimation: Local observability relaxed,” in Proc. IEEE Power & Energy
Society General Meeting, Detroit, MI, Jul. 2011.
[76] W. Xu, M. Wang, J. F. Cai, and A. Tang, “Sparse error correction
from nonlinear measurements with applications in bad data detection
for power networks,” IEEE Trans. Signal Process., vol. 61, no. 24, pp.
6175–6187, Dec. 2013.
[77] A. S. Zamzam, N. D. Sidiropoulos, and E. Dall’Anese, “Beyond relaxation and Newton-Raphson: Solving AC OPF for multi-phase systems
with renewables,” IEEE Trans. Smart Grid, 2017.
[78] 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).
[79] J. Zhao, M. Netto, and L. Mili, “A robust iterated extended Kalman
filter for power system dynamic state estimation,” IEEE Trans. Power
Syst., 2016 (to appear).
[80] L. Zhao and A. Abur, “Multiarea state estimation using synchronized
phasor measurements,” IEEE Trans. Power Syst., vol. 20, no. 2, pp.
611–617, May 2005.
[81] H. Zhu and G. B. Giannakis, “Robust power system state estimation
for the nonlinear AC flow model,” in Proc. North American Power
Symposium, Urbana-Champaign, IL, Sep. 2012.
[82] H. Zhu, G. B. Giannakis, and A. Cano, “Distributed in-network decoding,” IEEE Trans. Signal Process., vol. 57, no. 10, pp. 3970–3983, Oct.
2009.
[83] H. Zhu and G. B. Giannakis, “Estimating the state of AC power systems using semidefinite programming,” in Proc. North American Power
Symposium, Boston, MA, Aug. 2011.
[84] ——, “Power system nonlinear state estimation using distributed
semidefinite programming,” IEEE J. Sel. Topics Signal Process., vol. 8,
no. 6, pp. 1039–1050, Dec. 2014.

43

