Radiality Constraints for Resilient Reconfiguration of Distribution Systems: Formulation and Application to Microgrid Formation

TL;DR

This paper introduces a novel radiality constraint formulation for distribution system reconfiguration, enhancing flexibility and optimality, especially in microgrid formation for resilience, validated through theoretical and case study analyses.

Contribution

It proposes a new radiality constraints method that fully enables topological flexibilities in distribution systems, improving feasibility and optimality in reconfiguration problems.

Findings

Enhanced microgrid formation model with flexible merge/separation

The new constraints improve solution feasibility and optimality

Case studies demonstrate superiority over existing models

Abstract

Network reconfiguration is an effective strategy for different purposes of distribution systems (DSs), e.g., resilience enhancement. In particular, DS automation, distributed generation integration and microgrid (MG) technology development, etc., are empowering much more flexible reconfiguration and operation of the system, e.g., DSs or MGs with flexible boundaries. However, the formulation of DS reconfiguration-related optimization problems to include those new flexibilities is non-trivial, especially for the issue of topology, which has to be radial. That is, existing methods of formulating radiality constraints can cause underutilization of DS flexibilities. Thus, this work proposes a new method for radiality constraints formulation fully enabling the topological and some other related flexibilities of DSs, so that the reconfiguration-related optimization problems can have extended feasibility and enhanced optimality. Graph-theoretic supports are provided to certify its theoretical validity. As integer variables are involved, we also analyze the tightness and compactness issues. The proposed radiality constraints are specifically applied to post-disaster MG formation, which is involved in many DS resilience-oriented service restoration and/or infrastructure recovery problems. The resulting new MG formation model, which allows more flexible merge and/or separation of sub-grids, etc., establishes superiority over the models in the literature. Case studies are conducted on two test systems.

Figures

Tables

Table 1. TABLE I: Comparing Different Types of Spanning Tree Constraints

Type of formu- lation/model	Number of variables	Number of constraints	Spanning tree polytope ?
Subtour elimination	$\left|𝑳\right|$	$𝒪\left(2^{\left|𝑵\right|}\right)$	Yes
Directed cutset	$3$$\cdot$$\left|𝑳\right|$	$𝒪\left(2^{\left|𝑵\right|-1}\right)$	Yes (extended formulation)
Single-com- modity flow	$3$$\cdot$$\left|𝑳\right|$	$\left|𝑵\right|$$+$$2$$\cdot$$\left|𝑳\right|$	No
Primal & dual graphs-based	$4$$\cdot$$\left|𝑳\right|$	$2$$\cdot$$\left|𝑳\right|$	Yes (extended formulation)
Directed multi- commodity flow	$2$$\cdot$$\left|𝑵\right|$$\cdot$$\left|𝑳\right|$ $+$$\left|𝑳\right|$	$\left|𝑵\right|$2$+$$2$$\cdot$$\left|𝑵\right|$$\cdot$$\left|𝑳\right|$ $-$$\left|𝑵\right|$$-$$\left|𝑳\right|$$+$$1$	Yes (extended formulation)

Table 2. TABLE II: Notations of the Spanning Tree Constraints ( 3 )-( 9 )

Parameters
$𝑵/𝑳$	Set of all DS nodes/branches.
$i_r$	Index for the substation node.
Variables
$f_{ij}^k$	Flow of commodity $k$ from node $i$ to node $j$.
${\lambda }_{ij}$	Binary, $1$ if arc $(i,j)$ is included in the directed spanning tree, $0$ otherwise.
${\beta }_{ij}$	Binary, $1$ if the fictitious connection status of branch $(i,j)$ is closed, $0$ if open.

Table 3. TABLE III: Notations of the MG Formation Model ( 10 )-( 21 )

Parameters
${𝑵}_{𝒓}/{𝑵}_{𝒈}$	Set of substation nodes/DG nodes.
${𝑵}_{𝒐}/{𝑵}_{𝒄}$	Set of nodes with faulted open/closed load switches.
${𝑳}_{𝒐}/{𝑳}_{𝒄}$	Set of faulted open/closed branches.
$p_i^c/q_i^c$	Real/reactive power demand of the load at node $i$.
${\overline{p}}_i^g/{\overline{q}}_i^g$	Real/reactive power capacity of the power source at node $i$.
${\overline{v}}_i/{\underset{¯}{v}}_i$	Maximum/minimum squared voltage magnitude of node $i$.
$r_{ij}/x_{ij}/{\overline{S}}_{ij}$	Resistance/reactance/apparent power capacity of branch $(i,j)$.
${\omega }_i$	Priority weight of the load at node $i$.
${\vartheta }_i$	Number of branches starting or ending with node $i$.
$M$	A large enough positive number.
Variables
${\delta }_i$	Binary, $1$ if the load at node $i$ is picked up, $0$ otherwise.
${\epsilon }_i$	Binary, $1$ if node $i$ is energized, $0$ otherwise.
${\alpha }_{ij}$	Binary, $1$ is branch $(i,j)$ is closed, $0$ if open.
$p_i^g/q_i^g$	Real/reactive power output of the power source at node $i$.
$v_i$	Squared voltage magnitude of node $i$.
$P_{ij}/Q_{ij}$	Real/reactive power flow on branch $(i,j)$.

Table 4. TABLE IV: Comparing Different MG Formation Models

———	Model in [18]	Model in [19]	Proposed model
Applicable systems	Radial	Radial or meshed	Radial or meshed
Radiality constraints	——	Single-commo- dity flow-based	Proposed radi- ality constraints
Allocation of DGs into MGs	One DG in each MG	One DG in each MG	Flexible
Number of MGs	Fixed	Fixed	Flexible
Unenergized islands	Not allowed	Not allowed	Allowed
Loads w/ faulted closed switches	Forced to pick up	Forced to pick up	Flexible

