# A Two-Level Simulation-Assisted Sequential Distribution System Restoration
Model With Frequency Dynamics Constraints

Qianzhi Zhang, Zixiao Ma, Yongli Zhu, and Zhaoyu Wang  This work was supported
in part by the U.S. Department of Energy Wind Energy Technologies Office under
Grant DE-EE00008956 (_Corresponding author: Zhaoyu Wang_). The authors are
with the Department of Electrical and Computer Engineering, Iowa State
University, Ames, IA 50011 USA (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>wzy@iastate.edu).

###### Abstract

This paper proposes a service restoration model for unbalanced distribution
systems and inverter-dominated microgrids (MGs), in which frequency dynamics
constraints are developed to optimize the amount of load restoration and
guarantee the dynamic performance of system frequency response during the
restoration process. After extreme events, the damaged distribution systems
can be sectionalized into several isolated MGs to restore critical loads and
tripped non-black start distributed generations (DGs) by black start DGs.
However, the high penetration of inverter-based DGs reduces the system
inertia, which results in low-inertia issues and large frequency fluctuation
during the restoration process. To address this challenge, we propose a two-
level simulation-assisted sequential service restoration model, which includes
a mixed integer linear programming (MILP)-based optimization model and a
transient simulation model. The proposed MILP model explicitly incorporates
the frequency response into constraints, by interfacing with transient
simulation of inverter-dominated MGs. Numerical results on a modified IEEE
123-bus system have validated that the frequency dynamic performance of the
proposed service restoration model are indeed improved.

###### Index Terms:

Frequency dynamics, service restoration, network reconfiguration, inverter-
dominated microgrids, simulation-based optimization.

## Nomenclature

Sets

$\Omega_{\rm BK}$

Set of bus blocks.

$\Omega_{\rm G}$

Set of generators.

$\Omega_{\rm BS}$

Set of generators with black start capability.

$\Omega_{\rm NBS}$

Set of generators without black start capability.

$\Omega_{\rm K}$

Set of distribution lines.

$\Omega_{\rm SW_{K}}$

Set of switchable lines.

$\Omega_{\rm NSW_{K}}$

Set of non-switchable lines.

$\Omega_{\rm L}$

Set of loads.

$\Omega_{\rm SW_{L}}$

Set of switchable loads.

$\Omega_{\rm NSW_{L}}$

Set of non-switchable loads.

$\Omega_{\phi}$

Set of phases.

Indices

$BK$

Index of bus block.

$k$

Index of line.

$i,j$

Index of bus.

$t$

Index of time instant.

$\phi$

Index of three-phase $\phi_{a},\phi_{b},\phi_{c}$.

Parameters

$a_{\phi}$

Approximate relative phase unbalance.

$D_{\rm P},D_{\rm Q}$

$P-\omega$ and $Q-V$ droop gains.

$f_{0}$

Nominal steady-state frequency.

$f^{\rm min}$

Minimum allowable frequency during the transient simulation.

$M$

Big-M number.

$P_{i}^{\rm G,M},Q_{i}^{\rm G,M}$

Active and reactive power output maximum limits of generator at bus $i$.

$P_{k}^{\rm K,M},Q_{k}^{\rm K,M}$

Active and reactive power flow maximum limits of line $k$.

$p_{k,\phi}$

Phase identifier of line $k$.

$R,L$

Aggregate resistance and inductance of connections from the inverter
terminal’s point review.

$\hat{R}_{k},\hat{X}_{k}$

Matrices of resistance and reactance of line $k$.

$T$

Length of rolling horizon.

$U_{i}^{\rm m},U_{i}^{\rm M}$

Minimum and maximum limit for squared nodal voltage magnitude of bus $i$.

$V_{\rm bus}$

Bus voltage.

$Z_{k},\hat{Z}_{k}$

Matrices of original impedance and equivalent impedance of line $k$.

$\alpha$

Hyper-parameter in frequency dynamics constraints.

$\Delta f^{\rm max}$

User-defined maximum allowable frequency drop limit.

$\Delta f^{\rm meas}$

Measured maximum transient frequency drop.

$w_{i}^{\rm L}$

Priority weight factor for load of bus $i$.

$\omega_{\rm c}$

Cut-off frequency of the low pass filter.

$\omega_{\rm set},V_{\rm set}$

Set points of frequency and voltage controllers.

$\omega_{0}$

Nominal angular frequency.

Variables

$f^{\rm nadir}$

Frequency nadir during the transient simulation.

$I_{\rm d},I_{\rm q}$

$dq$-axis current.

$P,Q$

Filtered terminal output active and reactive power.

$P^{\rm L},Q^{\rm L}$

Restored active and reactive loads.

$P_{i,\phi,t}^{\rm G}$

Three-phase active power output of generator at bus $i$, phase $\phi$, time
$t$.

$P_{i,t}^{\rm G,MLS}$

Maximum load step at bus $i$, time $t$.

$P_{k,\phi,t}^{\rm K}$

Three-phase active power flow of line $k$, phase $\phi$, time $t$.

$P_{i,\phi,t}^{\rm L}$

Restored active load at bus $i$, phase $\phi$, time $t$.

$Q_{i,\phi,t}^{\rm G}$

Three-phase reactive power output of generator at bus $i$, phase $\phi$, time
$t$.

$Q_{k,\phi,t}^{\rm K}$

Three-phase reactive power flow of line $k$, phase $\phi$, time $t$.

$U_{i,\phi,t}$

Squared of three-phase voltage magnitude.

$V$

Output voltage of the inverter.

$x_{i,t}^{\rm B}$

Binary energizing status of bus, if $x_{i,t}^{\rm B}=1$ then the bus $i$ is
energized at time $t$.

$x_{B,t}^{\rm BK}$

Binary energizing status of bus block, if $x_{B,t}^{\rm BK}=1$ then the bus
block $B$ is energized at time $t$.

$x_{i,t}^{\rm G}$

Binary switch on/off status of grid-following generator, if $x_{i,t}^{\rm
G}=1$ then the grid-following generator at bus $i$ is switched on at time $t$.

$x_{k,t}^{\rm K}$

Binary connection status of line, if $x_{k,t}^{\rm K}=1$ then the line $k$ is
connected at time $t$.

$x_{i,t}^{\rm L}$

Binary restoration status of load, if $x_{i,t}^{\rm L}=1$ then the load $i$ is
restored at time $t$.

$\Delta P_{i,t-1}^{G,MLS}$

Change of the maximum load step.

$\theta$

Output phase angle of the inverter.

$\omega$

Output angular frequency of the inverter.

##  I Introduction

EXTREME events can cause severe damages to power distribution systems [1],
e.g. substation disconnection, line outage, generator tripping, load shedding,
and consequently large-scale system blackouts [2]. During the network and
service restoration, in order to isolate faults and restore critical loads, a
distribution system can be sectionalized into several isolated microgirds
(MGs) [3]. Through the MG formation, buses, lines and loads in outage areas
can be locally energized by distributed generations (DGs), where more outage
areas could be restored and the number of switching operations could be
minimized [4, 5, 6, 7, 8, 9]. In [4], the self-healing mode of MGs is
considered to provide reliable power supply for critical loads and restore the
outage areas. In [5], a networked MGs-aided approach is developed for service
restoration, which considers both dispatchable and non-dispatchable DGs. In
[6] and [7], the service restoration problem is formulated as a mixed integer
linear programming (MILP) to maximize the critical loads to be restored while
satisfying constraints for MG formation and remotely controlled devices. In
[8], the formation of adaptive multiple MGs is developed as part of the
critical service restoration strategy. In [9], a sequential service
restoration framework is proposed to generate restoration solutions for MGs in
the event of large-scale power outages. However, the previous methods mainly
use the conventional synchronous generators as the black start units, and only
consider steady-state constraints in the service restoration models, which
have limitations in the following aspects:

(1) An inverter-dominated MG can have low-inertia: With the increasing
penetration of inverter-based DGs (IBDGs) in distribution systems, such as
distributed wind and photovoltaics (PVs) generations, the system inertia
becomes lower [10, 11]. When sudden changes happen, such as DG output
changing, load reconnecting, and line switching, the dynamic frequency
performance of such low-inertia distribution systems can deteriorate [12].
This issue becomes even worse when restoring low-inertia inverter-dominated
MGs. Without considering frequency dynamics constraints, the load and service
restoration decisions may not be implemented in practice.

