JOURNAL OF LATEX CLASS FILES, VOL. 13, NO. 9, SEPTEMBER 2014

1

From State Estimation to Network Reconstruction

arXiv:1701.09084v1 [] 28 Jan 2017

Farnaz Basiri, Jose Casadiego, Marc Timme, Dirk Witthaut

Abstract—We develop methods to efficiently reconstruct the
topology and line parameters of a power grid from the measurement of nodal variables. We propose two compressed sensing
algorithms that minimize the amount of necessary measurement
resources by exploiting network sparsity, symmetry of connections and potential prior knowledge about the connectivity. The
algorithms are reciprocal to established state estimation methods,
where nodal variables are estimated from few measurements
given the network structure. Hence, they enable an advanced
grid monitoring where both state and structure of a grid are
subject to uncertainties or missing information.
Index Terms—State estimation, Power system control, SCADA
systems, Compressed sensing

I. I NTRODUCTION
The secure and reliable operation of complex power grids
requires a precise knowledge of the grid topology and the state
of connected generation and transmission elements. A central
tool is ‘state estimation’, i.e the estimation of the current
state of the nodes from few, possibly noisy measurements
and the grid data. Various methods and algorithms have been
developed in this field, see e.g. [1] for a review. A prime
example is the estimation of the voltage phase angles from
power flow measurements (see e.g. [2] for basic textbook
examples). State estimation algorithms typically assume that
the network structure, i.e. the location and parameters of
transmission lines and transformers, are known. However,
parts of the information might be lacking or uncertain, for
example for contingencies or attacks, such that the question
arises: Is it possible to reconstruct the network structure from
local measurements only?
In this article we develop two algorithms that enable the
faithful reconstruction of the network structure or parts of it
from the measurement of nodal variables. Special attention
is paid to the efficiency of these methods: We show how
to reduce the number of measurements using methods from
compressed sensing exploiting the structural properties of
power grids. We here focus on the DC approximation to
formulate the fundamental ideas of network reconstruction
most clearly.
The presented algorithms can find several applications in
the monitoring and operation of (smart) power grids. While
it is unlikely that the entire network structure is unknown,
situations regularly occur where parts of the information
F. Basiri is with the Forschungszentrum Jülich, Institute for Energy and
Climate Research - Systems Analysis and Technology Evaluation (IEK-STE),
52428 Jülich.
J. Casadiego and M. Timme are with the Max Planck Institute for Dynamics
and Self-Organization (MPIDS), 37077 Göttingen, Germany.
D. Witthaut is with the Forschungszentrum Jülich, Institute for Energy and
Climate Research - Systems Analysis and Technology Evaluation (IEK-STE),
52428 Jülich, Germany and the Institute for Theoretical Physics, University
of Cologne, 50937 Köln, Germany.

is lacking. Manual switches still exist in many distribution
grids, transmission lines can get lost in contingency cases and
targeted attacks can cause entire regions of a grid to collapse.
We present an example where the connectivity in two regions
of the grid is a priori unknown, but can be reconstructed from
a single snapshot of the nodal variables. In future power grids,
the network’s structure will become even more variable: The
transmission of high-voltage direct current (HVDC) lines in
hybrid power grids can be actively controlled and the effective
line parameters can be regulated with flexible alternating current transmission system (FACTS) devices. Then it becomes
essential to include lacking information or uncertainties of
the grid topologies in any grid monitoring tool. The methods
presented in this article shall contribute to the development of
such advanced monitoring algorithms [3]–[5].
II. T HE DC APPROXIMATION
The DC approximation describes the flow of real power
flow in AC power grids in a linearized way. It is derived
from the full nonlinear load flow equations under the following
additional simplifying assumptions [2], [6]–[8]:
• All bus voltage magnitudes |Vs | are close to 1 p.u. and
approximated to be exactly |Vs | = 1 p.u.
• Voltage angle differences across branches are small, such
that sin(ϕs − ϕr ) ≈ ϕs − ϕr and cos(ϕs − ϕr ) ≈ 1.
• Branches can be considered lossless. In particular, branch
resistances and charging capacitances are negligible.
The approximate real power flow from node s to node r over
a transmission line with reactance Xsr is then given by
|Vs ||Vr |
−1
sin(ϕs − ϕr ) ≈ Xsr
(ϕs − ϕr )
(1)
Xsr
where ϕs is the phase angle at bus s. The real power balance
at the node s reads
N
N
X
X
−1
Ps =
Psr =
Xsr
(ϕs − ϕr ).
(2)
Psr =

r=1

r=1

For notational convenience we summarize these equations
in matrix form. We define the vector of all voltage angles and
power injections,
ϕ = (ϕ1 , . . . , ϕN )T ∈ RN , P

= (P1 , . . . , PN )T ∈ RN ,

where N is the number of buses or nodes in the grid and the
superscript ‘T ’ denotes the transpose of a vector or matrix.
The DC approximation then yields the condition
Bϕ = P ,

(3)
N ×N

