A Two-Step Distribution System State Estimator
with Grid Constraints and Mixed Measurements
Miguel Picallo∗† , Adolfo Anta∗ and Ara Panosyan∗
∗ General

arXiv:1703.10815v1 [] 31 Mar 2017

Electric Global Research
Munich, Germany
Email: miguel.picallocruz|anta|panosyan@ge.com

Abstract—In this work we study the problem of State Estimation (SE) in large-scale, 3-phase coupled, unbalanced distribution
systems. More specifically, we address the problem of including
mixed real-time measurements, synchronized and unsynchronized, from phasor measurement units and smart meters, into
existing solutions. We propose a computationally efficient twostep method to update a prior solution using the measurements,
while taking into account physical constraint caused by buses
with no loads. We test the method on a benchmark test feeder
to illustrate the effectiveness of the approach.

I. I NTRODUCTION
The operation of power networks requires accurate monitoring of its state: bus voltages, line currents, consumption and
generation. This is specially relevant in transmission networks
where volatile and distributed generation and consumption
cause bidirectional power flow and voltage drops. Normally,
State Estimation (SE) consists in estimating the bus voltage
phasors and then computing the currents and loads using
the power flow equations derived from the structure of the
network, represented by the admittance matrix. SE is typically
performed by taking several measurements and solving a
weighted least squares problem using an iterative approach
like Newton-Raphson [1]–[4].
However, SE has not yet been studied and implemented
deeply in distribution networks, since they used to have a
simple radial structure with a single source bus connected
to the grid and injecting power, being the point of common
coupling (PCC). Nowadays, this is changing, since the increasing number of sources of distributed generation like PVpanels, electrical vehicles, etc. inject power in the network,
causing bidirectional power flow. Consequently, SE becomes
necessary in distribution networks. One of the major difficulties for its implementation lies in the structural difference
between transmission and distribution networks: in the latter,
unbalanced loads and lower X/R ratios cause coupled phases.
Therefore, the power flow equations need to be solved for the
3 phases simultaneously, and fast methods like fast decoupled
power flow [5] cannot be used. Another limitation is the lack
of sufficient real-time measurement units in distribution networks: despite Phasor Measurement Units (PMU) and smart
meters are being introduced, their high cost prevent installing
the required number of sensors to make the system fully
This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No 675318 (INCITE).

Bart De Schutter†
† Delft

Center for Systems and Control
Delft, The Netherlands
Email: b.deschutter@tudelft.nl

observable. Some methods have been proposed to perform
SE by exploiting the network structure and nature of the
measurements, such as branch-current-based SE [6]–[8]. Other
methods use the few measurements available to improve the
SE accuracy of previous methods based on load assumptions
only [9]–[11]. Yet, most of this work is limited to transmission
networks, small-scale distribution networks, or single-phase
networks. Some recent work extends it more general networks,
but assumes measurement availability at every bus [12].
Our main contributions can be summarized as follows:
First, we include equality constraints more efficiently by a
dimension reduction. Then, we extend some recently proposed
methods that split the problem [9], to large-scale, 3-phase
coupled, unbalanced and constrained distribution systems, with
linear as well as nonlinear real-time measurements considering
different sources of measurement noises.
The rest of the paper is structured as follows: Section II
presents some background about power networks. Section III
discusses the different types of measurements. Section IV
presents the standard methodologies for general SE. Section V
proposes the linear dimension reduction method and extends
recent methodologies. Section VI defines the test case simulation and shows its results. Finally, Section VII presents the
conclusions and proposed future work.
II. D ISTRIBUTION S YSTEM M ODEL
A distribution system can be modeled as a graph G =
(V, E, W) with nodes V = {1, ..., Nbus }, edges E = {(vi , vj ) |
vi , vj ∈ V} and edge weights W. The nodes represent the
buses in the network, where power is injected or consumed, the
edges the electrical connections in the grid, and their weights
express the impedance between the buses, which is determined
by the length and type of the line cables.
In 3-phase networks buses may have up to 3 phases, so that
the voltage Vi at bus i lives in CPi , where Pi is the number of
its phases. Moreover, Vsource ∈ C3 are the known voltages at
2π T
the source bus with reference angles Vsource = [0 −2π
3
3 ] .
The state of the network is then typically represented with
the vector of the remaining bus voltages V ∈ CN where N
depends on the number of buses and phases per bus.
The impedance of the connection between buses i and j is a
square symmetric matrix Zi,j ∈ CPi,j ×Pi,j where Pi,j is the
number of phases in the connection i → j. Since phases are
coupled, Zi,j has off-diagonal terms different from 0. These
impedance matrices are used to compute the admittance matrix

