# Chance-Constrained Optimal Distribution Network Partitioning to Enhance Grid
Resilience

Shuchismita Biswas, Manish K. Singh, and Virgilio A. Centeno The authors are
with Power and Energy Center (PEC), Department of Electrical and Computer
Engineering, Virginia Tech, Blacksburg, VA, U.S.A. Emails :
<EMAIL_ADDRESS>
###### Abstract

This paper formulates a chance-constrained optimal distribution network
partitioning (ODNP) problem addressing uncertainties in load and renewable
energy generation; and presents a solution methodology using sample average
approximation (SAA). The objective is to identify potential sub-networks in
the existing distribution grid; that are likely to survive as self-adequate
islands if supply from the main grid is lost. This constitutes a planning
problem. Practical constraints like ensuring network radiality and
availability of grid-forming generators are considered. Quality of the
obtained solution is evaluated by comparison with- a) an upper bound on the
probability that the identified islands are supply-deficient, and b) a lower
bound on the optimal value of the true problem. Performance of the ODNP
formulation is illustrated on a modified IEEE 37-bus feeder. It is shown that
the network flexibility is well utilized; the partitioning changes with risk
budget; and that the SAA method is able to yield good quality solutions with
modest computation cost.

###### Index Terms:

microgrids, renewable energy generation, radiality, chance-constrained
optimization, resilience, DER

##  I Introduction

In recent years, the adoption of renewable energy based distributed energy
resources (DERs) has increased due to the recognition of their economic and
environmental benefits. A primary advantage of DERs is their ability to
sustain local loads if the main grid is lost, possibly due to natural
disasters. DERs and loads may be clustered together to form _microgrids_ , a
resiliency resource, which supply essential loads and aid service restoration
during and after outages [1]. According to the IEEE 1547.4-2011 standard,
microgrids: _1)_ have DERs and load; _2)_ can operate in both grid-connected
and islanded modes; and _3)_ are intentionally planned[2]. For safe
operations, microgrids must have adequate control capabilities.

Utilities are interested in identifying parts of the existing distribution
network that can be converted to microgrids via economically viable
retrofitting. This is because DERs cannot supply local loads during an outage
if adequate control and protection schemes are not in place. Hence, optimally
splitting a network into microgrids constitutes a pertinent planning problem
[3, 4, 5, 6, 7, 8, 9, 10, 11]. This optimal distribution network partitioning
(ODNP) task seeks to identify potential _self-adequate_ sub-networks that can
survive the loss of the main grid as islands. Both exact [3, 5, 4, 7, 6] and
heuristic [9, 10, 8, 11] methods have been proposed for ODNP. Self-adequacy in
the objective function has been surrogated by either expected power flow on
microgrid boundary lines [10, 9, 7, 11], or expected load-generation imbalance
within microgrids [3, 5, 4, 7, 6, 8]. Moreover, dynamic identification of
boundary lines in response to faults have also been proposed [6, 7]. A method
for determining self-sufficient islands in transmission networks is described
in [12], but cannot be directly extended to ODNP without including
distribution system specific constraints.

Distribution networks are usually operated radially for protection
coordination, and this radiality needs to be maintained while separating into
microgrids. In [9, 5, 10, 7, 11, 6], ODNP is demonstrated on an already radial
feeder and radiality is not explicitly enforced. This approach ignores
normally open switches, and under-utilizes network flexibility. Radiality is
considered in [4], but another restrictive condition is imposed- each
microgrid is assigned exactly one DER. This single DER constraint is also
present in [6]. In this approach, the number of partitions are predetermined,
leading to sub-optimal solutions. A radiality constraint without specifying
the number of microgrids was recently presented in [3], and the formulation in
the current work builds upon this approach.

A critical aspect that has been overlooked in the existing microgrid planning
literature is the requirement of grid-forming generators in viable islands.
The 1547 standard mandates that an island should have at least one generator
that provides voltage and frequency support during a system disturbance, or
has black-start capabilities [2]. An exhaustive path search based method for
checking connectivity to black-start generators has been proposed in [13].
Another multiple commodity flow based approach outlined in [14] separately
checks nodes for their connectivity to black-start nodes. Both these
approaches become computationally prohibitive for large networks. The ODNP
formulation put forth in the present work guarantees that all nodes in each
microgrid will be connected to at least one grid-forming generator. The
formulation is somewhat similar to the single commodity flow model of [14] but
uses fewer constraints and shows faster performance (empirically observed to
be 10 to 20% faster).

ODNP is further complicated by the uncertainty in demands and available
generation capacity. In [9, 5, 10, 11, 7], the load-generation uncertainty is
addressed by constructing typical daily profiles, over which optimization is
performed. However, the quality of solution obtained is not evaluated. The
present work formulates a chance-constrained ODNP (cc-ODNP) to identify
optimal microgrids in the planning stage. This is computationally challenging
as the underlying deterministic formulation uses mixed integer linear program
(MILP) and is not convex by nature. Hence, sampling and integer programming
[15] has been used to solve an approximation of the cc-ODNP, and the quality
of the solution obtained is rigorously evaluated.

The main contributions of this work are as follows. _First_ , a deterministic
ODNP for identifying optimal microgrids, given real-time load-generation
values is formulated in Section III. Practical constraints are comprehensively
addressed, with the following novel aspects: _i)_ efficiently formulating a
generalized radiality condition, and _ii)_ ensuring every microgrid includes a
grid-forming generator, without any pre-assignment. Thus, optimality is not
compromised. _Second,_ a probabilistic ODNP problem is formulated, and solved
using a computationally tractable sample average approximation (SAA) based
MILP. While the SAA approach offers asymptotic equivalence to the original
probabilistic formulation, in practice, computational resources restrict the
number of scenarios that can be analyzed. Hence, rigorous stochastic tools
have been used in section IV to statistically assess the quality of an
obtained solution, in terms of confidence in feasibility and relevance of the
attained objective. _Third,_ in section V, performance of ODNP is demonstrated
through extensive numerical tests on a modified IEEE 37-bus feeder. It is
shown that the SAA approach is able to efficiently utilize network
flexibility, and outperforms a robust clustering based methodology in terms of
objective cost.

##  II Preliminaries

In this section, some mathematical preliminaries are revisited before
expounding on the problem formulation. Calligraphic symbols represent sets,
lower case bold letters represent column vectors, and upper case bold letters
denote matrices. All zero and all one vectors and matrices of appropriate size
are denoted by $\mathbf{0}$ and $\mathbf{1}$ respectively.

###  II-A Graph Theory

A graph $\mathcal{G}:=(\mathcal{V},\mathcal{E})$ consists of a _vertex set_
$\mathcal{V}$ and an _edge set_ $\mathcal{E}$, where an edge is an unordered
pair of distinct vertices of $\mathcal{G}$. Edge $e_{ij}\in\mathcal{E}$ is
denoted by its incident vertices $(i,j)$, such that $i,j\in\mathcal{V}$. If
$e_{ij}\in\mathcal{E}$, then vertices $i$ and $j$ are adjacent. Two edges are
adjacent if they have a common vertex. A subgraph of $\mathcal{G}$ is a graph
$\mathcal{H}:=(\mathcal{X,Y})$ such that $\mathcal{X}\subseteq\mathcal{V}$ and
$\mathcal{Y}\subseteq\mathcal{E}$. If $\mathcal{X}=\mathcal{V}$, then
$\mathcal{H}$ is a spanning subgraph of $\mathcal{G}$. $\mathcal{H}$ is an
induced subgraph of $\mathcal{G}$ if vertices in $\mathcal{X}$ are adjacent in
$\mathcal{H}$ if and only if they are adjacent in $\mathcal{G}$.