where the nodal susceptance matrix B ∈ R
has elements
 X
N

X −1 if k = n;
j=1 nj
Bnk =
(4)
 −X −1
if k 6= n.
nk

JOURNAL OF LATEX CLASS FILES, VOL. 13, NO. 9, SEPTEMBER 2014

2

The matrix B is a Laplacian matrix, which has one zero eigenvalue with eigenvector (1, 1, . . . , 1)T [9]. This eigenvector
represents a global shift of all voltage angles without physical
significance. To avoid this complication, one often fixes the
voltage angle at a reference bus and excludes this bus from
the calculation, i.e. the corresponding row and column will be
missing in B. We finally note that mathematically equivalent
models of flow are used to describe hydraulic networks [10]
or vascular networks of plants [11].
III. S TATE E STIMATION
State estimation is a mathematical procedure by which the
state of an electric power system is extracted from a set
of measurements [1]. In general, any measurement can be
expressed as a function of the system state. Let
zi = hi (x)

(5)

denote a measured quantity, where x is the system state and
hi is a function specific to the measured quantity zi . Assume
that M measurements are taken. Then, all measurements can
be written in compact form
z = h(x)

(6)

where
x ∈ RN ×1 is the system state,
z ∈ RM ×1 is a vector of measured quantities and
h is a vector function, i.e. a mapping RN → RM
Typically more measurements than the number of state variables to be determined are taken, i.e. M > N , such that
equation (6) is overdetermined. It is then solved in a leastsquares fashion, i.e. one calculates x such that kz − h(x)k22 is
minimized. For underdetermined equations (6) typically many
solutions exist.
It is often helpful to work with a simplified DC approximation model for the measurement equations in analyzing
the inherent limitations of various methods related solely to
the measurement configuration [12]. The state of the grid is
described by the voltage phase angles x = ϕ and shall be
estimated from measurements. As explicated above, the real
power flow measured from bus s to r can be approximated by
−1
Psr = Xsr
(ϕs − ϕr ) + esr

N
X

Psr + e0s

Let us come back to the defining equation for the DC
approximation
Bϕ = P .
(9)
In state estimation, we assume that the nodal phase angles
ϕ are difficult to measure and shall be estimated from the
knowledge of the two other quantities in this equation, B
and P . But we can also reverse the problem and ask what
happens if we don’t know the nodal susceptance matrix B.
Can we efficiently reconstruct its entries from measurements
of the two remaining quantities P and ϕ? This problem may
appear hypothetical at a first glance, as the network structure
and parameters are generally known. But there are situations
where at least parts of the network structure are unknown:
Switches can be open or closed, transmission lines can undergo
failures, or plans may be inaccurate. Even more, in times
of war or terrorism we may be extremely unsure about the
physical integrity of parts of the grid. Can we reconstruct the
structure of the grid from measurements only?
In the following we develop a mathematical theory of
network reconstruction. To introduce the method, we start
from the hypothetical situation that B is completely unknown,
whereas P and ϕ can be measured perfectly. We then study
more realistic cases, where only parts of the grid topology are
unknown.
We first notice that a single measurement of P and ϕ will
obviously not be enough to reconstruct the entire grid topology. However, we may repeat our measurement at different
times t1 , t2 , . . . , tM such that we obtain a large number of
conditions on the entries of B of the form
X
Brk ϕk (tm ) = Pr (tm )
(10)
k

which hold for all rows r = 1, . . . , N and all measurement
time steps m = 1, . . . , M . Furthermore we know that the rowsums of the nodal susceptance matrix (4) vanish,
X
Brk = 0.
(11)
k

In total we thus have N × (M + 1) conditions which we can
use to reconstruct the unknown entries of B.
A. Row-wise reconstruction

(7)

where ϕs is the phase angle at bus s and esr is a measurement
error. Similarly, a power injection measurement
Ps =

IV. N ETWORK R ECONSTRUCTION

(8)

r=1

at a given bus s can be expressed as a sum of flows along
incident branches to that bus. The measurement vector z
thus consists of a subset of the power flows Psr and power
injections Ps . In the DC approximation the vector function
h are all linear, which strongly simplifies the problem such
that it is especially suitable to introduce new concepts such as
network reconstruction.

In a first approach we aim to reconstruct the nodal admittance matrix row-by-row (cf. [13]–[16]). To simplify notation
we collect all measurement conditions (10) and the condition
(11) for a given row r and rewrite them in matrix form





Pr (t1 )
ϕ1 (t1 ) · · · ϕN (t1 ) 
B
r1
 . 
 ..

..
 .
  ..   .. 
.
 (12)

 .  = 
Pr (tM )
ϕ1 (tM ) · · · ϕN (tM )
BrN
0
1
···
1
| {z }
|
{z
} =:B r
{z
}
|
=:Φ

N ×1

=:P r

Here, B r ∈ R
denotes the transpose of the rth row of the
nodal admittance matrix B, which we want to reconstruct. The
matrix Φ ∈ R(M +1)×N summarizes all measurement results

