# Ensembles of Realistic Power Distribution Networks

Rounak Meyur, Anil Vullikanti, Samarth Swarup, Henning Mortveit, Virgilio
Centeno,
Arun Phadke, H. Vincent Poor and Madhav Marathe

## Preliminaries of Power Distribution System

The power distribution network connects the HV (greater than 33kV)
distribution substation to LV (208-480V) residential consumers. In most
practical networks, this is done through two sets of networks: (i) the MV
(usually 6-11kV) primary network connects a step-down transformer at the
substation to local pole-top transformers along the road network and (ii) the
LV secondary network connects the residences to the local pole-top
transformers. Additionally, most distribution networks are operated in a
_radial_ structure (with no cycles or loops) to facilitate protection
coordination [1]. Residences are primarily connected in chains (the degree of
a node is at most 2) to avoid branching and thereby maintaining a healthy
voltage level. We term this configuration a _star-like tree_ since the edges
emerge from a single root transformer node and connect residences without
further branching. The schematic of a typical power distribution networks is
shown in Fig. 1.

Figure 1: Schematic of power distribution network (left) and corresponding
graph structure (right). The proposed approach first creates the secondary
network (red) connecting the residences to local transformers. Thereafter, the
road network (black edges) is used as a proxy to connect the transformers
forming the primary network.

## Related Work

Over the past decade, researchers have developed methods to create synthetic
power networks of different sizes. These attempts include interdependent
networks created using simple logical models [2] to models created based on
statistical attributes of actual networks [3].The simple logical network
models allowed network scientists to carry out a linear analysis which reveals
big picture results for large networks. However, neglecting rich system data
(line parameters and residence load demand) often makes these simple models
less useful in certain kinds of applications [4]. More recently, models with
actual system parameters have been proposed in [5, 6, 7, 8]. However, most of
these networks are limited to the transmission system and provides little
importance to the distribution power networks.

For a long time, power system researchers were dependent on the standard IEEE
test systems [9, 10] to carry out experiments and validate
control/optimization algorithms. These test systems were considered a replica
of the actual power grid and therefore, used as a test-bed before deploying a
methodology on the actual power system. These models either consider the
residential customers as passive consumers of electricity or assume an
aggregated version of multiple consumers as system load demand. These provided
an almost accurate system analysis as long as we do not consider the role of
consumers. However, with the advent of distributed energy resources (DERs) and
electric vehicles (EVs), residential customer behavior has become an integral
part of power system analysis [11, 12, 13]. Therefore, modern power network
models need to include detailed residential consumer attributes for accurate
representation of the actual power grid. Table 1 lists the different power
network models present in the literature and a summary of the network
features.

The literature dealing with the creation synthetic distribution networks can
be broadly classified into two categories [14, 15]. The first category
comprises methods that use one or more heuristics to synthesize the networks.
The second category of methods uses mathematical programming to encode various
physical and structural constraints that real-world networks usually obey. The
methods in the second category involve solving one or more optimization
problems that often require large computation time [16]. The work by
Schweitzer et al. [17] is one of the earliest papers in the second category.
The author analyze a large dataset of actual distribution networks in
Netherlands. Then they statistically fit these parameters to create a class of
stochastic network models. Random networks are then sampled from using the
stochastic model. This method is reminiscent of “configuration models” used to
construct random social networks [18]. The generated synthetic networks have
random interconnections, yet the network attributes follow similar
distributions as the actual networks. In [19], the authors have used an
extensive dataset of actual distribution networks and trained a generative
adversarial neural network (GANN). Thereafter, random networks are generated
from the trained generative model and hence multiple synthetic networks can be
generated which resemble actual distribution networks. However, these methods
are dependent on the availability of an extensive amount of actual
distribution network dataset. The lack of availability of such information
uniformly from all geographic regions makes the statistical models ‘biased’
and hence cannot be used as a generalized framework to create synthetic
networks.

On the other end of the spectrum, we have multiple heuristic methods [20, 21,
16, 22, 23, 24] to create random feeder networks with or without using the
road network information. The heuristics are primarily used in the initial
steps to assign substation locations and construct a feeder network. The
heuristics include clustering residential loads to assign substation location
at the cluster centroids [21], distribute an aggregated residential load to
certain load points in a hierarchical fashion [23] and use minimum spanning
tree algorithms [20, 21, 22]. Some of these methods construct an imaginary
road network [20, 21], while others have used openly available information
[22]. Neither of these methods include the power engineering constraints in
the network synthesis phase, rather the cables and shunt compensators are
chosen such that power engineering constraints are satisfied.

The mathematical programming methods are mostly restricted to distribution
system expansion planning [25, 26, 27, 28] or distribution system
restructuring problem [29, 30, 31]. We have followed similar works to
construct the mathematical framework for our approach (Steps 1 and 2).
However, one of the main difference is that in all these works the feeder
nodes are known beforehand. In our approach, the feeder nodes are identified
as an output of the proposed optimization framework. The development of such a
framework allows us to use it in Step 3 of our algorithm to construct an
ensemble of synthetic networks which are statistically similar to each other.

Table 1: Summary of power system network data available in literature Authors/Datasets | Description/Comments
---|---
EPRI Test Circuits [10] | A set of test feeders, representative of small, medium, and large-scale distribution networks.
IEEE Test Feeders [9] | A repository for small and large-scale distribution feeders with no associated geographic attribute.
GridLAB-D Feeders [32] | A set of 24 actual radial distribution feeders from 5 climate regions in the US.
D’Souza et.al [2] | Random graph models to create synthetic interdependent infrastructure networks.
Zussman et.al [3] | Synthetic transmission networks with geographic attributes created using a statistical distribution of actual network attributes.
Birchfield et.al [5, 6] | Large scale synthetic transmission networks for the state of Texas, US with synthetic loads aggregated at the zip-code level.
Overbye et.al [7] | Stochastic geometry-based approach was used to place transformers in the distribution network without physical constraints.
Trpovski et.al [27] | Optimization framework to create distribution networks connecting zipcode level aggregated consumer loads to substations.
Kadavil et.al [23] | Creates a feeder network and populates it with stochastic load profiles with no associated geographic attribute.
Schweitzer et.al [17, 22] | Distribution networks populated with random loads created to mimic statistical distribution of actual network attributes.
Mateo et al. [20, 21, 16] | Heuristic-based approach to connect consumers with no explicit structural and operational constraints.
Li et.al [8] |  Combined synthetic transmission and distribution networks created for Texas, US by adding models of [5, 6] with RNM.
Liang et.al [19] | Generative adversarial networks (GAN) based framework to create an ensemble of synthetic networks.
Bidel et.al [24] | Synthetic networks were created after learning the statistical distribution of load and network published by network operators.

Our work differs from prior works in the following ways. The points are
further summarized in Table 2.

  1. 1.

We propose an optimization framework which creates a minimum length network
which satisfies standard structural and power flow constraints. Earlier
papers, e.g. [24, 22, 16, 19] often do not use such well-defined mathematical
models that ensure all the constraints are satisfied while creating the
synthetic networks.

  2. 2.

Although the heuristics used in the methods proposed in [33, 22, 20, 21, 16,
34] are able to create a distribution network as a minimum spanning tree, they
do not ensure that the power engineering voltage constraints are satisfied as
per ANSI standards [35]. Therefore, the synthesized networks might be
topologically realistic but need not be be realistic from the perspective of
power engineering. The power engineering constraints are satisfied only after
adding shunt capacitors and voltage regulators, which improves the voltage
profile.

  3. 3.

Most heuristics described in previous papers [16, 34, 19, 24] deal with the
issue of finding locations of substation feeders by clustering residences. In
contrast, we start with given locations of substations [36] and residences
[37] based on real-world data and construct the network connecting these
points. Further, in our primary network creation method (Step-2), we have
included the optimal substation feeder selection problem effectively in our
optimization framework, which has not been used in any prior work.

  4. 4.

The power flow constraints used in our methodology are similar to the
constraints used in distribution feeder planning problems, as in [27, 28, 31,
38]. These constraints define the fundamental physics of distribution network
operation and hence are similar to related works. In addition, we include two
novel components in our optimization framework: (i) Prior works included
either a commodity flow model [31] to ensure tree structure, or included
multiple constraints to avoid the occurrence of cycles [29]. Here, we have
theoretically proved in Proposition (1) that the power flow constraints are
sufficient to ensure a tree structure. (ii) Prior papers on distribution
network planning and expansion, which include a similar optimization framework
assume that the number of feeders are known beforehand. In our framework, we
do not make such assumptions; rather, the optimal set of feeders are
identified as an output of the primary network creation problem.

  5. 5.

We construct an ensemble of networks which indeed is an important contribution
of the paper. We do not consider the optimal network as the ‘sole’ output;
rather we consider it as ‘one’ random realization of the actual distribution
network and propose a framework to create multiple feasible and realistic (but
not optimal) networks. Using our well-defined optimization framework in Step-1
and Step-2, we are able to create an ensemble of feasible networks by solving
a restricted version of the optimization problem. Prior works have not
considered this aspect of an ensemble of synthetic power distribution
networks.

Table 2: Table showing comparison with other related works Previous Works | Includes | Excludes
---|---|---
D’Souza et.al [2] | I | II,III,IV,V,VI,VII
Zussman et.al [3] | I | II,III,IV,V,VI,VII
Birchfield et.al [5, 6] | I | II,III,IV,V,VI,VII
Overbye et.al [7] | I | II,III,IV,V,VI,VII
Trpovski et.al [27] | I,IV,V | II,III,VI,VII
Kadavil et.al [23] | I | II,III,IV,V,VI,VII
Schweitzer et.al [17, 22] | I,II | III,IV,V,VI,VII
Reference Network Model (RNM) [20, 21, 16] | I,II,V | III,IV,VI,VII
Liang et.al [19] | I,VII | II,III,IV,V,VI
Bidel et.al [24] | I,II | III,IV,V,VI,VII
I | Geographic information embedding
II | High resolution residential hourly demand profile
III | Avoids usage of actual distribution networks
IV | Well defined mathematical framework with a guaranteed solution
V | Networks satisfies ANSI voltage constraints
VI | Optimal choice of feeders
VII | Ensemble of networks

## Methods

### Details of the dataset used

We use open-source, publicly available information regarding several
infrastructures to generate the synthetic distribution networks: (i) road
network data from OpenStreetMap (OSM) [39], (ii) geographic locations of HV
(greater than 33kV) substations from data sets published by EIA [36], and
(iii) residential electric power demand information developed in earlier work
from our research group [37]. The details of each dataset are provided in
Table 3. The rightmost column provides the size of each dataset for Montgomery
County in southwest Virginia, US.

Table 3: Datasets and related attributes used to generate synthetic distribution network Dataset | Source | Attributes |  | Example for
---
Montgomery County
of Virginia, US
Substation |  | Electric substation data
---
published by the US
Department of
Homeland Security [36]
|  • substation ID
---
• longitude
• latitude
20 substations
| Road
---
network
| GIS and electronic
---
navigable maps
published by
OpenStreetMap [39]
|  • node ID
---
• node longitude
• node latitude
• link ID
• link geometry
| 33882 nodes
---
41261 edges
Residences |  | Synthetic population
---
and electric load
demand profiles[37]
|  • residence ID
---
• longitude
• latitude
• hourly load
profile
35629 homes

### Step 1a: Map residences to the nearest road network link

Figure 2: Plots showing the steps in mapping a residence to the nearest road
network link. Bounding regions are drawn around each link (blue boxes) and the
residence (red box). The intersections (green boxes) between these bounding
regions are identified to shortlist nearby road links. The nearest road
network link (red edge) is then selected from these shortlisted links.

This section details the proposed mapping
$\mathscr{F}_{M}\colon\mathscr{H}\rightarrow\mathscr{E}_{R}$ between the set
of residences $\mathscr{H}$ and links $\mathscr{E}_{R}$ of road network
$\mathscr{G}_{R}\left(\mathscr{V}_{R},\mathscr{E}_{R}\right)$. We map each
residence $h\in\mathscr{H}$ to the nearest road network link
$e\in\mathscr{E}_{R}$. Thereafter, we compute the inverse mapping
$\mathscr{F}^{-1}\left(e\right)$ to identify a set of residential buildings
near each road network link $e\in\mathscr{E}_{R}$. This information would be
used in the successive steps to generate the secondary distribution network.
We denote the spatial embedding of residence $h\in\mathscr{H}$ as
$\mathbf{p}_{h}\in\mathbb{R}^{2}$. A road network link $e\in\mathscr{E}_{R}$
is denoted by $e:\left(u,v\right)$, where the nodes $u,v\in\mathscr{V}_{R}$
have spatial embedding $\mathbf{p}_{u},\mathbf{p}_{v}\in\mathbb{R}^{2}$
respectively.