Y ∈ C(N +3)×(N +3) , which defines the structure of the system.
Assuming lines with no shunt admittance, the element of Y
corresponding to the connection between phase l of bus i and
phase k of j can be expressed as:
(
− (Zi,j1 )l,k
if i 6= j
P
(Y )il ,jk =
(1)
1
if i = j
m6=i (Zi,m )l,k
where il is the index of phase l of bus i, similarly for jk ,
(·)i,j denotes the i, j element of a matrix, and (Y )il ,jk = 0
if no connection exists between these phases of these buses.
Using the admittance matrix Y , the formulas of the power
flow equations to compute the current I and the load S at all
buses are expressed as:



 


Isource
Vsource
Y Y
Vsource
=Y
= a b
I
V
Yc Yd
V
(2)
¯
S = diag(I)V
¯ denotes the complex conjugate, diag(·) represents
where (·)
the diagonal operator, converting a vector into a diagonal
matrix. Separating Y in blocks according to the indices of the
source bus Vsource , see (2), the voltage V for the non-source
buses can be rewritten as:
V = Yd−1 I − Yd−1 Yc Vsource
(3)
= Yd−1 (diag(V̄ ))−1 S̄ − Yd−1 Yc Vsource
so that the voltage under no load (I = 0, S = 0) is:
V0 = −Yd−1 Yc Vsource
(4)
III. M EASUREMENTS
Several different sources of information are available to
solve the SE problem:
1) Pseudo-measurements: Load estimations for every hour,
based on predictions and known installed load capacity
at every bus, and represented by Spseudo . Since they are
estimations rather than actual measurements, their noise is
modeled with a relative large standard deviation (a typical
value is σpseudo ≈ 50%):
S = Spseudo + <{Spseudo }ωpseudo,< + j={Spseudo }ωpseudo,=
with ωpseudo,< , ωpseudo,= ∼ N (0, σpseudo Id,N ), where
<{·}, ={·} denote the real and imaginary part, and Id,N
the identity matrix with size N × N . The noise covariance
of these measurements is then:
ΣSpseudo = E[(S − Spseudo )(S − Spseudo )∗ ]
(5)
2
2
= σpseudo
diag( Spseudo )
where E[·] denotes expectation, (·)∗ complex conjugate
transpose, |·| element-wise magnitudes of a complex vector
and (·)2 element-wise square.
2) Virtual measurements: Buses with no loads connected,
acting only as a connection [13]. They can be modeled
as physical constraints for the voltage states. Let ε =
{i, · · · , j} be the set with their indices, then:
(S)ε = 0, (I)ε = 0
(6)
where (·)ε denotes the elements at indices in ε.
3) Real-time measurements: Voltages and currents measurements from PMU, smart meters, or conventional remote
terminal units. There are two kinds: GPS-synchronized

measuring magnitude and phase angle, and unsynchronized
for only magnitude. We model the noises with a low
standard deviation for magnitude and angle, σmag ≈ 1%
and σang ≈ 0.01rad respectively, according to the IEEE
standard for PMU [14].
Synchronized measurements can be expressed using a linear approximation (see Appendix) with magnitude and angle noise caused by the measurements and by an imperfect
synchronization. For a number of NmeasL measurements
zmeasL ∈ CNmeasL we have:
zmeasL ≈ CmeasL V + diag(CmeasL V )(ωmag + jωang ) (7)
where ωmag , ωang ∼ N (0, σmeas Id,NmeasL ) using σmeas =
σang = σmag (since they are similar [14]), and CmeasL is
the matrix mapping state voltages to measurements. For
element j of CmeasL V measuring at phase l of bus i:
(C
V ) = (CmeasL )j,• V =
 measL j

for a voltage measurement

 Vil
(8)
(Y )il ,• V
for a current measurement

(Y
)
(V
−
V
)
for
a
branch-current
i
→
m
il ,ml
il
ml