JOURNAL OF LATEX CLASS FILES, VOL. 13, NO. 9, SEPTEMBER 2014

3

for the voltage phase angles at all time steps, while the vector
P r ∈ R(M +1)×1 summarizes the values of the power injection
only for node r. The last row of the matrix equation represents
the condition (11) for the row sum.
Now we have to distinguish, whether equation (12) is
(over)determined or underdetermined. If the number of linearly independent equations is larger or equal to N we
generally have enough information to directly compute B r .
In particular there is exactly one solution if rank(Φ̄|P¯r ) =
rank(Φ̄) = N , which can be calculated using Gaussian
elimination. Nevertheless, matrix inversion can be numerically
ill-conditioned, especially for systems with large N [17].
If the data includes some measurement noise, the system
of equations (12) can be overdetermined when measuring
M + 1 > N times such that it is solvable in a least squares
fashion, i.e. we have to determine the vector B r ∈ RN ×1
minimizing the 2-norm
kΦB r − P r k22 =

M
+1
X

(ΦB r − P r )2k .

(13)

k=1

To efficiently reconstruct the grid topology, we want to rely on
as few measurements as possible. If M +1 < N , we generally
have an underdetermined set of equations which admits many
possible solutions. Is it possible to obtain the correct solution
also in this case? Do we have more information about B r
which we can exploit? Indeed, we know that a power grid
is typically very sparse – a single substation is connected to
only few other substations. Thus we choose the one solution to
equation (12) which minimizes the number of non-zero entries
of B r . Unfortunately, the direct minimization is computationally hard in general. In 2006 Candes, Romberg, Tao and
Donoho showed that an efficient reconstruction is nevertheless
possible using a convex surrogate for sparsity: the 1-norm [18],
[19]. Under weak conditions, the correct sparse solution can
be calculated efficiently by minimizing the 1-norm
kB r k1 =

N
X

|B rk |

(14)

k=1

subject to the constraint (12). This problem can be mapped to
a linear program which we explicate in the following section.
B. Partial reconstruction
In applications we will rarely encounter the situation that the
structure of a power grid is completely unknown. In a typical
application we know the value of Brs of all transmission lines
(s, r) in one part of the grid and have to reconstruct only the
remaining entries of the matrix B – for instance we might
want to monitor the position of manual switches at remote
places of the grid.
To keep track of our knowledge about the network, we
define the matrix K ∈ RN ×N with entries

1 if the value of Brs is known for line (r, s)
Krs :=
0 otherwise.
(15)
We again proceed row by row. For the rth row we reduce
equation (12) to
red
Φred
= P red
(16)
r Br
r ,

where



P red
r


ϕc (t1 )
N
 .. 
X


= Pr −
Brc Krc  . 

ϕc (tM )
c=1
1

(17)

and Φred
is the submatrix of Φr obtained by deleting all
r
columns c for which Krc = 1 and B red
is the submatrix of
r
B r obtained by deleting all rows r0 for which Kr,r0 = 1. The
dimension of the reduced linear system of equations (16) is
smaller than in the original problem (12), such that the correct
solution can typically be found from fewer measurements.
C. Reconstruction from power flow measurements
The reconstruction scheme introduced above requires the
knowledge of the voltage phase angle at all nodes of the grid.
This data can in principle be obtained using phasor measurement units (PMUs), but these are typically very expensive.
The measurement of other quantities such as real power flows
is typically much simpler.
Fortunately, we usually do not have to reconstruct the entire
grid in most cases. In a typical application we know a lot
about the grid and have to reconstruct only parts of the matrix
B. Then the measurement of power injections and real power
flows of a certain subset of transmission lines can be enough
to perform the reconstruction of the remaining data using
the methods described above. Technically, this amounts to
combining methods of classical state estimation and network
reconstruction.
So assume that we can measure the power injection Pn
for all nodes n ∈ {1, . . . , N } at the time steps t1 , . . . , tM .
In addition we have knowledge about the transmission line
parameters, in particular Xsr , and measure the real power flow

−1
Psr (tm ) = Xsr
ϕr (tm ) − ϕs (tm )
(18)
for a subset L of all transmission lines at all time steps
t1 , . . . , tM . If this set of equations is fully determined or even
overdetermined we reconstruct the state vector ϕ(tm ) as in the
classical state estimation problem described in section III. As
before the resulting estimates for the nodal phase angles for
all time steps t1 , . . . , tM are then summarized in the matrix
Φ. In addition we have partial knowledge about the nodal
susceptance matrix, in particular we already know the entries
−1
Brs = Bsr = −Xsr

(19)

for all (s, r) ∈ L. This information can then be used to
reconstruct the remaining entries of B as described in the
previous section.
V. I MPLEMENTATION
A. Least-squares solution
If we have many measurements or prior knowledge available, then the system of equations (12) can be be overdetermined. In this case equation must be solved in a least squares
fashion, i.e. we have to solve
min kΦB r − P r k22 ,
Br

