This article has been accepted for publication in a future issue of this journal, but has not been

fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TSG.2018.2848958, IEEE
Transactions on Smart Grid
1

Maximum Loadability of Islanded Microgrids

with Renewable Energy Generation
Xiaoyuan Xu, Member, IEEE, Zheng Yan, Mohammad Shahidehpour, Fellow, IEEE,
Sijie Chen, Member, IEEE, Han Wang, Zhiyi Li, Member, IEEE, and Quan Zhou, Student Member, IEEE

PGi,lb, QGi,lb Minimum active and reactive power of unit i

Abstract--A microgrid (MG) is a small-scale power system PLi, QLi Active and reactive power loads of bus i
which is fed by constrained distributed generation (DG) units PLi0, QLi0 Nominal values of loads
and its continuous operation is affected by the variability of αpi, βqi Exponents of loads
available generation resources. In this paper a global sensitivity
hpi, hqi Frequency sensitivity parameters
analysis (GSA) method is proposed to evaluate the impact of
variable energy resources on the maximum loadability of PLi,b, QLi,b Base active and reactive power loads
islanded MGs (IMGs). First, a probabilistic optimization bPi, bQi Active and reactive power load changes
problem is formulated to calculate the IMG load margin k, c Parameters of a Weibull distribution
considering the droop characteristics of DG units and rmax Maximum solar irradiance
uncertainties of renewable generation, loads and distribution a, b Parameters of a Beta distribution
feeder parameters. Then, the global sensitivity (GS) is introduced
Yij Element of a bus admittance matrix
that can identify the impact of independent and correlated
variables on the IMG loadability. Next, the sparse polynomial δij Phase angle of Yij
chaos expansion (SPCE) method is used to obtain the ε Small positive number
probabilistic models for IMG load margins, and an efficient GSA pi Probability density function (PDF) of xi
method is proposed to calculate the GS of IMG loadability to Fi Cumulative distribution function (CDF) of xi
prevailing variables. The probabilistic models are considered for I, Φ PDF and CDF of a normal variable
IMG input variables and the impact of variable correlations on
IMG loadability is analyzed. The proposed method for In, Φn PDF and CDF of a n-dimensional normal vector
calculating the maximum loadability is tested using a 33-bus μ Mean vector
IMG and the results are compared with those of other GSA and Σ Covariance matrix
local sensitivity analysis (LSA) methods. C Copula function of correlated variables
Index Terms—Microgrids, maximum loadability, global n Total number of input random variables
sensitivity analysis, polynomial chaos expansion s Number of variables in the subset y
N Sample size of Monte Carlo simulation
NOMENCLATURE Nc Number of collocation points
Sets fl1…lk Item of ANOVA decomposition of function f
ΩB Set of buses D, f0 Variance and mean of f
ΩGd Set of droop-controlled DG units Dl1…lk Variance of fl1,…,lk
ΩGpv Set of DG units at PV buses Sy Global sensitivity of y
Ωs Set of input random variables ai Coefficient of polynomial chaos expansion (PCE)
Hi ith order Hermite polynomial
Parameters m Order of PCE
ω0 Angular frequency at no-load
Vi0 Voltage magnitude of bus i at no-load Variables
mpi, nqi Active and reactive power droop gains ω Angular frequency
PGi,ub, QGi,ub Maximum active and reactive power of unit i λ Load margin of IMG
Vi, θi Voltage magnitude and angle of bus i
PGi, QGi Active and reactive power of a DG unit
This work was supported by the National Key Research and Development PRi, QRi Active/reactive power of a renewable DG unit
Program of China (No. 2017YFB0902800), National Natural Science
Foundation of China (No. 51707115) and China Postdoctoral Science Pij Active power flow from bus i to bus j
Foundation (No. 2017M610251). X. Xu, Z. Yan, S. Chen, and H. Wang are Vi,1, Vi,2, Vi,3 Items of voltage magnitude of bus i
with the School of Electronic Information and Electrical Engineering, x Input random variable vector
Shanghai Jiao Tong University, Shanghai, 200240, China (e-mail:
[email protected]; [email protected]; [email protected]; xi ith input random variable
[email protected]). M. Shahidehpour, Z. Li, and Q. Zhou are with y, z Subset of x and its complementary subset
the Galvin Center for Electricity Innovation, Illinois Institute of Technology, z Variable vector of conditional probability
Chicago, IL 60616 USA (e-mail: [email protected]; [email protected];
[email protected]). M. Shahidehpour is also a Research Professor with distribution
the Renewable Energy Research Group, King Abdulaziz University, Jeddah, γ, α, β Normal variable vectors
Saudi Arabia. ˆ
[i , [i Normal variables

1949-3053 (c) 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TSG.2018.2848958, IEEE
Transactions on Smart Grid
2

I. INTRODUCTION problems. In [17]-[18], an improved GSA method was

M ICROGRID (MG) is operated either in grid-connected proposed to evaluate the impact of input variables on DS
or islanded mode. In the grid-connected mode, the power flows. This method was improved by using sparse
frequency and voltages of MGs are determined by the utility polynomial chaos expansion (PCE) to perform the high-
grid which is represented as an infinite bus. In islanded MGs dimensional model representation (HDMR) decomposition,
(IMGs), the frequency and voltages are regulated by DG units thus presenting analytical GS solutions for independent
governed by the respective droop characteristics. The IMG variables. However, the HDMR decomposition is non-unique
operation could enhance the reliability and the resilience of for correlated variables [19]. Besides, the correlation of
utility grid when the power system is facing outages. In the variables, which affects the solution for uncertain variables,
literature, various methods have been proposed for the was not detailed in previous GSA applications.
planning and operation of IMGs [1]-[8]. In this paper, GSA is applied to IMGs with renewable
The maximum loadability of transmisison systems, which energy generation. First, a probabilistic nonlinear optimization
maintains the large-scale voltage stability, was studied problem is established to calculate the load margin of IMGs
extensively. For a distribution system (DS), the related work is considering DG droop characteristics and variabilities of
quite limited since the DS loadability was traditionally renewable energy generation, loads and DS network
managed by the connected main grid. Recently, the maximum parameters. Then, GSA is applied to IMG loadability with
IMG loadability has attracted additional attention since IMGs correlated input variables using conditional probability and
are fed in the island mode by small DG units which could Copula functions. Next, a surrogate model of load margin is
enhance the loadability of IMG. The existing loadability established to calculate GS values, and an efficient GSA
approach evaluates the distance from the current operating method is proposed to analyze the impact of input variables on
point to the critical operating point for maintaining reliable the IMG load margin. The proposed GSA method is tested
and stable operations of IMGs. Ref. [9] studied the maximum using a 33-bus IMG with various types of input variables, and
IMG loadability by calculating the equilibrium point using the the results are compared with those considering other
Fischer-Burmeister-based method. Ref. [10] evaluated the sensitivity analysis methods to validate the effectiveness of the
maximum loadability of droop-controlled IMGs using the proposed GSA method.
optimal power flow method in which a complementary
constraint was added to describe the capacity-limited DG unit II. MAXIMUM LOADABILITY ASSESSMENT OF IMG
operations. Ref. [11] proposed a probabilistic algorithm to In this section, the models are presented for DG units and
determine the optimal droop settings of IMGs. The algorithm loads. Then, probability models are developed for input
adopted a hierarchical approach to expand the IMG load variables. Finally, a probabilistic optimization problem is
margin while satisfying the prevailing operating constraints. established to evaluate the maximum IMG loadability.
IMGs can face various uncertainties including load and
A. DG Model
renewable energy forecast errors as well as random DS
network outages. Accordingly, probabilistic methods are There are three types of power system nodes, including PV,
applied when evaluating the maximum IMG loadability in PQ and reference nodes, for balancing the total generation and
which the scenarios with minor impacts on IMG operations load. In IMGs, DG units are responsible for regulating the
are disregarded. Therefore, a quick identification of critical frequency and voltages while satisfying load sharing. DG units
variables could accelerate IMG loadability analyses and in IMGs are usually modeled as droop, PV or PQ nodes.
improve its controllability. In general, there are three droop control strategies: P-f/Q-V,
Global sensitivity analysis (GSA) is an effective method for P-V/Q-f and generalized droop, with the assumption of
quantifying the impact of IMG input variables on the MG inductive, resistive and complex impedances, respectively. In
performance [12]. Different from local sensitivity analysis this paper, the P-f/Q-V droop is adopted, which decouples
(LSA), GSA considers system nonlinearities and interactions active and reactive power controls and mimics the behavior of
among pertinent variables. By decomposing the output a synchronous generator in power systems. Moreover, a
variance into segments attributed to a subset of inputs, we virtual impedance control is applied to satisfy the inductive
define the ratio of the segment to the total variance as the system requirement.
global sensitivity (GS) of input variables for evaluating IMG The proposed GSA method is applied to a general
variables. stochastic problem, which is not affected by the DG unit
GSA was generally applied to artificial neural network [13] model. For the P-f/Q-V droop, active and reactive power of
and hydrological modeling [14], with a limited applications to DG units are determined as
power systems. Recently, [15] introduced GSA to evaluate ­° PGi (Z0  Z ) mpi
uncertainties in DS operations. Ref. [16] compared GSA with ® , i  :Gd (1)
°̄QGi (Vi 0  Vi ) nqi
LSA for identifying critical uncertainties in small signal
stability analyses. However, GS was mostly calculated using a For droop-controlled DG units, the output power follows
crude Monte Carlo simulation (MCS) method, which required the droop characteristic before it reaches the DG limits.
large sample sizes to obtain accurate results, making GSA a However, beyond the active power generation limit PGi,ub, the
time-consuming and expensive approach for large-scale DG unit is transformed into constant active power PGi,ub.

1949-3053 (c) 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TSG.2018.2848958, IEEE
Transactions on Smart Grid
3

Similarly, beyond the reactive power generation limit QGi,ub, A Beta model is used to the describe the variability of solar
the DG unit is transformed into constant reactive power QGi,ub irradiance for PV panels which is stated as
[10]. The droop-controlled capacity-limited DG units are a 1 b 1
described by the following constraints: * ab § r · § r ·
p r ¨ ¸ ¨1  ¸ (7)
* a * b © rmax ¹ © rmax ¹
PGi ,ub  PGi t 0, (Z0  Z ) mpi  PGi t 0, i  :Gd (2.1)
2) Load: The variability of active power loads is described by
PGi ,ub  PGi (Z0  Z ) mpi  PGi 0, i  :Gd (2.2) a Normal distribution. The reactive power load will maintain a
QGi ,ub  QGi t 0, (Vi 0  Vi ) nqi  QGi t 0, i  :Gd (2.3) constant power factor.
3) Parameter: System parameters corresponding to generators,
QGi ,ub  QGi (Vi 0  Vi ) nqi  QGi 0, i  :Gd (2.4) network, and loads are subject to forecast and measurement
For DG units operating in a PV mode, we define Vi,1 as the errors and system operating condition uncertainties. In this
desired voltage magnitude and use two non-negative variables paper, parameter uncertainties are represented by a uniform
Vi,2 and Vi,3 to describe the switching between PV and PQ distribution over a certain range [21].
modes, as stated in (3.1)-(3.4) [20]. If the reactive power 4) Correlation: Copula functions are adopted to represent
output is between QGi,lb and QGi,ub, (3.3) and (3.4) force Vi,2=0 correlated variables [22].
and Vi,3=0, with a voltage magnitude that is equal to the D. Maximum loadability assessment model
desired Vi,1. When the reactive power output reaches its limit,
the PV mode is switched to the PQ mode with constant The load margin is an effective index to evaluate the
reactive power. Specifically, if the reactive power output maximum loadability of power systems. Considering the
reaches its maximum, we have QGi=QGi,ub and QGi>QGi,lb. models of DG units, loads and their respective uncertainties,
Therefore, Vi,2=0 according to (3.3) and Vi=Vi,1-Vi,3. Similarly, the following optimization problem is considered to calculate
if the reactive power output reaches its minimum, we have the IMG load margin:
QGi=QGi,lb and Vi=Vi,1+Vi,2. max O (8.1)
Vi Vi ,1  Vi ,2  Vi ,3 , Vi ,1 , Vi ,2 , Vi ,3 t 0, i  :Gpv (3.1) subject to:
QGi ,lb d QGi d QGi ,ub , i  :Gpv (3.2)
PGi  PRi  PLi  Vi ¦ V j Yij cos Tij  G ij 0, i  : B (8.2)
QGi  QGi ,lb Vi ,2 =0, i  : Gpv (3.3)
QGi  QRi  QLi  Vi ¦ V j Yij sin T ij  G ij 0, i  : B (8.3)
QGi ,ub  QGi Vi ,3 =0, i  :Gpv (3.4)
Vi ,lb d Vi d Vi ,ub , i  : B (8.4)
B. Load Model
Pij ,lb d Pij d Pij ,ub , (i, j )  : B (8.5)
IMG loads, stated as a function of voltage magnitude and
frequency, are given as PGi ,lb d PGi d PGi ,ub , i  : Gd , : Gpv (8.6)
D pi
­
° PLi PLi 0 Vi Vi 0 ¬ª1  hpi (Z  Z0 ) ¼º QGi ,lb d QGi d QGi ,ub , i  : Gd , :Gpv (8.7)
® E qi
, i  :B (4)
°¯QLi QLi 0 Vi Vi 0 ª¬1  hqi (Z  Z0 ) º¼ (Z0  Z ) mpi  PGi t 0, i  : Gd (8.8)
In order to evaluate the maximum IMG loadability, bPi and (Vi 0  Vi ) nqi  QGi t 0, i  :Gd (8.9)
bQi are used to describe incremental loads; thus, PLi0 and QLi0
are expressed as PGi ,ub  PGi (Z0  Z ) mpi  PGi d H , i  : Gd (8.10)
­° PLi 0 PLi ,b  O bPi QGi ,ub  QGi (Vi 0  Vi ) nqi  QGi d H , i  :Gd (8.11)
® , i  :B (5)
°̄QLi 0 QLi ,b  O bQi
QGi  QGi ,lb Vi ,2 d H , QGi ,ub  QGi Vi ,3 d H , i  :Gpv (8.12)
where PLi,b and QLi,b are the base active and reactive loads,
which represent the current operating point. λ is the system Vi Vi ,1  Vi ,2  Vi ,3 , i  : Gpv (8.13)
load margin which describes the distance between existing and
Vi ,1 t 0,Vi ,2 t 0, Vi ,3 t 0, i  :Gpv (8.14)
critical operating points.
x
C. Modeling of Uncertainty Fi x ³f
pi t dt , i  :s (8.15)
1) Renewable-based DG: The probability distribution models
F x1 ,..., xn C F1 ( x1 ),..., Fn ( xn ) (8.16)
for renewable energy are designed as follows. A Weibull
model is adopted to describe the probability distribution of In the above formulation, the objective function is to
wind speed stated as maximize the IMG load margin. Constraints (8.2) and (8.3)
represent power balance. Constraints (8.4)-(8.7) are limits on
§ k ·§ x ·
k 1
§ § x ·k ·
p v ¨ ¸¨ ¸ exp ¨  ¨ ¸ ¸ , v t 0 (6) voltage magnitudes, feeder power flows and DG power
© c ¹© c ¹ ¨ ©c¹ ¸
© ¹ outputs. Constraints (8.8)-(8.11) describe the behavior of
droop-controlled DG units. Constraints (8.12)-(8.14) describe

