> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)

1

Risk-limiting Load Restoration for Resilience
Enhancement with Intermittent Energy Resources
Zhiwen Wang, Student Member, IEEE, Chen Shen, Senior Member, IEEE, Yin Xu, Member, IEEE,
Feng Liu, Member, IEEE, Xiangyu Wu, Student Member, IEEE, Chen-Ching Liu, Fellow, IEEE
Abstract—Microgrids are resources that can be used to
restore critical loads after a natural disaster, enhancing
resilience of a distribution network. To deal with the
stochastic nature of intermittent energy resources, such as
wind turbines (WTs) and photovoltaics (PVs), many
methods rely on forecast information. However, some
microgrids may not be equipped with power forecasting
tools. To fill this gap, a risk-limiting strategy based on
measurements is proposed. Gaussian mixture model (GMM)
is used to represent a prior joint probability density
function (PDF) of power outputs of WTs and PVs over
multiple periods. As time rolls forward, the distribution of
WT/PV generation is updated based the latest measurement
data in a recursive manner. The updated distribution is used
as an input for the risk-limiting load restoration problem,
enabling an equivalent transformation of the original
chance constrained problem into a mixed integer linear
programming (MILP). Simulation cases on a distribution
system with three microgrids demonstrate the effectiveness
of the proposed method. Results also indicate that
networked microgrids have better uncertainty management
capabilities than stand-alone microgrids.
Index Terms— Resilience; load restoration; microgrids;
solar power; wind power, probabilistic distribution

NOMENCLATURE
Most symbols used in this paper are listed below for ease of
reference; others are explained following their first appearance.
Indices
k,t
Indices for periods. And, t=k+1,…,K
g
v
l
m
Parameters
K
G
W
S
L
M
ωm, μm, Σm
Ω

Indices for generators and energy
storage systems (ESSs)
Indices for WTs and PVs
Indices for loads
Indices for Gaussian components
Number of periods
Number of generators and ESSs
Number of WTs
Number of PVs
Number of loads
Number of Gaussian components
Parameters of the mth component of a
GMM
Ω={ωm, μm, Σm |m=1,…,M}

min
pmax
g , pg

Maximum and minimum power outputs
of generator g
cl
Priority weight of load l
τ
Length of each period. τ=1 hour
rupg , rdng
Ramp-up and ramp-down limits of
generator g
pt,l
Power of load l in period t
ENg
Generation resource of generator g
ECg
Capacity of ESS g
ρc, ρd
Charging and discharging efficiency of
ESS g
Random Variables
X vt
power output of WT/PV v in period t
Decision Variables
pt,g
Scheduled power output of generator g
in period t
ch
dch
Pt,g , Pt,g
Scheduled charging/discharging power
of ESS g in period t
ut,l
Status of load l to determine in period t.
ut,l = 1 if load l restored in period t; ut,l
= 0, otherwise
I.

INTRODUCTION

A. Motivation
xtreme weather events, such as flooding and hurricanes,
cause severe damages to power systems from time to time.
For example, hurricane Sandy hit the east coast of U.S. in 2012,
causing power service interruption for 7.5 million people [1]. In
2008, a snow storm in Southern China left over 14 million
households out of power for days. To help power systems
defend against major disasters, the concept of “resilience” is
proposed. Resilience refers to the ability of a power system that
can withstand low-probability high-impact events and quickly
recover customer service from major disasters [2]-[4].
Distributed energy resources (DERs), energy storage systems
(ESSs), as well as microgrids (MGs), are valuable resources
that can enhance resiliency of a distribution system, as they can
be used to restore critical loads when the utility power is not
available.

E

B. Literature Review
Much research has been conducted on resilience-oriented
load restoration problems [5]-[7]. A three-stage restoration
method for distribution networks is proposed in [5] that can be

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
used to maximize the restored loads with available DERs.
Authors in [6] study an automatic reconfiguration algorithm of
microgrids with dynamical boundaries, considering transients
of line faults and generators. In [7], a method is proposed to
transform a load restoration problem into a linear integer
programming “maximum coverage,” resulting in an efficient
solution. However, these studies do not fully take renewable
DERs into consideration. One prominent feature of renewables
is the stochastic nature of generation. If uncertainties are not
fully considered in the load restoration procedure, undesirable
events, e.g., load shedding, over/under-frequency issues, may
happen. To reduce these risks, there is a great need for a
systematic method to manage uncertainties introduced by
renewable DERs.
There are several popular methodologies in the literature that
address the renewable uncertainty issue. The first approach is
using a forecast value to replace the actual wind/solar power.
For example, authors in [8] propose an operation and selfhealing strategy, regarding the power outputs of DERs as
deterministic forecast values. In a relevant study [9], based on
short-term forecasts, a model predictive control approach is
proposed to enhance power system resilience with microgrids.
These studies only consider forecasted scenarios. However, as
the actual values of renewable generation may deviate far from
the forecasts, these methods may not be appropriate.
The second approach is called “scenario-based stochastic
optimization”. This method generates many possible scenarios
with different probabilities, and determines an optimal strategy
that can minimize an expected cost in those scenarios. In [10],
forecast errors of wind/solar power are modeled by Beta
distributions. By running Monte Carlo simulation, one can
generate a number of scenarios of DG outputs. To reduce the
computational burden, the back forward reduction method is
adopted to retain representative scenarios. In a relevant study
[11], a two-stage stochastic programming model is proposed, in
which uncertainty scenarios are generated by auto-regressivemoving-average (ARMA) technique based on weather forecast
information. Although the scenario-based stochastic
optimization is able to consider different possibilities of
uncertainties, there is a drawback. In order to keep the
information of uncertainties as intact as possible, the number of
scenarios has to be large. As a consequence, the computation of
the approach becomes intractable.
The third approach is robust optimization (RO). In [12],
robust optimization is proposed to handle renewable and load
uncertainties, in which uncertain parameters are modeled by
convex and bounded intervals, e.g., ± 20% of installed
capacities. In [13], a two-stage robust optimization model is
proposed to coordinate the hardening and DER allocation.
Based on weather forecast, the spatial and temporal dynamics
of an uncertain natural disaster, e.g., a hurricane, are captured
by a multi-stage and multi-zone uncertainty set. Because RO
can be transformed into a single-level min (or max) problem,
the computation of RO is relatively inexpensive compared with
the scenario-based approach. However, as the optimal solution
of RO is always obtained at an extreme point, i.e., the worstcase scenario, RO suffers from the inherent conservativeness,

2

which decreases the overall efficiency of an optimal strategy.
In recent years, an important methodology that is developed
to deal with uncertainties, called “risk-limiting dispatch (RLD)”,
has drawn more and more attention [14]-[19]. RLD aims at
finding an optimal strategy with a minimum cost to restrain
risks brought by uncertainties, e.g., the insufficient generation
risk. In RLD, power delivered at a particular period is
determined in multiple decision stages. In each decision stage,
distributions of renewable uncertainties are continuously
updated as conditional distributions of the latest forecasts, e.g.,
the forecasted wind speed from a numerical weather prediction
(NWP) system, wind power forecast values from a forecasting
tool. Thereafter, optimal decisions are made to limit risks based
on the newly collected information, i.e., decision variables
depend on the prediction information. If more precise
information of uncertainties is obtained in future decision
stages, operators have a chance to take recourse actions to
remedy previous decisions. At the moment the uncertainty is
realized, the risk is restrained within a tolerance. As the
accuracy of descriptions of uncertainties increases with newly
collected information, it is less conservative. The original idea
of RLD is proposed in [14]. In [15], an analytical solution of a
simple RLD is obtained based on dynamic programming
principles. The congestion issue is discussed in [16], while the
wind power ramping issue is addressed in [17]. The
applications of RLD in unit commitment and economical
dispatch can be found in [18], [19]. The scope of this paper is a
load restoration problem with microgrids in the distribution
system level. A major concern on the application of RLD in the
distribution level is the absence of power forecasting tools,
which may result in infeasibility of the forecast-based RLD.
Existing methods can be categorized into two groups:
(1) The first one assumes that there is continually updated
forecast information, e.g., [8]-[11], [13]-[19].
(2) The second one assumes that uncertainty models are fixed.
For example, in [12], the uncertain parameters belong to a fixed
interval [-0.2 +0.2], according to personal experience or
historical data analysis. In a relevant study [21], renewable
uncertainties are modeled by fixed distributions.
In fact, both assumptions might be invalid: (1) in the
distribution level, forecasting tools may not be available; (2) In
a multiple-period optimization problem, uncertainty models of
WTs/PVs in different periods are not the same. Using fixed
models may result in conservative results.
In addition to the challenge arising from the limited forecast
information, there is one more challenge for RLD: the
correlation of renewables. In RLD, probabilistic constraints
associated with risks require arithmetical addition operations of
random variables, e.g., the power balance equation. When
random variables are correlated, “convolution technique” is no
longer applicable [20]. Therefore, it is difficult to transform the
probabilistic constraints into tractable forms.
C. Contributions
To address these important issues, this paper proposes a risklimiting load restoration strategy based on observations of
renewable uncertainties. Relative to the state-of-the-art, the

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)

