# Frequency-Constrained Resilient Scheduling of Microgrid: A Distributionally
Robust Approach

Zhongda Chu, Ning Zhang, and Fei Teng  Zhongda Chu and Fei Teng (Corresponding
author) are with Department of Electrical and Electronic Engineering, Imperial
College London, U.K. Ning Zhang is with Department of Electrical Engineering,
Tsinghua University, China.

###### Abstract

In order to prevent the potential frequency instability due to the high Power
Electronics (PE) penetration under an unintentional islanding event, this
paper presents a novel microgrid scheduling model which explicitly models the
system frequency dynamics as well as the long/short term uncertainty
associated with renewable energy resources and load. Synthetic Inertia (SI)
control is applied to regulating the active power output of the Inverter-Based
Generators (IBGs) to support the post-islanding frequency evaluation. The
uncertainty associated with the noncritical load shedding is explicitly
modeled based on the distributionally robust formulation to ensure resilient
operation during islanding events. The resulted frequency constraints are
derived analytically and reformulated into Second-Order Cone (SOC) form, which
are further incorporated into the microgrid scheduling model, enabling optimal
SI provision of Renewable Energy Sources (RESs) from the micorgrid
perspective. With the SOC relaxation of the AC power flow constraints, the
overall problem is constructed as a mixed-integer SOC Programming (MISOCP).
The effectiveness of the proposed model is demonstrated based on modified IEEE
14-bus system.

###### Index Terms:

microgrid scheduling, frequency dynamics, synthetic inertia, distributionally
robust optimization

##  I Introduction

Microgrids, distribution systems integrated with large scale of RESs, storage
devices and controllable loads have been a promising concept for reliable and
flexible electricity supply in an environment-friendly manner [1]. They are
connected to the main grid at the Point of Common Coupling (PCC), providing
the capability of power transmission in both directions. Microgrids can
operate in islanded mode by disconnecting itself from the main grid when
subjected to external disturbances, making them highly beneficial to customers
and utilities [2].

Due to the resiliency benefits of microgrids, extensive research has been
conducted on microgrid scheduling aiming to achieve optimal microgrid
operation and management [3, 4]. A stochastic microgrid scheduling model is
proposed in [5] to address the intermittency and variability of the RESs.
Applying the chance constrained approach, the authors in [6] formulate the
grid-connected microgrid scheduling problem as linear programming. The studies
in [7] develop a distributionally robust chance-constrained energy management
for islanded microgrids with the uncertainty of renewable generation captured
through a novel ambiguity set. [8] presents a resiliency-oriented microgrid
optimal scheduling model in islanded operating mode to minimize the microgrid
load curtailment. The demand response from Electric Vehicles (EV) is
incorporated into the islanded microgrid scheduling problem to minimize the
operating and EV charging costs [9].

Most of the literature in this vein focuses on the microgrid scheduling in
either grid-connected or islanded mode with less attention being paid on the
influence of the transition between the two modes. The islanding events of
microgrids force its exchanged power with the main grid to zero, resulting in
power unbalance between generation and demand. Furthermore, due to the high PE
penetration in microgrids, this unbalanced power can lead to large frequency
deviations, even blackout and system collapse.

The authors in [10] consider multi-period islanding constraints in a
centralized microgrid optimal scheduling model. The solution is examined for
islanding to ensure the microgrid has sufficient online capacity for quickly
switching to the islanded mode on request. [11, 12] propose an optimal
scheduling strategy for microgrid operation considering constraints of
islanding capability. Probability of Successful Islanding (PSI) is introduced
as a new concept to ensure there is enough reserve to cover the load after
islanding events. Similarly, in [13], a new optimal strategy for scheduling of
reconfigurable microgrids is presented while maintaining the PSI above a
certain level. Considering the uncertainty of RESs and demand, the microgrid
scheduling is formulated as a chance-constrained global optimization problem.

However, the above methods only focus on determining the feasible region of
the spinning reserve to guarantee the PSI while neglecting the detailed
frequency dynamics and the frequency support from IBGs. It is improved by [14]
where a comprehensive optimization and real-time control framework for
maintaining frequency stability of multi-microgrid networks under islanding
events are proposed. The frequency dynamics and SI from IBGs are also
considered leading to highly nonlinear frequency constraints. An iterative
algorithm is developed based on the cutting plan approach to incorporate the
post-contingency frequency constraints, which increases its implementing
complexity. Deep learning is applied in [15] to approximate the nonlinear
nadir constraint using a neural network such that an MILP-based microgrid
scheduling problem can be formulated. Nevertheless, the detailed SI modeling
from IBGs is not considered and the uncertainty associated with the
noncritical load shedding due to the forecasting errors and the relays has not
been discussed.

Motivated by previous research, this paper proposes a novel optimal microgrid
scheduling model considering the state-of-art SI control scheme from PV-
storage system and Wind Turbines (WTs) to minimize the microgrid operation
cost while ensuring the frequency security after islanding events. The
contributions of this paper are summarized as follow:

  * •

A novel microgrid scheduling model is proposed, which optimizes microgrid
operating conditions, noncritical load shedding as well as the SI from IBGs
such that the frequency constraints after microgrid islanding events can be
maintained.

  * •

A distributionally robust approach is adopted to account for the uncertainty
associated with noncritical load shedding leading to distributionally robust
chance constraint formulation of the frequency metrics.

  * •

The highly nonlinear frequency constraints are further effectively
reformulated into SOC form and incorporated into the microgrid scheduling
model, resulting in an overall MISOCP together with the SOC approximation of
the AC power flow.

The rest of this paper is structured as follows. Section II introduces the SI
control modeling of RESs and the overall microgrid frequency dynamics, based
on which the analytical expressions of the frequency metrics are derived. The
uncertainty associated with the noncritical load shedding is discussed in
Section III, leading to distributionally robust frequency constraints and the
nonlinear nadir constraint is further reformulated into SOC form. The overall
MISOCP-based microgrid scheduling model is presented in Section IV, followed
by case studies (Section V) illustrating the value and performance of the
proposed model. Finally, section VI concludes the paper.

##  II Frequency Dynamics of Microgrid with SI Provision from RESs

In this section, the microgrid frequency dynamics are investigated. The Rate
of Change of Frequency (RoCoF), frequency nadir and steady-state value
subsequent to islanding events are derived considering the provision of SI
from RESs. Due to the low efficiency during normal operation, the deloading
control strategy is not considered in this paper. Instead, in order to
maximize the energy captured by PV systems and WTs, we assume that they are
controlled through MPPT strategy during normal operation. For PV systems,
energy storage devices are used to provide frequency support during the system
disturbance, whereas stored kinetic energy is utilized in the WTs. It is the
synthetic inertia rather than the conventional droop control that is applied
in the proposed model to provide higher power injection at the beginning of
islanding events. It should be noted that only loss of generation (or increase
of load) is considered in this paper since the opposite situation can be
easily dealt with by shifting the operating point of the RESs away from the
optimal power point.

###  II-A Modeling of Frequency Support from Energy Storage Devices

A state-of-the-art VSC control scheme previously described in [16] is adapted
to provide constant synthetic inertia during a disturbance. Furthermore, since
the power injection associated with synthetic inertia approaches zero at the
frequency nadir and has no impact on the steady-state frequency, constant
power can be injected into the grid after the frequency nadir to improve the
steady-state frequency if needed. The active power setpoint is adjusted
through the outer control loop to achieve the desired power injection to the
grid. Note that although it is possible to achieve more sophisticated control
laws, e.g., adaptive virtual synchronous machine or online model predictive
control, it would make the scheduling model much more complicated, thus not
being considered in this paper.

In order to determine an optimal and feasible frequency support provision
during the scheduling process, it is essential to consider the limitation of
both instantaneous power and total energy of the energy storage devices.

The instantaneous output power of the energy storage device $b\in\mathcal{B}$
during an entire frequency event, $\mathcal{T}_{0}$ is confined by the maximum
charging/discharging rate,
$\bar{P}_{b}^{\mathrm{ch}}$/$\bar{P}_{b}^{\mathrm{dch}}$:

$\bar{P}_{b}^{\mathrm{ch}}\leq
P_{b}-2H_{s_{b}}\Delta\dot{f}(t)\leq\bar{P}_{b}^{\mathrm{dch}},\,\,\;\forall
t\in\mathcal{T}_{0}$ (1)

where $P_{b}$ and $H_{s_{b}}$ are the normal output power and synthetic
inertia from the storage unit $b$. Equation (1) is equivalent to:

$\bar{P}_{b}^{\mathrm{ch}}\leq
P_{b}+\max_{t\in\mathcal{T}_{0}}\Big{\\{}2H_{s_{b}}\left|\Delta\dot{f}(t)\right|\Big{\\}}\leq\bar{P}_{b}^{\mathrm{dch}}.$
(2)

Since $\left|\Delta\dot{f}(t)\right|$ is constrained by the RoCoF limit:

$0\leq\left|\Delta\dot{f}(t)\right|\leq\Delta\dot{f}(0)\leq\Delta\dot{f}_{\mathrm{lim}},$
(3)

(2) can be rewritten in linear form as follows:

$\bar{P}_{b}^{\mathrm{ch}}\leq
P_{b}+2H_{s_{b}}\Delta\dot{f}_{\mathrm{lim}}\leq\bar{P}_{b}^{\mathrm{dch}}.$
(4)

The energy required by synthetic inertia provision is negligible due to the
small time scale of the inertial response. Hence, it is not considered as a
constraint here.

Similarly, the limitation for the constant power provision $\Delta P_{C}$
after the frequency nadir can be formulated as:

$\displaystyle\bar{P}_{b}^{\mathrm{ch}}\leq P_{b}+\Delta P_{C}$
$\displaystyle\leq\bar{P}_{b}^{\mathrm{dch}}$ (5a) $\displaystyle\Delta
P_{C}T_{s}$ $\displaystyle\leq\mathrm{SoC}_{b}\cdot E_{c,b}$ (5b)

where $T_{s}$ is the time interval of the constant power provision; $E_{c}$
and $\mathrm{SoC}$ are the energy capacity and the state of charge of the
storage device. Note that the above frequency support model and the
operational constraints can also be applied to other microgrid battery storage
units that are not connected with PV systems. In addition to the synthetic
inertia provision, the energy storage devices are responsible to smooth the
renewable generation fluctuation in shorter timescale and balance the
inconsistent profiles between the load and generation in longer timescale
(storage the energy when the generation is more than the demand and release
the energy otherwise). All these services are simultaneously optimized during
the microgrid scheduling process in our proposed model, such that the storage
capacity can be optimally allocated in real time.

###  II-B Modeling of Frequency Support from WTs

The control framework proposed in [17] is applied to provide optimal synthetic
inertia from WTs. In the proposed model, active power is extracted from the
stored kinetic energy of WTs to facilitate the frequency evaluation during the
disturbance. Due to the complexity caused by incorporating the WT mechanical
dynamics into the control design and the restriction of the stored kinetic
energy, only short-term inertial response is provided from WTs.

Figure 1: Block diagram of WT synthetic inertia control scheme.

Furthermore, the secondary frequency dip associated with the rotor speed
deviation from the optimal operation point before the disturbance can be
eliminated. This is achieved by adding a Mechanical Power Estimator (MPE) in
the active power control loop of the WT grid-side converter as shown in Fig.
1. As a result, the additional output power from WTs ($\Delta P_{w}$) during a
system disturbance is the sum of synthetic inertia power $P_{SI}$ and the
output of the MPE $\Delta\tilde{P}_{a}$:

$\Delta P_{w}=-2H_{s}\Delta\dot{f}+\Delta\tilde{P}_{a}.$ (6)

In order to incorporate (6) into the system frequency dynamics, the highly
nonlinear expression of $\Delta\tilde{P}_{a}$ is further approximated by a
negative system damping term:

$\Delta\tilde{P}_{a}=D_{s}\Delta f=\gamma H_{s}^{2}\Delta f.$ (7)

In addition, the total available synthetic inertia form WTs in the system can
be estimated given the wind speed distribution as proposed in [17], where the
feasibility of frequency support from WTs and detailed control performance can
be found as well.

###  II-C Frequency Evaluation under Islanding Events

Without the frequency support from the RESs, the frequency dynamics in a
multi-machine microgrid can be expressed in the form of a single swing
equation, under the premise of the Centre-of-Inertia (CoI) model [18]:

$2H_{c}\frac{\partial\Delta f(t)}{\partial t}=-D_{0}\Delta f(t)+\Delta
R(t)-\underbrace{(\Delta P_{L_{0}}-\Delta P_{D})}_{\Delta P_{L}},$ (8)

where $\Delta P_{L_{0}}$, the loss of generation due to the islanding event at
$t=0$ is a decision variable and can be viewed as a step disturbance. $\Delta
P_{D}$ is the noncritical load shedding in order to maintain the post-
contingency frequency within the limits, which is a common practice in
microgrid after islanding events. It can be deferred or curtailed in response
to economic incentives or islanding requirements. Furthermore, $\Delta P_{D}$
is modeled as a decision variable with uncertainty and its mean ($\Delta
P_{D_{\mu}}$) and standard deviation ($\sigma$) are assumed to be partially
known. More details regarding to the uncertainty of $\Delta P_{D}$ are
discussed in Section III. $\Delta P_{L}$ is the equivalent loss of generation
which is always positive, as there is no point to shed more load than the lost
power. Moreover, the PFR $\Delta R(t)$ from conventional Synchronous
Generators (SGs) can be represented according to the following scheme [19]:

$\Delta R(t)=\begin{cases}\frac{R}{T_{d}}t&,\;0\leq t<T_{d}\\\ R&,\;T_{d}\leq
t\end{cases}$ (9)

with $T_{d}$ being the PFR delivered time and R being the total PFR delivered
by time $T_{d}$; The total inertia of SGs is computed as:

$H_{c}=\frac{\sum_{g\in\mathcal{G}}H_{g}P_{g}^{\mathrm{max}}y_{g}}{f_{0}}.$
(10)

Incorporate the frequency support from RESs as described in Section II-A and
II-B into (8) leading to:

$\displaystyle
2{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\underbrace{\left(H_{c}+\sum_{b\in\mathcal{B}}H_{s_{b}}+\sum_{w\in\mathcal{W}}H_{s_{w}}\right)}_{H}}\frac{\partial\Delta
f(t)}{\partial t}$ (11)
$\displaystyle=-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\underbrace{\left(D_{0}-\sum_{w\in\mathcal{W}}\gamma_{w}H_{s_{w}}^{2}\right)}_{D}}\Delta
f(t)+\Delta R(t)-\Delta P_{L}$

where $H$ and $D$ are the overall system inertia and damping respectively;
$b\in\mathcal{B}$ and $w\in\mathcal{W}$ are the set of energy storage units
and wind generation units. Note that the inertia in the system is now the
combination of SGs’ ($H_{c}$) and the SI from the energy storage devices
($H_{s_{b}}$) and WTs ($H_{s_{w}}$). The system damping is decreased by
$\sum_{w\in\mathcal{W}}\gamma_{w}H_{s_{w}}^{2}$ due to the SI provision from
WTs. It represents the the side effect of SI provision from wind turbines
through overproduction scheme, i.e., the output power reduction due to the
deviation from the optimal operating point [17].

Based on the frequency dynamics, the analytical expression of the maximum
instantaneous RoCoF
$(\Delta\dot{f}_{\mathrm{max}}\equiv\Delta\dot{f}|_{t=0^{+}})$ is identified
as:

$\Delta\dot{f}|_{t=0^{+}}=-\frac{\Delta P_{L}}{2H}.$ (12)

It can be maintained within the RoCoF limits by choosing an appropriate system
inertia $H$. Solving the differential equation (11) gives the microgrid
frequency evaluation during an islanding event:

$\Delta f(t)=\left(\frac{\Delta
P_{L}}{D}+\frac{2HR}{T_{d}D^{2}}\right)\left(e^{-\frac{D}{2H}t}-1\right)+\frac{R}{T_{d}D}t,$
(13)

valid $\forall t\in[0,t_{n}]$. The time instant $t_{n}$ of frequency nadir is
derived by setting the derivative of (13) to zero:

$\Delta\dot{f}(t_{n})=0\longmapsto
t_{n}=\frac{2H}{D}\ln{\left(\frac{T_{d}D\Delta P_{L}}{2HR}+1\right)}.$ (14)

Substituting (14) into (13) leads to the expression for frequency nadir
$(\Delta f_{\mathrm{max}}\equiv\Delta f(t_{n}))$:

$\Delta f(t_{n})=\frac{2HR}{T_{d}D^{2}}\ln{\left(\frac{T_{d}D\Delta
P_{L}}{2HR}+1\right)}-\frac{\Delta P_{L}}{D}.$ (15)

It is understandable that the dependence of the frequency nadir on $H$, $D$,
and $\Delta P_{L}$ through a highly-nonlinear relationship makes it difficult
to be incorporated into the microgrid scheduling model. A SOC reformulation is
proposed to cope with this problem as demonstrated in Section III-A.

It should be noted that in order to derive the analytical expressions of
maximum RoCoF and frequency nadir, only the inertial response from RESs are
incorporated in (11). However, when deriving the steady-state frequency
$(\Delta f_{\mathrm{max}}^{\mathrm{ss}}\equiv\Delta f|_{t=\infty})$, the
constant power injection $\Delta P_{C}$ from energy storage devices needs to
be considered:

$\Delta f|_{t=\infty}=\frac{R+\Delta P_{C}-\Delta P_{L}}{D}.$ (16)

Note that the secondary frequency response from conventional generators and
the associated frequency restoration process after steady-state are not
considered in this paper.

Having obtained the analytical expressions of the frequency metrics during an
islanding event (12) (15) and (16), they should be bounded within predescribed
limits in the microgrid scheduling model by selecting proper $H$, $R$, $\Delta
P_{C}$ and $\Delta P_{L}$. These quantities can all be decided
deterministically except the equivalent loss of generation $\Delta P_{L}$ due
to the uncertainty associated with the noncritical load shedding $\Delta
P_{D}$. At the beginning of an islanding event, the noncritical load is
disconnected from the microgrid in order to support the frequency evaluation.
However, the exact value of the shed load is unknown during the scheduling
period, thus increasing the complexity of the scheduling problem as elaborated
in the next section.

##  III Distributionally Robust Chance Constraints of Frequency Metrics

The uncertainty associated with the noncritical load shedding $\Delta P_{D}$
stems from different aspects. On the one hand, forecasting error always exists
during the microgrid scheduling process in terms of the actual demand and the
noncritial load percentage. On the other hand, depending on the specific load
shedding strategies[20, 21], the uncertainty level of $\Delta P_{D}$ varies.
Traditionally, the practice setting of the RoCoF and frequency relays are
mainly based on the experts’ experiment [22]. More advanced noncritical load
control schemes have also been proposed to mitigate frequency variations
during islanding events. For instance, the emergence demand response has been
considered in the setting of the under frequency load shedding relays and the
design of the load shedding schemes [23, 24, 25, 26]. Therefore, it is
complicated to derive the detailed distribution of $\Delta P_{D}$ through
either model-based or data-driven approaches at microgrid scheduling stage. As
a result, the equivalent loss of generation $\Delta P_{L}$ as defined in (8)
also presents uncertainty of the same level. In order to account for the
uncertainty of $\Delta P_{L}$, the frequency constraints are reformulated
through distributionally robust optimization. Assume the first- and second-
order moments of $\Delta P_{L}$ are decision-dependent whereas the exact
probability distribution $\mathbf{D}$ is unknown. This is modeled by the
following ambiguity set:

$\displaystyle\mathcal{P}=\Big{\\{}\mathbf{D}\in\Phi(\Delta P_{L}):\;$
$\displaystyle\mathbb{E}^{\mathbf{D}}(\Delta P_{L})=\Delta P_{L_{\mu}},$
$\displaystyle\mathrm{Var}^{\mathbf{D}}(\Delta P_{L})=\sigma^{2}\Big{\\}}$
(17)

where $\Phi(\cdot)$ is the probability density function; $\Delta
P_{L_{\mu}}=\Delta P_{L_{0}}-\Delta P_{D_{\mu}}$ and $\sigma$ denote the mean
and standard deviation of the equivalent loss of generation given the
distribution $\mathbf{D}$. Since the more noncritical load is about to be
shed, the higher uncertainty level it presents, it is reasonable to assumed
that $\sigma$ depends on the decision variable $\Delta P_{D_{\mu}}$ with a
linear coefficient $\alpha$, i.e.,

$\sigma=\alpha\Delta P_{D_{\mu}}.$ (18)

Having defined the mean and standard deviation of $\Delta P_{D}$ and $\Delta
P_{L}$, it should be clear now that the only decision to be made associated
with noncritical load shedding is $\Delta P_{D_{\mu}}$, meaning that
statistically the mean of noncritical load shedding equals to the decision
made by the system operator. However, for the realization in a single
islanding event, it is very likely for $\Delta P_{D}$ to deviate from the
decision variable $\Delta P_{D_{\mu}}$ characterized by its standard deviation
$\sigma$.

###  III-A Nadir Constraint Reformulation

Based on the method proposed in [17], the nadir constraint $\Delta
f(t_{n})\leq\Delta f_{\mathrm{lim}}$ can be converted into the following
nonlinear from:

$HR\geq\frac{\Delta P_{L}^{2}T_{d}}{4\Delta f_{\mathrm{lim}}}-\frac{\Delta
P_{L}T_{d}D_{0}}{4}+\frac{\Delta
P_{L}T_{d}\sum_{w\in\mathcal{W}}\gamma_{w}H_{s_{w}}^{2}}{4}.$ (19)

Since $\sum_{w\in\mathcal{W}}\gamma_{w}H_{s_{w}}^{2}$ is much less than
$D_{0}$, the $\Delta P_{L}$ in the last term of (19) is set to be a constant
($\Delta P_{L}^{\mathrm{max}}$ for conservativeness). As a result, (19) can be
rewritten as follows:

$HR-\underbrace{\frac{\Delta
P_{L}^{\mathrm{max}}T_{d}\sum_{w\in\mathcal{W}}\gamma_{w}H_{s_{w}}^{2}}{4}}_{c}\geq\frac{T_{d}}{4\Delta
f_{\mathrm{lim}}}\Delta P_{L}^{2}-\frac{T_{d}D_{0}}{4}\Delta P_{L}.$ (20)

Equation (20) can be viewed as a quadratic inequality of $\Delta P_{L}$. Since
$T_{d}/(4\Delta f_{\mathrm{lim}})>0$, it is equivalent to:

$\displaystyle\Delta P_{L}\in\Bigg{[}$
$\displaystyle\underbrace{\frac{D_{0}\Delta f_{\mathrm{lim}}}{2}-\frac{\Delta
f_{\mathrm{lim}}}{T_{d}}\sqrt{\frac{T_{d}^{2}D_{0}^{2}}{4}+\frac{4T_{d}}{\Delta
f_{\mathrm{lim}}}(HR-c)}}_{\Delta\underline{P}_{L}},$
$\displaystyle\underbrace{\frac{D_{0}\Delta f_{\mathrm{lim}}}{2}+\frac{\Delta
f_{\mathrm{lim}}}{T_{d}}\sqrt{\frac{T_{d}^{2}D_{0}^{2}}{4}+\frac{4T_{d}}{\Delta
f_{\mathrm{lim}}}(HR-c)}}_{\Delta\bar{P}_{L}}\Bigg{]}$ (21)

Therefore, the distributionally robust nadir chance constraint is formulated
as follows:

$\min_{\mathbf{D}\in\mathcal{P}}\mathrm{Pr}\Big{\\{}\Delta
P_{L}\in\big{[}\Delta\underline{P}_{L},\Delta\bar{P}_{L}\big{]}\Big{\\}}\geq\eta.$
(22)