(2) Frequency responses need to be considered: Previous studies [13, 14, 15,
16] have considered the impact of disturbances on frequency responses in the
service restoration problem using different approaches. In [13], the amount of
load restored by DGs is limited by a fixed frequency response rate and maximum
allowable frequency deviation. However, because the frequency response rate is
pre-determined in an off-line manner, the impacts of significant load
restoration, topology change, and load variations may not be fully captured by
the off-line model. In [14], the stability and security constraints are
incorporated into the restoration model. However, this model has to be solved
by meta-heuristic methods due to the nonlinearity of the stability
constraints, which may lead to large optimality gaps. In [15], even though the
transient simulation results of voltage and frequency are considered to
evaluate the potential MG restoration paths in an online manner, it adopts a
relatively complicated four-stage procedure to obtain the optimal restoration
path. In [16], a control strategy of real-time frequency regulation for
network reconfiguration is developed, nonetheless, it is not co-optimized with
the switch operations.

(3) Grid-forming IBDGs need to be considered: In previous studies on optimal
service restoration, IBDGs are usually modeled as grid-following sources
(i.e., PQ sources) to simply supply active and reactive power based on the
control commands. However, during the service restoration after a network
blackout and loss of connection to the upstream feeder, a grid-forming IBDG
will be needed to setup voltage and frequency references for the blackout
network [17]. During outages, the grid-following IBDGs will be switched off.
After outages, the grid-forming IBDGs have the black start capability, which
can restore loads after the faults are isolated. Because IBDGs are connected
with power electronics converters and have no rotating mass, there is no
conventional concept of “inertia” for IBDGs. Thus, control techniques such as
droop control [18, 19] and virtual synchronous generator (VSG) control [20,
21] are usually adopted to emulate the inertia property in IBDGs.

To alleviate the frequency fluctuations caused by service restoration, we
establish a MILP-based optimization model with frequency dynamics constraints
for sequential service restoration to generate sequential actions for remotely
controlled switches, restoration status for buses, lines, loads, operation
actions for grid-forming and grid-following IBDGs, which interacts with the
transient simulation of inverter-dominated MGs. Inspired by recent advances in
simulation-assisted methods [15, 22] and to incorporate the frequency dynamics
constraints explicitly in the optimization formulation, we associate the
frequency nadir of the transient simulation with respect to the maximum load
that a MG can restore. Although some previous works have considered the
transient simulation as well in finding the optimal restoration solution, they
either adopts a heuristic framework, or merely using the transient simulation
to validate the feasibility of the obtained restoration solution after solving
an optimization problem. By contrast, the proposed two-level simulation-
assisted restoration model directly incorporates the transient simulation
module on top of a strict MILP optimization problem via explicit constraints,
thus its solving process is more tractable and straightforward.

The main contribution of this paper is two-folded:

  * •

We develop a two-level simulation-assisted sequential service restoration
model within a rolling horizon framework, which combines a MILP-based
optimization level of service restoration and a transient simulation level of
inverter-dominated MGs.

  * •

Frequency dynamics constraints are developed and explicitly incorporated in
the optimization model, to associate the simulated frequency responses with
the decision variables of maximum load step at each stage. These constraints
help restrict the system frequency drop during the transient periods of
restoration. Thus, the generated restoration solution can be more secure and
practical.

The reminder of the paper is organized as follows: Section II presents the
overall framework of the proposed service restoration model. Section III
introduces frequency dynamics constrained MILP-based sequential service
restoration. Section IV describes transient simulation of inverter-dominated
MGs. Numerical results and conclusions are given in Section V and Section VI,
respectively.

##  II Overview of the Proposed service restoration Model

The general framework of the proposed two-level simulation-assisted service
restoration is shown in Fig. 1, including an optimization level of MILP-based
sequential service restoration model and a transient simulation level of
$7$th-order electromagnetic inverter-dominated MG dynamic model. After
outages, the fault-affected areas of the distribution system will be isolated.
Consequently, each isolated sub-network can be considered as a MG [23], which
can be formed by the voltage and frequency supports from the grid-forming
IBDGs, and active and reactive power supplies from the grid-following IBDGs.
In the proposed optimization level, each MG will determine its restoration
solutions, including optimal service restoration status of loads, optimal
operation of remotely controlled switches and optimal active and reactive
power dispatches of IBDGs. To prevent large frequency fluctuation due to a
large load restoration, the maximum restorable load for a given period is
limited by the proposed frequency dynamics constraints. In this way, the whole
restoration process is divided into multiple stages. As shown in Fig. 1, the
information exchanged between the optimization level and the simulation level
are the restoration solution (obtained from optimization) and MG system
frequency nadir value (obtained from transient simulation): at each
restoration stage, the optimization level will obtain and send the optimal
restoration solution to the simulation level; then, after receiving the
restoration solution, the simulation level will begin to run transient
simulation by the proposed dynamic model of each inverter-dominated MG, and
send the frequency nadir value to the optimization level for next restoration
stage.

Figure 1: The overall framework of the proposed service restoration model with
optimization level and simulation level.

To accurately reflect the dynamic frequency-supporting capacities of grid-
forming IBDGs during the service restoration process, a rolling-horizon
framework is implemented in the proposed service restoration model, as shown
in Fig. 2. More specifically, we repeatedly run the MILP-based sequential
service restoration model by incorporating the network configuration from the
preceding stage as the initial condition, and then feedback the frequency
nadir value from the transient simulation to the frequency dynamics
constraints. For each stage: (1) the horizon length will be fixed; (2) then
only the restoration solution of first horizon of the current stage is
retained and transferred to the simulation level, while the remaining horizons
are discarded; (3) this process will keep going until the maximum restored
load is reached in each MG. More details about the principles of rolling
horizon can be found in [24].

Figure 2: Implementation of rolling-horizon in the proposed restoration model.

##  III Frequency Dynamics Constrained Service Restoration

This section presents the mathematical formulation for coordinating remotely
controlled switches, grid-forming and grid-following IBDGs, and the sequential
restoration status of buses, lines and loads. Here, we consider a unbalanced
three-phase radial distribution system. The three-phase
$\phi_{a},\phi_{b},\phi_{c}$ are simplified as $\phi$. Define the set
$\Omega_{\rm L}=\Omega_{\rm SW_{L}}\cup\Omega_{\rm NSW_{L}}$, where
$\Omega_{\rm SW_{L}}$ and $\Omega_{\rm NSW_{L}}$ represent the set of
switchable load and the set of non-switchable loads, respectively. Define the
set $\Omega_{\rm G}=\Omega_{\rm BS}\cup\Omega_{\rm NBS}$, where $\Omega_{\rm
BS}$ and $\Omega_{\rm NBS}$ represent the set of grid-forming IBDGs with black
start capability and the set of grid-following IBDGs without black start
capability, respectively. Define the set $\Omega_{\rm K}=\Omega_{\rm
SW_{K}}\cup\Omega_{\rm NSW_{K}}$, where $\Omega_{\rm SW}$ and $\Omega_{\rm
NSW}$ represent the set of switchable lines and the set of non-switchable
lines, respectively. Define $\Omega_{\rm BK}$ as the set of bus blocks, where
bus block [9] is a group of buses interconnected by non-switchable lines and
those bus blocks are interconnected by switchable lines. It is assumed that
bus block can be energized by grid-forming IBDGs. By forcing the related
binary variables of faulted lines to be zeros, each faulted area remains
isolated during the restoration process.

###  III-A MILP-based Sequential Service Restoration Formulation

The objective function (1) aims to maximize the total restored loads with
priority factor $w_{i}^{L}$ over a rolling horizon $[t,t+T]$ as shown below:

$\max\sum_{t\in[t,t+T]}\sum_{i\in\Omega_{L}}\sum_{\phi\in\Omega_{\phi}}(w_{i}^{\rm
L}x_{i,t}^{\rm L}P_{i,\phi,t}^{\rm L})$ (1)

where $P_{i,\phi,t}^{\rm L}$ and $x_{i,t}^{\rm L}$ are the restored load and
restoration status of load at $t$. If the load demand $P_{i,\phi,t}^{\rm L}$
is restored, then $x_{i,t}^{\rm L}=1$. $T$ is horizon length in the rolling
horizon optimization problem. In this work, the amount of restored load is
also bounded by frequency dynamics constraints with respect to frequency
response and maximum load step. More details of frequency dynamics constraints
are discussed in Section III-B.

