1

Learning with End-Users in Distribution Grids:
Topology and Parameter Estimation

arXiv:1803.04812v1 [cs.SY] 12 Mar 2018

Sejun Park∗,† , Deepjyoti Deka† , Scott Backhaus† , Michael Chertkov†,§
∗KAIST, Daejeon, Korea
†Los Alamos National Laboratory, Los Alamos, USA
§Skoltech, Moskow, Russia

Abstract—Efficient operation of distribution grids in the smartgrid era is hindered by the limited presence of real-time nodal
and line meters. In particular, this prevents easy estimation of
grid topology and associated line parameters that are necessary
for control and optimization efforts in the grid. This paper
studies the problems of topology and parameter estimation in the
limited observability regime where measurements are restricted
to only the terminal nodes of the grid and all intermediate nodes
are unobserved/hidden. To this end, we propose two algorithms
for exact topology (and impedances) estimation that provably
reconstruct topology using voltage and injection measured only
at the end-users of the distribution grid. The first algorithm
requires time stamped voltage samples, statistics of nodal power
injections and permissible line impedances to recover the true
topology. The second and improved algorithm requires only timestamped voltage/complex power samples to recover both the
true topology and impedances without any additional input (e.g.,
number of grid nodes, statistics of injections at hidden nodes,
permissible line impedances). We discuss the computational and
sample complexity of our proposed algorithms and demonstrate
that topology (and impedance) estimation by our algorithms are
optimal with respect to number of nodal observability required.
We illustrate their performance through numerical experiments
on the IEEE and custom power distribution models.
Index Terms—Distribution networks, Missing data, Power
flows, Sample complexity, Topology learning, Impedance estimation

I. I NTRODUCTION
Distribution grids include the low and medium voltage
transmission lines that help transfer power from the distribution substation to the final consumers. Structurally a majority
of distribution grids are radial in structure. However unlike
traditional passive distribution grids, modern ones have smart
controllable loads, household renewable generators (e.g., solar
panels) and battery storage devices (e.g., electric vehicles).
Presence of active devices have made distribution grids dynamic, reconfigurable and an important location for smart
grid operations like demand response, frequency regulation
S. Park is a Ph.D student at Electrical Engineering Department of KAIST,
Daejeon, 34141, Korea. Email: sejun.park@kaist.ac.kr
D. Deka is with the Theory Division and the Center for Nonlinear Systems
of LANL, Los Alamos, NM 87544. Email: deepjyoti@lanl.gov
S. Backhaus is with the A Division of LANL, Los Alamos, NM 87544.
Email: backhaus@lanl.gov
M. Chertkov is with the Theory Division and the Center for Nonlinear
Systems of LANL, Los Alamos, NM 87544 and with Center for Energy
Systems of Skoltech, Moskow, 143026, Russia. Email: chertkov@lanl.gov

and inter-household energy settlements/transactions. However,
optimal operations under different regimes require real-time
state estimation in the grid, in particular of the current radial
topology of current operational lines, and their impedances.
In addition, real or near real-time estimation of the distribution grid topology and corresponding line impedances is not
straightforward due to limited availability of real-time measurement devices, unlike in high voltage transmission grids.
In recent years, Phasor Measurement Unit (PMU) technology
and its alternatives (Eg. micro-PMUs [1], FNETs [2]) have
become available in distribution grids but their presence is
not ubiquitous [3]. Among others, presence of underground
lines in urban areas (e.g., New York City) makes meter
placement and parameter calibration challenging. Thus, there
is a greater need to develop low complexity algorithms that
can optimally estimate topology and line parameters that are
able to overcome the restrictions due to sparse meter presence
and infrequent calibration of line parameters.
In recent years, few nodal measurements devices (PMUs,
micro-PMUs, FNETs) have been placed in distribution grids.
More importantly, new loads like smart air-conditioners or
electric vehicles that connect to the grid at the end-user
level have the ability to measure and communicate nodal
quantities - voltages and injections. In this work, we consider
such scenarios and analyze the joint problem of topology
and parameter estimation in grid where measurements are
restricted to the end/terminal nodes of the smart distribution
grids.
A. Prior Work
Past research in topology or parameter estimation have
proposed different algorithms that differ primarily on the
availability of data and type of measurements (nodal or line
based). For available line measurements, [4] uses cycle basis
and maximum likelihood tests to estimate the topology. For
nodal voltage based measurements available at all nodes,
graphical model based formulations have been proposed to
estimate the operational lines for both radial grids [5], [6]
and loopy grids [7], [8]. In a similar measurement regime
including nodal voltages, [9], [10] present greedy topology
learning schemes based on trends in second moments of
voltage magnitudes. Real-data driven and model-free schemes,
that use signature based tests to reconstruct topology and line
parameters are presented in [11]–[13].

2

It is worth mentioning that majority of the prior work
rely on on availability of nodal measurements (voltage and/or
injection) at all nodes of the grid. In work involving missing
nodes [9], [10], injection statistics at all nodes are assumed to
be known. However this may be stronger assumption due to
unavailable meters and historical information. Further none
of the mentioned work provides guaranteed topology and
impedance estimation in the presence of missing nodes. In this
manuscript, we provide efficient algorithms for both topology
and impedance estimation on all operational lines in the grid in
a severely measurement deficient regime where voltage and/or
injection measurements at only end-users (terminal nodes in
radial grid) are assumed to be known while all other nodal
quantities are assumed to be unobserved.
B. Contribution
We present two related algorithms for topology and
impedance recovery with only end-user measurements in a
radial grid that are provably correct for linearized power flows
[9], [14]–[16]. The first algorithm takes time-stamped voltage
magnitude samples and complex injection statistics of endusers and identifies operational edges from among an overcomplete set of permissible edges with known impedances.
Operational edges and missing nodes are identified based on
novel relations between second order moments of voltage and
power injections at the observed leaf nodes. We show that the
algorithm has O(|V|3 ) computational complexity.
We improve upon this to present our second algorithm that
jointly estimates operational edges and their impedances using
time-stamped voltage magnitude and injection samples (not
just statistics) at the leaf nodes. Unlike the first algorithm,
the second algorithm does not require any knowledge on the
number of missing/unobserved nodes or the set of permissible
edges. In other words, no information other than the voltages
and injections at the terminal nodes are used. The second
algorithm first recovers the impedance distance between all
observed leaf nodes and then uses this to iteratively identify
each operational edge along with its impedance. The second
algorithm has a similar computational complexity as the first
algorithm. Further we prove that it has a O(|V| log |V|) sample
complexity. The step-by-step operation of our learning algorithms is demonstrated on toy examples. Simulations results
on IEEE test cases with ac power flow models demonstrate the
practical use of our algorithms. To the best of our knowledge,
this is the first work which provides guaranteed topology and
impedance reconstruction in the distribution grids with only
terminal node observability. Summary on both algorithms is
presented in Table I. Parts of the work have been presented
in IEEE SmartGridComm 2016 [17] and PSCC 2018 [18].
The rest of the manuscript is organized as follows. Section
II introduces nomenclature and power flow relations in the
distribution grid. The first algorithm for topology estimation is
described in Section III. Next, we present the second algorithm
for joint topology and impedance estimation in Section IV.
We also include a detailed derivation of its sample and
computational complexity. Numerical experiments on IEEE
test cases are presented in Section V. Finally, Section VI is
reserved for conclusions and discussion of future work.