Here only the cases of generation loss are considered, i.e., $\Delta P_{L}>0$,
therefore, the probability of $\Delta P_{L}$ being negative is zero. Moreover,
it can be derived from (III-A) that $\Delta\underline{P}_{L}<0$ always holds.
Hence, (22) is equivalent to:

$\min_{\mathbf{D}\in\mathcal{P}}\mathrm{Pr}\Big{\\{}\Delta
P_{L}\leq\Delta\bar{P}_{L}\Big{\\}}\geq\eta.$ (23)

By applying Chebyshev inequality [27], the nonconvex constraint (23) can be
reformulated as follows:

$\Delta\bar{P}_{L}\geq\Delta
P_{L_{\mu}}+\underbrace{\sqrt{\frac{\eta}{1-\eta}}}_{\xi}\sigma.$ (24)

Substituting the expression of $\Delta\bar{P}_{L}$ into (24) yields:

$\displaystyle HR\geq$ $\displaystyle\frac{T_{d}}{4}\left[\frac{(\Delta
P_{L_{\mu}}+\xi\sigma)^{2}}{\Delta f_{\mathrm{lim}}}-D_{0}(\Delta
P_{L_{\mu}}+\xi\sigma)\right]$ $\displaystyle+$ $\displaystyle\frac{\Delta
P_{L}^{\mathrm{max}}T_{d}\sum_{w\in\mathcal{W}}\gamma_{w}H_{s_{w}}^{2}}{4}.$
(25)

Introduce ancillary variables $x_{1}$ such that:

$\displaystyle x_{1}^{2}$ $\displaystyle=\frac{(\Delta
P_{L_{\mu}}+\xi\sigma)^{2}}{\Delta f_{\mathrm{lim}}}-D_{0}(\Delta
P_{L_{\mu}}+\xi\sigma)$ $\displaystyle=\underbrace{\frac{\Delta
P_{L_{\mu}}+\xi\sigma}{\sqrt{\Delta
f_{\mathrm{lim}}}}}_{x_{2}}\left(\frac{\Delta
P_{L_{\mu}}+\xi\sigma}{\sqrt{\Delta
f_{\mathrm{lim}}}}-\underbrace{\sqrt{\Delta
f_{\mathrm{lim}}}D_{0}}_{d}\right).$ (26)

It can be proved that $x_{2}\geq d$ always holds given a small system damping.
Hence, (26) is a well-defined real value constraint. The nadir constraint
(III-A) can be thus rewritten as a SOC form:

$HR\geq\frac{T_{d}}{4}x_{1}^{2}+\frac{\Delta
P_{L}^{\mathrm{max}}T_{d}\sum_{w\in\mathcal{W}}\gamma_{w}H_{s_{w}}^{2}}{4}$
(27)

However, the nonconvex constraint (26) can not be included in the MISOC
optimization directly. A set of linear constraints are used to approximate
this relationship conservatively. Depending on the ratio of $x_{2}$ to $d$,
the relationship between $x_{1}$ and $x_{2}$ described by (26) can be piece-
wise characterized by $N$ linear expression:

$\displaystyle x_{1}=a_{n}$ $\displaystyle x_{2}+b_{n}D_{0},\,\,\forall
n\in\mathcal{N}$ $\displaystyle if\,\;k(n-1)+1\leq\frac{x_{2}}{d}<kn+1$ (28a)
$\displaystyle x_{1}=a_{N}$ $\displaystyle x_{2}+b_{N}D_{0},$ $\displaystyle
if\,\;K+1\leq\frac{x_{2}}{d}$ (28b)

where $k$ defines the step size of $x_{2}$ in the first $N-1$ linear
constraints, i.e., $n\in\mathcal{N}=\\{1,2,...,N-1\\}$:

$k=\frac{K}{N-1}$ (29)

with $K$ being a predefined constant. Theoretically, the ratio of $x_{2}$ to
$d$ can be vary large. To reduce the value of $N$, a single linearized
constraint is applied if $x_{2}$ is large enough compared to $d$, i.e., $K\leq
x_{2}/d$ as shown in (28). The coefficient $a_{n}$ and $b_{n}$ can be
calculated as:

$\displaystyle a_{n}$
$\displaystyle=\derivative{x_{1}}{x_{2}}\Bigr{|}_{x_{2}=(kn+1)D_{0}}=\frac{2nk+1}{2\sqrt{n^{2}k^{2}+nk}},\,\,\forall
n\in\mathcal{N}$ (30a) $\displaystyle b_{n}$
$\displaystyle=\frac{-nk-1}{2\sqrt{n^{2}k^{2}+nk}},\,\,\forall
n\in\mathcal{N}$ (30b) $\displaystyle a_{N}$ $\displaystyle=1$ (30c)
$\displaystyle b_{N}$ $\displaystyle=-0.5$ (30d)

However, the $N$ linear constraints defined in (28) cannot be included in the
optimization simultaneously. Instead, only one of them needs to hold while
others should be relaxed depending on the relationship between $x_{2}$ and
$d$. Therefore, $N$ binary variables, $z_{n},\,\forall
n\in\mathcal{N}\cup\\{N\\}$ are introduced to indicate to which interval
$x_{2}$ belongs:

$\displaystyle z_{n\in\mathcal{N}}$
$\displaystyle=\begin{cases}1&if\,\,dk(n-1)\leq x_{2}<dkn\\\
0&\mathrm{otherwise}.\end{cases}$ (31a) $\displaystyle z_{N}$
$\displaystyle=\begin{cases}1&if\,\,Kd\leq x_{2}\\\
0&\mathrm{otherwise}.\end{cases}$ (31b)

Equation (31a) can be rewrite in the following form by defining ancillary
binary variables $z_{n_{1}},\,z_{n_{2}},\,\forall n\in\mathcal{N}$:

$\displaystyle z_{n_{1}}$ $\displaystyle=\begin{cases}1&if\,\,dk(n-1)\leq
x_{2}\\\ 0&\mathrm{otherwise}\end{cases}$ (32a) $\displaystyle z_{n_{2}}$
$\displaystyle=\begin{cases}1&if\,\,x_{2}<dkn\\\
0&\mathrm{otherwise}\end{cases}$ (32b) $\displaystyle z_{n}$
$\displaystyle=z_{n_{1}}+z_{n_{2}}-1.$ (32c)

As a result, the conditional constraints (32) can be transformed into linear
constraints $\forall n\in\mathcal{N}$:

$\displaystyle 0<x_{2}-dkn+Mz_{n_{1}}\leq M$ (33a) $\displaystyle 0\leq
dk(n-1)-x_{2}+Mz_{n_{2}}<M$ (33b) $\displaystyle z_{n}=z_{n_{1}}+z_{n_{2}}-1$
(33c) $\displaystyle z_{n_{1}},\,z_{n_{2}}\,\in\\{0,1\\}$ (33d)

where $M>0$ is a sufficiently large constant. Similarly, linear constraints
for $Z_{N}$ can be expressed as follows:

$\displaystyle 0<x_{2}-Kd+Mz_{N}\leq M$ (34a) $\displaystyle
z_{N}\,\in\\{0,1\\}$ (34b)

Based on the interval indicators $z_{n}$, the equality constraints (28) are
relaxed as follows to ensure feasibility:

$x_{1}\geq a_{n}x_{2}+b_{n}d+(z_{n}-1)M^{\prime},\;\;\;\;\forall
n\in\mathcal{N}\cup\\{N\\}$ (35)

where $M^{\prime}>0$ is a large enough constant. It should be noted that the
equality in (35) will be automatically obtained during the optimization
process since $HR$ term in the nadir constraint (27) is positively correlated
to the objective function. As a result, the original distributionally robust
nadir chance constraint (III-A) is now reformulated into the SOC form:
(27)(33)(34)(35).

###  III-B RoCoF and steady-state Constraints Reformulation

Based on the expressions derived in (12) and (16), the distributionally robust
frequency constraints of RoCoF and steady-state can be formulated as follows:

$2H\Delta\dot{f}_{\mathrm{lim}}\geq\Delta P_{L_{\mu}}+\xi\sigma$ (36)

$R+\Delta P_{C}+(D_{0}-\sum_{w\in\mathcal{W}}\gamma_{w}H_{s_{w}}^{2})\Delta
f^{ss}_{\mathrm{lim}}\geq\Delta P_{L_{\mu}}+\xi\sigma$ (37)

with $\dot{f}_{\mathrm{lim}}$ and $\Delta f^{ss}_{\mathrm{lim}}$ being the
pre-specified limits for maximum instantaneous RoCoF and steady-state
frequency deviation. The nonlinear term associated with $H_{s_{w}}^{2}$ can be
effectively linearized as demonstrated in [17].

##  IV MISOCP-based Microgrid Scheduling

In this section, a two-stage stochastic microgrid scheduling model is
introduced to determine the optimal generator dispatch, wind/PV curtailment
and load shedding with frequency security constraints. The relationship
between the microgrid scheduling problem and the frequency control is
demonstrated through Fig. 2. During normal operation (grid-connected mode),
the microgrid operates according to the results obtained from the scheduling
problem with optimal operating conditions and the optimal frequency responses
updated in each hour, as indicated by the blue and yellow areas respectively
in the figure below. Notably, the commands related to the frequency services
would only lead to the parameter updates in the associated controllers. Those
services would not be triggered unless an islanding event is detected, which
can be achieved by monitoring the main breaker at the PCC or measuring the
RoCoF. Due to the frequency constraints in the microgrid scheduling model, the
frequency security after an islanding event at any time can be guaranteed
through a most cost-efficient way.

Figure 2: Relationship between microgrid scheduling and frequency control.

Consider a microgrid with a set of SGs units $g\in\mathcal{G}$ and loads
$l\in\mathcal{L}$. The generation units are further categorized into two
groups, i.e. $\mathcal{G}=\mathcal{G_{\mathrm{1}}\cup G_{\mathrm{2}}}$
representing the sets of fast and slow generators. Wind, PV and storage units
are represented by $w\in\mathcal{W}$, $m\in\mathcal{M}$ and $b\in\mathcal{B}$
respectively. The uncertainties of renewable generation and demand in the
microgrid are managed by the two-stage decision process. The unit commitment
decisions are made in the first stage except for the fast-start generators,
before the uncertainty is realized. Once most uncertain inputs (demand and
renewable generation) are realized, the power outputs of committed units as
well as the fast-start generators are decided to meet the load [28, 29]. The
two-stage decision process in power system operations makes it natural to
formulate the scheduling problem as a multi-stage stochastic program. Based on
the stochastic multi-temporal method proposed in [30], the stochastic
scheduling problem can be formulated as follows.

###  IV-A Objective Function

The objective of the scheduling problem is to minimize the microgrid average
operation cost for all scenarios ($\forall s\in\mathcal{S}$) along with
considered time horizon $t\in\\{0,1,...,T\\}$:

$\begin{split}\min\sum_{s\in\mathcal{S}}&\sum_{t\in
T}\pi_{s}(\sum_{g\in\mathcal{G}}c_{g}^{SU}z_{t,s,g}+\Delta
t(\sum_{g\in\mathcal{G_{\mathrm{1}}}}c_{g}^{R1}y_{t,s,g}+\\\
&\sum_{g\in\mathcal{G_{\mathrm{2}}}}c_{g}^{R2}p_{t,s,g}+\sum_{l\in\mathcal{L}}c^{VOLL}(p_{t,s,l}^{c}+(q_{t,s,l}^{c})^{2})))\end{split}$
(38)

where $\pi_{s}$ is the probability associted with scenario $s$; $c_{g}^{SU}$,
$c_{g}^{R1}/c_{g}^{R2}$ and $c^{VOLL}$ refer to start-up costs, running costs
of fixed/flexible generators and the value of lost load (VOLL); $z_{t,s,g}$
and $y_{t,s,g}$ are binary variables of generator $g$ at time step $t$ in
scenario $s$ with $1/0$ indicating starting up/not and on/off; $p_{t,s,g}$ and
$p_{t,s,l}^{c}/q_{t,s,l}^{c}$ denote the active power produced by generators
and active/reactive load shedding.

###  IV-B Constraints

The traditional microgrid scheduling constraints that related to the generator
operation, wind/PV curtailment and load shedding are omitted here. [31] can be
referred for more details.

####  IV-B1 Constraints of battery storage system

$\displaystyle\eqref{Hv_lim},\,\eqref{Pc_lim},\;\;\;\;\;\forall t,s,b$ (39a)
$\displaystyle\mathrm{SoC}_{t,s,b}E_{c,b}=\mathrm{SoC}_{t-1,s,b}E_{c,b}+\eta_{b}p_{t,s,b}\Delta
t,\;\;\;\;\;\forall t,s,b$ (39b)
$\displaystyle\mathrm{SoC}_{\mathrm{min}}\leq\mathrm{SoC}_{t,s,b}\leq\mathrm{SoC}_{\mathrm{max}},\;\;\;\;\;\forall
t,s,b$ (39c)
$\displaystyle\mathrm{SoC}_{0,s,b}=\mathrm{SoC}_{T,s,b},\;\;\;\;\;\forall
t,s,b.$ (39d)

The power injection from the battery storage system to the microgrid is
confined in (39a) by the upper bound of the charging and discharging rate with
$p_{b}$ in the original equations being replaced by $p_{t,s,b}$. The battery
state of charge is quantified by (39b) with the charging/discharging
efficiency $\eta_{b}$. (39c) imposes the upper and lower limits on the SoC of
the storage devices. The SoC at the end of the considered time horizon is set
to be a pre-specified value being equal to its initial value as in (39d).

####  IV-B2 Constraints of AC power flow and power balance

$\displaystyle W_{t,s,ij}W_{t,s,ij}^{*}\leq
W_{t,s,ii}W_{t,s,jj},\;\;\;\;\;\forall t,s,i,j$ (40a) $\displaystyle
V_{\mathrm{min},i}^{2}\leq W_{t,s,ii}\leq
V_{\mathrm{max},i}^{2},\;\;\;\;\;\forall t,s,i$ (40b) $\displaystyle
p_{t,s,i}^{G}=\sum_{\Omega_{g-i}}p_{t,s,g}+\sum_{\Omega_{w-i}}p_{t,s,w}$
$\displaystyle\quad\quad\quad\quad\quad+\sum_{\Omega_{m-i}}p_{t,s,m}+\sum_{\Omega_{b-i}}p_{t,s,b},\;\;\;\;\;\forall
t,s,i$ (40c) $\displaystyle
q_{t,s,i}^{G}=\sum_{\Omega_{g-i}}q_{t,s,g}+\sum_{\Omega_{w-i}}q_{t,s,w}$
$\displaystyle\quad\quad\quad\quad\quad+\sum_{\Omega_{m-i}}q_{t,s,m}+\sum_{\Omega_{b-i}}q_{t,s,b},\;\;\;\;\;\forall
t,s,i$ (40d) $\displaystyle
p_{t,s,i}^{D}=\sum_{\Omega_{l-i}}p_{t,s,l}-\sum_{\Omega_{l-i}}p_{t,s,l}^{c},\;\;\;\;\;\forall
t,s,i$ (40e) $\displaystyle
q_{t,s,i}^{D}=\sum_{\Omega_{l-i}}q_{t,s,l}-\sum_{\Omega_{l-i}}q_{t,s,l}^{c},\;\;\;\;\;\forall
t,s,i$ (40f) $\displaystyle
p_{t,s,i}=p_{t,s,i}^{G}-p_{t,s,i}^{D},\;\;\;\;\;\forall t,s,i$ (40g)
$\displaystyle q_{t,s,i}=q_{t,s,i}^{G}-q_{t,s,i}^{D},\;\;\;\;\;\forall t,s,i$
(40h) $\displaystyle
p_{t,s,i}=\sum_{ij\in\mathcal{R}}p_{t,s,ij},\;\;\;\;\;\forall t,s,i$ (40i)
$\displaystyle
q_{t,s,i}=\sum_{ij\in\mathcal{R}}q_{t,s,ij}-\mathrm{Im}{(W_{t,s,ii}\mathrm{j}Y_{i,sh})},\;\;\;\;\;\forall
t,s,i$ (40j) $\displaystyle
p_{t,s,ij}+\mathrm{j}q_{t,s,ij}=W_{t,s,ii}{Y_{i,sh}^{*}}$
$\displaystyle\quad\quad\quad\quad\quad-(W_{t,s,ii}-W_{t,s,ij})y_{ij}^{*},\;\;\;\;\;\forall
ij\in\mathcal{R},t,s$ (40k) $\displaystyle p_{t,s,ij}^{2}+q_{t,s,ij}^{2}\leq
S_{\mathrm{max},ij}^{2},\;\;\;\;\;\forall ij\in\mathcal{R},t,s$ (40l)

