IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. , NO. , 2017

1

Robust Power System Dynamic State Estimator with
Non-Gaussian Measurement Noise: Part I–Theory

arXiv:1703.04790v1 [] 14 Mar 2017

Junbo Zhao, Student Member, IEEE, Lamine Mili, Fellow, IEEE

Abstract—This paper develops the theoretical framework and
the equations of a new robust Generalized Maximum-likelihoodtype Unscented Kalman Filter (GM-UKF) that is able to suppress
observation and innovation outliers while filtering out nonGaussian measurement noise. Because the errors of the real and
reactive power measurements calculated using Phasor Measurement Units (PMUs) follow long-tailed probability distributions,
the conventional UKF provides strongly biased state estimates
since it relies on the weighted least squares estimator. By
contrast, the state estimates and residuals of our GM-UKF are
proved to be roughly Gaussian, allowing the sigma points to
reliably approximate the mean and the covariance matrices of
the predicted and corrected state vectors. To develop our GMUKF, we first derive a batch-mode regression form by processing
the predictions and observations simultaneously, where the statistical linearization approach is used. We show that the set of
equations so derived are equivalent to those of the unscented
transformation. Then, a robust GM-estimator that minimizes a
convex Huber cost function while using weights calculated via
Projection Statistics (PS’s) is proposed. The PS’s are applied
to a two-dimensional matrix that consists of serially correlated
predicted state and innovation vectors to detect observation and
innovation outliers. These outliers are suppressed by the GMestimator using the iteratively reweighted least squares algorithm.
Finally, the asymptotic error covariance matrix of the GM-UKF
state estimates is derived from the total influence function. In the
companion paper, extensive simulation results will be shown to
verify the effectiveness and robustness of the proposed method.
Index Terms—Dynamic state estimation, robust estimation,
unscented Kalman filter, non-Gaussian noise, total influence
function, outliers, cyber attacks, power system dynamics.

I. I NTRODUCTION
HE widespread deployment of synchro-Phasor Measurement Units (PMUs) on power transmission grids has
made possible the real-time monitoring and control of power
system dynamics. However, these functions cannot be reliably achieved without the development of a fast and robust
Dynamic State Estimator (DSE). Indeed, the state variable
estimates of the synchronous machines can be utilized by
power system stabilizers, automatic voltage regulators, and
under-frequency relays to enhance small signal stability and to
initiate generation outages and load shedding during transient
instabilities, among other actions [1], [2].
To date, a variety of dynamic state estimators have been
proposed in the literature; they are based on the Extended
Kalman Filter (EKF) [3], [4], the Iterated EKF (IEKF) [5],
[6], the unscented Kalman filter (UKF) [7]–[9], to cite a
few. However, all these methods suffer from several important shortcomings, precluding them from being adopted by
power utilities for power system real-time applications. To

T

Junbo Zhao and Lamine Mili are with the Bradley Department of Electrical Computer Engineering, Virginia Polytechnic Institute and State University, Northern Virginia Center, Falls Church, VA 22043, USA (e-mail:
zjunbo@vt.edu, lmili@vt.edu).

be specific, they cannot handle i) non-Gaussian process and
observation noise of the system nonlinear dynamic models and
ii) innovation, observation and structural outliers.
There are several reasons for these shortcomings. Firstly,
the current DSE approaches assume that both the process and
the observation noise of the system nonlinear dynamic models
are Gaussian. However, two recent investigations conducted by
PNNL [10], [11] revealed that the PMU measurement errors of
the voltage and current magnitudes obey non-Gaussian probability distributions. This is demonstrated in Fig.1 using real
PMU data provided to us by PNNL. This figure displays histograms and parametric probability density estimates of PMU
errors on nodal voltage magnitudes and angles, line current
magnitudes and angles, and line real and reactive powers. As
observed in Fig.1, except for the measurement errors on nodal
voltage and line current angles, which are roughly Gaussian,
the measurement errors on both nodal voltage and line current
magnitudes obey a bimodal Gaussian mixture distribution. As
for the measurement errors of line real and reactive powers
calculated from voltage and current phasors, they follow a
thick tailed distribution that may be approximated by either the
Laplacian or the Cauchy distribution. Recall that in contrast to
the Gaussian distribution, which is a short-tailed distribution,
a thick-tailed distribution is the one that allows the associated
random variable to take, as compared to a scale parameter,
large values with a non-negligible probability. Evidently, the
presence of non-Gaussian noise calls for new research and
development in robust power system DSE based on robust
statistics.
Secondly, three types of outliers associated with a given
dynamical system model have been defined by Gandhi and
Mili [12], namely observation outliers, which affect the metered values; innovation outliers, which corrupt the predicted
state estimates; and structural outliers, which affect the system
dynamic states and the observation functions. Observation outliers may result from large biases in PMU measurements due
to infrequent calibration, or instrument failures, or impulsive
communication noise [13], [14]. As for innovation outliers,
they may occur in several different ways. For example, some
of the generator models may not be well calibrated, resulting
in highly inaccurate model outputs that are inconsistent with
the measurements. This was precisely the case in the 1996
blackout, where the model being used predicted system stability while in reality the system was undergoing numerous
cascading failures, which resulted in a rapid system collapse
that occurred within minutes [15], [16]. Innovation outliers
may also be induced by the approximations in the state
prediction model or by a system process impulsive noise.
By contrast, structural outliers are induced by wrong circuit

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. , NO. , 2017

2

!"#$%$&'&()*Q,-.&()*#J*(6,*P6%.#"*L#'(%3,*4%3-&(0Q,*,""#"

