1

Distributed Gauss-Newton Method for AC State
Estimation Using Belief Propagation

arXiv:1702.05781v2 [] 1 Oct 2017

Mirsad Cosovic, Student Member, IEEE, Dejan Vukobratovic, Member, IEEE

Abstract—We present a detailed study on application of
factor graphs and the belief propagation (BP) algorithm to the
power system state estimation (SE) problem. We start from the
BP solution for the linear DC model, for which we provide a
detailed convergence analysis. Using insights from the DC
model, we use two different approaches to derive the BP
algorithm for the non-linear AC model. The first method
directly applies BP methodology, however, providing only
approximate BP solution for the AC model. In the second
approach, we make a key further step by providing the
solution in which the BP is applied sequentially over the AC
model, akin to what is done by the Gauss-Newton method. The
resulting BP-based Gauss-Newton algorithm has the
interpretation of a fully distributed Gauss-Newton method with
the same accuracy as the centralized SE, preserving a number
of advantages of the BP framework. The paper provides
extensive numerical study of the proposed algorithms, provides
details on their convergence properties, and gives a number of
useful insights for their implementation.
Index Terms—State Estimation, Electric Power System, Factor
Graphs, Belief Propagation, Distributed Gauss-Newton Method

I. I NTRODUCTION
Motivation: Electric power systems consist of generation,
transmission and consumption spread over wide geographical
areas and operated from the control centers by the system
operators. Maintaining normal operation conditions is of the
central importance for the power system operators [1, Ch. 1].
Control centers are traditionally operated in centralized and
independent fashion. However, increase in the system size
and complexity, as well as external socio-economic factors,
lead to deregulation of power systems, resulting in
decentralized structure with distributed control centers.
Cooperation in control and monitoring across distributed
control centers is critical for efficient system operation.
Consequently, existing centralized algorithms have to be
redefined based on a new requirements for distributed
operation, scalability and computational efficiency [2].
The system monitoring is an essential part of the control
centers, providing control and optimization functionality
This paper was presented in part at the IEEE International Conference
on Smart Grid Communications (SmartGridComm), Sydney, Australia,
November 2016; in part at the IEEE International Workshop on Signal
Processing Advances in Wireless Communications (SPAWC), Edinburgh,
United Kingdom, August 2016.
M. Cosovic is with Schneider Electric DMS NS, Novi Sad, Serbia
(e-mail: mirsad.cosovic@schneider-electric-dms.com). D. Vukobratovic is
with Department of Power, Electronic and Communications Engineering,
University of Novi Sad, Novi Sad, Serbia (e-mail: dejanv@uns.ac.rs). Demo
source code available online at https://github.com/mcosovic.
This paper has received funding from the EU 7th Framework Programme
for research, technological development and demonstration under grant
agreement no. 607774.

whose efficiency relies on accurate state estimation (SE).
The core of the SE is the SE algorithm that provides an
estimate of the system state (i.e., the set of all complex bus
voltages) based on the network topology and available
measurements. Besides the SE algorithm, the SE includes
several additional routines such as network topology
processors, observability analysis and bad data analysis.
The centralized SE assumes that the measurements
collected across the system are available at the control
center, where the centralized SE algorithm provides the
system state estimate. Precisely, the centralized SE algorithm
typically uses the Gauss-Newton method to solve the
non-linear weighted least-squares (WLS) problem [3], [4]. In
contrast, decentralized SE distributes communication and
computational effort across multiple control centers to
provide the system state estimate. There are two main
approaches to distributed SE: i) algorithms which require a
global control center to exchange data with local control
centers, and ii) algorithms with local control centers only
[5]. Distributed SE algorithms target the same state estimate
accuracy as achievable using the centralized SE algorithms.
Literature review: The mainstream distributed SE
algorithms exploit matrix decomposition techniques applied
over the Gauss-Newton method. These algorithms usually
achieve the same accuracy as the centralized SE algorithm
and work with either global control center [6]–[8] or without
it [9]–[12]. Recently, the SE algorithms based on distributed
optimization [13], and in particular, the alternating direction
method of multipliers became very popular [14]. Examples
include algorithms based on convex relaxation [15], [16],
and without convex relaxation [17]. In [18], the gossip-based
Gauss-Newton algorithm is proposed which provides flexible
communication model, but suffers from slight performance
degradation compared to the centralized SE. The work in
[19] presents a fully distributed SE algorithm for wide-area
monitoring which provably converges to the centralized SE.
We refer the reader to [20] for a detailed survey of the
distributed multi-area SE.
Belief-Propagation Approach: In this paper, we solve the
SE problem using probabilistic graphical models [21], a
powerful tool for modeling the dependencies among the
systems of random variables. We represent the SE problem
using a popular class of probabilistic graphical models called
factor graphs and solve it using the belief propagation (BP)
algorithm. Applying the BP algorithm on probabilistic
graphical models without loops, one obtains exact marginal
distributions or a mode of the joint distribution of the system
of random variables [21], [22]. The BP algorithm can be

2

also applied to graphical models with loops (loopy BP) [23],
although in that case, the solution may not converge to the
correct marginals/modes of the joint distribution. BP is a
fully distributed algorithm suitable for accommodation of
distributed power sources and time-varying loads. Moreover,
placing the SE into the probabilistic graphical modelling
framework enables not only efficient inference, but also, a
rich collection of tools for learning parameters or structure
of the graphical model from observed data [24].
The work in [25] provides the first demonstration of BP
applied to the SE problem. Although this work is elaborate
in terms of using, e.g., environmental correlation via
historical data, it applies BP to a simple linearized DC
model. The AC model is recently addressed in [26], where
tree-reweighted BP is applied using preprocessed weights
obtained by randomly sampling the space of spanning trees.
The work in [27] investigates Gaussian BP convergence
performance for the DC model. Although the above results
provide initial insights on using BP for distributed SE, it is
fair to say that a systematic analysis of applying the BP
algorithm on the SE problem, and in particular for the AC
model, is still missing. This paper intends to fill this gap.
Contributions: In this paper, we present a novel
distributed BP-based Gauss-Newton algorithm. In the process
of deriving the proposed algorithm, we provide a
step-by-step guide for application of BP algorithm to the SE
problem, giving this paper strong tutorial flavor. Our
methodology is to start with the simplest linear DC model
and use insights obtained therein to derive the native BP
solution for the non-linear AC model. Unfortunately, as
closed-form expressions for certain classes of BP messages
cannot be obtained, we present suitable approximations that
lead us to propose the AC-BP algorithm as an (approximate)
BP solution for the AC model. Then, as a main contribution,
we make a key further step where we change the perspective
of our BP approach and, instead of applying the BP directly
onto the non-linear AC model, we present the solution where
the BP is applied sequentially over the AC model, akin to
what is done by the Gauss-Newton method. The resulting
Gauss-Newton BP (GN-BP) represents a BP counterpart of
the Gauss-Newton method achieving the same accuracy,
however, preserving a number of advantages brought in by
the BP framework. For example, the proposed GN-BP
algorithm is of considerably lower computational complexity,
can be easily designed to provide asynchronous operation,
could be integrated as part of real-time systems [28], and is
flexible and easy to distribute and parallelize. Finally, this
paper provides a novel and detailed convergence analysis of
the DC-BP algorithm and points to extension of this analysis
for the proposed GN-BP algorithm. We close this paper by
extensive numerical study where we compare the proposed
BP-based algorithms with the centralized SE, and provide a
number of useful insights for their implementation.
Paper Organization: In Section II, we present basic
concepts of factor graphs and BP. Section III reviews the
centralized DC and AC SE. Section IV presents closed form
BP expressions and detailed convergence analysis for the DC
model. In Section V, the AC model is solved via BP-based

Gauss-Newton method. Section VI presents extended
numerical convergence and performance results. Concluding
remarks are given in Section VII.

II. FACTOR G RAPHS AND BP A LGORITHM
Factor graphs and BP algorithm are widely used tools for
probabilistic inference [21], [22]. In the standard setup, the
goal of the BP algorithm is to efficiently evaluate the marginals
of a system of random variables x = (x1 , . . . , xn ) described
via the joint probability density function g(x). Assuming that
the function g(x) can be factorized proportionally (∝) to a
product of local functions:
g(x) ∝

k
Y

ψi (Xi ),

(1)

i=1

where Xi ⊆ {x1 , . . . , xn }, the marginalization problem can
be efficiently solved using BP algorithm.
The first step is forming a factor graph, which is a
bipartite graph that describes the structure of the
factorization (1). The factor graph structure comprises the set
of factor nodes F = {f1 , . . . , fk }, where each factor node fi
represents local function ψi (Xi ), and the set of variable
nodes X = {x1 , . . . , xn }. The factor node fi connects to the
variable node xs if and only if xs ∈ Xi [29].
The BP algorithm on factor graphs proceeds by passing
two types of messages along the edges of the factor graph: i)
a variable node to a factor node, and ii) a factor node to a
variable node messages. Both variable and factor nodes in a
factor graph process the incoming messages and calculate
outgoing messages, where an output message on any edge
depends on incoming messages from all other edges. BP
messages represent ”beliefs” about variable nodes, thus a
message that arrives or departs a certain variable node is a
function (distribution) of the random variable corresponding
to the variable node.
Message from a variable node to a factor node: Consider
a part of a factor graph shown in Fig. 1 with a group of factor
nodes Fs = {fi , fw , ..., fW } ⊆ F that are neighbours of the
variable node xs ∈ X . The message µxs →fi (xs ) from the
µ f w → xs ( x s )

fw

.
.
.
fW

x s µ xs → f i ( x s )

fi

µ f W → xs ( x s )

Fig. 1. Message µxs →fi (xs ) from variable node xs to factor node fi

variable node xs to the factor node fi is equal to the product
of all incoming factor node to variable node messages arriving
at all the other incident edges:
Y
µxs →fi (xs ) =
µfa →xs (xs ),
(2)
fa ∈Fs \fi

3

where Fs \ fi represents the set of factor nodes incident to the
variable node xs , excluding the factor node fi . Note that each
message is a function of the variable xs .
Message from a factor node to a variable node:
Consider a part of a factor graph shown in Fig. 2 that
consists of a group of variable nodes Xi = {xs , xl , ..., xL } ⊆
X that are neighbours of the factor node fi ∈ F. The
µ xl → f i ( x l )

xl

.
.
.

f i µ f i → xs ( x s )

xs

µ xL → f i ( x L )

xL

Fig. 2. Message µfi →xs (xs ) from factor node fi to variable node xs

message µfi →xs (xs ) from the factor node fi to the variable
node xs is defined as a product of all incoming variable
node to factor node messages arriving at other incident
edges, multiplied by the function ψi (Xi ) associated to the
factor node fi , and marginalized over all of the variables
associated with the incoming messages:
Z

Z
···

µfi →xs (xs ) =

xL

xl

Y

ψi (Xi )




µxb →fi (xb ) · dxb , (3)

xb ∈Xi \xs

where Xi \xs is the set of variable nodes incident to the factor
node fi , excluding the variable node xs .
Marginal inference: The marginal of the variable node xs ,
illustrated in Fig. 3, is obtained as the product of all incoming
messages into the variable node xs :
Y
p(xs ) =
µfc →xs (xs ),
(4)
fc ∈Fs

where Fs is the set of factor nodes incident to the variable
node xs .
µ f w → xs ( x s )

fw

.
.
.
fW

x s µ f i → xs ( x s )

fi

III. SE IN E LECTRIC P OWER S YSTEMS
The SE algorithm estimates the values of the state
variables based on the knowledge of network topology and
parameters, and measured values obtained from measurement
devices spread across the power system. The knowledge of
the network topology and parameters is provided by the
network topology processor in the form of the bus/branch
model, where branches of the grid are usually described
using the two-port π-model [1, Ch. 1,2]. The bus/branch
model can be represented using a graph G = (V, E), where
the set of nodes V = {1, . . . , N } represents the set of buses,
while the set of edges E ⊆ V × V represents the set of
branches of the power network.
As an input, the SE requires a set of measurements M of
different electrical quantities spread across the power
network. Using the bus/branch model and available
measurements, the observability analysis defines the
measurement model [1, Ch. 4] which can be described as the
system of equations [4]:
z = h(x) + u,

where x = (x1 , . . . , xn ) is the vector of the state variables,
h(x) = (h1 (x), . . . , hk (x)) is the vector of measurement
functions, z = (z1 , . . . , zk ) is the vector of measurement
values, and u = (u1 , . . . , uk ) is the vector of uncorrelated
measurement errors. The SE problem is commonly an
overdetermined system of equations (k > n) [31, Sec. 2.1].
Each measurement Mi ∈ M is associated with measured
value zi , measurement error ui , and measurement function
hi (x). Under the assumption that measurement errors ui
follow a zero-mean Gaussian distribution, the probability
density function associated with the i-th measurement equals:
(
)
[zi − hi (x)]2
1
exp
,
(6)
N (zi |x, vi ) = √
2vi
2πvi
where vi is the measurement variance defined by the
measurement error ui , and the measurement function hi (x)
connects the vector of state variables x to the value of the
i-th measurement.
The SE in electric power systems deals with the problem of
determining state variables x according to the noisy observed
data z and a prior knowledge:

µ f W → xs ( x s )

Fig. 3. Marginal inference of the variable node xs

In the following, we focus on the Gaussian BP (GBP)
over the set of real-valued random variables where all the
local functions ψi (Xi ) and all the inputs to the BP algorithm
represent Gaussian distributions. Due to the fact that variable
node and factor node processing preserves ”Gaussianity” of
the messages, the resulting GBP algorithm operates by
exchanging the messages that represent Gaussian functions.
Therefore, each message exchanged in GBP is completely
represented using only two values: the mean and the
variance [30].

(5)

p(x|z) =

p(z|x)p(x)
.
p(z)

(7)

Assuming that the prior probability distribution p(x) is
uniform, and given that p(z) does not depend on x, the
maximum a posteriori (MAP) solution of (7) reduces to the
maximum likelihood solution, as given below [32]:
x̂ = arg max p(x|z) = arg max p(z|x) = arg max L(z|x). (8)
x
x
x