II. D ISTRIBUTION G RID S TRUCTURE AND P OWER F LOWS

missing

Fig. 1: Illustration of a radial distribution grid. The large red
node denotes the substation bus. Leaf nodes that represent
households/end-users are observed and colored solid blue,
while unobserved and possibly unknown intermediate nodes
are denoted by dotted blue circles. The unknown operational
lines are colored grey.
Radial Structure: We consider radial distribution grids in this
paper. Mathematically we define the grid by graph G = (V, E),
where the set of buses/nodes is denoted by V and the set of
undirected operational lines/edges is denoted by E. We denote
the set of all edges/lines including non-operational ones by Ē.
The operational grid is ‘radial’ with the substation bus as a root
node. Fig. 1 shows an illustration of a distribution grid. We
use lower-case alphabets a, b, c, . . . to represent buses/nodes
and pairs (ab) to denote a line/edge between nodes a and b.
We denote t ∈ V as a root node (reference bus). We denote
Pab as the unique path from node a to node b in operational
radial grid G. For a node a, All nodes whose path to the root
contains a is called descendants of a and denoted by a set
Da . If (ab) is an edge b ∈ Da , then node b is called ‘child’
of ‘parent’ node b. Nodes that children of same parent are
called ‘siblings’. L ⊂ V denotes the set of leaf nodes that are
observed in our learning algorithms. The remaining nodes in
the grid (termed intermediate/missing nodes) are assumed to
be unobserved. Next we discuss the power flow models used
in this paper for algorithm design.
Power Flow Models: In a radial grid G = (V, E), the
following Kirchhoff’s law of ac power flow expresses complex
power injections at each node in terms of the node-voltages
and line-impedances:
X v 2 − va vb exp(iθa − iθb )
a
.
(1)
pa + iqa =
∗
zab
b:(ab)∈E

Here, zab , va , θa , pa , qa denote impedance of (ab) ∈ E, voltage
magnitude, voltage phase, active and reactive power at a ∈
V respectively. Since Eq. (1) is non-convex, we consider a
linearized approximation that neglects second order terms in
Eq. (1) termed Linear Coupled Power Flow (LC-PF) model
[9], [16]:
X 

pa =
βab (θa − θb ) + gab (va − vb )
b:(ab)∈E

qa =

X 
b:(ab)∈E


βab (va − vb ) − gab (θa − θb )

(2)

3

TABLE I: Summary of topology learning algorithms with missing nodes
Algorithm
1

Output
Topology

2

Topology
Impedances

Observations (available only at leaf nodes)
Time-stamped voltage magnitude samples
Complex power injection statistics
Time-stamped voltage magnitude samples
Time-stamped complex power injection samples

2
2
where gab = rab /(x2ab + rab
), βab = xab /(x2ab + rab
) and
rab , xab are resistance, reactance of line (ab) respectively, i.e.,
zab = rab +ixab . Following the standard notation, we consider
the substation/root node as a reference bus and measure
voltage magnitude and phase at each non-substation bus with
respect to it. Further, due to the lossless nature of linearized
power flow model, the injection at the reference bus is the
negative of the sum of injections at all other nodes. We can
thus ignore the reference bus from the power flow analysis and
consider a reduced model comprising of power flow equations
at the non-reference buses in the grid. Further, by considering
only deviations from the respective steady state values, we
model p, q, v, θ as random variables with zero mean. Note that
the steady state mean value of each observed random variable
can be computed and subtracted out.
Linear Coupled Power Flow Model: Note that LC-PF model
Eq. (2) is equivalent to a first order approximation of voltage
magnitudes in the LinDistFlow equations introduced in [14],
[15], [19] for distribution grids. The LC-PF model Eq. (2) can
also be stated in the following matrix form [9]
−1
−1
v = H1/r
p + H1/x
q

−1
−1
θ = H1/x
p − H1/r
q

(3)

where v, θ, p, q are respectively vectors of voltage magnitude,
voltage phase, active and reactive power at the non-substation
buses of the grid. H1/r , H1/x represents the reduced weight
Laplacian matrices for G \ {t} where 1/rab , 1/xab are used
as edge-weights of (ab) respectively.1 We remind the reader
that va , θa , pa , qa are normalized to have zero mean. We also
−1
−1
mention a structural property of H1/r
, H1/x
that arises due to
the radial topology.
Lemma 1 ( [9], [20]). Let H1/r be the reduce weighted
Laplacian matrix for grid G. Then its inverse satisfies
X
−1
H1/r
(a, b) =
rcd .
(4)
(cd)∈Pat

T

Pbt

−1
H1/r