Equations (40a) and (40b) are second-order cone constraints of voltages [32]
where $W_{t,s,ij}=V_{t,s,i}V_{t,s,j}^{*},\,\forall t,s,i,j$;
$V_{t,s,i}/V_{t,s,j}$ are voltages at bus $i/j$ and
$V_{\mathrm{min},i}/V_{\mathrm{max},i}$ are minimum/maximum voltage at bus
$i$. Total active/reactive power generation and load at each bus are defined
in (40c)/(40d) and (40e)/(40f) with $\Omega_{g/w/m/b-i}$ and $\Omega_{l-i}$
being the set of synchronous/wind/PV/storage units and loads connected to bus
$i$. Note that the imported power from the main grid, $P_{t,s,im}/Q_{t,s,im}$
is include in $\Omega_{g-i}$ for simplicity. Power balance at each bus is
given by (40g) to (40j) where $p_{t,s,ij}/q_{t,s,ij}$ are active/reactive
power flow from bus $i$ to $j$ and $ij\in\mathcal{R}$ is the set of branches;
$Y_{i,sh}=Y_{j,sh}$ denotes shunt susceptances at both ends of the line. (40k)
and (40l) are the power flow and line rating constraints.

####  IV-B3 Frequency security constraints subsequent to islanding events

According to the derivation in Section III, (27)(33)(34)(35), (36) and (37)
are incorporated into the microgrid scheduling model as the frequency nadir,
RoCoF and steady-state constraints. Therefore, the optimal microgrid inertia
which includes both conventional and synthetic one, the PFR and the equivalent
loss of generation will be determined in the microgrid scheduling model to
ensure the minimum operational cost while maintaining the frequency
constraints.

Different operations are coordinated before, during and after islanding events
to ensure the frequency security. Before the islanding event, the microgrid
operates in grid-connected mode. All the operating points of the dispatchable
units and the imported power from the maingrid are set optimally according to
the results from the scheduling model as indicated by the blue area of the
above figure. These settings help to ensure the frequency security by
allocating proper reserve, system inertia and power exchange with the
maingrid.

At the time instant of an islanding event, the power exchange with the
maingrid becomes zero almost instantaneously leading to a step power
disturbance to the microgrid. This event can be detected with negligible
delays by monitoring the main breaker at PCC and measuring the RoCoF of the
frequency [14, 33].

After that, different frequency services instructed by the scheduling results,
begin to react automatically including the frequency response from SGs, the SI
provision from IBGs and the noncritical load shedding in order to facilitate
the frequency regulation. After around tens of seconds, the frequency reaches
and remains at the steady-state value until the recovery and restoration
processes. Note that the fault repair and microgrid recovery and restoration
processes are out of the scope of the proposed model.

##  V Case Studies

Figure 3: Modified 14-bus microgrid test system.

In order to demonstrate the performance of the proposed distributionally
robust chance constrained microgrid scheduling model, case studies are carried
out through the modified IEEE 14-bus distribution system [34] as shown in Fig.
3. The optimization problem is solved in a horizon of 24 hours with the time
step being 1 hour. System parameters are set as follows: load demand
$P_{D}\in[160,300]\,\mathrm{MW}$, damping $D=0.5\%P_{D}/1\,\mathrm{Hz}$, PFR
delivery time $T_{d}=10\,\mathrm{s}$. The frequency limits of nadir, steady-
state value and RoCoF are set as: $\Delta f_{\mathrm{lim}}=0.8\,\mathrm{Hz}$,
$\Delta f_{\mathrm{lim}}^{\mathrm{ss}}=0.5\,\mathrm{Hz}$ and
$\Delta\dot{f}_{\mathrm{lim}}=0.5\,\mathrm{Hz/s}$. Dispatchable SGs are
installed at Bus 1,2 and 3 with a total capacity of $240\,\mathrm{MW}$. The
PV-storage system and wind turbines locate at Bus 6 and 8 respectively. The
parameters of battery devices are listed in Table I. The weather conditions
are obtained from online numerical weather prediction [35]. The MISOCP-base
optimization problem is solved by Gurobi (8.1.0) on a PC with Intel(R)
Core(TM) i7-7820X CPU @ 3.60GHz and RAM of 64 GB.

TABLE I: Parameters of Battery Storage Devices

$\boldsymbol{\mathrm{SoC}_{\mathrm{min}}}$ | $\boldsymbol{\mathrm{SoC}_{\mathrm{max}}}$ | $\boldsymbol{\eta_{b}}$ | $\boldsymbol{|\bar{P}^{(d)ch}|\,\mathrm{[MW]}}$ | $\boldsymbol{E_{c}\,\mathrm{[MWh]}}$
---|---|---|---|---
$15\%$ | $85\%$ | $0.9$ | $50$ | $150$

###  V-A Frequency nadir validation

Figure 4: Microgrid frequency evaluation after an islanding event.

Figure 5:  Histogram of frequency nadirs in 24-hour scheduling.

Figure 6: Frequency comparison: Analytical and EMT results.

The approximation of the frequency nadir constraint discussed in Section III-A
is assessed through dynamic simulation via Matlab/Simulink. The microgrid
frequency evaluation subsequent to an islanding event with and without (w/o)
the nadir constraint is illustrated in Fig. 4. The system operating conditions
are selected at an arbitrary hour based on the optimal scheduled results:
$P_{D}=162.7\,\mathrm{MW}$, $R=50.1\,\mathrm{MW}$, $\Delta
P_{L}=37.0\,\mathrm{MW}$, $H=86.0\,\mathrm{MWs/Hz}$. It can be observed that
the microgrid frequency decreases dramatically after an islanding event if the
nadir constraint is not implemented. Even though the steady-state frequency is
within the limit, the RoCoF and nadir constraints are violated. On the
contrary, once the nadir constraint is considered in the microgrid scheduling
model, all the frequency constraints can be maintained. The frequency nadir of
$-0.77\,\mathrm{Hz}$ shows a good approximation yet conservativeness. All the
other conditions present a similar frequency evaluation trend thus not being
covered. Instead, to demonstrate the robustness of the proposed method in
terms of effectiveness of the nadir constraints, the frequency nadir in each
hour of the one-day scheduling if an unintentional islanding event occur is
obtained through the dynamic simulation with the results depicted in Fig. 5.
It is observed from the histogram that all the frequency nadirs during the
24-hour scheduling are close to the boundary (-0.8 $\mathrm{Hz}$) with the
mean and standard deviation being $-0.7814\,\mathrm{Hz}$ and
$0.008\,\mathrm{Hz}$ respectively, indicating a good robustness of the
proposed method.