One can find the solution (8) via maximization of the
likelihood function L(z|x), which is defined via likelihoods
of k independent measurements:
x̂ = arg max L(z|x) = arg max
x

x

k
Y
i=1

N (zi |x, vi ).

(9)

4

It can be shown that the solution of the MAP problem can
be obtained by solving the following optimization problem,
known as the WLS problem [33, Sec. 9.3]:
x̂ = arg min
x

k
X
[zi − hi (x)]2
i=1

vi

.

(10)

The state estimate x̂ representing the solution of the
optimization problem (10) is known as the WLS estimator,
and is equivalent to the MAP estimator (9).
A. AC State Estimation

only bus voltage angles x ≡ θ as state variables.
Consequently, the number of state variables is n = N .
The set of DC model measurements M involves only active
power flow {MPij }, (i, j) ∈ E, active power injection {MPi },
i ∈ V, and bus voltage angle {Mθi }, i ∈ V, measurements and
the model is dealing with linear measurement functions h(x)
given Appendix A.
The DC state estimate x̂ ≡ θ̂, which is a solution to the
WLS problem (10), is obtained through the non-iterative
procedure by solving the system of linear equations:


HT WH x̂ = HT Wz,
(12)

The AC SE model is defined using the measurement
functions h(x) that precisely follow the physical laws that
connect the measured variables and the state variables. As a
result, the system (5) in general represents the system of
non-linear equations.
In a usual scenario, the AC SE model takes bus voltage
magnitudes and bus voltage angles, transformer magnitudes
of turns ratio and transformer angles of turns ratio as state
variables x. Without loss of generality, in the rest of the
paper, for the AC SE we observe bus voltage angles θ =
(θ1 , . . . , θN ) and bus voltage magnitudes V = (V1 , . . . , VN )
as state variables x ≡ (θ, V). Consequently, the number of
state variables is n = 2N .
The typical set of measurements M in electric power
systems includes: active and reactive power flow and current
magnitude, {MPij , MQij , MIij }, (i, j) ∈ E, respectively;
active and reactive power injection, {MPi , MQi }, i ∈ V; and
bus voltage angle and magnitude, {Mθi , MVi }, i ∈ V. For
completeness, we provide the set of non-linear measurement
functions h(x) of the AC SE model given in Appendix A,
equations (62) - (65).
Based on the available set of measurements, the WLS
estimator x̂ ≡ (θ̂, V̂), i.e., the solution of the WLS problem
(10), can be found using the Gauss-Newton method:
h
i
J(x(ν) )T WJ(x(ν) ) ∆x(ν) = J(x(ν) )T Wr(x(ν) ) (11a)

where H ∈ RkxN is the Jacobian matrix of measurement
functions given in Appendix A, equations (66) - (68).

x(ν+1) = x(ν) + ∆x(ν) ,

B. DC State Estimation

For the DC model, the set of variable nodes is defined by
the state variables θ, i.e., X = {θ1 , . . . , θN }. The set of
factor nodes F = {f1 , . . . , fk } is defined by the set of
measurements M, with measurement functions given in
Appendix A, equations (66) - (68). Measurements define
likelihood functions N (zi |x, vi ) that are in turn equal to
local functions ψi (Xi ) associated to factor nodes. More
precisely, a factor node fi connects to a variable node
xs ∈ X if and only if the state variable xs is an argument of
the corresponding measurement function hi (x).

The DC model is obtained by linearisation of the AC
model. In the typical operating conditions, the difference of
bus voltage angles between adjacent buses (i, j) ∈ E is very
small θi − θj ≈ 0, which implies cos θij ≈ 1 and
sin θij ≈ θij . Further, all bus voltage magnitudes are Vi ≈ 1,
i ∈ V, and all shunt elements and branch resistances can be
neglected. This implies that the DC model ignores the
reactive powers and transmission losses and takes into
account only the active powers. Therefore, the DC SE takes

Example 1 (Constructing factor graph). In this toy example,
using a simple 3-bus model presented in Fig. 4a, we
demonstrate the conversion from a bus/branch model with a
given measurement configuration into the corresponding
factor graph for the DC model. The variable nodes represent
state variables, i.e., X = {θ1 , θ2 , θ3 }. Factor nodes are
defined by corresponding measurements, where in our
example, measurements MP12 and MP3 are mapped into
factor nodes F = {fP12 , fP3 }.
4

(11b)

where ν = {0, 1, 2, . . . } is the iteration index, ∆x(ν) ∈ Rn is
the vector of increments of the state variables, J(x(ν) ) ∈ Rkxn
is the Jacobian matrix of measurement functions h(x(ν) ) at
x = x(ν) (see Appendix A, for details), W ∈ Rkxk is a
diagonal matrix containing inverses of measurement variances,
and r(x(ν) ) = z−h(x(ν) ) is the vector of residuals [3, Ch. 10].

IV. BP- BASED A LGORITHM FOR THE DC SE
In the SE problem defined above, each measurement
function hi (x) depends on a limited (typically small) subset
of state variables x. Hence, the likelihood function L(z|x)
can be factorized into factors (9) affecting small subsets of
state variables. This motivates solving the SE problem
scalably and efficiently using probabilistic graphical models.
The solution involves defining the factor graph
corresponding to (9), and deriving expressions for BP
messages exchanged over the factor graph, as detailed next.
For completeness of exposition, we first present the
solution of the DC SE problem using the BP algorithm; we
refer to the corresponding method as the DC-BP.
Furthermore, and as a novel contribution, we provide an
in-depth convergence analysis of the DC-BP algorithm, and
present methods to improve its convergence. The material in
this section sets the stage for the main contribution of this
paper - the BP-based Gauss-Newton method for AC model
presented in the following section.
A. The Factor Graph Construction

5

bus 1
θ1

2

1 MP12

fP12

The DC model contains only linear measurement functions
which we represent in a general form as:
X
hi (xs , xl , . . . , xL ) = Cxs xs +
Cxb xb ,
(18)

bus 2
θ2

xb ∈Xi \xs

f P3

3
MP 3
θ3
bus 3

(a)

(b)
Fig. 4.
Transformation of the bus/branch model and measurement
configuration (subfigure a) into the corresponding factor graph for the DC
model (subfigure b).

B. Derivation of BP Messages

where Xi \xs is the set of variable nodes incident to the factor
node fi , excluding the variable node xs . From the expression
(3), and using (16)-(18), it can be shown that the message
µfi →xs (xs ) from the factor node fi to the variable node xs is
represented by the Gaussian function (15), with mean zfi →xs
and variance vfi →xs obtained as:
!
X
1
zi −
Cxb zxb →fi
(19a)
zfi →xs =
Cxs
xb ∈Xi \xs
!
X
1
2
vfi →xs = 2 vi +
Cxb vxb →fi .
(19b)
Cxs
xb ∈Xi \xs

Message from a variable node to a factor node: Let us
assume that the incoming messages µfw →xs (xs ), . . . ,
µfW →xs (xs ) into the variable node xs are Gaussian and
represented by their mean-variance pairs (zfw →xs , vfw →xs ),
. . . , (zfW →xs , vfW →xs ). According to (2), it can be shown
that the message µxs →fi (xs ) from the variable node xs to
the factor node fi is proportional (∝) to:

To summarize, after the factor node fi receives the messages
from all of the neighbouring variable nodes from the set Xi \
xs , it evaluates the message µfi →xs (xs ) according to (19a)
and (19b), and sends it to the variable node xs .
Marginal inference: According to (4), it can be shown that
the marginal of the state variable xs is represented by:

µxs →fi (xs ) ∝ N (zxs →fi |xs , vxs →fi ),

p(xs ) ∝ N (x̂s |xs , vxs ),

with mean zxs →fi and variance vxs →fi obtained as:
!
X zfa →xs
vxs →fi
zxs →fi =
vfa →xs

(13)

(14a)

fa ∈Fs \fi

1
vxs →fi

X

=

fa ∈Fs \fi

1
.
vfa →xs

(14b)

(15)

(16)

µxL →fi (xL ) ∝ N (zxL →fi |xL , vxL →fi ).
The Gaussian function associated with the factor node fi is
given by (Section III, equation (6)):
N (zi |xs , xl , . . . , xL , vi )
(
∝ exp

X
1
1
=
.
vx s
vfc →xs

(21b)

fc ∈Fs

Let us assume that the messages into factor nodes are
Gaussian, denoted by:
µxl →fi (xl ) ∝ N (zxl →fi |xl , vxl →fi )
..
.

(21a)

fc ∈Fs

After the variable node xs receives the messages from all of
the neighbouring factor nodes from the set Fs \ fi , it evaluates
the message µxs →fi (xs ) according to (14) and sends it to the
factor node fi .
Message from a factor node to a variable node: Due to
linearity of measurement functions hi (·), closed form
expressions for these messages is easy to obtain and follow a
Gaussian form:
µfi →xs (xs ) ∝ N (zfi →xs |xs , vfi →xs ).

with the mean value x̂s and variance vxs :
!
X zfc →xs
vxs
x̂s =
vfc →xs

(20)

)
[zi − hi (xs , xl , . . . , xL )]2
. (17)
2vi

Finally, the mean-value x̂s is adopted as the estimated value
of the state variable xs .
C. Iterative DC-BP Algorithm
The SE scenario is in general an instance of loopy BP
since the corresponding factor graph usually contains cycles.
Loopy BP is an iterative algorithm, with an iteration index
τ = {0, 1, 2, . . . }, and requires a message-passing schedule.
The scheduling where messages from variable to factor
nodes, and messages from factor nodes to variable nodes, are
updated in parallel in respective half-iterations, is known as
synchronous scheduling. Synchronous scheduling updates all
messages in a given iteration using the output of the
previous iteration as an input [34].
To present the algorithm precisely, we need to introduce
different types of factor nodes. The indirect factor nodes
Find ⊂ F correspond to measurements that measure state
variables indirectly. In the DC scenario, this includes active
power flow and power injection measurements. The direct
factor nodes Fdir ⊂ F correspond to the measurements that
measure state variables directly. For our choice of state
variables for the DC scenario, an example includes
measurements of angles.

6

Besides direct and indirect factor nodes, we define three
additional types of singly-connected factor nodes. The slack
factor node corresponds to the slack or reference bus where
the voltage angle has a given value. The initialization factor
node is required to initialize the algorithm at certain variable
nodes. Finally, the virtual factor node is a is singly-connected
factor node used if the variable node is not directly measured.
Both initialization and virtual factor nodes take the value of
”flat start” with variance σx2i → ∞ or a priori given mean
value and variance of state variables. We refer to direct factor
nodes and three additional types of singly-connected factor
nodes1 as local factor nodes Floc ⊂ F.
Algorithm 1 The DC-BP
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:

procedure I NITIALIZATION τ = 0
for each fs ∈ Floc do
(0)
send µfs →xs to incident xs ∈ X
end for
for each xs ∈ X do
(0)
(0)
send µxs →f = µfs →xs , to incident fi ∈ Find
i
end for
end procedure
procedure I TERATION LOOP τ = 1, 2, . . .
while stopping criterion is not met do
for each fi ∈ Find do
(τ )
Compute µf →xs using (19a)*, (19b)*
i
end for
for each xs ∈ X do
(τ )
Compute µxs →f using (14)
i
end for
end while
end procedure
procedure O UTPUT
for each xs ∈ X do
Compute x̂s , vxs using (21)
end for
end procedure
*Incoming messages are obtained in previous iteration τ − 1

Example 2 (Different types of factor nodes). In this example,
we consider the bus/branch model with four measurements
illustrated in Fig. 5a that we use to describe different types
of factor nodes. The corresponding factor graph is given in
Fig. 5b, with indirect factor nodes (red squares), direct factor
nodes (orange squares), the initialization factor node (green
square) and the virtual factor node (blue square).
4
The DC-BP algorithm is presented in Algorithm 1, and it
is an instance of a loopy GBP applied over a linear model
defined by linear measurement functions h(x). It is well
known that, if loopy GBP applied over a linear model
converges, it will converge to a fixed point representing a
solution of an equivalent WLS problem (12) [35]. Unlike
means, the variances of GBP messages need not converge to
1 We note that each variable node in the initialization step of the algorithm
has the corresponding local factor node attached. Local factor nodes only send,
but do not receive, the message to/from incident variable nodes. Direct and
virtual factor nodes repeatedly transmit the same message to the corresponding
variable node throughout BP iterations.

f θ1

1 MP12

2

fθ2
fP12

f θ3
fP23

3

MP23

θ1

θ2
bus 2
fP4

bus 1

θ3
bus 3

Mθ 3

4

f θ4

MP 4

(a)

θ4
bus 4

(b)
Fig. 5.
Transformation of the bus/branch model and measurement
configuration (subfigure a) into the corresponding factor graph with different
types of factor nodes for the DC model (subfigure b)

correct values. In the following, it is sufficient to investigate
the conditions for the DC-BP convergence as in case of
convergence, the obtained mean values are the WLS solution.
D. Convergence of DC-BP Algorithm
In this part, we present convergence analysis of DC-BP
algorithm with synchronous scheduling. In the following, it
will be useful to consider a subgraph of the factor graph that
contains the set of variable nodes X = {x1 , . . . , xN }, the set
of indirect factor nodes Find = {f1 , . . . , fm } ⊂ F, and a set
of edges B ⊆ X × Find connecting them. The number of
edges in this subgraph is b = |B|. Within the subgraph, we
will consider a factor node fi ∈ Find connected to its
neighboring set of variable nodes Xi = {xq , . . . , xQ } ⊂ X
by a set of edges Bi = {bqi , . . . , bQ
i } ⊂ B, where di = |Xi | is
the degree of fi . Next, we provide results on convergence of
both variances and means of DC-BP messages, respectively.
Convergence of the Variances: From equations (14b) and
(19b), we note that the evolution of the variances is
independent of mean values of messages and measurements.
Let vs ∈ Rb denote a vector of variance values of messages
from indirect factor nodes Find to variable nodes X . Note
that this vector can be decomposed as:
(τ )