Algorithm 1 is used to compute the nearest road network link to a given
residence. First, a bounding region of suitable size is evaluated for each
road network link. This is done such that any point in the region is within a
radius $r$ from any internal point of the road network link $e$. The spatial
embedding of an interior point along link $e=(u,v)$ is the convex addition of
the spatial embedding of the nodes $u$ and $v$. The bounding region for link
$e=(u,v)\in\mathscr{E}_{R}$ denoted by $\mathcal{B}_{e}$ is,

$\mathcal{B}_{e}=\big{\\{}\mathbf{p};\lvert\lvert\mathbf{p}-\mathbf{p}_{e}\rvert\rvert_{2}~{}\leq
r,\forall\mathbf{p}_{e}=\theta\mathbf{p}_{u}+(1-\theta)\mathbf{p}_{v},\theta\in[0,1]\big{\\}}$
(1)

Similarly, a bounding region is considered for a residential building
$h\in\mathscr{H}$ and denoted it by $\mathcal{B}_{h}$,

$\mathcal{B}_{h}=\big{\\{}\mathbf{p}\big{|}||\mathbf{p}-\mathbf{p}_{h}||_{2}\leq
r\big{\\}}$ (2)

The intersections between the bounding region of the building and those for
the links are stored and indexed in a _quad-tree_ data structure [40]. We
identify links $e_{1},e_{2},\cdots,e_{k}$ with bounding regions
$\mathcal{B}_{e_{1}},\mathcal{B}_{e_{2}},\cdots,\mathcal{B}_{e_{k}}$ which
intersect with $\mathcal{B}_{h}$. These $k$ links are comparably nearer to the
residential building than the others. Thus, the algorithm reduces the
computational burden of evaluating the distance between all road links and
residential buildings. Finally, the nearest road network link is identified by
computing the minimum perpendicular distance from $\mathbf{p_{h}}$ as,

$\mathscr{F}_{M}\left(h\right)=\textrm{arg}~{}\min_{e=\left(u,v\right)\in\\{e_{i}\\}_{i=1}^{k}}\frac{\lvert\lvert\left(\mathbf{p}_{u}-\mathbf{p}_{h}\right)\times\left(\mathbf{p}_{u}-\mathbf{p}_{v}\right)\rvert\rvert}{\lvert\lvert\mathbf{p}_{u}-\mathbf{p}_{v}\rvert\rvert}$
(3)

Algorithm 1 Map residences to the nearest road network link.

Input Set of residences $\mathscr{H}$, road network
$\mathscr{G_{R}}\left(\mathscr{V}_{R},\mathscr{E}_{R}\right)$, radius for
bounding region $r$.

1:for each link $e\in\mathscr{E}_{R}$ do

2: Define bounding region $\mathcal{B}_{e}$ using (1).

3:end for

4:for each residence $h\in\mathscr{H}$ do

5: Define bounding region $\mathcal{B}_{h}$ using (2).

6: Find links $e_{1},\cdots,e_{k}$ with bounding regions
$\mathcal{B}_{e_{1}},\cdots,\mathcal{B}_{e_{k}}$ which intersect with
$\mathcal{B}_{h}$.

7: Compute the nearest link to residence $h$ using (3)

8:end for

Output Mapping $\mathscr{F}_{M}\colon\mathscr{H}\rightarrow\mathscr{E}_{R}$

Implementation examples. The steps involved in mapping a given to the nearest
road network link are shown in Fig. 2. The first figure shows the residence
(magenta), which is required to be mapped to the nearest road link. We draw
bounding regions around each road link as well as the residence. Thereafter,
we identify the bounding regions which intersect with the bounding region of
the residence. The nearest road link is identified from these short-listed
links.

### Step1b: Connect residences to local transformers

######  Problem 1 ($\mathcal{P}_{\textrm{sec}}$ construction).

Given a road link $e\in\mathscr{E}_{R}$ with a set of residences
$\mathscr{F}_{M}^{-1}\left(e\right)$ assigned to it, construct an optimal
forest of trees, $\mathscr{G}_{S}\left(e\right)$, rooted at points (local
transformers) along the link and connecting the residences.

Here we present the detailed optimization framework for the problem
$\mathcal{P}_{\textrm{sec}}\left(e\right)$ to create the secondary
distribution network
$\mathscr{G}_{S}\left(e\right)(\mathscr{V}_{S}\left(e\right),\mathscr{E}_{S}\left(e\right))$
connecting local transformers $\mathscr{V}_{T}\left(e\right)$ along a given
road network link $e\in\mathscr{E}_{R}$ to the set of residences
$\mathscr{V}_{H}=\mathscr{F}_{M}^{-1}\left(e\right)$ mapped to it. It has
already been presented in an earlier work [41]; yet included here for the sake
of completion. We drop the dependency on $e$ from every notation since we are
discussing the same problem for each road link $e$.

Note that we do not have the set of local transformers $\mathscr{V}_{T}$ to
start with. Hence, we begin with a set of probable local transformers
interpolated along the link and denoted by $\mathscr{V}_{\textrm{prob}}$. We
need to connect the set of residences $\mathscr{V}_{H}$ to actual transformer
nodes $\mathscr{V}_{T}\subseteq\mathscr{V}_{\textrm{prob}}$ in a forest of
_star-like_ trees with each tree rooted at a transformer node. A candidate set
of edges $\mathscr{E}_{D}$ between the residences and probable transformer
nodes is chosen. The secondary network edges are to be selected from this
candidate set. We construct an undirected graph
$\mathscr{G}_{D}:=(\mathscr{V}_{D},\mathscr{E}_{D})$ with node set
$\mathscr{V}_{D}=\mathscr{V}_{\textrm{prob}}\cup\mathscr{V}_{H}$ and edge set
$\mathscr{E}_{D}$. We provide the formal problem statement to construct the
secondary distribution network.

######  Formal Problem Statement 1 (Secondary network creation problem).

Given the undirected graph $\mathscr{G}_{D}(\mathscr{V}_{D},\mathscr{E}_{D})$,
find $\mathscr{E}_{S}\subseteq\mathscr{E}_{D}$ such that the induced subgraph
network $\mathscr{G}_{S}(\mathscr{V}_{S},\mathscr{E}_{S})$ with
$\mathscr{V}_{S}=\mathscr{V}_{T}\bigcup\mathscr{V}_{H}$ is a forest of
starlike trees, $\mathscr{V}_{T}\subseteq\mathscr{V}_{\textrm{prob}}$ is the
set of root nodes, and the overall length of the network is minimized.

Note that a complete graph composed of the residence and transformer nodes can
always be considered as the candidate set $\mathscr{E}_{D}$. In this paper, a
Delaunay triangulation [42] of the residential nodes is considered to reduce
the size of the problem.

Table 4: Sets of nodes and edges in secondary network creation problem for each road link Notation | Description
---|---
$\mathscr{E}_{D}$ | Set of all candidate edges for the network
$\mathscr{E}_{S}$ | Set of chosen edges in the network
$\mathscr{V}_{D}$ | Set of all possible nodes in the network
$\mathscr{V}_{\textrm{prob}}$ | Set of all probable transformer nodes along link
$\mathscr{V}_{T}$ | Set of actual transformer nodes along link
$\mathscr{V}_{H}$ | Set of residence nodes mapped to the link

Edge weight assignment An edge $l:(i,j)\in\mathscr{E}_{D}$ is assigned a
weight $w(i,j)$

$w(i,j)=\begin{cases}\infty,\quad\quad\quad\quad\quad\quad\quad\quad\textrm{if
}i,j\in\mathscr{V}_{\textrm{prob}}\\\
\mathsf{dist}(i,j)+\lambda\mathsf{C}(i,j),\quad\textrm{otherwise}\end{cases}$

where
$\mathsf{dist}:\mathscr{V}_{D}\times\mathscr{V}_{D}\rightarrow\mathbb{R}$
denotes the geodesic distance between the nodes $i,j$. The function
$\mathsf{C}(i,j)$ penalizes the cost function if the edge $l$ connecting nodes
$i,j$ crosses the road network link $e$ and is defined as

$\mathsf{C}(i,j)=\begin{cases}0,\quad\textrm{if }i,j\textrm{ are on the same
side of the road network link }e\\\ 2,\quad\textrm{if }i,j\textrm{ are on the
opposite side of the road network link }e\\\ 1,\quad\textrm{if
}i\in\mathscr{V}_{\textrm{prob}}\textrm{ or
}j\in\mathscr{V}_{\textrm{prob}}\end{cases}$

$\lambda$ is a weight factor to penalize multiple crossing of edges over the
road links. It also penalizes multiple edges emerging from the root node. The
weights are stacked in $|\mathscr{E}_{D}|$-length vector $\mathbf{w}$. Note
that an edge between two probable transformer nodes is assigned a weight of
_infinity_ which is equivalent to not considering them as candidate edges.

Edge variables We introduce binary variables $x_{l}\in\\{0,1\\}$ for each
$l\in\mathscr{E}_{D}$. Variable $x_{l}=1$ indicates that the edge $l$ is
present in the optimal topology and $x_{l}=0$ denotes otherwise. Each edge
$l:=(u,v)$ is assigned a flow variable $f_{l}$ (arbitrarily) directed from
node $u$ to node $v$. The binary variable and flows can be respectively
stacked in $|\mathscr{E}_{D}|$-length vectors $\mathbf{x}$ and $\mathbf{f}$.

Node variables The average hourly load demand at the $i^{th}$ residence node
is denoted by $p_{i}$ and is strictly positive. We stack these average hourly
load demands at all residence nodes in a $|\mathscr{V}_{H}|$-length vector
$\mathbf{p}$. We also consider the reactive power load at each household.
Since the models in [37] did not explicitly model the reactive power demand at
each residence, we consider a reactive power consumption of $q_{i}=\gamma
p_{i}$ for each residence $i$. Here, $\gamma=\tan\left(\phi\right)$ with
$\cos\phi$ denoting the power factor. For all cases, we consider a power
factor of $0.95$ which renders $\gamma=0.33$. This is the addition we have
made to the initial models which were proposed in [41].

Degree constraint. Statistical surveys on distribution networks in [38, 43]
show that residences along the secondary network are mostly connected in
series with at most two neighbors. This is ensured by (4) which limits the
degree of residence nodes to $2$.

$\sum_{l:(h,j)}x_{l}\leq 2,\quad\forall h\in\mathscr{V}_{H}$ (4)

Power flow constraints. For the connected graph
$\mathscr{G}_{D}(\mathscr{V}_{D},\mathscr{E}_{D})$, we define the
$|\mathscr{E}_{D}|\times|\mathscr{V}_{D}|$ branch-bus incidence matrix
$\mathbf{A}_{\mathscr{G}_{D}}$ with the entry along $l^{\textrm{th}}$ row and
$k^{\textrm{th}}$ column as

$\begin{matrix}\mathbf{A}_{\mathscr{G}_{D}}(l,k):=\begin{cases}~{}~{}1,\quad
k=i\\\ -1,\quad k=j\\\
~{}~{}0,\quad\textrm{otherwise}\end{cases}&\forall{l=(i,j)\in\mathscr{E}_{D}}\end{matrix}$

Since the order of rows and columns in $\mathbf{A}_{\mathscr{G}_{D}}$ is
arbitrary, we can partition the columns as
$\mathbf{A}_{\mathscr{G}_{D}}=\begin{bmatrix}\mathbf{A_{T}}&\mathbf{A_{H}}\end{bmatrix}$,
without loss of generality, where the partitions are the columns corresponding
to transformer and residence nodes respectively. We call $\mathbf{A_{H}}$ the
reduced branch-bus incidence matrix.