D. Organization
The remainder of the paper is organized as follows: Section
II provides the problem formulation; Section III discusses a
method to update distributions of uncertainties based on
observations. In Section IV, the original problem is transformed
into an MILP. Case study results are presented in Section V.
Conclusion and limitations are provided at the end of this paper.
II. PROBLEM FORMULATION
This section provides the formulation of the risk-limiting
load restoration problem. First, assumptions are provided and
justified. Then, the resilience metric used in this paper is
elaborated. A framework is given to clarify the computational
procedure, followed by a full description of the constraints of
the problem. Finally, challenges of the problem are discussed.
A. Assumptions
(1) There is a centralized microgrid controller that collects
information and optimizes the operation. Intelligent electronic
devices (IEDs) receive data and measurements from sensors
and pass them to the central controller. IEDs also transfer
control commands from the central controller to microgrid
devices, such as diesel generators, ESSs, and load breakers.
This communication structure can be implemented following
the standard IEC 61850 via TCP/IP [22]. In fact, it is commonly
used in microgrid testbeds [23], [24].
(2) In order to focus on volatile renewables, this paper
assumes that uncertainties come from solar and wind generation
only. Load uncertainty is not considered. This is justified as
follows. Note that only a small number of critical loads will be
restored after a major disaster. Since the critical loads are used
to maintain basic and minimum societal functions, e.g.,
hospitals and street lighting, their demand does not change
significantly during the restoration process [25]. Therefore,
load variations are much smaller than fluctuations of renewable
generation. Thus, they are neglected. In addition, if load varies
during the restoration process, the maximum load demand will
be used as the input to the restoration problem [25], [26], which
ensures that the load demand will be satisfied.
B. Resilience Metric
The conceptual resilience curve in Fig. 1 that is associated

with an extreme event is proposed in [27], and adopted in [25],
[28]-[30]. There are 7 well-defined states, in which different
measures can be taken to improve the resilience level, e.g.,
resource preparation in the pre-event state, and service
restoration using microgrids in the restorative state and postrestoration state.
F(t)

System Performance

contributions of this paper are threefold:
(1) The resilience-oriented customer service restoration with
intermittency renewables is formulated as a risk-limiting
decision-making problem.
(2) In modeling of uncertainties, a methodology is proposed
to recursively update distributions of renewables as latest
observations of WT/PV generation are available. Both wind and
solar power uncertainties are modeled in a universal manner.
The spatial-temporal correlation of uncertainties is considered.
(3) For the solution methodology, probabilistic constraints
consisting of correlated random variables are converted into
linear inequalities, based on a property of the updated
distribution. By doing so, the original problem can be solved
with commercial solvers.

3

Outage Duration
T0
Norm
State

Pre-Event
State

Event
Process

PostInfrastructure
Restoration Repair
State
Post-Event
Restorative
Degraded State State

t0

te

t pe

tr

t pr

t ir

tn

t

Fig. 1 A conceptual resilience curve

This paper is concerned with a load restoration strategy to
enhance resilience. Therefore, the period from tr to tir is of
primary interest, i.e., the restorative state and post-restoration
state. When an extreme event takes place, the utility power is
not available from tr to tir. Hence, the concerned period is
referred to as the “outage duration”. The resilience level is
evaluated by the integral of the system performance function
over the outage duration [27]-[30], i.e., the blue-shadowed area
in Fig.1:
tir
(1)
R
F  t  dt



tr

The system performance function adopted in this paper is
similar to those used in relevant resilience studies [25], [26],
[28], [29]: the weights of restored loads. This is reasonable
because more critical loads being restored indicates that the
distribution system can better restore service to customers after
an extreme event, i.e., the system is more resilient. Priority
weights of critical loads (hospital, police station, etc.,) should
be sufficiently greater than that of non-critical ones so that
critical loads can be served as a priority with limited available
generation capacities and resources.
Let the outage duration T0 be discretized into several subintervals τ. Thereafter, resilience index (1) can be computed as
follows:
(2)
R   cl  ut ,l 
t

l

That is to say, the resilience level is evaluated by the
cumulative service time to loads weighted by their priority over
the outage duration.
Remark 1: In some studies [26], [28], the active power of a
load is involved in the resilience metric. Then, equation (2)
becomes:
(3)
R   cl  Pt ,l  ut ,l 
t

l

Both (2) and (3) can well represent resilience of a system.
In fact, if one regards the weighted power in (3) as the weight
in (2), metrics (2) and (3) are equivalent. Regarding this, both
metrics can be applied to the proposed risk-limiting load
restoration method. This paper prefers metric (2) because it is

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
more intuitive than (3). For example, a hospital should be
identified as a critical load with a large weighted coefficient, as
it is a life-critical infrastructure after an extreme event. Suppose
that the hospital consumes a small amount of active power. If
metric (3) is adopted, the hospital might not be restored due to
a small value of the weighted power, while metric (2) avoids
the problem.
Remark 2: It should be noted that the resilience index is
quite different from the traditional “reliability” concept.
Reliability is associated with common events, most of them
making a low impact, e.g., a short circuit on a line. For such
events, most power sources stay connected. Transmission lines
and distribution feeders are available. Since the power system
facilities are intact, several scenarios will lead to high reliability,
e.g., scheduling sufficient reserves in preparation for any single
failure of equipment, and utilization of automatic reclosers to
re-energize a feeder after a temporary fault. By doing so,
reliability indices, e.g., the system average interruption
frequency index (SAIFI) and system average interruption
duration index (SAIDI), will be low.
In contrast, resilience is designed for extreme events that
rarely happen (low-probability high-impact events), during
which the power system infrastructure, e.g., transmission lines,
substations, automatic reclosers, and remote-control switches,
might be damaged. As a result, numerous power system
components are not available. Traditional measures that are
designed for reliability enhancement may not apply. For
example, if the upstream substation is damaged by a hurricane,
the automatic recloser cannot re-energize downstream loads. To
improve the resilience level, i.e., to quickly recover from a
major disaster, focus should be placed on how to pick up critical
loads as many as possible when the utility power is unavailable.
Measures for the resilience enhancement include utilization of
microgrids to restore critical service and hardening of
transmission and distribution facilities.
More discussions on the differences between resilience and
reliability can be found in [31].
C. Framework
RLD designed for transmission systems is based on two
assumptions. First, as time moves forward, the sensing system
can obtain more precise information on uncertainties. Second,
there is a market mechanism, where a power/energy transaction
enforced in earlier times is less costly. Therefore, system
operators using RLD scheme determine the power/energy
delivered at a certain period in multiple decision stages, trying
to achieve a tradeoff between making decisions earlier with
cheaper prices and making decisions later with more precise
information. In each decision stage, based on the latest
information, the risks brought by uncertainties are limited
within a tolerance, i.e., constraints associated with risk indices
are added into optimization models.
For a load restoration problem, the first assumption of RLD
still holds, while the second one may not be valid. An example
is given in Fig. 2. For the delivery period between 9:00 to 10:00,
suppose that one can make decisions in three stages, 7:00, 8:00,
and 9:00, which are called the 1st, 2nd, and final decisions,

4

respectively. Clearly, the uncertainty information of the
delivery period obtained at 9:00 is more precise than that
obtained at 8:00, which is more precise than that obtained at
7:00. However, unlike the market mechanism in the
transmission system level, the weighted coefficients of
resilience index (2) are independent of time, i.e., enforcing
decisions that are made earlier than 9:00 will not bring
additional reward. Therefore, a rational decision maker will not
enforce the decisions made at 7:00 and 8:00. That is to say, the
1st and 2nd decisions are for advisory purpose only; they are not
implemented.
7:00

8:00

9:00

...

10:00

16:00

17:00

Delivery Period

1st advisory
decision

2nd advisory
decision

t

Final
decision

Fig. 2 Illustration of advisory decisions and final ones for a delivery period.