(τ ) T
vs(τ ) = [vs,1 , . . . vs,m
] ,

(22)

where the i-th element vs,i ∈ Rdi is equal to:
(τ )

(τ )

(τ )

vs,i = [vfi →xq , . . . , vfi →xQ ].

(23)

Substituting (14b) in (19b), the evolution of variances vs is
equivalent to the following iterative equation:

(24)
vs(τ ) = f vs(τ −1) .
More precisely, using simple matrix algebra, one can obtain
the evolution of the variances vs in the following matrix form:
h
i


e −1 ΠC
e · D(A) −1 + Σa C
e −1 i
vs(τ ) = C
(25)
where
e = CCT
C
A=

T
ΓΣ−1
s Γ

(26a)
+ L.

(26b)

7

(τ −1)

Note that in (25), the dependance on vs
is hidden in matrix
A, or more precisely, in matrix Σs . Next, we briefly describe
both the matrices and matrix-operators involved in (25).
The operator D(A) ≡ diag(A11 , . . . , Abb ), where Aii is
the i-th diagonal entry of the matrix A. The unit vector i is
of dimension b and is equal to i = [1, . . . , 1]T . The diagonal
(τ −1)
matrix Σs is obtained as Σs = diag vs 
∈ Rb×b .
The matrix C = diag C1 , . . . , Cm ∈ Rb×b contains
diagonal entries of the Jacobian non-zero elements, where
i-th element Ci = [Cxq, . . . , CxQ ] ∈ Rdi . The matrix
Σa = diag Σa,1 , · · · Σa,m ∈ Rb×b contains indirect factor
node variances, with the i-th entry Σ
[vi , . . . , vi ] ∈ Rdi .
 a,i =b×b
The matrix L = diag L1 , · · · Lm ∈ R
contains inverse
variances from singly-connected factor nodes to a variable
node,  if such nodes exist, where the i-th element
Li = lxq , · · · lxQ ∈ Rdi . For example, lxq equals:


1

, if xq is incident to fd,q
(27)
lxq = vfd,q →xq

0,
otherwise.

The matrix Π = diag Π1 , . . . Πm ∈ Fb×b
2 , F2 = {0, 1},
is a block-diagonal matrix in which the i-th element is a block
matrix Πi = 1i − Ii ∈ Fd2i ×di , where the matrix 1i is di × di
block matrix of ones, and Ii is di × di identity matrix. The
matrix Γ ∈ Fb×b
is of the following block structure:
2


01,1 Γ1,2 . . . Γ1,m
 Γ2,1 02,2 . . . Γ2,m 


Γ= .
(28)
,
..
..
 ..

.
.
Γm,1 Γm,2 . . . 0m,m
d ×d