measurement
where (·)j,• denotes row j.
Unsynchronized measurements have a nonlinear relation
CmeasNL (V ) with the voltages V . For NmeasNL measurements zmeasNL ∈ RNmeasNL :
zmeasNL = CmeasNL (V ) + diag(CmeasNL (V ))ωmeasNL (9)
with ωmeasNL ∼ N (0, σmeas Id,NmeasNL ). For element j of
CmeasNL (V ) measuring at phase l of bus i:
(CmeasNL (V ))j =


for a voltage measurement

 Vil
for a current measurement
(Y )il ,• V

(Y
)
(V
−
V
)
for a branch-current i → m
i
,m
i
m

l
l
l
l

measurement
(10)
Since the measurement noises are small, their covariance
matrices can be approximated using the measurements:
2
2
2
2
diag(|zmeasL | )
diag(|CmeasL V | ) ≈ 2σmeas
ΣmeasL = 2σmeas
2
2
2
diag(zmeasNL
)
diag(CmeasNL (V )2 ) ≈ σmeas
ΣmeasNL = σmeas
IV. S TANDARD M ETHODOLOGY FOR S TATE E STIMATION
In a power network, SE consists in estimating the state of
the network, represented by the bus voltages vector V . The
standard methodology for SE computes the maximum likelihood estimation by solving a constrained nonlinear weighted
least squares problem with all measurements and estimations:
(Spseudo )εc , zmeasL , zmeasNL , their noise covariance and the system constraints (Spseudo )ε = 0, where εc denotes the complimentary of ε. Typically this is solved using the NewthonRaphson method [1]. However, since this requires the use of
gradients and the power flow formulas in (2) are not holomorphic, i.e. complex differentiable, the problem is solved in real
T
variables using a polar representation Vpolar = [|V | ∠V T ]T :
min z − h(Vpolar ) W −1 s.t. g(Vpolar ) = 0
(11)
Vpolar
2
kxkA =

where
x∗ Ax is the norm with respect to the real
positive-definite matrix of weights A, z is the vector of mea-

surements in rectangular coordinates and h(·) is a nonlinear
function mapping Vpolar to measurements:




<{Sεc (V (Vpolar ))}
<{Spseudo,εc }


 ={Sεc (V (Vpolar ))} 
 ={Spseudo,εc } 







z=
 <{zmeasL }  , h(Vpolar ) =  <{(CmeasL V (Vpolar ))} 

 ={(C
 ={zmeasL } 
measL V (Vpolar ))} 