The proposed risk-limiting load restoration is designed as
follows. The framework is illustrated in Fig. 3. For
demonstration, the time resolution is set to be 1 hour in this
paper. It can be changed to other frames, e.g. 10 min or 30 min,
if needed. The proposed method is an online tool, which is
executed periodically. In each implementation, the proposed
method aims to maximize resilience index (2), and determines
the set of loads to be restored, the generation schedule of
dispatchable generators, and the charging/discharging strategy
of ESSs. An optimal load restoration strategy is determined
over the optimization window, in which only decisions in the
first time interval are implemented to instruct
generators/ESSs/loads for the next period k+1 whereas the rest
are abandoned. There are multiple periods within a certain
optimization window. In each period, constraints associated
with risks are considered based on the latest uncertainty
information. Specifically, this paper uses observations of
uncertainties in the past periods to update distributions of
renewable generation in the future.
Out age Duration
7:00

8:00

9:00

10:00

...

16:00

17:00

7:00

8:00

9:00

10:00

...

16:00

17:00

...

16:00

17:00

k=0

k=1
Using observations to update
distributions of uncertainties
7:00

8:00

9:00

10:00

k=2
Using observations to update
distributions of uncertainties

Past intervals in
which unc erta inties
are observed

Next inter val in
which decisions
are implem ented

Future interva ls in
which decisions
are abandoned

Optimiz ation
window

Fig. 3 A framework of the risk-limiting load restoration scheme.

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
Remark 3: The periodic process in Fig. 3 is also known as
receding horizon control [32]. Compared with the method in
[32], the proposed method has two major differences:
(1) The method in [32] addresses critical issues of a load
restoration problem in deterministic scenarios. Renewable
uncertainties have not been fully considered [32]. This paper
mainly aims at handling uncertainties introduced by renewables
in the load restoration process.
(2) Authors in [32] point out that the information should be
updated periodically. One original contribution of this paper is
a recursive methodology in Section III that uses observations to
update distributions of renewable uncertainties.
D. Constraints
1) Risk-limiting constraints
● Since power outputs of WTs/PVs are stochastic, the
following adequacy requirement should be met in each
period:
W +S
L
G

Pr  Pt , g   X tv   ut ,l Pt ,l     t  k  1 (4)
v 1
l 1
 g 1


Equation (4) indicates that the probability of restored
load demands being satisfied with available resources
should be greater than a given confidence level. This is a
typical risk-limiting constraint [13]. In this paper, α is
90%.
● Similar to (4), an energy adequacy requirement is
formulated as follows:
T W S

v 
  EN g   1    X t  
 gESS

t  k 1 v 1
Pr 

T
L


   ut ,l Pt ,l


t  k 1 l 1


k

EN g  k  1  EN g  0    Pt , g

g  ESS

(5)

(6)

where ENg (0) is generation resource available at the
beginning; 𝑃̅𝑡̅,𝑔 is the power output of generator g in the
past k periods.
Note that when real time k rolls forward, distributions of
random variables Xt v in (4)(5) should be updated, which
are conditioned on the latest information.
2) Deterministic constraints
● If g is a dispatchable generator, e.g., a diesel, then the
following output limit, generation resource limit, and
ramping limit should be satisfied:
(7)
Pgmin  Pt ,g  Pgmax  t  k  1, g  ESS
T

t  k 1

  EN g  k  1

t ,g

rgdn  Pt ,g  Pt 1,g  rgup

 g  ESS

and maximum/minimum SOC are formulated as follows:
dch,max
0  Pt dch
, t  k  1, g  ESS (11)
, g   t , g Pg

 t , g Pgch,max  Pt ,chg  0 , t  k  1, g  ESS
SOCgmin  SOCt , g  SOCgmax

(12)
(13)

t  k  1, g  ESS

Note that the use of an ESS would involve decisions in
multiple periods, as the current SOC depends on the previous
SOC and the recent charging/discharging strategy. So, the SOC
transition of an ESS is formulated as follows:
1
ch
SOCt , g  SOCt 1, g   Pt dch
1, g  d +Pt 1, g  c / EC g
(14)
t  k  2, g  ESS
In the literature [8], the SOC at the end is required to be the
initial value. This paper does not include such a constraint,
because there are necessities to charge an ESS after the outage.





E. Challenges
To solve the problem (2)-(14), two critical issues need to be
addressed:
(1) Updating distribution of Xt v as time rolls forward;
(2) Dealing with probabilistic constraints (4) and (5), e.g.
computing the distribution of an aggregation of X t v from
multiple WTs/PVs in a single period, as well as over the whole
outage duration. Usually, one may use the “convolution
technique” [20] to compute the summation of Xt v . However,
Xt v is interdependent. Hence, the “convolution technique” does
not apply.
The two issues are addressed in Section III and IV,
respectively.
III. MODELING UNCERTAINTY

t 1

P

5

(8)

 t  k  2, g  ESS (9)

● If g is an ESS, then two auxiliary binary variables are
introduced: χt,g being 1 represents the discharging mode;
while γt,g being 1 for charging mode. Note that an ESS
cannot be in charging and discharging modes at the same
time [8]. That is :
(10)
t , g   t , g  1
t  k  1, g  ESS
Besides, limits on charging/discharging power output

In this section, features of intermittent energy resource
uncertainties are discussed, followed by details for modeling,
as well as updating, the uncertainties of WT/PV generation.
A. Features of uncertainties in microgrids
Wind and solar power are commonly used in microgrids.
Their uncertainties have three features:
(1) A distribution of wind power is different from that of
solar power. Usually, wind speed follows a Weibull distribution
[33]. Thereafter, the distribution of wind power can be obtained
by the “wind speed-power curve”, which is a modified Weibull
distribution with two spikes at 0.0 and 1.0 p.u. The PV
generation production is usually a function of solar irradiance
and temperature. Authors in [34] report that the irradiance in
Arizona can be described reasonably by a Beta distribution, and
the temperature is Gaussian. Then, the distribution of solar
power is obtained as a variant of a Beta distribution.
Alternatively, authors in [35] adopt the Saunier model to
evaluate the irradiation and obtain a different distribution of
solar generation. In the presence of both WTs and PVs, different
forms of distributions make it difficult to solve a load
restoration problem in a universal manner.
(2) Wind/solar power outputs have spatial-temporal
correlations. First, in the distribution system level, renewables
are geographically closely located. Hence, the spatial

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
dependence of power outputs from adjacent wind farms/solar
arrays is significant. Second, power outputs over adjacent
periods are also correlated. An illustration of the spatialtemporal correlation of two adjacent WTs over 10 hours is
shown in Fig. 4. The hourly data of wind power is from the
“wind integration data set” of National Renewable Energy
Laboratory (NREL) [36]. Different colors correspond to
different values of correlation coefficients. According to Fig. 4,
the correlation coefficients of the two WTs in one period
(spatial correlation), and two adjacent periods of one WT
(temporal correlation), are about 0.89-0.96, indicating strong
interdependences. Similar results can be obtained for PVs.

In this paper, a prior joint distribution of X is presented by a
GMM with an adjustable parameter set Ω={ωm, μm, σm ;
m=1,…,M} as follows:

f X  x    m N m  x ; μ m ,  m 
M

M


m 1

m

 1, m  0

(17)

9

10

(16)

m 1

1.00

N m  x ; μ m ,  m  :

7

8

0.95

5

6

0.90

3

4

0.85

2

0.80

1

power of WT2
over different periods

B. A prior distribution of uncertainties
Let a random variable Xk denote the actual power output of
an intermittent energy resource in period k. Then, a random
vector, X, for outputs over K periods is defined as:
T
(15)
X   X1
Xk
XK 

e



1
 x μm T m1  x μm 
2

 2π 

K

1/ 2

where ωm is the weight, Nm(·) is the mth multivariate Gaussian
component with a mean vector μm and a covariance matrix σm.

10

GMM
1st Component
2nd Component
3rd Component

3.5

9

0.70

7

8

3.0
0.65

2.5

1

2

0.55
0.50
1

2

3

4

5

6

7

8

9

power of WT1
over different periods

10

1

2

3

4

5

6

7

8

9

power of WT2
over different periods

10

Legend of
correlation
coefficient

Fig. 4 The spatial-temporal correlation matrix of two WTs over 10 hours