1949-3053 (c) 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TSG.2018.2848958, IEEE
Transactions on Smart Grid
4

the switching between PV and PQ modes. (8.15) is the where N is the MCS sample size, k indicates the sample
probability distribution of variable xi and (8.16) is the Copula number for the variable vector, and z(k) and z'(k) represent two
function of correlated variables. Besides, (2.2), (2.4), (3.3) and sample vectors of z.
(3.4) are transformed to inequations in (8.10)-(8.12), and a
C. GSA with dependent variables
small positive number (ε=10-7) is introduced to guarantee a
strictly feasible region of the optimization problem [20]. When variables are correlated, the tight correspondence
between functional and variance decompositions is lost as the
III. GLOBAL SENSITIVITY ANALYSIS variance cannot be decomposed into a sum of squares as given
in (11) [23]. For correlated variables, the total variance of f(x)
In this section, GSA is developed for independent and
is decomposed as
dependent random variables, and the GS calculation is
discussed. D Dy ¬ª Ez f ( y, z ) ¼º  E y ¬ª Dz f ( y, z ) ¼º (17)
A. GSA with independent variables where
Assume x=(x1,…,xn) are uniformly distributed in an n- ­ Ez f ( y, z )
dimensional unit hypercube. Let Y=f(x1,…,xn) be a function ° ³Rns f ( y, z ) p( y, z y) dz
° 2
with n input variables represented as ® Dy ¬ª Ez f ( y , z ) ¼º ³R s ¬ª Ez f ( y , z ) ¼º p ( y )dy  f 0 (18)
2