(20)

JOURNAL OF LATEX CLASS FILES, VOL. 13, NO. 9, SEPTEMBER 2014

4

which leads to [15]

Start

B r = (ΦT Φ)−1 ΦT P r

(21)
Measurement of Φ and 𝑃

provided that Φ is full rank such that the inverse exists. Leastsquares solutions are already implemented in many numerical
solvers, for instance in the MATLAB function mldivide.

Add constraint ∑𝑐 𝐵𝑟𝑐 = 0

B. Minimizing the 1-norm

no

𝑟𝑘(Φ) < 𝑁 ?

The key to an efficient reconstruction of B from an underdetermined system of equations is the minimization of the
1-norm [18], [19]. This problem can be mapped to a linear
program which can be solved efficiently, i.e. the computation
time scales at most polynomially with the input length.
So we consider the linear system of equations (12) and
assume that it is solvable and underdetermined. The solutions
span an affine subspace of of RN with dimension D =
N − rank(Φ). All solutions can be written as
+ W y,
B r = B (sp)
r

Linear Sy
underdet

yes

𝑟 = 1
(sp)

Calculate Specific solution 𝐵𝑟
Solve linear program

(22)

(sp)

min𝑦 ‖𝑊𝑦 + 𝐵𝑟

‖

1

N ×D

where the columns of the matrix W ∈ R
form a basis
for the kernel (nullspace) of Φ and y ∈ RD is a vector of
is a specific solution to equation (12); in
parameters. B (sp)
r
the implementation it is obtained using the MATLAB function
mldivide. We now search for the vector y such that the 1norm (14) assumes its minimum. This optimization problem
is rewritten as a linear program [20]
min 1T s
s

+ Wy  s
B (sp)
r
,
B (sp)
+ W y  −s
r

such that

min f T x
x

such that

Ax  b.

(sp) 𝑇

Set 𝐵(𝑟, : ) = [𝑊𝑦 + 𝐵𝑟

(24)

using the further auxiliary variables
 


y
W
−1l
x=
, A=
,
s
−W −1l
 (sp) 
 
Br
0
b=
,
f
=
,
1
−B (sp)
r
where 1l ∈ RN ×N is the identity matrix and 0 ∈ RD is a
vector of zeros.
C. Row-Wise Reconstruction
We first describe the implementation of the row-wise reconstruction of B from the measurement data without any prior
knowledge about the entries of B, whose basic idea has been
outlined in section IV-A. The algorithm to solve this problem
is illustrated by the flow-chart in Fig. 1. We implement this
algorithm in Matlab, the central part of the program code is
listed in Fig. 2.
First data is collected and the matrix Φ and the vectors
P r are formed. Then the problem is solved depending on the
rank of the matrix Φ. If the system is (over)determined, we
directly solve it for B using the Matlab-function mldivide.

]

𝐵 = Φ\𝑃𝑟

𝑟 = 𝑟+1
yes

𝑟≤𝑁
no
Return 𝐵

(23)

where s ∈ RN is an auxiliary variable, 1 ∈ RN is a vector
of ones and  and  denotes entry-wise comparison. In the
standard form of many commercial solvers the optimization
problem reads

underdet

Calculate null space 𝑊of Φ

Fig. 1. Flow chart for row-wise reconstruction. MATLAB code for the central
part is listed in Fig. 2.

If the system is underdetermined we proceed row-by-row and
reconstruct B r as described in section V-B using the function
null to compute the matrix W and the function linprog
to solve the linear program (24).
D. Iterative reconstruction
In practice we typically have some prior knowledge about
the entries of the matrix B as outlined in section IV-B. An
algorithm to perform the reconstruction exploiting the prior
knowledge is illustrated by the flow-chart in Fig. 3
Most interestingly, we can also use this algorithm to greatly
improve the convergence of the reconstruction algorithm. For
intermediate values of M we typically face the situation that
some rows are successfully reconstructed while others are
not. If we have not reconstructed B successfully, but gained
further knowledge, we may use it in the following. Having
successfully reconstructed the rth row of B, we also know the
entries Br0 ,r = Br,r0 for all other rows r0 due to the symmetry
of the matrix. We thus propose to perform the reconstruction
iteratively. In each step of the reconstruction process we use
the initial knowledge about the entries of the matrix B and
also the knowledge gained in previous steps.
The algorithm shown in Fig. 3 starts with the input of the
measurement results and the prior knowledge of the system

JOURNAL OF LATEX CLASS FILES, VOL. 13, NO. 9, SEPTEMBER 2014

5

Start

% Case 1: Overdetermined system:
if rank(Phi) >= N
Breconstruct = Phi \ Pt;

Input measurement results Φ, 𝑃𝑟
and prior knowledge 𝐵, 𝐾 new
Set 𝐾 = 𝐾 new , 𝑟 = 0

% Case 2: Underdetermined system:
else

Proceed to next row 𝑟 = 𝑟 + 1