(3) Unlike transmission systems, some microgrids may not
be equipped with power forecasting tools due to cost
considerations. Hence, in some cases, forecast information for
WTs/PVs is unavailable, increasing the difficulty in decision
making of load restoration.
To overcome this hurdle, an intuitive idea is to use latest
observations to infer future uncertainties. Note that power
outputs in period k and k+1 are usually positively correlated. If
the power in k is observed and it is at a high value, there is a
good chance that power in k+1 is also at a high value. By doing
so, operators can update the inference and make decisions
correspondingly for the next periods. Based on this idea, one
classic approach, called “persistence forecasting”, is developed
[37], which uses past measurements of renewable power
outputs as forecasted values for the next periods. “Persistence
forecasting” is low cost. However, it only obtains an estimated
point, not a distribution, of future uncertainties. Therefore, the
classic “persistence forecasting” is not suitable for risk-limiting
constraints (4)(5). If one wants to use persistence forecasting
in the load restoration problem, the probabilistic constraints (4)
(5) need modifications. A detailed discussion is provided in
Appendix D.
In the following, a method that uses observations to infer
distributions of renewable generation for the next periods is
detailed.

PDF

6
5
3

4

0.60

(18)

det   m 

0.75

power of WT1
over different periods

6

4th Component

2.0

1.5
1.0
0.5
0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Value
Fig. 5 An illustrative example using a GMM to represent the PDF of a
random variable. In this example, K=1 and M=4. Weights for components are
0.4, 0.1, 0.3, and 0.2, respectively.

In Fig. 5, the idea of using a GMM to represent the PDF of a
random variable is presented. For illustration, suppose that the
number of periods K =1 and the number of Gaussian
components M = 4. The PDF of a GMM (in black) is a weighted
summation of four Gaussian distributions (in red, blue, green,
and purple, respectively). By adjusting the parameter set Ω,
GMM is able to characterize different kinds of non-Gaussian
correlated random variables [38]-[40]. Thus, it is suitable for
modeling uncertainties of WTs and PVs.
With historical data of X, the parameter set Ω of a GMM can
be determined through the maximum likelihood estimation
technique [41]. Specifically, this paper adopts a solver
gmdistribution.fit in MATLAB to estimate the parameter set Ω.
C. Using observations to update the prior distribution
Suppose in the second period, the actual power output in the
first period is observed. That is,
(19)
X1  x1
Then, the distribution of the remaining entries, X2,…,K, can be
updated as a conditional distribution with respect to X1=x1. In
Appendix A, the computational procedure of a conditional
distribution is provided:

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
f X 2,

,T

x

|X 1

2,

| x1 

,K

  l  x1 ,   N l  x2,

; μ l  x1 ,   ,  l  x1 ,   

M

l 1

,K

(20)

Note that the updated distribution of X2,…,K has a GMM form.
This is an important property of GMM, namely, conditional
invariance. When the third period comes with an observed
realization of X2, the distribution of X3,…,K can be updated in a
similar way as (20) via the formulae in Appendix A. A more
general updating formula for the period k+1 based on an
observation of Xk is given as follows:
f X k 1,

, K |X1,

,k

x



 Nl xk 1,

M

| x1,

,k

     x

;μl  x1,

,k

,   ,  l  x1,

k 1, , K

,K

l

l 1

1, , k

,k

, 

, 

(21)



[X1k+1 … XW+S
… X1K … XW+S
]T. Note that (24) is a GMM, and
k+1
K
1
T
random variables [Xk+1 … XW+S
can be regarded as a linear
k+1 ]
T
transformation of [X1k+1 … XW+S
…
X1K … XW+S
k+1
K ] :
1
 X k 1 




 X kW1+S 
 X k11 
(25)






  1 0 
 X1 
 X kW1+S 
E


 K 


 X W +S 
 K 
Using g Lemma 1 in Appendix B, one can obtain a
T
distribution of [X1k+1 … XW+S
k+1 ] :
f X k 1 |X1,

Realization
of Xk is
observed

Period k

f Xk ,

xk

Period k+1



Nl x

Distribution
of Xk, ,K

, K |X1, ,k 1

x

k, ,K

| x1,

, k 1



Equation (21) and Appendix A
Distribution of Xk+1,

f X k 1,

, K |X1, ,k

x

,K

k 1, , K

| x1,

,k



Fig. 6 A recursive procedure to update the distribution of power uncertainties

A flow chart is given in Fig. 6. The distribution of power
uncertainties for the next periods with respect to observations
in the past can be updated in a recursive manner.
D. Multiple renewable case
If there are multiple PVs and WTs, the random vector X is
augmented as follows:
T
X   X 11
X 1W  S
X k1
X kW  S
X K1
X KW  S  (22)




In period k, an aggregated observation of W+S components
is obtained, i.e.,
(23)
 X k1
X kW +S    x1k
xkW +S 








Then, the distribution of [X1k+1 … XW+S
… X1K … XW+S
]T with
k+1
K
respect to the aggregated observation can be updated via the
formulae in Appendix A. The updated distribution is denoted as
follows:
f X k 1,

,K

|X1,



,k

 Nl x

x

1, , S W
k 1, , K

1, , S W
k 1, , K

| x1,1,

;μ l  x

, S W
,k

1, , S W
1, , k

     x
M

l 1

1, , S W
1, , k

l

,   , l  x

1, , S W
1, , k

,   (24)

, 



The updated distribution (24) has a GMM form.
E. Integrating the updated distribution into the risk-limiting
constraints
When the kth observation is obtained, the conditional joint
distribution of [X1k+1 … XW+S
… X1K … XW+S
]T shown in (24) is
k+1
K
incorporated into the risk-limiting constraints (4)(5) as follows:
a) In (4), when t=k+1, the distribution of random variables
T
[X1k+1 … XW+S
k+1 ] should be extracted from the joint distribution of

7

,k

x

1, , S W
k 1

1, , S W
k 1

| x1,1,

, S W
,k

;Eμ l  x

1, , S W
1, , k

     x
M

l 1

1, , S W
1, , k

l

,   , E l  x

1, , S W
1, , k

,  

,  E

T

(26)



Obviously, this is a GMM.
b) In (4), when t>k+1, the distribution of random variables
1
T
[X 1t+1 … X W+S
… X W+S
]T.
t+1 ] should be conditioned on [X t
t
1
W+S T
However, in the kth period, [X t … X t ] have not been
observed. Similar to RLD, the idea of the proposed risk-limiting
load restoration strategy is to utilize the conditional expectation
of [X1t … XW+S
]T | [X1k … XW+S
]T to replace [X1t … XW+S
]T. The
t
k
t
1
W+S T
1
W+S T
distribution of [Xt … Xt ] | [Xk … Xk ] can be obtained in
a similar way as (25)(26), being a GMM. Thereafter, one can
use Lemma 2 in Appendix C to compute the conditional
expectation of [X1t … XW+S
]T | [X1k … XW+S
]T. Finally, using the
t
k
formulae in Appendix A, one can obtain the conditional
1
T
distribution of [X 1t+1 … X W+S
… X W+S
]T . Due to the
t+1 ] | [X t
t
conditional invariance of GMM, the obtained distribution is
still a GMM.
According to a) and b), the distributions of [X1t … XW+S
]T for
t
t=k+1 and t>k+1 are GMMs. For brevity, they are denoted as
follows:
f X 1,
t

,S W

1, ,S W
|X1,
,k



x

 N l ,t xt1,

1,
t

, S W

| x1,1,

, S W
,k

     x

;μl ,t  x1,1,

, S W
,k

,   ,  l ,t  x1,1,

, S W

M

l 1

1, , S W
1, , k

l ,t

, S W
,k

, 

, 

 (27)

t  k 1