Table 5. TABLE V: Computation Time and the Number of Infeasible Cases (IEEE 33-Node Test System)

———	Model in [18]	Model in [19]	Proposed model
Avg. comput. time	0.56 s	0.75 s	0.49 s
Infeasible cases	212/10000	272/10000	2/10000

Table 6. TABLE VI: Summary Statistics for the Restored Loads (IEEE 33-Node Test System)

——	Restored loads (kW)					Box plots
	avg.	std.	max.	med.	min.
1: Model in [18]	2686	531	3715	2750	1175
2: Model in [19]	2849	495	3715	2915	1380
3: Proposed model	3088	393	3715	3165	1425

Table 7. TABLE VII: Summary Statistics for the DG Capacity Utilization Rate (IEEE 33-Node Test System)

——	DG capacity utilization rate					Box plots
	avg.	std.	max.	med.	min.
1: Model in [18]	64.4%	9.8%	88.4%	65.1%	29.7%
2: Model in [19]	60.7%	15.1%	94.8%	61.0%	15.3%
3: Proposed model	74.1%	15.7%	100%	74.4%	6.6%

Table 8. TABLE VIII: Computation Time and the Number of Explored Nodes in the B&C Search Tree (IEEE 33-Node Test System)

Avg. $T_s$	Avg. $T_m$	Number of cases with:
			$T_s=T_m$	$T_s>T_m$
0.49 s	0.56 s	8911/10000	0/10000	1089/10000
Avg. $N_s$	Avg. $N_m$	Number of cases with:
			$N_s=N_m$	$N_s>N_m$
407	143	402/10000	21/10000	9577/10000
$T_s/T_m$ (resp. $N_s/N_m$): Computation time (resp. the number of ex- plored nodes in the B&C search tree) of the proposed MG formation model using the single-commodity flow-based method$/$directed multi- commodity flow-based method to formulate constraint (1).

Table 9. TABLE IX: Computation Time and the Number of Infeasible Cases (IEEE 123-Node Test System)

———	Model in [18]	Model in [19]	Proposed model
Avg. comput. time	13.41 s	4.36 s	1.79 s
Infeasible cases	109/10000	196/10000	2/10000

Table 10. TABLE X: Summary Statistics for the Restored Loads (IEEE 123-Node Test System)

——	Restored loads (kW)					Box plots
	avg.	std.	max.	med.	min.
1: Model in [18]	2917	482	4150	2915	1555
2: Model in [19]	3008	524	4350	3010	1265
3: Proposed model	3189	494	4450	3195	1665

Table 11. TABLE XI: Summary Statistics for the DG Capacity Utilization Rate (IEEE 123-Node Test System)

——	DG capacity utilization rate					Box plots
	avg.	std.	max.	med.	min.
1: Model in [18]	63.0%	9.1%	90.5%	63.1%	30.8%
2: Model in [19]	63.3%	10.6%	97.4%	63.1%	30.1%
3: Proposed model	70.9%	12.4%	100%	70.7%	24.9%

Table 12. TABLE XII: Computation Time and the Number of Explored Nodes in the B&C Search Tree (IEEE 123-Node Test System)

Avg. $T_s$	Avg. $T_m$	Number of cases with:
			$T_s=T_m$	$T_s>T_m$
1.79 s	1.54 s	7348/10000	0/10000	2652/10000
Avg. $N_s$	Avg. $N_m$	Number of cases with:
			$N_s=N_m$	$N_s>N_m$
496	71	108/10000	134/10000	9758/10000
$T_s/T_m$ (resp. $N_s/N_m$): Computation time (resp. the number of ex- plored nodes in the B&C search tree) of the proposed MG formation model using the single-commodity flow-based method$/$directed multi- commodity flow-based method to formulate constraint (1).

Equations

$$\mathbf{\beta }\in \mathbf{\Omega }.$$

$${\alpha }_{ij}\le {\beta }_{ij},\forall (i,j)\in \mathbf{L}.$$

$$\sum _{(j,i_r)\in \mathbf{L}}{f}_{j{i}_r}^k-\sum _{(i_r,j)\in \mathbf{L}}{f}_{i_{r}j}^k=-1,\forall k\in \mathbf{N}\setminus i_r.$$

$$\sum _{(j,i)\in \mathbf{L}}{f}_{ji}^k-\sum _{(i,j)\in \mathbf{L}}{f}_{ij}^k=0,\forall k\in \mathbf{N}\setminus i_r,\forall i\in \mathbf{N}\setminus \{i_r,k\}.$$

$$\sum _{(j,k)\in \mathbf{L}}{f}_{jk}^k-\sum _{(k,j)\in \mathbf{L}}{f}_{kj}^k=1,\forall k\in \mathbf{N}\setminus i_r.$$

$$0\le f_{ij}^k\le {\lambda }_{ij},0\le f_{ji}^k\le {\lambda }_{ji},\forall k\in \mathbf{N}\setminus i_r,\forall (i,j)\in \mathbf{L}.$$

$$\sum _{(i,j)\in \mathbf{L}}({\lambda }_{ij}+{\lambda }_{ji})=|\mathbf{N}|-1.$$

$${\lambda }_{ij}+{\lambda }_{ji}={\beta }_{ij},\forall (i,j)\in \mathbf{L}.$$

$${\lambda }_{ij},{\lambda }_{ji}\in \{0,1\},\forall (i,j)\in \mathbf{L}.$$

$$\max \sum _{i\in \mathbf{N}}{\delta }_i\cdot {\omega }_i\cdot p_i^c$$

$$s.t.\text{eqref}SF1-\text{eqref}SF2.$$