Assuming no network losses, (5a) represents the power balance equations at all
residence nodes. Note that the optimal network is obtained from
$\mathscr{G}_{D}$ after removing the edges for which $x_{l}=0$. Therefore, we
need to enforce zero flows $f_{l}$ for non-existing edges. The constraint (5b)
performs this task along with constraining the flows $f_{l}$ for existing
edges to be within pre-specified capacities $\overline{f}$.

$\displaystyle\mathbf{A_{H}}^{T}\mathbf{f}=\mathbf{p}$ (5a) $\displaystyle-$
$\displaystyle\overline{f}\mathbf{x}\leq\mathbf{f}\leq\overline{f}\mathbf{x}$
(5b)

Ensuring radial topology The radiality requirement of the secondary network
$\mathscr{G}_{S}$ can be enforced from a known graph theory property: _a
forest with $n$ nodes and $m$ root nodes has $n-m$ edges_. In our case,
$|\mathscr{V}_{H}|+|\mathscr{V}_{T}|$ nodes need to be covered in a forest of
trees with $|\mathscr{V}_{T}|$ root nodes, which leads us to the following
constraint.

$\sum_{l\in\mathscr{E}_{D}}x_{l}=|\mathscr{V}_{H}|$ (6)

However, we need to ensure that there are no disconnected cycles in the
optimal network. This can be done by ensuring that the residence points are
connected to a transformer node [30, 31]. In our case, this condition is
satisfied by the node power flow condition in (5a) if all the residential
nodes consume non-zero power.

######  Proposition 1.

The graph $\mathscr{G}_{S}(\mathscr{V}_{S},\mathscr{E}_{S})$ with reduced
branch-bus incidence matrix $\mathbf{A_{H}}$ (corresponding to columns of
$\mathscr{V}_{H}$) and node power demand vector
$\mathbf{p}\in\mathbb{R}^{|\mathscr{V}_{H}|}$, with strictly positive entries,
has exactly $|\mathscr{V}_{T}|$ connected components if and only if there
exists $\mathbf{f}\in\mathbb{R}^{|\mathscr{E}_{D}|}$, such that (5a) is
satisfied.

###### Proof.

Proving by contradiction, suppose
$\mathscr{G}_{S}(\mathscr{V}_{S},\mathscr{E}_{S})$ has more than
$|\mathscr{V}_{T}|$ connected components and there exists
$\mathbf{f}\in\mathbb{R}^{|\mathscr{E}_{D}|}$ satisfying the proposed
equality. Therefore, there exists a connected component
$\mathscr{G}_{C}(\mathscr{V}_{C},\mathscr{E}_{C})$ which is a maximal
connected subgraph with $\mathscr{V}_{C}\subset\mathscr{V}_{H}$ and
$\mathscr{V}_{C}\bigcap\mathscr{V}_{T}=\emptyset$. Let $\mathbf{A}_{C}$ denote
the bus incidence matrix of ${\mathscr{G}}_{C}$. By definition, it holds that
$\mathbf{A}_{C}\mathbf{1}=\mathbf{0}$.

Since graph $\mathscr{G}_{C}(\mathscr{V}_{C},\mathscr{E}_{C})$ is a maximal
connected subgraph of $\mathscr{G}_{S}$, there exists no edge $(i,j)$ with
$i\in\mathscr{V}_{C}$ and $j\in\mathscr{V}_{\overline{C}}$, where
$\mathscr{V}_{\overline{C}}=\mathscr{V}_{S}\setminus\mathscr{V}_{C}$. Since
the order of rows and columns of $\mathbf{A_{H}}$ are arbitrary, we can
partition without loss of generality as