where 0i,i is a block matrix di ×di of zeros, and Γi,j ∈ F2i j
with the (i, j)-th entry:
(
1, if both bqi and bqj are incident to xq
Γi,j (i, j) =
(29)
0, otherwise.
Note that the following holds: Γj,i =

ΓT
i,j .

Theorem 1. The variances vs from indirect factor nodes to
variable nodes always converge to a unique fixed point
(τ )
(τ =0)
limτ →∞ vs = v̂s for any initial point vs
> 0.

Proof. The theorem can be proved by showing that f vs
satisfies the conditions of the so-called standard function [36],
following similar steps as in the proof of Lemma 1 in [37].

(31d)

Note that the vector of means zs ∈ Rb can be decomposed as:
(τ )

)
) T
z(τ
= [zs,1 , . . . , z(τ
s
s,m ] ,

(32)

where the i-th element zs,i ∈ Rdi is equal to:
(τ )
(τ )
(τ )
(33)
zs,i = [zfi →xk , . . . , zfi →xK ].

T
The vector za = za,1 , · · · za,m ∈ Rb contains means of
indirect factor nodes, where za,i = [zi , . . . , zi ] ∈ Rdi . The
diagonal matrix Σ̂s ∈ Rb×b is obtained as  Σ̂s =
(τ )
limτ →∞ Σs . The vector zb = zb,1 , · · · , zb,m ∈ Rb
contains means from direct and virtual factor nodes to a
variable node, if such
 nodes exist, where the i-th element
zb,i = zxk , · · · zxK ∈ Rdi . For example, the element zxk of
zb,i is equal to:
(
zfd,k →xk , if xk is incident to fd,k
zxk =
(34)
0,
otherwise.

Theorem 2. The means zs from indirect factor nodes to
variable nodes converge to a unique fixed point
(τ )
limτ →∞ zs = ẑs :
−1
(35)
e
ẑs = I + Ω
z,
(τ =0)

for any initial point zs
ρ(Ω) < 1.

if and only if the spectral radius

Proof. The proof steps follow the proof of Theorem 5.2, [36].
To summarize, the convergence of the DC-BP algorithm
depends on the spectral radius of the matrix:
−1

 
T
Ω = C−1 ΠC · D(ΓΣ̂−1
· ΓΣ̂−1
. (36)
s Γ + L)
s
If the spectral radius ρ(Ω) < 1, the DC-BP algorithm will
converge and the resulting vector of mean values will be equal
to the solution of the MAP estimator.
E. Convergence of DC-BP with Randomized Damping

(31a)
(31b)

)
(τ −1)
)
zd = z(τ
+ α2 z(τ
q + α1 zr
r ,

(30)

where
−1
e
z = C−1 za − D · D(Â)
· Lzb
−1
−1
Ω = D · D(Â)
· ΓΣ̂s

(31c)

D = C−1 ΠC.

In this section, we propose an improved DC-BP algorithm
that applies synchronous scheduling with randomized
damping. Several previous works reported that damping the
BP messages improves the convergence of BP [37], [38].
Here, we propose a different randomized damping approach,
where each mean value message from indirect factor node to
a variable node is damped independently with probability p,
otherwise, the message is calculated as in the standard
DC-BP algorithm. The damped message is evaluated as a
linear combination of the message from the previous and the
current iteration, with weights α1 and 1 − α1 , respectively.
In Section VI, we demonstrate that the DC-BP with
randomized damping dramatically improves convergence as
compared to the standard DC-BP.
In the proposed damping, the equation (30) is redefined as:

Convergence of the Means: Equations (14a) and (19a)
show that the evolution of the mean values depends on the
variance values. Due to Theorem 1, it is possible to simplify
evaluation of mean values zs from indirect factor nodes Find
to variable nodes X by using the fixed-point values of v̂s .
The evolution of means zs becomes a set of linear equations:
)
−1)
z(τ
=e
z − Ωz(τ
,
s
s

T
Â = ΓΣ̂−1
s Γ +L

(τ )

(37)

8

where 0 < α1 < 1 is the weighting coefficient, and α2 =
(τ )
(τ )
1 − α1 . In the above expression, zq and zr are obtained
as:
)
−1)
z(τ
z − QΩz(τ
q = Qe
s
)
z(τ
r

= Re
z−

−1)
RΩz(τ
,
s

(38a)
(38b)

where diagonal matrices Q ∈ Fb×b
and R ∈ Fb×b
are
2
2
defined as Q = diag(1 − q1 , ..., 1 − qb ), qi ∼ Ber(p), and
R = diag(q1 , ..., qb ), respectively, where Ber(p) ∈ {0, 1} is
a Bernoulli random variable with probability p independently
sampled for each mean value message.
Substituting (38a) and (38b) in (37), we obtain:

 −1)
(τ )
zd = Q + α2 R e
z − QΩ + α2 RΩ − α1 R z(τ
.
s
(39)
(τ −1)
(τ −1)
Note that zr
= Rzs
. In a more compact form,
equation (39) can be written as follows:
(τ )

−1)
zd = z̄ − Ω̄z(τ
,
s

(40)

where

z̄ = Q + α2 R e
z

(41a)

Ω̄ = QΩ + α2 RΩ − α1 R.

(41b)

Theorem 3. The means zd from indirect factor nodes to
variable nodes converge to a unique fixed point
(τ )
(τ =0)
ẑd = limτ →∞ zd for any initial point zd
if and only if
the spectral radius ρ(Ω̄) < 1. For the resulting fixed point, it
holds that ẑd = ẑs .
Proof. To prove theorem it is sufficient to show that equation
(40) converges to the fixed point defined in (35). We can write:
zr

(τ −1)

= Re
z−

−2)
RΩz(τ
.
s

(42)

Substituting (38a), (38b) and (42) in (37):

(τ )
z
zd = Q + α2 R + α1 R e
 −1)
− QΩ + α2 RΩ z(τ
− α1 RΩzs(τ −2) . (43)
s
(τ )

The fixed point ẑd = limτ →∞ zd

is equal to:

−1
ẑd = I + QΩ + α2 RΩ + α1 RΩ

· Q + α2 R + α1 R e
z. (44)
From definitions of Q, R and α2 , we have QΩ + α2 RΩ +
α1 RΩ = Ω and Q + α2 R + α1 R = I, thus (44) becomes:
−1
(45)
e
ẑd = I + Ω
z.
This concludes the proof.
V. BP-BASED D ISTRIBUTED G AUSS -N EWTON M ETHOD
FOR AC SE
The rest of the paper focuses on the more challenging AC
model, where our goal is to derive the BP-based AC SE
algorithm. However, due to non-linearity of measurement
functions, the closed-form expressions for certain classes of
BP messages cannot be obtained. Using approximations, we

derive the AC-BP algorithm as an (approximate) BP solution
for the AC SE problem. Due to space constraints, we refer
the reader to Appendix C, for a detailed derivation of the
AC-BP algorithm. Unfortunately, due to approximations, the
AC-BP algorithm does not match the performance of the
centralized AC SE based on Gauss-Newton method (as we
show in Section VI).
As the main contribution of this paper, we adopt different
methodology to derive efficient BP-based AC SE method.
Namely, instead of using the AC-BP method, we present the
solution where BP is applied sequentially over the AC
model, akin to what is done by the Gauss-Newton method.
The resulting GN-BP algorithm represents a BP counterpart
of the Gauss-Newton method, achieving exactly the same
accuracy as the Gauss Newton method.

A. Gauss-Newton Method as a Sequential MAP Problem
Consider the Gauss-Newton method (11) where, at each
iteration step ν, the algorithm returns a new estimate of x
denoted as x(ν) . Note that, after a given iteration, an
estimate x(ν) is a vector of known (constant) values. If the
Jacobian matrix J(x(ν) ) has a full column rank, the equation
(11a) represents the linear WLS solution of the minimization
problem [39]:
min ||W1/2 [r(x(ν) ) − J(x(ν) )∆x(ν) ]||22 .

∆x(ν)

(46)

Hence, at each iteration ν, the Gauss-Newton method produces
WLS solution of the following system of linear equations:
r(x(ν) ) = g(∆x(ν) ) + u,

(47)

where g(∆x(ν) ) = J(x(ν) )∆x(ν) comprises linear functions,
while u is the vector of measurement errors. The equation
(11a) is the weighted normal equation for the minimization
problem defined in (46), or alternatively, equation (11a) is a
WLS solution of (47). Consequently, the probability density
function associated with the i-th measurement (i.e., the i-th
residual component ri ) at any iteration step ν:
N (ri (x(ν) )|∆x(ν) , vi )
1
exp
=√
2πvi

(

)
[ri (x(ν) ) − gi (∆x(ν) )]2
. (48)
2vi

The MAP solution of (9) can be redefined as an iterative
optimization problem where, instead of solving (11), we solve
the MAP (sub)problem:


∆x̂(ν) = arg max L r(x(ν) )|∆x(ν)
∆x(ν)

= arg max

∆x(ν)

k
Y



N ri (x(ν) )|∆x(ν) , vi

i=1
(ν)

x(ν+1) = x(ν) + ∆x̂

.

(49a)
(49b)

Next, we show that the solution of the above MAP
sub-problem (49a) over increment variables ∆x(ν) can be
efficiently obtained using the BP algorithm applied over the
underlying factor graph.

9

Note that, in general, the resulting BP algorithm represents
an instance of a loopy GBP over a linear model defined by
linear functions g(∆x(ν) ). Thus if it converges, it provides a
solution equal to the linear WLS solution ∆x(ν) of (11a).
The BP solution of ∆x̂(ν) in each iteration ν (outer iteration
loop) is obtained via iterative BP algorithm (inner iteration
loops). Every inner BP iteration τ = {0, 1, 2, . . . η(ν)}
outputs ∆x̂(ν,τ ) , where η(ν) is the number of inner BP
iterations within the outer iteration ν.

Marginal inference: It can be shown that the marginal of
the state variable ∆xs , according to (4), is represented by the
Gaussian function:
p(∆xs ) ∝ N (∆x̂s |∆xs , v∆xs ),

(55)

with mean ∆x̂s which represents the estimated value of the
state variable increment ∆xs and variance v∆xs :
!
X rfc →∆xs
v∆xs
(56a)
∆x̂s =
vfc →∆xs
fc ∈Fs

1

B. The Factor Graph Construction
From the factorization of the likelihood expression (49a),
one easily obtains the factor graph corresponding to the
GN-BP method as follows. The increments ∆x of state
variables x determine the set of variable nodes
X = {(∆θ1 , ∆V1 ), . . . , (∆θN , ∆VN )} and each likelihood
function N (ri (xν )|∆xν , vi ) represents the local function
ψ(Xi ). Since the residual equals ri (xν ) = zi − hi (xν ), in
general, the set of factor nodes F = {f1 , . . . , fk } is defined
by the set of measurements M. The factor node fi connects
to the variable node ∆xs ∈ {∆θs , ∆Vs } iff the increment of
the state variable ∆xs is an argument of the corresponding
function gi (∆x), i.e., if the state variable xs ∈ {θs , Vs } is an
argument of the measurement function hi (x).

C. Derivation of BP Messages
Message from a variable node to a factor node: It is easy
to show that the message µ∆xs →fi (∆xs ) is proportional to:
µ∆xs →fi (∆xs ) ∝ N (r∆xs →fi |∆xs , v∆xs →fi ),

(50)

with mean r∆xs →fi and variance v∆xs →fi obtained as:
!
X rfa →∆xs
v∆xs →fi
(51a)
r∆xs →fi =
vfa →∆xs
fa ∈Fs \fi

1
v∆xs →fi

1

X

=

fa ∈Fs \fi

vfa →∆xs

.

(51b)

Message from a factor node to a variable node: Similarly,
the message µfi →∆xs (∆xs ) is proportional to:
µfi →∆xs (∆xs ) ∝ N (rfi →∆xs |∆xs , vfi →∆xs ),

(52)

with mean rfi →∆xs and variance vfi →∆xs obtained as:
rfi →∆xs =

vfi →∆xs =

1
C∆xs
1
2
C∆x
s

!
ri −

X

C∆xb · r∆xb →fi

(53a)

∆xb ∈Xi \∆xs

!
vi +

X

2
C∆x
b

· v∆xb →fi .

(53b)

∆xb ∈Xi \∆xs

The coefficients C∆xp , ∆xp ∈ Xi , are Jacobian elements
of the measurement function associated with the factor node
fi (see Appendix A and D for details):
C∆xp =

∂h(xs , xl , . . . , xL )
.
∂xp

(54)

v∆xs

=

1

X
fc ∈Fs

vfc →∆xs

.

(56b)

To summarize, the MAP sub-problem defined in (49a) can
be efficiently solved using (51), (53) and (56).
D. Iterative GN-BP Algorithm
To introduce the GN-BP algorithm, we use the same
classification of factor nodes and the same message-passing
schedule as in Section IV-C for the DC-BP algorithm. The
slack factor node is specified as a node where the voltage
angle has a given value, therefore, the residual of the
corresponding state variable is equal to zero, and its variance
tends to zero. Residuals of initialization and virtual factor
nodes approach zero, while their variances tend to infinity. It
is important to note that local factor nodes send messages
represented by a triplet: mean (of the residual), variance and
the state variable value (see Appendix D, for details).
The GN-BP algorithm is presented in Algorithm 2. After
the initialization, the outer loop starts by computing residuals
for direct and indirect factor nodes, as well as the Jacobian
elements, and passes them to the inner iteration loop. The
inner iteration loop represents the main algorithm routine
which includes BP-based message inference described in the
previous subsection. The output of the inner iteration loop is
the estimate of the state variable increments. Finally, the
outer loop updates the set of state variables. The outer loop
iterations are repeated until the stopping criteria is met.
E. Convergence of the GN-BP Algorithm
The GN-BP in the inner iteration loop solves the system
of linear equations (47). By directly applying the results from
Section IV-D of the (linear) DC-BP convergence analysis, the
convergence of the GN-BP with synchronous scheduling in
each outer iteration loop ν depends on the spectral radius of
the matrix:


Ω(x(ν) ) = C(x(ν) )−1 ΠC(x(ν) )
−1


T
· D(ΓΣ̂−1
· ΓΣ̂−1
, (57)
s Γ + L)
s
where C(x(ν) ) contains diagonal entries of the Jacobian
J(x(ν) ) non-zero elements in every outer iteration ν.
According to Theorem 2, the GN-BP with synchronous
scheduling converges to a unique fixed point if
ρ(Ω(x(ν) )) < 1, ν = {0, 1, 2, . . . }. We note that this

10

Algorithm 2 The GN-BP
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:

Although similar results can be obtained for the GN-BP
algorithm (see Section V-E), for brevity, we present the
convergence results only for the DC-BP case.
For evaluation of accuracy, we focus on the proposed BPbased AC SE methods and compare them to the Gauss-Newton
method (the centralized SE algorithm). More precisely, we
first consider the IEEE 14 bus test case and provide a detailed
performance comparison between the AC-BP and the GN-BP
algorithms. Then, due to favorable performance, we provide
an in-depth evaluation of the GN-BP algorithm using the IEEE
14 bus test case, providing the results for IEEE 30 and IEEE
118 bus test cases in Appendix F.

procedure I NITIALIZATION ν = 0
for Each xs ∈ X do
(0)
Initialize value of xs , typically as a ”flat start”
end for
end procedure
procedure O UTER ITERATION LOOP ν = 0, 1, 2, . . . ; τ = 0
while stopping criterion for the outer loop is not met do
for Each fs ∈ Fdir do
(ν)
(ν)
compute residual rs = zs − xs
end for
for Each fs ∈ Floc do
(ν)
(ν)
send µfs →∆xs , xs to incident ∆xs ∈ X
end for
for Each ∆xs ∈ X do
(ν)(τ =0)
(ν)
(ν)
send µ∆xs →f = µfs →∆xs , xs to incident fi ∈ Fid
i
end for
for Each fi ∈ Find do
(ν)
compute residual ri = zi − hi (x(ν) )
(ν)
compute Jacobian elements Ci,∆xp ; ∆xp ∈ Xi
end for
procedure I NNER I TERATION LOOP τ = 1, 2, . . .
while stopping criterion for the inner loop is not met do
for Each fi ∈ Find do
(τ )
Compute µf →∆xs using (53)*
i
end for
for Each ∆xs ∈ X do
(τ )
Compute µ∆xs →f using (51)
i
end for
end while
end procedure
for Each ∆xs ∈ X do
(ν)
compute ∆x̂s using (56)
(ν)
(ν)
(ν+1)
update xs
= xs + ∆x̂s
end for
end while
end procedure
*Incoming messages are obtained in previous iteration τ − 1

A. Simulation Setup

convergence condition is required to hold in every outer
iteration ν.
Similarly, the convergence of the GN-BP with randomized
damping in every outer iteration loop ν is governed by the
spectral radius of the matrix:
Ω̄(x(ν) ) = QΩ(x(ν) ) + α2 RΩ(x(ν) ) − α1 R.

In the simulated model, we start with a given IEEE bus
test case and apply the power flow analysis to generate the
exact solution. In case of the DC model, we apply the DC
power flow analysis to calculate voltage angles and active
powers, while for the AC model, we apply the AC power
flow analysis to generate currents, voltages and powers
across the network. Further, we corrupt the exact solution by
the additive white Gaussian noise of variance v and we
observe the set of measurements.
In the DC model, the measurements contain active power
flows and power injections, and bus voltage angles, while in
the AC model, they contain active and reactive power flows,
active and reactive power injections, bus voltage magnitudes
and bus voltage angles. The set of measurements is selected
in such a way that the system is observable. More precisely,
for each scenario, we generate 1000 random measurement
configurations with the number of measurements equal either
to double or triple the size of the number of state variables
(i.e., we consider the redundancy to be equal 2 or 3).
To evaluate the performance of the BP-based SE algorithms,
we convert each of the above randomly generated IEEE bus
test cases with a given measurement configuration into the
corresponding factor graph and we run the DC-BP, AC-BP or
GN-BP algorithms over the factor graph. Note that, in order
to initialize the AC-BP, the GN-BP, and the Gauss-Newton
method, we use the ”flat start” assumption2 (Vi = 1, θi = 0,
i = 1, . . . , N ).
B. Convergence Analysis of DC-BP Algorithm

(58)

According to Theorem 3, in every outer iteration ν, the
algorithm will converge to a unique fixed point if
ρ(Ω̄(x(ν) )) < 1, ν = {0, 1, 2, . . . }, and the resulting fixed
point is equal to the fixed point obtained by the GN-BP with
synchronous scheduling.
VI. N UMERICAL R ESULTS
In this section, using numerical simulations, we analyze the
convergence and the accuracy of the BP-based SE methods.
For convergence evaluation, we focus on the DC-BP
algorithm whose convergence in analyzed in Sections IV-D
and IV-E. We verify these results via numerical simulations.

As detailed in Sec. IV-D and IV-E, the DC-BP with
synchronous scheduling with or without randomized
damping will converge if ρ(Ω) < 1 and ρ(Ω̄) < 1,
respectively. This condition is verified in our simulations,
thus we present the convergence performance by comparing
spectral radii of matrices Ω and Ω̄.
Figure 6 shows spectral radii for different redundancies for
IEEE 14 and IEEE 118 bus test case. For each scenario, the
randomized damping case behaves superior in terms of the
2 Usually, IEEE bus test cases include in-phase transformers with zero series
resistances. Consequently, the ”flat start” will cause undefined messages due to
the corresponding Jacobian elements being equal to zero. This problem can be
alleviated by initializing bus voltage angles with a small random perturbation
at transformers’ buses.

11

1.0
Synchronous
Scheduling
0.9
Randomized
0.8
Damping
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4
Spectral Radius ρ

Empirical CDF

Empirical CDF

spectral radius. As an interesting and somewhat extreme
case, for the IEEE 118 bus test case, the DC-BP algorithm
with synchronous scheduling could not converge at all, while
with randomized damping3 , we recorded convergence with
probability above 0.9. As expected, the algorithm with
randomized damping performs better for larger redundancy.

1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0.8

Synchronous
Scheduling
Randomized
Damping

0.9
1.0
1.1
Spectral Radius ρ

1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

(b)

0.8950

Empirical CDF

Empirical CDF

(a)

0.8945
0.8940
0.8935

1.0

Synchronous
Scheduling
Randomized
Damping

1.0

1.2

1.1 1.2 1.3 1.4
Spectral Radius ρ

1.5

1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

0.9935

0.9930

0.9925

0.9920

1.0

Note that WRSS is the value of the objective function of the
optimization problem (10) we are solving, thus it is suitable
metric for the SE accuracy. Finally, we normalize the obtained
WRSS by WRSSWLS of the centralized SE obtained using the
Gauss-Newton method after 12 iterations (which we adopt as
a normalization constant). This way, we compare the accuracy
of BP-based algorithms to the one of the centralized SE.
Besides WRSS, we use the mean absolute difference (MAD)
between the state variables in two consecutive iterations:
n
1X
|∆xi |.
MAD =
(60)
n i=1
The MAD value represents average component-wise shift of
the state estimate over the iterations, thus it may be used to
quantify the rate of convergence.
The IEEE 14 bus test case with fixed measurement
configuration containing 61 measurement devices, as shown
in Fig. 7, is used to compare the performance of the
BP-based SE algorithms. For each value of noise variance
v = {v1 , v2 , v3 , v4 } = {10−10 , 10−8 , 10−6 , 10−4 } p.u.,
using Monte Carlo approach, we generate 1000 random sets
of measurement values and feed them to the proposed
BP-based SE algorithms.
13

Synchronous
Scheduling
Randomized
Damping

1.0

1.1 1.2 1.3 1.4
Spectral Radius ρ

(c)

14

1.5

12
11

(d)

10

Fig. 6. The spectral radius of matrices Ω for synchronous scheduling and
Ω̄ for randomized damping for redundancy equal 2 for IEEE 14 (subfigure
a) and IEEE 118 (subfigure c) bus test case and for redundancy equal 3 for
IEEE 14 (subfigure b) and IEEE 118 (subfigure d) bus test case.

In scenarios when the standard DC-BP algorithm did not
converge (i.e., ρ(Ω) > 1), the algorithm with randomized
damping had smaller spectral radius, i.e., ρ(Ω̄) < ρ(Ω), and
moreover, the spectral radius was usually smaller than one
(ρ(Ω̄) < 1). On the other hand, when DC-BP converged
(ρ(Ω) < 1), then in most cases the randomized damping
method also converged, albeit with slightly higher spectral
radius (ρ(Ω) < ρ(Ω̄) < 1). These results demonstrate the
benefit of randomized damping for improved DC-BP
stability and convergence4 .
C. Performance Evaluation of AC-BP and GN-BP
In the following, compare the performance of the AC-BP
and GN-BP algorithms via extensive numerical study (the ACBP is presented in detail in Appendix B and C). We use the
weighted residual sum of squares (WRSS) as a metric:
WRSS =

k
X
[zi − hi (x)]2
i=1

σi2

.

(59)

3 Note that randomized damping parameters, which are optimized via
numerical experiments in Sec. VI-D, are set to p = 0.6 and α1 = 0.5.
4 Another trick for improving the convergence behavior is to connect virtual
factor nodes to all variable nodes that do not already have direct factor nodes.

6
9
1
7

5
4

8

2
3
Active and Reactive Power Flow
Injection Active and Reactive Power
Voltage Magnitude and Angle

Fig. 7. The IEEE 14 bus test case with given measurement configuration.

For both algorithms, we apply randomized damping due to
improved convergence, and initially, we set the randomized
damping parameters to p = 0.6 and α1 = 0.5. The GN-BP is
configured so that, for every outer iteration ν, the number of
inner iterations grows exponentially as η(ν) = ν e , where e is
the inner iteration exponent which we set to e = 4 [40]5 .
In order to present average performance results, we use box
plots in which we: i) set the minimum and the maximum value
to 0% and 90% of the obtained simulation results, ii) the box
is defined by the first and the third quartile, iii) the middle line
is the median value, and iv) the results above the maximum
value are represented as outliers.
Fig. 8 shows the weighted residual sum of squares of the
AC-BP WRSSτBP over the iterations τ , normalized by
5 Note that we perform detailed performance analysis as a function of
damping parameters in [40], Appendix E, and as a function of the number of
inner iterations in Section VI-D.

WRSSBP/WRSSWLS

WRSSτBP/WRSSWLS

1.25
1.20
1.15
1.10
1.05

WRSSBP/WRSSWLS

12

1.00003
1.00002
1.00001
1.00000

100

300

500

700

1.0003
1.0002
1.0001
1.0000

AC - BP

1.00

1.0005
1.0004

GN - BP

AC - BP

(a)

900 1100 1300 1500 1700 1900 2100 2275
Number of iterations τ

GN - BP

(b)

1.4
1.3

1.002

1.3
1.2
1.1
1.0
AC - BP

GN - BP

500

700

900 1100 1300 1500 1700 1900 2100 2275
Number of iterations τ

Fig. 8. The AC-BP normalized WRSS (i.e., WRSSτBP /WRSSWLS ) for the
low noise level v1 (subfigure a) and the high noise level v4 (subfigure b).