$$p_i^g-{\delta }_i\cdot p_i^c+\sum _{(j,i)\in \mathbf{L}}{P}_{ji}-\sum _{(i,j)\in \mathbf{L}}{P}_{ij}=0,\forall i\in \mathbf{N}.$$

$$q_i^g-{\delta }_i\cdot q_i^c+\sum _{(j,i)\in \mathbf{L}}{Q}_{ji}-\sum _{(i,j)\in \mathbf{L}}{Q}_{ij}=0,\forall i\in \mathbf{N}.$$

$$p_i^g=q_i^g=0,\forall i\in \mathbf{N}\setminus \{{\mathbf{N}}_{\mathbf{r}},{\mathbf{N}}_{\mathbf{g}}\}.$$

$$0\le p_i^g\le {\overline{p}}_i^g,0\le q_i^g\le {\overline{q}}_i^g,\forall i\in \{{\mathbf{N}}_{\mathbf{r}},{\mathbf{N}}_{\mathbf{g}}\}.$$

$$v_i-v_j\ge 2\cdot (P_{ij}\cdot r_{ij}+Q_{ij}\cdot x_{ij})+({\alpha }_{ij}-1)\cdot M,$$

$$v_i-v_j\le 2\cdot (P_{ij}\cdot r_{ij}+Q_{ij}\cdot x_{ij})+(1-{\alpha }_{ij})\cdot M,$$

$$\forall (i,j)\in \mathbf{L}.$$

$${\underline{v}}_i\le v_i\le {\overline{v}}_i,\forall i\in \mathbf{N}.$$

$$P_{ij}^2+Q_{ij}^2\le {\alpha }_{ij}\cdot {\overline{S}}_{ij}^2,\forall (i,j)\in \mathbf{L}.$$

$${\alpha }_{ij}=0,\forall (i,j)\in {\mathbf{L}}_{\mathbf{o}};{\alpha }_{ij}=1,\forall (i,j)\in {\mathbf{L}}_{\mathbf{c}}.$$

$${\delta }_i=0,\forall i\in {\mathbf{N}}_{\mathbf{o}};{\delta }_i\ge {\epsilon }_i,\forall i\in {\mathbf{N}}_{\mathbf{c}}.$$

$${\epsilon }_i=1,\forall i\in \{{\mathbf{N}}_{\mathbf{r}},{\mathbf{N}}_{\mathbf{g}}\}.$$

$$(\sum _{(i,j)\in \mathbf{L}}{\epsilon }_j\cdot {\alpha }_{ij}+\sum _{(j,i)\in \mathbf{L}}{\epsilon }_j\cdot {\alpha }_{ji})/{\vartheta }_i\le {\epsilon }_i$$

$$\le \sum _{(i,j)\in \mathbf{L}}{\epsilon }_j\cdot {\alpha }_{ij}+\sum _{(j,i)\in \mathbf{L}}{\epsilon }_j\cdot {\alpha }_{ji},\forall i\in \mathbf{N}\setminus \{{\mathbf{N}}_{\mathbf{r}},{\mathbf{N}}_{\mathbf{g}}\}.$$

Full text

I Introduction

Distribution system (DS) reconfiguration is an effective and multi-function strategy [1, 2, 3]. Its optimization thus has been extensively studied. As most DSs have to operate with a radial topology, the mathematical formulation of radiality constraints has been specifically investigated [4, 5, 6, 7, 8] (see more details in Section III). This issue is resolved for conventional DS reconfiguration. Nevertheless, now DSs can adopt much more adaptive reconfiguration and operation owing to the added flexibilities of automation equipment, distributed generations (DGs), and microgrid (MG) components, etc. In particular, regarding the topology issue, DSs and MGs now can have flexible boundaries [9, 10, 11]. For example, the DS is split into a to-be-optimized number of MGs in [11]. Such added flexibilities are actually empowering more resilient reconfiguration of DSs.

However, existing methods to formulate radiality constraints cannot fully include those new flexibilities in optimizing DS reconfiguration [12]. That is, the feasible region given by these formulations only refers to a subset of the actual one. DS flexibilities thus will be underutilized and less coordinated, which is especially adverse for resilience-oriented reconfiguration [13]. For example, post-disaster reconfiguration is faced with quite limited flexibilities due to many faults caused by an extreme weather event [14], etc. In such cases, the faults also complicates decision-making and even requires co-optimization with other recovery efforts [15], etc. Actually, to co-optimize with the repairing sequence of the damaged parts involves reconfiguring a DS network with a physical structure evolving with the repairing variables. Existing methods to formulate radiality constraints cannot properly handle these situations.

This work proposes a new formulation of radiality constraints that fully enables topological and some other related flexibilities in DS reconfiguration-related optimization problems. It is superior to the literature’s other attempts with the same or similar aims [12] [16] (see comparisons in Section II). Graph-theoretic justifications are provided to affirm the analytical validity of it. As integer variables are involved, the tightness and compactness issues are also analyzed [17]. Generally, adopting the proposed radiality constraints in DS reconfiguration optimization problems can attain extended feasibility and enhanced optimality.

For verification, the proposed radiality constraints are applied to construct a new optimization model for post-disaster MG formation, which reconfigures the DS to form MGs energized by DGs and/or other power sources. It is an essential strategy for many resilience-oriented DS restoration and/or recovery problems [15, 16, 14, 12, 13]. Our proposed model again establishes superiority over the literature’s two groups of MG formation models [18] [19]. As compared in Fig. 1, while their models exclude some DS flexibilities, our model allows more flexible merge and separation of sub-grids, etc. (See more details in Section IV.)