Additionally, the developed model is incorporated into the detailed EMT
simulation and analyzed for the test case with $P_{D}=199.6\,\mathrm{MW}$,
$R=57.0\,\mathrm{MW}$, $\Delta P_{L}=30.2\,\mathrm{MW}$,
$H=48.7\,\mathrm{MWs/Hz}$. As shown in Fig. 6, an unintentional islanding
event occurs at $t=0\,\mathrm{s}$. The analytical result represent the Center-
of-Inertia (CoI) frequency of the microgrid, whereas the 4 trajectories in the
EMT result represent the local frequencies at the generation buses (Bus 1, 2,
6 and 8). Note that the SG at Bus 3 is not online in this hour. High
oscillations depicted in the figure reflect the complexity of the EMT model at
hand, as well as the level of controller interaction characteristic of low
inertia system. It is observed that the frequency constraints in both cases
can be maintained and the EMT result stays close to the analytical one despite
a little mismatch after the frequency nadir, which is due to the approximation
of the SG model in the analytical derivation. For more detailed SG, VSM and WT
models, [16, 17] can be referred.

###  V-B Impact of islanding frequency constraints and SI from RESs

In this section, the influence of the frequency constraints subsequent to
microgrid islanding events as well as the value of SI are investigated. System
operation cost at different scheduling conditions is presented in Fig. 7 with
the cases defined as follows.

  * •

Base Case: Do not consider frequency dynamic constraints.

  * •

Case I: Consider frequency dynamic constraints, and SI is not allowed.

  * •

Case II: Consider frequency dynamic constraints, and SI is allowed.

Figure 7: Averaged cost at different operating conditions.

Figure 8: Microgrid imported power at various IBGs’ capacity.

Notably, the IBGs’ capacity refers to the total capacity of wind turbines
($P_{W}^{C}$) and PV systems ($P_{M}^{C}$) with $P_{W}^{C}/P_{M}^{C}=3/5$; To
avoid the PV and wind power curtailment due to the battery storage saturation,
the total battery capacity (${\bar{P}^{(d)ch}},\,E_{c}$) also varies with the
PV capacity, i.e., $\bar{P}^{(d)ch}:E_{c}:P_{M}^{C}=1:3:2$. It is observed
that the averaged system operation cost over 24 hours decreases along with the
increase of IBGs’ capacity in the system for all the three cases as more
energy is supplied by the RESs. In the base case, the system operation cost
always has the smallest value since the frequency dynamic constraints are not
considered. As a consequence, violations of RoCoF and nadir constraints would
be inevitable, should islanding events happen. For instance,
$\Delta\dot{f}_{\mathrm{max}}\in[-3.81,-1.70]\,\mathrm{Hz/s}$ with an average
of $-2.53\,\mathrm{Hz/s}$ and $\Delta
f_{\mathrm{max}}\in[-13.08,-6.89]\,\mathrm{Hz}$ with an average of
$-9.44\,\mathrm{Hz}$ are observed at IBGs’ capacity of 320 $\mathrm{MW}$. Once
the frequency dynamic constraints are included and SI is not allowed, the
averaged cost (blue curve) grows to maintain the frequency limits by dispatch
more partially-loaded SGs in the system for inertia provision only. The SI
provision from RESs (Case II) reduces the operational cost significantly. This
cost saving (the difference between Case I and II) becomes more obvious at
high IGB’s capacity where the cost of Case II is almost the same as that of
the base case, which highlights the effectiveness and value of SI provision
especially in high PE-penetrated microgrids. The computational time of each
optimization in different cases varies between $[59.14,169.84]\,\mathrm{s}$
with an average of $100.92\,\mathrm{s}$.

The averaged imported power from the main grid, $P_{i,avg}$ is also depicted
in Fig. 8. In the Base Case, the imported power starts to decrease after IBGs’
capacity becomes higher than about $110\,\mathrm{MW}$ since less energy is
needed from the main grid. If the frequency dynamic constraints are considered
(Case I), the microgrid cannot deal with the large disturbance without the SI.
Therefore, $P_{i,avg}$ is reduced by around a half compared to Base Case in
order to decrease the system disturbance level. In Case II, the system
available SI becomes higher as the IBGs’ capacity rises. Therefore, the
microgrid can withstand larger disturbance without violating the frequency
constraints, thus enabling more imported power from the main grid compared to
Case I, which also justifies the cost saving in Fig. 7.

The effects on the system SGs’ dispatch in different cases are investigated as
well with the results of 24 hours shown in Fig. 9 together with the demand and
total SI ($H_{r}$) profile. The IBGs’ capacity in the system is
$160\,\mathrm{MW}$. The implementation of the frequency dynamic constraints
(Case I) induces more power dispatched from SGs in almost all hours compared
to Base Case such that the power from the main grid could be reduced. With the
SI from RESs (Case II), more power can be supplied by the main grid and RESs
leading to a declined SG power. In addition, the total SI from RESs, varying
in the range of $[23-94]\,\mathrm{MWs/Hz}$ is also plotted, where its inverse
relationship with SG power in Case II is observed. In particular, during the
time of low net demand (i.e., $t\in[12,16]\,\mathrm{h}$), a significant amount
of SI is scheduled from RESs due to the lower inertia from SGs and vice versa.

Figure 9: Microgrid scheduling results within one day.

###  V-C Impact of uncertainty level of demand shedding during islanding
events

Figure 10: Averaged operation cost at different uncertainty levels.

In order to maintain the frequency constraints during unintentional islanding
events, noncritical load shedding is implemented to reduce the disturbance
magnitude. The IBGs’ total capacity is set to be $160\mathrm{MW}$. The
uncertainty level $\alpha$ associated with the noncritical load shedding as
defined in (18) is evaluated in this section. Its influence on the averaged
microgrid operation cost during 24 hours is depicted in Fig. 10 with different
confidence levels ($\eta=0.95$ and $0.90$). As expected, a higher uncertainty
level generates more operational cost since its effect on reducing the
disturbance becomes less reliable. Moreover, one can find that as the
confidence level is reduced from $0.95$ to $0.90$, the cost decreases by
approximate $10\%$, which highlights that the trade-off between the risk level
and microgrid operation cost needs to be well balanced. It is also worth
noticing that as the uncertainty level increases above some value (0.27 or
0.38), the microgrid operation cost becomes the same as the case where
noncritical load shedding is not implemented, represented by the dashed yellow
curve. Since the Chebyshev inequality is used in the deviation of nadir
constraints (24), which gives the lower bound regardless of the actual
distribution of the noncritical load shedding, conservative results are
obtained. Therefore, if this is the case in practice, system operators should
either pursue more knowledge of the load shedding distribution or decrease the
confidence level to achieve benefits in terms of microgrid operation cost
saving.

##  VI Conclusion and Future Work

This paper proposes a novel microgrid scheduling model enabling optimal
selection of the SI from IBGs while maintaining a minimum operational cost and
frequency constraints subsequent to an islanding event. Based on detailed
microgrid frequency dynamics and the state-of-art SI control schemes of the
IBGs, the frequency metrics subjected to an islanding event, which is modeled
as a step disturbance, are derived analytically. The uncertainty level
associated with the noncritical load shedding is modeled via an ambiguity set
without the knowledge of its specific distribution, leading to a
distributionally robust reformulation of the frequency constraints given a
certain confidence level. The nonlinear and nonconvex nadir constraint is
approximated by SOC relationship with the conservativeness and accuracy being
illustrated. An overall MISOCP-based optimization problem is formulated and
can be solved effectively using commercial solvers. Case studies demonstrate
the importance and necessity of islanding event consideration and the value of
SI provision from IBGs in terms of microgrid operation cost saving. The impact
of the uncertainty level of the noncritical load shedding is also
investigated.

The proposed model can be enhanced in several directions. Instead of assuming
the noncritical load can be shedding simultaneously, at the time instant when
the islanding event occurs, time delays in different types of load due to the
event detection and communications can be modeled in more detail. In addition,
the emulated damping from RESs can also be considered, providing more degrees
of freedom in the frequency regulation and microgrid scheduling.