A path from $i$ to $j$ is a sequence of distinct vertices starting at $i$ and
ending at $j$ such that consecutive vertices are adjacent. If there is a path
between all pairs of vertices of a graph $\mathcal{G}$, then $\mathcal{G}$ is
connected; else $\mathcal{G}$ is disconnected. An induced subgraph of
$\mathcal{G}$ that is maximal, subject to being connected, is called a
connected component of $\mathcal{G}$. A cycle is a sequence of adjacent edges
without repetition that starts and ends at the same node. A graph with no
cycles is called acyclic. A connected and acyclic graph is a tree. A spanning
tree subgraph of $\mathcal{G}$ is a tree that covers all vertices in
$\mathcal{G}$. An acyclic graph with multiple connected components is a
forest. A spanning forest subgraph of $\mathcal{G}$ is a forest that covers
all vertices in $\mathcal{G}$. Spanning forests may include connected
components with a single node. A review of graph theory fundamentals is
available in [16].

###  II-B Chance-Constrained Optimization

Stochastic optimization refers to a collection of methods for solving an
optimization problem with uncertain parameters. For many real-world
applications operating in uncertain environments, ensuring 100% reliability is
physically and economically impractical. This difficulty is often dealt with
by designing systems that assure a minimum reliability level with high
probability. Mathematical models of such reliability-constrained systems
involve the use of probabilistic or chance constraints [15]. A generic chance-
constrained optimization (CCO) problem is of the form

$\displaystyle\min_{\boldsymbol{x}\in\mathcal{X}}\quad$ $\displaystyle
f(\boldsymbol{x})$ ($P_{1}$) s. to
$\displaystyle\boldsymbol{h}(\boldsymbol{x})\leq\boldsymbol{0}$ ($C_{1}$)
$\displaystyle
Pr\\{\boldsymbol{g}(\boldsymbol{x},\boldsymbol{\xi})\leq\boldsymbol{0}\\}\geq
1-\varepsilon$ ($C_{2}$)

Here, $\boldsymbol{x}$ is the vector of decision variables, whose feasible
region is given by $\mathcal{X}\subset\mathds{R}^{n}$. The objective function
to be minimized is $f:\mathds{R}^{n}\rightarrow\mathds{R}$. Vector
$\boldsymbol{\xi}$ stacks the uncertain parameters with known probability
distribution, and $\varepsilon\in(0,1)$ is a tunable risk parameter. Problem
$P_{1}$ seeks to find an optimal decision vector $\boldsymbol{x}^{*}$ that
minimizes $f(\boldsymbol{x})$, such that the hard constraints $C_{1}$ are
always satisfied, while the chance constraint $C_{2}$ is satisfied with
probability at least $1-\varepsilon$.

In power systems literature, CCO has been previously used to address security
constrained economic dispatch and unit commitment problems [17]. This class of
problems is difficult to solve, due to two main reasons:

• Given a candidate solution $\boldsymbol{\bar{x}}\in\mathcal{X}$, accurately
computing
$Pr\\{\boldsymbol{g}(\boldsymbol{\bar{x}},\boldsymbol{\xi})\leq\mathbf{0}\\}$
can be very difficult, making it hard to check if constraint $C_{2}$ is
satisfied.

• The feasibility region defined by a chance constraint is usually not convex
[15]. This makes finding an optimal solution difficult even when the
feasibility of $\boldsymbol{\bar{x}}$ can be checked.

These difficulties may be overcome by considering a sample average
approximation (SAA) of the original problem where the true distribution of
$\boldsymbol{\xi}$ is replaced by an empirical distribution with discrete
support. The SAA is still a chance-constrained stochastic problem, but with a
different distribution for $\boldsymbol{\xi}$, and may be solved via integer
programming [15]. This method has been shown to yield good candidate solutions
if the sampling is ample and rich. In this work, the SAA approach will be
incorporated to solve a probabilistic ODNP and the solution obtained will be
further analyzed to verify how well it solves the original chance constrained
problem.

##  III Problem Formulation

Given a distribution network with DERs, planners would like to optimally
construct microgrids, such that DERs sustain internal loads if supply from the
main grid is lost. Load served is to be maximized. Both load and generation
vary with weather and assuring self-adequacy for the worst case may lead to
very conservative solutions. Hence, a solution that works well for most
operating conditions might be preferred. Thus microgrids may be designed to be
self-adequate with probability at least $(1-\varepsilon)$ across all possible
operating scenarios, where $\varepsilon$ is a tunable risk parameter. The
value of $\varepsilon$ may be chosen based on available storage resources.
Once optimal microgrids are identified, they need to be equipped with control
capabilities and boundary line switches. It must be noted that depending on
the generation capacity of installed DERs, all load may not be served by
microgrids.

Our mathematical formulation is put forth in three steps. First, a
deterministic version of the problem, d-ODNP is presented where the load
served is maximized for a given scenario of demands and generation. Next, the
chance constraints arising from the randomness in generation and demands are
added. DERs are assumed to be dispatchable subject to stochastic generation
capacity. Such units in practice could be photo-voltaic (PV) generators,
diesel generators (DGs) or combined heat and power plants (CHPs) that are
plausible in a low/medium voltage network setup. Non-dispatchable generators
may be seamlessly incorporated in the formulation as negative stochastic
demands. Finally, a SAA based algorithm is proposed that can tractably solve
the probabilistic ODNP.

###  III-A Distribution Network Model

A single-phase distribution network may be represented by a connected directed
graph $\mathcal{G}_{N}:=(\mathcal{V}_{N},\mathcal{E}_{N})$, where vertices
denote buses and edges denote lines. The substation node is indexed by $0$;
and the set of all other nodes is denoted by
$\mathcal{V}:=\mathcal{V}_{N}\setminus\\{0\\}$. Each edge
$e_{i,j}\in\mathcal{E}_{N}$ is assigned an arbitrary direction from node $i$
to $j$. If $e_{i,j}\in\mathcal{E}_{N}$, then $e_{j,i}\notin\mathcal{E}_{N}$.
The task at hand considers that the main grid is unavailable, hence
partitioning needs to be carried out on
$\mathcal{G}:=(\mathcal{V},\mathcal{E})$, the induced subgraph of
$\mathcal{G}_{N}$ on vertex set $\mathcal{V}$. In the present setup all lines
are considered switchable. Any non-switchable edge coinciding with a microgrid
boundary would need to be retrofitted with a switch. Moreover, edges include
lines with existing normally open and normally closed switches, and hence
$\mathcal{G}$ is not necessarily radial.

Each node has an associated demand $(\xi^{d_{p}}_{i}+j\xi^{d_{q}}_{i})$ and
generation capacity $(\xi^{g_{p}}_{i}+j\xi^{g_{q}}_{i})$. The demand and
generation capacities are not precisely known at the planning stage and only a
probability distribution, possibly empirical, may be available. Let $v_{i}$ be
the voltage magnitude at bus $i$ and $(p_{i}+jq_{i})$ be the complex power
injection. Bus voltages, demand, generation capacity, and complex power
injections are respectively stacked into vectors $\mathbf{v}$,
$\boldsymbol{\xi^{d_{p}}}+j\boldsymbol{\xi^{d_{p}}},\boldsymbol{\xi^{g_{p}}}+j\boldsymbol{\xi^{g_{q}}},\mathbf{p}+j\mathbf{q}$.
All quantities are in per units.

Let us introduce two binary decision variables $b^{n}_{i}\in\\{0,1\\}$ and
$b^{e}_{ij}\in\\{0,1\\}$ that respectively dictate if vertex $i\in\mathcal{V}$
and edge $e_{ij}\in\mathcal{E}$ are energized. If $e_{ij}$ is energized, then
adjacent vertices $i$ and $j$ need to be energized as well. Mathematically,