In the following, Section II details the proposed method of formulating radiality constraints. Section III analyzes its tightness and compactness issues. Sections IV and V apply it to resilient MG formation. Section VI provides the conclusion.

II Proposed Radiality Constraints

This work proposes to construct radiality constraints based on two simple graph-theoretic concepts and their relationships.

First, the definition of spanning tree, which may be already well-known, is still given here to make this paper self-contained:

Definition 1

A spanning tree is a graph that connects all the vertices and contains no cycles [20].

Second, spanning forest, the other involved concept which is less common, is defined as follows:

Definition 2

A spanning forest is a graph with no cycles [20].

As in Fig. 2, a spanning forest is a graph whose connected components are spanning trees. Actually, a spanning forest with $\kappa$ components can be called as a $\kappa$-tree. For example, Fig. 2(b) is a $3$-tree. And, a $1$-tree is a spanning tree, e.g., Fig. 2(a). Obviously, both spanning trees and spanning forests are radial.

Let $\kappa_{1}$ be the number of substations in the DS. In light of Definitions 1 and 2, we can observe that normal DS reconfiguration for loss reduction [1] and supply capacity improvement [2], etc., is essentially forming a $\kappa_{1}$-tree. That is, a spanning tree is formed if $\kappa_{1}=1$, while a spanning forest is formed if $\kappa_{1}>1$. It is also required that each component has a substation node.

Regarding resilient DS reconfiguration (e.g., MG formation for service restoration), it is to form a $\kappa_{2}$-tree with $\kappa_{2}\geq\kappa_{1}$. In this case, it is required that each component has at most one substation node, if any. Actually, the value of $\kappa_{2}$ should be optimized in resilient reconfiguration. If it is predefined, the optimality can be impacted. In all relevant cases, we will have $\kappa_{2}>1$, which indicates that a spanning forest is formed. That is, the topology issues in resilient DS reconfiguration can be resolved by requiring the network to be a spanning forest. Thus, we can formulate the radiality constraints as equations (1)-(2), which are inspired by the remark as below:

Remark 1

An arbitrary subgraph of a spanning tree is a spanning forest.

Specifically, a subgraph of a graph consists of a subset of the vertices and edges in the graph. For example, Fig. 2(b) is a subgraph of the graph in Fig. 2(a). According to Definitions 1 and 2, Remark 1 is naturally true. Based on that, the proposed radiality constraints are formulated as follows:

[TABLE]

[TABLE]

In (1)-(2), $\boldsymbol{L}$ is the set of DS branches; $\boldsymbol{\alpha}\triangleq\{\alpha_{ij},\forall(i,j)\in\boldsymbol{L}\}$, where $\alpha_{ij}$ is the connection status of branch $(i,j)$ ($1$ if closed, [math] if open); $\boldsymbol{\beta}\triangleq\{\beta_{ij},\forall(i,j)\in\boldsymbol{L}\}$, where $\beta_{ij}$ is the fictitious connection status of branch $(i,j)$ ($1$ if closed, [math] if open). By “fictitious”, it indicates that $\boldsymbol{\beta}$ are just auxiliary variables, which do not actually determine the network topology of the DS. The topology is still determined by variables $\boldsymbol{\alpha}$. The symbol $\boldsymbol{\Omega}$ is the set of all incidence vectors of spanning tree topologies that the network can form via reconfiguration (see Section III).

Thus, constraint (1) enforces $\boldsymbol{\beta}$ to form a fictitious spanning tree. Constraint (2) then restricts the DS to close a subset of the closed branches in the spanning tree determined by $\boldsymbol{\beta}$. That is, $\boldsymbol{\alpha}$ form a subgraph of the fictitious spanning tree. Remark 1 indicates that the resulting network topology determined by $\boldsymbol{\alpha}$ is a spanning forest. Note that constraint (1) is just expressed conceptually here. We will elaborate on its explicit formulations in Section III. Besides, constraints (1)-(2) implicitly assume that the DS has only one substation node. Section III will briefly explain the reason for this assumption, and introduce a simple method that enables the proposed radiality constraints to fully handle a DS with multiple substation nodes. For uncontrollable branches without switches or with faulted switches, constraints can be added to specify their connection status (see Section IV).

Remark 1 only specifies that any $\boldsymbol{\alpha}$ satisfying constraints (1)-(2) is topologically feasible for the DS. The following theorem indicates that any topologically feasible $\boldsymbol{\alpha}$ satisfies constraints (1)-(2). As it is somewhat less straightforward, its proof is also provided. Remark 1 and Theorem 1 together certify the validity of the proposed radiality constraints.

Theorem 1

A spanning forest subgraph of a connected graph is also the subgraph of a spanning tree subgraph of the connected graph.

Proof: Assume a $\kappa$-tree denoted as $\mathcal{G}_{\kappa}$ to be a subgraph of a graph $\mathcal{G}$. As $\mathcal{G}$ is connected, there exists at least one edge in $\mathcal{G}$ that can link the $\kappa$th component of $\mathcal{G}_{\kappa}$ to another component of $\mathcal{G}_{\kappa}$. Adding this edge to $\mathcal{G}_{\kappa}$, which does not create cycles, $\mathcal{G}_{\kappa}$ becomes a $(\kappa-1)$-tree denoted as $\mathcal{G}_{\kappa-1}$. Repeating this process, $\mathcal{G}_{\kappa}$ ends up to be a $1$-tree denoted as $\mathcal{G}_{1}$. Obviously, the spanning forest $\mathcal{G}_{\kappa}$ is a subgraph of the spanning tree $\mathcal{G}_{1}$, which is a subgraph of $\mathcal{G}$. This completes the proof.        $\square$