Thus, the (a, b)th entry in
is equal to the sum of line
resistances on edges common to paths from node a and b to
the root. As Pat ⊂ Pbt for parent child pair a, b, Eq. (4) gives
the following for parent node a, child node b and ∀c,
(
rab if c is a descendant of b
−1
−1
H1/r (a, c) − H1/r (b, c) =
.
0
otherwise
(5)
A. Assumptions on Distribution Grid Structure
We now present two assumptions on distribution grid structure and statistics of power injections which is required for
our topology/impedance learning algorithms.
1G

\ {t} denotes a subgraph of G induced by V \ {t}.

Prior Information
Line impedances of all permissible lines
None

Assumptions
Uncorrelated power injections
Missing nodes have degree ≥ 3
Uncorrelated power injections
Missing nodes have degree ≥ 3

Assumption 1. All missing nodes have a degree at least 3.
Assumption 1 implies that each missing node has at least
two children and that all leaf nodes are observed. Note that
Assumption 1 ensures that the total number of missing nodes
is at least one less than the number of observed nodes. In
its absence, the system is thus under-determined and multiple
configurations satisfy the available measurements (see [17]).
We discuss later topology learning without Assumption 1
learns a Kron-reduced model [21] of the grid where nodes
of degree 2 are reduced. Assumption 1 is akin to similar
assumptions for recovery in graphical models [22]. In addition,
we assume that the complex power injections at different nodes
are uncorrelated
Assumption 2. E[pa pb ] = E[qa qb ] = E[pa qb ] = 0 ∀a 6= b.
As considered in prior studies [5], [9], Assumption 2 is welljustified over sufficiently short time intervals while considering
deviations of injections at end-users. For intermediate nodes
that are involved in separation of power into downstream lines,
leakage or device losses contribute to the net power injection
and are independent from other nodes. In particular, note that
Assumption 2 does not restrict the class of distributions that
can be used to model individual nodes power injection and
applies for both positive and negative nodal injections. We
discuss techniques to extend our work to cases with correlated
user injection profiles and multi-phase systems in future work.
III. T OPOLOGY L EARNING A LGORITHM WITH VOLTAGE
S AMPLES
In this section, we discuss properties of voltages and injection at leaf nodes and use them to design the first topology
learning algorithm, Algorithm 1, introduced in [17]. The
algorithm uses voltage samples and injection statistics at all
leaf nodes and identifies all operational edges from a overcomplete set of permissible edges Ē with known impedances.
Using the LC-PF Eqs. (3), we can write the second moments
of nodal voltages with that of nodal injections as:
−1
−1
−1
−1
E[vv T ] = H1/r
E[ppT ]H1/r
+ H1/x
E[qq T ]H1/x
−1
−1
−1
−1
+ H1/r
E[pq T ]H1/x
+ H1/x
E[qpT ]H1/r

(6)

Note that E[ppT ], E[qq T ], E[qpT ], E[pq T ] are diagonal matrices from Assumption 2. For the notational convenience,
we first define the variance φab of the difference of voltage
measurements of a, b ∈ V as follow:
φab = E[(va − vb )2 ] = E[va2 ] + E[vb2 ] − 2E[va vb ].

(7)

Note that φab can be computed for all leaf pairs in L using the
observed voltage samples. Under LC-PF model, the following
result holds:

4

Theorem 1. Let a, b ∈ L have a common parent k1 ∈ V \ L.
Then

1:

2
E[p2a ] + x2ak1 E[qa2 ] + 2rak1 xak1 E[pa qa ]
φab = rak
1

+ rk21 b E[p2b ] + x2k1 b E[qb2 ] + 2rk1 b xk1 b E[pb qb ].

(8)

The derivation follows by expanding φab using Eq. 6 and
using Lemma 1. Note that aside from pathological cases,
Theorem 1 is satisfied only by true parent k1 of nodes a, b.
Thus, the equality can be used to identify the true parent of
sibling leaves (see Fig. 2a, 2b for an example). The next result
involves φ values at three leaf nodes in L.
Theorem 2 (Theorem 2, [17]). Let leaf nodes a, b ∈ L have
common parent node k1 . Consider leaf c ∈ L such that c, k1 ∈
Dk2 (descendants of node k2 ) and Pk1 t ∩Pct = Pk2 t for some
intermediate node k2 (see Fig. 2f for an example). Then
φac −φbc = E[p2a ]((rak2 )2 − (rkk12 )2 ) + E[qa2 ]((xka2 )2 − (xkk21 )2 )
+ 2E[pa qa ](rak2 xka2 − rkk12 xkk21 ) − E[p2b ]((rbk2 )2 − (rkk12 )2 )
− E[qb2 ]((xkb 2 )2 − (xkk21 )2 ) − 2E[pb qb ](rbk2 xkb 2 − rkk12 xkk21 )
(9)
P
P
e
e
where rd = (f g)∈Pde rf g , xd = (f g)∈Pde xf g .
The proof for Theorem 2 uses algebraic expansions of the
expression for φac − φbc and application of Lemma 1. we
refer the reader to check [17] for details. Crucially, Theorem
2 enables us to identify edges between missing intermediate
nodes. For example, consider the case in Fig. 2c. If edges from
node k1 to leaves a, b have already been discovered, then we
can assert the existence of edge (k1 k2 ) to k1 ’s parent k2 by
checking if Eq. 9 holds. We use this result in our algorithm
to learn edges iteratively from parents of leaves to the root.
However it needs to be mentioned that the right side of Eq. 9
does not depend on the path to node c. Thus c can only be
identified as a descendant of k2 . Its true location cannot be
identified in particular if leaf node c do not have another leaf
node as sibling. The location of such leaf nodes are determined
once the rest of the network is recovered. We arrange the
identified intermediate edges in reverse order and check for
Eq. (9) to arrive at the true parent of unidentified leaf node c.
The post-order node traversal [23] is necessary to ensure that
the true parent of c is checked before other intermediate nodes
on Pct (path to the root). The steps are outlined in Algorithm
1.
Now, we briefly explain how Algorithm 1 works. Voltage
samples and injection statistics at leaf nodes are taken as input.
The output is the set of operational edges E from among the
set of input permissible edges Ē with known impedances. In
Steps 5-11, Algorithm 1 identifies sibling relationships of leaf
nodes and find their parent using Theorem 1. In Steps 12-24,
the algorithm identifies the edge between missing nodes using
Theorem 2, as explained in the example in the last paragraph.
Similarly, in Steps 26-33, the algorithm finds the parent of leaf
nodes without sibling leaf nodes. Fig. 2 illustrates Algorithm
1 on a test-case step by step.
2 Eq.

(9) is checked under the assumption that (kk2 ) exists.

Algorithm 1 Topology Learning Algorithm with Voltage
Magnitude Samples

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:

Input: L, M = V \ L, {E[p2a ], E[pa qa ], E[qa2 ] : a ∈ L},
v 1 , . . . , v m , {rab , xab : (ab) ∈ Ē}
Output: (V, E)
Initialization: para ← ∅, desa ← ∅ for all a ∈ V, E ← ∅
Compute φab = E[(va − vb )2 ] for all a, b ∈ L
for a ∈ L do
if para = ∅, ∃b ∈ L, ∃k1 ∈ M s.t. a, b, k1 satisfy
Eq. (8) with tolerance τ1 then
E ← E ∪ {(ak1 ), (bk1 )}
para ← {k1 }, parb ← {k1 }, desk1 ← {a, b}
end if
end for
L ← {a : a ∈ L, para = ∅}
do
M1 ← {k : k ∈ M, park = ∅, desk 6= ∅}
M2 ← {k : k ∈ M, desk = ∅}
for k ∈ M1 with a, b ∈ desk do
k1 ← para
if ∃k2 ∈ M − M2 , φac − φbc satisfies Eq. (9) with
tolerance τ2 for some c ∈ desk2 then2
E ← E ∪ {(kk2 )}, park ← {k2 }
else if ∃k2 ∈ M2 s.t. φac − φbc satisfies Eq. (9)
with tolerance τ2 for some c ∈ L then
E ← E ∪ {(kk2 )}, park ← {k2 }
desk2 ← desk
end if
end for
while |{k1 : k1 ∈ M1 , park1 6= ∅}| > 0
Form a post-order traversal node set W using para for all
a ∈ M such that desa 6= ∅
for c ∈ L do
for j = 1 to |W| do
k2 ← W(j) with a, b ∈ desk2 , k1 ← para
if φac −φbc satisfies Eq. (9) with tolerance τ2 then
E ← E ∪ {(ck2 )}, j ← |W|, W ← W \ {k2 }
end if
end for
end for
if |M1 | = 1 then
Join k ∈ M1 to root
end if

Computational Complexity of Algorithm 1: Algorithm 1
has three major parts, Steps 5-11, Steps 12-24 and Steps
26-33 where the rest part has computational complexity
O(|V|2 ) which arises from Step 4. Steps 5-11 iterate over
a set (a, b, k1 ) ∈ L × L × M where each iteration takes
O(1) computations. Therefore, the complexity for steps 511 is O(|V|3 |). Steps 12-24 and steps 26-33 iterate over sets
(c, k1 , k2 ) ∈ L × M1 × M and (c, k2 ) ∈ L × M respectively,
where each iteration takes O(1) computations. Therefore,
the complexity for Steps 12-24 and 29-24 are O(|V|3 |) and
O(|V|2 |) respectively. Hence, the overall computational complexity of Algorithm 1 is O(|V|3 ).

5

𝑘2
𝑘2

𝑘1
𝑎

(a)

(b)

𝑐
𝑎

𝑏

𝑘2

𝑘

𝑘1 = 𝑘
𝑏

(c)

𝑘1
𝑎

𝑘1

𝑐
𝑏

(d)

𝑘2

𝑎

(e)

𝑘1

𝑐
𝑏

𝑎

𝑏

(f)

Fig. 2: Illustration of Algorithm 1 (a) an original topology with missing blue nodes and observed green nodes (b & c) finding
parent k1 of sibling leaf nodes a, c (Steps 5 − 11) (d & e) iterative recovery of missing parent and grandparent of k1 (Steps
12 − 24) (f) finding parents of leaf nodes without sibling leaves (Steps 25 − 36) to get recovered topology

Modification for Finite Samples: Note that in reality due to
finite samples, the equality relations (8, 9) will not hold with
equality. In that setting, we compute the relative difference
between the left and right sides for either relation (8) or (9).
We consider the relations to be satisfied in Algorithm 1 if
the relative differences are respectively less than user defined
tolerances τ1 , τ2 .
The main bottleneck of Algorithm 1 is that it requires
knowledge about permissible lines and impedances which
might not be available in real distribution grids. In particular,
the line impedances of operational lines may need to be
estimated. To overcome these drawbacks, we propose a new
algorithm in the next section which additionally requires
power injection samples on leaf nodes instead of covariance
measurements but is able to recover the true topology and line
impedances on operational lines. More importantly, it does not
require any information on missing nodes or permissible lines,
i.e., the algorithm only requires voltage and power samples of
observed leaf nodes.
IV. T OPOLOGY AND I MPEDANCE L EARNING A LGORITHM
WITH VOLTAGE M AGNITUDE AND P OWER S AMPLES
Our algorithm, termed Algorithm 3 uses time-stamped observations of voltage magnitudes and complex injections at the
end-nodes as input. Our algorithm uses the notion of additive
‘distance’ defined as a distance over the graph
P which satisfies
the weighted metric property, d(a, b) =
(cd)∈Pab d(c, d).
We first estimate the distances between all leaf node pairs,
and then utilize the recursive grouping algorithm [24] to learn
operational topology of the grid. Under Assumption 2 and
using the LC-PF Eq. (3) in the case of observed nodes a, b,
one derives the following identity
−1
−1
E[va pb ] = H1/r
(a, b)E[p2b ] + H1/x
(a, b)E[pb qb ]
−1
−1
E[va qb ] = H1/r
(a, b)E[pb qb ] + H1/x
(a, b)E[qb2 ]

(10)

where E[va pb ], E[va qb ], E[p2b ], E[pb qb ], E[qb2 ] can be computed
from measurements at observed nodes a and b. Using Eq. (10),
−1
−1
we can estimate the value of H1/r
(a, b) and H1/x
(a, b) for
2
2
any observed a, b ∈ V unless E[pb ]E[qb ] = (E[pb qb ])2 . To
avoid such pathological situations, we make the following
assumption.
Assumption 3. There exists a constant M such that for all
node a ∈ V, E[p2a ]E[qa2 ] − (E[pa qa ])2 ≥ M.

−1
Once H1/r
(a, b) is estimated, we can derive the resistance
distance (effective resistance) between observed nodes a, b as
X
−1
−1
−1
dr (a, b) =
rcd = H1/r
(a, a) + H1/r
(b, b) − 2H1/r
(a, b)
(cd)∈Pab

(11)
Note that for tree networks, effective resistance is an additive
distance metric between nodes a and b in the grid. Similarly,
one can estimate the additive reactance distance dx (a, b)
between observed nodes a, b. Following estimation of dr (a, b)
for all pairs of observed nodes, we can utilize the recursive
grouping algorithm (RG) [24] which, under Assumption 1,
leads to consistent topology and impedance estimation of the
power grid G.
A. Recursive Grouping Algorithm
Here we introduce the recursive grouping (RG) algorithm
that recovers the true radial topology given any additive
distance d(·, ·) between all leaf nodes (distances to non-leaf
nodes observed can be incorporated if available). Let exact
values of d(·, ·) be known for every pair of observed nodes.
RG algorithm uses the following lemma [24]. We note that
the notation ‘parent’ and ‘child’ in algorithms and lemmas in
this section is not related with the substation node as defined
in Section II.
Lemma 2. For Φabc := d(a, c)−d(b, c), the following relation
holds:
a) Φabc = d(a, b) for all c ∈ V \ {a, b} if and only if a is a
leaf node and b is its parent.
b) −d(a, b) ≤ Φabc = Φabc0 ≤ d(a, b) for all c, c0 ∈ V \
{a, b} if and only if a, b are leaf nodes with common
parent, i.e., they belong to the same group of siblings.
Using Lemma 2 a), the parent-child relationships for a set of
observed nodes O can be identified. Similarly, using Lemma
2 b), sibling groups can be identified.
The formal description of RG is given in Algorithm 2. The
input of RG is a set of observed nodes O ⊂ V and the additive
distance d(a, b) for all a, b ∈ O. We discuss the working
of RG steps through an illustrative example given in Fig. 3,
where green nodes represent O. First, RG finds groups of
siblings and their parents using Lemma 2, as shown in Fig.
3b. Edges are added between all identified parent-child pairs.
For identified siblings without observed parent, a parent node
is inserted and connected to its children as shown in Fig. 3c.