$b^{n}_{i}+b^{n}_{j}\geq 2b^{e}_{ij}\quad\forall{e_{ij}\in\mathcal{E}}$ (1)

For safe operations, ANSI standards mandate that voltages at active buses
should be within $\pm 5\%$ p.u. of the nominal value [18]. Mathematically,

$0.95~{}\mathbf{b^{n}}\leq\mathbf{v}\leq 1.05~{}\mathbf{b^{n}}.$ (2)

Let the power flow on line $e_{ij}\in\mathcal{E}$ be $P_{ij}+jQ_{ij}$. The
line capacity constraints may be formulated as follows.

$\displaystyle-b^{e}_{ij}P^{max}_{ij}\leq P_{ij}\leq
b^{e}_{ij}P^{max}_{ij}\quad\forall{e_{ij}\in\mathcal{E}}$ (3a)
$\displaystyle-b^{e}_{ij}Q^{max}_{ij}\leq Q_{ij}\leq
b^{e}_{ij}Q^{max}_{ij}\quad\forall{e_{ij}\in\mathcal{E}}$ (3b)

Flow constraints of the form $P_{ij}^{2}+Q_{ij}^{2}\leq S_{ij}^{2}$ are not
used here to avoid quadratic constraints. A polytopic approximation of this
constraint proposed in [19] could also be used.

###  III-B Power Flow Model

Power injection at energized buses is assumed to be controllable subject to
the maximum capacity, while loads are assumed to be inelastic. Some generators
can absorb reactive power, and this flexibility is included in the formulation
as negative reactive power generation. Let the maximum reactive power
absorption at bus $i$ be denoted as $q_{i}^{min}$. The value of $q_{i}^{min}$
can be set to zero to indicate the absence of reactive power absorption
capabilities. The nodal power injections are governed by the following
expressions:

$\displaystyle 0\leq p^{g}_{i}\leq
b^{n}_{i}\xi^{g_{p}}_{i}\quad\forall{i\in\mathcal{V}}$ (4a)
$\displaystyle-b^{n}_{i}q_{i}^{min}\leq q^{g}_{i}\leq
b^{n}_{i}\xi^{g_{q}}_{i}\quad\forall{i\in\mathcal{V}}$ (4b)

Here, constraints (4a)-(4b) establish that generation output
$p^{g}_{i},q^{g}_{i}$ at bus $i$ is bounded. The maximum capacity of DGs may
be fixed based on the machine rating. The active power capacity of PV
generators is limited by solar irradiance levels. Similarly, the maximum
active power generation by CHPs is affected by the local heating demand. These
factors make $\xi^{g_{p}}$ stochastic, in general. Furthermore, given the
limit on active power, a corresponding stochastic limit on reactive power may
be obtained based on the operating curves of the respective generators. For
inelastic loads, the constraints for power consumption may be formulated as:

$\displaystyle
p^{d}_{i}=b^{n}_{i}\xi^{d_{p}}_{i}\quad\forall{i\in\mathcal{V}}$ (5a)
$\displaystyle
q^{d}_{i}=b^{n}_{i}\xi^{d_{q}}_{i}\quad\forall{i\in\mathcal{V}}$ (5b)

Constraints (5a)-(5b) state that consumption at bus $i$ is equal to its
demand; if energized. Thus the net power injections are:

$\displaystyle p_{i}=p^{g}_{i}-p^{d}_{i}\quad\forall{i\in\mathcal{V}}$ (6a)
$\displaystyle q_{i}=q^{g}_{i}-q^{d}_{i}\quad\forall{i\in\mathcal{V}}$ (6b)

For the power flow, the linearized distribution flow (LDF) model proposed in
[20] is followed. Despite being an approximation for the full AC power flow
model, LDF has been used extensively and shown to perform well in literature
[21]. Thus, ignoring line losses, the power balance at each node entails:

$\displaystyle\sum_{e_{ij}\in\mathcal{E}}P_{ij}-\sum_{e_{jk}\in\mathcal{E}}P_{jk}=p_{j}\quad\forall{j\in\mathcal{V}}$
(7a)
$\displaystyle\sum_{e_{ij}\in\mathcal{E}}Q_{ij}-\sum_{e_{jk}\in\mathcal{E}}Q_{jk}=q_{j}\quad\forall{j\in\mathcal{V}}$
(7b)

Let $r_{ij}+jx_{ij}$ be the impedance of line $e_{ij}\in\mathcal{E}$. Then,
the relationship between voltages and power injections may be linearized as:
$v^{2}_{i}-v^{2}_{j}=2(r_{ij}P_{ij}+x_{ij}Q{ij})$. Assuming small voltage
deviations, the squared terms may be approximated as $v^{2}_{i}\simeq
2v_{i}-1$. Using these results,

$\displaystyle b^{e}_{ij}(v_{i}-v_{j}-r_{ij}P_{ij}-x_{ij}Q_{ij})=0,~{}\forall
e_{ij}\in\mathcal{E}$ (8)

Here, the indicator $b^{e}_{ij}$ is multiplied to enforce the voltage drop
relation only for the energized lines. Bilinear terms like $b^{e}_{ij}v_{i}$
in (8) can be handled by McCormick linearization, where the product terms are
replaced by their linear convex envelopes to yield a relaxation of the
original non-convex feasible set [22]. If at most one of the variables is
continuous and the rest are binary, this relaxation is exact. For
illustration, let us consider a term $z=xy$, where $x$ is binary and $y$ is a
continuous variable bounded in $y\in[\underline{y},\overline{y}]$. Here,
$z=xy$ may be equivalently expressed as four linear inequality constraints.

$\displaystyle x\underline{y}$ $\displaystyle\leq z\leq x\overline{y}$ (9a)
$\displaystyle y+(x-1)\overline{y}$ $\displaystyle\leq z\leq
y+(x-1)\underline{y}$ (9b)

Note that putting $x=0$ in (9a)-(9b) yields $z=0$. Similarly, putting $x=1$
yields $z=y$. All such bilinear terms appearing henceforth in this work will
be treated in a similar manner.

###  III-C Radiality Constraint

Network radiality is essential for distribution system operations. Some
approaches proposed for enforcing radiality are cycle elimination [13] and
virtual commodity flow [3, 23]. The ODNP task needs to identify
$\mathcal{G}^{\prime}=(\mathcal{V},\mathcal{E}^{\prime})$, a spanning subgraph
of $\mathcal{G}$, such that every connected component, or simply component, of
$\mathcal{G}^{\prime}$ is a tree, i.e. $\mathcal{G}^{\prime}$ is a forest. A
spanning forest may include isolated nodes, i.e. it may have components with a
single node. The radiality constraints in this context will be formulated
using the condition stated next.

Figure 1: Components $(\\{A,B,C\\},\\{AB,BC\\})$,$(\\{D,E,F\\},\\{DE,DF\\})$
and $(\\{G\\},\\{\\})$ form a spanning forest. Adding edges $BD$ and $DG$
creates a spanning tree.

Proposition 1.([3]) Given a spanning forest subgraph
$\mathcal{F}:=(\mathcal{V},\mathcal{E}_{F})$ of a connected graph
$\mathcal{G}:=(\mathcal{V},\mathcal{E})$, there exists at least one spanning
tree subgraph of $\mathcal{G}$, expressed as
$\mathcal{T}:=(\mathcal{V},\mathcal{E}_{T})$, such that
$\mathcal{E}_{F}\subseteq\mathcal{E}_{T}\subseteq\mathcal{E}$.