The proposed radiality constraints (1)-(2) can be interpreted as a two-step method to regulate DS topology in reconfiguration. In the first step, constraint (1) ensures network radiality by having a spanning tree to be the DS topology’s supergraph, which is determined by the fictitious connection status of branches. (If $\mathcal{G}_{A}$ is a subgraph of $\mathcal{G}_{B}$, then $\mathcal{G}_{B}$ is said to be a supergraph of $\mathcal{G}_{A}$.) In the second step, constraint (2) enables more flexible reconfiguration by allowing the DS to select a subgraph of the fictitious spanning tree to be the actual network topology.

By contrast, common methods of formulating radiality constraints in the literature (e.g., [6] [8]) can be seen as a one-step process that directly enforces the DS topology to be a spanning tree or a spanning forest. Their models (the single-commodity flow model [6], etc.) can also be used to formulate constraint (1) here. Nevertheless, a new formulation will be presented in Section III. Above all, the proposed two-step method for radiality constraints formulation enables many more flexibilities in DS reconfiguration, such as more adaptive merge or separation of sub-grids, and more flexible allocation of power sources into sub-grids. Such benefits will be detailed in Sections IV and V applying the proposed radiality constraints to resilient MG formation.

Due to the critical importance of network reconfiguration in DS resilience enhancement, etc., many publications also propose to develop new methods of ensuring radiality to allow undiminished flexibilities and adaptivities of DSs in optimization. Specifically, references [12][16] have presented two applicable methods, which are compared to our proposed method as below:

II-1 Validity

Remark 1 and Theorem 1 theoretically prove the validity of constraints (1)-(2) here. References [12] [16] however lack such analytical proofs for their proposed formulations.

II-2 Tightness

With constraint (1) properly formulated, our proposed model attains the tightest formulation of the spanning forest polytope of the DS, i.e., the convex hull of incidence vectors of possible spanning forest topologies (see Section III). The methods in [12] [16] produce relatively less tight formulations.

II-3 Compactness

Constraint (1) can also be formulated in less tight but more compact manners (see Section III). In this regard, the resulting topology constraints (1)-(2) are generally more compact than (i.e., with fewer variables and constraints), and still as tight as or even tighter than those of [12] [16].

II-4 Application convenience

The application of our proposed method can be quite straightforward, i.e., simply adding constraint (2) to the commonly-used single-commodity flow-based radiality constraints. The methods in [12] [16] however involve the introduction of a virtual source node and virtual branches, and the modeling of a virtual DC optimal power flow subproblem and its Karush-Kuhn-Tucker conditions etc., respectively.

II-5 Applicability

The proposed radiality constraints (1)-(2) can be adopted in different DS optimization problems involving reconfiguration. Topological and some other related flexibilities thus can be fully enabled in optimization. The methods in [12] [16] have some limitations in this regard. For example, the radiality constraints proposed in [12] allow the merge of sub-grids in optimization, but do not enable their possible separation.

III Tightness and Compactness Issues

Constraint (1) is only expressed conceptually in Section II. As aforementioned, common methods or models for representing radiality constraints in DS reconfiguration-related publications can be used for its explicit formulation. Those methods and models are revisited here. A new model is also presented. The tightness and compactness issues are specifically discussed.

DS reconfiguration is essentially a mixed-integer programming (MIP) problem generally solved by the branch-and-cut (B&C) method, which needs to solve its linear programming (LP) relaxations (i.e., relaxing integer$/$binary constraints). The computational complexity of a MIP depends on its tightness and compactness, etc. A MIP formulation is tight if its feasible region is similar to that of its LP relaxation, contributing to a smaller gap between the optimal values (solutions) of the MIP and its LP relaxation, and less explored nodes in the B&C search tree (i.e., fewer iterations and faster convergence). A compact MIP formulation has a small number of variables and constraints, leading to shorter computation time for each explored node in the B&C search tree. Tightness and compactness are often conflicting objectives in formulating a MIP [17].

Several relevant concepts are further introduced. Here, incidence vectors of spanning trees are values of $\boldsymbol{\beta}$ defining fictitious spanning tree topologies of the DS; and spanning tree polytope is the convex hull of such incidence vectors. For example, the system in Fig. 3(a) has $\boldsymbol{\beta}=[\beta_{12},\beta_{13},\beta_{23}]$. Its incidence vectors of spanning trees include $[1,1,0]$, $[1,0,1]$ and $[0,1,1]$, and its spanning tree polytope is the shaded region with those vectors as the vertices in Fig. 3(b). Incidence vectors of spanning forests and spanning forest polytope are defined likewise.

That is, constraint (1) actually requires $\boldsymbol{\beta}$ to be an incidence vector of a spanning tree. A tight and compact formulation is desired to explicitly represent this requirement. LP relaxations of the tightest formulations define the spanning tree polytope.

Some methods or models to express radiality constraints in DS reconfiguration-related publications can be used to formulate the spanning tree constraint (1). They are revisited here:

III-1 Loop-eliminating method [4]

It enumerates all loops of the DS, and enforces each one to be open. However, finding all loops in a graph is NP-hard. This method is essentially equivalent to the subtour elimination formulation of spanning tree constraints in graph theory. As one of the tightest formulations, its LP relaxation defines the spanning tree polytope and thus has integer extreme points. It has only $|\boldsymbol{L}|$ variables but also an exponential number of constraints, which limits its application.

III-2 Path-based model [5]

It enumerates all paths to the substation for each node, and activates only one of them. Still, finding all paths between two graph nodes is NP-hard. This model also has a counterpart in graph theory, i.e., the directed cutset formulation of spanning tree constraints with 3$$\cdot$$|\boldsymbol{L}| variables and an exponential number of constraints. It is one of the tightest formulations with integer vertices of its LP relaxation, too.