Constraints (2)-(11) are defined by the unbalanced three-phase version of
linearized DistFlow model [25, 26] in each formed MG during the service
restoration process. Constraints (2) and (3) are the nodal active and reactive
power balance constraints, where $P_{k,\phi,t}^{\rm K}$ and $Q_{k,\phi,t}^{\rm
K}$ are the active and reactive power flows along line $k$, and $P^{\rm
G}_{i,\phi,t}$ and $Q^{\rm G}_{i,\phi,t}$ are the power outputs of the
generators. Constraints (4) and (5) represent the active and reactive power
limits of the lines, where the limits ($P_{k}^{\rm K,M}$ and $Q_{k}^{\rm
K,M}$) are multiplied by the line status binary variable $x_{k,t}^{\rm K}$.
Therefore, if a line is disconnected or damaged $x_{k,t}^{\rm K}=0$, then
constraints (4) and (5) will be relaxed, which means that power cannot flow
through this line. In the proposed model, there are two types of IBDGs, grid-
forming IBDGs with black start capability and grid-following IBDGs without
black start capability. On the one side, the grid-forming IBDGs can provide
voltage and frequency references in the MG during the restoration process,
which can energize the bus and restore the part of the network that is not
damaged if the fault is isolated. Therefore, the grid-forming IBDGs are
considered to be connected to the network at the beginning of restoration. On
the other side, the grid-following IBDGs are switched off at the beginning of
restoration. If the grid-following IBDGs are connected to an energized bus
during the restoration process, then they can be switched on to supply active
and reactive powers. In constraints (6) and (7), the active and reactive power
outputs of the grid-forming IBDGs are limited by the maximum active and
reactive capacities $P_{i}^{\rm G,M}$ and $Q_{i}^{\rm G,M}$, respectively.
Constraints (8) and (9) limit the active and reactive outputs of the grid-
following IBDGs. Note that the constraints (8) and (9) of grid-following IBDGs
are multiplied by binary variable $x_{i,t}^{\rm G}$. Consequently, if one
grid-following IBDG is not energized ($x_{i,t}^{\rm G}=0$) during the
restoration process, then constraints (8) and (9) of this grid-following IBDG
will be relaxed.

$\displaystyle\sum_{k\in\Omega_{\rm K}(i,.)}P_{k,\phi,t}^{\rm
K}-\sum_{k\in\Omega_{\rm K}(.,i)}P_{k,\phi,t}^{\rm K}=P^{\rm
G}_{i,\phi,t}-x_{i,t}^{\rm L}P_{i,\phi,t}^{\rm L},\forall i,\phi,t$ (2)
$\displaystyle\sum_{k\in\Omega_{\rm K}(i,.)}Q_{k,\phi,t}^{\rm
K}-\sum_{k\in\Omega_{\rm K}(.,i)}Q_{k,\phi,t}^{\rm K}$ $\displaystyle=Q^{\rm
G}_{i,\phi,t}-x_{i,t}^{\rm L}Q_{i,\phi,t}^{\rm L},\forall i,\phi,t$ (3)
$-x_{k,t}^{\rm K}P_{k}^{\rm K,M}\leq P_{k,\phi,t}^{\rm K}\leq x_{k,t}^{\rm
K}P_{k}^{\rm K,M},\forall k\in\Omega_{\rm K},\phi,t$ (4) $-x_{k,t}^{\rm
K}Q_{k}^{\rm K,M}\leq Q_{k,\phi,t}^{\rm K}\leq x_{k,t}^{\rm K}Q_{k}^{\rm
K,M},\forall k\in\Omega_{\rm K},\phi,t$ (5) $0\leq P_{i,\phi,t}^{\rm G}\leq
P_{i}^{\rm G,M},\forall i\in\Omega_{\rm BS},\phi,t$ (6) $0\leq
Q_{i,\phi,t}^{\rm G}\leq Q_{i}^{\rm G,M},\forall i\in\Omega_{\rm BS},\phi,t$
(7) $0\leq P_{i,\phi,t}^{\rm G}\leq x_{i,t}^{\rm G}P_{i}^{\rm G,M},\forall
i\in\Omega_{\rm NBS},\phi,t$ (8) $0\leq Q_{i,\phi,t}^{\rm G}\leq x_{i,t}^{\rm
G}Q_{i}^{\rm G,M},\forall i\in\Omega_{\rm NBS},\phi,t$ (9)

Constraints (10) and (11) calculate the voltage difference along line $k$
between bus $i$ and bus $j$, where $U_{i,\phi,t}$ is the square of voltage
magnitude of bus $i$. We use the big-M method [9] to relax constraints (10)
and (11), if lines are damaged or disconnected, then $x_{k,t}^{\rm K}=0$. The
$p_{k,\phi}$ represents the phase identifier for phase $\phi$ of line $k$. For
example, if line $k$ is a single-phase line on phase a, then
$p_{k,\phi_{a}}=1$, $p_{k,\phi_{b}}=0$ and $p_{k,\phi_{c}}=0$. Constraint (12)
guarantees that the voltage is limited within a specified region [$U^{\rm
m}_{i}$,$U^{\rm M}_{i}$], and will be set to 0 if the bus is in an outage area
$x_{i,t}^{\rm B}=0$.

$\begin{split}U_{i,\phi,t}-U_{j,\phi,t}\geq&2(\hat{R}_{k}P_{k,\phi,t}^{\rm
K}+\hat{X}_{k}Q_{k,\phi,t}^{\rm K})\\\ &+(x_{k,t}^{\rm
K}+p_{k,\phi}-2)M,\forall k,ij\in\Omega_{\rm K},\phi,t\end{split}$ (10)
$\begin{split}U_{i,\phi,t}-U_{j,\phi,t}\leq&2(\hat{R}_{k}P_{k,\phi,t}^{\rm
K}+\hat{X}_{k}Q_{k,\phi,t}^{\rm K})\\\ &+(2-x_{k,t}^{\rm
K}-p_{k,\phi})M,\forall k,ij\in\Omega_{\rm K},\phi,t\end{split}$ (11)
$x_{i,t}^{\rm B}U^{\rm m}_{i}\leq U_{i,\phi,t}\leq x_{i,t}^{\rm B}U^{\rm
M}_{i},\forall i,\phi,t$ (12)

where $\hat{R}_{k}$ and $\hat{X}_{k}$ are the unbalanced three-phase
resistance matrix and reactance matrix of line $k$. To model the unbalanced
three-phase network, we assume that the distribution network is not too
severely unbalanced and operates around the nominal voltage, then the relative
phase unbalance can be approximated as $a_{\phi}=[1,{e}^{-{\bf
i}2\pi/3},{e}^{{\bf i}2\pi/3}]^{T}$ [25]. Therefore, the equivalent unbalanced
three-phase system line impedance matrix $\hat{Z}_{k}$ can be calculated based
on the original line impedance matrix $Z_{k}$ and $a_{\phi}$ in (13).
$\hat{R}_{k}$ and $\hat{X}_{k}$ are the real and imaginary parts of
$\hat{Z}_{k}$, as shown in (14). Note that the loads and IBDGs are also
modelled in a three-phase form. More details about the model of unbalance
three-phase distribution system can be found in [26].

$\hat{Z}_{k}=a_{\phi}a_{\phi}^{H}\odot Z_{k}$ (13)
$\hat{R}_{k}=real(\hat{Z}_{k}),\hat{X}_{k}=imag(\hat{Z}_{k})$ (14)

Constraints (15)-(22) ensure the physical connections among buses, lines,
IBDGs and loads during restoration process. In constraint (15), the grid-
following IBDGs will be switched on $x_{i,t}^{\rm G}=1$, if the connected bus
is energized $x_{i,t}^{\rm B}=1$; otherwise, $x_{i,t}^{\rm G}=0$. Constraint
(16) implies a switchable line can only be energized when both end buses are
energized. Constraint (17) presents that a non-switchable line can be
energized once one of two end buses is energized. Constraint (18) ensures that
a switchable load can be energized $x_{i,t}^{\rm L}=1$, if the connected bus
is energized $x_{i,t}^{\rm B}=1$; otherwise, $x_{i,t}^{\rm L}=0$. Constraint
(19) allows that a non-switchable load can be immediately energized once the
connected bus is energized. Constraints (20)-(22) ensure that the grid-
following IBDGs, switchable lines and loads cannot be tripped again, if they
have been energized at the previous time $t-1$.