c) Since the risk-limiting constraint (5) is associated with
random variables [X 1k+1 … X W+S
… X 1K … X W+S
]T, the joint
k+1
K
distribution in (24) is used for the risk-limiting constraint (5)
without modifications.
F. Comparison with a probabilistic forecasting tool
Usually, a probabilistic forecasting tool is designed to obtain
quantitative information of renewable generation random
variables conditioned on explanatory variables [38]. Take the
wind power as an example. The quantitative information could
be confidence intervals for any quantile and/or full PDFs of
wind generation. Typical explanatory variables include the
forecasted wind speed from a NWP system as well as forecasted
wind generation from a point forecasting tool. If one regards the
past observations [X11 … XW+S
… X1t … XW+S
]T as explanatory
1
t
variables, the proposed method becomes a probabilistic

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
forecasting tool for the future wind generation [X1t+1 … XW+S
…
t+1
T
X 1K … X W+S
]
.
Compared
with
the
traditional
probabilistic
K
forecasting tools [38], the proposed method has the following
features and advantages:
(1) The traditional methods rely on auxiliary forecasted
information of explanatory variables. In contrast, the proposed
method only relies on observations, which simplifies the
implementation and makes the method appropriate in cases
where the auxiliary forecasted information is not available.
(2) Usually, the traditional methods provide univariate PDFs
1
and/or quantiles of random variables X1t+1, …, XW+S
t+1 , …, XK , …,
W+S
XK , respectively. If so, it is not easy to handle the arithmetical
addition operations of random variables in the probabilistic
1
W+S
constraints (4)(5), as X 1t+1 , …, X W+S
are
t+1 , …, X K , …, X K
correlated random variables. The proposed method obtains a
1
W+S T
joint distribution of a random vector [X1t+1… XW+S
t+1 … XK … XK ] .
The joint distribution has a GMM form, which makes the
computation of the probabilistic constraints (4)(5) feasible.
G. Advantages
The proposed method modeling uncertainties has the
following advantages:
(1) Even if distributions of power outputs of WTs and PVs
are different from one another, the proposed method is able to
provide a universal model with satisfactory estimation by
adjusting the parameter set Ω. Also, the spatial-temporal
correlation can be taken into account.
(2) Based on latest collected observations, the proposed
method is able to infer distributions recursively as time rolls
forward, increasing the accuracy of descriptions of future
uncertainties.
(3) The updated distribution in each period is a GMM,
facilitating the solution of the original problem (2)-(14).This
point will be detailed in the next section.

L

u

t ,l

l 1

G

Pt ,l   Pt , g  CDFX1sum 1   

SOLUTION METHODOLOGY

where CDF-1(·) denotes the quantile of X t sum , which can be
computed through its distribution function (29).
B. Equivalent transformation of energy constraint (5)
Adopting a similar idea as (28) and (29), one can compute
the distribution of the aggregation of Xt v over multiple periods
from k+1 to T as follows:
X sum
1,T 
f X sum

 1,T

x

     x
l 1



T
1,
N l xsum
1,T ; e μ l  x1,

, S W
,

W +S

X
v 1

t   1

v
t

M

t



l 1

Nl xtsum ; eT μl ,t  x1,1,

, S W
,k

, S W
,k

,  

,   , eT  l ,t  x1,1,

,  

,   , eT  l  x1,1,

, S W
,k

,  e

(32)
, S W
,

,  e



L


u


Pt ,l 

1



EN g   1  CDFX1sum

 gESS
t  1 l 1
sum
-1
where CDF (·) of X κ+1,T
distribution function (32).
t ,l

 1,T

1    (33)

can be computed through its

C. Equivalent MILP
At this point, the original problem in Section III is converted
to an equivalent MILP as follows:
Objective :

Max (2)

Subject to:

(6)-(14), (30), and (33)

The MILP can be solved with available commercial solvers,
such as intlinprog in MATLAB [42].
D. Implementation procedure
The implementation procedure of the proposed scheme for
the kth period is shown in Fig. 7.
MILP

Optimal Results
ESS

Objective: max (2)
Realization of Xk
is observed
,K

Section
III-E

Subject to: (6)(14), (30) and (33)

, K |X1, ,k

Dispatchable
Generators
Load
Breakers

Fig. 7 Implementation of the proposed scheme

(28)

Based Lemma 1 in Appendix B, the distribution of Xt sum is
obtained as follows,
f X sum  xtsum    l,t  x1,1,

, S W
,

(31)

v
t

Thereafter, the equivalent transformation of (5) is:
T

f X k 1,

Equivalent transformation of power constraint (4)
Define the aggregation of Xt v in period t by
X tsum 

1,
1,

l

W S

t  1 v 1

M

sum

 1,T

Distribution of Xk+1,

A.

T


X


Renewable Generation

This section transforms probabilistic constraints (4)(5) into
equivalent deterministic linear inequailities.

 t    1 (30)

t

g 1

Modeling Uncertainty

IV.

8

 (29)

t  k 1

where e is a unit vector with a proper dimension.
Note that Lemma 1 holds even if Xt v is correlated with others.
Hence, the correlation of random variables can be handled.
Thereafter, moving the deterministic terms of (4) into left
side and random variables to the right side, one can convert the
probabilistic constraint (4) into (30) as follows:

First, measurements of [X1k … XW+S
]T are obtained, based on
k
which the joint distribution of future uncertainties is updated.
Then, the joint distribution is fed into the risk-limiting
constraints (4)(5), enabling the constraints to be converted into
linear inequalities. Thereafter, MILP is solved. Finally, the
optimal results, i.e., loads to be restored, the scheduled power
of dispatchable generators, and the charging/discharging
strategy of ESSs are sent to equipment for implementation.
E. Optimal control dynamic dispatch (OCDD) formulation
In the literature [43], the model (2)-(14) is categorized as
“dynamic economic dispatch (DED)”, because Pt,g are decision
variables. There is also another methodology for power system
operation called “optimal control dynamic dispatch”, in which
the ramp rate rt,g are decision variables [43].

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
These two models are subject to similar constraints and
implemented periodically. However, authors in [43] indicate
that there are two differences between the two models:
(1) The optimal decision of the OCDD depends on initial
values of P0,g, while the DED does not.
(2) The OCDD considers the ramp limit on P1,g -P0,g, while
the DED neglects such a constraint.
The load restoration problem in this paper can be formulated
as OCDD. The full OCDD formulation is provided in Appendix
E. The proposed method to update distributions of uncertainties
and the solution methodology to handle probabilistic
constraints (4)(5) can be readily used in the OCDD formulation.
V. CASE STUDY
A. Test system and data
In the following tests, there are three microgrids, each with a
diesel generator and an ESS. Besides, MG1 and MG2 both have
a WT; MG3 has a PV. The duration of an outage from 07:00 to
17:00 is 10 hour with 1-hour time resolution. As a result, the
total number of random variables is 30 (=3 × 10 ). Hourly
historical data of WTs/PVs comes from publicly available
datasets, i.e., “solar integration data set” and “wind integration
data set”, of NREL [36]. The data set consists of actual power
outputs of wind/solar generation for 3 years. This paper uses
data records of the first 1000 days as the “training set”, and the
rest as “test set”. The “training set” is used to estimate the
parameters of GMM using the maximum likelihood estimation
technique, while the “test set” is used to test the performance of
the proposed method. The proposed method is applicable to
other time resolution of data, e.g., 30 min, 10 min, if needed.
Load information of the IEEE 342-node system comes from
Pacific Northwest National Laboratory (PNNL) [44]. Important
information of the test system is listed in Table I.
TABLE I
IMPORTANT SYSTEM INFORMATION

MG1

Similar results can be obtained for WT2 and PV3.
To quantify the fitting performance, this paper uses three
probabilistic metrics.
(1) The likelihood function value is used [41]. This is a
positively-oriented score: the higher, the better. Note that the
likelihood function value is not an “absolute value metric”,
which means there is not a universal threshold to distinguish
between accuracy and inaccuracy. To overcome this hurdle, this
paper compares the likelihood function of GMM with three
widely-used statistical models [45]. They are the Gaussian
model, Gaussian Copula model, and t Copula model. By doing
so, the fitting performance of GMM is evaluated. The
likelihood function values of GMM and the other three models
are listed in Table II. It can be seen that GMM has the largest
value, remarkably outperforming the other three.
TABLE II
FITTING TEST OF MODELING [XT YT]T

: 2.0 MW, 5.0 MWh
: 0.5 MW, 2.0 MWh, SOC70%
: 2.0 MW
: 32 loads, 6.48MW in total

Diesel2
ESS2
WT2
Loads

: 3.0 MW, 8.0 MWh
: 1.5 MW, 3.0 MWh, SOC60%
: 2.0 MW
: 30 loads, 8.79 MW in total

Log-likelihood function values(103)

Gaussian

1.8345

Gaussian Copula

3.7658

t Copula

4.9025

GMM

7.1324

Historical data

GMM

0.12

Information
Diesel1
ESS1
WT1
Loads

Methods

Fig. 8 Joint distributions of X at 500 points

0.10
0.08
PDF

MG No.

9

0.06
0.04

MG2

0.02

0

0

5

10

15

20

Value of eTX

MG3

Diesel3
ESS3
PV3
Loads

: 2.5 MW, 10.0 MWh
: 1.0 MW, 4.0 MWh, SOC70%
: 2.0 MW
: 34 loads, 6.54 MW in total