zmeasNL
CmeasNL (V (Vpolar ))
where S(V ) represents S as a function of V , V (Vpolar ) is
defined similarly. The function g(·) indicates the null loads:
g(Vpolar ) = [<{[Sε (V (Vpolar ))}T ={[Sε (V (Vpolar ))}T ]T
and the weight matrix W −1 is the inverse of the measurement
noises in rectangular variables



!
2
c
<{(S
)
}
pseudo ε
2
diag
0
0 
 σpseudo
={(Spseudo )εc }2


W =


0
Σrect,measL
0 
0
0
ΣmeasNL
where 0 denotes a matrix with 0 entries, and
2
Σrect,measL ≈ σmeas
"
#
2
diag(|zmeasL | )
2diag(<{zmeasL }={zmeasL })
·
2
2diag(<{zmeasL }={zmeasL })
diag(|zmeasL | )
This problem (11) is typically solved iteratively using Lagrangian multipliers λ [15], so that the size of the optimization
T
variables is given by [Vpolar
λ]T ∈ R2N +2|ε| .
V. P ROPOSED M ETHODS FOR S TATE E STIMATION
In this section we propose new SE methods that significantly
improve the standard methodology in (11). First, we present a
method that simplifies the problem and reduces its complexity,
while obtaining an equal solution; then, we propose another
method that significantly reduces the computation cost while
achieving approximately the same accuracy.
A. Weighted least squares in subspace
Previous methods either do not consider constraints, or
do it by introducing Lagrange multipliers as in (11), which
increases the complexity of the algorithm. Here we propose
an alternative method to include the constraints using a dimension reduction through a linear transformation that produces
feasible solutions. The space of feasible solutions is:
{V | Iε (V ) = (Yd )ε V + (Yc )ε Vsource = 0}
(12)
= {V | ∃x ∈ CN −|ε| s.t. V = F x + Vp }
where Vp satisfies Iε (Vp ) = 0, and F is such that:
(Yd )ε F = 0
(13)
Therefore, V0 is chosen for Vp , and F is built using an
orthonormal basis of the kernel of (Yd )ε e.g. by computing
the Singular Value Decomposition of (Yd )ε = USV ∗ , and
taking the columns of V with null singular value. Then,
F = ker((Yd )ε ) and F ∗ F = Id , where ker(·) denotes the
kernel subspace. If we consider the rectangular representation
of voltages Vrect = [<{V }T ={V }T ]T instead of the polar

representation Vpolar , the subspace equations can still be represented linearly in real variables:
{Vrect | <{Iε (Vrect )} = 0, ={Iε (Vrect )} = 0}
(



 )
<{(Yd )ε } −={(Yd )ε }
<{(Yc )ε Vsource }
V +
=0
= Vrect
={(Yd )ε } <{(Yd )ε } rect
={(Yc )ε Vsource }
= {Vrect | ∃x̃ ∈ C2N −2|ε| s.t. Vrect = F̃ x̃ + Ṽp }
(14)
where now



!
<{V0 }
<{(Yd )ε } −={(Yd )ε }
Ṽp =
, F̃ = ker
(15)
={V0 }
={(Yd )ε } <{(Yd )ε }
Then the problem in (11) becomes:

Vrect = F̃ arg min z − h̃(F̃ x̃ + Ṽp )
x


W −1

+ Ṽp

(16)

where h̃(·) is h(·) in rectangular coordinates, using V (Vrect )
instead of V (Vpolar ), so that h̃(Vrect ) = h(Vpolar ). This again
can be solved using the Newton-Raphson method:
x̃k+1 = x̃k + ∆x̃k
with
∆x̃k = (HkT W −1 Hk )−1 HkT W −1 (z − h̃(F̃ x̃k + Ṽp ))
Hk = ∇x̃ h̃(F̃ x̃ + Ṽp ) |x̃k = ∇Vrect h̃(Vrect ) |Vrect,k F̃
Note that now no Lagrangian multipliers are required. We
eliminate the equality constraints by embedding the solution
into a smaller subspace of feasible solutions. This simplifies
the method and reduces its computational cost, since we
reduce the size of the optimization variables x̃ ∈ R2N −2|ε| .
B. Power flow solution plus optimal linear update
The methods (11) and (16) are iterative, and thus can
become computationally expensive as the size of the network
N increases. This may be a limitation in large distribution
systems when fast calculations are required every few seconds.
Therefore, we propose an alternative method splitting the
problem into two parts as in [9]: first, we use only the
constraints and the load estimates Spseudo available beforehand
to solve the power flow problem offline and to obtain a prior
solution defined as Vprior ; then, we update in real-time the
solution using the measurements zmeasL , zmeasNL to compute a
better estimate, represented by Vpost . We also prove that several
alternatives for the update are equal under certain conditions.
1) Obtaining a prior solution:
a) The prior estimate can be computed using the previous
method in (16) without real-time measurements, so that
z = [<{Spseudo,εc }T ={Spseudo,εc }T ]T and h̃(Vrect ) =
[<{Sεc (V (Vrect ))}T ={Sεc (V (Vrect ))}T ]T . This method
will estimate real and imaginary parts Vrect,prior , the estimation error covariance of which can be approximated by:
−1
Σrect,prior ≈ F̃ (HkT Wrect
Hk )−1 F̃ T
(17)
The covariance matrix in complex variables is:
Σprior = Σrect,prior,<< + Σrect,prior,==
(18)
−jΣrect,prior,<= + jΣTrect,prior,<=
where Σrect,prior,<< , Σrect,prior,== indicate the covariance of
real and imaginary parts respectively, and Σrect,prior,<= the
covariance between real and imaginary parts.

b) An alternative approach using complex variables would be
the method proposed in [9]. This method can be extended
to a 3-phase coupled, unbalanced system using the voltage
V0 under no loads defined in (4). There, V0 is used as
initial value in this iterative method, where the iterations
are inspired by (3), i.e. for the k-iteration we have:
Vk+1 = Yd−1 (diag(V̄k ))−1 S̄ + V0
(19)
The estimation error covariance can be approximated by
the covariance after the first iteration, which corresponds
to the linear approximation in [16]:
Σprior ≈ Yd−1 (diag(V̄0 ))−1 ΣSpseudo (diag(V0 ))−1 (Yd−1 )∗
Remark 1. The estimated current at any iteration k is:
Ik+1 = (diag(V̄k ))−1 S̄
(20)
If the loads satisfy (S)ε = 0, then (Ik+1 )ε = 0. This means
that there exists a vector x ∈ CN−|ε| s.t. Vk = F x + V0
and thus Vk is feasible ∀k w.r.t. (12). Note that we use
F ∈ CN ×(N −|ε|) instead of F̃ since Vk ∈ CN .
Remark 2. This method can be extended to real variables
to get Σprior in rectangular coordinates:


<{Spseudo }
Vrect,k+1 = Bk
+ Vrect,0
(21)
={Spseudo }
with
# !!−1
"
#
"
2
<{Yd−1 } −={Yd−1 }
|Vk |
Bk =
diag
2
={Yd−1 } <{Yd−1 }
|Vk |


diag(<{Vk }) diag(={Vk })
·
diag(={Vk }) −diag(<{Vk })
and estimation error covariance:
!

2
<{S
}
pseudo
2
B0T (22)
Σrect,prior = B0 σpseudo
diag
={Spseudo }2
The advantage of approach b) w.r.t. a) is that b) can be solved
in complex variables. In both cases, since the load estimates
Spseudo are known beforehand, this part of the problem can
be computed offline, taking this iterative computational cost
outside the real-time operation of the system.
2) Updating with only synchronized measurements:
If there are not any unsynchronized measurements, the mapping from the state variables to the real-time measurements (8)
is linear and thus holomorphic in x ∈ CN −|ε| . Consequently,
the update can be solved using complex variables to obtain a
linear expression. Different alternative methods are possible,
here we present two of them and later show their equivalence:
a) If the prior solution Vprior is feasible, then there exists a
vector xprior s.t. Vprior = F xprior +V0 . It can be computed as
xprior = F ∗ (Vprior − V0 ), with estimation error covariance
F ∗ Σprior F , or xrect,prior = F̃ T (Vrect,prior − Vrect,0 ) for real
variables. Then the maximum likelihood update can be
computed using weighted least squares:

 2
x − xprior
Vpost = V0 +F arg min
CmeasL (F x + V0 ) − zmeasL
x
−1
Wpost

(23)

−1
where the optimal weights Wpost
are:
 ∗

F Σprior F 0
Wpost =
(24)
0
ΣmeasL
Taking derivatives at (23) we obtain the solution:
∗
−1
Vpost = F ((F ∗ Σprior F )−1 + F ∗ CmeasL
Σ−1
measL CmeasL F )
∗
·((F ∗ Σprior F )−1 xprior + F ∗ CmeasL
Σ−1
measL
·(zmeasL − CmeasL V0 )) + V0

(25)
Remark 3. Considering the first order (and thus linear)
approximation of the measurement function h̃(·), in [10] is
proven that solving the weighted least squares problem in
one step (16) and in two steps: weighted least squares for
power flow plus update linear (23), produce the same exact
solutions. Consequently, we can expect a similar accuracy
between the one step and the two steps methods.
b) An alternative proposed [9] is to improve the solution using
an minimum variance linear update:

Vpost = Vprior + K zmeasL − CmeasL Vprior
(26)
Then the error covariance of estimation Vpost is:
∗
Σpost = Σprior + K(ΣmeasL + CmeasL Σprior CmeasL
)K ∗
(27)
∗
∗
−KCmeasL Σprior − Σprior CmeasL K
and the optimal gain K is computed minimizing the
expected error E[(Vpost − V )∗ (Vpost − V ) | Vprior , zmeasL ] =
E[tr(Σpost ) | Vprior , zmeasL ]:
K = arg minK E[tr(Σpost ) | Vprior , zmeasL ]
(28)
∗
∗
= Σprior CmeasL
(CmeasL Σprior CmeasL
+ ΣmeasL )−1
Proposition 1. If Vprior is a feasible solution to (12), as the
solutions in Section V-B1, then the solution (26) satisfies (12).
Proof. For Vprior : (Yd )ε Vprior +(Yc )ε Vsource = 0. Also, since the
first term in K is F , because K starts with Σprior = F Σxprior F ∗ ,
we have (Yd )ε K = 0 due to (13). Therefore, Vpost from (26)
satisfies (12):
(Yd )ε Vpost + (Yc )ε Vsource
= (Yd )ε Vprior + (Yc )ε Vsource + (Yd )ε K(zmeasL − CmeasL V )
=0+0=0