$x_{i,t}^{\rm G}\leq x_{i,t}^{\rm B},\forall i\in\Omega_{\rm NBS},t$ (15)
$x_{k,t}^{\rm K}\leq x_{i,t}^{\rm B},x_{k,t}^{\rm K}\leq x_{j,t}^{\rm
B},\forall k,ij\in\Omega_{\rm SW_{K}},t$ (16) $x_{k,t}^{\rm K}=x_{i,t}^{\rm
B},x_{k,t}^{\rm K}=x_{j,t}^{\rm B},\forall k,ij\in\Omega_{\rm NSW_{K}},t$ (17)
$x_{i,t}^{\rm L}\leq x_{i,t}^{\rm B},\forall i\in\Omega_{\rm SW_{L}},t$ (18)
$x_{i,t}^{\rm L}=x_{i,t}^{\rm B},\forall i\in\Omega_{\rm NSW_{L}},t$ (19)
$x_{i,t}^{\rm G}-x_{i,t-1}^{\rm G}\geq 0,\forall i\in\Omega_{\rm NBS},t$ (20)
$x_{k,t}^{\rm K}-x_{k,t-1}^{\rm K}\geq 0,\forall k\in\Omega_{\rm SW_{k}},t$
(21) $x_{i,t}^{\rm L}-x_{i,t-1}^{\rm L}\geq 0,\forall i\in\Omega_{\rm
SW_{L}},t$ (22)

Constraints (23)-(25) ensure that each formed MG remains isolated from each
other and each MG can maintain a tree topology during the restoration process.
Constraint (23) implies that if one bus $i$ is located in one bus block,
$i\in\Omega_{\rm BK}$, then the energization status of bus and the
corresponding bus block keep the same. Here $x_{B,t}^{\rm BK}$ represents the
energization status of bus block $BK$. To avoid forming loop topology,
constraint (24) guarantees that a switchable line cannot be closed at time $t$
if its both end bus blocks are already energized at previous time $t-1$. Note
that the DistFlow model is valid for radial distribution network, therefore,
loop topology is not considered in this work. If one bus block is not
energized at previous time $t-1$, then constraint (25) makes sure that this
bus block can only be energized at time $t$ by at most one of the connected
switchable lines. Constraints (26) and (27) ensure that each formed MG has a
reasonable restoration and energization sequence of switchable lines and bus
blocks. Constraints (26) implies that energized switchable lines can energize
the connected bus block. Constraints (27) requires that a switchable line can
only be energized at time $t$, if at least one of the connected bus block is
energized at previous time $t-1$.

$x_{i,t}^{\rm B}=x_{i,t}^{\rm BK},\forall i\in\Omega_{\rm BK},t$ (23)
$\begin{split}(x_{i,t}^{\rm BK}-x_{i,t-1}^{\rm BK})&+(x_{j,t}^{\rm
BK}-x_{j,t-1}^{\rm BK})\\\ &\geq x^{\rm K}_{k,t}-x^{\rm K}_{k,t-1},\forall
k,ij\in\Omega_{\rm SW_{K}},t\geq 2\end{split}$ (24)
$\begin{split}\sum_{ki,k\in\Omega_{i}}(x^{\rm K}_{ki,t}-x^{\rm
K}_{ki,t-1})&+\sum_{ij,j\in\Omega_{i}}(x^{\rm K}_{ij,t}-x^{\rm K}_{ij,t-1})\\\
&\leq 1+x_{i,t-1}^{\rm BK}M,\forall k,ij\in\Omega_{\rm SW_{K}},t\geq
2\end{split}$ (25) $x_{i,t-1}^{\rm BK}\leq\sum_{ki,k\in\Omega_{i}}(x^{\rm
K}_{ki,t})+\sum_{ij,j\in\Omega_{i}}(x^{\rm K}_{ij,t}),\forall
k,ij\in\Omega_{\rm SW_{K}},t\geq 2$ (26) $x^{\rm K}_{ij,t}\leq x_{i,t-1}^{\rm
BK}+x_{j,t-1}^{\rm BK},\forall ij\in\Omega_{\rm SW_{K}},t\geq 2$ (27)

###  III-B Simulation-based Frequency Dynamics Constraints

By considering the frequency dynamics of each isolated inverter-dominated MG
during the transitions of network reconfiguration and service restoration,
constraints (28) and (30) have been added here to avoid the potential large
frequency deviations caused by MG formation and oversized load restoration.
The variable of maximum load step $P_{i,t}^{G,MLS}$ has been applied in
constraint (28) to ensure that the restored load is limited by an upper bound
for each restoration stage, as follows:

$\begin{split}0\leq P_{i,t}^{\rm G,MLS}&\leq P_{i,t-1}^{\rm G,MLS}\\\
&+\alpha(\Delta f^{\rm max}-\Delta f^{\rm meas}),\forall i\in\Omega_{\rm
BS},t\geq 2\end{split}$ (28)

In constraint (28), the variable $P_{i,t}^{G,MLS}$ is restricted by three
items: a hyper-parameter $\alpha$ representing the virtual frequency-power
characteristic of IBDGs, a user-defined maximum allowable frequency drop limit
$\Delta f^{\rm max}$ and the measured maximum transient frequency drop from
the results of simulation level $\Delta f^{\rm meas}$. The hyper-parameter
$\alpha$ is used to curb the frequency nadir during transients from too low.
This can be shown by the following expressions:

$\displaystyle\alpha(\Delta f^{\rm max}-\Delta f^{\rm meas})$
$\displaystyle=\alpha(f_{0}-f^{\rm min}-(f_{0}-f^{\rm nadir}))$
$\displaystyle=\alpha(f^{\rm nadir}-f^{\rm min})$
$\displaystyle\triangleq\Delta P_{i,t-1}^{\rm G,MLS}$ (29)

where $f_{0}$ is the nominal steady-state frequency, e.g. 60Hz. $f^{\rm
nadir}$ is the lowest frequency reached during the transient simulation.
$f^{\rm min}$ is the minimum allowable frequency. $\Delta P_{i,t-1}^{G,MLS}$
is the incremental change of the maximum load step for the next step $t$
(estimated at step $t-1$). Finally, constraint (30) ensures the restored load
and frequency response of the IBDGs do not exceed the user-defined thresholds.

$\begin{split}-x_{i,t}^{\rm G}P_{i,t}^{\rm G,MLS}\leq P_{i,\phi,t}^{\rm
G}-P_{i,\phi,t-1}^{\rm G}\leq&x_{i,t}^{\rm G}P_{i,t}^{\rm G,MLS}\\\
&,i\in\Omega_{\rm BS},\phi,t\geq 2\end{split}$ (30)

Note that the generator ramp rate is not a constant number anymore as in
previous literature, but is varying with the value of $P_{i,t}^{\rm G,MLS}$
from (28) during the optimization process combining with transient simulation
information of frequency deviation. When $f^{\rm nadir}$ is approaching
$f^{\rm min}$, that implies a necessity to reduce the potential amount of
restored load in the next step. Thus the incremental change of maximum load
step $\Delta P_{i,t}^{\rm G,MLS}$ is reduced to reflect the above purpose.
During the restoration process, the restored load in each restoration stage is
determined by maximum load step and available DG power output through power
balance constraints (2), (3) and constraints (28), (30) in optimization level;
then, the frequency deviation in each restoration stage is determined by
restored load through transient model in simulation level, which is introduced
in the next section.

##  IV Transient Simulation of Inverter-Dominated MG Formation

In optimization level, our target is to maximize the amount of restored load
while satisfying a series of constraints. One of these constraints should be
frequency dynamics constraint which is derived from simulation level. However,
due to the different time scales and nonlinearity, the conventional dynamic
security constraints cannot be directly solved in optimization problem, such
as Lyapunov theory, LaSalle’s theorem and so on. Therefore, we need a
connection variable between the two levels.

For this purpose, we assume that the changes of typologies between each two
sequential stages can be represented by the change of restored loads $P^{\rm
L}$. The sudden load change of $P^{\rm L}$ results in a disturbance in MGs in
the time-scale of simulation level. During the transience to the new
equilibrium (operation point), the system states such as frequency will
deviate from their nominal values. Therefore, it is natural to estimate the
dynamic security margin with the allowed maximum range of deviations.