## References

  * [1] M. Agrawal and A. Mittal, “Micro grid technological activities across the globe : A review,” 2011.
  * [2] A. Gholami, F. Aminifar, and M. Shahidehpour, “Front lines against the darkness: Enhancing the resilience of the electricity grid through microgrid facilities,” _IEEE Electrification Magazine_ , vol. 4, no. 1, pp. 18–24, 2016.
  * [3] S. M. Nosratabadi, R.-A. Hooshmand, and E. Gholipour, “A comprehensive review on microgrid and virtual power plant concepts employed for distributed energy resources scheduling in power systems,” _Renewable and Sustainable Energy Reviews_ , vol. 67, pp. 341 – 363, 2017.
  * [4] S. Parhizi, H. Lotfi, A. Khodaei, and S. Bahramirad, “State of the art in research on microgrids: A review,” _IEEE Access_ , vol. 3, pp. 890–925, 2015.
  * [5] W. Su, J. Wang, and J. Roh, “Stochastic energy scheduling in microgrids with intermittent renewable energy resources,” _IEEE Trans. Smart Grid_ , vol. 5, no. 4, pp. 1876–1883, 2014.
  * [6] J. Liu, H. Chen, W. Zhang, B. Yurkovich, and G. Rizzoni, “Energy management problems under uncertainties for grid-connected microgrids: A chance constrained programming approach,” _IEEE Trans. Smart Grid_ , vol. 8, no. 6, pp. 2585–2596, 2017.
  * [7] Z. Shi, H. Liang, S. Huang, and V. Dinavahi, “Distributionally robust chance-constrained energy management for islanded microgrids,” _IEEE Trans. Smart Grid_ , vol. 10, no. 2, pp. 2234–2244, 2019.
  * [8] A. Khodaei, “Resiliency-oriented microgrid optimal scheduling,” _IEEE Trans. Smart Grid_ , vol. 5, no. 4, pp. 1584–1591, 2014.
  * [9] Y. Li and K. Li, “Incorporating demand response of electric vehicles in scheduling of isolated microgrids with renewables using a bi-level programming approach,” _IEEE Access_ , vol. 7, pp. 116 256–116 266, 2019\.
  * [10] A. Khodaei, “Microgrid optimal scheduling with multi-period islanding constraints,” _IEEE Trans. Power Syst._ , vol. 29, no. 3, pp. 1383–1392, 2014.
  * [11] G. Liu, M. Starke, B. Xiao, X. Zhang, and K. Tomsovic, “Microgrid optimal scheduling with chance-constrained islanding capability,” _Electric Power Systems Research_ , vol. 145, pp. 197 – 206, 2017.
  * [12] B. Zhao, Y. Shi, X. Dong, W. Luan, and J. Bornemann, “Short-term operation scheduling in renewable-powered microgrids: A duality-based approach,” _IEEE Trans. Sustain. Energy_ , vol. 5, no. 1, pp. 209–217, 2014\.
  * [13] M. Hemmati, B. Mohammadi-Ivatloo, M. Abapour, and A. Anvari-Moghaddam, “Optimal chance-constrained scheduling of reconfigurable microgrids considering islanding operation constraints,” _IEEE Systems Journal_ , pp. 1–10, 2020.
  * [14] A. Gholami and X. A. Sun, “Towards resilient operation of multimicrogrids: An misocp-based frequency-constrained approach,” _IEEE Control Netw. Syst_ , vol. 6, no. 3, pp. 925–936, 2019.
  * [15] Y. Zhang, C. Chen, G. Liu, T. Hong, and F. Qiu, “Approximating trajectory constraints with machine learning microgrid islanding with frequency constraints,” _IEEE Trans. Power Syst._ , pp. 1–1, 2020.
  * [16] U. Markovic, Z. Chu, P. Aristidou, and G. Hug, “Lqr-based adaptive virtual synchronous machine for power systems with high inverter penetration,” _IEEE Trans. Sustain. Energy_ , vol. 10, no. 3, pp. 1501–1512, 2019.
  * [17] Z. Chu, U. Markovic, G. Hug, and F. Teng, “Towards optimal system scheduling with synthetic inertia provision from wind turbines,” _IEEE Trans. Power Syst._ , vol. 35, no. 5, pp. 4056–4066, 2020.
  * [18] F. Teng and G. Strbac, “Full stochastic scheduling for low-carbon electricity systems,” _IEEE Trans. Autom. Sci. Eng._ , vol. 14, no. 2, pp. 461–470, April 2017.
  * [19] H. Chávez, R. Baldick, and S. Sharma, “Governor rate-constrained opf for primary frequency control adequacy,” _IEEE Trans. Power Syst._ , vol. 29, no. 3, pp. 1473–1480, May 2014.
  * [20] B. Delfino, S. Massucco, A. Morini, P. Scalera, and F. Silvestro, “Implementation and comparison of different under frequency load-shedding schemes,” in _2001 Power Engineering Society Summer Meeting. Conference Proceedings_ , vol. 1, 2001, pp. 307–312 vol.1.
  * [21] N. N. A. Bakar, M. Y. Hassan, M. F. Sulaima, M. N. Mohd Nasir, and A. Khamis, “Microgrid and load shedding scheme during islanded mode: A review,” _Renewable and Sustainable Energy Reviews_ , vol. 71, pp. 161 – 169, 2017\.
  * [22] Y. R. Omar, I. Z. Abidin, S. Yusof, H. Hashim, and H. A. A. Rashid, “Under frequency load shedding (ufls): Principles and implementation,” in _2010 IEEE International Conference on Power and Energy_ , 2010, pp. 414–419.
  * [23] A. Rafinia, J. Moshtagh, and N. Rezaei, “Towards an enhanced power system sustainability: An milp under-frequency load shedding scheme considering demand response resources,” _Sustainable Cities and Society_ , vol. 59, p. 102168, 2020.
  * [24] L. Chang-Chien, L. N. An, T. Lin, and W. Lee, “Incorporating demand response with spinning reserve to realize an adaptive frequency restoration plan for system contingencies,” _IEEE Trans. Smart Grid_ , vol. 3, no. 3, pp. 1145–1153, 2012.
  * [25] Y. Dong, X. Xie, K. Wang, B. Zhou, and Q. Jiang, “An emergency-demand-response based under speed load shedding scheme to improve short-term voltage stability,” _IEEE Trans. Power Syst._ , vol. 32, no. 5, pp. 3726–3735, 2017.
  * [26] H. Mortaji, S. H. Ow, M. Moghavvemi, and H. A. F. Almurib, “Load shedding and smart-direct load control using internet of things in smart grid demand response management,” _IEEE Trans. Ind. Appl._ , vol. 53, no. 6, pp. 5155–5163, 2017.
  * [27] G. Calafiore and L. Ghaoui, “On distributionally robust chance-constrained linear programs,” _Journal of Optimization Theory and Applications_ , vol. 130, pp. 1–22, 01 2006.
  * [28] A. J. Conejo, M. Carrión, and J. M. Morales, _Decision making under uncertainty in electricity markets_. Springer, Boston, MA, 2010.
  * [29] P. A. Ruiz, C. R. Philbrick, E. Zak, K. W. Cheung, and P. W. Sauer, “Uncertainty management in the unit commitment problem,” _IEEE Trans. Power Syst._ , vol. 24, no. 2, pp. 642–651, 2009.
  * [30] F. Teng and G. Strbac, “Assessment of the role and value of frequency response support from wind plants,” _IEEE Trans. Sustain. Energy_ , vol. 7, no. 2, pp. 586–595, April 2016.
  * [31] Z. Chu, M. Zhao, and F. Teng, “Modelling of dynamic line rating in system scheduling: A misocp formulation,” in _2020 IEEE Power Energy Society General Meeting (PESGM)_ , 2020, pp. 1–5.
  * [32] J. A. Taylor, _Convex Optimization of Power Systems_. Cambridge University Press, 2015.
  * [33] A. Sharma and R. Sunitha, “Unintentional islanding detection in microgrid,” in _2017 International Conference on Energy, Communication, Data Analytics and Soft Computing (ICECDS)_ , 2017, pp. 2519–2523.
  * [34] “Power Systems Test Case Archive,” 1993, accessed: 2020-10-23. [Online]. Available: https://labs.ece.uw.edu/pstca/pf14/pg\\_tca14bus.htm
  * [35] “Met Office. Weather Forecasts,” 2020, accessed: 2020-10-23. [Online]. Available: https://www.metoffice.gov.uk/public/weather