WRSSWLS (i.e., WRSSτBP /WRSSWLS ). We observe that the
AC-BP converges for both the low and the high noise level,
however, for the high noise level, the solution of the AC-BP
algorithm does not correspond to the solution of the
centralized SE. This is expected, since as the noise variance
increases, the accuracy of Gaussian approximation of the BP
messages is decreasing, which affects the accuracy of the
AC-BP solution.
Fig. 9 shows WRSSνBP /WRSSWLS of the GN-BP over the
outer iterations ν. As shown, the GN-BP algorithm converges
to the solution of the centralized SE for both the low and the
high noise level.

Number of inner iterations τ
979 2275 4676 8722 15333 25333 39974 60710

WRSSνBP/WRSSWLS

Number of inner iterations τ
979 2275 4676 8722 15333 25333 39974 60710

2.0
1.5

6
7
8
9 10 11 12
Number of outer iterations ν

Fig. 10. The normalized WRSS for the AC-BP and GN-BP after 2275
iteration for variances v1 (subfigure a), v2 (subfigure b), v3 (subfigure c)
and v4 (subfigure d).

addition, in terms of implementation, the AC-BP suffers two
drawbacks as compared to the GN-BP: i) the AC-BP
messages have considerably more complex form, and ii) the
AC-BP requires prior knowledge (e.g., historical data).
Consequently, in the remaining part of this section, we focus
on the performance of the GN-BP.
D. Design of Inner Iteration Schemes
In order to reduce the number of non-converging
simulations due to an insufficient number of inner iterations,
in this section we test the GN-BP convergence performance
with different inner iteration schemes. The damping
parameters, investigated in detail in Appendix E, are set to
p = 0.4 and α1 = 0.3. As a reference point, the MAD values
of the centralized SE for variances v1 and v4 are shown in
Fig. 11.

1.004
1.003
1.002
1.001
1.000

1.0

GN - BP

(d)

5

6
7
8
9 10 11 12
Number of outer iterations ν

(b)

(a)
WRSSν

Fig. 9. The GN-BP normalized WRSS (i.e.,
BP /WRSSWLS ) for the
low noise level v1 (subfigure a) and the high noise level v4 (subfigure b).

To compare BP-based algorithms, we consider WRSSBP
of the BP-based algorithms after the same total number of
iterations. More precisely, the AC-BP algorithm is terminated
after τ = 2275 iterations, while the GN-BP is terminated
after ν = 5 outer iterations, which corresponds exactly to the
total of 2275 inner iterations. Fig. 10 demonstrates that, for
the low noise level, the AC-BP converges faster then the
GN-BP. However, as the noise level increases, the GN-BP
converges faster than the AC-BP.
To summarize, the GN-BP is able to reach exactly the
same solution as the centralized SE, while the AC-BP
generally fails to achieve it for higher noise levels. In

10−2

10−2

10−4

10−4

10−6

10−6
MAD

300

(b)

WRSSνBP/WRSSWLS

1.004

1.4

(c)

1.1

100

5

1.006

AC - BP

1.2

1.0

2.5

1.008

1.000

MAD

WRSSτBP/WRSSWLS

1.5

WRSSBP/WRSSWLS

WRSSBP/WRSSWLS

(a)

10−8
10−10

10−8
10−10

10−12

10−12

10−14

10−14

10−16

10−16
1 2 3 4 5 6 7 8 9 10 11 12
Number of iterations ν

1 2 3 4 5 6 7 8 9 10 11 12
Number of iterations ν

(a)

(b)

Fig. 11. The MAD values of the centralized SE for the low noise level v1
(subfigure a) and high noise level v4 (subfigure b).

For comparison, the MAD values of the GN-BP with the
exponential inner iteration scheme (η(ν) = ν e , e = 4) is
shown in Fig. 12. Note that the GN-BP exhibits comparable
convergence performance to the centralized SE algorithm.
Note also that it is difficult to directly compare the two, due
to a large difference in computational loads of a single
(outer) iteration. For example, the complexity of a single
iteration remains constant but significant (due to matrix

13

3 3 4 0
4 9 75 76 22 33 33 97 71
17 98 35 97 22 46 87 15 25 39 60

1

10−2

10−2

10−4

10−4

10−6

10−6

10−10

10−4

10−6

10−6

10−8
10−10
10−12

10−12
10−14

10−16

10−16
1 2 3 4 5 6 7 8 9 10 11 12
Number of outer iterations ν

1 2 3 4 5 6 7 8 9 10 11 12
Number of outer iterations ν

(a)

(b)

10−10

10−12

10−12
10−14

10−16

10−16
1 2 3 4 5 6 7 8 9 10 11 12
Number of outer iterations ν

1 2 3 4 5 6 7 8 9 10 11 12
Number of outer iterations ν

(a)

(b)

Fig. 12. The MAD values of the GN-BP with the exponential inner iteration
scheme for the low noise level v1 (subfigure a) and high noise level v4
(subfigure b).

In the following, we define an alternative accuracy-based
inner iteration scheme. Namely, instead of prescribing the
number of inner iterations in advance, we let the number of
inner iterations evolve until a certain accuracy-based
criterion is met. More precisely, the algorithm in the inner
iteration loop is running until the following convergence
criterion is reached:
(ν,τ −1

|rf →∆x − rf →∆x )| < (ν),

(61)

where rf →∆x represents the vector of mean-value messages
from factor nodes to variable nodes, and (ν) is the
threshold at iteration ν. Note that, due the fact that the BP is
a distributed algorithm, imposing this condition in real
systems would require additional communication overhead.
Fig. 13 shows MAD values and the total number of inner
iterations for the threshold (ν) = [10−3 , 10−6 , . . . , 10−15 ,
10−16 , . . . , 10−16 ], where the number of outer iterations is
ν = 1, . . . , 12. Comparing the results with the centralized SE
(Fig. 11), we note that the GN-BP has almost identical
convergence rate. We note that the proposed accuracy-based
inner iteration scheme dramatically reduces the number of
non-converging simulations.
To summarize, we have demonstrated that by using
different inner iteration schemes, one can influence and
control the trade-off between the GN-BP convergence rate
and the total number of inner iterations. For additional
results on scalability and complexity of GN-BP, we refer the
reader to Appendix F.
VII. C ONCLUSIONS
In this paper, we presented an in-depth study of the
application of the BP algorithm to the SE problem in power
systems. We provided detailed derivation, convergence and
performance analysis of BP-based SE algorithms for both
DC and AC model. The main contribution of our study is
the GN-BP algorithm, which is shown to represent a
BP-based implementation of the iterative Gauss-Newton
method. In a thorough numerical study, we have shown how

10−8
10−10

10−14

10−8

10−14

(ν,τ )

10−2

10−4

200000

200000

100000
60710
39974
25333
15333
8722
4676

100000
60710
39974
25333
15333
8722
4676

Number of inner iterations τ

10−8

3 3 4 0
4 9 75 76 22 33 33 97 71
17 98 35 97 22 46 87 15 25 39 60

Number of inner iterations τ

MAD

MAD

1

Number of inner iterations τ

10−2

MAD

Number of inner iterations τ

MAD

inversion) over iterations for the centralized SE algorithm,
while it gradually increases for the GN-BP starting from an
extremely low complexity at initial outer iterations.

2275
979
354
98

2275
979
354
98

1 2 3 4 5 6 7 8 9 10 11 12
Number of outer iterations ν

1 2 3 4 5 6 7 8 9 10 11 12
Number of outer iterations ν

(c)

(d)

Fig. 13. The MAD values and the number of inner iterations for the low noise
level v1 (subfigure a and subfigure c) and the high noise level v4 (subfigure
b and subfigure d) of the GN-BP with threshold (ν) = [10−3 , 10−6 , . . . ,
10−15 , 10−16 , . . . , 10−16 ].

to fine-tune the GN-BP implementation and balance between
the accuracy and the convergence rate. In our future work,
we will investigate extensions to asynchronous BP in
dynamic and real-time state estimation, including suitable
models for communication constraints.

APPENDIX A
M EASUREMENT F UNCTIONS AND JACOBIAN E LEMENTS
In the following, we present measurement functions and
corresponding Jacobian elements for both, the AC and DC SE
model. Note that, for simplicity, we assume that the model
does not contain phase-shifting transformers.

AC SE Model
The active and reactive power flow at the branch (i, j) ∈ E
that connects buses i and j:
hPij (·) = Vi2 (gij + gsi ) − Vi Vj (gij cos θij + bij sin θij )

(62a)

hQij (·) = −Vi2 (bij + bsi ) − Vi Vj (gij sin θij − bij cos θij ),(62b)

where Vi and Vj are bus voltage magnitudes, while θij =
θi −θj is the bus voltage angle difference between bus voltage
angles at buses i and j. The parameters in above equations
include the conductance gij and susceptance bij of the branch,
as well as the conductance gsi and susceptance bsi of the
branch shunt element connected at the bus i.

14

The Jacobian expressions corresponding to hPij (·) are:
∂hPij (·)
= Vi Vj (gij sinθij − bij cosθij )
∂θi
∂hPij (·)
= −Vi Vj (gij sinθij − bij cosθij )
∂θj
∂hPij (·)
= −Vj (gij cosθij + bij sinθij ) + 2Vi (gij + gsi )
∂Vi
∂hPij (·)
= −Vi (gij cosθij + bij sinθij ).
∂Vj
The Jacobian expressions corresponding to hQij (·) are as
follows:
∂hQij (·)
= −Vi Vj (gij cosθij + bij sinθij )
∂θi
∂hQij (·)
= Vi Vj (gij cosθij + bij sinθij )
∂θj
∂hQij (·)
= −Vj (gij sinθij − bij cosθij ) − 2Vi (bij + bsi )
∂Vi
∂hQij (·)
= −Vi (gij sinθij − bij cosθij ).
∂Vj

Jacobian expressions corresponding to hPi (·) are:
X
∂hPi (·)
= Vi
Vj (−Gij sinθij + Bij cosθij )
∂θi
j∈Hi \i

∂hPi (·)
= Vi Vj (Gij sinθij − Bij cosθij )
∂θj
X
∂hPi (·)
=
Vj (Gij cosθij + Bij sinθij ) + 2Vi Gii
∂Vi
j∈Hi \i

∂hPi (·)
= Vi (Gij cosθij + Bij sinθij ),
∂Vj
where Hi \ i is the set of buses adjacent to the bus i. The
Jacobian expressions corresponding to hQi (·) are:
X
∂hQi (·)
= Vi
Vj (Gij cosθij + Bij sinθij )
∂θi
j∈Hi \i

∂hQi (·)
= Vi Vj (−Gij cosθij − Bij sinθij )
∂θj
X
∂hQi (·)
=
Vj (Gij sniθij − Bij cosθij ) − 2Vi Bii
∂Vi
j∈Hi \i

∂hQi (·)
= Vi (Gij sinθij − Bij cosθij ).
∂Vj
The current magnitude at the branch (i, j) ∈ E that
connects buses i and j:
hIij (·) = [Ac Vi2 + Bc Vj2
− 2Vi Vj (Cc cos θij − Dc sin θij )]1/2 . (63)
Ac = (gij + gsi )2 + (bij + bsi )2 ;

2
Bc = gij
+ b2ij

Cc = gij (gij + gsi ) + bij (bij + bsi ); Dc = gij bsi − bij gsi .
The Jacobian expressions corresponding to current magnitude
measurement function hIij (·) are:
∂hIij (·) Vi Vj (Dc cos θij + Cc sin θij )
=
∂θi
Iij
∂hIij (·)
Vi Vj (Dc cos θij + Cc sin θij )
=−
∂θj
Iij
∂hIij (·) Vj (Dc sin θij − Cc cos θij ) + Ac Vi
=
∂Vi
Iij
∂hIij (·) Vi (Dc sin θij − Cc cos θij ) + Bc Vj
=
.
∂Vj
Iij

The voltage magnitude and angle on the bus i ∈ V:
hVi (·) = Vi ; hθi (·) = θi .

(65)

The Jacobian expressions corresponding to measurements of
voltage magnitude and voltage angle are as follows:
∂hVi (·)
= 1;
∂Vi
∂hθi (·)
= 1;
∂θi

∂hVi (·)
=0
∂Vj
∂hθi (·)
= 0.
∂θj

DC SE Model
The active power flow at the branch (i, j) ∈ E that connects
buses i and j:
hPij (·) = −bij (θi − θj ).

(66)

The Jacobian expressions corresponding to hPij (·) are as
follows:
∂hPij (·)
∂hPij (·)
= −bij ;
= bij .
∂θi
∂θj

The active and reactive power injection into the bus i ∈ V:

The active power injection into the bus i ∈ V is:
X
hPi (·) = −
bij (θi − θj ),

(67)

j∈Hi \i

hPi (·) = Vi

X

Vj (Gij cos θij + Bij sin θij )

(64a)

Vj (Gij sin θij − Bij cos θij ),

(64b)

j∈Hi

hQi (·) = Vi

X
j∈Hi

where Hi is the set of buses adjacent to the bus i, including
the bus i. The parameters Gij and Bij are conductance and
susceptance of the complex bus matrix [1, Sec. 2.3]. The

where Hi \ i is the set of buses adjacent to the bus i. The
Jacobian expressions corresponding to hPi (·) are as follows:
X
∂hPi (·)
=−
bij ;
∂θi
j∈Hi \i

X
∂hPi (·)
=
bij .
∂θj
j∈Hi \i

The voltage angle on the bus i ∈ V:
hθi (·) = θi .

(68)

15

The Jacobian expressions corresponding to hθi (·) are as
follows:
∂hθi (·)
= 1;
∂θi

∂hθi (·)
= 0.
∂θj