Figure 3: Diagram of studied MG control system.

Since the frequency of each inverter-dominated MG is mainly controlled by the
grid-forming IBDGs, we can approximate the maximum frequency deviation during
the transience by observing the dynamic response of the grid-forming IBDGs
under sudden load change. In this paper, the standard outer droop control
together with inner double-loop control structure is adopted for each IBDGs
unit. As shown in Fig. 3, the three-phase output voltage $V_{0,abc}$ and
current $I_{0,abc}$ are measured from the terminal bus of the inverter and
transformed into $dq$ axis firstly. Then, the filtered terminal output active
and reactive power $P$ and $Q$ are obtained by filtering the calculated power
measurements $P^{\rm meas}$ and $Q^{\rm meas}$ with cut-off frequency
$\omega_{\rm c}$. Finally, the voltage and frequency references for the inner
control loop are calculated with droop controller. Since the references can be
accurately tracked by inner control loop with properly tuned PID parameters in
the much faster time-scale, the output voltage $V$ and frequency $\omega$ can
be considered equivalently as the references generated by the droop
controller. Thus, the inverter can be modelled effectively modelled by using
the terminal states and line states of the inverter [18, 19]. In this work,
the transient simulation is conducted with the detailed mathematical MG model
(31)–(37) adopted from [18], where the droop equations (34) and (35) are
replaced by the ones proposed in [19] to consider the restored loads.

$\displaystyle\dot{P}$ $\displaystyle=\omega_{\rm c}(V\cos{\theta}I_{\rm
d}+V\sin{\theta}I_{\rm q}-P),$ (31) $\displaystyle\dot{Q}$
$\displaystyle=\omega_{\rm c}(V\sin{\theta}I_{\rm d}-V\cos{\theta}I_{\rm
q}-Q),$ (32) $\displaystyle\dot{\theta}$ $\displaystyle=\omega-\omega_{0},$
(33) $\displaystyle\dot{\omega}$ $\displaystyle=\omega_{\rm c}(\omega_{\rm
set}-\omega+D_{\rm P}(P-P^{\rm L})),$ (34) $\displaystyle\dot{V}$
$\displaystyle=\omega_{\rm c}(V_{\rm set}-V+D_{\rm Q}(Q-Q^{\rm L})),$ (35)
$\displaystyle\dot{I}_{\rm d}$ $\displaystyle=(V\cos{\theta}-V_{\rm
bus}-RI_{\rm d})/L+\omega_{o}I_{\rm q},$ (36) $\displaystyle\dot{I}_{\rm q}$
$\displaystyle=(V\sin{\theta}-RI_{\rm q})/L-\omega_{o}I_{\rm d},$ (37)

where $\omega_{\rm set}$ and $V_{\rm set}$ are the set points of frequency and
voltage controllers, respectively; $\omega_{\rm c}$ is cut-off frequency;
$D_{\rm P}$ and $D_{\rm Q}$ are $P-\omega$ and $Q-V$ droop gains,
respectively; $P^{\rm L}$ and $Q^{\rm L}$ are the restored active and reactive
loads, respectively; $\theta$ is phase angle; $\omega$ is angular frequency in
$rad/s$; $\omega_{0}$ is a fixed angular frequency; $V_{\rm bus}$ is bus
voltage; $I_{\rm d}$ and $I_{\rm q}$ are $dq$-axis currents; $R$ and $L$ are
aggregate resistance and inductance of connections from the inverter
terminal’s point view, respectively. In (34), it can be observed that, the
equilibrium can be achieved when $\omega=\omega_{\rm set}$ and $P=P^{\rm L}$,
which means that the output frequency tracks the frequency reference when the
output power of the simulation level tracks the obtained restored load of the
optimization level.

Note that constraint (28) is the connection between the optimization level and
simulation level in our proposed two-level simulation-assisted restoration
model, which incorporates the frequency response of inverter-dominated MG from
the simulation level into the optimization level. The variable $P_{i,t}^{\rm
G,MLS}$ is restricted by frequency response in constraint (28). Meanwhile,
$P_{i,t}^{\rm G,MLS}$ also limits the IBDG power output in constraint (30). In
constraints (2) and (3), the power balance is met between restored load and
power supply of IBDGs. Therefore, we associate the frequency nadir of the
transient simulation with respect to the restored load by incorporating the
frequency dynamics constraints explicitly in the optimization level.

After the process of fault detection [27] and sub-grids isolation are
finished, the proposed service restoration model will begin to work. Each
isolated network will begin to form a MG depending on the location of the
nearest grid-forming IBDG with black start capability. The flowchart of the
proposed restoration method is shown in Fig. 4 and the interaction between the
proposed transient simulation and the established optimization problem of
service restoration is described as follows:

Figure 4: Flowchart of the proposed two-level simulation-assisted restoration
method.

(a) Solving the optimal service restoration problem: Given horizon length $T$
in each restoration stage, the MILP-based sequential service restoration
problem (1)–(28) and (30) is solved, and the restoration solution is obtained
for each formed MG.

(b) Transient simulation of inverter-dominated MGs: After receiving
restoration solutions of current stage from optimization level, the frequency
response is simulated by (31)–(37) and the frequency nadir is calculated for
each inverter-dominated MG.

(c) Check the progress of service restoration and stopping criteria: If the
maximum service restoration level is reached for all the MGs, then stop the
restoration process; otherwise, go back to (a) to generate the restoration
solution with newly obtained frequency responses of all MGs for next
restoration stage.

##  V Numerical Results

###  V-A Simulation Setup

A modified IEEE 123-bus test system [28] in Fig. 5 is used to test the
performance of the proposed frequency dynamics constrained service restoration
model. In Fig. 5, blue dotted line and blue dot stand for single-phase line
and bus, orange dashed line and orange dot stand for two-phase line and bus,
black line and black dot stand for three-phase line and bus, respectively. The
modified test system has been equipped with multiple remotely controlled
switches, as shown in Fig. 5. In Table I, the locations and capacities of
grid-following and grid-forming IBDGs are shown. Four line faults on lines
between substation and bus 1, bus 14 and bus 19, bus 14 and bus 54 and bus 62
and bus 70 are detected, as shown in red dotted lines of Fig. 5. They are
assumed to be persisting during the restoration process until the faulty areas
are cleared to maintain the radial topology and isolate the faulty areas.
Consequently, four MGs can be formed for service restoration with grid-forming
IBDGs and switches. For the sake of simplicity, we assume that the weight
factors for all loads are set to 1 during the restoration process. We
demonstrate the effectiveness of our proposed service restoration model
through numerical evaluations on the following experiments: (i) Comparison
between a base case (i.e., without the proposed frequency dynamics
constraints) and the case with the proposed restoration model. (ii) Cases with
the proposed restoration model under different values of hyper-parameters. All
the case studies are implemented using a PC with Intel Core i7-4790 3.6 GHz
CPU and 16 GB RAM hardware. The simulations are performed in MATLAB R2019b,
which integrates YALMIP Toolbox with IBM ILOG CPLEX 12.9 solver and ordinary
differential equation solver.

Figure 5: Modified IEEE 123 node test feeder. TABLE I: Locations and Capacities of Grid-following and Grid-forming IBDGs in modified IEEE 123 Node Test Feeder. Type | Locations | Capacities
---|---|---
| Grid-following
---
IBDG (1-$\phi$)
| 5, 11, 16, 28, 40, 42,
---
47, 81, 83, 90, 97, 107
110, 116
| 80 kW for single-phase
---
40 kVAr for single-phase
| Grid-following
---
IBDG (3-$\phi$)
| 24, 33, 41, 48, 52,
---
59, 69, 91, 105, 109
| 100 kW per $\phi_{a},\phi_{b},\phi_{c}$
---
50 kVAr per $\phi_{a},\phi_{b},\phi_{c}$
| Grid-forming
---
IBDG (3-$\phi$)
14, 19, 62, 72 |  | 100 kW per $\phi_{a},\phi_{b},\phi_{c}$
---
50 kVAr per $\phi_{a},\phi_{b},\phi_{c}$

###  V-B Sequential Service Restoration Results

As shown in (28), the relationship between the maximum load step and the
frequency nadir is influenced by the value of hyper-parameter $\alpha$ in the
frequency-dynamics constraints. Therefore, different $\alpha$ values may lead
to different service restoration results. In this case, the horizon length $T$
and the hyper-parameter $\alpha$ are set to 4 and 0.1, respectively.