6

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 3: Illustration of Algorithm 2 (a) an original topology where blue nodes are missing and green nodes are observed (O)
(b) partitions Πs (marked by the dashed line) of O generated by a node, its siblings and parent in the first iteration of RG
(c) addition of edges and missing nodes and updated O after the first iteration of RG and partitions in the second iteration
(d) updated O after the second iteration and partition of O in the third iteration (e) result after the third iteration of RG (f)
recovered topology
d(·, ·) are updated for the newly added parents using the fact
that distances are additive. For siblings a, b ∈ O and their
newly added parent h, the distance d(a, h) and d(c, h) for any
c ∈ O are calculated by
1
(d(a, b) + Φabc ), any c ∈ O
2
d(c, h) = d(a, c) − d(a, h)

d(a, h) =

(12)
(13)

Finally, RG updates the set O with newly added parents and
nodes without established parent or child relations illustrated
by green nodes in Fig. 3c. The process is iterated and new
edges are added unless |O| ≤ 2, which applies when an edge
can be added to remaining vertices or when a single vertex is
left. Fig. 3d-3f illustrates iterations of the RG after the first
one.
For topology estimation in radial grids, we propose the
following two stage algorithm with missing modes:
1. For all pair of observed nodes a, b ∈ O, calculate
dr (a, b) and dx (a, b) using Eqs. (10, 11) and second order
moments.
2. Recover missing nodes and lines using the recursive
grouping algorithm.
The formal statement of the algorithm is presented in Algorithm 3. Note that by learning the impedances based distances,
the resistance and reactance of each operational edge is jointly
estimated along with the topology. This is possible due to the
availability of injection samples that enable computation of
the right side of Eqs. (10). The previous Algorithm 1 used
only injection statistics which are not sufficient for impedance
estimation. Next, we briefly discuss how RG works in the
absence of Assumption 1, i.e. when missing intermediate
nodes have degree 2 or some leaf node is unobserved.
Recursive Grouping without Assumption 1: Note that RG
estimates additive distance (here resistive/reactive distance)
to identify sibling nodes and then recovers their parent. If
some internal node a has degree 2 (its child has no sibling),
then a cannot be identified using Lemma 2 b). Instead if a’s
parent k1 has degree > 2, then a’s child will get connected to
k1 . In other words, RG outputs a topology without degree
2 nodes by adding edges between their parent and child.
This reduced graph is exactly the Kron-reduced model [21]
is a coarsest partition if for any Π0 and for any S 0 ∈ Π0 , there exists
S ∈ Π such that S 0 ⊂ S. The coarsest partition Π in Algorithm 2 represents
a collection of sets of siblings and their parent.
3Π