APPENDIX B
BP- BASED A LGORITHM FOR THE AC SE
The AC model is characterized by the set of state
variables x ≡ (θ, V), while measurement functions are
defined in Appendix A, equations (62)-(65).
The factor graph construction: According to (9), in the
AC scenario, the set of state variables θ and V determines the
set of variable nodes X = {(θ1 , V1 ), . . . , (θN , VN )}, while
the set of factor nodes F = {f1 , . . . , fk } is defined by the set
of measurements M. A factor node fi connects to a variable
node xs ∈ X if and only if the state variable xs is an argument
of the corresponding measurement function hi (x).
Example 3 (Constructing factor graph). In this example,
using a simple 3-bus model presented in Fig. 14a, we
demonstrate the conversion from a bus/branch model with a
given measurement configuration into the corresponding
factor graph for the AC model. The variable nodes represent

bus 1
θ1
2

1 MP12

bus 2
fP12

V1

θ2
V2

fP3
3

MP 3

(a)

θ3

V3

bus 3

(b)
Fig. 14.
Transformation of the bus/branch model and measurement
configuration (subfigure a) into the corresponding factor graph for the DC
model (subfigure b) and the AC model (subfigure c)

state variables X = {θ1 , V1 , θ2 , V2 , θ3 , V3 }. Factor nodes
are defined by corresponding measurements, where in our
example, measurements MP12 and MP3 are mapped into
factor nodes F = {fP12 , fP3 }.
4
Message from a variable node to a factor node: Due to
the fact that variable node output messages do not depend on
measurement functions (see equation (2)), relations (13) and
(14) hold for the AC-BP.
Message from a factor node to a variable node: Due to
non-linear measurement functions hi (·), the integral in (3)
for the AC-BP cannot be evaluated in closed form.
Consequently, the message from a factor node to a variable
node will not be Gaussian. In the following, as an
approximation, we assume that for the AC-BP, the message
µfi →xs (xs ) also has the Gaussian form (15). According to
DC-BP we provide arguments that lead us to approximations
used to derive messages for the AC-BP.

Mean value evaluation: The expression for the mean of
the DC-BP zfi →xs is exact and equals (19a). Although the
expression (19a) is obtained by directly evaluating (3) for the
linear DC model, we note that it has a useful interpretation
via conditional expectation. For that purpose, let us define a
vector xb = Xi \ xs , and let zxb →fi denote a vector of mean
values of messages from variable nodes Xi \ xs to the factor
node
fi .
Then,
the
conditional
expectation
E[hi (xs , xb )|xb = zxb →fi ] can be calculated as:
E[hi (xs , xb )|xb = zxb →fi ]
= Cxs E[xs |xb = zxb →fi ] +

X

Cxb zxb →fi = zi . (70)

xb ∈Xi \xs

From the BP perspective, the conditional expected value
E[xs |xb = zxb →fi ] represents the mean zfi →xs . Hence, it is
possible to define the conditional expectation of non-linear
measurement function hi (·):
E[hi (xs , xb )|xb = zxb →fi ] = zi .

(71)

Due different forms of non-linear measurement functions
hi (·), see equations (62)-(64), the equation (71) will produce
different forms of conditional expectation E[xs |zxb →fi ] ≡
zfi →xs :
aE[x2s |xb

aE[xs |zxb →fi ] + b = 0

(72a)

= zxb →fi ] + bE[xs |xb = zxb →fi ] + c = 0

(72b)

aE[sin2 xs |xb = zxb →fi ] + bE[sin xs |xb = zxb →fi ] + c = 0,

(72c)

where a, b and c are coefficients derived from non-linear
measurement functions (see Appendix C for details).
Due to quadratic form of (72b) and (72c), we may obtain
two possible values for the mean value zfi →xs . Thus in order
to unambiguously define zfi →xs , we assume that certain a
e V),
e
e ≡ (θ,
priori knowledge of state variables, denoted as x
is available (e.g., historical data). Given the prior data, we
evaluate the mean value as:
 (1)

, if ∆ > 0 and d1 < d2
z

 fi →xs
(2)
(73)
zfi →xs = zfi →xs , if ∆ > 0 and d1 > d2



x
es ,
if ∆ < 0
where ∆ is the discriminant of the quadratic polynomial, and
(1)
(2)
d1 = |zfi →xs − x
es |, d2 = |zfi →xs − x
es |, (see Appendix C for
details).
The variance evaluation: The expression for the variance
of the DC-BP vfi →xs is equal (19b). Let us provide another
interpretation of the variance vfi →xs . For this purpose, we
observe the factor graph presented in Fig. 15. Consider the
set of messages µxb →fi = N (zxb →fi |xb , vxb →fi ) arriving to
the factor node fi from any variable node neighbour
xb ∈ Xi . Informally, we note that this message carries a
“belief” about itself that the variable node xb sends to the
factor node fi , representing collective evidence the rest of
the factor graph provides about the variable node xb . Let us
represent this belief by an equivalent factor node attached to
each variable node. Thus for a set of variable nodes Xi , we
introduce a set of factor nodes Feq = {fs , fl , . . . , fL },

16

It can be shown, using (78), that the variance vfi →xs is
governed by (19b) where the coefficients Cxp , xp ∈ Xi are
defined by Jacobian elements (see Appendix A and C for
details):
∂hi (·)
.
C xp =
(79)
∂xp xs =zfi →xs

xl

fl

µ xl → f i ( x l )
x
fi µxs →fi (xs ) s

.
.
.

.
.
.
fL

fs

µ xL → f i ( x L )

xb =zxb →fi

xL

Fig. 15. Factor graph which provides interpretation of the variance σf2

i →xs

.

where for each xb ∈ Xi , the corresponding factor node
fb
∈
Feq is singly-connected to xb and by
N (zxb →fi |xb , vxb →fi ). Note that, from the perspective of
SE, this factor node can observed as a measurement defined
by the value zxb →fi , variance vxb →fi , and measurement
function hb (xb ) = xb .
Let us now solve the system illustrated in Fig. 15 using
the WLS method. It is easy to show that the corresponding
Jacobian matrix6 H and weighted matrix W have the
following form:


Cxs Cxl . . . CxL
 1
0
...
0 


 0
1
.
.
.
0 
H=
(74)

 ..
..
.. 
 .

.
.
0
0
...
1
W = diag(1/vi , 1/vxs →fi , 1/vxl →fi , . . . , 1/vxL →fi ). (75)
A variance-covariance matrix of WLS method is defined as:
V(xi ) = (HT WH)−1




var(xs ) cov(xs , xl ) . . . cov(xs , xL )
 cov(xl , xs ) var(xl ) . . . cov(xl , xL ) 



=


..
.

..
.

.


..
.

cov(xL , xs ) cov(xL , xl ) . . .

(76)

var(xL )

According to (76), and using (74) and (75), the variance
var(xs ) is:
1
var(xs )

=

"

1
vxs →fi

+

1 
Cx2s

vi +

X

Cx2b vxb →fi



#−1
.

(77)

xb ∈Xi \xs

Consider the second term on the right-hand side of (77). Recall
that it represents the inverse of the variance vfi →xs of the
message from the factor node fi to the variable node xs , as
defined by (19b). Therefore, we have demonstrated that by
applying WLS on the factor graph in Fig. 15, one can obtain
the expression for the variance of the message from the factor
node fi to the variable node xs .
For the AC SE that deals with non-linear measurement
functions, it is possible to define a linear approximation of
the variance-covariance matrix at a given point xi using the
Gauss-Newton method (11a):
V(xi ) = [J(xi )T WJ(xi )]−1 .

(78)

6 Note that the measurement function of the factor node f is given by (18),
i
while for all other factor nodes fb ∈ Feq , it is equal to hb (xb ) = xb .

Note that the coefficients above are evaluated at the point xi =
(xs , xb ), where the values in xi represent the mean-values of
the corresponding messages.
To summarize, the message evaluation for the AC-BP is
governed by (72) and (19b), where coefficients are obtained
using (79).
Marginal inference: The marginal of the state variable xs
is governed by (21).
Iterative AC-BP algorithms Here, the indirect factor
nodes Find ⊂ F include measurements of power flows,
power injections and currents. The direct factor nodes
Fdir ⊂ F include measurements of voltage magnitudes and
angles.
The AC-BP algorithms are presented in Algorithm 3. Note
that, the initialization step for the DC-BP and AC-BP is
different. This is due to the fact that the variance of the
message from a factor node to a variable node for the
AC-BP depends not only on the mean values of incoming
messages, but also on the mean value of the message whose
variance is being calculated.
Algorithm 3 The AC-BP
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:

procedure I NITIALIZATION τ = 0
for Each fs ∈ Floc do
(0)
send µfs →xs to incidence xs ∈ X
end for
for Each xs ∈ X do
(0)
(0)
send µxs →f = µfs →xs , to incidence fi ∈ Fid
i
end for
for Each fi ∈ Find do
(0)
(0)
send µf →xs = µxs →f to incidence xs ∈ X
i
i
end for
end procedure
procedure I TERATION LOOP τ = 1, 2, . . .
while stopping criterion is not met do
for Each fi ∈ Find do
(τ )
Compute µf →xs using (72)*, (19b)*
i
end for
for Each xs ∈ X do
(τ )
Compute µxs →f using (14)
i
end for
end while
end procedure
procedure O UTPUT
for Each xs ∈ X do
Compute x̂s , vxs using (21)
end for
end procedure
*Incomming messages are obtained in previous iteration τ − 1

17

The mean zfi →Vi follows from the quadratic equation, where
we selected a solution using (73).

APPENDIX C
T HE AC-BP A LGORITHM : D ETAILED M ESSAGE
D ERIVATION
Here we present an example of evaluation of the message
from a factor node to a variable node for the AC-BP algorithm.
We consider a simple model containing buses i and j, with
the active power flow measurement Mi ≡ MPij at the branch
(i, j). The mean zi , variance vi and the measurement function
hi (θi , Vi , θj , Vj ) defined as (62a) is associated with the active
power flow measurement Mi . The corresponding factor graph
is shown in Fig. 16.
Further, all incoming messages from variable nodes to the
factor node fi have Gaussian form. Therefore, these
messages, denoted as µθi →fi (θi ), µVi →fi (Vi ), µθj →fi (θj )
and µVj →fi (Vj ), are represented by their mean-variance pair
(zθi →fi , vθi →fi ), (zVi →fi , vVi →fi ), (zθj →fi , vθj →fi ) and
(zVj →fi , vVj →fi ), respectively (Fig. 16a - Fig. 16d).

θi

µ θ i →fi ( θ i )

µ θ j →fi ( θ j )
fi

Vi
µ fi → V i ( Vi )

θj

θi

Vj

Vi

µ fi →θ i ( θ i )

θj

µ θ j →fi ( θ j )

Vi
µ V j → fi ( V j )

(c)

θj

θi

Vj

Vi

σf2i →Vi =

Cθi =

CVi =

∂hi (Vi , xb )
∂θi
∂hi (Vi , xb )
∂Vi

Vi =zfi →Vi
xb =zxb →fi

CVj =
Vi =zfi →Vi
xb =zxb →fi

∂hi (Vi , xb )
∂θj
∂hi (Vi , xb )
∂Vj

Vi =zfi →Vi
xb =zxb →fi

Vi =zfi →Vi
xb =zxb →fi

aE[Vj |xb = zxb →fi ] + b = 0,
where: xb = (θi , Vi , θj ) and zxb →fi = (zθi →fi , zVi →fi ,
zθj →fi ), with coefficients:

µ fi → θ j ( θ j )
fi

a = zi − zV2i →fi (gij + gsi )

θj

b = zVi →fi (gij cos zθij →fi + bij sin zθij →fi ).

Vj
µ V i → fi ( Vi )

Cθj =

• The message µfi →Vj (Fig. 16b): The mean zfi →Vj is defined
according to (72a) as:

Vj

µ θ i →fi ( θ i )

1
(vi + Cθ2i vθi →fi + Cθ2j vθj →fi + CV2j vVj →fi ),
CV2i

where coefficients are defined according to Jacobian elements
of the measurement function hi (·):

µ fi → V j ( V j )

(b)

fi

µ V i → fi ( Vi )

fi

µ V i → fi ( Vi )

µ V j → fi ( V j )

(a)

θi

µ θ j →fi ( θ j )

µ θ i →fi ( θ i )

The variance vfi →Vi is determined using (19b) as:

µ V j → fi ( V j )

(d)

Due the fact that the conditional expected value E[Vj |xb =
zxb →fi ] represents the mean zfi →Vj , we obtain:

Fig. 16. Messages from from factor node fi to variable nodes: Vi (subfigure
a), Vj (subfigure b), θi (subfigure c) and θj (subfigure d).

According to assumption (see Section Appendix B), the
messages from the factor node fi to variable nodes have
Gaussian form: µfi →Vi (Vi ), µfi →Vj (Vj ), µfi →θi (θi ) and
µfi →θj (θj ) with their mean-variance pair (zfi →Vi , vfi →Vi ),
(zfi →Vj , vfi →Vj ), (zfi →θi , vfi →θi ) and (zfi →θj , vfi →θj ). In
the following, we consider calculation of each of these
messages.
• The message µfi →Vi (Fig. 16a): Let us first consider the
mean zfi →Vi . The equation (19a) for the active power flow
measurement boils down to (72b):

azfi →Vj + b = 0.
The variance vfi →Vj is determined using (19b) as:
vfi →Vj =

1
(vi + Cθ2i vθi →fi + CV2i vVi →fi + Cθ2j vθj →fi ),
CV2j

where coefficient are defined according to Jacobian elements
of the measurement function hi (·).

The messages µfi →θi and µfi →θj (Fig. 16c and Fig. 16d):
Means zfi →θi and zfi →θj are defined according to (72c):
•

aE[Vi2 |xb = zxb →fi ] + bE[Vi |xb = zxb →fi ] + c = 0,

aE[sin2 xs |xb = zxb →fi ] + bE[sin xs |xb = zxb →fi ] + c = 0,

where: xb = (θi , θj , Vj ) and zxb →fi = (zθi →fi , zθj →fi ,
zVj →fi ), with coefficients:

where: xb = (Vi , θj , Vj ) and zxb →fi = (zVi →fi , zθj →fi ,
zVj →fi ) for the message µfi →θi , xb = (θi , Vi , Vj ) and
zxb →fi = (zθi →fi , zVi →fi , zVj →fi ) for the message µfi →θj ,
and xs ∈ {θi , θj }. Due the fact that the all variables and
messages preserve Gaussian distribution, the conditional
expectations of sine functions are equal to E[sin2 xs |xb =
zxb →fi ] = sin2 zfi →xs and E[sin xs |xb = zxb →fi ] =
sin zfi →xs , which allows us to compute the mean:

a = gij + gsi
b = −zVj →fi (gij cos zθij →fi + bij sin zθij →fi )
c = −zi ,
where zθij →fi is determined as zθi →fi −zθj →fi . Due the
fact that the conditional expected value E[Vi |xb = zxb →fi ]
represents the mean zfi →Vi , we can write:
a(zf2i →Vi + vfi →Vi ) + bzfi →Vi + c = 0.

a sin2 zfi →xs + b sin zfi →xs + c = 0.
To simplify expressions, we introduce coefficients a = A2 +

18

B 2 , b = −2BC and c = −A2 + C 2 :
A = gij cos zθj →fi − bij sin zθj →fi ,

xs ≡ θ i

A = gij cos zθi →fi + bij sin zθi →fi ,

xs ≡ θj

B = gij sin zθj →fi + bij cos zθj →fi ,

xs ≡ θ i

B = gij sin zθi →fi − bij cos zθi →fi ,

xs ≡ θi

C=

zV2i →fi (gij

+ gsi ) − zi
,
zVi →fi zVj →fi

xs ∈ {θi , θj }

The variance vfi →θi is determined using (19b) as:
vfi →θi =

1
(vi + CV2i vVi →fi + Cθ2j vθj →fi + CV2j vVj →fi ),
Cθ2i

where coefficients are defined, as above, by calculating
Jacobian elements of the measurement function hi (·).
The variance vfi →θj is determined according to (19b) as:
vfi →θj =

1
(vi + Cθ2i vθi →fi + CV2i vVi →fi + CV2j vVj →fi ),
Cθ2j

where coefficient follow Jacobian elements of the
measurement function hPi (·).
Using the same methodology, it is possible to define
corresponding equations for means and variances for every
type of measurement functions.
APPENDIX D
T HE GN-BP A LGORITHM : T OY E XAMPLE
An illustrative example presented in Fig. 5a will be used to
provide a step-by-step presentation of the proposed algorithm.
Input data for SE from measurement devices are
Gaussian-type functions represented by means and variances:
{zV1 , zθ2 , zθ3 , zP12 , zP3 } and {vV1 , vθ2 , vθ3 , vP12 , vP3 }.
The corresponding factor graph is given in Fig. 14, where
indirect factor nodes (red squares) are frP12 and frP3 , while
direct factor nodes (orange squares) are frV1 , frθ2 and frθ3 .
The slack factor node frθ1 (yellow square) corresponds to
the slack or reference bus. The initialization (green square)
and virtual (blue square) factor nodes are frV2 and frV3 ,
respectively.
Each variable node in the initialization step of the
algorithm has the corresponding local factor node. Local
factor nodes only send, but do not receive, the messages to
the incident variable nodes. Direct and virtual factor nodes
always repeat the same message to their corresponding
variable nodes through iterations.
f r θ1

f r θ2

Algorithm Initialization
1) The AC SE in electric power systems assumes ”flat start”
or a priori given values of state variables:
x(ν=0) = [θ1 θ2 θ3 V1 V2 V3 ](ν=0) .
2) The residual of the slack factor node is set to rθ1 = 0 with
variance vθ1 → 0.
3) The value of initialization factor nodes and virtual factor
nodes are set to rV2 → 0 and rV3 → 0, with variances
vV2 → ∞ and vV3 → ∞.
Iterate - Outer Loop: ν= 0, 1, 2, . . .; τ = 0
4) Each direct factor node computes residual, e.g.:
(ν)