As shown in Fig. 6, the system is partitioned into four MGs by energizing the
switchable lines sequentially, and the radial structure of each MG is
maintained at each stage. Inside each formed MG, the power balance is achieved
between the restored load and power outputs of IBDGs. The value in brackets
nearby each line switch in Fig. 6 represents the number of restoration stage
when it closes. In Table II, the restoration sequences for switchable IBDGs
and loads are shown, where the subscript and superscript are the bus index and
the MG index of grid-following IBDGs and loads, respectively. It can be
observed that MG2 only needs 3 stages to be fully restored, while MG1 and MG3
can restore in 4 stages. However, due to the heavy loading situation, MG4 is
gradually restored in 5 stages to ensure a relatively smooth frequency
dynamics.

Figure 6: Restoration solutions for the formed MG1-MG4, where the restoration stage when line switch closes is shown in red. TABLE II: Restored Grid-following IBDGs and Loads at Each Restoration Stage. | Restoration
---
stage
| Restored
---
grid-following IBDGs
| Restored
---
loads
1 |  | $G_{11}^{1},G_{5}^{1},G_{24}^{2}$
---
$G_{28}^{2},G_{40}^{2},G_{41}^{2}$
$G_{42}^{2},G_{97}^{4}$
| $L_{14}^{1},L_{8}^{1},L_{9}^{1},L_{10}^{1},L_{11}^{1},L_{12}^{1}$
---
$L_{13}^{1},L_{1}^{1},L_{2}^{1},L_{3}^{1},L_{4}^{1},L_{5}^{1},L_{6}^{1}$
$L_{7}^{1},L_{19}^{2},L_{20}^{2},L_{21}^{2},L_{22}^{2},L_{23}^{2}$
$L_{24}^{2},L_{25}^{2},L_{26}^{2},L_{27}^{2}L_{28}^{2},L_{29}^{2}$
$L_{30}^{2},L_{31}^{2},L_{36}^{2},L_{37}^{2}L_{38}^{2},L_{39}^{2}$
$L_{40}^{2},L_{41}^{2},L_{42}^{2},L_{43}^{2},L_{44}^{2},L_{62}^{3}$
$L_{65}^{3},L_{63}^{3},L_{64}^{3},L_{72}^{4},L_{71}^{4},L_{93}^{4}$
$L_{94}^{4},L_{95}^{4},L_{70}^{4},L_{96}^{4},L_{97}^{4},L_{98}^{4}$
$L_{99}^{4},L_{100}^{4},L_{101}^{4},L_{102}^{4},L_{103}^{4}$
$L_{104}^{4}$
2 |  | $G_{33}^{2},G_{47}^{2},G_{48}^{2}$
---
$G_{69}^{3},G_{105}^{4},G_{107}^{4}$
$G_{109}^{4},G_{110}^{4},G_{116}^{4}$
| $L_{32}^{2},L_{33}^{2},L_{34}^{2},L_{35}^{2},L_{45}^{2},L_{46}^{2}$
---
$L_{47}^{2},L_{48}^{2},L_{49}^{2},L_{66}^{3},L_{67}^{3},L_{68}^{3}$
$L_{69}^{3},L_{105}^{4},L_{106}^{4},L_{107}^{4},L_{108}^{4}$
$L_{109}^{4},L_{110}^{4},L_{111}^{4},L_{112}^{4},L_{113}^{4}$
$L_{114}^{4},L_{115}^{4},L_{116}^{4},L_{117}^{4},L_{118}^{4}$
$L_{119}^{4}$
3 |  | $G_{52}^{2}$
---
| $L_{50}^{2},L_{51}^{2},L_{52}^{2},L_{53}^{2},L_{84}^{4},L_{85}^{4}$
---
$L_{86}^{4},L_{87}^{4}$
4 |  | $G_{16}^{1},G_{59}^{3},G_{90}^{4}$
---
$G_{91}^{4}$
| $L_{15}^{1},L_{16}^{1},L_{17}^{1},L_{18}^{1},L_{54}^{3},L_{55}^{3}$
---
$L_{56}^{3},L_{57}^{3},L_{58}^{3},L_{59}^{3},L_{60}^{3},L_{60}^{3}$
$L_{61}^{3},L_{73}^{4},L_{74}^{4},L_{75}^{4},L_{76}^{4},L_{77}^{4}$
$L_{78}^{4},L_{79}^{4},L_{88}^{4},L_{89}^{4},L_{90}^{4},L_{91}^{4}$
$L_{92}^{4}$
5 |  | $G_{81}^{4},G_{83}^{4}$
---
| $L_{80}^{4},L_{81}^{4},L_{82}^{4},L_{83}^{4}$
---

For each restoration stage, the restored loads and frequency nadir in MG1-MG4
are shown in Table III. Total 1773 kW of load are restored at the end of the 5
stages. It can be observed the service restoration actions happened in certain
stages rather than in all stages. For example, MG1 restores 280.5 kW of load
in Stage 1, but it restores no more load until Stage 4. While MG4 takes action
on service restoration in each stage. It is because the sequential service
restoration is limited by operational constraints, among which the maximum
load step in each stage is again limited by the proposed frequency-dynamics
constraints. Note that a larger amount of restored load in the optimization
level will typically cause a lower frequency nadir in the simulation level,
then a low frequency nadir will be considered in constraint (28) and help the
optimization level to restrict a larger amount of restored load in next
restoration stage. Because the first stage is the entry point of the
restoration process, there is no prior frequency nadir information to be used
in constraint (28), therefore, the restored load in the first stage is
typically the largest among all stages, which leads to a corresponding lowest
frequency nadir among all stages.

TABLE III: Restored Loads, Frequency Nadir and Computation Time for MG1-MG4. Cases |  | Restored load
---
(kW)
| Frequency nadir
---
(Hz)
MG1 ($T=4$ and $\alpha=0.1$) | Stage 1 | 280.5 | 59.7044
Stage 2 | 280.5 | 59.9992
Stage 3 | 280.5 | 59.9992
Stage 4 | 346.5 | 59.9200
Stage 5 | 346.5 | 59.9989
MG2 ($T=4$ and $\alpha=0.1$) | Stage 1 | 230.0 | 59.7079
Stage 2 | 360.0 | 59.8201
Stage 3 | 420.0 | 59.9146
Stage 4 | 420.0 | 59.9984
Stage 5 | 420.0 | 59.9984
MG3 ($T=4$ and $\alpha=0.1$) | Stage 1 | 212.5 | 59.7116
Stage 2 | 212.5 | 59.9990
Stage 3 | 212.5 | 59.9990
Stage 4 | 382.5 | 59.7656
Stage 5 | 382.5 | 59.9985
MG4 ($T=4$ and $\alpha=0.1$) | Stage 1 | 192.0 | 59.7910
Stage 2 | 324.0 | 59.8541
Stage 3 | 414.0 | 59.9003
Stage 4 | 570.0 | 59.8230
Stage 5 | 624.0 | 59.9364

The comparison of total restored loads with and without considering the
proposed frequency dynamics constraints is shown in Fig. 7. Note that the
total amount of restorable load of the base case model (i.e., without the
frequency dynamics constraints) is the same as that of the proposed model with
the frequency dynamics constraints. That is because the total load of the test
system is fixed and less than the total DG generation capacity in both models.
However, the base case needs 6 stages to fully restore the all the loads,
while the proposed model can achieve that goal in the first 5 stages (as it is
observed, no more loads between Stage 5 and Stage 6 are restored). While In
the early stages 1 to 3, the restored load of the proposed model is a little
bit less than the base case. A further analysis is that: during the early
restoration stages, the proposed model generated a restoration solution that
prevents too low frequency nadir during transients. The base case restores
more loads at Stage 1 to Stage 3 without considering such limitation on the
frequency nadir. However, Stage 4 is a turning point when the proposed model
restores more loads than the base case. Therefore, the proposed model restores
less loads than the base case during early stages (here, Stage 1 to Stage 3),
while it restores more loads than the base case during later stages (from
Stage 4). Such restoration pattern (restored load at each stage) of the base
case model and the proposed model may vary case by case if the system topology
or other operational constraints are changed. Therefore, if we implement the
base case model and the proposed model in another test system with different
topology or constraint settings, the base case model may restore fewer loads
than the proposed model in the early stages and the turning point stage may
change as well.