Algorithm 2 Recursive Grouping Algorithm (RG)
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:

Input: O, {d(a, b) : a, b ∈ O}
Output: (V, E), {d(a, b) : a, b ∈ V}
Initialization: V = O, E = ∅
while |O| > 2 do
ON EW ← ∅.
Compute Φabc = d(a, c) − d(b, c) for all a, b, c ∈ O.
Find a coarsest partition Π of O such that any two
nodes in S ∈ Π are either leaves and sibling, or a parent
and a leaf child.3
for S ∈ Π do
if |S| = 1 then
ON EW ← ON EW ∪ S.
else if a parent
 pS ∈ S exists then
E ← E ∪ (pS a) : a ∈ S \ {pS }
ON EW ← ON EW ∪ {pS }
else
Add a parent hS of S as follows
V ← V ∪ {hS }
E ← E ∪ (hS a) : a ∈ S \ {hS }
ON EW ← ON EW ∪ {hS }
end if
end for
Update d(·, ·) for ON EW using Eqs. (12, 13).
O ← ON EW .
end while
if |O| = 2 then
E ← E ∪ {(ab) : a, b ∈ O, a 6= b}
end if

derived by removing degree 2 nodes from the grid graph.
Crucially, the line impedances of new edges will be the sums
of impedances of the two lines connected to corresponding
unidentified nodes. This also preserves the impedance in the
original graph between the parent and child of each degree
2 node. Similarly, if some leaf nodes are not observed, their
resistive/reactive distances to other leaves are not computed
and hence RG outputs the Kron-reduced model of the topology
without them.
B. Recursive Grouping with Finite Samples
In a practical scenario, due to finite number of injection
and voltage samples we can compute only the approximated

7

Algorithm 3 Topology/Impedance Learning Algorithm with
Voltage and Power Samples
Input: O, {E[va pb ], E[va qb ], E[p2a ], E[qa2 ], E[pa qa ] : a, b ∈
O}
2: Output: (V, E), {rab , xab : (ab) ∈ E}
3: for "
a, b ∈ O do#

−1 

−1
H1/r
(a, b)
E[p2b ] E[pb qb ]
E[va pb ]
4:
←
−1
E[pb qb ] E[qb2 ]
E[va qb ]
H1/x
(a, b)
1:

5:
6:
7:
8:
9:
10:
11:
12:
13:

end for
for a, b ∈ O do
−1
−1
−1
dr (a, b) ← H1/r
(a, a) + H1/r
(b, b) − 2H1/r
(a, b)
−1
−1
−1
dx (a, b) ← H1/x (a, a) + H1/x (b, b) − 2H1/x (a, b)
end for
(V, E), {dr (a, b) : a, b ∈ V} ← RG(O, {dr (a, b) : a, b ∈
O})
for (ab) ∈ E do
rab ← dr (a, b), xab ← dx (a, b) where dx (a, b) is
obtained using (V, E)
end for

value dbr of dr rather than the exact value. In other words
the variance of the distance is nonzero. To account for it, we
allow some tolerance ε for finding parent-child and sibling
relationships in Lemma 2. In addition, to test the relationship
of a, b we only use nodes that are close enough to both a and
b, i.e., nodes in Kab where Kab satisfies
Kab = {c ∈ O \ {a, b} : dbr (a, c), dbr (b, c) < τ }
for some constant τ . Let us now present rules which guide the
relationships of nodes using samples.
b abc | ≤ ε for all
a) Set a as a parent of b if |dbr (a, b) − Φ
c ∈ Kab .
b) Set a, b as siblings if
b abc − min Φ
b abc ≤ ε.
max Φ
c∈Kab

c∈Kab

Update of the distance is done in a similar manner where we
use averaging to mitigate the variability from finite sample
sizes. For a ∈ O and its newly added parent h, we set
dbr (a, h) =
1
2(|C(h)| − 1)

X
b∈C(h)\a

X
1
b abc
dbr (a, b) +
Φ
|Kab |

!

c∈Kab

where C(h) denotes the children set of h. Likewise, for c ∈
/
C(h),

X 
1
dbr (c, h) =
dbr (a, c) − dbr (a, h) .
|C(h)|
a∈C(h)

C. Sample and Computational Complexity
Algorithm 3 terminates in O(d|V|3 ) steps where d is the
depth of the grid (number of iterations) and |V|3 follows
from the calculation of Φabc . Here we analyze and show that
under some mild assumptions, the algorithm requires only
O(|V| log |V|) samples to correctly recover the operational

topology and line impedances, as formalized in the next
theorem.
Theorem 3. Suppose that a radial graph/grid (V, E) has a
constant depth. Under Assumptions 1-3 and assuming the LCPF model, if line impedances are bounded by nonzero value,
nodal power injections are sub-Gaussian with constantly
bounded sub-Gaussian parameters, and the number of samples
is greater than C|V| log(|V|/η) for some constant C, then
Algorithm 3 recovers the true topology and impedances with
probability 1 − η.
The definition of sub-Gaussian and the proof of Theorem 3
are presented in the supplementary material.

D. Recursive Grouping in non-linear Power flows
In the last section, we introduce the topology and impedance
learning algorithm, Algorithm 3. As discussed, Algorithm 3
provably outputs exact topology and impedance with low
sample/computational complexity in linearized LC-PF model.
However, as we are interested in simulations over samples
generated by non-linear ac power flow model, there are some
limitation for directly applying Algorithm 3 into real examples.
In this section, we address these bottlenecks and propose a
simple variant of Algorithm 3 for practical implementations.
First, Algorithm 2 is very sensitive to the tolerances used for
finite sample lengths. If tolerance is too small, the algorithm
outputs an error as it cannot find sibling relationships. In
contrast, if the tolerance is too large, the algorithm outputs
a loose topology with small number of missing nodes which
in the worst case can result in a star topology. Second, even
with infinite number of samples, as the real model is not
linear, the approximated distance does not converge to the
real distance. This causes serious problem in large grids as
different distances may have different errors that may not be
handled by a fixed tolerance. Third, the algorithm does not
utilize dr and dx at once.
To resolve these issues, we dynamically vary the tolerance
value in our experiments as follows:
1. In Algorithm 2, if no parent-child relationship is updated,
set d ← dx instead of dr and iterate. If parent-child
relationship is updated, set d ← dr and reset ε to the
input value.
2. If the algorithm does not find parent-child relationship
after step 1, increase the tolerance value ε ← αε (α > 1)
and set d ← dr . If parent-child relationship is updated,
reset ε to the input value. Otherwise, go to Step 1.
Note that this procedure uses both resistive and reactive distances to determine edges. If both fail due to a small tolerance
value, the algorithm dynamically increases the tolerance to
find the appropriate relationships. The possibility to consider
several values of tolerance can also help in handling the nonlinearity due ac power flow models. We note that this modified
algorithm is at least good as Algorithm 3 under the linearized
LC-PF model.