In other words, some edges may be removed from a spanning tree to obtain a
spanning forest. This idea is illustrated in fig. 1. The solid lines show
edges in a spanning forest and the addition of dashed edges creates a spanning
tree. Hence, the radiality of
$\mathcal{G}^{\prime}:=(\mathcal{V},\mathcal{E}^{\prime})$ holds true if there
exists a fictitious spanning tree subgraph
$\mathcal{T}:=(\mathcal{V},\mathcal{E}_{T})$ of $\mathcal{G}$ such that
$\mathcal{E}^{\prime}\subseteq\mathcal{E}_{T}$. Let us first establish a
condition to select a spanning tree, and then extract the required spanning
forest from it. The base topology of $\mathcal{G}$ may be captured by a
branch-bus incidence matrix $\mathbf{\tilde{A}}$ of dimension
$|\mathcal{E}|\times(|\mathcal{V}|)$, with the following entries:

$\mathbf{\tilde{A}}_{e_{ij},k}:=\begin{cases}1\quad&,\quad k=i\\\
-1\quad&,\quad k=j\\\ 0&,\quad\text{otherwise}\end{cases}$ (10)

The first column $\mathbf{a}_{1}$ of $\mathbf{\tilde{A}}$ may be separated as
$\mathbf{\tilde{A}}=[\mathbf{a}_{1}\hskip 3.0pt\mathbf{A}]$. This yields the
reduced branch-bus incidence matrix $\mathbf{A}$ of $\mathcal{G}$. An
efficient model for imposing graph connectivity put forth in [23] posits that
a graph with vertex set $\mathcal{V}$ and reduced branch-bus incidence matrix
$\mathbf{A}$, is connected if and only if there exists a vector
$\mathbf{f}\in\mathds{R}^{|\mathcal{E}|}$, such that
$\mathbf{A}^{T}\mathbf{f}=\mathbf{1}$. For proof, see [23]. For a physical
interpretation, consider every vertex in $\mathcal{V}\setminus\\{1\\}$ injects
unit virtual commodity into the network represented by the graph. Then,
$\mathbf{f}$ denotes the flow of commodities on the edges. If this flow setup
is feasible, then there must be a withdrawal of $|\mathcal{V}|-1$ units at
vertex $1$, and every vertex in $\mathcal{V}\setminus\\{1\\}$ must have a path
to reach vertex $1$. Vertex $1$ may be arbitrarily chosen. Stating a well-
known result from graph theory, a tree with $n$ vertices has exactly $n-1$
edges. Hence, the radiality constraints become:

$\displaystyle\mathbf{A}^{T}\mathbf{f}$ $\displaystyle=\mathbf{1}$ (11a)
$\displaystyle-(|\mathcal{V}|-1)\boldsymbol{\theta}\leq$ $\displaystyle\text{
}\mathbf{f}\leq(|\mathcal{V}|-1)\boldsymbol{\theta},$ (11b)
$\displaystyle\boldsymbol{\theta}$ $\displaystyle\in\\{0,1\\}^{|\mathcal{E}|}$
(11c) $\displaystyle\mathbf{1}^{T}\boldsymbol{\theta}$
$\displaystyle=|\mathcal{V}|-1$ (11d) $\displaystyle\mathbf{b^{e}}$
$\displaystyle\leq\boldsymbol{\theta}$ (11e)

Constraints (11a)-(11d) ensure that auxiliary binary indicator variables
$\boldsymbol{\theta}$ on the edge-set of the base graph describe a spanning
tree. Then, (11e) states that the edges selected via the binary variables
$\mathbf{b^{e}}$ are a subset of this spanning tree; and hence form a spanning
forest structure as per Proposition 1. Since the radiality constraints are
thus posed, the number of components in $\mathcal{G}^{\prime}$ need not be
pre-assigned.

###  III-D Connection to grid-forming generators

As stated in Section I, each microgrid should have a grid-forming generator.
Let $\mathcal{V}_{s}$ be the set of buses with grid-forming capabilities.

Then, the connectivity constraint becomes:

$\displaystyle\sum_{e_{jk}\in\mathcal{E}}f^{\prime}_{jk}-\sum_{e_{ij}\in\mathcal{E}}f^{\prime}_{ij}$
$\displaystyle=b^{n}_{j},\forall j\in\mathcal{V}\setminus\mathcal{V}_{s}$
(12a)
$\displaystyle-(|\mathcal{V}|-|\mathcal{V}_{s}|)\mathbf{b^{e}}\leq\mathbf{f^{\prime}}$
$\displaystyle\leq(|\mathcal{V}|-|\mathcal{V}_{s}|)\mathbf{b^{e}}$ (12b)

Constraint (12a) states that every energized non grid-forming node in
$\mathcal{G}$ injects unit virtual commodity into the network. Constraint
(12b) bounds flows on energized lines and fixes flows on deenergized ones at
0. Here, $\mathbf{f^{\prime}}\in\mathds{R}^{|\mathcal{E}|}$ is a vector
representing virtual line flows. It must be emphasized that
$\mathbf{f^{\prime}}$ is different from $\mathbf{f}$ in (11a). Both these
vectors are used to impose connectivity conditions, and have no physical
significance related to the actual power flow. Again, this setup is feasible
only when all units injected by energized non grid-forming nodes can be
withdrawn at grid-forming nodes. Some grid-forming nodes may not be energized
in the optimal topology. This is implicitly considered through the topological
constraint in (1), that ensures all edges connected to a deenergized node are
deenergized as well. Therefore, no path exists from an energized non grid-
forming to a deenergized grid-forming node. If multiple buses in a microgrid
host generators with grid-forming capabilities, only one must be assigned as
the reference bus that determines the system operating point. Power sharing
strategies among multiple dispatchable generators in a microgrid have been
widely studied in literature, see [24] and references therein.

###  III-E Deterministic ODNP

The d-ODNP is solved for one realisation of the power generation capacity
$(\boldsymbol{\xi^{g_{p}}},\boldsymbol{\xi^{g_{q}}})$, and demands
$(\boldsymbol{\xi^{d_{p}}},\boldsymbol{\xi^{d_{q}}})$. The central idea is to
sustain maximum load through microgrids if supply from the main grid is lost,
thereby minimizing service interruption. Therefore, the objective for d-ODNP
becomes maximizing load served. The entire deterministic optimization setup
may be mathematically expressed as follows.

$\displaystyle\min\quad$ $\displaystyle{-\mathbf{1}^{T}\mathbf{p^{d}}}$
(ODNP-1) s. to
$\displaystyle\eqref{constr:topo}-\eqref{LDF},\eqref{eq:connectivity}-\eqref{eq:forest},\eqref{eq:flow_bs}-\eqref{eq:flowbound}$

The relative priority of loads has not been considered in (ODNP-1). However,
this cost may be modified by assigning weights to loads in proportion to their
criticality.

###  III-F Probabilistic ODNP

The problem (ODNP-1) applies to one realization of the generation-demand
scenario
$(\boldsymbol{\xi^{g_{p}}},\boldsymbol{\xi^{g_{q}}},\boldsymbol{\xi^{d_{p}}},\boldsymbol{\xi^{d_{q}}})$.
However, a more realistic goal would be to identified microgrids that are
optimal in some sense for a large set of realizations of the generation-demand
scenarios. In the latter setup, the decision variables
$\boldsymbol{\psi}_{1}:=\\{\mathbf{b^{n},~{}b^{e},~{}f,~{}f^{\prime}}\boldsymbol{,\theta}\\}$
shall remain fixed for all realizations of the uncertainties. The realization
dependent variables would be
$\boldsymbol{\psi}_{2}:=\\{\mathbf{p^{d},~{}q^{d},~{}p^{g},~{}q^{g},~{}v,~{}P,~{}Q}\\}$.
Collecting all the uncertainties in
$\boldsymbol{\xi}:=\\{\boldsymbol{\xi^{g_{p}}},\boldsymbol{\xi^{g_{q}}},\boldsymbol{\xi^{d_{p}}},\boldsymbol{\xi^{d_{q}}}\\}$,
the probabilistic ODNP may seek to solve-