B. Validation of modeling uncertainties
This subsection provides test results using GMM to model
uncertainties and updating distributions based on observations.
First, a GMM with 40 components is adopted to fit the
historical data of actual power outputs of WTs/PVs over 10
periods. The results of WT1 are detailed in the following.

Fig. 9 GMM fits historical data of eTX

(2) The joint distribution of X based on GMM is compared
with the empirical joint distribution of X. According to
Glivenko–Cantelli theorem [46], the empirical distribution
converges to the true underlying distribution with probability
one. Therefore, the empirical distribution can be used as a
benchmark. The values of GMM based joint distribution and
the empirical distribution at 500 points are shown in Fig. 8. The
maximum absolute error (MAE) [47] and root mean square
error (RMSE) [47] over the 500 points are computed. Both

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
MAE and RMSE are small (0.0324, 0.0089, respectively). The
results of Fig. 8 verify that GMM based joint distribution
coincides with the empirical distribution, i.e., GMM well
represents the historical data.

With
update

determined. Then, as the next period arrives, power outputs of
WTs/PVs are realized. Since power outputs of WTs/PVs are
random, an adjustment is activated: if the power supply is less
than the restored load demand, the dispatchable generators
should increase outputs and consume resources to meet the load
demand; otherwise, there should be wind and/or solar spillage.
For brevity, it is assumed that there is no way to regulate power
outputs of ESSs once they are determined. By doing so, it is
straightforward and convenient to use the restored loads,
generator regulations, and wind/solar spillage to evaluate the
effectiveness of the load restoration strategy.
Distribution with update

Resilience value

80

Without
update

40
20
0

C. Results of the load restoration plan
1) Individual microgrid case
With the validated model of uncertainties, this subsection is
concerned with the load restoration using the updated
distributions. In this test, each microgrid operates in a standalone mode with no connection to others, instructs its own
energy resources, and supplies its own loads. In order to show
the advantage of the recursively updated distribution, the
distribution without update is provided as a comparison.
The load restoration strategy is tested with the “test set” in an
operation simulation. First, a generation/load instruction for the
next arriving time is obtained. At this point, restored loads are

1

1.5

Power (MW)

Fig. 10 distributions of WT1 in MG1 over 10 hours with/without update

Distribution without update

60

Confidence
level

(3) Note that when X is represented by a GMM, the PDF of
an aggregated random variables eTX can be obtained based on
Lemma 1. Obviously, if a GMM well represents X, the PDF of
eTX, which is derived from GMM, should match well with the
historical data of eTX. Based on this principle, this paper
compares the PDF of eTX derived from GMM with the
histogram of eTX. The results are shown in Fig. 9. It can be seen
that GMM matches well with the histogram. MAE and RMSE
are also small (0.0112, and 0.0045, respectively). The results
further indicate that GMM is a good representation for X.
Secondly, this paper computes the recursively updated
distributions. In Fig. 10, WT1 in MG1 is taken as an example.
The width of confidence levels of the updated distribution
becomes narrower as time rolls forward. Meanwhile, the actual
power is well bounded within confidence intervals. As a
comparison, with the same confidence level, the width of the
distribution without update is relatively broad, indicating a
rough description of wind power uncertainties. The results
reveal that the updated distribution provides a more elaborate
description of wind power uncertainties than the one without
updating. Similar results are obtained for WT2 in MG2, and
PV3 in MG3.

10

2

3

4

5
6
Time

7

8

9

10

Wind spillage with update

Generator regulation with update

Wind spillage without update

Generator regulation without update

1.0
0.5
0

1

2

3

4

5
6
Time

7

8

9

10

Fig. 11 Resilience values, generator regulations, and wind spillage of MG1

Fig. 11 shows resilience values, wind spillage, and generator
regulations of MG1 (MG2 and MG3 are omitted). Three results
can be obtained from Fig. 11:
(1) At the beginning, the resilience values of distributions
with/without update are the same, since there has not been any
observation to update the distribution.
TABLE III
COMPARISON OF DISTRIBUTION WITH/WITHOUT UPDATE

Case

Resilience values

Distribution with update

513.3

Distribution without update

342.6

(2) As time rolls forward, the recursively updated distribution
takes advantages with more restored loads and less wind
spillage over the distribution without update. This can be
explained as follows: based on the distribution without update
which provides a rough description of wind uncertainties,
operators have to make conservative decisions to ensure that the
probabilistic power/energy adequacy constraints (4)(5) are
satisfied. That is, one uses more energy from the diesel
generator of limited generation resources to serve loads,
resulting in more wind spillage. As a comparison, the updated
distribution provides an accurate estimation for the future
uncertainties. Therefore, even with more wind power utilized,

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
operators are confident that the adequacy requirements (4)(5)
can be satisfied. The resilience values over the outage duration
are listed in Table III. The results also verify that the updated
distribution outperforms the distribution without update.
(3) An interesting phenomenon in Fig. 11 is that the wind
spillage appears more frequently than the generator regulation.
Here is an explanation: the confidence levels α in probabilistic
constraints (4)(5) are 90%. That is, once the restored loads are
determined, there is only a 10% chance that the load demand
cannot be met, i.e., the generator regulation is thereafter
activated. Therefore, in this test, the generator regulation is
infrequently observed.
2) Networked microgrid case
In this test, the three microgrids connect with each other,
become networked microgrids, and share resources to restore
loads. This subsection discusses the advantage of networked
microgrids over stand-alone ones. The resilience values of
networked microgrids are computed. In Fig. 12, the networked
resilience values are greater than the summation of three
resilience values in stand-alone microgrids. This is reasonable:
in networked microgrids, generation resources from one
microgrid can be transferred to another microgrid to restore
critical loads if needed.

VI. CONCLUSION

100

To enhance the resilience of a distribution system with
microgrids, this paper proposes a method using observations to
infer distributions of uncertainties of WTs/PVs based on the
conditional invariance property of a GMM. The transformation
of the original risk-limiting load restoration problem into an
MILP benefits from the recursively updated distributions being
a GMM.
This paper is focused on uncertainties. The power/energy
adequacy requirements are considered, laying a foundation for
load restoration with intermittent energy resources for practical
use. Furthermore, concerning the operational and dynamic
constraints, there are two on-going studies:
(1) Transmission limits should be considered. If the power
flow is linear, e.g., the linearized DistFlow, a method proposed
in previous authors’ work [48] is readily applicable to deal with
linear chance-constrained line power limits without major
change. For the AC power flow in which the line power is a
nonlinear implicit function of random power injections, the
tractable computation of the risks pertaining to line limits
remains an issue to resolve. The sample average approximation
(SAA) is a potential solution.
(2) Dynamic constraints should be considered. Usually,
distributed generators and ESSs have limited capabilities to
withstand large transient shocks. Therefore, it is necessary to
run dynamic simulations to ensure that the load restoration
strategy does not cause instability or damage to equipment.
Further discussions that are useful for the future work can be
found in [6], [7].

50

APPENDIX A

Summation of stand-alone microgrids

Networked microgrids

250

Resilience value

11

200
150

0

1

2

3

4

5

6

7

8

9

10

Time

CONDITIONAL INVARIANCE OF A GMM
Let the random vector X be decoupled into two parts:

Fig. 12 Comparison of restored loads between networked microgrids and standalone microgrids

D. Computation time
As far as the computation time is concerned, all tests are
implemented on a Core-i5 PC, with a 2.39-GHz processor and
8 GB of RAM. The coding environment is Matlab. In the rolling
plan, the optimization problem at the first period (k=0) has the
heaviest computational burden since it determines decisions
with the longest optimization window (see Fig. 3). It takes
11.649s to solve the problem for the individual case and
34.889s for the networked case. Note that the time resolution of
the problem is set to be 1 hour in this paper. Therefore, the
proposed method meets the time requirement. Even if the time
resolution is reduced to 10 min, the proposed method is still
applicable.
As far as communication latency is concerned, this paper
assumes that the communication structure follows the standard
IEC 61850 via TCP/IP [22]. Since a microgrid is not a large
system, the communication delay is not significant. Suppose
that the communication latency is 1 second. The computation
and communication time together is within a minute. Therefore,
the communication latency does not cause a major hurdle.

X :  Y T

Z T 

T

(34)

where Y represent observations, and Z represent the rest entries
of X.
The distribution of X is a GMM. That is:

f X  x   fYZ  y, z    m N m  y, z ; μ m ,  m 

(35)

μ y 
 yy  myz 