III-3 Single-commodity flow-based model [6]

It closes |\boldsymbol{N}|$$-$$1 branches ($\boldsymbol{N}$: set of DS nodes), and ensures connectivity by imposing unit fictitious flow from the substation to each node. With 3$$\cdot$$|\boldsymbol{L}| (or reduced to 2$$\cdot$$|\boldsymbol{L}|) variables and a linear number of constraints, it is one of the most compact formulations. Due to its simplicity, it is also the most commonly used DS radiality model. Still, it is less tight. The projection of its LP relaxation into the $\boldsymbol{\beta}$-space defines a region larger than the spanning tree polytope, and therefore has fractional extreme points.

III-4 Parent-child relation-based method [7]

It instructs every node except the substation to have one parent. However, this method can produce a disconnected graph with loops [8].

III-5 Primal and dual graphs-based model [8]

It froces each primal (dual) node to be connected to another primal (dual) node, and forbids the primal and dual spanning trees to be intersected. With 4$$\cdot$$|\boldsymbol{L}| variables and a linear number of constraints, it is one of the most compact formulations. Originally proposed in [22], it actually has a totally unimodular constraint matrix to define a polyhedron with integer vertices. Thus, it is also one of the tightest formulations with its LP relaxation defining the spanning tree polytope in the $\boldsymbol{\beta}$-space. Still, it is only applicable to planar graphs, and needs to construct dual graphs. It also has a flaw to be avoided by appropriate selection of roots [23]. Generally, its applicability is limited and its use is not convenient.

Using notations in Table II, we further introduce a new flow-based formulation of the spanning tree constraint (1) as below:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The above formulation is called as the directed multicommodity flow-based model of the spanning tree constraints. It defines a fictitious commodity for each node $k\neq i_{r}$, and enforces 1 unit of commodity $k$ delivered from the substation node $i_{r}$ to node $k$. Constraint (6) implies that each commodity can flow on an arc only if the arc is included in the directed spanning tree defined by variables $\lambda_{ij}$. Other equations are self-explanatory.

The LP relaxation of the above formulation (3)-(9) defines the spanning tree polytope (i.e., $\mathit{conv}(\boldsymbol{\Omega})$) in the $\boldsymbol{\beta}$-space by a polynomial number of variables and constraints. As compared in Table I, it is generally the most compact formulation among the ones that are the tightest and applicable for both planar and non-planar graphs. Now the models listed in Table I actually cover a wide spectrum of tightness and compactness levels for the formulation of spanning tree constraints. Researchers and practitioners may choose an appropriate model based on their preferences and needs, etc.

Next, the tightness and compactness of constraints (1)-(2), rather than solely constraint (1), are further discussed briefly.

Proposition 1

If the LP relaxation of the explicit formulation of constraint (1) defines the spanning tree polytope in the $\boldsymbol{\beta}$-space, the LP relaxation of the proposed radiality constraints (1)-(2) defines the spanning forest polytope in the $\boldsymbol{\alpha}$-space.

Proof: Assume that the LP relaxation of constraints (1)-(2) has a fractional vertex $(\boldsymbol{\alpha^{*}},\boldsymbol{\beta^{*}})$ in the $(\boldsymbol{\alpha},\boldsymbol{\beta})$-space. Further assume that $\boldsymbol{\beta^{*}}$ has fractional entries. Thus, $\boldsymbol{\beta^{*}}$ can be represented by an integer vertex $\boldsymbol{\beta^{m}}$ and another point $\boldsymbol{\beta^{n}}$ of the spanning tree polytope: $\boldsymbol{\beta^{*}}=\xi\cdot\boldsymbol{\beta^{m}}+(1-\xi)\cdot\boldsymbol{\beta^{n}}$ with $0<\xi<1$. For $(i,j)$ with $\beta_{ij}^{m}=0$ or $\beta_{ij}^{m}=1$, we have $\beta_{ij}^{*}=(1-\xi)\cdot\beta_{ij}^{n}$ and $\beta_{ij}^{*}=\xi+(1-\xi)\cdot\beta_{ij}^{n}$, respectively. As $0\leq\beta_{ij}^{n}\leq 1$ and $0\leq\alpha_{ij}^{*}$ $\leq\beta_{ij}^{*}\leq 1$, we have $0\leq\frac{\alpha_{ij}^{*}}{1-\xi}\leq\beta_{ij}^{n}$ for $(i,j)$ with $\beta_{ij}^{m}=0$, and have $0\leq\frac{\alpha_{ij}^{*}-\xi\cdot\chi_{ij}}{1-\xi}\leq\beta_{ij}^{n}$ with $\chi_{ij}\in[\frac{\alpha_{ij}^{*}-\beta_{ij}^{*}+\xi}{\xi},\frac{\alpha_{ij}^{*}}{\xi}]\cap[0,1]$ for $(i,j)$ with $\beta_{ij}^{m}=1$. Then, we can construct $\boldsymbol{\alpha^{m}}$ and $\boldsymbol{\alpha^{n}}$ to have $\alpha_{ij}^{m}=0$ and $\alpha_{ij}^{n}=\frac{\alpha_{ij}^{*}}{1-\xi}$ for $(i,j)$ with $\beta_{ij}^{m}=0$, and have $\alpha_{ij}^{m}=\chi_{ij}$ and $\alpha_{ij}^{n}=\frac{\alpha_{ij}^{*}-\xi\cdot\chi_{ij}}{1-\xi}$ for $(i,j)$ with $\beta_{ij}^{m}=1$, so that $(\boldsymbol{\alpha^{*}},\boldsymbol{\beta^{*}})=\xi\cdot(\boldsymbol{\alpha^{m}},\boldsymbol{\beta^{m}})+(1-\xi)\cdot(\boldsymbol{\alpha^{n}},\boldsymbol{\beta^{n}})$. This contradicts the assumption of $(\boldsymbol{\alpha^{*}},\boldsymbol{\beta^{*}})$ being a vertex. The case with $\boldsymbol{\beta^{*}}$ as a vertex of $\mathit{conv}(\boldsymbol{\Omega})$ can be analyzed similarly. Thus, vertices of the projection of constraints (1)-(2)’s LP relaxation into $\boldsymbol{\alpha}$-space are 0-1 incidence vectors of spanning forests of the DS. $\square$