Proposition 2. If Vprior is a feasible solution to (12), then both
alternatives (26) and (25) are equal.
Proof. We start from expression (25), using Woodbury’s identity [17] on the inverse term and (28) for F we get:
∗
Vpost = F ((F ∗ Σprior F ) − (F ∗ Σprior F )F ∗ CmeasL
∗
·(CmeasL Σprior CmeasL
+ ΣmeasL )−1 CmeasL F (F ∗ Σprior F ))
∗
·((F ∗ Σprior F )−1 xprior + F ∗ CmeasL
Σ−1
measL
·(zmeasL − CmeasL V0 )) + V0
= (F (F ∗ Σprior F ) − KCmeasL F (F ∗ Σprior F ))
∗
·((F ∗ Σprior F )−1 xprior + F ∗ CmeasL
Σ−1
measL
·(zmeasL − CmeasL V0 )) + V0
= V0 + (Id − KCmeasL )F xprior
∗
+(Id − KCmeasL )F (F ∗ Σprior F )F ∗ CmeasL
Σ−1
measL
·(zmeasL − CmeasL V0 ))

(29)

Doing some manipulation to the formula of K in (28) we get:
∗
(Id − KCmeasL )Σprior CmeasL
Σ−1
(30)
measL = K
So that finally we can convert the expression in (25) to (26):
Vpost = V0 + (Id − KCmeasL )F xprior + K(zmeasL − CmeasL V0 )
= Vprior + K(zmeasL − CmeasL Vprior )

The maximum likelihood estimator of weighted least square
coincides with the minimum variance estimator of the linear
update, since we consider gaussian distributions for the noise.
3) Updating with unsynchronized measurements:
If there are nonlinear measurements, since the magnitude of a
complex number as in (10) is not an holomorphic function, the
update step needs to be solved using real variables. It becomes:


<{zmeasL − CmeasL V (Vrect,prior )}
Vrect,post = Vrect,prior + K̃  ={zmeasL − CmeasL V (Vrect,prior )} 
zmeasNL − CmeasNL (V (Vrect,prior ))
(31)
Given that

 

<{CmeasL V (Vrect,prior )}
<{CmeasL } −={CmeasL }
=
V
={CmeasL V (Vrect,prior )}
={CmeasL } <{CmeasL } rect,prior
and defining
 
 
<{CmeasL } −={CmeasL }

CmeasLNL =  ={CmeasL } <{CmeasL }
∇Vrect CmeasNL (V (Vrect )) |Vrect,prior


Σrect,measL
0
ΣmeasLNL =
0
ΣmeasNL
The first-order approximation of the estimation error covariance of this nonlinear filter is:
T
Σrect,post ≈ K̃(ΣmeasLNL + CmeasLNL Σrect,prior CmeasLNL
)K̃ T
T
+Σrect,prior − K̃CmeasLNL Σrect,prior − Σrect,prior CmeasLNL
K̃ T
(32)
where now the optimal gain K̃ is:
T
T
K̃ = Σrect,prior CmeasLNL
(CmeasLNL Σrect,prior CmeasLNL
+ΣmeasLNL )−1

Remark 4. As in Proposition 2, it can be proven that the
first-order approximation of the estimation error covariance
(32) of (31) is equal to the one approximating the estimation
error covariance of (23) including nonlinear measurements,
computed similarly as in (17).
VI. S IMULATION R ESULTS
A simulation of 24 hours with 15 min intervals is run on
a test case to compare the different methods for SE. Here we
describe the settings of this simulation and analyze the results.
A. Settings
• System: We use the 123-bus test feeder [18] available online
[19], see Fig. 1. This is a challenging example, since it is 3phase coupled, unbalanced, and larger than other examples
in the literature [7], [9].
• Measurements (see Fig. 1): Voltage measurements (red
solid, circle for phasor, square for magnitude only) are
placed at buses 79, 95, 83 and 300, current measurements
(blue dashed, circle and square) at bus 65 and 48, and branch
current measurements (blue dashed arrow) at branch 150
(after the regulator) → 149. The standard deviation used to
simulate noise in the measurements is σmeas = 0.01 [14].

Fig. 1: 123-bus test feeder with measurements location. The
network image has been taken from [19].

•

•

•