!"#$%$&'&()*Q,-.&()*#J*(6,*P6%.#"*10"",-(*4%3-&(0Q,*,""#"

<

IC
<A

?@A

!"#$%$&'&()*Q,-.&()*L%'0,

!"#$%$&'&()*Q,-.&()*L%'0,

<C

?

B@A

B

?A
?C
BA
BC

C@A

A

C
DC@H

DC@I

DC@?

C

C@?

C@I

C@H

C
DC@CA

C@F

DC@CI

DC@C<

DC@C?

DC@CB

C

C@CB

C@C?

C@C<

C@CI

>

>
!"#$%$&'&()*+,-.&()*/0-1(&#-

!"#$%$&'&()*+,-.&()*/0-1(&#-

C@BB
C@CFF

!""#" $%&'"%()'%#*
#+ ',- .,/&#"
0#1'/2- /*21-

C@B
C@CE
C@CF

C@CF
C@CG?

!""#" $%&'"%()'%#*
#+ ',- .,/&#"
3)""-*' /*21-

C@CHI

C@CG
C@CAH

J;>=

C@CH
C@CIF

C@CA
C@CI

C@CI
C@C<?

C@C<
C@C?I

C@C?
C@CBH

C@CB
C@CCF

C
D?@A

D?

DB@A

DB

DC@A

C

C@A

B

B@A

?

?@A

C

<

D?@A

>
2&.(#3"%4

5#6-.#-*78

D?

DB@A

DB

DC@A

C

C@A

B

B@A

?

?@A

>

9#"4%'

:#3-#"4%'*;<!=

2&.(#3"%4

5#6-.#-*78

9#"4%'

!"#$%$&'&()*+,-.&()*/0-1(&#-*#J*K,%1(&L,*!#M,"*N""#"

!"#$%$&'&()*+,-.&()*/0-1(&#-*#J*K,%'*!#M,"*N""#"
??
C@?

C@?

4/)3,5 $%&'"%()'%#*

C@BF

BF

4/)3,5 $%&'"%()'%#*
C@BH

BH

6/.1/3- $%&'"%()'%#*

BI

C@BI

B?

C@B?

J;>=

6/.1/3- $%&'"%()'%#*

C@B

C@B

CF

C@CF

CH

C@CH

CI

C@CI

C?

C@C?

C

C
DC@CCH

DC@CCA

DC@CCI

DC@CC<

DC@CC?

DC@CCB

C

C@CCB

C@CC?

C@CC<

C@CCI

C@CCA

C@CCH

DC@CCA

DC@CCI

DC@CC<

DC@CC?

DC@CCB

2&.(#3"%4

O%016)

C

C@CCB

C@CC?

C@CC<

C@CCI

>

>
:%P'%1,

2&.(#3"%4

O%016)

:%P'%1,

Fig. 1. Error distributions of the phasor voltage angle and magnitude, the phasor current angle and magnitude, the real and reactive power using field PMU
data.

breaker statuses or gross errors in the model parameters of
the transmission lines, or of the automatic voltage regulators,
or of the synchronous machines. In [17], it is reported that
wrong estimates of the parameters of the synchronous machine
models may result from the use of erroneous metered values.
It turns out that the conventional filters, namely the EKF, the
IEKF, the UKF, and the Particle Filter (PF) are not robust
to any type of outliers. For instance, it is demonstrated in
[18] that their performances are significantly degraded in
the presence of observation outliers. To address this issue,
Rouhani and Abur [9] developed a robust UKF-based DSE
using the Least-Absolute-Value (LAV) estimator. However,
the authors do not address the vulnerability of the DSE to
innovation outliers. In [6], a robust IEKF was proposed to
handle observation and innovation outliers, but it may suffer

from divergence problems if the nonlinearity of the system
model is strong. In addition, both [6], [9] do not address the
non-Gaussianity of the measurement noise.
In this paper, a robust Generalized Maximum-Likelihoodtype UKF (GM-UKF) method is proposed to suppress observation and innovation outliers while filtering out non-Gaussian
measurement noise. Our choice of the UKF is motivated
by the fact that, considering the real-time implementation
requirements for nonlinear DSE, it achieves a more balanced
performance between computational efficiency and ability to
cope with strong system nonlinearities than the EKF, or the
IEKF, or the PF. However, the UKF is based on the sigma
points, which reliably approximate the mean and the covariance matrices of the state estimates only under the Gaussian
assumption of the process and observation noises. We show

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. , NO. , 2017

3

that this assumption is further stressed by the reliance of the
UKF on the weighted least squares estimator. Interestingly,
the state estimates calculated by our GM-UKF are shown
to be asymptotically Gaussian even when the noises obey
thick-tailed distributions, which is precisely the case when
using PMU measurements. Furthermore, we show that the
state estimates obtained from the application of statistical
linearization to the nonlinear discrete-time state space system
model are equivalent to those of the unscented transformation.
Therefore, our filter allows the sigma points to provide good
results.
It is developed according to the following steps. We first
derive a redundancy batch-mode regression form by processing
the predictions and observations simultaneously; this overdetermined system of equations provides the data redundancy
needed for the detection and suppression of the innovation
and observation outliers. This is achieved by means of a
robust GM-estimator defined as the minimum of the Huber
convex cost function while using weights calculated via the
Projection Statistics (PS’s). The latter are applied to a twodimensional matrix consisting of serially correlated predicted
state and innovation vectors. Then, a statistical test is applied
to them to flag the outliers. Finally, the GM-estimator is solved
via the iteratively reweighted least squares algorithm and the
asymptotic error covariance matrix of the state estimates is
calculated from the total influence function.
The rest of the paper is organized as follows. Section II
presents the problem formulation. Section III develops the
theory of the proposed GM-UKF and finally Section IV
concludes the paper.
II. P ROBLEM F ORMULATION
A. Nonlinear Discrete-Time Dynamical System Model
A discrete-time state space representation of a general
nonlinear dynamical system is expressed as
xk = f (xk−1 , uk ) + wk ,

(1)

zk = h (xk , uk ) + vk ,

(2)

n×1

matrices is required, which can be by itself a difficult task to
achieve in some cases or computationally costly.
To be specific, given a state estimate at time step k-1,
bk−1|k−1 ∈ Rn×1 , having a covariance matrix given by
x
xx
Pk−1|k−1
, its statistics are captured by 2n weighted sigma
points defined as

q
xx
bk−1|k−1 ±
,
(3)
nPk−1|k−1
χik−1|k−1 = x
i

1
with weights wi = 2n
, i = 1, ..., 2n. Then, each sigma point
is propagated through the nonlinear system process model (1),
yielding a set of transformed samples expressed as


(4)
χik|k−1 = f χik−1|k−1 .

Next, the predicted sample mean and sample covariance matrix
of the state vector are calculated by

xx
=
Pk|k−1

2n
X
i=1

B. Dynamic State Estimation using UKF
The main idea underlying the UKF is the application of
a deterministic sampling technique known as the unscented
transformation, which allows us, under the Gaussian noise
assumption, to choose a set of sample points, termed sigma
points, that have the same mean and covariance matrix as those
of the a priori state vector [19]. These sigma points are then
propagated through the non-linear functions f and h, yielding
an estimation of the a posteriori state statistics by using the
Kalman filter approach, i.e., the sample mean and the sample
covariance matrix. Consequently, no calculation of Jacobian

2n
X

wi χik|k−1 ,

(5)

i=1

bk|k−1 )(χik|k−1 − x
bk|k−1 )T +Qk .
wi (χik|k−1 − x

(6)
Finally, the measurement updating is performed and the filxx
bk|k with the covariance matrix Pk|k
tered state x
are calculated
by


zz
xz
Pk|k−1
Kk = Pk|k−1

−1

,


bk|k = x
b k|k−1 + Kk zk − zbk|k−1 ,
x
zz
xx
xx
KkT ,
− Kk Pk|k−1
= Pk|k−1
Pk|k

where zbk|k−1 =

2n
P

i=1

(7)

(8)
(9)

i
wi zk|k−1
is the predicted measurement

i
vector and zk|k−1
= h(χik|k−1 ); the self and cross-covariance
zz
xz
matrices, Pk|k−1 and Pk|k−1
, are respectively calculated by
zz

Pk|k−1 =

2n
X
i=1

m×1

where xk ∈ R
and zk ∈ R
are the state vector
and the measurement/observation vector at time sample k,
respectively; f and h are vector-valued nonlinear functions;
wk and vk are the system process and observation noise, respectively; they are assumed to be independent and identically
distributed with zero mean and covariance matrices Qk and
Rk , respectively; uk is the system input vector.

bk|k−1 =
x

xz
Pk|k−1
=

i
i
− zbk|k−1 )T + Rk ,
− zbk|k−1 )(zk|k−1
wi (zk|k−1

2n
X
i=1

(10)

i
bk|k−1 )(zk|k−1
wi (χik|k−1 − x
− zbk|k−1 )T . (11)

C. Motivation of the Use of a Robust UKF
If the system process and measurement noises obey a
bk−1|k−1 ,
Gaussian probability distribution, the filtered state, x
will follow a Gaussian distribution as well. In that case, the
bk−1|k−1
sample mean and the sample covariance matrix of x
will be captured by the sigma points and the UKF will produce
reliable state estimates. However, the Gaussianity assumption
may not hold true in practice. This is precisely the case in
power systems; for instances, impulsive process noise may
occur due to system model inaccuracy at a certain time window
and the PMU measurement noise may not follow a Gaussian
distribution as shown in Fig.1. Consequently, the sigma points
may not capture the complete statistics of the state vector,
resulting in poor or even diverged estimations. Furthermore,
since the UKF lacks statistical robustness, it is sensitive to any

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. , NO. , 2017

4

type of outliers, including observation, innovation and structural outliers. In power system DSE, observation outliers refer
to the phase biases and gross errors in PMU measurements
[14]; innovation outliers may be induced by incorrect generator
parameter values, failure of brushless exciter rotating diodes,
or impulsive system process noise; and structural outliers may
be caused by transmission parameter errors or topology errors.
In the following section, we will propose a robust GM-UKF
that is able to suppress observation and innovation outliers and
to filter out various types of thick-tailed measurement noises.
Note that the problem of the identification and suppression
of structural outliers is outside the scope of this paper since
it requires a different formulation; it will be addressed in a
future work.
III. T HE P ROPOSED GM-UKF
Our GM-UKF consists of four major steps, namely a batchmode regression form step, a robust pre-whitening step, a
robust regression state estimation step, and a robust error
covariance matrix updating step. They are described next.
A. Derivation of the Batch-Mode Regression Form
In this subsection, we first show the equivalence of statistical
linearization and the unscented transformation using sigma
points. We then derive the proposed batch-mode regression
form. The former claim is presented in the following Theorem:
bk−1|k−1 and
Theorem 1. Given the state estimate vector x
xx
its associated covariance matrix Pk−1|k−1
, statistical linear
regression applied to an arbitrary nonlinear function g(x)
yields results that are equivalent to those of the unscented
transformation using the sigma points generated according to
(3).
Proof. Consider a nonlinear function y = g(x) evaluated in
2n points, i.e., (χi , γi ), where γi = g (χi ) for i= 1,..., 2n.
Assuming that the nonlinear function is statistically linearized
b and b
as y = Ax + b + ζ, the objective is to find A
b so that
the point-wise linearization error ζi is minimized, i.e.,
2n
n
o
X
b b
(12)
wi ζi T ζi ,
A,
b = arg min
i=1

where ζi = γi − (Aχi + b). By taking the derivative of the
objective function with respect to A and b and let them equal
to zero, respectively, we obtain
b
b = y − Ax,
(13)
where x =
2n
P

i=1

2n
P

T
−1
b = Pxy
A
Pxx
,

wi χi ; y =

i=1

2n
P

wi g (χi ) =

i=1
T

(14)

wi (χi − x)(χi − x) ; Pxy =

2n
P

i=1

2n
P

wi γi ; Pxx =

i=1

wi (χi − x)(γi − y)T .

Then, the estimation error covariance matrix is calculated as
2n
P
wi ζbi ζbiT
Pζζ =
i=1


T
2n
P
b (χi − x) γi − y − A
b (χi − x)
wi γi − y − A
=
i=1

bT = Pyy − P T P −1 Pxy ,
b xx A
= Pyy − AP
xy xx

(15)

where Pyy =

2n
P

i=1

T

wi (γi − y)(γi − y) . Now, by taking the

expectation and the outer product of the statistical linearized
model, respectively, we obtain the posterior statistics given by
b +
yb = Ax

2n
X
i=1

b =
wi χi − Ax

2n
X

wi χi ,

(16)

i=1

bT + Pζζ
b xx A
Pyy = AP
2n
P
T
T
−1
T
−1
wi (γi − y)(γi − y) − Pxy
= Pxy
Pxx
Pxy +
Pxx
Pxy
=

2n
P

i=1

i=1

T

wi (γi − y)(γi − y) ,

(17)
which are the same expressions as those obtained by applying
the unscented transformation to the nonlinear function y =
g(x). Thus, the proof is completed.
b is no longer the
Remark. In statistical linearization, A
Jacobian matrix of g(x) at a given point. The error covariance
matrix Pζζ is used to compensate the linearization errors
of the higher order Taylor series expansion terms. This is
however explicitly contained in the unscented transformation
process.
By applying statistical linearization to the nonlinear system
bk|k−1
process model, we obtain the predicted state vector x
xx
bk|k−1 =
along with its covariance matrix Pk|k−1
. We define x
xk − δk , where
state vector; δk is the prediction
 xk is the true
xx
error; and E δk δkT = Pk|k−1
. Then, statistical linearization
can be applied to the nonlinear observation equation, yielding


bk|k−1 + h x
bk|k−1 + νk + εk , (18)
zk = Hk xk − x

xz
xx
where Hk = (Pk|k−1
)T (Pk|k−1
)−1 , which is no longer
a Jacobian matrix. Here, the covariance
of the statistical

e k = E νk ν T = P zz
linearization error term is R
−
k
k|k−1
xz
xx
xz
zz
xz
(Pk|k−1
)T Pk|k−1
Pk|k−1
, where Pk|k−1
and Pk|k−1
are two
covariance matrices that are calculated by following the same
steps as those of the UKF. By processing the predictions and
the observations simultaneously, we get the following batchmode regression form:




 
bk|k−1 − h(b
xk|k−1 )
zk + Hk x
νk + εk
Hk
xk +
=
bk|k−1
−δk
I
x
(19)
which can be rewritten in a compact form as

fk xk + eek ,
zek = H

(20)

and the error covariance matrix is given by




Σk|k−1
0
T
= Sk SkT , (21)
Wk = E eek eek =
xx
0
Pk|k−1

T
ek;
where Σk|k−1 = E[(νk + εk )(νk + εk ) ] = Rk + R
I is an identity matrix; Sk is calculated by the Cholesky
decomposition technique.

Theorem 2. The weighted least squares estimator of the
batch-mode regression form (20) yields an estimated state
xx
vector x
bk|k and its associated covariance matrix Pk|k
that
are equivalent to those of the UKF.

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. , NO. , 2017

5

lj2

lj2
medj(lj2)

Confidence ellipse

Outliers

Projection

M
lj
medj(lj1)

Projected values

lj1

Median value

lj1

Fig. 2. Application of the projection statistics to the matrix Zk for detecting outliers in a two-dimensional dataset that consists of the innovation vectors and
the predicted state vectors.

Proof. It is well-known that the state estimate of (20) using
the weighted least squares estimator is given by
−1

fT Wk zek ,
fk
fT Wk H
bk|k = H
(22)
x
H
k
k


−1

xx
fk
fT Wk H
. By apwith the covariance matrix Pk|k
= H
k
plying an algebraic substitution and using the matrix inversion
lemma, we get

−1 −1

xx
T −1
xx
Pk|k = Hk Rk Hk + Pk|k−1
−1

xx
xx
xx
xx
Hk Pk|k−1
HkT + Rk
= Pk|k−1
− Pk|k−1
HkT Hk Pk|k−1
xx
xx
zz
T
= (I − Kk Hk ) Pk|k−1 = Pk|k−1 − Kk Pk|k−1 Kk ,
(23)
where the gain matrix is expressed as
xx
xx
xz
zz
Kk = Pk|k−1
HkT (Hk Pk|k−1
HkT + Rk )−1 = Pk|k−1
(Pk|k−1
)−1
(24)
Thus, we can conclude that the estimation error covariance
is identical to that of the UKF in (9). By applying similar
substitutions and using the matrix inversion lemma, we can
also show that the estimated state vector is given by

bk|k = x
bk|k−1 + Kk zk − zbk|k−1 ,
(25)
x

which completes the proof.

Remark. In the literature, a few Huber estimator-based robust
UKF methods have been proposed and applied to various
applications in signal processing, target tracking, to name a
few [20]–[22]. However, in their developed regression models,
νk that compensates higher order Taylor series expansion
error terms is neglected completely. As a consequence, the
estimation results are biased. In addition, they are unable
to handle innovation outliers and filter out non-Gaussian
measurement noise.
B. Robust Prewhitening
Before carrying out a robust regression, we uncorrelate the
state prediction errors of the batch-mode regression form. This
is achieved by pre-multiplying Sk−1 on both sides of (20),
yielding
fk xk + S −1 eek ,
Sk−1 zek = Sk−1 H
(26)
k

which can be further organized to the compact form
yk = Ck xk + ξk ,

(27)

T

where E[ξk ξk ] = I. However, if outliers occur, the application of Sk−1 will corrupt the prewhitening [12]. To overcome
this problem, we first detect the outliers and calculate the
weights using the projection statistics (PS) [12], [23]. Those
weights will be incorporated in the objective function that
is defined in the proposed GM-estimator shown in Section
III-C. Now, we describe the procedures used to calculate the
weights. We apply the PS to a 2-dimensional matrix Zk that
contains serially correlated samples of the innovations and of
the predicted state variables. Note that the innovation vector
is defined as the difference between the observations and their
associated predicted values at the previous step. Formally, we
have


zk−1 − h(b
xk−1|k−2 ) zk − h(b
xk|k−1 )
Zk =
, (28)
bk−1|k−2
bk|k−1
x
x

where zk−1 −h(b
xk−1|k−2 ) and zk −h(b
xk|k−1 ) are the innovation vectors while x
bk−1|k−2 and x
bk|k−1 are the predicted state
vectors at time instants k-1 and k, respectively. We may also
apply the PS to higher dimensional samples, but we found
that 2 dimensions are enough to identify outliers. The PS
values of the predictions and of the innovations are separately
calculated because the values taken by the former and the latter
are centered around different points.
The PS of the jth row vector, ljT , of the predictions (respectively the innovations) in Zk is defined as the maximum of
the standardized projections of all the ljT ’s on every direction
ℓ that originates from the coordinatewise medians of the predictions (respectively the innovations) and that passes through
every data point, and where the standardized projections are
based on the sample median and the median-absolute-deviation
[23]. The implementation of the PS to detect outliers in matrix
Zk is displayed in Fig. 2, while its mathematical expression
is defined as [23].

ljT ℓ − medi liT ℓ
 ,
(29)
P Sj = max
T ℓ − med lT ℓ
kℓk=1 1.4826 medκ lκ
i i
where i, j, κ = 1, 2, ..., m + n.
Once the PS values are calculated, they are compared to
a statistical threshold to identify outliers. Extensive Monte
Carlo simulations and Q-Q plots reveal that the probability

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. , NO. , 2017

6

25

12
20

Sample quantile

Sample quantile

10

8

6

15

10

4
5

2

0
0

2

4
6
8
2
Theoretical quantile−χ distribution

10

12

0

0

5

10

15

20

25

2
4

Theoretical quantile−χ distribution

2

Fig. 3. Q-Q plots of the sample quantiles of the PS vs. the corresponding quantiles of the χ22 and χ24 distributions, where (a) and (b) represent Q-Q plots of
PS with Gaussian and Laplace noise, respectively.

distributions of the PS applied to Zk , whose data points
obey bivariate Gaussian and Laplace probability distributions,
follow chi-square distributions with degree of freedom 2 and
4, respectively (See Fig. 3). This investigation allows us to
apply statistical tests to the PS and to flag all the data points
that satisfy PSi > η as outliers. The latter are downweighted
via


(30)
̟i = min 1, d2 P Si2 ,
where the parameter d is set equal to 1.5 to yield good statistical efficiency at different distributions without increasing too
much the bias induced by outliers. As an example, when the
noise is assumed to be Laplacian, the PS obeys a chi-square
distribution with 4 degrees of freedom. In that case, we can
choose the statistical detection threshold η as χ24,0.975 at a
significance level of 97.5%.
C. Robust Filtering and Solution
To suppress the outliers and filter out thick-tailed nonGaussian measurement noise, we develop a robust GMestimator that minimizes the following objective function:
J (xk ) =

m+n
X

̟i2 ρ (rSi ) ,

(31)

i=1

where ̟i is calculated by (30); rSi = ri /s̟i is the standardb is the residual, where cTi is the ith
ized residual; ri = yi −cTi x
row vector of the matrix Ck ; s = 1.4826 · bm ·mediani |ri | is
the robust scale estimate; bm is a correction factor to achieve
unbiasedness for a finite sample of size m + n at a given
probability distribution; ρ(·) is the nonlinear function of rSi .
In this paper, the convex Huber-ρ function [24] is adopted,
that is
 1 2
for |rSi | < λ
2 rSi ,

ρ (rSi ) =
,
(32)
λ |rSi | − λ2 2, elsewhere
where the parameter λ between the quadratic and the linear
segment of ρ(·) is typically chosen between 1.5 to 3 in the
literature.
To minimize (31), one takes its partial derivative with
respect to xk and sets it equal to zero, yielding
m+n
X ̟i c i
∂J (xk )
−
=
ψ (rSi ) = 0,
∂xk
s
i=1

(33)

where ψ (rSi ) = ∂ρ (rSi )/∂rSi is the so-called ψ-function. By
dividing and multiplying the standardized residual rSi to both
sides of (33) and putting it in a matrix form, we get
b (yk − Ck xk ) = 0,
CkT Q

(34)

b =diag(q (rS )) and q (rS ) = ψ (rS )/rS .
where Q
i
i
i
i
By using the IRLS algorithm [25], [26], the state vector
correction at the j iteration is calculated through
−1

(j+1)
b (j) yk ,
b (j) Ck
CkT Q
(35)
∆b
xk|k = CkT Q
(j+1)

where ∆b
xk|k

(j+1)

when ∆b
xk|k

(j+1)

bk|k
= x
∞

(j)

bk|k . The algorithm converges
−x

≤ 10−2 .

D. Asymptotic Error Covariance Matrix of the GM-UKF State
Estimates
Upon convergence of the iterative algorithm, the error
xx
covariance matrix Pk|k
is updated so that the state prediction
for the next step can be performed. To this end, consider the
ǫ-contamination model G = (1 − ǫ) Φ + ǫ∆r , where G and
Φ are the contaminated and the true cumulative probability
distribution function of the residual, respectively; ∆r is the
point mass to model outliers or unknown non-Gaussian distributions. The error covariance matrix is updated based on the
following theorem:
Theorem 3. Let T (·) be the functional form of the GMestimator with a bounded ψ(·) function and Φα be the
empirical cumulative probability distribution function, then
√
d
xx
),
(36)
α(T (Φα ) − T (Φ)) → N (0, Pk|k
d

where α = m + n; → means convergence in probability
xx
distribution; Pk|k
= E[IF (x; Φ, T ) · IF (x; Φ, T )T ] with the
influence function IF (x; Φ, T ) evaluated at Φ.
Proof. By taking a first-order Taylor series expansion of the
functional form of the estimator T with respective to Φ, we
get
T (Φα ) = T (Φ) + T ′ (Φα − Φ) + Rem(Φα − Φ),

(37)

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. , NO. , 2017

7

which can
√ be reorganized into the following form by multiplying α on both sides of the equality:
√
α (T (Φα ) − T (Φ))
√
√
(38)
= αT ′ (Φα − Φ) + αRem (Φα − Φ)
Z
√
√
= α IF (x; Φ, T )d(Φα − Φ) + αRem(Φα − Φ)(39)
Z
√
√
(40)
= α IF (x; Φ, T )dΦα + αRem(Φα − Φ)

Let’s now derive the IF (x; Φ, T ) of our GM-UKF at the
cumulative probability distribution Φ.
Corollary 3.2. The total influence function of the GM-UKF
defined by (33) using the regression model (27) is expressed
as
−1
Z
1 ′
ψ (rSi ) CC T T (Φ) dΦ
IF (x; Φ, T ) =
̟Cψ (rSi ) .
s
(45)

α

√
1 X
IF (xi ; Φ, T ) + αRem (Φα − Φ) ,
=√
α i=1

(41)

where the definition of the influence function is applied to
yield (38) to (39); byR virtue of Fisher consistency at the
distribution Φ, that is, IF (x; Φ, T )dΦ = 0, (39) reduces to
(40); finally, by using the property of the empirical cumulative
probability distribution function, we have
Z
α
1X
IF (xi ; Φ, T ),
(42)
IF (x; Φ, T )dΦα =
α i=1
yielding (40) to (41).
Following the work of Fernholz [27], we can show that
√
p
αRem (Φα − Φ) → 0,
(43)

Proof. In our previous work [6], the total influence function of
a GM-estimator based on a nonlinear regression model given
by y = ϕ(x) + ξ is expressed as
IF (x; Φ, T )

  ′
R ψ (rSi ) ∂ϕ(x) ∂ϕ(x) T
− ̟ψ (rSi ) D
=
s
∂x
∂x

T (Φ) dΦ

−1

· ̟ ∂ϕ(x)
∂x ψ (rSi ) ,

(46)
2
where D = ∂∂xϕ(x)
is
the
Hessian
matrix
of
ϕ(x).
Since
we
i ∂xj
have ϕ(x) = Cx for the GM-UKF, (46) reduces to
IF (x; Φ, T ) =

Z

1 ′
ψ (rSi ) CC T
s

dΦ
T (Φ)

−1

̟Cψ (rSi ) .
(47)

p

where → means probability convergence. Therefore, by applying the central limit theorem and Slutsky’s lemma to (41),
it follows that
√
d
xx
),
(44)
α(T (Φα ) − T (Φ)) → N (0, Pk|k
where

xx
Pk|k

T

= E[IF (x; Φ, T ) · IF (x; Φ, T ) ].

Discussion: The UKF is able to provide good results only
when the process and observation noises obey a Gaussian
distribution [19]. Indeed, in that case the filtered state vector
bk|k is Gaussian and the mean and covariance matrix of
x
bk|k can be accurately estimated by the sample mean and the
x
sample covariance matrix of the sigma points. However, this
property no longer holds true if the Gaussianity assumption
of the noises is violated. In that case, the state estimate vector
bk|k obtained from the UKF is significantly biased due the
x
filter lack of statistical robustness to thick-tailed non-Gaussian
noise. By contrast, our GM-UKF guarantees the asymptotic
bk|k for thick-tailed non-Gaussian noises and
Gaussianity of x
yields reliable state estimates with good statistical efficiency.
Corollary 3.1. Assume that the system process noise is
contaminated about a Gaussian distribution. Then, the data
points defined by the row vectors of the matrix Zk follow
asymptotically a Gaussian distribution.

Proof. From the definition of the matrix Zk given by (28)
and Theorem 3, we can see that the predicted state vector
is roughly Gaussian. Furthermore, under the assumption that
the minority of the measurements obey a thick-tailed nonGaussian distribution, the innovation vectors can be shown to
be approximately Gaussian. From this, we conclude that Zk
is asymptotically Gaussian.

Now, we are in a position to derive the covariance matrix
xx
Pk|k
from (44). First, let us prove the following theorem:
Theorem 4. The sample variance of the robust scale estimator
s of the GM-standardized residuals tends to one as the number
of observation tends to infinity.
Proof. By the law of large numbers, the distribution
of the residuals tends to the Gaussian distribution, i.e.,
Φ ∼ N (µ, σ 2 ). Since the median absolute deviation
s
(MAD= 1.4826·b
) is a consistent estimator for the standard
m
deviation σ of a Gaussian distribution, we get




|X−µ|
s
s
1
=P
|X
−
µ|
≤
≤
=
P
2
1.4826·b
σ
1.4826·bm ·σ
m


s
= 2Φ 1.4826·b
−
1.
m ·σ
(48)

Therefore, we obtain s/σ = 1.4826 · bm · Φ−1 43 = bm → 1
as m tends to infinity, where Φ is the cumulative probability
function of the standard Gaussian distribution. On the other
hand, from the equation (27) and the fact that E[ξk ξk T ] = I,
the residuals can be shown to actually follow the standard
Gaussian distribution. Therefore, EF [s2 ] = s2 → 1.
Finally, by Theorems 3 and 4, the asymptotic error covariance matrix of our GM-UKF at time sample k is updated by
xx
Pk|k
= E[IF (x; Φ, T ) · IF (x; Φ, T )T ]
−1

−1 T
EΦ [ψ 2 (rSi )]
T
Ck Q̟ Ck CkT Ck
=
2 Ck Ck
′ r
E
ψ
{ Φ [ ( Si )]}
(49)

where Q̟ = diag ̟i2 .

IEEE TRANSACTIONS ON POWER SYSTEMS, VOL. , NO. , 2017

E. Discussions on the Statistical Efficiency of the GM-UKF
In this section, we discuss the statistical efficiency of our
proposed GM-UKF under various probability distributions of
the noise. Firstly, under Gaussian measurement noise, the
outliers detected by the PS will be downweighted by the linear
segment of the ρ-function while all the good measurements
will be assigned weights equal or close to one since most
of them will be processed by the quadratic segment of the
ρ-function. As a result, the state estimator exhibits a high
statistical efficiency. Secondly, if the measurement noise obeys
a Laplace distribution, those measurements associated with the
thick tails of that distribution will have standardized residuals
corresponding to the linear segment of the ρ-function. This
means that for them, the GM-estimator behaves like the least
absolute value estimator; since the latter is the maximumlikelihood estimator at that distribution, it will have a high
asymptotic statistical efficiency. On the other hand, when the
estimation error covariance matrix is updated, all the outliers
with respect to the Gaussian distribution, which include the
measurements associated with the tails of the Laplacian or the
Cauchy distribution, will be heavily downweighted through
the matrix Q̟ , yielding bounded biases and variances in the
state estimates.
IV. C ONCLUSION
In this first part of a two-part series paper, we present
the fundamental theory of the proposed GM-UKF. We show
first that the UKF estimates the state vector via a weighted
least squares estimator under the Gaussianity assumption of
the system process or measurement noises; consequently, it
yields strongly biased state estimates when the noises follow
non-Gaussian probability distributions, which is precisely the
case when processing PMU measurements. By contrast, the
state estimates and residuals of our GM-UKF are proved
to be asymptotically Gaussian, allowing the sigma points to
reliably approximate the mean and the covariance matrices
of the predicted and corrected state vectors. Furthermore,
by relying on the projection statistics and the GM-estimator,
the proposed GM-UKF is able to suppress observation and
innovation outliers while exhibiting high statistical efficiency
of the state estimates. In addition, we derive the expression
of the asymptotic error covariance matrix of the GM-UKF
state estimates from the total influence function of the GMestimator. In the companion paper, we will discuss the implementation of our GM-UKF in power systems and analyze
its performance by carrying out extensive simulations under
various scenarios.
R EFERENCES
[1] I. Kamwa, R. Grondin, and Y. Hebert, “Wide-area measurement based
stabilizing control of large power systems-a decentralized/hierarchical
approach,” IEEE Trans. Power Syst., vol. 16, no. 1, pp. 136–153, Feb.
2001.
[2] H. Ni, G. T. Heydt, L. Mili, “Power system stability using robust wide
area control,” IEEE Trans. Power Syst., Vol. 17, No. 4, pp. 1123–1131,
Nov. 2002.
[3] E. Ghahremani, I. Kamwa, “Dynamic state estimation in power system
by applying the extended Kalman filter with unknown inputs to phasor
measurements,” IEEE Trans. Power Syst., vol. 26, no. 4, pp. 2556–2566,
Nov. 2011.

8

[4] Z. Huang, K. Schneider, and J. Nieplocha, “Feasibility studies of applying Kalman filter techniques to power system dynamic state estimation,”
in Proc. 8th Int. Power Eng. Conf., Singapore, Dec. 2007, pp. 376–382.
[5] L. Fan and Y. Wehbe, “Extended Kalman filtering based real-time
dynamic state and parameter estimation using PMU data,” Elect. Power
Syst. Res., vol. 103, pp. 168–177, Oct. 2013.
[6] J. B. Zhao, M. Netto, L. Mili, “A robust iterated extended Kalman filter
for power system dynamic state estimation,” IEEE Trans. Power Syst.,
DOI: 10.1109/TPWRS.2016.2628344, 2016.
[7] A. K. Singh, B. C. Pal, “Decentralized dynamic state estimation in power
systems using unscented transformation,” IEEE Trans. Power Syst., vol.
29, no. 2, pp. 794–804, Sep 2014.
[8] S. Wang, W. Gao, A. P. S. Meliopoulos, “An alternativemethod for power
system dynamic state estimation based on unscented transform,” IEEE
Trans. Power Syst., vol. 27, no. 2, pp. 942–950, May 2012.
[9] A. Rouhani, A. Abur, “Linear phasor estimator assisted dynamic state
estimation,” IEEE Trans. Smart Grid, 2016.
[10] Z. Huang, N. Zhou, R. Diao, S. Wang, S. Elbert, D. Meng, S. Lu, “Capturing real-time power system dynamics: opportunities and challenges,”
Proc. IEEE Power Eng. Soc. General Meeting, July 2015.
[11] N. Zhou, Z. Huang, D. Meng, “Capturing dynamics in the power
grid: formulation of dynamic state estimation through data assimilation.
Technical Report PNNL-23213, Pacific Northwest National Laboratory,
2014.
[12] M. Gandhi, L. Mili, “Robust Kalman filter based on a generalized
maximum-likelihood-type estimator, IEEE Trans. Signal Processing, vol.
58, no. 5, pp. 2509–2520, 2010.
[13] K. Martin, J. Hauer, T. Faris, “PMU testing and installation considerations at the Bonneville power administration,” Proc. IEEE Power Eng.
Soc. General Meeting, Jun. 2007, pp. 1–6.
[14] K. D. Jones, A. Pal, J. S. Thorp, “Methodology for performing synchrophasor data conditioning and validation,” IEEE Trans. Power Syst.,
vol. 30, no. 3, pp. 1121–1130, May. 2015.
[15] D. N. Kosterev, C. W. Taylor, W. A. Mittelstadt, ”Model Validation for
the August 10, 1996 WSCC System Outage,” IEEE Trans. Power Syst.,
vol. 14, no. 3, pp. 967–979, 1999.
[16] Z. Huang, P. Du, D. Kosterev, S. Yang, “Generator dynamic model
validation and parameter calibration using phasor measurements at the
point of connection,” IEEE Trans. Power Syst., Vol. 28, no. 2, pp. 1939–
1949, 2013.
[17] M. A. M. Ariff, B. C. Pal, A. K. Singh, “Estimating Dynamic Model
Parameters for Adaptive Protection and Control in Power System,” IEEE
Trans. Power Syst., vol. 30, no. 2, pp. 829–839, Mar. 2015.
[18] N. Zhou, D. Meng, Z. Huang, G. Welch, “Dynamic state estimation of
a synchronous machine using PMU data: A comparative study,” IEEE
Trans. Smart Grid., Vol. 6, no. 1, pp: 450–460, 2015.
[19] S. Julier, J. Uhlmann, and H. F. Durrant-Whyte, “A new method for
the nonlinear transformation of means and covariances in filters and
estimators,” IEEE Trans. Autom. Control, vol. 45, no. 3, pp. 477–482,
Mar. 2000.
[20] L. Y. Chang, B. Hu, G. K. Chang, et al, “Huber-based novel robust
unscented Kalman filter,” IET Science, Measurement & Technology, vol.
6, no. 6, pp. 502–509, 2012.
[21] L. B. Chang, B. Hu, G. K. Chang, et al, “Robust derivative-free Kalman
filter based on Huber’s M-estimation methodology,” Journal of Process
Control, Vol. 23, no. 10, pp. 1555-1561, 2013.
[22] L. Y. Chang, K. Li, B. Hu, “Huber’s M-estimation based process
uncertainty robust filter for integrated INS/GPS,” IEEE Sensors Journal,
Vol. 15, no. 6, 2015.
[23] L. Mili, M. Cheniae, N. Vichare, and P. Rousseeuw, “Robust state
estimation based on projection statistics,” IEEE Trans. Power Syst., vol.
11, no. 2, p. 1118–1127, 1996.
[24] P. J. Huber, Robust Statistics. New York: Wiley, 1981.
[25] L. Thomas, L. Mili, “A robust GM-estimator for the automated detection
of external defects on barked hardwood logs and stems,” IEEE Trans.
Signal Processing, vol. 55, no. 7, pp. 3568–3576, 2007.
[26] F. R. Hampel, E. M. Ronchetti, P. J. Rousseeuw, and W. A. Stahel,
Robust Statistics: The Approach Based on Influence Functions. New
York: John Wiley & Sons, Inc., 1986.
[27] L. Fernholz, “Von mises calculus for statistical functionals,” Lecture
Notes in Statistics. New York: Springer-Verlag, 1983, vol. 19.