eps = 0.05
eps = 0.07
eps = 0.1

0.8
0.6
0.4
0.2

0.8
0.6
0.4
0.2
0

2000
5000 10000
number of samples

10%

eps = 0.05
eps = 0.07
eps = 0.1

8%
6%
4%
2%
1000

2000 5000 10000
number of samples

(a) 10 vertices

0.8
0.6

eps = 0.05
eps = 0.07
eps = 0.1

0.4
0.2

2000
5000 10000
number of samples

10%

eps = 0.05
eps = 0.07
eps = 0.1

8%
6%
4%
2%
0%

0%

1

0
1000

impedance error

1000

impedance error

eps = 0.05
eps = 0.07
eps = 0.1

1000

2000
5000 10000
number of samples

10%

impedance error

0

1

correct recovery ratio

1

correct recovery ratio

correct recovery ratio

8

eps = 0.05
eps = 0.07
eps = 0.1

8%
6%
4%
2%
0%

1000

2000 5000 10000
number of samples

(b) 50 vertices

2000
5000
10000
number of samples

(c) 100 vertices

Fig. 4: Experimental results for synthetic grids with 10, 50 and 100 vertices averaged (each case) over 100 random
radial grids.
P
|rab −r̂ab |
1
Upper graphs show average accuracy of our algorithm while lower graphs represent the impedance error 2|E|
(ab)∈E
|rab | +
|xab −x̂ab |
|xab |

observed within the correctly recovered topologies. “eps” in the graphs denotes ε.

V. E XPERIMENTS
In this section, we present experimental results of Algorithm
3 on custom examples and IEEE test cases for both LC-PF
(approximate) and non-linear ac power flow models.
Custom Examples with LC-PF samples: We first run simulations on randomly designed networks with voltages generated
by LC-PF model. For space constraints, we only simulate
Algorithm 3 for the random grids and postpone discussion
of Algorithm 1 to real grids.
In each simulation run we construct a random radial grid,
generate complex power injections and associate voltage measurements computed from the LC-PF model. For the topology
of the grid, we generate a random tree with maximum degree
5. The line resistance and reactance are independently sampled
from the uniform distribution over the interval [.1, .2]. For
complex power injections, we sample these from the independent normal distribution, i.e., pa , qa ∼ N (0, 1) and produce
nodal voltage magnitudes and phases using LC-PF Eq. (3).
The input of the algorithm is the complex power injections,
voltage magnitudes of the leaf/end-user nodes.
Under this setting, we run experiments changing number of
vertices from 10 to 100, the number of samples from 1000 to
10000, and also changing tolerance ε. To quantify performance
of our algorithm, we record accuracy in the recovered topology
and errors in estimated impedances, averaged over 100 random
radial grids for each instance of network size and number
of samples. Fig. 4 shows the correct recovery ratio and the
average error in estimating line P
impedances. The average error
|rab −r̂ab |
1
ab |
+ |xab|x−x̂
.
is defined for the grid as 2|E|
(ab)∈E
|rab |
ab |
One observes that our algorithm recover line impedances with
small error even in the demanding case of 1000 samples. We

also observe that larger ε results in a higher accuracy for the
small number of samples but it becomes less accurate for the
large number of samples (compare ε = 0.1 to ε = 0.07).
However, if the threshold is too small (ε = 0.05), the algorithm
performance decreases for all samples sizes. Note that similar
results (thus not shown) are derived when changing variance
of the complex power injections.

Fig. 5: Illustration of 56 bus distribution grid with 33 leaf/enduser nodes, 22 internal nodes and 1 substation (colored red).
IEEE test cases with non-linear ac samples: Here we discuss
more realistic simulations of Algorithm 1 and 3 on samples
generated by ac power flow model in test cases. We use a
matpower test case with 33 nodes [25], and a modified case
with 56 nodes [5] derived from the IEEE 123 test feeder [26].
Note that modifications were made to ensure radial topology
and that the internal nodes all have minimum degree ≥ 3 (see
Assumption 1). The modified grids are illustrated in Fig. 5

9

Fig. 6: Illustration of 33 bus distribution grid with 20 leaf/enduser nodes, 12 internal nodes and 1 substation (colored red).

algorithm works similarly for both MATPOWER samples and
the LC-PF samples. In Fig. 9, in accordance with observations
for custom model experiments, the algorithm performance
decreases as the threshold ε increases. For 33 bus model, we
perform similar experiments and report results in Fig. 10. One
can observe that the algorithm performs better for the 33 bus
network at low number of samples compare to the 56 bus test
case. This is because the depth of the grid is smaller and the
number of vertices is lesser.

and Fig. 6. We generate the complex power injections from
the independent normal distribution as in the case of the
custom models. However, we obtain the corresponding voltage
magnitude by using ac power flow equations in MATPOWER
[27]. We also compare the performance of our algorithm
on LC-PF samples generated with the same complex power
injections to see the effect of non-linear power flows.
We first show simulation results of Algorithm 1 where
the input includes complex nodal injection statistics, voltage
magnitude at all leaf/end-user nodes, and a set of permissible
edges with known impedances from which the operational
edges are determined. We consider the 33 bus network with
nodal injection variance of 10−4 p.u. per node. We include
50 additional edges of comparable impedances along with the
true operational edges to create the input permissible edge set
of 82 edges. Fig. 7 shows the errors in topology estimation
for different values of tolerances with increasing number of
samples. Note that for LC-PF and ac power flow voltage
magnitude samples, the errors are comparable.

number of errors

6

var = 1e-4, eps = 1e-4, MP
var = 1e-4, eps = 1e-4, LC
var = 1e-5, eps = 1e-4, MP
var = 1e-5, eps = 1e-4, LC
var = 1e-6, eps = 1e-4, MP
var = 1e-6, eps = 1e-4, LC

4

2

0
10000

6
4
2

number of errors

number of errors

8

tau1 = 8e-3, tau2 = 2e-2, LC
tau1 = 9e-3, tau2 = 2e-2, LC
tau1 = 9e-3, tau2 = 4e-2, LC
tau1 = 2e-2, tau2 = 6e-2, MP

50000

Fig. 8: Errors in topology estimation of Algorithm 3 observed
as the variance of the complex power injection changes in
56 bus model. “var” denotes the variance of the complex
power injections and MP, LC implies that the samples are
from MATPOWER and LC-PF respectively.

10
10

20000 30000 40000
number of samples

var = 1e-4, eps = 5e-5, MP
var = 1e-4, eps = 1e-4, MP
var = 1e-4, eps = 5e-4, MP
var = 1e-4, eps = 1e-3, MP

8
6
4
2
0
10000

0
40000 50000 60000 70000 80000 90000
number of samples