$\displaystyle\min_{\boldsymbol{\psi_{1}}}~{}~{}~{}$
$\displaystyle-\mathds{E}_{\boldsymbol{\xi}}[\mathbf{1}^{T}\mathbf{p^{d}}]$
(ODNP-2) s. to
$\displaystyle\textrm{Pr}\left(\exists~{}\boldsymbol{\psi_{2}}|\mathds{1}(\eqref{constr:topo}-\eqref{LDF},\eqref{eq:connectivity}-\eqref{eq:flowbound})=1\right)\geq
1-\epsilon$

The probabilistic constraint in (ODNP-2) is very difficult to enforce in
practice. However, we will next discuss some reformulations that simplify the
setup without loss of generality. First, note that if the power demands at all
nodes are zero, then for a feasible $\boldsymbol{\psi}_{1}$ satisfying (1) and
(11a)-(12b), there always exist a $\boldsymbol{\psi}_{2}$ that satisfy all
other constraints. Therefore, the probabilistic constraint may be equivalently
posed by enforcing all constraints other than (5a)-(5b) as hard constraints,
and putting the probabilistic requirement on (5a)-(5b). Setting aside (5b) for
expository convenience, notice that the equality constraints in (5a) may be
decomposed into the following inequality constraints.

$\displaystyle p^{d}_{i}-b^{n}_{i}\xi_{i}^{d_{p}}\leq 0,\quad\forall i$ (14a)
$\displaystyle-p^{d}_{i}+b^{n}_{i}\xi_{i}^{d_{p}}\leq 0,\quad\forall i$ (14b)

Now, (14a) can be posed as a hard constraint, leaving (14b) as the main chance
constraint. To reiterate, the ODNP task seeks to identify potential microgrids
within an existing distribution network, such that load served is maximized,
and the microgrids are self-adequate with probability at least
$(1-\varepsilon)$, if the main grid is lost. Islanded microgrids are called
self-adequate when their internal load can be met by their internal
generation. Mathematically,

$\displaystyle Pr(-p^{d}_{i}+b^{n}_{i}\xi_{i}^{d_{p}}\leq 0,\forall i)\geq
1-\varepsilon$ (15)

The self-adequacy condition of microgrids may be thus posed at the node level
since constraints (11a)-(11e) ensure that the network topology is a spanning
forest, and hence loads within a microgrid can be supplied only from
generators within the same microgrid. If needed, one may relax the the
reliability requirement by modifying (15) slightly. For instance, the
condition, a microgrid should be able to meet 90% of its internal load could
be written as $Pr(-p^{d}_{i}+0.9\times b^{n}_{i}\xi_{i}^{d_{p}}\leq 0,\forall
i)\geq 1-\varepsilon$.

###  III-G Sample Average Approximation

Recall from Section II-B that a chance constraint needs to be satisfied with a
probability specified by a risk parameter $\varepsilon$. The chance-constraint
$(C_{2})$ may also be rewritten as $q(\boldsymbol{x})\leq\varepsilon$, where
$q(\boldsymbol{x})=Pr(\boldsymbol{g}(\boldsymbol{x},\boldsymbol{\xi})>\boldsymbol{0})$.
Let $\boldsymbol{\xi^{1}},\boldsymbol{\xi^{2}},...,\boldsymbol{\xi^{N}}$ be
$N$ independent and identically distributed (iid) samples of the uncertainty
vector $\boldsymbol{\xi}$; then $\hat{q_{N}}(\boldsymbol{x})$, an estimator of
$q(\boldsymbol{x})$ is equal to the proportion of realizations in the sample
where
$\boldsymbol{g}(\boldsymbol{x},\boldsymbol{\xi^{i}})>\boldsymbol{0},i=1,..,N$.
This is a sample average approximation of the chance-constrained problem
($P_{1}$) for the samples
$\boldsymbol{\xi^{1}},\boldsymbol{\xi^{2}},\dots,\boldsymbol{\xi^{N}}$:

$\min_{\boldsymbol{x}\in\mathcal{X}}\quad
f(\boldsymbol{x})\quad\textrm{s.to.}\quad\hat{q_{N}}(\boldsymbol{x})\leq\gamma$

Here, $\gamma\in(0,1)$ and is the risk level for the SAA problem. Assuming
that the SAA can be solved, a) if $\gamma<\varepsilon,$ and $N$ is
sufficiently large, SAA is a restriction on the true problem and a feasible
solution of SAA is likely to be feasible for the true problem as well, b) if
$\gamma>\varepsilon$, SAA is a relaxation of the true problem and the optimal
value of SAA is likely to be a lower bound to the optimum for true problem. It
can be shown that for $\gamma=\varepsilon$, the SAA optimum approaches its
true counterpart with probability one as $N$ approaches infinity [15].

The chance-constrained SAA problem $P_{2}$ can be solved using MILP for $N$
iid samples of $\boldsymbol{\xi}$ as shown below[15].

$\displaystyle\min_{\boldsymbol{x}\in\mathcal{X}}\quad$ $\displaystyle
f(\boldsymbol{x})$ ($P_{3}$) s. to
$\displaystyle\boldsymbol{h}(\boldsymbol{x})\leq\boldsymbol{0}$ (16a)
$\displaystyle\boldsymbol{g}(\boldsymbol{x},\boldsymbol{\xi}^{\alpha})\leq
M(1-z_{\alpha}),\quad$ $\displaystyle\alpha=1,2,\dots,N$ (16b)
$\displaystyle\mathbf{1}^{T}\mathbf{z}\geq(1-\gamma)N$ (16c)
$\displaystyle\mathbf{z}\in\\{0,1\\}^{N},\quad$
$\displaystyle\alpha=1,2,\dots,N$ (16d)

Here, $\alpha$ is used to index samples of $\boldsymbol{\xi}$, $z_{\alpha}$ is
a binary variable and $M$ is a large number such that
$M>\max_{\boldsymbol{x}\in\mathcal{X}}\boldsymbol{g}(\boldsymbol{x},\boldsymbol{\xi}^{\alpha})$
for all $\alpha=1,2,\dots,N$. Vector $\mathbf{z}$ stacks all $z_{\alpha}$
values. In constraint (16b), if $z_{\alpha}$=1, then the chance constraint is
not violated. If $z_{\alpha}=0$, then no bound is imposed. The cardinality
constraint in (16c) bounds the proportion of constraint violations.

For ODNP the hard constraints are given by
$\\{\eqref{constr:topo}-\eqref{eq:qgen},\eqref{eq:qload}-\eqref{LDF},\eqref{eq:connectivity}-\eqref{eq:flowbound},\eqref{eq:p_hard}\hskip
2.0pt\forall\alpha\\}$. Probability of violating the chance-constraint
$\\{\eqref{eq:p_chance}\hskip 2.0pt\forall\alpha\\}$ is to be bounded.
Equation (16b) becomes:

$\displaystyle-\mathbf{p^{d}}_{\alpha}+diag(\mathbf{b^{n}})\times\boldsymbol{\xi}^{\mathbf{d_{p}},\alpha}\leq
M(1-z_{\alpha})\times\mathbf{1},\quad\forall\alpha$ (17)

Putting everything together, the problem becomes:

$\displaystyle\min\quad$
$\displaystyle-\frac{1}{N}\sum_{\alpha=1}^{N}\mathbf{1}^{T}\mathbf{p^{d}}_{\alpha}$
(ODNP-3) s. to $\displaystyle\mathbf{p^{d}}_{\alpha}\leq
z_{\alpha}\times\boldsymbol{\xi^{d_{p}}},\quad\forall\alpha$ (18)
$\displaystyle\eqref{constr:topo}-\eqref{eq:qgen},\eqref{eq:connectivity}-\eqref{eq:flowbound},\eqref{eq:p_hard},\eqref{eq:cardinality}-\eqref{eq:G}$

The objective function in (ODNP-3) is the sample-based estimator of the
objective in (ODNP-2) designed to maximize average load served across
considered scenarios. Constraint (18) fixes bus consumptions at zero when
constraint (14b) is not satisfied. This motivates the optimal solution for the
ODNP to be one that also increases $\mathbf{1}^{T}\mathbf{z}$, lower bounded
by $(1-\gamma)N$.

The optimal topology obtained by solving (ODNP-3) is determined by vectors
$\mathbf{b^{n}}^{*}$ and $\mathbf{b^{e}}^{*}$. In practice, only lines
connecting an energized node to a deenergized node will need to be
disconnected to isolate the microgrids.

##  IV Solution Validation

Consider a candidate solution $\boldsymbol{\bar{x}}$ found by the SAA approach
of (ODNP-3). To adjudge its quality, two aspects need to be analyzed: _a)_ Can
it be said with some desired confidence that $\boldsymbol{\bar{x}}$ a feasible
solution for the true problem (ODNP-2)? _b)_ If yes, then how far is
$f(\boldsymbol{\bar{x}})$ from the optimal value $f(\boldsymbol{x}^{*})$? A
method for checking an upper bound of
$Pr\\{\boldsymbol{g}(\boldsymbol{\bar{x}},\boldsymbol{\xi})>\mathbf{0}\\}$ and
lower bound on $f(\boldsymbol{x}^{*})$ is shown in [15] and references
therein.

• Upper bound on violation probability : Consider $N^{\prime}$ iid
realizations of $\boldsymbol{\xi}$, such that $N^{\prime}>>N$, where $N$ is
the number of $\boldsymbol{\xi}$ samples considered for solving the SAA
problem. Here, $N^{\prime}$ may be large as the samples will not be used in
solving an optimization problem and hence not pose computational issues. Let
$\hat{q}_{N^{\prime}}(\boldsymbol{\bar{x}})$ be an estimator of
$q(\boldsymbol{\bar{x}})$; equal to the proportion of times the event
$\mathds{1}(\boldsymbol{g}(\boldsymbol{\bar{x}},\boldsymbol{\xi^{j}})>\boldsymbol{0})=1$
is observed in $N^{\prime}$ trials. Estimator
$\hat{q}_{N^{\prime}}(\boldsymbol{\bar{x}})$ of $q(\boldsymbol{\bar{x}})$ is
unbiased, implying
$\mathds{E}(\hat{q}_{N^{\prime}}(\boldsymbol{\bar{x}}))=q(\boldsymbol{\bar{x}})$.
Also, for a large $N^{\prime}$, its distribution may be approximated by a
normal distribution with mean $q(\boldsymbol{\bar{x}})$ and variance
$q(\boldsymbol{\bar{x}})(1-q(\boldsymbol{\bar{x}}))/N^{\prime}$ [15]. This
yields an approximate $(1-\beta)$-confidence upper bound on
$q(\boldsymbol{\bar{x}})$:

$U_{\beta,N^{\prime}}(\boldsymbol{\bar{x}}):=\hat{q}_{N^{\prime}}(\boldsymbol{\bar{x}})+z_{\beta}\sqrt{\hat{q}_{N^{\prime}}(\boldsymbol{\bar{x}})(1-\hat{q}_{N^{\prime}}(\boldsymbol{\bar{x}}))/N^{\prime}}$

Here, $z_{\beta}=\Phi^{-1}(1-\beta)$, where $\Phi$ is the cumulative
distribution function for the standard normal distribution, $\beta\in(0,1)$.
We compare $U_{\beta,N^{\prime}}(\boldsymbol{\bar{x}})$ to $\varepsilon$ to
check if $\boldsymbol{\bar{x}}$ is a feasible solution.

• Lower bound on optimal value : A procedure for deriving a lower bound for
$f(\boldsymbol{x}^{*})$ is shown in [15]. Let the SAA problem (ODNP-3) be
solved for $N^{\prime\prime}$ iid samples of $\boldsymbol{\xi}$ and risk level
$\gamma\geq 0$; and denote this problem by $P^{N^{\prime\prime}}_{\gamma}$.
Let the true problem (ODNP-2) with risk $\varepsilon$ be denoted as
$P_{\varepsilon}$. Now, the probability that at most $\lfloor\gamma
N^{\prime\prime}\rfloor$ constraint violations are observed in
$N^{\prime\prime}$ trials while solving $P^{N^{\prime\prime}}_{\gamma}$, when
the true violation probability is $\varepsilon$, becomes:

$\Theta_{N^{\prime\prime}}:=B(\lfloor\gamma
N^{\prime\prime}\rfloor;\varepsilon,N^{\prime\prime})$

where,

$B(k;q,N):=\sum_{r=0}^{k-1}\binom{N}{r}q^{r}(1-q)^{N-r}$

is the cumulative density function of the binomial distribution. Say, solving
$P^{N^{\prime\prime}}_{\gamma}$ yields an objective value
$f(\boldsymbol{\bar{x}})$. Assuming $P_{\varepsilon}$ has an optimal solution
$f(\boldsymbol{x}^{*})$, $Pr\\{f(\boldsymbol{\bar{x}})\leq
f(\boldsymbol{x}^{*})\\}\geq\Theta_{N}^{\prime\prime}$. This result yields a
method for obtaining lower bounds with a specified confidence level
$(1-\beta)$. Consider two positive integers $M$ and $N^{\prime\prime}$, such
that $M>N^{\prime\prime}$. Generate $M$ independent sets of $N^{\prime\prime}$
iid samples of $\boldsymbol{\xi}$, and solve the SAA problem for each of the
$M$ sets to obtain values $f(\boldsymbol{\bar{x}}_{j}),j=1,2,\dots,M$. These
can be viewed as iid samples of the random variable
$f(\boldsymbol{\bar{x}}))$. Let $L$ be the largest integer such that
$B(L-1;\Theta_{N}^{\prime\prime},M)\leq\beta$. If the optimal values are
arranged in a non-decreasing order $f(\boldsymbol{\bar{x}}_{1})\leq
f(\boldsymbol{\bar{x}}_{2})\leq\dots f(\boldsymbol{\bar{x}}_{M})$, it can be
shown that with probability at least $(1-\beta)$,
$f(\boldsymbol{\bar{x}}_{L})$ is lower than the true optimal value
$f(\boldsymbol{x^{*}})$.

Note that $f(\boldsymbol{\bar{x}}_{L})>f(\boldsymbol{x^{*}})$ if and only if
more than $L$ of the observed $f(\boldsymbol{\bar{x}}_{j})$ values are greater
than $f(\boldsymbol{x^{*}})$. Considering event
$f(\boldsymbol{\bar{x}}_{j})\leq f(\boldsymbol{x^{*}})$ as a success,
$f(\boldsymbol{\bar{x}}_{L})>f(\boldsymbol{x^{*}})$ if and only if there are
fewer than $L$ successes in $M$ trials, with success probability
$\Theta_{N^{\prime\prime}}$. Probability of fewer than $L$ successes in $M$
trials is $B(L-1,\Theta_{N^{\prime\prime}},M)$, and the bounding procedure
described in this section restricts this probability value to $\beta$.

##  V Numerical Results