% Calculate nullspace of Phi
W = null(Phi);
D = size(W,2);

Set Φ𝑟red , 𝑃𝑟red , 𝑁𝑟red = 𝑁 − Σ𝑐 𝐾𝑟𝑐
yes

𝑁𝑟red = 0 ?
no

% treat all rows separately:
for r=1:N

𝑟(Φred ) < 𝑁red

Linear Sy
underdet

underdet

no

yes
Solve linear program
min ‖𝐵𝑟red ‖1 s.t. Φred 𝐵𝑟red = 𝑃𝑟red

% one specific solution
Br_sp = Phi \ Pt(:,r);
no

%
f
A
b
y
y

Choose sparse solution
= [zeros(D,1); ones(N,1)];
= [W, -eye(N); - W, -eye(N)];
= [-Br_sp; +Br_sp];
= linprog(f,A,b);
= y(1:D);

𝐵𝑟red = Φ𝑟red \𝑃𝑟red

Successful?
new
Set 𝐾𝑟𝑐
= 1 ∀𝑐 = 1, … 𝑁
Add results 𝐵𝑟red to matrix 𝐵

no

𝑟 = 𝑁?

yes
𝐾

% Write result to matrix B
Breconstruct(r,:) = (W*y+Br_sp)’;

new

=𝐾?

yes

no
no

new
𝐾𝑟𝑐

= 1 ∀ 𝑟, 𝑐no
?

yes

Output
Warning

end
Return 𝐵

end
Fig. 2.

MATLAB code for row-wise reconstruction

which is encoded in the matrices B and K. It solves the
reconstruction problem iteratively, processing row-by-row in
each iteration. For each row, we first calculate the reduced
as defined in section IV-B. Then we
and P red
matrices Φred
r
r
attempt to reconstruct the vector B red
as described above: If
r
the problem is (over)determined we solve the linear system
of equations using the Matlab-function mldivide, if it is
underdetermined we minimize the 1-norm kB red
r k1 subject to
the constraint (16). If the reconstruction was successful, we
add the results to the matrix B and update the knowledge
matrix K. To facilitate the bookkeeping we define two matrices K and K new : We set K = K new at the beginning
of each iteration and only modify K new during the step. The
iteration stops when the reconstruction has been successfully
new
completed, i.e. Krc
= 1 for all r, c = 1 . . . , N , or when
no further progress has been made, i.e. K new has not been
modified during the last step.
One unsolved problem remains: If we do not know the
matrix B a priori, how can we know if the reconstruction of
the rth row has been successful? In the underdetermined case
we cannot decide whether the reconstructed values Brc are
definitely correct – but we can decide if they are reasonable in
terms of the connectivity of the grid. Power grids are generally
very sparse: A single substation is connected to only few other

Fig. 3.

Flow chart for the iterative reconstruction algorithm.

substations. Hence we seek for a solution which is most sparse,
i.e. which minimizes
kB r k0 = number of non-zero entries of B r .

(25)

However, we cannot minimize kB r k0 directly as this is a
computationally hard problem. The ingenious contribution of
Candes, Romberg, Tao and Donoho [18], [19] was the proof
that the most sparse solution can be obtained by minimizing kB r k1 with high probability if a sufficient amount of
measurement data is available. Otherwise we find a solution
which minimizes the 1-norm, but has a large number of nonzero entries kB r k0 . In algorithm in Fig. 3 we thus adopt the
definition that the reconstruction of the rth row is assumed
to be successful if (a) the linear system is overdetermined
or (b) the reconstructed row B r is sufficiently sparse, i.e.
kB r k0 ≤ dmax for some upper limit dmax . A typical value
used in the following is dmax = 15.
VI. A PPLICATIONS AND P ERFORMANCE
A. Reconstruction from time series
We demonstrate the applicability of network reconstruction
from time series measurements for a test grid taken from [21]
illustrated in Figure 4 (a). This data set includes hourly data
for demand and generation Pk (t) for one year and all N =

JOURNAL OF LATEX CLASS FILES, VOL. 13, NO. 9, SEPTEMBER 2014

6

(a)

φk(tm)
[degree]

Pk(tm)
[MW]

(b)

2000

0

1000

−100

0

−200

10

10

−10

−10

0

number of errors

(c)

minimum number of measurements

300

20
time t

m

40
[hour]

60

0

20
time t

m

40
[hour]

60

400

200
100
0
10 0
60

200
100

25

30

35
40
45
50
number of measurements M

55

60

Fig. 4.
Example of the complete grid reconstruction for a test grid with
N = 117 nodes. (a) Test grid from [21]. It is assumed that the network
topology, i.e. the location and parameters of the transmission lines (grey) are
completely unknown. (b) The reconstruction is based on the time series of
the nodal power injections Pk (t) and voltage angles ϕk (t). Two time series
are shown. (c) Performance of the iterative algorithm shown in Fig. 3. Perfect
reconstruction is obtained with M ≥ 55 measurements. We use the test grid
and the load and generation time series of the scenario 2013 from [21]. For
the sake of simplicity we have removed the import/export nodes and rescaled
the generation to exactly match the load.