Proposition 1 implies that the tightness and compactness features of constraints (1)-(2) essentially follow the explicit formulation of constraint (1). Thus, the proposed radiality constraints (1)-(2) also cover a wide spectrum of tightness and compactness levels for the formulation of spanning forest polytope. An appropriate model can be selected based on one’s preferences, etc.

Specifically, to exploit the aforementioned properties of the formulations, the DS is assumed to have only one substation node. For a DS with multiple substation nodes, one can merge them into one node in modeling constraints (1)-(2), but still treat them as separate nodes in modeling DS operational constraints.

IV Application to Resilient MG Formation

The proposed radiality model (1)-(2) can be applied in different optimization problems involving DSs and/or MGs with flexible boundaries [9, 10, 11], etc. Here, we apply it to regulate the DS topology in resilient MG formation to verify its advantages.

MG formation has recently been extensively studied in [13] [18] [19] [24] [25], etc. It is to reconfigure the DS into multiple MGs energized by DGs, so as to restore critical loads, etc. The literature currently has two major types of optimization models for this problem [18] [19]. Using the proposed radiality model, a new formulation for resilient MG formation is constructed:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Notations are listed in Table III. The objective function (10) maximizes the weighted sum of restored loads. Constraints (1)-(2) are used to ensure radiality and enable topological flexibilities, etc. Constraint (1) is formulated using the methods in Section III. Equations (11)-(12) enforce real and reactive power balance, respectively. Equation (13) imposes zero power output for nodes without power sources. Constraint (14) indicates real and reactive power capacities of substations or DGs. Constraint (15) represents the DistFlow model with the much smaller quadratic terms ignored [1] [26]. It is relaxed for open branches. Constraint (16) expresses the voltage magnitude limits. Constraint (17) is convex though non-linear, and can be linearized by the technique in [3], etc. It limits the apparent power on a branch by its capacity. Equation (18) restricts the connection status of faulted open or faulted closed branches. Equation (19) prohibits picking up the loads with faulted open switches, and enforces picking up the loads with faulted closed switches as long as their located nodes are energized. Equation (20) specifies that substation and DG nodes are energized. Constraint (21) derives the energization status of other nodes by examining if they are connected to an energized node. The non-linear and non-convex terms, i.e., $\varepsilon_{j}\cdot\alpha_{ij}$ and $\varepsilon_{j}\cdot\alpha_{ji}$, can be equivalently linearized and convexified by the McCormick envelopes [27].

Note that the areas with surviving access to the main grid power via substation nodes are also considered. Operational constraints may prohibit those areas to be fully restored by the substations, and forming MGs with DGs can help achieve better restoration of them [14]. For statement simplicity, a sub-grid powered by a substation is also counted as a MG here.

Table IV compares the proposed model with the two existing MG formation models in the literature. The summaries in the table are self-explanatory. Brief explanations are given as below:

1. The model in [18] is designed for radial systems. It is extended in [25] to deal with meshed systems. Still, the extended model does not eliminate loops. Thus, it may be used for the bulk power system, but does not apply to meshed DSs, which have to open some branches to be operated in radial topologies. Both model in [19] and our model can handle meshed DSs, as radiality constraints are included to avoid loops when reconfiguring the network. Single-commodity flow-based radiality constraints are used in[19], while our model adopts the proposed radiality constraints to fully enable topological flexibilities, etc.

2. Both models in [18] and [19] allocate one DG to each MG. Therefore, the number of MGs is actually fixed and equal to the number of DGs and/or other power sources. As for our proposed model, the number of DGs designated to different MGs, and the resulting number of MGs, are flexible. Thus, it can also be used for dynamic MG formation, which involves adaptive merge and/or separation of MGs when damaged parts of the DS are sequentially repaired, etc. Such flexibilities introduce many benefits. For example, larger MGs can be formed to better match DGs with different-sized loads, so as to enhance capacity utilization rates of DGs and the restoration of critical loads.

3. Both models in [18] [19] energize all nodes included in $\boldsymbol{N}$. Thus, their models cannot consider the nodes in load islands without power sources and isolated by faulted open branches. Discrete/non-dispatchable loads at nodes $\boldsymbol{N_{c}}$ are also forced to be picked up in their models. Our model can consider such load islands, and can be easily extended to optimize the allocation of mobile power sources in such islands and their merge with other islands after faulted open branches are repaired. Our model can also intentionally form unenergized islands in optimization, enabling more flexible pick-up of the loads at nodes $\boldsymbol{N_{c}}$. This flexibility can be critical, as some DSs may have a large portion of loads not equipped with switches, causing $\boldsymbol{N_{c}}$ to be a large set. In [13], with additional binary variables, the model in [19] is modified to permit more flexible allocation of DGs into MGs. However, it does not fully enable topological and some related flexibilities. For example, it still requires energizing all nodes.

In general, as indicated in Table IV, the proposed MG formation model is more adaptive and allows more flexibilities.

V Case Studies

In this section, the proposed resilient MG formation model is demonstrated on two test systems. We use a computer with an Intel i5-4278U processor and 8GB of memory. Involved MIP problems are solved by Gurobi 7.5.2 with the default settings.