$\mathbf{A_{H}}=\begin{bmatrix}\mathbf{A}_{\overline{C}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{A}_{C}\end{bmatrix}$

We can partition vectors $\mathbf{f}$ and $\mathbf{p}$ conformably to
$\mathbf{A_{H}}$ to get the following equality

$\begin{bmatrix}\mathbf{A}_{\overline{C}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{A}_{C}\end{bmatrix}^{T}\begin{bmatrix}\mathbf{f}_{\overline{C}}\\\
\mathbf{f}_{C}\end{bmatrix}=\begin{bmatrix}\mathbf{p}_{{\overline{C}}}\\\
\mathbf{p}_{C}\end{bmatrix}$

where $\mathbf{p}_{C},\mathbf{p}_{{\overline{C}}}$ are respectively the
vectorized power demands for node sets $\mathscr{V}_{C}$ and
$i\in\mathscr{V}_{\overline{C}}$. From the second block, it implies that

$\mathbf{p}_{C}=\mathbf{A}_{C}^{T}\mathbf{f}_{C}\Rightarrow\mathbf{1}^{T}\mathbf{p}_{C}=\mathbf{1}^{T}\mathbf{A}_{C}^{T}\mathbf{f}_{C}=\mathbf{0}$
(7)

Since all entries of $\mathbf{p}$ are strictly positive from the initial
assumption, we have $\mathbf{1}^{T}\mathbf{p}_{C}\neq 0$ which contradicts (7)
and completes the proof. ∎

Generating optimal network topology The optimal secondary network is obtained
after solving the optimization problem.

$\displaystyle\min_{\mathbf{x}}$ $\displaystyle~{}\mathbf{w}^{T}\mathbf{x}$
(8) s.to. $\displaystyle~{}(\ref{eq:sec-degree}),(\ref{eq:pf-
sec}),(\ref{eq:radial})$

Though the mixed integer optimization problem (MILP) in (8) has been relaxed,
the time complexity depends on the size of the edge set $\mathscr{E}_{D}$. The
computation time can be significantly reduced by considering the edges of a
Delaunay triangulation of nodes $\mathscr{V}_{H}$ as the edge set instead of
choosing $\binom{|\mathscr{V}_{H}|}{2}$ combinatorial edges between the nodes.
The MILP is solved individually for each road link $e\in\mathscr{E}_{R}$ to
construct the secondary network $\mathscr{G}_{S}\left(e\right)$ corresponding
to the road link. The overall secondary network is a combination of such
generated networks.

Implementation examples. Here we present an example of the secondary
distribution network created for the state of Virginia. The steps involved in
constructing the synthetic secondary network connecting residences along a
road network link are shown in Fig. 3. The first figure shows the residences
(red) mapped to the road link and the second figure shows the probable
transformers (green points) along the link. A Delaunay triangulation is
considered to connect the points and obtain the set of possible edges.
Thereafter, the MILP (8) is solved which identifies the optimal set of edges
as shown in the third figure. We observe that all residences are connected in
star-like trees with roots as the transformers along the road link.

Figure 3: Plots showing steps in creating the secondary network for a road
network link. First the probable transformer locations are identified along
the link. Thereafter, the network is generated by solving the optimization
problem (8). The network originates from the road link and connects residences
mapped to it in a forest of star-like trees rooted at the transformers.

### Step 2a: Map local transformers to the nearest substation

The secondary network results in local transformer nodes
$\mathscr{V}_{T}=\bigcup_{e\in\mathscr{E}_{R}}\mathscr{V}_{T}\left(e\right)$
along the road network links of the geographic region (in our case, throughout
the state of Virginia). The task of primary network creation is to connect
these transformer nodes to substation nodes in the region. In order to tackle
the size of the problem, we partition the task into a number of subtasks. We
perform this partition such that each substation is required to connect only
the transformer nodes near to it. Therefore, we produce a map
$\mathscr{F}_{V}$ from the set of local transformers $\mathscr{V}_{T}$ to the
set of substations $\mathscr{S}$. Since we use the road network as a proxy for
creating primary networks, we map each local transformer to the nearest
substation along the road network.

Algorithm 2 Map local transformer to nearest substation.

Input Set of substations $\mathscr{S}$, set of local transformers
$\mathscr{V}_{T}$, road network
$\mathscr{G_{R}}\left(\mathscr{V}_{R},\mathscr{E}_{R}\right)$.

1:for each substation $s\in\mathscr{S}$ do

2: Find the nearest road network node $v_{s}\in\mathscr{V}_{R}$.

3:end for

4:for each local transformer $t\in\mathscr{V}_{T}$ do

5: Compute the nearest substation using (9).

6:end for

Output Mapping $\mathscr{F}_{V}\colon\mathscr{V}_{T}\rightarrow\mathscr{S}$

Algorithm 2 details the steps involved in creating this mapping. First, we
find the geographically located nearest road network node
$v_{s}\in\mathscr{V}_{R}$ for each substation $s\in\mathscr{S}$. Let
$\mathsf{NetDist}(u,v)$ denotes the shortest path distance between nodes
$u,v\in\mathscr{V}_{R}$ along road network $\mathscr{G}_{R}$. We identify the
nearest substation to each local transformer using,

$\mathscr{F}_{V}\left(v\right)=\textrm{arg}\min_{s\in\mathscr{S}}\mathsf{NetDist}\left(v_{s},v\right)$
(9)

Figure 4: Plot showing the mapping between local transformer nodes and
substations in Montgomery County of southwest Virginia. Each color depicts a
particular set of transformer nodes mapped to a substation. Note that same
color has been used to denote more than one mapping. The primary network
creation problem (in the succeeding step) would generate a primary network
connecting the substation to the mapped transformer nodes.

The inverse mapping $\mathscr{F}_{V}^{-1}\left(s\right)$ assigns a set of
local transformer nodes to each substation $s\in\mathscr{S}$. These sets of
local transformer nodes are mutually exclusive and are exhaustive. We perform
the primary network creation for each of these sets individually.

Implementation example. The secondary network creation step generates local
transformer nodes all over the footprint of the region under consideration,
which in our case, is the Montgomery County of southwest Virginia. Next, we
perform the mapping of the transformer nodes to the nearest substation along
the road network. We show the partitions in Fig. 4 where each color denotes a
partition of the transformer nodes.

### Step 2b: Connect local transformers to the substation

######  Problem 2 ($\mathcal{P}_{\textrm{prim}}$ construction).

Given a substation $s\in\mathscr{S}$ with an assigned set of local transformer
nodes $\mathscr{F}_{V}^{-1}\left(s\right)$, construct a tree network
$\mathscr{G}_{P}\left(s\right)$ using the road network $\mathscr{G}_{R}$ as a
proxy which connects all local transformers while ensuring acceptable node
voltages by power engineering standards.

Here we present the detailed optimization framework for the problem
$\mathcal{P}_{\textrm{prim}}\left(s\right)$ to create the primary distribution
network
$\mathscr{G}_{P}\left(s\right)(\mathscr{V}_{P}\left(s\right),\mathscr{E}_{P}\left(s\right))$
connecting substation $s$ to local transformers
$\mathscr{V}_{T}=\mathscr{F}_{V}^{-1}\left(s\right)$ mapped to it. We drop the
dependency on $s$ from every notation since we are discussing the same problem
for each substation $s$.

We start with an undirected graph
$\mathscr{G}_{R}(\mathscr{V},\mathscr{E}_{R})$ which is the road network
subgraph induced from the mapped transformer nodes. Here
$\mathscr{V}=\mathscr{V}_{R}\bigcup\mathscr{V}_{T}$ comprises of transformer
nodes as well as road nodes. The goal of the primary network creation problem
is to select edge set $\mathscr{E}_{P}\subseteq\mathscr{E}_{R}$ to generate
the optimal primary network
$\mathscr{G}_{P}(\mathscr{V}_{P},\mathscr{E}_{P})$. Note that all transformer
nodes are required to be included. On the contrary, road nodes are dummy
points with no load and are only used to connect the local transformers. For
example, we cannot have a road node as a leaf node since it does not connect
transformer nodes. The node set
$\mathscr{V}_{P}=\mathscr{V}_{R}^{\star}\bigcup\mathscr{V}_{T}$ comprises road
and transformer nodes respectively where
$\mathscr{V}_{R}^{\star}\subseteq\mathscr{V}_{R}$ are the selected road nodes.

In all practical distribution networks, multiple feeder lines originate from a
substation and connect the local transformers. To this end, we construct the
primary network $\mathscr{G}_{P}$ as a forest of trees with the root of each
tree connected to the substation through high voltage feeder lines. These root
nodes are road network nodes with no loads connected to them directly. Note
that some nodes in $\mathscr{V}_{R}^{\star}$ which are the root nodes in the
constructed primary network. Here, we present the formal problem statement for
constructing the primary distribution network.

######  Formal Problem Statement 2 (Primary network creation problem).

Given a connected network $\mathscr{G}_{R}(\mathscr{V},\mathscr{E}_{R})$ where
$\mathscr{V}=\mathscr{V}_{R}\bigcup\mathscr{V}_{T}$, find
$\mathscr{E}_{P}\subseteq\mathscr{E}_{R}$ such that the induced subgraph
network $\mathscr{G}_{P}(\mathscr{V}_{P},\mathscr{E}_{P})$ is a forest of
trees with each tree rooted at some $r\in\mathscr{V}_{R}^{\star}$ where
$\mathscr{V}_{P}=\mathscr{V}_{R}^{\star}\bigcup\mathscr{V}_{T}$ and
$\mathscr{V}_{R}^{\star}\subseteq\mathscr{V}_{R}$.

Table 5: Sets of nodes and edges in primary network creation problem Notation | Description
---|---
$\mathscr{E}_{R}$ | Set of all candidate edges in the network
$\mathscr{V}$ | Set of all nodes in the network
$\mathscr{V}_{T}$ | Set of all transformer nodes in the network
$\mathscr{V}_{R}$ | Set of all road nodes in the network
$\mathscr{E}_{P}$ | Set of all chosen edges in the optimal primary network
$\mathscr{V}_{P}$ | Set of all chosen nodes in the optimal primary network
$\mathscr{V}_{R}^{\star}$ | Set of all chosen road nodes in the optimal primary network
Table 6: Binary variables in primary network creation problem Notation | Description
---|---
$x_{e}$ | 1 if edge $e\in\mathscr{E}_{R}$ is included in optimal network
$y_{r}$ | 1 if road node $r\in\mathscr{V}_{R}$ is included in optimal network
$z_{r}$ | 0 if road node $r\in\mathscr{V}_{R}$ is a root node

Binary variables. We assign a binary variable $x_{e}\in\\{0,1\\}$ for each
edge $e\in\mathscr{E}_{R}$. Variable $x_{e}=1$ indicates that the edge $e$ is
included in the network, and $x_{e}=0$ otherwise. Further, we introduce binary
variables $y_{r},z_{r}\in\\{0,1\\}$ for each road network node
$r\in\mathscr{V}_{R}$. Variable $y_{r}=1$ indicates that road node $r$ is part
of the primary network and vice versa. A selected road node ($y_{r}=1$) may be
chosen to be a root node or otherwise. Binary variable $z_{r}=0$ indicates
that road node $r$ is a root node, and $z_{r}=1$ implies that $r$ is not a
root node. Note that an unselected road node ($y_{r}=0$) needs to be treated
as a non-root node ($z_{r}=1$) and hence

$1-z_{r}\leq y_{r},\quad\forall r\in\mathscr{V}_{R}$ (10)

Topology constraints. We need to ensure the following: (i) a selected non-root
road node should not be a terminal node, (ii) the network does not consist of
any cycles, and (iii) the network covers all transformer nodes.

Let $e=(r,j)\in\mathscr{E}_{R}$ denote an edge which is incident on the road
node $r\in\mathscr{V}_{R}$. The degree of a road node $r$ in graph
$\mathscr{G}_{R}$ is denoted by $\sum_{e:(r,j)}x_{e}$. For the first condition
to hold true, we need to ensure that the degree of a selected non-root road
node (with $y_{r}=1,z_{r}=1$) is at least $2$, and the degree of an unselected
road node (with $y_{r}=0$) is $0$. Finally, the degree of a selected road node
which is also a root node (with $y_{r}=1,z_{r}=0$) is positive. This is
ensured through the following inequalities.

$\displaystyle\sum_{e=(r,j)}$ $\displaystyle x_{e}\leq|\mathscr{E}_{R}|y_{r},$
$\displaystyle\forall r\in\mathscr{V}_{R}$ (11a) $\displaystyle\sum_{e:(r,j)}$
$\displaystyle x_{e}\geq 2(y_{r}+z_{r}-1),$ $\displaystyle\forall
r\in\mathscr{V}_{R}$ (11b) $\displaystyle\sum_{e=(r,j)}$ $\displaystyle
x_{e}\geq y_{r},$ $\displaystyle\forall r\in\mathscr{V}_{R}$ (11c)

To enforce the second condition (no cycles in the network) or ‘radiality’
condition, we use results from graph theory. We know that a forest with $n$
nodes and $m$ components has $n-m$ edges. In our case, the total number of
nodes is $|\mathscr{V}_{T}|+\sum_{r\in\mathscr{V}_{R}}y_{r}$, while the number
of components is the number of root nodes, i.e,
$\sum_{r\in\mathscr{V}_{R}}(1-z_{r})$. Therefore, the radiality constraint is
given below

$\sum_{e\in\mathscr{E}_{R}}x_{e}=|\mathscr{V}_{T}|+\sum_{r\in\mathscr{V}_{R}}y_{r}-\sum_{r\in\mathscr{V}_{R}}(1-z_{r})$
(12)

However, this is not a sufficient condition for radiality since it does not
avoid the formation of disconnected components with cycles. We can extend
Proposition 1 in [41] and ensure radiality with the power balance equations at
each node in the network.

Table 7: Power flow variables and parameters in primary network creation problem Notation | Description
---|---
$V_{j}=v_{j}e^{\mathsf{j}\theta_{j}}$ | complex bus voltage phasor at bus $j\in\mathscr{V}$; magnitude: $v_{j}$, angle: $\theta_{j}$
$I_{e}$ | complex current flowing through edge $e\in\mathscr{E}_{R}$
$S_{ij}=P_{ij}+\mathsf{j}Q_{ij}$ | complex power flowing through edge $e=(i,j)\in\mathscr{E}_{R}$ from node $i$ to node $j$; real power flow: $P_{ij}$, reactive power flow: $Q_{ij}$
$s_{j}=p_{j}+\mathsf{j}q_{j}$ | complex power injection at bus $j\in\mathscr{V}$; real power injection: $p_{j}$, reactive power injection: $q_{j}$
$Z_{e}=R_{e}+\mathsf{j}X_{e}$ | complex impedance of edge $e\in\mathscr{E}_{R}$; resistance: $R_{e}$, reactance: $X_{e}$
$f_{e}$ | power flow limit in edge $e\in\mathscr{E}_{R}$
$\overline{s}_{j}$ | substation feeder rating
$\underline{v},\overline{v}$ | lower and upper limits of bus voltage

Power balance constraints. The complex power injection at bus
$j\in\mathscr{V}$ is denoted by $s_{j}=p_{j}+\mathsf{j}q_{j}$. Note that the
power injections are the difference between generation and load demand at each
bus. For transformer nodes, the power injection is the negative load demand;
for non-root road nodes, the power injection is zero; and for root nodes, the
power injection is bounded by a defined limit. It is assumed that the root
nodes are connected to the substation through high voltage feeders. Therefore,
the feeder capacity can be considered as the defined limit.

$\displaystyle\sum_{e:(i,j)}x_{e}S_{ij}-\sum_{e:(j,k)}x_{e}S_{jk}=-s_{j},\quad$
$\displaystyle\forall j\in\mathscr{V}_{T}$ (13a)
$\displaystyle\sum_{e:(i,j)}x_{e}S_{ij}-\sum_{e:(j,k)}x_{e}S_{jk}=0,\quad$
$\displaystyle\forall j\in\mathscr{V}_{R},z_{j}=1$ (13b)
$\displaystyle\bigg{|}\sum_{e:(i,j)}x_{e}S_{ij}-\sum_{e:(j,k)}x_{e}S_{jk}\bigg{|}\leq\overline{s}_{j},\quad$
$\displaystyle\forall j\in\mathscr{V}_{R},z_{j}=0$ (13c)

In the secondary network creation method, we discuss about the residential
reactive power model. We assumed a constant power factor of $0.95$ for all
residences. This leads to a linear relationship between the real and reactive
power throughout the network with $Q_{ij}=\gamma P_{ij}$ for all edges. We can
drop either of the real or reactive power terms in (13), since they would lead
to the same equality constraints. Such assumptions are standard for problems
related to network planning [38, 30, 31]. Therefore, we can rewrite (13) by
dropping the reactive power terms (imaginary parts) and combining (13b) and
(13c) to get:

$\displaystyle\sum_{e:(i,j)}P_{ij}-\sum_{e:(j,k)}P_{jk}=-p_{j},\quad\forall
j\in\mathscr{V}_{T}$ (14a)
$\displaystyle-\overline{p}_{j}(1-z_{j})\leq\sum_{e:(i,j)}P_{ij}-\sum_{e:(j,k)}P_{jk}\leq\overline{p}_{j}(1-z_{j}),\forall
j\in\mathscr{V}_{R}$ (14b)

Power flow constraints. We define complex node voltage
$V_{j}=v_{j}e^{\mathsf{j}\theta_{j}}$ (magnitude $v_{j}$ and angle
$\theta_{j}$) for each node $j\in\mathscr{V}$ and complex power
$S_{i,j}=P_{i,j}+\mathsf{j}Q_{i,j}$ for each edge $e=(i,j)\in\mathscr{E}_{R}$
flowing from node $i$ to $j$, where $P_{i,j}$ and $Q_{i,j}$ respectively
denote real and reactive power flowing along edge $e\in\mathscr{E}_{R}$. The
current flowing through the edge is defined by $I_{e}$ and the complex
impedance of the edge is denoted by $Z_{e}=R_{e}+\mathsf{j}X_{e}$. The
following equations/inequalities describe the relationship between the
variables and the associated constraints. Note that the constraints are only
activated for those edges which are selected in the optimal network.
Therefore, we have introduced the binary variable $x_{e}$ in (15a). We assume
that HV feeder lines from the substation to the root nodes in the optimal
primary distribution network end in voltage regulators, which ensure that the
root nodes have a voltage of $1.0$pu.

$\displaystyle x_{e}(V_{i}-V_{j}-Z_{e}I_{e})=0,\quad$ $\displaystyle\forall
e\in\mathscr{E}_{R}$ (15a) $\displaystyle S_{ij}=V_{j}I_{e}^{\star},\quad$
$\displaystyle\forall e\in\mathscr{E}_{R}$ (15b)
$\displaystyle|S_{ij}|\leq\overline{f}_{e},\quad$ $\displaystyle\forall
e\in\mathscr{E}_{R}$ (15c) $\displaystyle\underline{v}\leq
v_{j}\leq\overline{v},\quad$ $\displaystyle\forall j\in\mathscr{V}$ (15d)
$\displaystyle v_{j}=1,\quad$ $\displaystyle\forall j\in\mathscr{V},z_{j}=0$
(15e)

Generating optimal primary network. Each edge $e=(i,j)\in\mathscr{E}_{R}$ is
assigned a weight $w_{e}=w(i,j)=\mathsf{dist}(i,j)$ which is the geodesic
distance between the nodes. Additionally, for every road node
$j\in\mathscr{V}_{R}$, we compute its geodesic distance from the substation
$s$, denoted by $\mathsf{d}_{j}$. The optimal primary network topology is
obtained by solving the optimization problem:

$\displaystyle\min_{\mathbf{x},\mathbf{y},\mathbf{z}}$
$\displaystyle~{}\sum_{e\in\mathscr{E}_{R}}x_{e}w_{e}+\sum_{j\in\mathscr{V}_{R}}(1-z_{j})\mathsf{d}_{j}$
(16) s.to. $\displaystyle~{}(\ref{eq:sup-non-root}),(\ref{eq:sup-prim-
connectivity}),(\ref{eq:sup_prim_radial}),(\ref{eq:sup-prim-pbal-
relax}),(\ref{eq:sup-prim-pf})$

Relaxing non-linear power flow constraints. Note that power flow constraints
in (15) are non-convex because of the quadratic equality constraint in (15b)
and the bilinear terms in (15a). First, we consider relaxing the quadratic
equality constraint. Multiplying by the complex conjugates in (15a) and (15b),
we get the following result:

$v_{j}^{2}=v_{i}^{2}-2R_{e}P_{ij}-2X_{e}Q_{ij}+(R_{e}^{2}+X_{e}^{2})|I_{e}|^{2}$
(17)

Assuming small $R_{e},X_{e}$, the current magnitudes can be eliminated. Note
that (17) is still non-linear in $v_{i},v_{j}$; however, it is linear in
squared magnitude. This relaxed model is known in the literature as the
Linearized Distribution Flow (LDF) model. The squared voltage is often
approximated as $|v_{i}|^{2}\approx 2v_{i}-1$ which leads to the relation
$v_{j}=v_{i}-R_{e}P_{i,j}-X_{e}Q_{i,j}$ [44]. Notice that the constraint in
(17) is required to be enforced on only those edges $e\in\mathscr{E}_{R}$ for
which $x_{e}=1$, which leads us to the relaxed version of (15a) as

$x_{e}(v_{i}-v_{j}-R_{e}P_{ij}-X_{e}Q_{ij})=0\;.$ (18)

Here (18) is non-convex due to the bi-linear terms in the equality. In order
to deal with this non-convexity, McCormick relaxation has been used widely in
several previous works [29, 30]. In general, McCormick relaxation replaces the
non-convex equality constraint with its convex envelope [45]. However, in the
case of bilinear variables with at least one binary variable, this relaxation
becomes exact. The convex relaxed version of (15) is:

$\displaystyle-$ $\displaystyle(1-x_{e})M\leq
v_{i}-v_{j}-R_{e}P_{ij}-X_{e}Q_{ij}\leq(1-x_{e})M,~{}\forall
e\in\mathscr{E}_{R}$ (19a) $\displaystyle-$
$\displaystyle\overline{f}_{e}{x_{e}}\leq{P_{ij}}\leq\overline{f}_{e}\mathbf{x_{e}},~{}\quad\forall
e\in\mathscr{E}_{R}$ (19b) $\displaystyle(1-v_{r})\leq z_{r},~{}\quad\forall
r\in\mathscr{V}_{R}$ (19c)
$\displaystyle\underline{v}\leq{v_{j}}\leq\overline{v},~{}\quad\forall
j\in\mathscr{V}$ (19d)

The overall relaxed optimization problem for creating the primary network
topology is:

$\displaystyle\min_{\mathbf{x},\mathbf{y},\mathbf{z}}$
$\displaystyle~{}\sum_{e\in\mathscr{E}_{R}}x_{e}w_{e}+\sum_{j\in\mathscr{V}_{R}}(1-z_{j})d_{j}$
(20) s.to. $\displaystyle~{}(\ref{eq:sup-non-root}),(\ref{eq:sup-prim-
connectivity}),(\ref{eq:sup_prim_radial}),(\ref{eq:sup-prim-pf-
relax}),(\ref{eq:sup-prim-pbal-relax})$

Implementation example. Next, we show the steps involved in creating the
primary distribution network within the boundary of Montgomery County of
Virginia, US in Fig. 5. We use the road network obtained from Open Street Maps
as a proxy for the primary network, i.e., the edges in the primary network are
chosen from the road network edges by solving the optimization problem in
(20). The first figure shows the road network which is used to construct the
primary network. The zoomed-in inset figures show two examples of road network
subgraphs assigned to two substations. Note that this assignment is computed
from the inverse mapping $\mathscr{F}_{V}^{-1}$ in the preceding step. The
second figure shows the solution of the optimization problem where the primary
network is created from the set of edges in the road network. The primary
network is a tree originating from a substation and connects the local
transformers through multiple feeders. Note that the road network in the first
figure had multiple cycles, which are no longer present in the primary network
in the second figure. Finally, the secondary network is appended to the local
transformers to obtain the final distribution network, which is shown in the
third figure.

Figure 5: Plots showing the steps involved in creating the primary network of
Montgomery County in southwest Virginia. The edges of the primary network
(middle figure) are chosen from the road network (left figure) by solving the
optimization problem in (20). The network originates from the substations and
follows the road links to connect the local transformers in a tree structure
with no loops. The secondary network is appended to the local transformers to
obtain the final synthetic distribution network (right figure).

### Step 3: Constructing an ensemble of networks

Figure 6: Plots summarizing the optimal distribution network constructed using
Steps 1 and 2. The overall network (left figure) consists of the primary
network (blue) and secondary network (red). The optimal primary network
(middle figure) is constructed from the underlying road network (right
figure). The edges in the primary network are a subset of the road network
graph edges.

Using the optimization frameworks in (8) and (20), we create the optimal
(sometimes near-optimal) distribution network as shown in the left figure of
Fig. 6. It consists of the primary network which connects the substation to
the local transformers in a tree structure, and these transformers are
connected to the residences in chains, forming the secondary network. The
middle figure of Fig. 6 shows the near-optimal primary network denoted by
$\mathscr{G}_{P}^{0}(\mathscr{V}_{P}^{0},\mathscr{E}_{P}^{0})$. Note that the
MILP (20) identified
$\mathscr{G}_{P}^{0}(\mathscr{V}_{P}^{0},\mathscr{E}_{P}^{0})$ as the solution
by choosing edges from the underlying road network
$\mathscr{G}_{R}(\mathscr{V}_{R},\mathscr{E}_{R})$. The right figure of Fig. 6
shows the underlying road network corresponding to the primary network. We
note that all primary network edges are selected from the road network graph.
The resulting primary network is a tree, and it connects all the transformer
nodes. The road network nodes $\\{3,5,6\\}$ are used to construct the network,
whereas some other road network nodes such as $\\{18,19\\}$ are not included
in the generated primary network.

We model the variant network creation process as a Markov chain $\mathcal{M}$
where each state denotes a valid realization of the primary network. The
transitions
$\mathscr{G}_{P}^{t}(\mathscr{V}_{P}^{t},\mathscr{E}_{P}^{t})\rightarrow\mathscr{G}_{P}^{t+1}(\mathscr{V}_{P}^{t+1},\mathscr{E}_{P}^{t+1})$
of $\mathcal{M}$ are specified in the following manner. Let $s$ denote the
substation node in the network
$\mathscr{G}_{P}^{t}(\mathscr{V}_{P}^{t},\mathscr{E}_{P}^{t})$ and
$\mathscr{N}(s)$ denote the neighbors of $s$. We choose a random set of edges
$\mathscr{E}_{D}=\\{e_{d}\in\mathscr{E}_{P}^{t}|e_{d}=(u,v),u\neq v\neq s\\}$
to be deleted and thereafter solve a restricted version of (20) with the
following additional constraints.

$\displaystyle x_{e}=0,\textrm{ for }e\in\mathscr{E}_{D};\quad x_{e}=1,\forall
e\in\mathscr{E}_{P}^{t}\setminus\mathscr{E}_{D}$ (21a) $\displaystyle
y_{r}=1,\forall r\in\mathscr{V}_{R}\cap\mathscr{V}_{P}^{t}$ (21b)
$\displaystyle z_{r}=0,\forall r\in\mathscr{N}(s);\quad z_{r}=1,\forall
r\in\mathscr{V}_{R}\setminus\mathscr{N}(s)$ (21c)

Algorithm 3 Constructing ensemble of networks.

Input Near-optimal primary network
$\mathscr{G}_{P}^{0}(\mathscr{V}_{P}^{0},\mathscr{E}_{P}^{0})$, optimal
secondary network $\mathscr{G}_{S}$, underlying road network
$\mathscr{G}_{R}(\mathscr{V}_{R},\mathscr{E}_{R})$, size of ensemble $N$.

1:Initialize network count, $t\leftarrow 0$

2:while $t<N$ do

3: Choose a random set of edges $\mathscr{E}_{D}\subset\mathscr{E}_{P}^{t}$ to
be deleted.

4: Formulate optimization problem $\mathcal{P}_{t}$ by combining (20) with
constraints in (21).

5: if $\mathcal{P}_{t}$ is feasible then

6: Solve $\mathcal{P}_{t}$ to get the variant primary network
$\mathscr{G}_{P}^{t+1}(\mathscr{V}_{P}^{t+1},\mathscr{E}_{P}^{t+1})$

7: Increment network count $t\leftarrow t+1$.

8: end if

9: Augment the secondary network to the variant primary network;
$\mathscr{G}^{t}\leftarrow\mathscr{G}_{S}\bigcup\mathscr{G}_{P}^{t}(\mathscr{V}_{P}^{t},\mathscr{E}_{P}^{t})$

10:end while

Output Ensemble of networks
$\mathscr{G}^{1},\mathscr{G}^{2},\cdots,\mathscr{G}^{N}$

Figure 7: Schematic of the Markov chain process to generate an ensemble of
primary networks from the optimal primary distribution network. Every
transition in the Markov chain involves deletion of a random edge from the
primary network followed by solving a restricted version of the primary
network creation problem. If no feasible solution exists for the restricted
version of the problem, the chain transitions to itself.

The restrictions reduce the order of the optimization problem (20) to a
significant extent and thereby requires less time to reach an optimal
solution. (21a) forces all edges in the set
$\mathscr{E}_{P}^{t}\setminus\mathscr{E}_{D}$ to be selected in the new
network $\mathscr{G}_{P}^{t+1}(\mathscr{V}_{P}^{t+1},\mathscr{E}_{P}^{t+1})$.
The edge set $\mathscr{E}_{D}$ is not selected and the optimization problem
chooses edges from $\mathscr{E}_{R}\setminus\mathscr{E}_{P}^{t}$. Further,
(21b) forces all road network nodes in $\mathscr{G}_{P}^{t}$ to be selected in
$\mathscr{G}_{P}^{t+1}$; the other road nodes are left free to be selected.
Finally (21c) keeps the feeder nodes (which are connected to substation node
through HV feeders) similar in $\mathscr{G}_{P}^{t}$ and
$\mathscr{G}_{P}^{t+1}$. Therefore, the number of binary variables reduces
from $|\mathscr{E}_{R}|+2|\mathscr{V}_{R}|$ to
$\left|\mathscr{E}_{R}\setminus\mathscr{E}_{t}\right|+\left|\mathscr{V}_{R}\setminus\mathscr{V}_{t}\right|$,
which is a significant improvement. Here $|\cdot|$ denotes the cardinality of
a set. Algorithm 3 lists the steps involved in constructing the ensemble of
$N$ networks. We start with the near-optimal primary network and create the
variant primary networks. Finally, we augment the variant primary networks to
the optimal secondary network to create the networks in the ensemble.

The schematic of the process is shown in Fig. 7 where each transition involves
deletion of a random edge followed by solving the restricted version of the
optimization problem. For example, in case of the first transition, the edge
$(4,10)$ is chosen at random to be deleted. The restricted version of the
optimization problem forces all the edges in the optimal primary network
except $(4,10)$ to be selected. The remaining edges are selected from the
candidate set of edges (shown in dotted green lines) so that the output result
is a feasible network.

### Post-processing step

The final step of the synthetic network creation process is to assign _line
types_ to the different distribution lines in the network. We use the catalog
of distribution lines [46] to assign the line types and parameters to
secondary and primary distribution lines. A simultaneity factor of $0.8$ has
been used to rate the transformers and distribution lines similar to the steps
used in [20]. This means that when all residences are consuming their
respective peak hourly load, the distribution lines are loaded to at most
$80\%$ of their rating. Table 8 lists the node and edge attributes in the
generated synthetic power distribution networks.

Figure 8: Plots showing histogram of node voltages for three synthetic
distribution networks operating at peak load. We observe that majority of the
nodes are within acceptable voltage limits set by the ANSI C.84 A standard.
The substation voltage is set at 1.02 p.u. to avoid low voltages at the leaf
residence nodes.

Figure 9: Plots showing variation of voltage with distance from substation
(bottom three figures) for three synthetic distribution networks operating at
peak load. The HV feeders are used to connect distant residences while
maintaining a healthy voltage profile. The primary network maintains the
voltage level within acceptable engineering standard for distribution; while
significant voltage drop occurs at in the LV secondary network.

A key aspect of creating the synthetic networks is regulating the voltage at
different nodes in the network. ANSI C.84 Range A limits the acceptable node
voltage within $0.95-1.05$ p.u.; however, in certain sections, we observe
undervoltage. Since such occurrences of undervoltage are common in real
distribution networks as well, we relax the acceptable limits as follows:
accept networks where voltage at $>99.5\%$ of residences are within the limits
set by ANSI C.84 Range A ($0.95-1.05$ p.u.). In the rejected cases, we perform
either of the following post-processing steps: (i) increase the load tap
changer (LTC) setting at the distribution setting to 1.02 p.u. (ii) add line
regulators at the end of HV feeders, which are used to connect remote loads.
The three figures in Fig. 8 show the histogram of node voltages after
performing these post-processing steps. The three plots in Fig. 9 show the
variation of node voltage with distance from the substation. The blue lines
show the HV feeder lines, which results in the minimum voltage drop. The black
lines denote the MV primary network lines, which causes voltage drop as we
tend to move away from the substation. Finally, the red lines are the LV
secondary lines resulting in a major voltage drop at the residences. However,
in all cases, we notice that $>99.5\%$ of customer node voltages are within
the acceptable ANSI C.84 A limits. The details of the output synthetic power
distribution networks are provided in Table 8. We have uploaded our data to
the GitHub repository [47]. Table 9 provides a list of lines used in the
synthetic secondary and primary distribution networks. A summary of the
percentage of these lines in the entire network of Virginia is also appended
in the table.

Table 8: Node and edge attributes in created synthetic power distribution networks | Attribute Name | Attribute Details |  | Size of Output Network for
---
Montgomery County, Virginia (US)
| Node ID | integer ID of node | 63220 nodes
| Node geometry |  $(x,y)$: longitude and latitude of node | 20 substation nodes
Nodes | Node label | ‘S’: substation, ‘T’: local transformer, ‘H’: residence, ‘R’: auxiliary node | 18759 transformer nodes
| Node average load | Average hourly load demand in Watts (only for residences) | 35629 residence nodes
| Node peak load | Peak hourly demand in Watts (only for residences) | 8812 auxiliary nodes
| Node phase | Phase assigned to node (A, B or C) |
| Edge ID | integer IDs of connecting nodes |
| Edge geometry | shapely LineString geometry of edge | 63200 edges
| Edge label | ‘E’: HV feeder line, ‘P’: MV primary network, ‘S’: LV secondary network | 156 HV feeder lines
Edges | Edge name | Name of conductor type used for the edge | 27415 primary edges
| Length | Length of the edge in meters | 35629 secondary edges
| Resistance | Resistance of the edge conductor in p.u. |
| Reactance | Reactance of the edge conductor in p.u. |
| Edge phase | Phase assigned to edge (A, B, C or mixed) |
Table 9: Catalog of LV and MV distribution network lines Line Name | Line Type |  | Resistance
---
(ohms/$1000$ ft)
| Reactance
---
(ohms/$1000$ ft)
| Current
---
(A)
| Voltage
---
(kV)
| Percentage
---
of lines
in Virginia (%)
OH Voluta | TRPLX #6 | 0.661 | 0.033 | 95 | 0.24 | 73.40
OH Periwinkle | TRPLX #4 | 0.416 | 0.031 | 125 | 0.24 | 8.70
OH Conch | TRPLX #2 | 0.261 | 0.030 | 165 | 0.24 | 7.00
OH Neritina | TRPLX 1/0 | 0.164 | 0.030 | 220 | 0.24 | 5.00
OH Runcina | TRPLX 2/0 | 0.130 | 0.029 | 265 | 0.24 | 2.20
OH Zuzara | TRPLX 4/0 | 0.082 | 0.027 | 350 | 0.24 | 3.70
OH Swanate | ACSR #4 | 0.407 | 0.113 | 140 | 12.47 | 90.80
OH Sparrow | ACSR #2 | 0.259 | 0.110 | 185 | 12.47 | 3.90
OH Raven | ACSR 1/0 | 0.163 | 0.104 | 240 | 12.47 | 2.70
OH Pigeon | ACSR 3/0 | 0.103 | 0.099 | 315 | 12.47 | 1.70
OH Penguin | ACSR 4/0 | 0.082 | 0.096 | 365 | 12.47 | 0.80

### Three phase networks

We also created three phase synthetic distribution networks as a post-
processing step. The synthetic networks created thus far are positive sequence
networks, which are suitable for planning studies. We create networks
comprising of three phases (A, B, and C) which would be useful to run studies
on the operation of the network involving phase unbalances.We aim to create a
balanced three phase network where the power delivered at each substation
feeder is balanced in the three phases. For the sake of simplicity, we
consider branching of three different phases only in secondary networks. This
implies that each residence is assigned one of the three phases, and the
branching occurs at the local pole top transformers. Each branch originating
from a transformer feeds residences with the same phase. Thus, we have a
problem of three way partitioning the set of residences, such that the total
load supplied to residences in each partition is balanced.

######  Problem 3 (Three phase network creation problem).

Given a set of residences $\mathscr{H}$ connected to a substation feeder with
respective load demand $p_{h}$ for every $h\in\mathscr{H}$, find a partition
into three sets $\mathscr{H}_{A},\mathscr{H}_{B},\mathscr{H}_{C}$ such that

$\sum_{h\in\mathscr{H}_{A}}p_{h}=\sum_{h\in\mathscr{H}_{B}}p_{h}=\sum_{h\in\mathscr{H}_{C}}p_{h}$

We define an optimization problem similar to the one proposed in [22], which
assigns a phase to every residence connected to a substation feeder and
minimizes the deviation between the total load fed by each pair of phases. We
provide the details of the optimization framework, which we solve to assign
one of the three phases to each residence connected to a substation feeder. We
assign binary variables $u^{A}_{h},u^{B}_{h},u^{C}_{h}\in\\{0,1\\}$ to every
residence $h\in\mathscr{H}$ in order to denote the phase assigned to the
residence. For example, $u^{A}_{h}=1$ denotes that residence $h$ is assigned
phase A. We stack the binary variables to construct vectors
$\mathbf{u}^{A},\mathbf{u}^{B},\mathbf{u}^{C}\in\\{0,1\\}^{|\mathscr{H}|}$
respectively. Since we assign exactly one of the phases (A,B,C) to each
feeder, we enforce the following constraint.

$\mathbf{u}^{A}+\mathbf{u}^{B}+\mathbf{u}^{C}=\mathbf{1}$ (22)

Let $p_{n}$ denote the average hourly power demand at residence $n$. We stack
the power demand to obtain the vector
$\mathbf{p}\in\mathbb{R}^{|\mathscr{H}|}$. Hence, we can compute the total
power fed by three phases at the substation feeder as follows.

$\displaystyle p^{A}=\mathbf{p}^{T}\mathbf{u}^{A}$ (23a) $\displaystyle
p^{B}=\mathbf{p}^{T}\mathbf{u}^{B}$ (23b) $\displaystyle
p^{C}=\mathbf{p}^{T}\mathbf{u}^{C}$ (23c)

The phases are assigned to each residence by solving the optimization problem
in (24).

$\displaystyle\min_{\mathbf{u}^{A},\mathbf{u}^{B},\mathbf{u}^{C}}$
$\displaystyle~{}|p^{A}-p^{B}|+|p^{B}-p^{C}|+|p^{C}-p^{A}|$ (24) s.to.
$\displaystyle~{}(\ref{eq:3phase-sum}),(\ref{eq:phase-load})$

Figure 10: Plots showing two 3-phase networks generated from the optimal
positive sequence networks.

## Results

### Motifs in the network

Figure 11: Plots showing number of 4-node paths (top left), 4-node star motifs
(top right), 5-node paths (bottom left), and 5-node star motifs (bottom right)
as a function of network size (measured as number of nodes in the network).
Colors depict motif numbers in urban versus rural areas. Urban distribution
networks have a larger number of star motifs than rural networks. In contrast,
the path motif count does not differ significantly across rural and urban
areas. Urban networks are often larger than rural networks as measured by
number of nodes due to larger population size.

We consider subgraphs of sizes $4$ and $5$ in the generated synthetic power
distribution networks. Since the distribution networks are tree graphs, the
relevant subgraphs are paths and stars. We define a $k$-path motif as a
subgraph of $k$ nodes that form a path. A $k$-star motif is a subgraph of $k$
nodes which form a star, i.e., it consists of a single central node with
degree $k-1$, and the remaining $k-1$ nodes are connected to the central node.
We compare the path and star motifs in the created synthetic distribution
networks in Fig. 11. The two colors denote the motif counts for rural and
urban networks.

### Structural Comparison of Network Ensembles

We consider the structural comparison of the ensemble of networks here. We
create an ensemble of $20$ networks for each distribution network created for
Montgomery County in southwest Virginia. Fig. 12 shows the results of
structural comparison among $4$ of these networks. Note that the networks
consist of the same residences, transformers and substations and hence the
number of nodes and edges are similar to each other. However, the geometries
of each network differ from the others. We do not see many structural
differences between the networks at a visual level. However, a comparison of
the geometry of the networks in the rectangular grid cells shows some
differences, but these are not too large, as that would create networks where
the power engineering constraints do not hold true.

Figure 12: Plots comparing the Hausdorff distance between pairs of networks in
the ensemble. The comparison is performed for each rectangular grid cell and
color depicts the magnitude of deviation in meters. Figure 13: Plot comparing
the structure of the actual distribution network (red) with the synthetic
distribution network (blue) for the town of Blacksburg in southwest Virginia.
The black dotted network denotes the underlying road network, which is used as
a proxy to create the synthetic network.

## Validation

### Variation in network structure

In this section, we perform a visual comparison of the generated synthetic
network to the actual network covering the same geographical region. Fig. 13
shows the actual distribution network of the town of Blacksburg in southwest
Virginia along with the synthetic network generated for the same. At first
sight, the two networks are adjacent to and almost overlap each other. The
inset figure confirms this hypothesis. A zoomed-in view of the figure shows
that the two networks are similar to each other. While the synthetic network
retraces the road network, as expected based on our assumption, the actual
network is also adjacent to it. This validates our primary assumption that the
primary distribution network follows the road network.

Feeder selection. The significant structural difference between the two
networks is the substation feeder to which it is connected. Fig 14 depicts the
structural differences between the two networks. The actual network is
connected to the nearest geographically located substation, where as the
synthetic network is connected to a distant feeder. It is to be noted that the
synthetic distribution network is generated using first principles with the
assumption that it follows the road network to the maximum extent (it is
easier to place distribution poles along the road network). It is seen that
though the nearest substation (from which the actual network is fed) has no
road links connecting itself to the neighborhood. On the contrary, the other
substation, though located at a further distance, has road links connecting
itself to the neighborhood under consideration. Therefore, the synthetic
network generation algorithm prefers the latter substation over the former.

Figure 14: Plot explaining the visual structural difference between actual
distribution network (red) and synthetic distribution network (blue) generated
using first principles from the knowledge of road network (black).

### Statistical Validation

This section contains additional results related to the statistical validation
of created synthetic networks. The created networks of the town of Blacksburg
are compared to the actual physical counterpart obtained from a power company.

Degree Distribution. The degree $k_{i}$ of a node $i$ in a graph with $n$
nodes and internode adjacency matrix $\mathbf{A}$ with elements $a_{ij}$ is
defined by $k_{i}=\sum_{j=1}^{n}a_{ij}$. In general, the node degree computes
the number of lines connected to it in the network. We consider three separate
sections of the two networks and compare the degree distribution for each
section. We observe that the two distributions match each other very closely.
We observe that the degree of the nodes ranges between $1$ and $4$ for either
network and the degree distribution shows that majority of nodes have a degree
of $2$. The networks are created in a manner so that when operating with
average hourly load demand, the network maintains an acceptable voltage
profile. Every branch originating from a node leads to a small voltage drop
caused by the load demand of the children nodes. This, in turn, causes a
voltage drop in the parent nodes. Hence, for a feasible network to operate
within acceptable voltage limits, it is expected that multiple branching from
a single node is avoided as much as possible. This is the same observation we
notice in the actual and synthetic networks, where the node degree hardly
exceeds a degree of $3$.

Figure 15: Plots comparing degree distribution of the actual and synthetic
networks for three sections of region.

Hop Distribution. A _path_ $\mathscr{P}_{ij}$ in a graph between two nodes $i$
and $j$ is defined as a sequence of adjacent edges starting from node $i$ and
terminating at node $j$. A _tree_ network (devoid of any cycles) has a unique
path between any pair of nodes. The hops between nodes $i$ and $j$ in a tree
network is defined as $|\mathscr{P}_{ij}|$ where $|\cdot|$ denotes the
cardinality of a set. Essentially, it denotes the number of edges between the
nodes $i$ and $j$. Here we are interested in the number of edges between any
node $i$ and the substation node, which is the ‘root’ node in the synthetic
networks. Therefore, we define the number of hops between any node $i$ and
substation (root) node by $h_{i}=|\mathscr{P}_{ri}|$. We consider the
empirical distribution for the hops $h_{i}$’s of all nodes and term it as the
_hop distribution_ of the networks. In this analysis, we study how the nodes
in each network are distributed around the root node. We observe that the hop
distributions for the synthetic and actual networks are significantly
different. This is primarily because of the dissimilarity in the spatial
distribution of nodes in each network which is studied later on in the later
section.

Figure 16: Plots comparing hop distribution of the actual and synthetic
networks for three sections of the geographical region.

Figure 17: Plots comparing distribution of distance of each node to root in
the actual and synthetic networks for three sections.

Reach Distribution. Since the synthetic networks have a geographic attribute
associated with them, we can use the geographic distance of each node from the
substation node to study how the nodes in the network are physically present
around the root node. Each edge $e$ has an attribute $w_{e}$ which denotes the
length of the edge in meters. We define _reach_ of a node $i$ in the synthetic
network as $l_{i}=\sum_{e\in\mathscr{P}_{ri}}w_{e}$. Essentially, the reach of
a node denotes the physical distance between the node and the substation
(root) node. We consider the empirical distribution for the reach $l_{i}$’s of
all nodes and term it as the _reach distribution_ of the networks. Fig. 17
compares the reach distribution of synthetic and actual network for three
different sections. An interesting observation is that the actual network
consists of a large number of long edges. This is primarily because of the
fact that the network has been built over multiple years as the population
grew in the geographical location. However, the synthetic networks are
generated using first principles and as an output of an optimization problem
where the generated network has the minimum length (or requires minimal
installation and maintenance cost).

Figure 18: Plots comparing motif count to network size ratio of actual and
synthetic distribution networks over three different sections.

Motif counts. We also compare the 4-node and 5-node path and star motifs in
actual and synthetic distribution networks. We consider three different
sections in the network. Since the size of the actual and synthetic networks
are very different, we compare the ratio of motif counts to the network size
in Fig. 18. We observe that the ratios are very similar for all three sections
of the networks. This validates the structural resemblance of the networks
through statistical attributes.

### Structural Validation

This section contains additional results related to the structural validation
of created synthetic networks. The created networks of the town of Blacksburg
are compared to the actual physical counterpart.

Figure 19: Plots comparing spatial distribution of nodes in the actual and
synthetic networks. The comparison is made for rectangular grid cells where
data regarding actual network is available. Two different grid resolutions are
considered for the spatial distribution comparison. Also we compare the node
distribution by moving the horizontal grid demarcations and observe varied
results.

Problems related to structural validation. The actual distribution network has
been obtained from a power distribution company as images. The associated
geographic topology is obtained by overlaying these images on maps from
OpenStreetMap. The distribution network in the considered region is owned by
multiple power companies and we have obtained the network for only one of
them. For this reason, we focus our comparison on only those portions where we
have data pertaining to both networks. Further, the above-mentioned overlaying
process has been performed manually and thereby has introduced some errors in
assigning the geographical attributes to the actual network. To this end, the
effective way to compare the structural attributes of the networks is to
divide the entire geographic region into multiple rectangular grid cells and
perform a comparison in each cell separately.

Spatial distributions of nodes. We compare the spatial distribution of nodes
in the actual and synthetic networks as follows. Let the actual and synthetic
networks for a region be denoted by
$\mathscr{G}_{\textrm{act}}(\mathscr{V}_{\textrm{act}},\mathscr{E}_{\textrm{act}})$
and
$\mathscr{G}_{\textrm{syn}}(\mathscr{V}_{\textrm{syn}},\mathscr{E}_{\textrm{syn}})$,
and let
$\mathscr{V}^{\textrm{CELL}}_{\textrm{act}}\subseteq\mathscr{V}_{\textrm{act}}$
and
$\mathscr{V}^{\textrm{CELL}}_{\textrm{syn}}\subseteq\mathscr{V}_{\textrm{syn}}$
denote the set of nodes in the actual and synthetic networks lying within a
rectangular grid cell. We use the following metric to compare the percentage
deviation in the node distribution for grid cells where
$\lvert\mathscr{V}_{\textrm{act}^{\textrm{CELL}}}\rvert\neq 0$:

$\mathsf{D}_{\textrm{N}}^{\textrm{CELL}}=\dfrac{\dfrac{\lvert\mathscr{V}^{\textrm{CELL}}_{\textrm{act}}\rvert}{\lvert\mathscr{V}_{\textrm{act}}\rvert}-\dfrac{\lvert\mathscr{V}^{\textrm{CELL}}_{\textrm{syn}}\rvert}{\lvert\mathscr{V}_{\textrm{syn}}\rvert}}{\dfrac{\lvert\mathscr{V}^{\textrm{CELL}}_{\textrm{act}}\rvert}{\lvert\mathscr{V}_{\textrm{act}}\rvert}}\times
100$

Fig. 19 shows the spatial comparison of node distribution for uniform
rectangular cell partitions of two different resolutions followed by shifting
the demarcations horizontally.. The color code denotes the intensity of the
percentage deviation in the distribution. Note that some grid cells are shaded
with black dots to denote the unavailability of network data. It is easy to
note that the spatial distribution of nodes is a function of the size of the
rectangular grid cells as well as the location of the demarcations between
cells. To address this issue, one might consider comparing networks through
multiple grid sizes (varying the resolution) or by considering a single-sized
rectangular grid cells followed by shifting the horizontal and vertical
demarcations of the grid to account for different sized cells.

Geometry comparison We consider a metric space $(\mathscr{M},d)$ where
$d(x,y)$ is the measure between two elements $x,y\in\mathscr{M}$. Given two
non-empty subsets $\mathscr{X},\mathscr{Y}\subseteq\mathscr{M}$ of the metric
space, the Hausdorff distance $d_{\textrm{H}}(\mathscr{X},\mathscr{Y})$
between the sets is defined as

$d_{\textrm{H}}(\mathscr{X},\mathscr{Y})=\sup_{x\in\mathsf{X}}~{}\inf_{y\in\mathsf{Y}}d(x,y)$
(25)

Figure 20: Plots comparing Hausdorff distance between edge geometries of
actual and synthetic networks. The comparison is made for rectangular grid
cells where data regarding the actual network is available. Two different grid
resolutions are considered for the spatial distribution comparison. Also we
compare the Hausdorff distance by moving the horizontal grid demarcations and
observe varied results.

In the context of networks associated with geographic attribute, we define the
metric space $\mathscr{M}$ to be the set of all points along edges in the
network and the metric $d(x,y):=\mathsf{dist}(x,y)$ which is the geodesic
distance between a pair of points $x,y\in\mathscr{M}$. We interpolate points
along each edge of either network in order to evaluate an exact comparison of
the network geometries. Let the set of interpolated points along the actual
and synthetic networks be represented by $\mathscr{P}_{\textrm{act}}$ and
$\mathscr{P}_{\textrm{syn}}$ respectively. Note that
$\mathscr{P}_{\textrm{act}},\mathscr{P}_{\textrm{syn}}\subseteq\mathscr{M}$.
We define the Hausdorff distance between networks $\mathscr{G}_{\textrm{act}}$
and $\mathscr{G}_{\textrm{syn}}$ as

$\mathsf{D}_{\textrm{H}}\left(\mathscr{G}_{\textrm{act}},\mathscr{G}_{\textrm{syn}}\right):=\max_{x\in\mathscr{P}_{\textrm{act}}}\min_{y\in\mathscr{P}_{\textrm{syn}}}\mathsf{dist}(x,y)$
(26)

Fig. 20 shows a Hausdorff distance-based network geometry comparison between
actual and synthetic networks for uniform rectangular grid partitions of two
different resolutions. This is accompanied by shifting the grid cell
demarcations horizontally. We observe that as the cell demarcations are
shifted horizontally, the Hausdorff distance between the edge geometries in
each grid cell alters. Therefore, the grid resolution and grid demarcations
play an important role in the geometry comparison as well.

## Impact of photovoltaic penetration

Installing photovoltaic (PV) generation on residence rooftops is an easy and
viable option to improve the resilience of distribution systems. However, the
continuous increase in PV penetration has led to serious operational issues
such as overvoltage and unbalances in the distribution grid [48]. Therefore,
it is necessary to identify the maximum allowable PV penetration without
violating operational and performance constraints.

It has been observed that a penetration as low as $2.5\%$ is capable of
causing voltage violation when a large PV generation is installed at a single
point in the medium voltage (MV) network. On the contrary, the low voltage
(LV) network can withstand multiple PV generators with penetration levels up
to $110\%$ [49]. In this work, we focus on the impact of PV penetration in the
LV networks through the installation of PV generation on residence rooftops.
We aim to address the following question: does PV penetration provide any
advantage to grid operation, and how much PV penetration is acceptable without
violating operational constraints.

While planning a practical network, the LV feeders are designed with a
suitable rating to supply a certain load to certain groups of residences
without considering the impact of distributed generation (DGs). Most grid
connected DGs needs to confer with the CAN/CSA-C22.2 No.257-06 [50] standard,
which specifies the electrical requirements for the interconnection of
inverter-based micro-distributed resource systems with grid-connected low-
voltage systems. For single phase connection, it is considered normal
operating conditions (NR-Normal Range) when the voltage level is between 0.917
and 1.042 pu. On extreme operation conditions (ER-Emergency Range), the under
and over steady state voltage limits are 0.88 pu and 1.058 pu, respectively.
It is worth mentioning that although networks are allowed to operate under
extreme conditions, the utility would need to take a corrective action.

Figure 21: Comparison of voltage profile without and with PV penetration in
the distribution network. With the introduction of rooftop PVs, the voltage
profile is alleviated to 1 pu, which is recommended by engineering standards.

If a distribution network is operated without any PV, the substation feeder is
solely responsible for delivering power to all residences in the network.
Therefore, the farthest connected residences suffer the maximum voltage drop
and experience a voltage sag. This is denoted in the left figure of Fig. 21
where we observe that a significant percentage of nodes experience
undervoltage (less than 0.9 pu). The introduction of rooftop PVs installed on
randomly chosen $50\%$ of residences alleviates the voltage profile of the
entire grid as shown in the right figure of Fig. 21. In this example, we have
considered a $30\%$ PV penetration , that is the total PV generation in the
distribution network amounts to $30\%$ of the total load. Therefore, the
substation feeders are responsible for delivering the remaining $70\%$ of the
power.

We perform the comparison on two different synthetic networks. Fig. 22 shows
the histogram of node voltages for comparing the impact of PV penetration in
an urban feeder (top three figures) and a rural feeder (bottom three figures).
The urban feeder network consists of shorter length lines. We observe that LV
level PV generation is less likely to cause overvoltage issues as compared to
a single node MV level PV integration. In the case of a rural feeder with
longer lines, we observe a similar trend. However, the percentage of nodes
experiencing severe overvoltage (around 1.05 pu which is the extreme limit of
an acceptable overvoltage) is higher for rural networks as compared to urban
feeder networks for all penetration levels when we consider a single node PV
integration at the MV level. Therefore an optimal placement of PV generators
is required for the rural feeders so that they do not suffer from overvoltage
issues.

Figure 22: Comparing the impact of PV integration in an urban feeder (top
figures) and rural feeder (bottom figures) for three different penetration
levels. Colors depict PV addition at multiple locations in LV network (blue)
and single node PV addition in MV network (red). The error bar shows the
variation in impact when analyzed over an ensemble of $20$ networks. LV level
penetration is less likely to cause overvoltage as compared to similar
penetration in the MV network. Rural feeders are more prone to overvoltage
issues.

## Short Circuit Analysis

We can use the created networks to perform short circuit analysis. Since we
have created a positive sequence network for the distribution system, we can
carry out a short circuit analysis for three phase symmetrical faults. To this
end, we construct the bus impedance matrix $\mathbf{Z}_{\textrm{bus}}$ using
the traditional methodology [51].Thereafter, we compute the post-fault voltage
at bus $j$ for a three phase symmetrical fault at bus $k$ using (27).

$V^{\textrm{post-fault}}_{j}=V^{\textrm{pre-
fault}}_{j}-\frac{Z_{jk}}{Z_{kk}}V^{\textrm{pre-fault}}_{k}$ (27)

Note that the pre-fault voltages are computed by solving the power flow
problem using the LDF model [44] with each residence consuming average hourly
demand. Fig. 23 shows the percentage of residences experiencing undervoltage
instantaneously after the occurrence of a symmetrical three phase fault at
different locations in the network. We perform the short circuit analysis for
two different feeders in a single synthetic network. We simulate symmetrical
three phase fault at various locations along each feeder and observe the post-
fault voltages at each residence nodes.

Though we have created positive sequence networks, we are able to perform
short circuit analysis for unsymmetrical faults (line-ground and line-line).
This can be performed under the assumption that all lines in the distribution
network are overhead lines. Under such an assumption, the zero sequence
impedance is three times the positive sequence impedance of the lines.
Thereafter, the fault currents can be computed accordingly [52, 53, 51].

Figure 23: Plots showing short circuit analysis on two different feeders in
one of the created synthetic networks. Three phase symmetrical faults are
simulated at multiple locations along each feeder and post-fault voltages at
residences are computed.

## References

  * [1] A. Zamani, T. Sidhu, and A. Yazdani.  A strategy for protection coordination in radial distribution networks with distributed generators.  In IEEE PES General Meeting, pages 1–8, July 2010.
  * [2] R. M. DŚouza, C. D. Brummitt, and E. A. Leicht.  Modeling Interdependent Networks as Random Graphs: Connectivity and Systemic Risk, pages 73–94.  Springer International Publishing, 2014.
  * [3] S. Soltan and G. Zussman.  Generation of synthetic spatially embedded power grid networks.  In 2016 IEEE Power and Energy Society General Meeting (PESGM), pages 1–5, 2016.
  * [4] C. D. Brummitt, P. D. H. Hines, I. Dobson, C. Moore, and R. M. DŚouza.  Transdisciplinary electric power grid science.  Proceedings of the National Academy of Sciences, 110(30):12159–12159, 2013.
  * [5] K. M. Gegner, A. B. Birchfield, Ti Xu, K. S. Shetye, and T. J. Overbye.  A methodology for the creation of geographically realistic synthetic power flow models.  In 2016 IEEE Power and Energy Conference at Illinois (PECI), pages 1–6, Feb 2016.
  * [6] A. B. Birchfield, K. M. Gegner, T. Xu, K. S. Shetye, and T. J. Overbye.  Statistical considerations in the creation of realistic synthetic power grids for geomagnetic disturbance studies.  IEEE Transactions on Power Systems, 32(2):1502–1510, Mar 2017.
  * [7] R. Atat, M. Ismail, M. F. Shaaban, E. Serpedin, and T. Overbye.  Stochastic geometry-based model for dynamic allocation of metering equipment in spatio-temporal expanding power grids.  IEEE Transactions on Smart Grid (Early Access), 11(3):1–12, Oct 2019.
  * [8] H. Li, J.L. Wert, A. Birchfield, T.J. Overbye, T.G.S. Roman, C. M. Domingo, F.E.P. Marcos, P. D. Martinez, T. Elgindy, and B. Palmintier.  Building highly detailed synthetic electric grid data sets for combined transmission and distribution systems.  IEEE Open Access Journal of Power and Energy, 7:478–488, 2020.
  * [9] IEEE test feeders, 2014.  Last accessed 20 Oct 2021.
  * [10] EPRI test circuits, 2019.  Last accessed 20 Oct 2021.
  * [11] J.E.T. Bistline and G.J. Blanford.  The role of the power sector in net-zero energy systems.  Energy and Climate Change, 2, Dec 2021.
  * [12] N.A. Sepulveda, J.D. Jenkins, A. Edington, D.S. Mallapragada, and R.K. Lester.  The design space for long-duration energy storage in decarbonized power systems.  Nature Energy, 6(5):506–516, May 2021.
  * [13] S. Joshi, S. Mittal, P. Holloway, P. R. Shukla, B. Ó Gallachóir, and J. Glynn.  High resolution global spatiotemporal assessment of rooftop solar photovoltaics potential for renewable electricity generation.  Nature Communications, 12(1):5738, Oct 2021.
  * [14] P. S. Georgilakis and N. D. Hatziargyriou.  A review of power distribution planning in the modern power systems era: Models, methods and future research.  Electric Power Systems Research, 121:89–100, 2015.
  * [15] M. Resener, S. Haffner, L. A. Pereira, and P. M. Pardalos.  Optimization techniques applied to planning of electric power distribution systems: a bibliographic survey.  Energy Systems, 9(3):473–509, Aug 2018.
  * [16] C. Mateo, F. Postigo, F. de Cuadra, T. G. S. Roman, T. Elgindy, Pablo Dueñas, B. Hodge, V. Krishnan, and B. Palmintier.  Building large-scale u.s. synthetic electric distribution system models.  IEEE Transactions on Smart Grid, 11(6):5301–5313, 2020.
  * [17] E. Schweitzer, A. Scaglione, A. Monti, and G. A. Pagani.  Automated generation algorithm for synthetic medium voltage radial distribution systems.  IEEE Journal on Emerging and Selected Topics in Circuits and Systems, 7(2):271–284, Jun 2017.
  * [18] E.A. Bender and E.R. Canfield.  The asymptotic number of labeled graphs with given degree sequences.  Journal of Combinatorial Theory, Series A, 24(3):296–307, 1978\.
  * [19] M. Liang, Y. Meng, J. Wang, D.L. Lubkeman, and N. Lu.  Feedergan: Synthetic feeder generation via deep graph adversarial nets.  IEEE Transactions on Smart Grid, 12(2):1163–1173, 2021.
  * [20] C.M. Domingo, T.G.S. Roman, A. Sanchez-Miralles, J. P. Gonzalez, and A.C. Martinez.  A reference network model for large-scale distribution planning with automatic street map generation.  IEEE Transactions on Power Systems, 26(1):190–197, Feb 2011.
  * [21] L. Gonzalez-Sotres, C.M. Domingo, A. Sanchez-Miralles, and M.A. Miro.  Large-scale MV/LV transformer substation planning considering network costs and flexible area decomposition.  IEEE Transactions on Power Delivery, 28(4):2245–2253, Oct 2013\.
  * [22] S. S. Saha, E. Schweitzer, A. Scaglione, and N.G. Johnson.  A framework for generating synthetic distribution feeders using openstreetmap.  In 2019 North American Power Symposium (NAPS), pages 1–6, 2019\.
  * [23] R. Kadavil, T. M. Hansen, and S. Suryanarayanan.  An algorithmic approach for creating diverse stochastic feeder datasets for power systems co-simulations.  In 2016 IEEE PES General Meeting (PESGM), pages 1–5, Jul 2016.
  * [24] A. Bidel, T. Schelo, and T. Hamacher.  Synthetic distribution grid generation based on high resolution spatial data.  In 2021 IEEE International Conference on Environment and Electrical Engineering and 2021 IEEE Industrial and Commercial Power Systems Europe (EEEIC / ICPS Europe), pages 1–6, 2021.
  * [25] D. Quiroga, E. Sauma, and D. Pozo.  Power system expansion planning under global and local emission mitigation policies.  Applied Energy, 239:1250–1264, 2019.
  * [26] S. You, S. W. Hadley, M. Shankar, and Y. Liu.  Co-optimizing generation and transmission expansion with wind power in large-scale power grids—implementation in the us eastern interconnection.  Electric Power Systems Research, 133:209–218, 2016.
  * [27] A. Trpovski, D. Recalde, and T. Hamacher.  Synthetic distribution grid generation using power system planning: Case study of Singapore.  In 2018 53rd International Universities Power Engineering Conference (UPEC), pages 1–6, Sep 2018.
  * [28] G. Byeon, P.V. Hentenryck, R. Bent, and H. Nagarajan.  Communication-constrained expansion planning for resilient distribution systems.  INFORMS Journal on Computing, 32(4):968–985, Apr 2020.
  * [29] M. K. Singh, V. Kekatos, and C.C. Liu.  Optimal distribution system restoration with microgrids and distributed generators.  In 2019 IEEE PES General Meeting, 2019.
  * [30] M. K. Singh, V. Kekatos, S. Taheri, K. P. Schneider, and C.C. Liu.  Enforcing radiality constraints for der-aided power distribution grid reconfiguration, 2019.
  * [31] S. Lei, C. Chen, Y. Song, and Y. Hou.  Radiality constraints for resilient reconfiguration of distribution systems: Formulation and application to microgrid formation, 2019.
  * [32] GridLAB-D feeders, 2015.  Last accessed 20 Oct 2021.
  * [33] Z. Wang, A. Scaglione, and R. J. Thomas.  Generating statistically correct random topologies for testing smart grid communication and control networks.  IEEE Transactions on Smart Grid, 1(1):28–39, 2010.
  * [34] V. Krishnan, B. Bugbee, T. Elgindy, C. Mateo, P. Duenas, F. Postigo, J. S. Lacroix, T. G. S. Roman, and B. Palmintier.  Validation of synthetic u.s. electric power distribution system data sets.  IEEE Transactions on Smart Grid, 11(5):4477–4489, Sept 2020.
  * [35] American National Standards Institute ANSI.  ANSI C84.1-2020: Electric Power Systems Voltage Ratings (60 Hz), 2020.
  * [36] U.S. Department of Homeland Security.  Electric substations.  https://hifld-geoplatform.opendata.arcgis.com/datasets/electric-substations, 2019\.
  * [37] S. Thorve, S. Swarup, A. Marathe, Y. Chungbaek, E. K. Nordberg, and M. V. Marathe.  Simulating residential energy demand in urban and rural areas.  In 2018 Winter Simulation Conference (WSC), pages 548–559, Dec 2018\.
  * [38] N. Rotering, C. Schröders, J. Kellermann, and A. Moser.  Medium-voltage network planning with optimized power factor control of distributed generators.  In 2011 IEEE Power and Energy Society General Meeting, pages 1–8, Jul 2011.
  * [39] Open Street Maps, 2021.  Last accessed 20 Apr 2021.
  * [40] Raphael A. F. and J. L. Bentley.  Quad trees a data structure for retrieval on composite keys.  Acta Informatica, 4:1–9, 1974.
  * [41] R. Meyur, M. Marathe, A. Vullikanti, H. Mortveit, V. Centeno, and A. Phadke.  Creating realistic power distribution networks using interdependent road infrastructure.  In IEEE Big Data Conference, Dec 2020.
  * [42] F. P. Preparata and M. I. Shamos.  Computational Geometry: An Introduction.  Springer-Verlag, Berlin, Heidelberg, 1985.
  * [43] F. Postigo, C. Mateo, T. Gómez, B. Palmintier, B.M. Hodge, V. Krishnan, F. de Cuadra, and B. Mather.  A review of power distribution test feeders in the united states and the need for synthetic representative networks.  Energies, 10(11):1896.1–14, Nov 2017.
  * [44] S. Bolognani and F. Dörfler.  Fast power system analysis via implicit linearization of the power flow manifold.  In 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 402–409, Sep 2015.
  * [45] G. P. McCormick.  Computability of global solutions to factorable nonconvex programs: Part i – convex underestimating problems.  Mathematical Programming, 10(1):147–175, Dec 1976.
  * [46] C. Mateo, F. Postigo, F. De Cuadra, T. Gómez, T. Elgindy, P. Dueñas, B. Hodge, V. Krishnan, and B. Palmintier.  RNM-US Catalog of Equipment.  https://dx.doi.org/10.21227/4vne-pd37, 2020.
  * [47] R. Meyur.  Synthetic distribution.  https://github.com/rounak-meyur/synthetic-distribution/tree/master/output, 2021\.  Last accessed 5 Apr 2022.
  * [48] S. Fatima, V. Püvi, and M. Lehtonen.  Review on the pv hosting capacity in distribution networks.  Energies, 13(18), 2020.
  * [49] T. Aziz and N. Ketjoy.  Pv penetration limits in low voltage networks and voltage variations.  IEEE Access, 5:16784–16792, 2017.
  * [50] CAN/CSA-C22.2 No.257-06. interconnecting inverter-based micro-distributed resources to distribution systems, 2006.
  * [51] J. J. Grainger and W. D. Stevenson.  Power System Analysis.  McGraw-Hill, 1994.
  * [52] S. H. Horowitz and A. G. Phadke.  Power System Relaying.  Research Studies Press, 2nd Edition, Taunton, UK, 1995.
  * [53] C. L. Wadhwa.  Electrical Power Systems.  New Age International, 2009.