Fig. 7: Errors in topology estimation of Algorithm 1 for
different values of tolerances τ1 ,τ2 . MP, LC implies that the
samples are from MATPOWER and LC-PF respectively.
Next we discuss Algorithm 3 where the input comprises of
voltage and injection samples. Under this setting, we measure
performance of our algorithm by varying the number of
voltage and injection samples available, the variance of the
complex power injection and the threshold value, ε used in
Algorithm 3. To quantify errors in topology estimation, we
count the number of edge difference between the recovered
topology and the true topology. Fig.8 and Fig. 9 show our 56
bus model experimental results. In Fig. 8, we observe that the

20000 30000 40000
number of samples

50000

Fig. 9: Errors in topology estimation of Algorithm 3 by varying
ε in 56 bus model.

VI. C ONCLUSION AND F UTURE W ORK
Real-time operation and control in the distribution grid
relies on accurate and fast estimation of its topology and line
impedances from sparse data. In this manuscript, we present
two algorithms that recovery topology (and line impedances)
using voltage and injection measurements collected only from
the end-users/leaf nodes in the radial distribution grid, while
all intermediate nodes are unobserved. The first algorithm
uses injection statistics at end users and learns the topology.

10

number of errors

10

var = 1e-4, eps = 5e-5, MP
var = 1e-4, eps = 1e-4, MP
var = 1e-4, eps = 5e-4, MP
var = 1e-4, eps = 1e-3, MP

8
6

4
2
0

1000

2000
3000
4000
number of samples

5000

Fig. 10: Error in topology estimation of Algorithm 3 by
varying ε in 33 bus model.

On the other hand, the second algorithm owing to the presence of injection samples is able to do joint topology and
impedance estimation. We show that either algorithm has a
computational complexity which scales as |V|3 . Further we
show that under some mild technical conditions, the second
algorithm guarantees to output the correct topology with only
O(|V| log |V|) samples. The theoretical analysis and exactness
of both algorithms is discusses using a linearized ac power
flow model and related properties of radial networks. However,
we demonstrate the performance of our algorithms through
numerical simulations with samples generated from non-linear
ac power flow model in MATPOWER.
This work opens up several directions for possible extension. We plan to analyze such learning algorithms for
three phase power distribution grids under a similar linearized
scheme. Further extension to the case with correlated injections at the end-users will be analyzed to make the algorithms
function under more generalized conditions. Finally we plan
to pursue theoretical extensions of this work from radial grids
to the case of loopy grids with large girth, where effective
distances are approximately additive.
R EFERENCES
[1] A. von Meier, D. Culler, A. McEachern, and R. Arghandeh, “Microsynchrophasors for distribution systems,” Innovative Smart Grid Technologies Conference (ISGT), 2014 IEEE PES, pp. 1–5, 2014.
[2] Z. Zhong, C. Xu, B. J. Billian, L. Zhang, S.-J. S. Tsai, R. W.
Conners, V. A. Centeno, A. G. Phadke, and Y. Liu, “Power system
frequency monitoring network (fnet) implementation,” Power Systems,
IEEE Transactions on, vol. 20, no. 4, pp. 1914–1921, 2005.
[3] R. Hoffman, “Practical state estimation for electric distribution networks,” in IEEE PES Power Systems Conference and Exposition. IEEE,
2006, pp. 510–517.
[4] R. Sevlian and R. Rajagopal, “Feeder topology identification,” arXiv
preprint arXiv:1503.07224, 2015.
[5] S. Bolognani, N. Bof, D. Michelotti, R. Muraro, and L. Schenato, “Identification of power distribution network topology via voltage correlation
analysis,” in Decision and Control (CDC), 2013 IEEE 52nd Annual
Conference on. IEEE, 2013, pp. 1659–1664.
[6] D. Deka, S. Backhaus, and M. Chertkov, “Estimating distribution grid
topologies: A graphical learning based approach,” in Power Systems
Computation Conference (PSCC), 2016. IEEE, 2016, pp. 1–7.
[7] Y. Liao, Y. Weng, G. Liu, and R. Rajagopal, “Urban distribution grid
topology estimation via group lasso,” arXiv preprint arXiv:1611.01845,
2016.

[8] D. Deka, M.Chertkov, S. Talukdar, and M. V. Salapaka, “Topology
estimation in bulk power grids: Theoretical guarantees and limits,” in
accepted in the Bulk Power Systems Dynamics and Control SymposiumIREP, 2017.
[9] D. Deka, M. Chertkov, and S. Backhaus, “Structure learning in power
distribution networks,” IEEE Transactions on Control of Network Systems, 2017.
[10] D. Deka, S. Backhaus, and M. Chertkov, “Learning topology of the
power distribution grid with and without missing data,” in Control
Conference (ECC), 2016 European. IEEE, 2016, pp. 313–320.
[11] G. Cavraro, R. Arghandeh, A. von Meier, and K. Poolla, “Data-driven
approach for distribution network topology detection,” arXiv preprint
arXiv:1504.00724, 2015.
[12] V. Arya, T. Jayram, S. Pal, and S. Kalyanaraman, “Inferring connectivity
model from meter measurements in distribution networks,” in Proceedings of the fourth international conference on Future energy systems.
ACM, 2013, pp. 173–182.
[13] J. Peppanen, J. Grimaldo, M. J. Reno, S. Grijalva, and R. G. Harley, “Increasing distribution system model accuracy with extensive deployment
of smart meters,” in PES General Meeting— Conference & Exposition,
2014 IEEE. IEEE, 2014, pp. 1–5.
[14] M. Baran and F. Wu, “Optimal sizing of capacitors placed on a radial
distribution system,” Power Delivery, IEEE Transactions on, vol. 4,
no. 1, pp. 735–743, Jan 1989.
[15] ——, “Optimal capacitor placement on radial distribution systems,”
Power Delivery, IEEE Transactions on, vol. 4, no. 1, pp. 725–734, Jan
1989.
[16] S. Bolognani and S. Zampieri, “On the existence and linear approximation of the power flow solution in power distribution networks,” Power
Systems, IEEE Transactions on, vol. 31, no. 1, pp. 163–172, 2016.
[17] D. Deka, S. Backhaus, and M. Chertkov, “Learning topology of
distribution grids using only terminal node measurements,” in IEEE
Smartgridcomm, 2016.
[18] S. Park, D. Deka, and M. Chertkov, “Exact topology and parameter
estimation in distribution grids with minimal observability,” in Power
Systems Computation Conference (PSCC), 2018 (accepted to appear).
IEEE, 2018.
[19] M. E. Baran and F. F. Wu, “Network reconfiguration in distribution
systems for loss reduction and load balancing,” IEEE Transactions on
Power delivery, vol. 4, no. 2, pp. 1401–1407, 1989.
[20] J. Resh, “The inverse of a nonsingular submatrix of an incident matrix,”
IEEE Transactions on Circuit Theory, vol. 10, pp. 131–132, 1963.
[21] F. Dorfler and F. Bullo, “Kron reduction of graphs with applications to
electrical networks,” IEEE Trans. on Circuits and Systems I: Regular
Papers, vol. 60, no. 1, pp. 150–163, 2013.
[22] J. Pearl, Probabilistic reasoning in intelligent systems: networks of
plausible inference. Morgan Kaufmann, 2014.
[23] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, Introduction
to Algorithms. The MIT Press, 2001.
[24] M. J. Choi, V. Y. Tan, A. Anandkumar, and A. S. Willsky, “Learning latent tree graphical models,” The Journal of Machine Learning Research,
vol. 12, pp. 1771–1812, 2011.
[25] “IEEE 1547 Standard for Interconnecting Distributed Resources with
Electric Power Systems.” [Online]. Available: http://grouper.ieee.org/
groups/scc21/1547/1547 index.html
[26] W. H. Kersting, “Radial distribution test feeders,” in Power Engineering
Society Winter Meeting, 2001. IEEE, vol. 2. IEEE, 2001, pp. 908–912.
[27] R. D. Zimmerman, C. E. Murillo-Sánchez, and R. J. Thomas, “Matpower: Steady-state operations, planning, and analysis tools for power
systems research and education,” IEEE Transactions on power systems,
vol. 26, no. 1, pp. 12–19, 2011.
[28] R. Vershynin, “Introduction to the non-asymptotic analysis of random
matrices,” arXiv preprint arXiv:1011.3027, 2010.
[29] O. Rivasplata, “Subgaussian random variables: An expository note,”
2012.