n n
°
f ( x) f 0  ¦ fi ( xi )  ¦ f ij ( xi , x j )    f1...n ( x1 ,..., xn ) (9) ° f 0 ³ n f ( y, z ) p ( y, z ) dydz
i 1 i j ¯ R

Which is defined as the analysis of variation (ANOVA) where z is a random vector produced from the joint PDF p(y,z),
decomposition of f if the following condition is satisfied: and z is a random vector produced from the conditional
³f l1 ...lk xl1 ,..., xlk dxt 0, t l1 ,..., lk (10) probability distribution p( y, z | y ) .
The ANOVA decomposition is unique only if variables are Based on the definition of variance-based sensitivity, the
independent. Moreover, using (10) we have, GS of subset y is given as
n n Dy ª¬ Ez f ( y, z ) º¼
D ¦D ¦D
i 1
i
i j
ij  ...  D1...n (11) Sy
D
(19)

where D is the total variance of f and Dl1…lk is the variance of Substituting (18) into (19), the full expression of Sy is
fl1…lk(xl1,…,xlk), which are calculated as 1ª ª º p( y )dy  f 02 º (20)
2
Sy ³ ³
D «¬ R ¬ R
s n  s
f ( y , z ) p ( y , z y ) dz
¼ »¼
­D
³ f dx  f0
2 2
° MCS can be used to estimate Sy, but the direct application
® (12)
°̄ Dl1 ...lk ³ fl1 ...lk dxl1 ...dxlk of (20) requires costly computations of a double loop [24].
2