Figure 7: Total restored load with and without considering frequency dynamics
constraints.

In Fig. 8a and Fig. 8b, a zoom in view of the frequency response of MG4 and
the frequency response of MG4 in Stage 1 are shown for better observation of
the frequency dynamic performance. The frequency responses with and without
the frequency dynamics constraints are represented by blue and red lines,
respectively. By this comparison, it can be observed that both the rate of
change of frequency and frequency nadir are significantly improved by
considering frequency dynamics constraints in the proposed restoration model.
However, if the frequency dynamics constraints are not considered to prevent a
large frequency drop, unstable frequency oscillation may happen. The reason of
the oscillation phenomenon in Fig. 8b is the too large $P^{L}$, which deviates
the initial state of MG in the current stage out of the region of attraction
of the original stable equilibrium. This in turn demonstrates the necessity to
incorporate that frequency dynamics constraint in the optimization level. Note
that $\omega_{\rm set}$ is set to 60 Hz in the droop equation (34), the
equilibrium can be achieved when $\omega=\omega_{\rm set}$ and $P=P^{L}$,
which means that the output frequency tracks the frequency reference when the
output power of the simulation level tracks the target restored load
calculated from the optimization level.

Figure 8: Frequency responses of MG4 with and without frequency dynamics
constraints: (a) Subplot of frequency response of MG4 during 5.0 s to 5.8 s;
(b) Frequency responses of MG4 in Stage 1.

Fig. 9 shows the frequency responses of each inverter-dominated MG based on
the proposed restoration model. The results show that the MG frequency drops
when the load is restored. Because the maximum load step is constrained in the
proposed MILP-based sequential service restoration model, the frequency nadir
is also constrained. When load is restored as the frequency drops, the
frequency nadir can be effectively maintained above the $f^{\rm min}$
threshold.

Figure 9: Frequency responses of inverter-dominated MGs: (a) MG1; (b) MG2; (c)
MG3; (d) MG4.

###  V-C Impact of Hyper-parameters in Frequency Dynamics Constraints

Compared to other MGs, MG4 is heavily loaded with the largest number of nodes.
Based on the results of Fig. 6, MG4 needs more stages to be fully restored
compared to other MGs. Therefore, MG4 is chosen to test the effect of
different $\alpha$ values. In Fig. 10a and Fig. 10b, the frequency responses
of MG4 during the period of 3.1 s to 5.1 s, the period of 9.3 s to 11.3 s and
the whole restoration process are shown, where the frequency with
$\alpha=0.1$, $\alpha=0.2$ and $\alpha=1.0$ are represented by blue solid
line, red dashed line and yellow dotted line, respectively. It can be observed
that 5 stages are required to fully restore all the loads when $\alpha=0.1$;
while only 4 restoration stages are needed when $\alpha=0.2$ or $\alpha=1.0$.
During the period of 3.1 s to 5.1 s in left of Fig. 10a, the frequency nadirs
with $\alpha=0.2$ or $\alpha=1.0$ are lower than the frequency nadir with
$\alpha=0.1$, which means more loads can be restored with larger value of
$\alpha$. During the period of 9.3 s to 11.3 s in right of Fig. 10b, the
frequency nadir with $\alpha=0.1$ is lower than the frequency nadirs with
$\alpha=0.2$ and $\alpha=1.0$, it is because the total restored loads for
different $\alpha$ values are same, with $\alpha=0.2$ or $\alpha=1.0$, it can
restore more loads in the early restoration stage, therefore they just need
less loads to be restored in the late restoration stage. However, $\alpha=0.1$
restores less loads in the early restoration stage, it has to restore more
loads in the late restoration stage. As shown in Fig. 10c, the overall dynamic
frequency performance with $\alpha=0.1$ is still better than the cases with
$\alpha=0.2$ and $\alpha=1.0$. Hence, there is a trade-off between dynamic
frequency performance and restoration performance regarding the choice of
$\alpha$: too small $\alpha$ may lead to too slow restoration and the
frequency nadir may be high in the early restoration stage and the frequency
nadir may be low in the late restoration stage; in turn, a large $\alpha$ may
lead to less number of restoration stages, too large $\alpha$ may cause too
low frequency in early stages and deteriorate the dynamic performance of the
system frequency in a practical restoration process.

Figure 10: Frequency responses of MG4 with different $\alpha$: (a) Frequency
responses during 3.1 s to 5.1 s; (b) Frequency responses during 9.3 s to 11.3
s; (c) Frequency responses during the whole restoration process.

We also shows that different values of the horizon length $T$ may cause
different service restoration results. Table IV summarizes the total restored
loads and computation time using different horizon lengths in the proposed
service restoration model. On the one side, the restored loads of case with
$T=2$ and $T=3$ are less than that of the cases with $T\geq 4$, where the
total restored load can reach the maximum level. Therefore, the results with
small number of horizon length $T=2$ and $T=3$ are sub-optimal restoration
solutions. On the other side, the longer horizon length also leads to heavy
computation burden and increase the computation time. Similar to the impact of
$\alpha$, there can be a trade-off between the computation time and the
quality of solution when determining the value of $T$.

TABLE IV: Restored Loads, Frequency Nadir and Computation Time with Different Horizon Lengths. | Total restored load (kW) | Computation time (s)
---|---|---
$T=2$ | 1362.5 | 26.8870
$T=3$ | 1410.5 | 32.6725
$T=4$ | 1773.0 | 48.5629
$T=5$ | 1773.0 | 61.9968
$T=6$ | 1773.0 | 88.0216

In Fig. 11, the frequency responses of MG1 to MG4 are depicted during the
restoration process with different values of droop gain $D_{p}$. In the test
case, the original setting of $D_{p}$ is $1\times 10^{-5}$. It can be observed
that the different values of $D_{p}$ will cause different restoration
solutions and frequency responses. As indicated by the arrow in Fig. 11a, MG1
can be fully restored in four stages when $D_{p}=1\times 10^{-5}$ or $2\times
10^{-5}$, however, if the $D_{p}=3\times 10^{-5}$, MG1 needs five stages to be
fully restored. Similar observation can be found for restoration stage in Fig.
11c for MG3, it needs five stages to be fully restored when $D_{p}$ equals
larger values (such as $2\times 10^{-5}$ or $3\times 10^{-5}$), while it only
needs four stages when $D_{p}$ equals smaller values (such as $=1\times
10^{-5}$). As shown in Fig. 11b and Fig. 11d, larger value of $D_{p}$ will
also lead to larger frequency drop during restoration process.

Figure 11: Frequency responses of inverter-dominated MGs with different values
of $D_{p}$ during restoration process: (a) MG1; (b) MG2; (c) MG3; (d) MG4.

##  VI Conclusion

To improve the dynamic performance of the system frequency during service
restoration of a unbalanced distribution systems in an inverter-dominated
environment, we propose a simulation-assisted optimization model considering
frequency dynamics constraints with clear physical meanings. Results
demonstrate that: (i) The proposed frequency dynamics constrained service
restoration model can significantly reduce the transient frequency drop during
MGs forming and service restoration. (ii) Other steady-state performance
indicators of our proposed method can rival that of the conventional methods,
in terms of the final restored total load and the required number of
restoration stages. Investigating on how to choose the best hyper-parameters,
such as $\alpha$, horizon length $T$ and droop gain $D_{p}$ will be the next
research direction.