Performance of the proposed methodology is illustrated through computational
experiments on a 3.6 GHz Intel Core i7-4790 CPU with 32 GB RAM. Optimization
tasks are solved using YALMIP and Gurobi [25, 26].

###  V-A Experiment Set-up

Figure 2: Modified IEEE 37-bus benchmark feeder with DER locations and
normally open switches

The ODNP problem is solved for a modified version of the IEEE 37-bus benchmark
feeder (fig. 2), converted to its single-phase equivalent by: a) assigning
average three-phase load as bus spot-loads, and b) assigning average three-
phase impedances as line impedances. Four normally open switches are added
(shown with dotted edges in fig. 2). Grid-forming generators are placed at
nodes 742, 718 and 710. PV generators of equal rated capacity are added at
nodes 702, 705, 707, 709 and 737. There are 22 buses with non-zero load. Total
rated capacity of grid-forming and PV generators are 13% and 29% of the rated
system load respectively. Such a contrived feeder model has been intentionally
chosen to capture potential flexibilities of a much larger network. In
practice, a feeder with fewer generators and tie-lines, would have fewer load-
generation scenarios and would be faster to solve for.

Load-generation scenarios were constructed as described next. Data
corresponding to hourly solar generation in California from NREL’s solar power
dataset were used to synthesize five annual generation profiles [27]. The
first 50 generators in the dataset were used; every 10 generators were
aggregated to obtain one profile. The normalized profiles were then scaled to
match the rated capacity of the generators. It is further assumed that the PV
generators are set to work at unity power factor, implying that they do not
participate in reactive power support. This is without loss of generality
since PV generators with reactive power support may be indicated with non-zero
entries in the $\boldsymbol{\xi^{g_{q}}}$ vector. In a similar manner, hourly
load profiles were constructed for residential and commercial buildings in
California with data available from OpenEI [28]. The normalized profiles were
scaled such that the 75th percentile of load data coincided with the nominal
spot load of the corresponding bus. Thus, a total of 8760 scenarios were
constructed for a year; denoted as set $\mathcal{S}$.

###  V-B Chance-Constrained ODNP

As stated previously, the original chance-constrained problem and its SAA
counterpart become equivalent in limit as the number of scenarios considered
$N$ increases. However, a higher $N$ value also results in high computation
time. This increasing trend is illustrated in fig. 5, the markers show median
time for 10 runs conducted over the same scenario sets. For computational
tractability, let the SAA problem be solved on a smaller sample set
$\mathcal{S^{\prime}}\subseteq\mathcal{S}$; if $\mathcal{S^{\prime}}$ is
sufficiently representative of $\mathcal{S}$, then the candidate solution
obtained will be close to the true solution for the original CCO problem.

Performance of the SAA approach is compared to a clustering based methodology,
wherein set $\mathcal{S}$ is divided into clusters and the ODNP task is
designed to yield a solution that holds for some representative samples drawn
from these clusters. Let us call these two Method 1 and Method 2 respectively.

$\bullet$ Method 1: Scenarios are sampled from $\mathcal{S}$ at random with
uniform probability and used to solve (ODNP-3).

$\bullet$ Method 2: Using principal component analysis followed by
hierarchical clustering, set $\mathcal{S}$ is divided into 10 clusters [29].
The scenario clusters are visualized in fig. 5. Once the clusters are
determined, equal number of samples are drawn from each cluster at random.
Evidently samples can only be drawn in multiples of 10. The ODNP is solved
such that the optimal topology is feasible for all selected samples, i.e.
$\gamma=0$.

Figure 3: Increase in median computation time with number of scenarios

Figure 4: 2-D visualization of scenario clustering

(a) Average load served (b) Constraint violation probability

Figure 5: Variation in performance when number of scenarios is varied from 10
to 150. Proposed cc-ODNP is compared to a clustering based method.

The performance of the two solution methods is compared in fig. 5. For SAA,
the value of $\gamma$ used is 0.1. A 95% confidence lower bound on the
objective value is found using the methodology described in Section IV. With
50 runs of independently generated sets of 20 scenarios each, and
$\gamma=0.7$, this lower bound is determined to be -0.18536. The parameters
$M,N^{\prime\prime}$ and $\gamma$ here were chosen following the
recommendations outlined in [30]. The 95% confidence upper bound on
feasibility of the candidate solution $U_{0.05,1000}(\mathbf{\bar{x}})$ is
estimated using a set of 1000 scenarios. Median computation time for the
feasibility checking process was 1.438 seconds.

Evidently, as more scenarios are considered, both average load served and
$U_{0.05,1000}(\mathbf{\bar{x}})$ decrease. The trends are not strict as
additional scenarios can introduce favorable cases with lower cost. Observe
that Method 2 yields a more robust solution (low violation probability) in
lieu of a higher objective cost. Method 1 achieves a cost close to the
theoretical lower bound while bounding supply-deficiency probability to an
acceptable level. Table I summarizes observations when both methods are run
with 100 scenarios.

(a) Average load served (b) Median time taken for 10 runs

Figure 6: Variation in cc-ODNP performance with $\gamma$ considering 50
scenarios.

(a) Average load served (b) Median time taken for 10 runs

Figure 7: Variation in performance with number of normally open switches. Number of scenarios considered is 50, $\gamma=0.1$.  TABLE I: Performance comparison | Method 1 | Method 2
---|---|---
$\mathbf{|f(\bar{x})-f(x^{*})|}$ | 0.00814 | 0.34416
$\mathbf{U_{0.05,1000}(\bar{x})}$ | 0.08 | 0.0025

Figure 8: Optimal microgrid topology considering all switches, 100 scenarios
and $\gamma=0.1$. Average load served is 0.1935 p.u..
$U_{0.05,1000}(\mathbf{\bar{x}})=0.08005$.

Figure 9: Optimal microgrid topology for the base radial network, 100
scenarios and $\gamma=0.1$. Average load served is 0.1483 p.u..
$U_{0.05,1000}(\mathbf{\bar{x}})=0.06465$.

###  V-C Choice of risk parameter

The optimal topology depends highly on the risk parameter. Of course, if a
utility has a high risk budget, they may plan to cover a larger amount of
loads with the microgrids. The risk appetite may be dictated by a number of
factors, like the installed storage capacity and criticality of loads. The
intuition of higher load served with higher risk values is experimentally
verified and shown in fig. 7. It can also be seen that computation time
increases with $\gamma$; possibly because for higher values of $\gamma$, the
feasibility space that the solver has to search for an optimal solution to
ODNP grows in size.

###  V-D Number of Normally Open Switches

Any topology determination problem is combinatorial in nature, and hence the
search space and solution time increases with the number of graph edges. In
the ODNP task, network flexibility may be better utilized to serve more load
by leveraging normally open switches. However, addition of extra edges
introduces additional binary decision variables, thereby increasing solution
time. In fig. 7, it is shown that as more switches are added to the base
radial 37-bus network, ODNP yields higher average load served (i.e. lower
objective values), but the computation time goes up. These data points were
determined by solving the ODNP problem over 50 scenarios sampled with method 1
and using $\gamma=0.1$. For each of these cases, multiple combinations of
normally open switches are possible. However, switches were added one at a
time in a random sequence to the base network for illustration.

###  V-E Microgrid Topology

Optimal microgrids determined for the base radial network with and without
normally open switches are shown in fig. 9 and 9 respectively. Microgrid
components are indicated in color while external elements are in grey. When
all switches are considered, higher load can be served. When only the base
radial network is considered, load served by microgrids is lower, and so is
the supply-deficiency violation probability.