V-A IEEE 33-Node Test System [1]

Fig. 4 depicts an illustrative case based on a scenario with 27

faults in total. The model in [18] does not apply to meshed DSs. Thus, as indicated in Fig. 4(a), it does not consider the normally open branches. It forms 6 MGs energized by the 6 DGs, respectively. As the model in [19] can handle meshed DSs, it gets better results via further reconfiguration involving those normally open switches, etc. Fig. 4(b) shows that it restores more loads by forming 6 larger MGs, which is essentially a 6-tree. By contrast, as in Fig. 4(c), our proposed model forms a 7-tree, i.e., 4 MGs and 3 load islands. Specifically, MG 4 and MG 6 in Fig. 4(b) are merged into a single MG in Fig. 4(c), so that the loads at nodes 31 and 33 can also be restored. MG 2 in Fig. 4(b) is separated into a MG and a load island in Fig. 4(c), so that node 23 is not energized and its load is not forced to be picked up. For space limit, we do not detail on other differences. The models in [18] [19] and our proposed model have DG capacity utilization rates of 55.6%, 67.4% and 95.4%, respectively. They restore 1500kW, 1820kW and 2575kW loads, respectively. Our model achieves more coordinated matching among the different-sized DGs and loads, and thus attains better service restoration. Generally, as the proposed radiality constraints can fully enable topological and some other related flexibilities of the DS, our MG formation model has extended feasibility and enhanced optimality.

To establish superiority of the proposed radiality constraints and MG formation model, we further run 10000 cases based on randomly generated scenarios of DS faults. Table V shows that our model has far less infeasible cases. In many scenarios, it finds a feasible operating point of the DS, while the models in [18] [19] return infeasibility. Such results verify that our model has extended feasibility as the proposed radiality constraints fully enable topological and some related flexibilities. Besides, while its search space is the largest, its average computation time is the shortest. Thus, its enlarged feasible set may possess more computationally tractable characteritics. As in Table VI, our model has the highest average/median/minimum restored loads and the smallest standard deviation. The same maximum value corresponds to the cases with all loads restored. On average, the proposed model restores 15.0% and 8.4% more loads than models in [18] and [19], respectively. Actually, our model performs equally well or better in all cases. The histograms in Fig. 5(Left) depict the outperformance. Such results validate the enhanced optimality of the proposed model. Here, one of the major reasons for its superiority is more coordinated matching among the DGs and loads. Table VII indicates that the average DG capacity utilization rate of our model is much higher than those of the models in [18] [19]. The smaller minimum and larger standard deviation are due to the cases with the substation contributing much power injection for service restoration.

The single-commodity flow-based model is used to formulate constraint (1) in all previous cases. Here, for comparison, the directed multicommodity flow-based model is used instead. Consequently, the proposed radiality constraints and MG formation model become tighter, though less compact. The revised model is run on the same 10000 scenarios. As indicated in Table VIII, $N_{m}$ is much smaller than $N_{s}$, both on average and in 95.8% of the cases. The histogram in Fig. 5(Right) also details the extra explored nodes in the B&C search tree of the less tight model. Although $T_{s}$ is slightly shorter than $T_{m}$ on average, $T_{m}$ is much shorter than $T_{s}$ in 10.9% cases. Specifically, in those cases, $T_{s}$ and $T_{m}$ are 1.92 s and 0.69 s on average, respectively; $N_{s}$ and $N_{m}$ are 815 and 186 on average, respectively. In general, the computation time of the revised tighter model is more consistent. It also reduces the computation time in many cases that require the less tight model to explore much more nodes in the B&C search tree. That is, the proposed radiality constraints can satisfy the need for tighter formulations of DS reconfiguration-related optimization problems.

V-B IEEE 123-Node Test System [28]

Fig. 6 provides an illustrative case based on a 60-fault sce-

nario using the proposed MG formation model. It forms 7 MGs powered by 8 DGs, and a sub-grid powered by the substation. The DG capacity utilization rate is 75.0%, and 3730kW loads are restored. We again run 10000 cases on random scenarios of DS faults. Table IX-XII and Fig. 7 show results similar to the previous system. That is, the proposed radiality constraints and MG formation model’s superiority is also established on this larger system. For example, our model matches DGs and loads in a more coordinated manner, thus achieving better service restoration. For space limit, we do not go into the details. In short, our MG formation model has extended feasibility and enhanced optimality due to the proposed radiality constraints’ fully enabling topological and some related flexibilities of DSs.

It is worth mentioning that, as in Table XII, $T_{m}$ is averagely shorter than $T_{s}$ for this larger system. Specifically, in the cases with $T_{m}<T_{s}$, $T_{m}$ and $T_{s}$ are 1.26 s and 5.59 s on average, respectively; $N_{m}$ and $N_{s}$ are 89 and 937 on average, respectively. That is, the computation time of the tighter version of our MG formation model is not only shorter on average, but also more consistent. It greatly reduces the computation time especially for cases requiring the less tight model to explore much more B&C search tree nodes. Thus, the tighter version of the proposed radiality constraints can allow both more flexible DS operation and more efficient computation for reconfiguration-related optimization problems, e.g., MG formation here.

VI Conclusion

This work proposes a new method for formulating radiality constraints that fully enables topological and some other related flexibilities in DS reconfiguration optimization problems. It is specifically applied to resilient post-disaster MG formation to attain extended feasibility and enhanced optimality. As verified in case studies, compared to the existing MG formation models, our model based on the proposed radiality constraints achieves a higher resilience enhancement via more coordinated utilization of DS flexibilities, etc., and reduces the computational complexity. Future work includes exploring the effects of increased topological flexibilities on other DS performance metrics, etc.