Therefore, a more efficient derivation [19] is stated as:

Consider an arbitrary set of s variables y=(xi1,…,xis), 1≤s<n, 1ª
p ( y )dy ª ³ ns f ( y, z ) p ( y, z y )dz
D ¬ ³R
and a complementary subset z=(xj1,…,xjn-s), such that x=(y,z). Sy s
¬ R
The variance corresponding to y is defined as (21)
º 2º
s u ³ ns f ( y , z ') p( y, z ' y )dz '  f 0
¼ ¼
Dy ¦ ¦
k 1 l1 ,..., lk I s
Dl1 ...lk , I s i1 ,..., is (13) R

D. GS calculation for correlated normal variables

The variance-based GS index Sy is defined to evaluate the
For n normal variables γ=(ξ1,…,ξn), consider a subset
impact of input variables y on the system output:
α=(ξi1,…,ξis), 1≤s<n and its complementary subset β. The
Dy mean vector and the covariance matrix of γ are partitioned as
Sy (14)
D ­P T
ªP P E º¼
°° J ¬ D
B. GS calculation for independent variables (22)
® ª 6D 6DE º
For independent uniformly distributed variables, the °6 J « »
subset’s variance Dy is stated as [12]: ¯° ¬ 6 ED 6 E ¼
The PDF of the n-dimensional normal variable vector γ is
Dy ³f x f y, z ' dxdz '  f 02 (15)
1 § 1 ·
In (J ) exp ¨  (J  P r )T 6J1 (J  P r ) ¸ (23)
Therefore, Sy can be estimated by MCS as (2S ) 6J
n © 2 ¹
­ 1 ª1 N
º
°S y
D «¬ N
¦ f y k ,z k f y k ,z' k  f 02 »
¼
For correlated normal variables, the conditional distribution
° k 1 In - s (D , E | D ) is also a normal distribution with the mean vector
°° 1 N μβc and the covariance matrix Σβc stated as
® f0 ¦ f y k ,z k
N k1
(16)
° ­° P E c PE  6 ED 6D1 D  PD
° 1 N 2
® (24)
°D
°̄
¦ f y k ,z k
N k1
 f 02 °̄6 E c 6 E  6 ED 6D16DE

1949-3053 (c) 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TSG.2018.2848958, IEEE
Transactions on Smart Grid
5

Therefore, for correlated normal variables, (21) becomes where a0, ai1, ai1i2 and ai1i2i3 are PCE coefficients, and Hi is the
1ª ith order Hermite polynomial (i=1,…,m).
Is (D )dD ª ³ ns f (D , E )In  s (D , E D )d E
D ¬ ³R
SD The number of coefficients in (31) is (n+m)!/(n!m!), which
s
¬ R
(25) grows by increasing n or m. Moreover, a sufficient number of
u ³ ns f (D , E ')In  s (D , E ' D )d E 'º  f 02 º samples, like 2-3 times the number of coefficients, is used to
R ¼ ¼
estimate these coefficients, making PCE computationally
E. GS calculation for correlated non-normal variables intractable for large uncertain problems. Sparse PCE (SPCE)
Let the marginal probability distribution of variable xi be Fi was proposed by Blatman and Sudret to address the PCE
(i=1,…,n). The Gaussian Copula is used to establish the joint problem when the output is dominated by lower order
probability distribution of x1,…,xn with the following formula: interactions between variables. The main idea of SPCE is to
apply a least angle regression-based truncation scheme to
C F1 ( x1 ),..., Fn ( xn ) ) n ) 1 ( F1 ( x1 )),..., ) 1 ( Fn ( xn )); 6 (26) detect significant PCE coefficients and retain a smaller
The PDF of the Gaussian Copula satisfies the following number of PCE terms. The theory of truncation scheme is
presented in [25] as we use the UQLab toolbox to construct
In ) 1 ( F1 ( x1 )),..., ) 1 ( Fn ( xn )); 6 the SPCE model in this paper [26].
p F1 ( x1 ),..., Fn ( xn ) (27)
–
n
i 1
I ) 1 ( Fi ( xi )) B. Procedure of GSA
According to the transformation law of probabilities (27), The probabilistic nonlinear optimization problem (8) is
the GS for correlated non-normally distributed variables is expressed in a compact form as Y=f(x), where x are n random
presented as [19] variables, and Y is the IMG load margin. The GSA procedure
1ª for evaluating the impact of variables on the load margin is
ª
Sy ³
D ¬« R
I ss
(D ) d D ³
¬ Rn s
f + s D , + n  s E I n  s (D , E D ) d E given as follows.
(28) Step 1: Obtain the IMG parameters and the probability models
u ³ f + s D , + n  s E ' In  s (D , E ' D )d E 'º  f 02 º» for random variables
R n s
¼ ¼
where Step 2: Select a set of variables, defined as y, from n input
variables x
­+ s D F 1
( ) ([ )),..., F 1
( ) ([ )) Step 3: Establish surrogate models of load margins
° i1 i1 is is
® (29) a. Select Nc samples of independent standard normal
°̄+ n  s E Fj1 1 () ([ j1 )),..., F jn1s () ([ jns ))
variables as collocation points of input variables
Based on (28), Sy is estimated by MCS as b. Transform collocation points into samples with a joint
1 ª1 N distribution p(y,z) or a conditional distribution p ( y, z |y )
2º
Sy
D «¬ N k 1
¦ f y k
, z k
f y k
, z ' k
 f 0 »