The determined microgrids do not violate the self-adequacy and power systems
constraints for more than $\gamma$ fraction of cases, are radial and contain
at least one grid-forming generator. Notice that despite hosting a solar
generator, bus 707 is not included in any of the microgrids. This may be
because there are no possible ways to connect bus 707 to a grid-forming
generator without violating one of the prescribed constraints.

##  VI Conclusion

The power grid is critical for maintaining essential sectors like healthcare,
transportation and emergency services. This has motivated research towards
boosting grid resilience. Efficiently planned microgrids can help minimize
load interruptions and aid restoration during and after outages. To this end,
this work proposes a chance-constrained optimal network partitioning problem
and presents a computationally tractable solution methodology. Uncertainty in
load and renewable energy generation as well as constraints like maintaining
network radiality and availability of grid-forming generators are addressed.
Experiments on a modified version of the IEEE 37-bus feeder show that good
quality candidate solutions can be found with modest computation cost. Future
work will focus on extending the present planning-stage formulations to multi-
phase topologies and near real-time applications.

## References

  * [1] K. P. Schneider, F. K. Tuffner, M. A. Elizondo, C. Liu, Y. Xu, and D. Ton, “Evaluating the feasibility to use microgrids as a resiliency resource,” _IEEE Trans. Smart Grid_ , vol. 8, no. 2, pp. 687–696, Mar. 2017.
  * [2] _IEEE Guide for Design, Operation, and Integration of Distributed Resource Island Systems with Electric Power Systems_ , IEEE Std. 1547.4-2011.
  * [3] S. Lei, C. Chen, Y. Song, and Y. Hou, “Radiality constraints for resilient reconfiguration of distribution systems: formulation and application to microgrid formation,” _IEEE Trans. Smart Grid_ , vol. 11, no. 5, pp. 3944–3956, 2020.
  * [4] T. Ding, Y. Lin, G. Li, and Z. Bie, “A new model for resilient distribution systems by microgrids formation,” _IEEE Trans. Power Syst._ , vol. 32, no. 5, pp. 4145–4147, Sep. 2017.
  * [5] Z. Wang and J. Wang, “Self-healing resilient distribution systems based on sectionalization into microgrids,” _IEEE Trans. Power Syst._ , vol. 30, no. 6, pp. 3139–3149, Nov. 2015.
  * [6] C. Chen, J. Wang, F. Qiu, and D. Zhao, “Resilient distribution system by microgrids formation after natural disasters,” _IEEE Trans. Smart Grid_ , vol. 7, no. 2, pp. 958–966, Mar. 2016.
  * [7] M. Barani, J. Aghaei, M. Akbari, T. Niknam, H. Farahmand, and M. Korpas, “Optimal partitioning of smart distribution systems into supply-sufficient microgrids,” _IEEE Trans. Smart Grid_ , vol. 10, no. 3, pp. 2523–2533, May 2019.
  * [8] R. Jovanovic, A. Bousselham, and S. Voss, “Partitioning of supply/demand graphs with capacity limitations: an ant colony approach,” _Journal of Combinatorial Optimization_ , vol. 35, pp. 224–249, Jan. 2018.
  * [9] R. A. Osama, A. F. Zobaa, and A. Y. Abdelaziz, “A planning framework for optimal partitioning of distribution networks into microgrids,” _IEEE Systems Journal_ , vol. 14, no. 1, pp. 916–926, Mar. 2020.
  * [10] S. A. Arefifar, Y. Mohamed, and T. H. El-Fouly, “Supply-adequacy-based optimal construction of microgrids in smart distribution systems,” _IEEE Trans. Smart Grid_ , vol. 3, no. 3, pp. 1491–1502, May 2012.
  * [11] M. E. Nassar and M. M. Salama, “Adaptive self-adequate microgrids using dynamic boundaries,” _IEEE Trans. Smart Grid_ , vol. 7, no. 1, pp. 105–113, Jan. 2016.
  * [12] S. Biswas, E. Bernabeu, and D. Picarelli, “Optimal distribution system restoration with microgrids and distributed generators,” in _IEEE Power Energy Society Innovative Smart Grid Technologies Conference (ISGT)_ , Washington DC, Feb. 2020, pp. 1–5.
  * [13] M. K. Singh, V. Kekatos, and C. Liu, “Optimal distribution system restoration with microgrids and distributed generators,” in _IEEE PES General Meeting_ , Atlanta, GA, Aug. 2019, pp. 1–5.
  * [14] G. Patsakis, D. Rajan, I. Aravena, and S. Oren, “Strong mixed-integer formulations for power system islanding and restoration,” _IEEE Trans. Power Syst._ , vol. 34, no. 6, pp. 4880–4888, Nov. 2019.
  * [15] S. Ahmed and A. Shapiro, “Solving chance-constrained stochastic programs via sampling and integer programming,” in _State-of-the-Art Decision-Making Tools in the Information-Intensive Age_. Informs, 2008, pp. 261–269.
  * [16] C. Godsil and G. Royle, _Algebraic Graph Theory_. New York, NY: Springer, 2001.
  * [17] X. Geng and L. Xie, “Data-driven decision making in power systems with probabilistic guarantees: Theory and applications of chance-constrained optimization,” _Annual Reviews in Control_ , vol. 47, pp. 341 – 363, 2019\.
  * [18] _American National Standard for Electric Power Systems and Equipment—Voltage Ratings (60 Hertz)_ , ANSI Std. C84.1-1995, 2016.
  * [19] R. A. Jabr, “Linear decision rules for control of reactive power by distributed photovoltaic generators,” _IEEE Trans. Power Syst._ , vol. 33, no. 2, pp. 2165–2174, Mar. 2018.
  * [20] M. E. Baran and F. F. Wu, “Network reconfiguration in distribution systems for loss reduction and load balancing,” _IEEE Trans. Power Del._ , vol. 4, no. 2, pp. 1401–1407, Apr. 1989.
  * [21] S. Bolognani and F. Dorfler, “Fast power system analysis via implicit linearization of the power flow manifold,” in _Annual Allerton Conf. on Communication, Control, and Computing_ , Allerton, IL, Sep. 2015, pp. 402–409.
  * [22] G. P. McCormick, “Computability of global solutions to factorable nonconvex programs: Part I - convex underestimating problems,” _Mathematical Programming_ , vol. 10, no. 1, pp. 147–175, Dec. 1976.
  * [23] M. K. Singh, V. Kekatos, S. Taheri, K. P. Schneider, and C.-C. Liu, “Enforcing radiality constraints for DER-aided power distribution grid reconfiguration,” _arXiv preprint 1910.03020_ , 2019.
  * [24] S. Ahn, J. Park, I. Chung, S. Moon, S. Kang, and S. Nam, “Power-sharing method of multiple distributed generators considering control modes and configurations of a microgrid,” _IEEE Trans. Power Del._ , vol. 25, no. 3, pp. 2007–2016, 2010.
  * [25] J. Löfberg, “YALMIP : A toolbox for modeling and optimization in MATLAB,” in _Proc. of the CACSD Conference_ , Taipei, Taiwan, 2004.
  * [26] Gurobi Optimization LLC, “Gurobi Optimizer Ref. Manual,” 2020.
  * [27] NREL, “Solar power data for integration studies.” [Online]. Available: https://www.nrel.gov/grid/solar-power-data.html
  * [28] OpenEI, “Commercial and residential hourly load profiles for all tmy3 locations in the United States.” [Online]. Available: https://openei.org/datasets/files/961/pub/
  * [29] P.-N. Tan, M. Steinbach, and V. Kumar, _Introduction to Data Mining_. Addison-Wesley Longman Publishing Co., Inc., 2005.
  * [30] A. Nemirovski and A. Shapiro, “Convex approximations of chance constrained programs,” _SIAM J. on Optimization_ , vol. 17, p. 969–996, Dec. 2006.