117 nodes. The nodal voltage angles ϕk (t) are then obtained
by solving the DC approximation (3) for all time steps. The
resulting time series of the power injections and angles are
shown in Figure 4 (b) for two nodes as an example.
The algorithm shown in Figure 3 can now reconstruct
the entire network topology - i.e. all entries of the nodal
susceptance matrix B - from the time series data. To evaluate
the performance of the algorithm we vary the amount of input
data. For each value of M we run the algorithm and count the
number of successfully reconstructed entries and the number
of reconstruction errors. To account for small numerical errors,
we say that an entry (r, k) of the nodal susceptance matrix is
successfully reconstructed if
reconstructed
true
|Br,k
− Br,k
| < ,

(26)

where the numerical tolerance is chosen as  = 10−3 ×
maxij |Bij |.
Figure 4 (c) shows that Mmin = 55 measurements are
sufficient to faithfully reconstruct the entire nodal susceptance matrix B. Most importantly, a faithful reconstruction is
possible in the strongly underdetermined case: The minimum
number of measurements Mmin = 55 is less than half of the
matrix dimension N = 117. This is possible because we can

10 1

10 2

10 3

10 2

10 3

(b)

40
20
0
10 0

300

0
20

(a)

10 1

grid size N

Fig. 5.
Performance of (a) the row-wise reconstruction algorithm shown
in Fig. 1 and (b) the iterative algorithm shown in Fig. 3. We have tested the
algorithms for several test grids from [22], [23] for synthetic power injection
data, choosing either the ϕk (blue circles) or the Pk (red crosses) randomly
(see text for details). Shown is the minimum number of measurements M
needed to obtain a faithful reconstruction as a function of the grid size N .
Results have been averaged over 10 random realizations.

exploit the sparsity of the matrix B as an additional structural
information.
B. Minimum Measurement Requirements and Scaling
The algorithms presented in figures 1 and 3 allow for a
reconstruction of the network topology also in the underdetermined case. But how many measurements Mmin are needed
for a faithful reconstruction and how does this number scale
with the grid size N ?
To systematically study the efficiency of the reconstruction
algorithms, we consider various test grids of different size
taken from [22], [23] and use randomized synthetic data for
the time series. We consider two different types of time series
data. First we draw the voltage phase angles ϕk (t) uniformly
at random from the interval [−π/8, +π/8] (type I). This
scenario is not realistic but close to the original mathematical
work [18], [19], where random Gaussian sampling vectors are
considered. Second, we choose the power injections Pk (t)
at random from a normal distribution with mean zero and
standard deviation as in the original test grid (type II). For
both types and both reconstruction algorithms (row-wise and
iterative), we run the reconstruction algorithm as a function
of the number of measurements M taken into account and
determine the minimum number of measurements Mmin for
a faithful reconstruction, i.e. zero reconstruction errors. Each
numerical experiment is repeated 10 times.
First, we observe that the iterative algorithm allows for a
faithful reconstruction from underdetermined input data for all
grids under consideration and both data types. The minimum
number of measurements Mmin is much smaller than the grid
size N in all cases. The row-wise algorithm performs less well
as it does not exploit the symmetry of the matrix B.

JOURNAL OF LATEX CLASS FILES, VOL. 13, NO. 9, SEPTEMBER 2014
power
flow
[MW]

(a)

7

prior knowledge

(b)
1

?

29.6
22.2

bus number

37
10

20

14.8

?

7.4

30
1

10

0

(c)

phase
angles
φk [°]

20
bus number

30

reconstructed network

(d)
1

4.8
3.6

bus number

6
10

20

2.4
1.2

30
1

10

0

20
bus number

30

Fig. 6. Reconstruction of missing information from flow measurements (a)
The IEEE 30-bus test grid. Real power flows are measured on the colored
lines. In two regions the grid topology is assumed to be unknown. (b) The
prior knowledge matrix K. Grey shading indicates that the respective line
parameter Bjk is known. The white region is completely unknown. (c) The
voltage phase angles ϕk are reconstructed from the power flow measurements
using standard state estimation methods. (d) The completely grid topology is
perfectly reconstructed using the algorithm shown in Fig. 3.

Second, efficient reconstruction is possible with both algorithms for the type I input data. The theory of compressed
sensing shows that a reconstruction of sparse vectors is possible with high probability from a number of measurements
scaling only logarithmically with the problem dimension N
if some conditions are satisfied [18], [19]. This logarithmic
scaling is well confirmed by the results of our numerical
experiments. For the iterative algorithm the scaling seems to
be even slower.
However, the proofs for the efficiency of compressed sensing depend on two features of the sampling vectors: isotropy
and incoherence. Roughly speaking, isotropy means that all
directions in the vector space are sampled equally well.
Incoherence guarantees that the redundancy in the information
gained by two measurements is small. Both features are no
longer guaranteed for type II input data. In fact, we observe
that the row-wise reconstruction algorithm does not perform
well in this case. The minimum number of measurements
Mmin is only slightly smaller than the problem dimension N
such that the iterative algorithm performs much better.
C. Reconstruction from flow measurements
The iterative algorithm allows to take into account prior
knowledge to efficiently reconstruct the missing information
about the grid topology. Such a problem can arise in practice
when a larger part of the grid is subject to damages or attacks.
An example of such a situation is depicted in figure 6 (a). It is
assumed that we have no reliable information about the grid
in two areas indicated by the question marks. However, we
do have information about the connections of all remaining
nodes and we assume that we can measure the real power
flow along the colored solid lines. The prior knowledge of the