This sensors locations are chosen since they correspond to
nodes with big loads and/or at the end of a feeder.
Load Profiles: Historical data from a German local utility is
scaled so that the mean load is the base load of the 123-bus
test feeder and the relative standard deviation is σpseudo =
50% at every hour. The values of Spseudo are an average of
the bus loads for every hour.
Metric: For every simulation time step t, a normalized root
mean square error (nRMSE) is determined for the voltages
estimated by every method j:
v
u
N
1 u
1 X
2
t
nRMSEt,Vmethodj =
Vmethodj ,t,i − Vt,i
|Vbase | N i=1
where Vt,i , Vmethodj ,t,i denote the i-element at time t of
the actual voltage and the voltage estimated by method j
respectively. This method quantifies the percentage deviation
of the solutions relative to the network’s base voltage Vbase .
Nomenclature: The solutions compared are denoted as:
VWLS , VWLSNL , which correspond to the solutions of (16)
without and with nonlinear measurements respectively; Vprior
corresponding the prior solution (19) before the update; and
Vpost , VpostNL corresponding the updated solutions (23), (31).

B. Results
Fig. 2a shows how a solution with real-time measurements
like Vpost clearly outperforms the prior solution Vprior using
only load estimations. This shows that even when few measurements are available, in this case only 4 measurements in
a 123-bus network, using them can dramatically increase the
accuracy of the prediction. Moreover, Fig. 2b shows that the
two-step solutions Vpost and VpostNL perform statistically as well
as the respective weighted least squares solutions VWLS and
VWLSNL . However, the methods using weighted least squares
are computationally much more expensive than the others,
specially the one with nonlinear measurements VWLSNL , see
Fig. 3 for the comparison of different computation times.
These algorithms are coded in Python and run on an Intel
Core i7-5600U CPU at 2.60GHz with 16GB of RAM.
To summarize, the solutions including real-time measurements
improve the accuracy of the estimation, specially if adding extra nonlinear measurements. However, the two-step approach

8

0.9

7

0.8
0.7

nRMSE [%]

nRMSE [%]

6
5
4
3
2

0.6
0.5
0.4
0.3
0.2

1
0

1
where ωmag ∼ N (0, 1%|u|) (or |u|
ωmag ∼ N (0, 0.01)),
ωθ ∼ N (0, 0.01rad), according to the standards in [14]. Using
trigonometric identities for cos(θu + ωθ ) and sin(θu + ωθ );
using approximations sin(ωθ ) ≈ ωθ and cos(ωθ ) ≈ 1,
since |ωθ |  1; and neglecting second order terms of noise
O(ωθ ωmag ), then we have:
ũ ≈|u| (cos(θu ) + j sin(θu )) + ωmag (cos(θu ) + j sin(θu ))
+ωθ |u| (− sin(θu ) + j cos(θu ))

1
1
= u + ωmag |u|
u + jωθ u = u + u |u|
ωmag + jωθ

0.1

Vprior

Vpost

0.0

VWLS

VWLSNL

(a)

Vpost

VpostNL

(b)

Fig. 2: Box-plots of MAPE errors for the different solutions.
The red line indicates the median, the red square the mean.
0.35
0.30

time [s]

0.25
0.20
0.15
0.10
0.05
0.00

VWLS

VWLSNL

Vpost

VpostNL

Fig. 3: Box-plots of execution times for the different solutions.
The red line indicates the median, the red square the mean.
has a smaller computational cost, and remains small after
adding nonlinear measurements.
VII. C ONCLUSIONS
We have proposed a method for state estimation in largescale, 3-phase coupled, unbalanced and constrained distribution systems, with mixed measurements: phasor-synchronized
(linear) and magnitude-unsynchronized (nonlinear), and different sources of noise. We have shown that the proposed method
is as accurate as standard methods, but computationally more
efficient in real-time operation, since it is a non-iterative
method. The computational cost reduction is achieved by
splitting the problem into an offline and an online problem.
The iterative offline method estimates the prior solution, taking
most of the computational cost, while the online problem
updates this solution in a single step.
Future work could include adding different sources of distributed generation and investigating how they affect the
solutions; developing methodologies to compute the minimum
and optimal allocation of measurement units to guarantee a
given accuracy level; and adding dynamic state equations to
exploit the historical information provided by the sensors.
A PPENDIX
A PMU measure of a complex number u = |u| ejθu with
magnitude and angle noise can be expressed as:
ũ = (|u| + ωmag )ej(θu +ωθ )
= (|u| + ωmag )(cos(θu + ωθ ) + j sin(θu + ωθ ))