μ m :  mz  ,  m   mzy
zz 
μ m 
m m 

(36)

M

m 1

Then, the conditional distribution of Z with respect to Y=y is
given as follows [49]:
f Z |Y  z | y    l  y,  l N l  z ; μ l  y,   ,  l  y,    (37)
M

l 1

l  y,    l

Nl  y ; μly ,  ly 
M


m 1

μl

m

N m  y ; μ my ,  my 

 y,    μlz   lzy   lyy 
l

(38)

1

 y,     lzz   lzy   lyy

 y -μ 
 
y
l

1

yz
l

(39)
(40)

The distribution of Z with respect to Y=y is a GMM, enabling
a recursive way to update distributions.

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
APPENDIX B

Distribution with update

12

Persistence forecasting

80

Lemma 1 [49] : If the distribution of X is a GMM as in (16),
and XLT is a defined as a linear transformation of X:
(41)
X LT = AX + C
Then, the distribution of XLT is computed as follows:

Resilience value

LINEAR INVARIANCE OF A GMM

f X LT  xLT    m N m  xLT ; Aμ m  C, A m A T  (42)

60
40
20

M

0
1

m 1

2

3

4

5

Lemma 1 holds even if entries of X are correlated.
Wind spillage with update
Wind spillage of the
persistence forecasting

APPENDIX C
Lemma 2: If the distribution of X is a GMM as in (16), the
expectation of X is
M

E  X    m μ m

(43)

m 1

Power (MW)

0.8

EXPECTATION OF A GMM

8

9

10

Generator regulation with update
Generator regulation of the
persistence forecasting

0.4
0.2
0
2

3

4

5

6

7

8

9

10

APPENDIX D

Time

COMPARISON WITH PERSISTENCE FORECASTING

Fig. D-1 Resilience value, generator regulations, and wind spillage of MG1

A. Modification of probabilistic constraints (4)(5)
Since “persistence forecasting” is a point forecasting tool of
uncertainties, the probabilistic constraints (4)(5) should be
modified as follows:
G
W +S
L
P 
X v  fore  
u P
 t    1 (44)


g 1

gESS

7

0.6

1



6

Time



t,g

v 1



t

EN g   1  

l 1

T

W S

t ,l

t ,l

T

t  1 v 1

t  1 l 1

  t   t ,l t ,l
 v 1
 l 1
T
W S
T
L


EN g   1  E    X tv     Pt ,l

gESS
 t  1 v 1
 t  1 l 1
g 1

Resilience value

Proposed updated expectation

568.2

Distribution without update

408.3

APPENDIX E

t ,l

where X vr (fore) represents the forecasted value from
“persistence forecasting” technique.
If one wants to compare the proposed method with
“persistence forecasting” technique, the updated distribution
should be modified as a point prediction. An intuitive idea is
using the expectation of the updated distribution as a point
prediction. Thereafter, probabilistic constraints (4)(5) should
be formulated as expectations. That is:
G
L
W +S v 
P E
X

u P
 t    1 (46)



Case

L


 X  fore   
 P  (45)


v
t

TABLE D-I
RESILIENCE VALUE COMPARISON

t,g

(47)

where E[·] is the expectation operation.
The distributions of the aggregation of Xt v have been derived
in (29)(32), which are GMMs. Therefore, the expectations can
be computed based on Lemma 2 in Appendix C.
At this point, one can compare “persistence forecasting” with
the proposed updated expectation.
B. Test results
The test procedure is the same as those in Section V-C. The
results are shown in Fig. D-1 and Table D-I. It can be seen that
the proposed method with updated expectations obtains a
higher resilience value, surpassing the “persistence forecasting”
method.

AN EXTENSION TO THE OCDD MODEL
C. Full formulation
● Risk limiting constraints:
t
L
G 

 W +S
Pr   P0, g   ri , g    X tv   ut ,l Pt ,l    (48)
i 1
l 1
 v 1
 g 1 

T W S

v 
  EN g   1    X t  
(49)
 gESS

t  1 v 1
Pr 

T
L


   ut ,l Pt ,l


t  1 l 1



● Deterministic constraints
t

Pgmin  P0, g   ri , g  Pgmax
i 1

(50)

t    1,  g  ESS
T




P
 

t  1


  ri , g   EN g   1
i 1

 g  ESS
t

0, g

t , g   t , g  1

t    1, g  ESS

(51)
(52)

t

dch
dch,max
0  P0,dch
g   ri , g   t , g Pg
i 1

(53)

t    1, g  ESS
t

 t , g Pgch,max  P0,chg   ri ch
,g  0
i 1

t    1, g  ESS

(54)

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
[3]

t

SOCgmin  SOC0, g   rsoci , g  SOCgmax
i 1

(55)

t    1, g  ESS



dch 1
ch
ch
rsoct , g    P0,dch
g   ri , g  d +P0, g   ri , g  c  / ECg (56)
i 1
i 1


t    2, g  ESS
t 1

[4]

t 1

[5]

[6]

D. Test results
This paper takes the optimization problem of the first period
(k=0) as an example to show the difference between DED and
OCDD models. As discussed in Section IV-E, the OCDD
considers an initial value of P0,g and the ramp limit on P1,g -P0,g,
while the DED does not. To clearly illustrate this phenomenon,
P0,g in OCDD is set to be 0 . That is to say, the maximum value
of P1,g equals to the ramping rate of generator g. In this paper,
the ramping rate is set to be 0.3MW/h.

[7]

[8]

[9]

[10]

[11]

[12]
[13]

[14]
Fig. E-1 Scheduled power of the diesel generator in MG1

[15]

[16]

[17]

[18]

[19]
Fig. E-2 Objective functions of the two models
[20]

Tests are conducted on MG1. Results of the two models are
shown in Fig. E-1 and Fig. E-2. In Fig. E-1, it is observed that
the scheduled power of the diesel in OCDD model and that in
DED model are different at the beginning, but they coincide
later. Similar conclusions are obtained for the objective
functions in Fig. E-2. Because the OCDD model has more
constraints than DED, the optimal solution of OCDD is no
better than that of DED. According to Fig. E-2, the objective
functions of the two models are close to each other (342.57,
341.58, respectively). The results of the two figures are
consistent with those reported in [43].

[21]

[22]

[23]

[24]

[25]

REFERENCES
[1]

[2]

L. Che, M. Khodayar, and M. Shahidehpour, “Only connect: Microgrids
for distribution system restoration,” IEEE Power Energy Mag., vol. 12,
no. 1, pp. 70–81, Jan./Feb. 2014.
Office of the press secretary of the white house, presidential policy
directive – critical infrastructure security and resilience [Online].
Available:http://www.whitehouse.gov/the-press-office/2013/02/12/presi
dential-policy-directive-critical-infrastructure-security-and-resil

[26]

[27]

13