grid topology is encoded in the matrix K, which is illustrated
in figure 6 (b). White entries indicate where we do not know
the entries of the matrix B.
The missing information is found by combining classical
state estimation and network reconstruction. First, the voltage
phase angles for all nodes of the network are reconstructed as
described in section III, the results being shown in figure 6 (c).
Second, the iterative reconstruction algorithm shown in figure
3 is applied using the measured power injections, the prior
knowledge B and K as well as the estimated voltage phase
angles as input data. The algorithm then faithfully reconstructs
all the missing information about the grid topology. The
−1
reconstructed values of the line susceptances Xsr
shown in
figure 6 (d) exactly match the true values. Most importantly,
the full reconstruction is possible already for M = 1 in this
case, i.e. with a single snapshot of the power injections P ,
making use of the prior knowledge.
VII. C ONCLUSION AND O UTLOOK
We have introduced two algorithms to reconstruct the structure of a power grid from nodal measurements only. One
algorithm is straightforward and row-based, reconstructing the
lines of the grid node by node. The second is iterative and
takes into account both prior knowledge about the presence of
absence of lines as well as the knowlegde generated by that
algorithm during previous steps. We have demonstrated how
these algorithms can be used to reconstruct the entire network
structure from time series or missing information on the grid
topology from a single snapshot.
The presented algorithms exploit several structural properties of power grids to reduce the number of necessary
measurements. Power grids are typically very sparse, i.e.
each substation is connected to only few other substations.
Hence, methods from compressed sensing can be used which
allow for a faithful reconstruction also in the underdetermined
case. The iterative reconstruction algorithm makes use of the
symmetry of the nodal susceptance matrix. Loosely speaking
the algorithm solves the simple parts of the reconstruction
problem first and than uses the gained information for the
remaining parts. We have shown that this trick leads to a vast
reduction of the measurement resources.
The present work presents a step towards hybrid state
estimation/network reconstruction algorithms. In future smart
grids one can face the situation where both nodal variables
and the grid topology are subject to uncertainties or lacking
information. A pure state estimation starting from perfect
knowledge of the network structure is then no longer sufficient.
Advanced grid monitoring algorithms must be able of coping
with missing information in both state and structure.
ACKNOWLEDGEMENTS
We thank D. Gross and R. Kueng for inspiring discussions. We gratefully acknowledge support from the Helmholtz
Association (via the joint initiative “Energy System 2050
- A Contribution of the Research Field Energy” and the
grant no. VH-NG-1025 to D.W.) and the Federal Ministry of
Education and Research (BMBF grant nos. 03SF0472B and
03SF0472E to D.W. and M.T.).

JOURNAL OF LATEX CLASS FILES, VOL. 13, NO. 9, SEPTEMBER 2014