## References

  * [1] E. O. of the President, “Economic benefits of increasing electric grid resilience to weather outages,” White House, Tech. Rep., 2020.
  * [2] A. M. Salman, Y. Li, and M. G. Stewart, “Evaluating system reliability and targeted hardening strategies of power distribution systems subjected to hurricanes,” _Reliab. Eng. Syst. Saf._ , vol. 144, pp. 319–333, Dec. 2015\.
  * [3] H. Haggi, R. R. nejad, M. Song, and W. Sun, “A review of smart grid restoration to enhance cyber-physical system resilience,” in _2019 IEEE Innovative Smart Grid Technologies - Asia (ISGT Asia)_ , 2019, pp. 4008–4013.
  * [4] Z. Wang and J. Wang, “Self-healing resilient distribution systems based on sectionalization into microgrids,” _IEEE Trans. Power Syst._ , vol. 30, pp. 3139–3149, Nov. 2015.
  * [5] A. Arif and Z. Wang, “Networked microgrids for service restoration in resilient distribution systems,” _IET Gener. Transm. Distrib._ , vol. 11, no. 14, pp. 3612–3619, Aug. 2017.
  * [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] S. Yao, P. Wang, and T. Zhao, “Transportable energy storage for more resilient distribution systems with multiple microgrids,” _IEEE Trans. on Smart Grid_ , vol. 10, pp. 3331–3341, May 2019.
  * [8] L. Che and M. Shahidehpour, “Adaptive formation of microgrids with mobile emergency resources for critical service restoration in extreme conditions,” _IEEE Trans. Power Syst._ , vol. 34, no. 1, pp. 742–753, Jan. 2019.
  * [9] B. Chen, C. Chen, J. Wang, and K. L. Butler-Purry, “Sequential service restoration for unbalanced distribution systems and microgrids,” _IEEE Trans. Power Syst._ , vol. 33, pp. 1507–1520, Mar. 2018.
  * [10] Y. Wen, W. Li, G. Huang, and X. Liu, “Frequency dynamics constrained unit commitment with battery energy storage,” _IEEE Trans. Power Syst._ , vol. 31, no. 6, pp. 5115–5125, Nov. 2016.
  * [11] H. Gu, R. Yan, T. K. Saha, E. Muljadi, J. Tan, and Y. Zhang, “Zonal inertia constrained generator dispatch considering load frequency relief,” _IEEE Trans. Power Syst._ , vol. 35, no. 4, pp. 3065–3077, Jul. 2020.
  * [12] Y. Wen, C. Y. Chung, X. Liu, and L. Che, “Microgrid dispatch with frequency-aware islanding constraints,” _IEEE Trans. Power Syst._ , vol. 34, no. 3, pp. 2465–2468, May 2019.
  * [13] O. Bassey, K. L. Butler-Purry, and B. Chen, “Dynamic modeling of sequential service restoration in islanded single master microgrids,” _IEEE Trans. Power Syst._ , vol. 35, no. 1, pp. 202–214, Jan. 2020.
  * [14] B. Qin, H. Gao, J. Ma, W. Li, and A. Y. Zomaya, “An input-to-state stability-based load restoration approach for isolated power systems,” _Energies_ , vol. 11, pp. 597–614, Mar. 2018.
  * [15] Y. Xu, C. Liu, K. P. Schneider, F. K. Tuffner, and D. T. Ton, “Microgrids for service restoration to critical load in a resilient distribution system,” _IEEE Trans. Smart Grid_ , vol. 9, pp. 426–437, Jan. 2018.
  * [16] Y. Du, X. Lu, J. Wang, and S. Lukic, “Distributed secondary control strategy for microgrid operation with dynamic boundaries,” _IEEE Trans. Smart Grid_ , vol. 10, no. 5, pp. 5269–5285, Sept. 2019.
  * [17] B. K. Poolla, D. Grob, and F. Dorfler, “Placement and implementation of grid-forming and grid-following virtual inertia and fast frequency response,” _IEEE Trans. Power Syst._ , vol. 34, pp. 3035–3046, Jul. 2019\.
  * [18] P. Vorobev, P. Huang, M. A. Hosani, J. L. Kirtley, and K. Turitsyn, “High-fidelity model order reduction for microgrids stability assessment,” _IEEE Trans. Power Syst._ , vol. 33, pp. 874–887, Jan. 2018.
  * [19] J. M. Guerrero, L. Hang, and J. Uceda, “Control of distributed uninterruptible power supply systems,” _IEEE Trans. Ind. Electron._ , vol. 55, no. 8, pp. 2845–2859, 2008.
  * [20] K. Y. Yap, C. R. Sarimuthu, and J. M.-Y. Lim, “Virtual inertia-based inverters for mitigating frequency instability in grid-connected renewable energy system: A review,” _Appl. Sci._ , vol. 9, no. 24, p. 5300, Dec. 2019.
  * [21] H. Bevrani, T. Ise, and Y. Miura, “Virtual synchronous generators: A survey and new perspectives,” _Int. J. Electr. Power Energy Syst._ , vol. 54, pp. 244–254, Jan. 2014.
  * [22] Y. Zhu, C. Liu, K. Sun, D. Shi, and Z. Wang, “Optimization of battery energy storage to improve power system oscillation damping,” _IEEE Trans. Sustain. Energy_ , vol. 10, no. 3, pp. 1015–1024, 2019.
  * [23] Y. Kim, J. Wang, and X. Lu, “A framework for load service restoration using dynamic change in boundaries of advanced microgrids with synchronous-machine dgs,” _IEEE Trans. Smart Grid_ , vol. 9, no. 4, pp. 3676–3690, Jul. 2018.
  * [24] Z. Wang, J. Wang, B. Chen, M. M. Begovic, and Y. He, “Mpc-based voltage/var optimization for distribution circuits with distributed generators and exponential load models,” _IEEE Trans. Smart Grid_ , vol. 5, no. 5, pp. 2412–2420, 2014.
  * [25] B. A. Robbins and A. D. Domínguez-García, “Optimal reactive power dispatch for voltage regulation in unbalanced distribution systems,” _IEEE Trans. Power Syst._ , vol. 31, no. 4, pp. 2903–2913, 2016.
  * [26] Q. Zhang, K. Dehghanpour, and Z. Wang, “Distributed CVR in unbalanced distribution systems with PV penetration,” _IEEE Trans. Smart Grid_ , vol. 10, no. 5, pp. 5308–5319, Sept. 2019.
  * [27] Y. Yuan, K. Dehghanpour, F. Bu, and Z. Wang, “Outage detection in partially observable distribution systems using smart meters and generative adversarial networks,” _IEEE Trans. Smart Grid_ , vol. 11, no. 6, pp. 5418–5430, Nov. 2020.
  * [28] 123-bus feeder. [Online]. Available: https://site.ieee.org/pes-testfeeders/resources/

|  Qianzhi Zhang (S’17) is currently pursuing his Ph.D. in the Department of
Electrical and Computer Engineering, Iowa State University, Ames, IA. He
received his M.S. in electrical and computer engineering from Arizona State
University in 2015. He has worked with Huadian Electric Power Research
Institute from 2015 to 2016 as a research engineer. His research interests
include the applications of machine learning and advanced optimization
techniques in power system operation and control.
---|---
|  Zixiao Ma (S’18) is currently a Ph.D. student in the Department of
Electrical and Computer Engineering at the Iowa State University, Ames, IA,
USA. He received his B.S. degree in Automation and M.S. degree in Control
theory and Control Engineering from Northeastern University in 2014 and 2017
respectively. His research interests are focused on the power system load
modeling, microgrids, nonlinear control and model reduction.
---|---
|  Yongli Zhu (S’12) received his B.S. degree from Huazhong University of
Science and Technology in 2009, M.S. degree from State Grid Electric Power
Research Institute in 2012, and Ph.D. degree from the University of Tennessee,
Knoxville in 2018. He joined Iowa State University in the position of postdoc
researcher in 2020. His research interests include power system stability,
microgrid, and machine learning applications in power systems.
---|---
|  Zhaoyu Wang (S’13–M’15–SM’20) is the Harpole-Pentair Assistant Professor
with Iowa State University. He received the B.S. and M.S. degrees in
electrical engineering from Shanghai Jiaotong University, and the M.S. and
Ph.D. degrees in electrical and computer engineering from Georgia Institute of
Technology. His research interests include optimization and data analytics in
power distribution systems and microgrids. He is the Principal Investigator
for a multitude of projects focused on these topics and funded by the National
Science Foundation, the Department of Energy, National Laboratories, PSERC,
and Iowa Economic Development Authority. Dr. Wang is the Chair of IEEE Power
and Energy Society (PES) PSOPE Award Subcommittee, Co-Vice Chair of PES
Distribution System Operation and Planning Subcommittee, and Vice Chair of PES
Task Force on Advances in Natural Disaster Mitigation Methods. He is an editor
of IEEE Transactions on Power Systems, IEEE Transactions on Smart Grid, IEEE
Open Access Journal of Power and Energy, IEEE Power Engineering Letters, and
IET Smart Grid. Dr. Wang was the recipient of the National Science Foundation
(NSF) CAREER Award, the IEEE PES Outstanding Young Engineer Award, and the
Harpole-Pentair Young Faculty Award Endowment.
---|---