(ν)

rV1 = zV1 − V1 .
5) Local factor nodes send messages represented by a triplet
(residual, variance, state variable), to incident variable
nodes, e.g.:
(ν)
(ν) 
µfr →∆θ1 := rθ1 , vθ1 , θ1
θ1
(ν)
(ν) 
(ν)
µfr →∆V1 := rV1 , vV1 , V1 .
V1

6) Variable nodes forward the incoming messages received
from local factor nodes along remaining edges, e.g.:
(ν)
(ν,τ )
(ν)
(ν) 
= µfr →∆θ1 := rθ1 , vθ1 , θ1
µ∆θ1 →fr
P12

θ1

(ν,τ )
r∆θ1 →fr
P

:=

P3

= µfr

θ1

→∆θ1

frP12

∆V1

∆V2

(ν,τ )

P3

(ν,τ ) 

P3

.

7) Indirect factor nodes compute residuals, e.g.:
(ν)

(ν)

(ν)

(ν)

(ν)

rP12 = zP12 − hP12 (θ1 , θ2 , V1 , V2 ).
8) Indirect factor nodes compute appropriate Jacobian
elements associated with state variables, e.g.:
(ν)

∂hP12 (·)
(ν)
(ν)
(ν) (ν)
= V1 V2 (g12 sinθ12 − b12 cosθ12 )
∂θ1
∂hP12 (·)
(ν)
(ν)
(ν)
=
= −V1 (g12 cosθ12 + b12 sinθ12 ).
∂V2

CP12 ,∆θ1 =
(ν)

CP12 ,∆V2

Iterate - Inner Loop: τ = 1, 2, . . . , η(ν)
9) Indirect factor nodes send messages as pairs along incident
edges according to (53), e.g.:

(τ )
(τ )
(τ )
µfP →∆θ2 := rfP →∆θ2 , vfP →∆θ2
(τ )

rfr

f rV2

P12

→∆θ2

12

=

1
(ν)
CP12 ,∆θ2

12

h
(ν)
(ν)
(ν,τ −1)
rP12 − CP12 ,∆θ1 · r∆θ1 →fr

(ν,τ −1)

(ν)

∆ θ 3 ∆ V3

(ν) 

:= rθ1 , vθ1 , θ1

(ν,τ )

−CP12 ,∆V1 · r∆V1 →fr
f r θ3

, θ1

∆θ2

f rP 3

f rV 1

P12

:= r∆θ1 →fr , v∆θ1 →fr , θ1

12

∆θ 1

, v∆θ1 →fr
(ν)

(ν)

(ν,τ )

µ∆θ1 →fr

(ν,τ ) 

(ν,τ )

12

f rV3

Fig. 17. Factor graph of the bus/branch model with given measurement
configuration shown in Fig. 5a.

(τ )

vfr

P12

→∆θ2
(ν)

=

1
(ν)
(CP12 ,∆θ2 )2

(ν)

P12

i

(ν,τ −1)

− CP12 ,∆V2 · r∆V2 →fr

P12

h
(ν)
(ν,τ −1)
vP12 + (CP12 ,∆θ1 )2 · v∆θ1 →fr

(ν,τ −1)

+(CP12 ,∆V1 )2 · v∆V1 →fr

P12

(ν)

P12

P12

i
.

(ν,τ −1)

+ (CP12 ,∆V2 )2 · v∆V2 →fr

P12

19

=

1
1
+ (τ )
vθ2 v
fr →∆θ2
P3

12

(ν)

(ν,τ )

r∆θ2 →fr

P12

P12

rθ2
+
vθ 2

=

!
(τ )
rfr →∆θ2
P3
(τ )
v∆θ2 →fr .
(τ )
P12
vfr →∆θ2
P
3

Iterate - Outer Loop: ν= 0, 1, 2, . . .; τ = η(ν)
11) Variable nodes compute marginals according to (21), e.g.:
(ν)

833
500
300
100
50
30
15
7

(ν)
∆θ̂2

P12

=

(ν)
rθ2

vθ 2

+

P12

(τ )

→∆θ2

(τ )

vfr

P12

→∆θ2

+

rfr

P3

→∆θ2

(τ )

vfr

P3

(ν+1)

(ν)

50
30

p = 0. 2

p = 0. 4

p = 0. 5
p = 0. 8

p = 0. 6
p = 1. 0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Weighting coefficient α1

(b)

!
(ν)

v̂∆θ2 .

→∆θ2

12) Variable nodes update the state variables, e.g.:
θ2

100

Fig. 18. The number of non-converging simulations (out of 1000 simulations
in total) for the low noise level v1 (subfigure a) and the high noise level v4
(subfigure b).

P3

(τ )

rfr

p = 0. 4
p = 0. 6
p = 1.0

300

(a)

1
1
1
=
+ (τ )
+ (τ )
vθ 2 v
vfr →∆θ2
fr
→∆θ2

(ν)

= θ2 + ∆θ̂2 .

13) Repeat steps 4-13 until convergence.
Note that, after step 8, initialization factor nodes are
removed from the factor graph. Also in each iteration, virtual
factor nodes repeat the same message as in the initial step
and messages from a virtual factor node to a variable node
should not be included in calculation of the marginals.
APPENDIX E
S ELECTION OF R ANDOMIZED DAMPING PARAMETERS
In this appendix, we analyze the GN-BP convergence
performance in order to select appropriate values for
randomized damping parameters p and α1 introduced in
equation (37). We generate 1000 random measurement
configurations with the redundancy equals to 3. The set of
measurement values are corrupted by additive noise of both
the low noise level v1 and the high noise level v4 . The
number of inner iterations is the same as in Section VI-B.
We record the number of non-converging simulations and
evaluate the MAD value, both of which will depend on the
convergence parameters p and α1 . To recall, the probability
p defines a fraction of a factor node to variable node
messages from the current iteration that are combined with
the corresponding messages from the previous iteration. The
weighting coefficient α1 defines the ratio that determines
how messages from the current and the previous iteration are
combined. For example, p = 0.2 specifies that 20% of the
messages from the current iteration will be combined with
their values in the previous iteration, while 80% of messages
are keeping the values calculated in the current iteration.
Furthermore, if α1 = 0.1, then for the 20% of messages, the
new value is obtained as a linear combination of the values
calculated in the current and the previous iteration with

Fig. 18 shows the number of non-converging simulations
for low noise level v1 and high noise level v4 as a function
of p and α1 . In general, for the selection of p and α1 for
which only a small fraction of messages are combined with
their values in the previous iteration, and that is the case for
p close to zero or α1 close to one, we observe a large
number of non-converging simulations. This clearly
demonstrates the necessity of using (37) to “slow down” the
BP progress, thus increasing the algorithm stability and
providing improved convergence.
For a fixed value of α1 , we investigate the convergence
rate of GN-BP for different values of p. We expect that, for
any selected α1 , the GN-BP will converge faster for smaller
values of p, as lower p leads to a reduced “slow down”
effect. For example, Fig. 19 shows median MAD values for
α1 = 0.5 over the range of p values (where the median is
calculated over convergent simulation outcomes only). Figure
clearly shows that lower p leads to faster convergence.
However, one needs to be careful with selection of p in
order to avoid the combinations of p and α1 that lead to
large number of non-converging outcomes.

Number of inner iterations τ
1

Number of inner iterations τ

3 3 4 0
4 9 75 76 22 33 33 97 71
17 98 35 97 22 46 87 15 25 39 60

10−1

1

10−5
10−7
10−9

3 3 4 0
4 9 75 76 22 33 33 97 71
17 98 35 97 22 46 87 15 25 39 60

10−1

p = 0.2
p = 0.5
p = 0.8
p = 1.0

10−3

MAD

(ν)

v̂∆θ2

p = 0. 2
p = 0. 5
p = 0. 8

500

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Weighting coefficient α1

(ν)

p(∆θ2 ) ∝ N (∆θ̂2 |∆θ2 , v̂θ2 )
1

841
Number of non - converging
simulations

1
(ν,τ )
v∆θ2 →fr
P

P12

p = 0.2
p = 0.5
p = 0.8
p = 1.0

10−3
10−5
MAD

P12

coefficients 0.1 and 0.9, respectively. In the following, we
numerically investigate the randomized damping parameter
pairs that lead to acceptable trade-off between the number of
non-converging simulations and the rate of convergence.

Number of non - converging
simulations

10) Variable nodes send messages as pairs along incident edges
to indirect factor nodes according to (14), e.g.:

(ν,τ )
(ν,τ )
(ν,τ )
µ∆θ2 →fr
:= r∆θ2 →fr , v∆θ2 →fr

10−7
10−9

10−11

10−11

10−13

10−13

10−15

10−15

1 2 3 4 5 6 7 8 9 10 11 12
Number of outer iterations ν

1 2 3 4 5 6 7 8 9 10 11 12
Number of outer iterations ν

(a)

(b)

Fig. 19. The median MAD values for α1 = 0.5 over the range of p values
for the low noise level v1 (subfigure a) and the high noise level v4 (subfigure
b).

20

Number of inner iterations τ
1

Number of inner iterations τ

3 3 4 0
4 9 75 76 22 33 33 97 71
17 98 35 97 22 46 87 15 25 39 60

1

10−2

10−2

10−4

10−4

10−6

10−6
MAD

MAD

To select convergence parameters, we observed the median
MAD values and the number of non-converging simulations.
Table I lists the p and α1 pairs that demonstrate good
performance of both metrics of interest. From the table, we
can select e.g., p = 0.4 and α1 = 0.3, to achieve good
balance between the number of non-converging simulations
and the convergence rate.

10−8
10−10

10−12

C ONVERGENCE PERFORMANCES

10−14

10−14

0.4
0.4
0.4
0.5
0.5
0.5
0.5
0.6
0.8
0.8
1.0
1.0

α1
0.1
0.2
0.3
0.2
0.3
0.4
0.5
0.5
0.5
0.6
0.6
0.7

Median MAD values per iteration (i.e. 7,8,9)
Variance v1
Variance v2
7
8
9
7
8
9
−8
−9
−11
−5
[10 ] [10 ] [10 ] [10 ] [10−6 ] [10−8 ]
5.547
5.218
4.522
5.929
5.351
4.928
4.264
4.762
5.850
4.887
5.590
4.593

5.123
4.206
3.271
5.340
4.543
3.754
2.835
4.034
5.179
3.823
5.231
3.638