R EFERENCES
[1] A. Abur and A. G. Exposito, Power System State Estimation: Theory
and Implementation. CRC Press, 2004.
[2] Y.-F. Huang, S. Werner, J. Huang, N. Kashyap, and V. Gupta, “State
estimation in electric power grids: Meeting new challenges presented by
the requirements of the future grid,” IEEE Signal Processing Magazine,
vol. 29, no. 5, pp. 33–43, 2012.
[3] A. Monticelli, “Electric power system state estimation,” Proceedings of
the IEEE, vol. 88, no. 2, pp. 262–282, 2000.
[4] B. Hayes and M. Prodanovic, “State estimation techniques for electric
power distribution systems,” in 2014 IEEE European Modelling Symposium (EMS), 2014, pp. 303–308.
[5] B. Stott and O. Alsaç, “Fast decoupled load flow,” IEEE Transactions
on Power Apparatus and Systems, no. 3, pp. 859–869, 1974.
[6] M. E. Baran and A. W. Kelley, “State estimation for real-time monitoring
of distribution systems,” IEEE Transactions on Power Systems, vol. 9,
no. 3, pp. 1601–1609, 1994.
[7] M. B. A. Kelley, “A branch current based state estimation method for
distribution systems,” IEEE Transactions on Power Systems, vol. 10, pp.
483–491, 1995.
[8] W.-M. Lin, J.-H. Teng, and S.-J. Chen, “A highly efficient algorithm in
treating current measurements for the branch-current-based distribution
state estimation,” IEEE Transactions on Power Delivery, vol. 16, no. 3,
pp. 433–439, 2001.
[9] L. Schenato, G. Barchi, D. Macii, R. Arghandeh, K. Poolla, and
A. Von Meier, “Bayesian linear state estimation using smart meters and
pmus measurements in distribution grids,” in 2014 IEEE International
Conference on Smart Grid Communications, 2014, pp. 572–577.
[10] M. Zhou, V. A. Centeno, J. S. Thorp, and A. G. Phadke, “An alternative
for including phasor measurements in state estimators,” IEEE Transactions on Power Systems, vol. 21, no. 4, pp. 1930–1937, 2006.
[11] Y. Hu, A. Kuh, T. Yang, and A. Kavcic, “A belief propagation based
power distribution system state estimator,” IEEE Computational Intelligence Magazine, vol. 6, no. 3, pp. 36–46, 2011.
[12] S. Bolognani, R. Carli, and M. Todescato, “State estimation in power
distribution networks with poorly synchronized measurements,” in 2014
IEEE 53rd Annual Conference on Decision and Control. IEEE, 2014,
pp. 2579–2584.
[13] G. N. Korres, N. D. Hatziargyriou, and P. J. Katsikas, “State estimation in multi-microgrids,” European Transactions on Electrical Power,
vol. 21, no. 2, pp. 1178–1199, 2011.
[14] K. Martin, D. Hamai, M. Adamiak, S. Anderson, M. Begovic, G. Benmouyal, G. Brunello, J. Burger, J. Cai, B. Dickerson et al., “Exploring
the IEEE standard C37. 118–2005 synchrophasors for power systems,”
IEEE Transactions on Power Delivery, vol. 23, no. 4, pp. 1805–1811,
2008.
[15] G. N. Korres, “A robust algorithm for power system state estimation
with equality constraints,” IEEE Transactions on Power Systems, vol. 25,
no. 3, pp. 1531–1541, 2010.
[16] S. Bolognani and S. Zampieri, “On the existence and linear approximation of the power flow solution in power distribution networks,” IEEE
Transactions on Power Systems, vol. 31, no. 1, pp. 163–172, 2016.
[17] M. A. Woodbury, “Inverting modified matrices,” Statistical Research
Group Memorandum Reports. Princeton, NJ: Princeton University,
no. 42, 1950.
[18] W. Kersting, “Radial distribution test feeders,” IEEE Transactions on
Power Systems, vol. 6, no. 3, pp. 975–985, 1991.
[19] DistributionTestFeeders, “IEEE PES Distribution Systems Analysis
Subcommittee’s Radial Test Feeders,” 1991, accessed: 2016-10-20.
[Online]. Available: ewh.ieee.org/soc/pes/dsacom/testfeeders.html