U.S. Department of Energy and National Energy Technology Laboratory,
Operates resiliently against attack and natural disaster [Online]. Available:
http://www.smartgridinformation.info/pdf/1438_doc_1.pdf
V. Madani, et al., “Distribution automation strategies challenges and
opportunities in a changing landscape,” IEEE Trans. Smart Grid, vol. 6,
no. 4, pp. 2157-2165, July 2015.
T. T. H. Pham, Y. Besanger, and N. Hadjsaid, “New challenges in power
system restoration with large scale of dispersed generation insertion,”
IEEE Trans. Power Syst., vol. 24, no. 1, pp. 398-406, Feb. 2009.
Y. J. 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, early access article.
Y. Xu, C. C. Liu, K. Schneider, F. Tuffner, and D. Ton, “Microgrids for
service restoration to critical load in a resilient distribution system,” IEEE
Trans. Smart Grid, early access article.
Z. Wang, B. Chen, J. Wang, and C. Chen, “Networked microgrids for selfhealing power systems,” IEEE Trans. Smart Grid, vol. 7, no. 1, pp. 310319, Jan. 2016.
H. Farzin, M. F. Firuzabad, and M. M. Aghtaie, “Enhancing power system
resilience through hierarchical outage management in multi-microgrids,”
IEEE Trans. Smart Grid, vol. 7, no. 6, pp. 2869-2879, Nov. 2016.
Z. Wang and J. Wang, “Self-healing resilient distribution systems based
on sectionalization into microgrids,” IEEE Trans. Power Syst., vol. 30, no.
6, pp. 3139-3149, Nov. 2015.
A. Gholami, T. Shekari, F. Aminifar, and M. Shahidehpour, “Microgrid
scheduling with uncertainty: the quest for resilience,” IEEE Trans. Smart
Grid, vol. 7, no. 6, pp. 2849-2858, Nov. 2016.
A. Khodaei, “Resiliency-oriented microgrid optimal scheduling,” IEEE
Trans. Smart Grid, vol. 5, no. 4, pp. 1584–1591, Jul. 2014.
W. Yuan, J. Wang, F. Qiu, C. Chen, C. Kang, and B. Zeng, “Robust
optimization-based resilient distribution network planning against natural
disasters,” IEEE Trans. Smart Grid, vol. 7, no. 6, pp. 2817-2826, Nov.
2016.
P. P. Varaiya, F. F. Wu, and J. W. Bialek, “Smart operation of smart grid:
risk-limiting dispatch,” Proc. IEEE, vol. 99, no. 1, pp. 40-57, Jan. 2011.
R. Rajagopal, E. Bitar, P. Varaiya, and F. Wu, “Risk-limiting dispatch for
integrating renewable power,” Int. J. Elec. Power, vol. 44, no. 1, pp. 615628, Jan. 2013.
B. Zhang, R. Rajagopal, and D. Tse, “Network risk-limiting dispatch:
optimal control and price of uncertainty,” IEEE Trans. Autom. Control,
vol. 59, no .9, pp. 2442-2456, Sep. 2014.
C. Wu, G. Hug, and S. Kar, “Risk-limiting economic dispatch for
electricity markets with flexible ramping products,” IEEE Trans. Power
Syst., vol. 31, no. 3, pp. 1990-2003, May 2016.
C. Peng, Y. Hou, N. Yu, and W. Wang, “Risk-limiting unit commitment
in smart grid with intelligent periphery,” IEEE Trans. Power Syst., early
access article.
C. Peng, Y. Hou, N. Yu, J. Yan, S. Lei, and W. Wang, “Multi-period risklimiting dispatch in power systems with renewables integration,” IEEE
Trans. Industrial Informatics, vol. 13, no. 4, pp. 1843-1854, Aug. 2017.
S. M. Ross, Introduction to Probability Models 11th Edition. San Diego,
CA, USA: Elsevier, 2015.
A. R. Malekpour and A. Pahwa, “Stochastic networked microgrid energy
management with correlated wind generators,” IEEE Trans. Power Syst.,
vol. 32, no. 5, pp. 3681-3693, Sep. 2017.
IEEE Standard for Exchanging Information Between Networks
Implementing IEC 61850 and IEEE Std 1815(TM), IEEE Std 1815.12015, Dec. 2016
C. Wang, et al., “A highly integrated and reconfigurable microgrid testbed
with hybrid distributed energy sources,” IEEE Trans. Smart Grid, vol. 7,
no. 1, pp. 451-459, Jan. 2016.
A. Ruiz-Alvarez, A. Colet-Subirachs, F. A. Figuerola, O. GomisBellmunt, and A. Sudria-Andreu, “Operation of a utility connected
microgrid using an IEC 61850-based multi-level management system,”
IEEE Trans. Smart Grid, vol. 3, no. 2, pp. 858-865, June 2012.
H. Gao, Y. Chen, S. Mei, S. Huang, and Y. Xu, “Resilience-oriented prehurricane resource allocation in distribution systems considering electric
buses,” Proc. IEEE, vol. 105, no. 7, pp. 1214-1233, July 2017.
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, March 2016.
M. Panteli and P. Mancarella, “The grid: stronger, bigger, smarter?
presenting a conceptual framework of power system resilience,” IEEE
Power Energy Mag., vol. 13, no. 3, pp. 58–66, May/Jun. 2015.

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT)
[28] H. Gao, Y. Chen, Y. Xu, and C. C. Liu, “Resilience-oriented critical load
restoration using microgrids in distribution systems,” IEEE Trans. Smart
Grid, vol. 7, no. 6, pp. 2837-2848, Nov. 2016.
[29] Z. Bie, Y. Lin, G. Li, and F. Li, “Battling the extreme: a study on the
power system resilience,” Proc. IEEE, vol. 105, no. 7, pp. 1253-1266,
July 2017.
[30] Z. Li, M. Shahidehpour, F. Aminifar, A. Alabdulwahab, and Y. Al-Turki,
“Networked microgrids for enhancing the power system resilience,” Proc.
IEEE, vol. 105, no. 7, pp. 1289-1310, July 2017.
[31] C. C. Liu, “Distribution systems: reliable but not resilient? [In My View],”
IEEE Power and Energy Magazine, vol. 13, no. 3, pp. 93-96, May-June
2015.
[32] A. Gholami and F. Aminifar, “A hierarchical response-based approach to
the load restoration problem,” IEEE Trans. Smart Grid, vol. 8, no. 4, pp.
1700-1709, July 2017.
[33] J. V. Seguro and T. W. Lambert, “Modern estimation of the parameters of
the Weibull wind speed distribution for wind energy analysis,” Journal of
Wind Engineering and Industrial Aerodynamics, vol. 85, no. 1, pp, 75-88,
2000.
[34] M. Fan, V. Vittal, G. T. Heydt, and R. Ayyanar, “Probabilistic power flow
studies for transmission systems with photovoltaic generation using
cumulants,” IEEE Trans. Power Syst., vol. 27, no. 4, pp. 2251-2261, Nov.
2012.
[35] F. J. Ruiz-Rodriguez, J. C. Hernandez, and F. Jurado, “Probabilistic load
flow for radial distribution networks with photovoltaic generators,” IET
Renewable Power Generation, vol. 6, no. 2, pp. 110-121, March 2012.
[36] C. Draxl, B. M. Hodge, A. Clifton, and J. McCaa. “Overview and
meteorological validation of the wind integration national dataset toolkit,”
Technical Report, NREL/TP-5000-61740. National Renewable Energy
Laboratory. Golden, CO. 2015.
[37] M. Milligan, M. Schwartz, and Y. Wan. “Statistical wind power
forecasting models: results for us wind farms,” Preprint. NREL/CP-50033956. National Renewable Energy Laboratory, Golden, CO., 2003.
[38] M. Nijhuis, M. Gibescu, and J. F. G. Cobben, “Gaussian mixture based
probabilistic load flow for lv-network planning,” IEEE Trans. Power Syst.,
early access article.
[39] R. Singh, B. Pal, and R. Jabr. “Statistical representation of distribution
system loads using gaussian mixture model,” IEEE Trans. Power Syst.,
vol. 25, no. 1, pp. 29-37, Feb. 2010.
[40] G. Valverde and A. T. Saric, “Probabilistic load flow with non-Gaussian
correlated random variables using Gaussian mixture models,” IET Gen.,
Trans. & Distri., vol. 6, no. 7, pp. 701–709, 2012.
[41] A. Hald, “On the history of maximum likelihood in relation to inverse
probability and least squares,” Statistical Science, vol. 14, no. 2, pp. 21422, 1999.
[42] Mathworks Support Documentation on Mixed-integer Linear
Programming (MILP). [Online]. Available: https://www.mathworks.com/
help/ optim/ug/intlinprog.html
[43] X. Xia, J. Zhang, and A. Elaiw, “An application of model predictive
control to the dynamic economic dispatch of power generation,” Control
Engineering Practice, vol. 9, no. 6, pp. 638-648, 2011.
[44] K. P. Schneider, P. Phanivong, and J. S. Lacroix, “IEEE 342-node low
voltage networked test system,” in Proc. IEEE PES General Meeting,
National Harbor, MD, 2014.
[45] G. Papaefthymiou and D. Kurowicka, “Using copulas for modeling
stochastic dependence in power system uncertainty analysis,” IEEE Trans.
Power Syst., vol. 24, no. 2, pp. 40-49, Feb. 2009.
[46] A. W. Van der Vaart, Asymptotic Statistics. Cambridge, England:
Cambridge University Press, 1998
[47] J. Zhang, et al., “A suite of metrics for assessing the performance of solar
power forecasting,” Sol. Energy, vol. 111, pp.701-709, 2012
[48] Z. W. Wang, C. Shen, F. Liu, X. Y. Wu, C-C Liu, and F. Gao. “Chanceconstrained economic dispatch with non-Gaussian correlated wind power
uncertainty,” IEEE Trans. Power Syst. early access article.
[49] J. T. Flam, “The linear model under Gaussian mixture inputs,” Ph.D.
dissertation, Dept. of Elect. & Telecom. Norwegian Univ. of Science and
Technology, Trondheim, Norway, 2013.

14