9.942
9.427
6.697
13.782
10.061
7.316
5.309
7.354
14.215
8.595
13.882
7.097

5.402
4.848
4.663
5.634
5.160
4.826
3.971
4.528
5.503
4.644
5.456
4.500

4.776
3.851
3.134
5.177
4.847
3.605
2.887
3.179
5.140
3.948
4.963
3.613

10.411
8.775
6.593
14.042
9.514
8.060
5.096
6.902
14.454
8.800
13.086
7.084

Number of
non-converging
simulations
v1
v2
8
8
9
7
10
11
10
9
9
10
7
8

33
35
31
31
33
30
35
33
32
30
34
32

There are two main causes of non-converging behavior.
The main cause is that the number of inner iterations in each
outer iteration follows the exponential rule (τ = ν q ). In
general, the GN-BP performs very well with such an
exponential increase of inner iterations, however, it suffers
from insufficient number of inner iterations in the first few
outer iterations. Our experiments show that increasing the
number of inner iteration only in the first few outer iterations
can dramatically reduce the number of non-converging
simulations (see Section VI-D). The second cause is a
combination of “weakly” connected variable nodes and
measurements, especially if those nodes represent the buses
incident to transformers. In that case, to prevent the
non-converging behavior, it is sufficient to connect virtual
factor nodes to such variable nodes. Therefore, with some
adjustments, one can eliminate the non-converging cases and
make the GN-BP algorithm fully reliable.
APPENDIX F
S CALABILITY AND C OMPLEXITY OF GN-BP A LGORITHM
Scalability of the GN-BP algorithm: To demonstrate
scalability of the GN-BP algorithm, we provide MAD values
for the IEEE 30 and IEEE 118 bus test case. As before, we
generate 1000 random measurement configurations with the
redundancy equal to 3. The set of measurement values are
corrupted by additive noise of both the low noise level v1
and the high noise level v4 . The convergence parameters are
set to p = 0.4 and α1 = 0.3 and we use the exponential
inner iteration scheme. From Fig. 20 and Fig. 21, we note
that the GN-BP converges, however, the rate of convergence
slightly decreases as the size of the system increases. As
described in the previous subsection, one can use different
inner iteration schemes to fine-tune the rate of convergence.
Complexity of the GN-BP algorithm: Let us first
consider the complexity of the centralized SE based on the

10−16
1 2 3 4 5 6 7 8 9 10 11 12
Number of outer iterations ν

1 2 3 4 5 6 7 8 9 10 11 12
Number of outer iterations ν

(a)

(b)

Fig. 20. The MAD values of the GN-BP for the IEEE 30 bus test case with
the exponential inner iteration scheme (η(ν) = ν e ) for the low noise level
v1 (subfigure a) and high noise level v4 (subfigure b).

Number of inner iterations τ
1

Number of inner iterations τ

3 3 4 0
4 9 75 76 22 33 33 97 71
17 98 35 97 22 46 87 15 25 39 60

1

10−2

3 3 4 0
4 9 75 76 22 33 33 97 71
17 98 35 97 22 46 87 15 25 39 60

10−2

10−4

10−4

10−6

MAD

p

10−10

TABLE I

MAD

Convergence
parameters

10−8

10−12

10−16

3 3 4 0
4 9 75 76 22 33 33 97 71
17 98 35 97 22 46 87 15 25 39 60

10−8
10−10

10−6
10−8

10−12

10−10

1 2 3 4 5 6 7 8 9 10 11 12
Number of outer iterations ν

1 2 3 4 5 6 7 8 9 10 11 12
Number of outer iterations ν

(a)

(b)

Fig. 21. The MAD values of the GN-BP for the IEEE 118 bus test case with
the exponential inner iteration scheme (η(ν) = ν e ) for the low noise level
v1 (subfigure a) and high noise level v4 (subfigure b).

Gauss-Newton method. For n state variables and k
measurements, each iteration of the centralized SE involves a
constant number of matrix multiplications and a matrix
inversion whose complexity is of the order of O(kn2 ) and
O(n3 ), respectively. Due to the fact that k is proportional to
n for observable systems, and that the centralized SE takes
constant number of iterations to converge, the overall
complexity of the centralized SE scales as O(n3 ). However,
this can be reduced by employing matrix inversion
techniques that exploit the sparsity of involved matrices. The
resulting sparsity-aware centralized SE methods operate with
complexity that scales as O(n2 ) [41], [42].
The complexity of BP depends on the sparsity of the
underlying factor graph, as the computational effort per
iteration is proportional to the number of edges in the factor
graph. For each of the k measurements, the degree (the
number of incident edges) of the corresponding factor node
is limited by a (typically small) constant. Indeed, for any
type of measurements, the corresponding measurement
function depends only on a few state variables corresponding
to the buses in the local neighbourhood of the bus/branch
where the measurement is taken. Thus as n and k grow
large, the number of edges in the factor graph scales as

21

O(n), which corresponds to the computation complexity of
BP per iteration. The scaling of the number of BP iterations
ρ as n grows large is a challenging problem. Based on
discussion in [43] for full matrices, the number of iterations
is likely to scale with condition number of the underlying
matrix, which for well-conditioned matrices may scale as
low as O(1). However, we leave the more detailed analysis
on the scaling of the number of GN-BP iterations for our
future work.
To summarize, BP approach builds upon the factor graph
structure that directly exploits the underlying system sparsity,
thus achieving minimal complexity of O(n) per iteration,
while the scaling of the number of iterations needs further
study. In contrast to the optimized centralized methods
whose complexity scales as O(n2 ), the BP method can be
flexibly distributed by arbitrarily segmenting the underlying
factor graph into disjoint areas. In the extreme case of the
fully-distributed BP algorithm (that we focus on in this
paper), each factor graph node operates locally and
independently. Thus, the SE problem is distributed across
O(n) nodes, and if implemented to run in parallel, can be
O(n) times faster than the centralized solution. In addition,
for fully-distributed BP, none of the nodes need to store the
system-level matrices (whose storage-size typically scales as
O(n2 )), and storing only constant-size set of local
parameters is sufficient.
[2] F. F. Wu, K. Moslehi, and A. Bose, “Power system control centers: Past,
present, and future,” Proc. of the IEEE, vol. 93, pp. 1890–1908, Nov 2005.
[3] A. Monticelli, “Electric power system state estimation,” Proceedings of the
IEEE, vol. 88, no. 2, pp. 262–282, 2000.
[4] F. C. Schweppe and D. B. Rom, “Power system static-state estimation, part ii:
Approximate model,” IEEE Transactions on Power Apparatus and Systems,
vol. PAS-89, no. 1, pp. 125–130, Jan 1970.
[5] 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, 2012.
[6] G. N. Korres, “A distributed multiarea state estimation,” IEEE Transactions on
Power Systems, vol. 26, no. 1, pp. 73–84, Feb 2011.
[7] W. Jiang, V. Vittal, and G. T. Heydt, “Diakoptic state estimation using phasor
measurement units,” IEEE Transactions on Power Systems, vol. 23, no. 4, pp.
1580–1589, Nov 2008.
[8] L. Zhao and A. Abur, “Multi area state estimation using synchronized phasor
measurements,” IEEE Transactions on Power Systems, vol. 20, no. 2, 2005.
[9] A. Minot, Y. Lu, and N. Li, “A distributed gauss-newton method for power
system state estimation,” in 2016 IEEE Power and Energy Society General
Meeting (PESGM), July 2016, pp. 1–1.
[10] D. Marelli, B. Ninness, and M. Fu, “Distributed weighted least-squares
estimation for power networks,” IFAC-PapersOnLine, vol. 48, no. 28, 2015.
[11] X. Tai, Z. Lin, M. Fu, and Y. Sun, “A new distributed state estimation technique
for power networks,” in 2013 American Control Conference, June 2013.
[12] R. Ebrahimian and R. Baldick, “State estimation distributed processing [for
power systems],” IEEE Transactions on Power Systems, vol. 15, no. 4, pp.
1240–1246, 2000.
[13] A. J. Conejo, S. de la Torre, and M. Canas, “An optimization approach to
multiarea state estimation,” IEEE Transactions on Power Systems, vol. 22,
no. 1, pp. 213–221, 2007.
[14] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed
optimization and statistical learning via the alternating direction method of
multipliers,” Foundations and Trends R in Machine Learning, vol. 3, no. 1,
pp. 1–122, 2011.
[15] H. Zhu and G. Giannakis, “Power system nonlinear state estimation using
distributed semidefinite programming,” IEEE Journal on Selected Topics in
Signal Processing, vol. 8, no. 6, pp. 1039–1050, 12 2014.
[16] J. Matamoros, A. Tsitsimelis, M. Gregori, and C. Antn-Haro, “Multiarea
state estimation with legacy and synchronized measurements,” in 2016 IEEE
International Conference on Communications, May 2016, pp. 1–6.

R EFERENCES
[1] A. Abur and A. Expósito, Power System State Estimation: Theory and
Implementation, ser. Power Engineering. Taylor & Francis, 2004.
[17] V. Kekatos and G. B. Giannakis, “Distributed robust power system state
estimation,” IEEE Transactions on Power Systems, vol. 28, no. 2, 2013.
[18] X. Li and A. Scaglione, “Robust decentralized state estimation and tracking
for power systems via network gossiping,” IEEE Journal on Selected Areas in
Communications, vol. 31, no. 7, pp. 1184–1194, 2013.
[19] L. Xie, D. H. Choi, S. Kar, and H. V. Poor, “Fully distributed state estimation
for wide-area monitoring systems,” IEEE Transactions on Smart Grid, vol. 3,
no. 3, pp. 1154–1169, Sept 2012.
[20] A. Gómez-Expósito, A. de la Villa Jaén, C. Gómez-Quiles, P. Rousseaux, and
T. Van Cutsem, “A taxonomy of multi-area state estimation methods,” Electric
Power Systems Research, vol. 81, no. 4, pp. 1060–1069, 2011.
[21] J. Pearl, Probabilistic Reasoning in Intelligent Systems: Networks of Plausible
Inference. San Francisco, USA: Morgan Kaufmann Publishers Inc., 1988.
[22] C. M. Bishop, Pattern Recognition and Machine Learning. Springer, 2006.
[23] Y. Weiss and W. T. Freeman, “On the optimality of solutions of the maxproduct belief-propagation algorithm in arbitrary graphs,” IEEE Transactions
on Information Theory, vol. 47, no. 2, pp. 736–744, 2001.
[24] D. Koller and N. Friedman, Probabilistic Graphical Models: Principles and
Techniques. MIT Press, 2009.
[25] 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.
[26] Y. Weng, R. Negi, and M. Ilic, “Graphical model for state estimation in electric
power systems,” in Smart Grid Communications, 2013 IEEE International
Conference on, Oct 2013, pp. 103–108.
[27] T. Sui, D. E. Marelli, and M. Fu, “Convergence analysis of gaussian belief
propagation for distributed state estimation,” in Decision and Control (CDC),
2015 IEEE 54th Annual Conference on. IEEE, 2015.
[28] M. Cosovic and D. Vukobratovic, “Fast real-time DC state estimation in
electric power systems using belief propagation,” 2017. [Online]. Available:
http://arxiv.org/abs/1705.01376
[29] F. R. Kschischang, B. J. Frey, and H. A. Loeliger, “Factor graphs and the sumproduct algorithm,” IEEE Transactions on Information Theory, vol. 47, no. 2,
pp. 498–519, Feb 2001.
[30] H. A. Loeliger, J. Dauwels, J. Hu, S. Korl, L. Ping, and F. R. Kschischang,
“The factor graph approach to model-based signal processing,” Proc. of the
IEEE, vol. 95, no. 6, pp. 1295–1322, 2007.
[31] A. Monticelli, State Estimation in Electric Power Systems: A Generalized
Approach, ser. Kluwer international series in engineering and computer
science. Springer US, 1999.
[32] D. Barber, Bayesian Reasoning and Machine Learning.
Cambridge
University Press, 2012.
[33] A. Wood and B. Wollenberg, Power Generation, Operation, and Control, ser.
A Wiley-Interscience publication. Wiley, 1996.
[34] G. Elidan, I. McGraw, and D. Koller, “Residual Belief Propagation: Informed
Scheduling for Asynchronous Message Passing,” ArXiv e-prints, Jun. 2012.
[35] Y. Weiss and W. T. Freeman, “Correctness of belief propagation in gaussian
graphical models of arbitrary topology,” in Advances in neural information
processing systems, 2000, pp. 673–679.
[36] B. L. Ng, J. Evans, and S. Hanly, “Distributed downlink beamforming in
cellular networks,” in 2007 IEEE International Symposium on Information
Theory, June 2007, pp. 6–10.
[37] C. Fan, X. Yuan, and Y. J. Zhang, “Scalable uplink signal detection in c-rans
via randomized gaussian message passing,” IEEE Transactions on Wireless
Communications, vol. PP, no. 99, pp. 1–1, 2017.
[38] M. Pretti, “A message-passing algorithm with damping,” Journal of Statistical
Mechanics: Theory and Experiment, vol. 2005, no. 11, 2005.
[39] P. Hansen, V. Pereyra, and G. Scherer, “Least squares data fitting with
applications.” Johns Hopkins University Press, 2012, ch. 9, p. 166.
[40] M. Cosovic and D. Vukobratovic, “Distributed gauss-newton method for ac
state estimation: A belief propagation approach,” in 2016 IEEE International
Conference on Smart Grid Communications (SmartGridComm), Nov 2016,
pp. 643–649.
[41] M. K. Enns, W. F. Tinney, and F. L. Alvarado, “Sparse matrix inverse factors
[power systems],” IEEE Transactions on Power Systems, vol. 5, no. 2, 1990.
[42] F. L. Alvarado, “Computational complexity in power systems,” IEEE
Transactions on Power Apparatus and Systems, vol. 95, no. 4, July 1976.
[43] D. Bickson, “Gaussian Belief Propagation: Theory and Aplication,” ArXiv eprints, Nov. 2008.