¼
(30) c. Solve the deterministic optimization problem (8) and
obtain load margin samples
The application of (30) requires: 1) samples of joint and
d. Calculate SPCE coefficients based on the collocation
conditional probability distributions, 2) evaluation of system
points and load margin samples
output based on inputs. The procedure to perform GSA for the
e. Establish the surrogate models of f ( y, z ) or f ( y, z )
maximum IMG loadability is given in the following section.
Step 4: Obtain surrogate model output using MCS
IV. GSA FOR MAXIMUM IMG LOADABILITY In order to improve the MCS efficiency, the Sobol sequence
is used to obtain the surrogate model output using the
MCS is usually used to perform GSA which requires many
following steps:
samples to obtain accurate results, which limits the practical
a. Generate 2n-dimensional Sobol sequences with a sample
GSA applications. In this paper, an efficient GSA method is
size of N
proposed to analyze the impacts of variables on the loadability
b. Transform the first n-dimensional Sobol sequences into
of IMGs. Instead of directly calculating GSs with MCS, a
samples of standard normal variables, and obtain the
surrogate model for the load margin is first established with
surrogate model output of f(y,z) using the samples
PCE, and then GS is calculated using the surrogate model.
c. Combine the first s-dimensional with the last n-s-
A. Polynomial Chaos Expansion dimensional Sobol sequences, transform them into
PCE is used to represent the random model output by a set samples of standard normal variables, and obtain the
of coefficients in polynomial chaos bases. For n independent surrogate model output of f ( y , z ) using the samples
standard normal variables [ˆ1 ,..., [ˆn , PCE is stated as Step 5: Calculate the GS using (30) and assess the importance
n n i1
of variables y.
Y a0  ¦ ai1 H1 [ˆi1  ¦¦ ai1i2 H 2 [ˆi1 , [ˆi2
i1 1 i1 1 i2 1 V. CASE STUDIES
i1 i2
(31)
n
In this section, the proposed GSA method is tested using a
 ¦¦¦ ai1i2i3 H 3 [ˆi1 , [ˆi2 , [ˆi3  ...
i1 1 i2 1 i3 1
33-bus IMG. The variabilities of loads, feeder parameters and
renewable power are modeled as random variables, and the
correlation of renewable power outputs is considered. The

1949-3053 (c) 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TSG.2018.2848958, IEEE
Transactions on Smart Grid
6

proposed method is compared with other GSA and LSA 3-4 PV Beta a=2.06, b=2.5
methods for evaluating the impacts of variables. Programs are 5-36 Load Normal σ/μ=0.2
37-73 Feeder parameter Uniform lb=0.9, ub=1.1
developed with Matlab 8.1 on a PC with Intel Core i7 3.6GHz
CPU and 8GB of RAM. Sobol sequences are used to generate
samples of random variables. The OPTI toolbox [27] is used B. Surrogate Model
to solve the nonlinear optimization problem (8). The UQLab In this section, the surrogate models of load margins are
toolbox is used to calculate SPCE coefficients. established by PCE and SPCE. The 2nd-order PCE with
A. IMG System Description accurate results in power system simulations [30] is used in
our simulations. Note that higher order PCEs can also be used.
Fig. 1 shows the diagram of the 33-bus IMG system. The
There are 73 variables in the test system; thus the number of
feeder parameters and nominal load power of the system are
coefficients of the 2nd-order PCE is 2,775. Accordingly, 8,000
given in [28]. Three dispatchable DG units and four
collocation points are used to estimate the coefficients using a
renewable-based DG units are allocated to feed the IMG
least square method [17]. The 2nd-order PCE is also applied
system, and the parameters of DG units are given in Tables I
and 1,000 collocation points are used to estimate the SPCE
and II [11]. The load parameters are designed as αp=0.92,
coefficients. For comparison, MCS with a sample size of
βq=4.04, hp=0.8 and hq=-2.0 [29]. The upper and lower limits
50,000 is used to solve the probabilistic nonlinear optimization
of voltage magnitude are 1.05 and 0.95 p.u.. The upper and
problem (8).
lower limits of the system frequency are 1.005 and 0.995 p.u..
The probability distributions of load margins are given in
The probability models of input random variables are given
Fig. 2. The deterministic load margin is 0.49. However, if
in Table III. The probability distribution models of wind speed
input random variables are considered, based on SPCE results,
and solar irradiance are described by Weibull and Beta
the probability that the load margin is larger than 0.7 is 11.4%
distributions, respectively [30]. The active power load has a
and the probability that the load margin is smaller than 0.3 is
normal distribution with a standard deviation which is 20% of
16.6%. Therefore, input variables affect the variability of load
its mean. The reactance of a feeder l is in [0.9Rl, 1.1Rl], where
margins significantly. Moreover, CDF of load margins
Rl is the base reactance value. The resistance is modeled to
obtained by SPCE, PCE and MCS are nearly the same, which
maintain a constant R/X. Moreover, the power outputs of the
validates the accuracy of surrogate models.
same type of renewable-based DG units have a rank
correlation coefficient of 0.5 for the Gaussian Copula. There
are 73 variables in IMG including independent and correlated
variables.

Fig. 2 CDFs of load margin.

Fig. 1 33-bus IMG test system.
TABLE I
Parameters of droop-controlled DG units (Sbase = 1MVA)
PG,ub QG,ub
Bus mp nq V0 ω0
(p.u.) (p.u.)
4 3 1.5 7.29E-4 1.65E-5 1.0009 1.016
8 0.5 1 4.03E-3 6.17E-2 1.0008 1.016
22 1.5 1 1.22E-3 3.21E-2 1.0009 1.040

TABLE II
Parameters of renewable-based DG units (Sbase = 1MVA)
Bus Type Capacity (p.u.) Power factor
8 PV 0.75 0.9
14 PV 0.75 0.9
28 Wind 0.75 0.9
30 Wind 0.75 0.9 Fig. 3 Coefficients of PCE and SPCE.
TABLE III The coefficients of PCE and SPCE are given in Fig. 3.
Models of input random variables Although there are 2,775 coefficients in PCE and only 140
Index Variable Model Parameter
1-2 Wind Weibull c=7.41, k=2.06
non-zero coefficients in SPCE, the coefficients with absolute

1949-3053 (c) 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TSG.2018.2848958, IEEE
Transactions on Smart Grid
7

values larger than 0.01 are almost the same, making the results variability with a GS that is 0.656 for GSA-PCE and 0.659 for
of PCE and SPCE consistent. Compared with PCE, SPCE GSA-SPCE. When the renewable power uncertainty is
obtains accurate results with much fewer coefficients for high- considered, the load impact decreases. The wind power is the
dimensional and uncertain problems and thus fewer most influencing factor, followed by the load at bus 30 and PV
collocation points are needed to estimate coefficients. power. The GSs of feeder parameters are smaller than 0.01.
GSA is also compared with the commonly used LSA [15]
C. Sensitivity Analysis
in assessing the importance of variables and 6 variables with
In this section, the proposed method is used to calculate largest sensitivities are given in Table V. The critical variables
sensitivities of random variables with respect to IMG load among independent variables identified by LSA and GSA are
margins. The performances of GSA-MCS, GSA-PCE and the same, but the rankings of variables are different. This is
GSA-SPCE are compared. The GS of wind power at bus 28 is because local sensitivities (LSs) are based on the partial
given in Fig. 4. The three methods obtain similar GSs with derivative of the output with respect to the input while GSs
large sample sizes. The result of GSA-PCE is closer to that of consider nonlinear systems and apply the probability models
GSA-MCS, since full polynomial chaos bases are used in PCE of variables. Hence GSA is more accurate for evaluating the
while only significant polynomial chaos bases are retained in impact of variables. Besides, LSA cannot be used for
SPCE. Nevertheless, GSA-SPCE also obtains relatively correlated variables and LSs do not provide any credible
accurate GSs. It is found that GSA obtains stable results when information on the segment of output uncertainty which is
more than 10,000 samples are considered. Therefore, GSA attributed to a subset of inputs.
suffers from heavy computational burdens if GSs are
calculated based on the original problem. The sample size and
computation time of corresponding methods are given in
Table IV. Although GSA-PCE and GSA-SPCE require
additional time for establishing surrogate models, they are
much more efficient than GSA-MCS since GSs are calculated
based on surrogate models. GSA-SPCE is more efficient than
GSA-PCE because fewer samples are used to estimate
surrogate model coefficients.

Fig. 5 GSs of random variables.

TABLE V
Sensitivities of variables obtained by LSA and GSA
Fig. 4 Relation between GS and sample size. GSA-SPCE
LSA
Ranking Independent Correlated
TABLE IV Variable LS Variable GS Variable GS
Sample sizes and computation time of different GSA methods 1 2 0.225 2 0.312 2 0.474
GSA- GSA- GSA- 2 33 0.170 1 0.213 1 0.403
MCS PCE SPCE 3 1 0.166 33 0.210 33 0.157
Establishing Sample size - 8000 1000 4 35 0.075 4 0.087 4 0.131
surrogate models Time (s) - 684.7 89.8 5 3 0.057 35 0.034 3 0.093
Sample size 50000 50000 50000 6 4 0.055 3 0.032 35 0.026
Calculating GSs
Time (s) 8796.4 1.5 0.6
Total computation time (s) 8796.4 686.2 90.4 Finally, the impact of the correlations of input variables is
investigated. The rank correlation of Gaussian Copula varies
The GS of each input variable is given in Fig. 5. There are
between -1 and 1, and GSs for various types of uncertainties
73 variables in total and only those with GSs larger than 0.01
are given in Fig. 6. The GS of wind and PV power increases
are presented. In Fig. 5, the results obtained by GSA-SPCE
and the GS of loads decreases as we increase the correlation
and GSA-PCE are nearly the same for all variables. Moreover,
coefficient of renewable power. The GS of feeder parameters
only a small set of variables have relatively large GSs, and
remains small in all simulation cases.
most of the variables have minor impacts on load margins.
The GS of wind power with various correlations are given
When the renewable generation output is assumed to be
in Fig. 7. For independent wind power, the sum of the GS of
deterministic, the 33rd input variable, which is the load
each wind power is almost equal to the GS of the total wind
demand at bus 30, accounts for most of the load margin

1949-3053 (c) 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TSG.2018.2848958, IEEE
Transactions on Smart Grid
8

power, which indicates that the wind power interactions nearly margin. The DG unit outputs in cases 1 and 2 are similar. The
have no impact on load margins. For positively or negatively reactive power of unit 1 at critical operating points is always
correlated wind power, the sum of the GS of each wind power equal to its maximum, and the active or reactive power of
is larger or smaller than the GS of total wind power. This is other units varies within certain ranges.
because of the variable attributes and its correlations with
other variables are considered for calculating GS. The sum of
the GS of each wind power will similarly consider the impact
of correlations.

Fig. 8 Cumulative distributions of load margin when different variables are

considered.

Fig. 6 GS of different types of correlated uncertainties.

Fig. 7 GS of each wind power and the total wind power. Fig. 9 Bus voltage magnitudes at critical operating points when different
variables are considered.
D. Maximum Loadability Analysis
In order to validate the GSA results, we consider the
following three cases and evaluate the probabilistic value of
the maximum IMG loadability.
Case 1: All 73 input random variables are considered.
Case 2: 7 variables with GSs larger than 0.01 are considered.
Case 3: 66 variables with a GS smaller than 0.01 are
considered.
Fig. 8 shows that if we merely consider the 7 variables with
large GSs, the load margin CDFs will almost be the same.
However, if we only consider the 66 variables with small GSs, Fig. 10 DG unit outputs at critical operating points when different variables
the load margin CDFs will be quite different although most are considered.
variables are considered. This result indicates that most E. Influence of PHEV
variables have a minute impact on the load margin, and the
proposed method accurately identifies the most effective As the number of plug-in hybrid electric vehicles (PHEVs)
variables. increases rapidly, the charging demand has become an
The bus voltage magnitudes and droop-controlled DG unit important part of the electricity load. In the following, we use
outputs at critical operating points are given in Figs. 9 and 10, the PHEV model in [31] to evaluate the load demand
respectively. The minimum voltage magnitude is 0.95 which variability of charging stations. The method establishes a
reaches the lower limit. The maximum voltage magnitude is single PHEV charging demand model and employs queuing
1.03 at bus 22 which is smaller than the upper limit. For some theory to describe the behavior of multiple PHEVs. Finally,
buses, voltage magnitude variations are different for cases 1 the probability model of the total demand of a PHEV charging
and 2, which indicates that although some variables have an station is obtained with the maximum likelihood estimation.
impact on voltage magnitudes, they do not affect the load We consider two cases to analyze the GS of PHEVs.
S1: All loads are described by normal distributions.

1949-3053 (c) 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TSG.2018.2848958, IEEE
Transactions on Smart Grid
9

S2: The loads at buses 22 and 29 are replaced by PHEV reduce the number of variables in IMG analyses and help
charging stations with models shown in Table VI. manage the impact of variabilities with proper control
The GS of each variable is calculated and load demands strategies, such as the installation of energy storage systems.
with GSs larger than 0.01 are given in Fig. 11. When load
demands are described by normal distributions, the load VII. REFERENCES
demand at bus 30 accounts for most of the load margin [1] S. Parhizi, H. Lotfi, A. Khodaei, and S. Bahramirad, “State of the art in
variability with a GS equal to 0.672. When the integration of research on microgrids: A Review,” IEEE Access, vol. 3, pp. 890-925,
2015.
PHEV charging stations is considered, the GS of load demand
[2] Z. Wang, B. Chen, J. Wang, and J. Kim, “Decentralized energy
at bus 30 decreases and the GS of load demand at bus 29 management system for networked microgrids in grid-connected and
increases. In addition, the GS of the PHEV charging demand islanded modes,” IEEE Trans. Smart Grid, vol. 7, no. 2, pp. 1097-1105,
at bus 22 remains small, although the PHEV charging demand 2016.
at bus 22 has the same probability model as that at bus 29. [3] L. Che, X. Zhang, M. Shahidehpour, A. Alabdulwahab, and Y. Al-Turki,
“Optimal planning of loop-based microgrid topology,” IEEE Trans.
Therefore, compared with load demands described by normal Smart Grid, vol. 8, no. 4, pp. 1771-1781, 2017.
distributions, the load demand of PHEV charging stations in [4] J. Li, Y. Liu, and L. Wu, “Optimal operation for community based multi-
this case has a higher impact on the maximum IMG party microgrid in grid-connected and islanded modes,” IEEE Trans.
loadability. The influence of charging demand will decrease Smart Grid, 2017.
by selecting suitable locations for PHEV charging stations. [5] M. A. Allam, A. A. Said, M. Kazerani, and E. F. El-Saadany, “A novel
dynamic power routing scheme to maximize loadability of islanded
TABLE VI hybrid ac/dc microgrids under unbalanced ac loading,” IEEE Trans.
Load demands of PHEV charging stations Smart Grid, 2017.
Weibull distribution [6] H. R. Baghaee and M. Abedi, “Calculation of weighting factors of static
Mean power security indices used in contingency ranking of power systems based on
Bus Scale Shape Power factor
(MW) fuzzy logic and analytical hierarchy process, ” Int. J. Elect. Power Energy
parameter parameter
22 0.192329 1.87103 0.17076 0.95 Syst., vol. 33, no. 4, pp. 855-860, 2011.
29 0.192329 1.87103 0.17076 0.95 [7] H. R. Baghaee, M. Mirsalim, G. B. Gharehpetian, and H. A. Talebi,
“Reliability/cost-based multi-objective Pareto optimal design of stand-
alone wind/PV/FC generation microgrid system,” Energy, vol. 115, part
1, no. 1, pp. 1022-1041, 2016.
[8] H. R. Baghaee, M. Mirsalim, G. B. Gharehpetian, and A. K. Kaviani,
“Security/cost-based optimal allocation of multi-type FACTS devices
using multi-objective particle swarm optimization,” Simulation: Int.
Trans. Society for Modeling and Simul., vol. 88, no. 8, pp. 999-1010,
2012.
[9] G. Diaz and C. Gonzalez-Moran, “Fischer-burmeister-based method for
calculating equilibrium points of droop-regulated microgrids," IEEE
Trans. Power Syst., vol. 27, no. 2, pp. 959-967, 2012.
[10] M. M. A. Abdelaziz and E. F. El-Saadany, “Maximum loadability
Fig. 11 GS of load demand at each bus. consideration in droop-controlled islanded microgrids optimal power
flow,” Electric Power Syst. Res., vol. 106, pp. 168-179, 2014.
[11] M. M. A. Abdelaziz, H. E. Farag, and E. F. El-Saadany, “Optimum droop
VI. CONCLUSIONS parameter settings of islanded microgrids with renewable energy
In this paper, a GSA method is proposed to evaluate the resources,” IEEE Trans. Sustain. Energy, vol. 5, no. 2, pp. 434-445, 2014.
impact of variables on the maximum IMG loadability. [12] I. M. Sobol, “Global sensitivity indices for nonlinear mathematical
models and their Monte Carlo estimates,” Mathematics and computers in
Considering MG characteristics operating in the islanded simulation, vol. 55, no. 1, pp. 271-280, 2001.
mode and the probability distribution models of uncertainties, [13] F. Fernández-Navarro, M. Carbonero-Ruz, D. Becerra Alonso, and M.
a probabilistic optimization problem is established to evaluate Torres-Jiménez, “Global Sensitivity Estimates for Neural Network
IMG load margins. Then, GSs are presented to evaluate the Classifiers,” IEEE Transactions on Neural Networks and Learning
Systems, vol. PP, no. 99, pp. 1-13.
importance of independent and correlated variables. In order
[14] D. M. King and B. J. C. Perera, “Morris method of sensitivity analysis
to improve the efficiency of GSA, SPCE is used to establish applied to assess the importance of input variables on urban water supply
probability models of load margins, and GSs are calculated yield – A case study,” Journal of Hydrology, vol. 477, pp. 17-32, 2013.
based on surrogate models. [15] R. Preece and J. V. Milanović, “Assessing the applicability of uncertainty
The proposed method is tested using a 33-bus IMG system importance measures for power system studies,” IEEE Trans. Power
Syst., vol. 31, no. 3, pp. 2076-2084, 2016.
which is compared with other SA methods. The simulation [16] K. N. Hasan, R. Preece, and J. V. Milanović, “Priority ranking of critical
results indicate that SPCE overcomes the PCE deficiency in uncertainties affecting small-disturbance stability using sensitivity
high-dimensional problems by reducing the number of analysis techniques,” IEEE Trans. Power Syst., vol. 32, no. 4, pp. 2629-
2639, 2017.
polynomial chaos bases, thus much fewer coefficients are
[17] F. Ni, P. H. Nguyen, and J. F. G. Cobben, “Basis-adaptive sparse
estimated in establishing surrogate models. Compared with polynomial chaos expansion for probabilistic power flow,” IEEE Trans.
LSs, GSs consider the power system nonlinearity and describe Power Syst., vol. 32, no. 1, pp. 694-704, 2017.
the impact of variable dependence on the output. GSA-SPCE [18] F. Ni, M. Nijhuis, P. H. Nguyen, and J. F. G. Cobben, “Variance-based
obtains accurate GSs and is much more efficient than GSA- global sensitivity analysis for power systems," IEEE Trans. Power Syst.,
vol. PP, no. 99, pp. 1-1.
MCS and GSA-PCE. For IMGs with various uncertainties,
accurate identification of the most influencing variables will

1949-3053 (c) 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TSG.2018.2848958, IEEE
Transactions on Smart Grid
10

[19] S. Kucherenko, S. Tarantola, and P. Annoni, “Estimation of global Jiao Tong University, Shanghai, China. His research interests include demand
sensitivity indices for models with dependent variables,” Comput. Phys. response, transactive energy system, and blockchain and machine learning
Commun., vol. 183, no. 4, pp. 937-946, 2012. applications in power systems. He is a recipient of Shanghai Pujiang Talent
[20] W. Rosehart, C. Roman, and A. Schellenberg, “Optimal power flow with Program.
complementarity constraints,” IEEE Trans. Power Syst., vol. 20, no. 2, pp.
843–822, 2005.
Han Wang received the B.S. degree from the College of Electrical
[21] B. Das, “Consideration of input parameter uncertainties in load flow
Engineering, Zhejiang University, Hangzhou, China, in 2015. He is currently
solution of three-phase unbalanced radial distribution system,” IEEE
working toward the Ph.D. degree at the Department of Electrical Engineering,
Trans. Power Syst., vol. 21, no. 3, pp. 1088-1095, 2006.
Shanghai Jiao Tong University, Shanghai, China. His research interests
[22] G. Papaefthymiou and D. Kurowicka, “Using copulas for modeling include power system uncertainty quantification and microgrid control.
stochastic dependence in power system uncertainty analysis,” IEEE
Trans. Power Syst., vol. 24, no. 1, pp. 40-49, 2009.
[23] E. Borgonovo and E. Plischke, “Sensitivity analysis: a review of recent Zhiyi Li (GSM’14-M’17) received the B.S. degree in electrical engineering
advances,” European Journal of Operational Research, vol. 248, no. 3, from Xi’an Jiaotong University, Xi’an, China, in 2011 and the M.S. degree in
pp. 869-887, 2016. electrical engineering from Zhejiang University, China, in 2014. He
[24] A. Saltelli and S. Tarantola, “On the relative importance of input factors completed his Ph.D. degree in 2017 in the Electrical and Computer
in mathematical models: safety assessment for nuclear waste disposal,” Engineering Department at Illinois Institute of Technology. He is a visiting
Journal of the American Statistical Association, vol. 97, no. 459, pp. 702- faculty in the Robert W. Galvin Center for Electricity Innovation at Illinois
709, 2002. Institute of Technology.
[25] G. Blatman and B. Sudret, “Adaptive sparse polynomial chaos expansion
based on least angle regression,” Journal of Computational Physics, vol.
230, no. 6, pp. 2345-2367, 2011. Quan Zhou (GSM’16) received the B.S. and M.S. degrees from the
Department of Electrical Engineering, Shanghai Jiao Tong University,
[26] S. Marelli and B. Sudret, “UQLab: A framework for uncertainty
Shanghai, China, in 2011 and 2016, respectively. He is currently pursuing the
quantification in Matlab,” in Proc. 2nd Int. Conf. Liverpool.
Ph.D. degree in the Electrical and Computer Engineering Department, Illinois
Vulnerability, Risk Anal. Manag., 2014, 2554-2563.
Institute of Technology. His research interests include microgrid control and
[27] J. Currie and D. I. Wilson, "OPTI: Lowering the barrier between open optimization.
source optimizers and the industrial MATLAB user," Foundations of
Computer-Aided Process Operations, Georgia, USA, 2012.
[28] M. E. Baran and F. F. Wu, “Network reconfiguration in distribution
systems for loss reduction and load balancing,” IEEE Trans. Power Del.,
vol. 4, no. 2, pp. 1401-1407, 1989.
[29] P. Li, Y. Shi, Z. Wu, M. Hu, X. Wang, and W. Shi, “Power flow
calculation method similar to Benders decomposition for islanded
microgrid,” Automation of Electric Power Systems, vol. 41, no. 14, pp.
119-125, 2017. (Chinese)
[30] Z. Ren, W. Li, R. Billinton, and W. Yan, “Probabilistic power flow
analysis based on the stochastic response surface method,” IEEE Trans.
Power Syst., vol. 31, no. 3, pp. 2307-2315, 2016.
[31] G. Li and X. P. Zhang, “Modeling of plug-in hybrid electric vehicle
charging demand in probabilistic power flow calculations”, IEEE Trans.
Smart Grid, vol. 3, no. 1, pp. 492-499, 2012.

BIOGRAPHIES
Xiaoyuan Xu (S’15-M’16) received the B.S. and Ph.D. degrees in electrical
engineering from Shanghai Jiao Tong University, Shanghai, China, in 2010
and 2016, respectively. He is currently a Postdoctoral Researcher with
Shanghai Jiao Tong University. He is also a visiting scholar with the Illinois
Institute of Technology, Chicago, IL, USA. His research interests include
power system uncertainty quantification and power system optimization.

Zheng Yan received the B.S. degree from Shanghai Jiao Tong University,
Shanghai, China, in 1984 and the M.S. and Ph.D. degrees from Tsinghua
University, Beijing, China, in 1987 and 1991, respectively, all in electrical
engineering. He is currently a Professor and the Director of Electrical
Engineering with Shanghai Jiao Tong University. His current research
interests include application of optimization theory to power systems, power
markets, and dynamic security assessment.

Mohammad Shahidehpour (F’01) received the Honorary Doctorate degree

in electrical engineering from the Polytechnic University of Bucharest,
Bucharest, Romania. He is the Bodine Chair Professor and the Director with
the Robert W. Galvin Center for Electricity Innovation, Illinois Institute of
Technology, Chicago, IL, USA. He is a member of the US National Academy
of Engineering and a Fellow of the American Association for the
Advancement of Science (AAAS).

Sijie Chen (S’10-M’15) received the B.E. and Ph.D. degrees from Tsinghua
University, Beijing, China, in 2009 and 2014, respectively. He is currently an
Assistant Professor with the Department of Electrical Engineering, Shanghai

1949-3053 (c) 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.