11

A PPENDIX
Definition 1. A zero mean random variable X is sub-Gaussian
2
if there exists a constant A ≥ 0 such that P(|X| > t) ≤ e1−At
for all t ≥ 0.
A. Proof of Theorem 3
Proof. We first provide the following key lemma that the
estimated resistance distances are uniformly bounded from the
true distances if Ω(|V| log(|V|/η)) samples are given.
Lemma 3. For any constant δ, there exists a constant C such
that if the number of samples is greater than C|V| log(|V|/η),
then
b abc − Φabc | ≤ δ
|Φ
∀a, b, c ∈ V
with probability 1 − η.
b abc − Φabc | ≤ δ holds
One can easily observe that once |Φ
for sufficiently small constant δ for all a, b, c ∈ V, RG recovers
the true topology. Now we prove Lemma 3 and completes the
proof.
Proof of Lemma 3. Before starting the proof, we first introduce the following definition and lemmas about sub-Gaussian
and sub-exponential random variables.
Lemma 4 (Lemma 5.5 of [28], Theorem 2.7 of [29]). For subGaussian random variables X, Y , the following statements
hold.
a) X + Y is sub-Gaussian
2 2
b) There exists a constant σ such that E[eλX ] ≤ eλ σ /2 .
c) There exists a constant C such that E[X 2 ] ≤ C.
Definition 2. A random variable X is sub-exponential if there
exists a constant A such that P(|X| > t) ≤ e1−At .
Lemma 5 (Lemma 5.14, Lemma 5.15, Proposition 5.16 of
[28], Theorem 2.7 of [29]). The following statements hold:
a) A random variable X is sub-Gaussian if and only if
X 2 is sub-exponential. In addition, X 2 − E[X 2 ] is subexponential
b) For a zero-mean sub-exponential random variable X,
there exists constant σ, K ≥ 0 such that E[eλX ] ≤
2 2
eλ σ /2 for all |λ| ≤ K.
c) For zero-mean random variables X1P
, . . . , Xn satisfying
2 2
n
E[eλXi ] ≤ eλ σi /2 for λ ≤ Ki , X = i=1 αi Xi satisfies
that
2 Pn
2
E[eλX ] ≤ eλ ( i=1 αi σi )
for λ ≤ mini Ki .
d) Let X1 , . . . , Xn be independent zero-mean random vari2 2
ables satisfying E[eλXi ] ≤ eλ σ /2 for λ ≤ K, for some
σ, K ≥ 0 and for all i. Then, for all t ≥ 0, we have
!

 2

n
X
t
tK
P
Xi ≥ t ≤ 2 exp − min
,
.
4nσ 2 2
i=1
We note that the later statement of Lemma 5 a) directly
follows from the definition and Lemma 5 c) is an analogue of
Theorem 2.7 of [29]. To obtain the result of Lemma 3, we first
note that (pb + qc )2 − E[(pb + qc )2 ], p2b − E[p2b ], qb2 − E[qb2 ] are

zero-mean sub-exponential due to Lemma 4 a) and Lemma 5
a). Using this fact, one can show that
E[eλ(pb qc −E[pb qc ]) ]
h 1
i
2
2
2
2
2
2
= E e 2 λ(((pb +qc ) −E[(pb +qc ) ])−(pb −E[pb ])−(qc −E[qc ]))
≤ eλ

2

σ 02 /2

(14)
0

0

0

0

for some constants K , σ where K , σ are chosen to satisfy
same inequality for p2b , qb2 . The inequality of Eq. (14) follows
from Lemma 5 c).
2 2
Now we show that E[eλ(va pb −E[va pb ]) ] ≤ eλ σ /2 for |λ| ≤
1/2
K/|V|
for some constants σ, K ≥ 0.
E[eλ(va pb −E[va pb ]) ]
 P 

−1
λ
H −1 (a,c)(pb pc −E[pb pc ])+H1/x
(a,c)(pb qc −E[pb qc ])
= E e c∈V 1/r
(15)
≤e

λ2 σ 2 |V|/2

for λ ≤ K

(16)

where Eq. (15) follows from Eq. (2) and Eq. (16) follows from
−1
−1
Lemma 5 c), Eq. (14) and the fact that H1/r
(a, c), H1/x
(a, c)
are constantly bounded.
Using Eqs. (14, 16), we address to bound the difference beb a pb ]−E[va pb ]|
tween the sample mean and the expectation |E[v
b
where E denotes the sample mean. To this end, we define the
error event
b a pb ] − E[va pb ]| ≥ ε}.
Eva pb (ε) := {|E[v
Lemma 5 d) directly leads us to obtain the following inequality
for any constant ε:


 2
ε n Kεn
,
P(Eva pb (ε)) ≤ 2 exp − min
4|V|σ 2 2
η
(17)
≤
4|V|2
for n ≥ C|V| log(|V|/η) for some constant C where n denotes
the number of samples. One can observe that same inequality
holds for va qb Similarly, for small enough constant ε we have

 2

ε n K 0 εn
P(Ep2b (ε)) ≤ 2 exp − min
,
2σ 02
2
η
≤
(18)
6|V|
for n ≥ C 0 log(|V|/η) for some constant C 0 . Same inequality
holds for pb qb and qb2 .
Now, we define the global error event
(
[
P(E(ε)) :=
(Eva pb (ε) ∪ Eva qb (ε))∪
a,b∈V

)
[

(Ep2b (ε) ∪ Epb qb (ε) ∪ Eqb2 (ε)) .

b∈V

Using two inequalities Eqs. (17, 18), we apply the union bound
to obtain to bound the error probability as follows:
η
η
P (E(ε)) ≤ 2|V|2 ×
+ 3|V| ×
= η.
2
4|V|
6|V|
Once we obtain the above union bound, if E[va pb ], E[va qb ],

12

E[p2b ],E[pb qb ],E[qb2 ] are constantly bounded, one can choose a
small enough constant ε so that if E(ε) does not occur using
b abc −Φabc | ≤ δ for all a, b, c ∈ V. Since
Assumption 3, then |Φ
Lemma 4 c) shows that E[va pb ], E[va qb ], E[p2b ], E[pb qb ], E[qb2 ]
are bounded, this completes the proof of Lemma 3.