R EFERENCES
[1] A. Meliopoulos, B. Fardanesh, and S. Zelingher, “Power system state
estimation: modeling error effects and impact on system operation,” in
Proceedings of the 34th Annual Hawaii International Conference on
System Sciences (HICSS-34). IEEE, 2001, p. 2012.
[2] A. J. Wood, B. F. Wollenberg, and G. B. Sheblé, Power Generation,
Operation and Control. New York: John Wiley & Sons, 2014.
[3] Y. Yuan, G.-B. Stan, S. Warnick, and J. Goncalves, “Robust dynamical
network structure reconstruction,” Automatica, vol. 47, no. 6, pp. 1230
– 1235, 2011.
[4] J. Yu, G. G. Karady, and L. Gu, “Applications of embedded hvdc in
power system transmission,” in 2012 Power Engineering and Automation
Conference, 2012, pp. 1–6.
[5] S. Bolognani, R. Carli, and M. Todescato, “State estimation in power
distribution networks with poorly synchronized measurements,” in 53rd
IEEE Conference on Decision and Control, 2014, pp. 2579–2584.
[6] J. J. Grainger and W. D. Stevenson Jr., Power System Analysis. New
York: McGraw-Hill, 1994.
[7] K. Purchala, L. Meeus, D. V. Dommelen, and R. Belmans, “Usefulness
of dc power flow for active power flow analysis,” in IEEE Power
Engineering Society General Meeting, 2005, 2005, pp. 454–459 Vol.
1.
[8] D. Van Hertem, J. Verboomen, K. Purchala, R. Belmans, and W. L.
Kling, “Usefulness of DC power flow for active power flow analysis
with flow controlling devices,” in 8th IEE International Conference on
AC DC Power Transmission (ACDC 2006). IEE, 2006, pp. 58–62.
[9] M. Newman, Networks: An Introduction. New York: Oxford University
Press, 2010.
[10] N. Hwang and R. Houghtalen, Fundamentals of Hydraulic Engineering
Systems. Upper Saddle River, NJ: Prentice Hall, 1996.
[11] E. Katifori, G. J. Szöllösi, and M. O. Magnasco, “Damage and Fluctuations Induce Loops in Optimal Transport Networks,” Phys. Rev. Lett.,
vol. 104, p. 48704, 2010.
[12] A. Abur and A. G. Exposito, Power system state estimation: theory and
implementation. Boca Raton: CRC Press, 2004.
[13] M. Timme, “Revealing Network Connectivity from Response Dynamics,” Phys. Rev. Lett., vol. 98, no. 22, p. 224101, 2007.
[14] S. G. Shandilya and M. Timme, “Inferring network topology from
complex dynamics,” New J. Phys., vol. 13, no. 1, p. 013004, 2011.
[15] M. Timme and J. Casadiego, “Revealing networks from dynamics: an
introduction,” J. Phys. A Math. Theor., vol. 47, no. 34, p. 343001, 2014.
[16] X. Han, Z. Shen, W.-X. Wang, and Z. Di, “Robust Reconstruction of
Complex Networks from Sparse Data,” Phys. Rev. Lett., vol. 114, no. 2,
p. 028701, 2015.
[17] J. Stoer and R. Bulirsch, Introduction to Numerical Analysis. New
York: Springer, 1993.
[18] D. L. Donoho, “Compressed sensing,” Inf. Theory, IEEE Trans., vol. 52,
no. 4, pp. 1289–1306, 2006.
[19] E. J. Candes, J. K. Romberg, and T. Tao, “Stable signal recovery from
incomplete and inaccurate measurements,” Comm. Pure Appl. Math.,
vol. 59, no. 8, pp. 1207–1223, 2006.
[20] S. Boyd and L. Vandenberghe, Convex optimization.
Cambridge:
Cambridge University Press, 2004.
[21] H. Barrios, A. Roehder, H. Natemeyer, and A. Schnettler, “A benchmark
case for network expansion methods,” in Proc. IEEE PowerTech 2015
Eindhoven, 2015, pp. 1–6.
[22] R. D. Zimmerman, C. E. Murillo-Sanchez, and R. J. Thomas, “MATPOWER: Steady-State Operations, Planning and Analysis Tools for
Power Systems Research and Education,” IEEE Trans. Power Syst.,
vol. 26, p. 12, 2011.
[23] S. Fliscounakis, P. Panciatici, F. Capitanescu, and L. Wehenkel, “Contingency ranking with respect to overloads in very large power systems
taking into account uncertainty, preventive and corrective actions,” IEEE
Trans. Power Syst., vol. 28, p. 4909, 2013.

Farnaz Basiri Farnaz Basiri received the B.Sc. and M.Sc. degree in electrical
power engineering from the RWTH Aachen University, Aachen, Germany in
2013 and 2015, respectively. Since 2015 she joined the Reaserch Group at
Reaserch Center Jülich, Germany as a Research Associate and is currently
pursuing the Ph.D. degree.

8

Jose Casadiego received his Licentiate degree in Physics from the University
of Carabobo (Venezuela) in 2009. He completed his PhD studies in the
International Max Planck Research School for Physics of Biological and
Complex systems (Göttingen) in 2016 working on model-free approaches for
revealing interactions in complex networks from collective dynamics. He now
works as a Postdoctoral Researcher at the Network Dynamics Group of the
Max Planck Institute for Dynamics and Self-Organization.

Marc Timme studied physics and applied mathematics at the Universities
of Würzburg (Germany), Stony Brook (New York, USA) and Göttingen
(Germany). He holds an M.A. in Physics from Stony Brook and a Doctorate
in Theoretical Physics (Göttingen). After postdoctoral and visiting stays at
Cornell University (New York, USA) and the National Research Center of
Italy (Sesto Fiorentino) he is heading the Max Planck Research Group on
Network Dynamics at the Max Planck Institute for Dynamics and SelfOrganization. He is Adjunct Professor at the University of Göttingen.

Dirk Witthaut received his Diploma (MSc) and PhD in Physics from the
Technical University of Kaiserslautern, Germany, in 2004 and 2007, respectively. He has been a Postodoctoral Researcher at the Niels Bohr Institute
in Copenhagen, Denmark and at the Max Planck Institute for Dynamics and
Self-Organization in Göttingen, Germany and a Guest Lecturer at the Kigali
Institue for Science and Technology in Rwanda. Since 2014 he is heading
a Research Group at Forschungszentrum Jülich, Germany and he is Junior
Professor at the University of Cologne, Germany.

