A novel probabilistic small signal stability analysis (PSSSA) method considering load correlation is proposed in this paper. The superiority Latin hypercube sampling (LHS) technique combined with Monte Carlo simulation (MCS) is utilized to investigate the probabilistic small signal stability of power system in presence of load correlation. LHS helps to reduce the sampling size, meanwhile guarantees the accuracy and robustness of the solutions. The correlation coefficient matrix is adopted to represent the correlations between loads. Simulation results of the twoarea, fourmachine system prove that the proposed method is an efficient and robust sampling method. Simulation results of the 16machine, 68bus test system indicate that load correlation has a significant impact on the probabilistic analysis result of the critical oscillation mode under a certain degree of load uncertainty.
1. Introduction
Power system always faces various small disturbances. Investigating the stability characteristics of dynamic systems under the action of small disturbance events, plays an important role in improving the overall security and reliability of the power grid. Traditional small signal stability analysis (SSSA) methods focus on some special or the worst operating condition case. The main limitation of those methods is that the impact of various uncertainties on the stability of power system cannot be demonstrated comprehensively.
Uncertainty, however, is the inherent nature in a real world physical system. In complex power system, uncertainties come from many aspects, such as volatility of load, uncertainty of line parameter, change of generator power, zero drift of controller parameter, model approximation, and so on. These uncertainty factors, under the proper conditions, may gently push the system into an unstable operating mode and affect the power system small signal stability in a certain extent. Therefore, probabilistic small signal stability analysis (PSSSA) considering various uncertainties attracts the interest of the researchers in the past few decades. Furthermore, with the rapid development of smart grid in recent years, to cope with the uncertainties, the PSSSA is a key issue
[1]
.
Reference
[2]
has introduced the probability concept into the power system dynamic stability analysis. The system’s joint instability probability was evaluated by using the eigenvalue ensitivity to calculate the statistic characteristics of eigenvalue locations in the complex plane. Up to now, the PSSSA methods of power system can be classified into analytical methods and numerical methods. The analytical methods include point estimate method
[3
,
4]
, cumulant method
[5]
, and probabilistic collocation method (PCM)
[6
,
7]
. These methods usually need to approximately simplify the complex nonlinear functional relationship between the eigenvalues and input random variables, and the calculation error is inevitable. In addition, these methods, based on the known probability distribution of random variables, have difficulty to handle the unknown empirical distribution situation. Monte Carlo simulation (MCS) based numerical method is a generic way to deal with all kinds of uncertainties, but it requires a large number of repetitive simulations to obtain accurate solutions. In
[8

10]
, the computation time of MCSbased PSSSA is studied and numerical results prove that a large amount of time is required to calculate the probability distribution function (PDF) of the system critical oscillation mode. Especially in a largescale power system, this method tends to be unacceptable.
On the other hand, from the viewpoint of modeling, the existing PSSSA research considers the nodal power injections as independent random variables normally
[11]
. In fact, different nodal power injections always exhibit correlation in a certain degree due to the impact of the consumer rest rules, geographical distribution, climate change and other factors. As a result, the correlations between different nodal loads should be duly taken into account in order to approximate to the actual situation. Billinton and Huang
[12]
study the system reliability index considering the influence of nodal load correlation. The effects of nodal load uncertainty and correlation in system adequacy assessment were illustrated in
[13]
. These studies show that the results are more accurate and credible when the nodal load correlation is considered in the probabilistic analysis of power system. Therefore, for the same reason, the SSSA study has to consider load uncertainty and correlation simultaneously.
Compared to the traditional simple random sampling (SRS) technique, Latin hypercube sampling (LHS) technique can greatly improve sampling efficiency, and ensure the sample values covering the entire distribution range of the input random variable. This technique has been applied in probabilistic load flow calculation
[14
,
15]
, and reliability evaluation in complex power system
[16

18]
. The results show that this method has efficient sampling property to deal with unknown empirical distribution random variables and can handle the correlation between multivariate input random variables conveniently.
In this paper, a LHSbased PSSSA method considering load correlation is proposed. The rest of this paper is organized as follows. Section 2 gives the fundamental concept of the PSSSA and the modeling of system uncertainties and load correlation. Section 3 presents the basic steps of LHSbased MCS. The processing approach of correlation between different loads is also discussed. The advantageous performance of the proposed method, taking load correlation into account, is investigated by two typical test systems in Section 4. Section 5 summarizes the main conclusion of this paper and points out the future work.
2. Fundamental Theory
The basic principle, modeling, and method of PSSSA are described in this section.
 2.1. Mathematical Model of SSSA
Multimachine power system dynamic behavior is represented by a set of nonlinear differentialalgebraicequations (DAE):
where
x
are the dynamic state variables,
y
are the algebraic variables. The differential equations
f
consist of the dynamic equations of the equipment such as generators, exciters, and governors, while the algebraic equations
g
include the network power flow equations.
When the system operates in a given initial equilibrium point (
x
_{0}
,
y
_{0}
), the full system linearization algebraicdifferentiation equations (2) can be obtained by Taylor series expansion by ignoring the higherorder terms.
where A = ▽
_{x}f
, B = ▽
_{y}f
, C = ▽
_{x}g
, and D = ▽
_{y}g
where. Then, the system state Eq. (3) is obtained by eliminating the algebraic variable Δ
y
where A
_{s}
= A−BD
^{−1}
C is the system state matrix, thus implicitly assuming that D is nonsingular.
The eigenvalues
λ
= [
λ
_{1}
,
λ
_{2}
, …,
λ_{n}
] of the system state matrix A
_{s}
are calculated by solving the Eq. (4).
The eigenvalues may be real or complex. If A
_{s}
is real matrix, complex eigenvalues always occur in conjugate pairs. Each pair of the complex eigenvalues represents one oscillation mode of the power system.
According to the Liapunov’s indirect method, a dynamical system with small signal stability corresponding every eigenvalue’s real part is less than zero, and vice versa. As a consequence, the distribution character of oscillation modes located in the complex plane reveals the small signal stability of the power system.
 2.2. PSSSA and probabilistic index
Probabilistic assessment of power system small signal stability is to evaluate the instability probability of all oscillation modes under all kinds of input random variables. When the probability density distribution of input random variable is known, the probability distribution of every oscillation mode can be obtained by repetitive modal analysis.
In order to quantificationally evaluate the instability probability of an oscillation mode, a probabilistic instability index (PIstI) is introduced. For a specified oscillation mode corresponding
λ_{i}
=
α_{i}
±j
β_{i}
, the PIstI can be calculated as
where
f
(
α_{i}
) is the PDF of the output variable
α_{i}
.
Another commonly used indicator in SSSA is mode damping ratio
In the practical operation of power system, even if all of the oscillation modes are stable, a too small mode damping ratio is unallowed, simply because this operating condition is considered to be insecure. A mode damping ratio greater than 5% (even 3% in a largescale interconnected power grid) is accepted to achieve preferable dynamic performance after a disturbance event
[19]
.
Apparently, damping ratio of an oscillation mode is influenced by input random variables. To estimate probabilistic insecurity margin, the probabilistic insecurity index (PIseI) is introduced and can be calculated as
where
f
(
ζ_{i}
) is the PDF of mode damping ratio
ζ_{i}
.
 2.3. System uncertainties modeling
The assessment and presentation of the effects of nodal power injection uncertainties, including the system generators and loads power injection, on small signal stability are now widely recognized as important parts of power system probabilistic security analysis
[11]
. From the practical point of view, nodal power possesses random volatility all the time. And the short time volatility characteristics of nodal power mainly consist of many stochastic factors. Thus, it is inappropriate to adopt the deterministic nodal power injection model.
System uncertainties are represented for all the PV generator output powers
P_{e}
and the system PQ load powers
P_{L}
,
Q_{L}
in the this paper. In the simplest case, the uncertainty of each nodal generator output and load power are both described as a normal distribution with whoes mean values of corresponding normal distributions are set as the nominal values. Different levels of standard deviation are considered here to denote the levels of system parameter uncertainty. For nodal power
p_{i}
which belongs to generator or load power vector P=[
p
_{1}
,
p
_{2}
,…,
p
_{n}
]
^{T}
, the PDF
f
(
p
_{i}
) can be expressed as:
where
μ_{i}
and
σ_{i}
^{2}
are mean and variance of nodal power
p_{i}
respectively. This description of uncertainty can adapt to many kinds of system models.
It is important to note that although the uncertainties considered in this work are assumed to be normal distributions, the proposed sampling technique in the next section is also adapted to any other distribution.
 2.4. Load correlation modeling
Load correlation comes from the intrinsic properties of load demand. The significant impact of load correlation on power system probabilistic load flow calculation has been proved in
[14
,
15
,
20]
. Since different initial conditions of the power flow influence the eigenvalues of linearization power system, especially when the system state matrix has a large condition number, the result of PSSSA is also affected by taking load correlation into account.
However, the relationships between multiple input variables are complicated. In most situations, different nodal load behavior reflects a certain degree of correlation. It is not reasonable to consider the nodal load to be completely independent or completely correlated corresponding correlation coefficient 0 or 1.
In this paper, a general simplified approach named correlation coefficient matrix is used to describe the nodal load correlation in a short time scale under certain load uncertainty level. Such an approach can easily describe the linear correlation between the different nodal load powers.
Assume
C
_{Z}
is the correlation coefficient matrix of input random variable vector
Z
=[
z
_{1}
,
z
_{2}
, … ,
z
_{n}
]
^{T}
. The correlation coefficient
ρ_{wij}
can be calculated by (10).
where
σ_{zi}
and
σ_{zj}
are the standard deviations of
z
_{i}
and
z
_{j}
,
Cov
(
z_{i}, z_{j}
) is the covariance of
z
_{i}
and
z
_{j}
, and
ρ_{zy}
is the correlation coefficient of w
_{i}
and w
_{j}
.
For a series of ideal independent random variables, the correlation coefficient matrix is a standard unit matrix. But in most situations, C
_{Z}
is a symmetric matrix with
ρ_{zy}
∈(−1, 1), where
ρ_{zy}
∈(−1, 0) represents negative correlation. In this paper, we assume that the correlation coefficient matrix C
_{Z}
is positive definite.
 2.5. SRSMCS based PSSSA
Simple random sampling Monte Carlo simulation (SRSMCS) is a widely used computational technique to characterize the random behavior of a physical system. This technique, which can provide comprehensive information, has been used in PSSSA as a reference standard for other methods. Although SRSMCS can obtain high accuracy results, this method is timeconsuming. The computational burden of PSSSA based on SRSMCS is extremely heavy, especially when applied in a largescale power system.
In PSSSA problems, the nodal injection power distributions of generators and loads are the input random variables. The output random variables are system eigenvalues and damping ratios. When performing SRSMCS based PSSSA, the following several assumptions are default in the case study:

1) All PV generator nodal output active powers are input random variables with independent normal distributions and the voltage amplitudes are considered to be constant at the nominal values.

2) All PQ loads are input random variables with normal distribution and the load power factors are considered to be constant at the nominal values.

3) All PQ loads are converted to constant PQ load model in the linearization procedure.

4) The correlation between load demands is described by a correlation coefficient matrix.

5) Generator and branch fault are not considered.
With these simplified conditions, the probabilistic analysis problem of small signal stability can be better expressed in the presence of system uncertainties and load correlation.
3. Proposed Method
In this section, a LHS technique is proposed first. Then how to deal with the correlation between input random variables is described in detail. The computational procedure of LHSMCS based PSSSA method is also given.
 3.1 LHS technique
LHS is a stratified sampling technique, proposed by Mckay et al in 1979
[21]
. This technique has been applied in many engineering areas. Compared with the traditional SRS, LHS has the following advantages
[14
,
21]
:

1) Under the same sample size, LHS covers a larger random variable sampling space and has better convergence property;

2) LHS is a sampling technique with more robustness than SRS.
LHS consists of two major steps, namely, sampling and permutation. The sampling’s purpose is to produce a series of presentation samples, which describe the distribution of each input random variable, and to ensure that the distribution area can be covered by sampling points completely. Permutation is aimed at making the correlation of sampled values minimized by changing the arrange order of sampling values of input random variables.
In the sampling process, a sampling technology, named lattice sampling, is adopted in this paper. For
M
input random variables
X_{1}
,
X_{2}
, …
X_{M}
in a probabilistic problem, the probability cumulative function of
X_{m}
is expressed as
Y_{m}
=
F_{m}
(
X_{m}
). For the sample size
N
, the range of
Y_{m}
is divided into
N
equal interval areas, and the width of each interval area is 1/
N
. Every sampling value is chosen from the midpoint of each interval area. Then the sampling values of
X_{m}
are calculated by the inverse function of
where
N
is the number of maximum sampling.
A row of sampled values
X
_{m1}
,
X
_{m2}
, …
X_{mN}
represent presentation samples. Once the
M
input random variables are sampled, a
M
×
N
primary sampling matrix X can be obtained.
In the random permutation process, the sampled value of every input variable is chosen randomly from the corresponding row in matrix X. Under this circumstance, the undesired correlations between different input random variables may affect the accuracy of the solutions. In order to reduce the undesired correlation of sampling values of multiple input random variables, various permutation skills, such as Cholesky decomposition
[14]
and Nataf transformation
[15
,
22]
, have been proposed in LHS. The Cholesky decomposition is employed in this paper to eliminate the correlation among rows of the sampling matrix X since the calculation burden of the Cholesky decomposition method is small when handling the problems with large numbers of input random variables
[14]
.
In the Cholesky decomposition based permutation process, after the primary sampling matrix X is obtained, a
M
×
N
ordering matrix
L
_{0}
is generated to indicate the order number for permutation of X, then the samples in every row of X are replaced with the order according to the order numbers in the corresponding row of ordering matrix
L
_{0}
.
The original ordering matrix
L
_{0}
is generated randomly, and X arranged according to this random generated
L
_{0}
is not completely independent in each row. The associated
M
×
M
correlation coefficient matrix of original ordering matrix
L
_{0}
can be represented by
C_{X}
, which is a positive definite and symmetric matrix and can be decomposed by Cholesky decomposition into
where
G
is a lower triangular matrix. Then a new
M
×
N
matrix
F
, whose correlation coefficient matrix is an identity matrix, can be obtained by the following
Generally, the elements in matrix
F
are not necessarily integer or positive number, and cannot be used to indicate orders in the sampling matrix. So the rows of new ordering matrix
L
_{1}
are updated to the orders of the element in the corresponding row of
F
. In the Cholesky decomposition based permutation process, sampling matrix X is arranged according to this updated
L
_{1}
instead of randomly generated
L
_{0}
in random permutation. Since
F
has an identity correlation matrix, the updated order matrix
L
_{1}
, which indicates the orders of the element in a row of
F
, has a much smaller correlation than
L
_{0}
. And each row in the sampling matrix X is considered to approximate independence.
Assume
C_{X′}
is the correlation coefficient matrix of arranged matrix X. It should be noted that after the permutation process, the correlation coefficient matrix
C_{X′}
is not accurately equal to
C_{X}
. There are small differences between correlation coefficient matrix
C_{X′}
and
C_{X}
, because
C_{X}
is just the order correlation matrix of X.
Through the above two steps, LHSbased MCS generates more uniform, independent random sample values and thus improves the convergence performance, reduces the sampling size and improves the calculation speed with the same calculation accuracy.
 3.2 Processing load correlation
The LHS analyzed above can make the sampling values of input random variables to be independent. Although the sampling values can’t be transformed to be completely independent, the correlation between them is quite small, and sampled values can be regarded as approximate independent variables.
In order to take load correlation into account, we need to generate the normal distribution samples of random variables with a correlation coefficient matrix
C_{Z}
. The basic idea is similar to the inverse process of removing correlation in the permutation process of LHS.
Let Y=[
y
_{1}
,
y
_{2}
,…,
y
_{m}
] be the
m
independent standard normal variables, a lower triangular matrix
E
is obtained by Cholesky decomposition of
C_{Z}
. Then
Z^{′}
are the standard normal random variables with correlation matrix
C_{Z}
.
According to the definition of the correlation coefficient matrix, the correlation coefficient matrix is the same as the covariance matrix. Then,
Z^{′}
are the standard normal random variables with correlation coefficient matrix
C_{Z}^{′}
.
For nonstandard normal distribution
P
, the desired correlations between the samples of different random variables can also be induced by the above Cholesky decomposition method. Detailed process is referred to
[15]
.
It’s worth to note that the permutation of the LHS samplings according to the Cholesky Decomposition can not obtain the exact desired correlation coefficient matrix between the random input variables. But the final achieved correlation coefficient matrix
C_{Z}^{′}
is very close to the desired correlation coefficient matrix
C_{Z}
. Assume ∆
ρ_{rms}
denotes the root mean square correlation error index of
C_{Z}^{′}
and
C_{Z}
, and its squared value can be expressed in (15) as follows to evaluate the error degree of correlation.
where
ρ_{kj}
and
ρ_{kj}^{′}
are the offdiagonal elements in correlation coefficient matrix
C_{Z}
and
C_{Z}^{′}
, respectively.
M
is the number of input random variables.
 3.3 Procedure of LHSMCS method
The procedure of solving PSSSA problems with LHSMCS is as follows.

1) Read the study power system basic data for PSSSA, the sample sizeN, the numberKof PV generators, the numberMof loads, the PDF of input random variables and the desired correlation coefficient matrixCZ. In this procedure, only the uncertainties of PV generators and nodal loads are considered, thus the number of input random variables is equal to the number of PV generatorsKplus loadsM.

2) Generate two normal distribution primary sampling matrix X and Y with dimensionM×NandK×Nby LHS. The same process of X and Y will be omitted for space saving.

3) Randomly generate anM×Noriginal ordering matrixL0of primary sampling matrix X to indicate rank number for permutation.

4) Calculate the undesired correlation coefficient matrix CXof random sampling matrix X with ordering matrixL0.

5) Decompose CXand CZby the Cholesky decomposition into CX=G·GTand CZ=E · ET.

6) Compute the new ordering matrix ofF=G−1L0.

7) Form the new independent random sampling matrix X' with new ordering matrixF.

8) Process correlation procedures. Z′=EX′is the correlated sampling matrix with correlation coefficient matrix CZ.

9) Setn= 1.

10) Perform deterministic SSSA with the ith column of the correlated sampling matrix Z′to obtain the eigenvalues and damping ratios.

11) Setn=n+ 1 ; ifn＜N, go to Step 10 to perform deterministic SSSA with the next column of the sampling matrix. Otherwise, proceed to Step 12.

12) Identify the electromechanical oscillation mode from all of the eigenvalues.

13) Analyze the statistical features of the critical oscillation mode and calculate its PIstI or PIseI.
4. Case Study
In this section, the proposed method is applied to two standard test systems. The twoarea, fourmachine system is used to confirm the efficiency of the recommended LHSbased method. The 16machine, 68bus test system is used to investigate the influence of load correlation on the probabilistic insecurity index of the critical oscillation mode and check the robustness of the proposed procedure in identifying critical oscillation modes. The error degree of correlation between multiple input random variables is also evaluated. The programs are developed in Matlab2012a, and PSAT
[23]
is adopted to solve the deterministic SSSA. This work is implemented on a PC with Intel Core2 2.93GHz CPU and 3GB of RAM.
 4.1 The twoarea, Fourmachine system case
The twoarea, fourmachine system consists of four machines and two large loads located at node 7 and 9, as shown in
Fig.1
. Each machine has been represented by a sixthorder detailed dynamic model. All four machines are equipped with the IEEE type ST1 excitation control. With and without power system stabilizers (PSS) are considered in the case simulations. System loads are modeled as constant PQ load models. There are total 5 input random variables in this system, 3 for PV generators and 2 for correlated loads. Detailed system parameters refer to
[24]
.
Twoarea, Fourmachine test system.
SRSMCS with 20,000 samples is used as a reference to evaluate the accuracy of the proposed method. The system has three electromechanical modes, one interarea mode and two local modes. The interarea mode is unstable when no PSS is equipped at the nominal system condition. The mean values and standard deviations based on SRSMC without and with PSS are shown in
Tables 1
and
2
.
Electromechanical oscillation modes with SRSMCS accurate calculation (without PSS)
Electromechanical oscillation modes with SRSMCS accurate calculation (without PSS)
Electromechanical oscillation modes with SRSMCS accurate calculation (with PSS)
Electromechanical oscillation modes with SRSMCS accurate calculation (with PSS)
In order to verify the effectiveness of the presented method comprehensively, the relative error indexes (REIs) (
ε_{μ}
and
ε_{σ}
) of the output variables’ mean values and standard deviations are used to assess the calculation accuracy from statistical aspects.
where
μ_{a}
and
σ_{a}
are the exact mean values and standard deviations of the output random variables obtained by SRSMCS, while
μ_{L}
and
σ_{L}
are those obtained by LHSMCS with sample size N.
In this simulation, the sample size changes from 100 to 900 with a step size of 200. In each sample size, the REIs of the three oscillation modes with PSS are listed in
Table 3
. In order to avoid the influence of the numerical calculation uncertainty, the listed values are the average of 50 times.
REIs of there electromechanical oscillation modes with LHSMCS method
REIs of there electromechanical oscillation modes with LHSMCS method
From the above table, the average REIs of all the three oscillation modes are very small and they are decreasing with the increase of sample size. Moreover, the REIs with sample size 500 achieve the required calculation accuracy. For example, the relative error of the mean value is 6.35 × 10
^{−3}
, while the standard deviation is 4.11 × 10
^{−2}
, corresponding to oscillation mode 1. That means the calculated values of the two methods are very close.
The PDF curve of interarea oscillation mode’s real part is shown in
Fig. 2
corresponding to LHS with N=500 and SRS with N=20,000. It can be seen that the two curves are very proximity.
Comparison between the PDF of interarea mode’s real part by LHS with N=500 and SRS with N= 20,000.
The REIs in
Table 2
and the PDF curve in
Fig. 2
have verified that the LHSMCS method has high calculation accuracy with a small sample size, and the sample size of 500 can meet a high accuracy.
In order to compare the convergence, for each sample size, the LHSMCS and SRSMCS methods run 50 times, respectively. The REIs of all oscillation modes are calculated, where the mean and standard deviation of REIs corresponding to oscillation mode 1 are given in
Table 3
. The REIs of other oscillation modes have the similar results.
The simulation results in
Table 4
indicate that the average mean value and standard deviation of REI, obtained by LHSMCS method, are much smaller than those obtained by SRSMCS, especially in small sample sizes. For instance, the mean value and standard deviation of
obtained by LHSMCS are 6.35×10
^{−3}
and 4.11×10
^{−3}
, while those obtained by SRSMCS are 11.7×10
^{−3}
and 9.33×10
^{−3}
. The smaller values of the average mean value and standard deviation prove that the LHS is more robust and converges much faster than SRS.
REIs comparisons of two methods with 50 times
REIs comparisons of two methods with 50 times
Moreover, the computing time for the LHSMCS method with sample size 500 is 78.2s while SRSMCS method is 76.6s. The process of sampling and permutation increases a little time. When the sample size and the number of input variables become larger, the time difference is not significant. However, to achieve the same calculation accuracy as LHSMCS method, the SRSMCS takes much longer time. That means with the same small sample size, LHSMCS is a more effective method compared with LHSMCS.
From all of the simulation results, we can conclude that LHSMCS is very efficient in PSSSA problems under small sample size.
 4.2 16machine, 68bus test system case
The effectiveness of the LHSMCS is discussed in the previous subsection with a simple power system. In this subsection, the influences of load correlation on probabilistic small signal stability are studied through the 16machine, 68bus test system which contains multiple correlated loads. This system is the reducedorder equivalent model of New England test system (NYTS ) and New York power system (NYPS) interconnected networks
[19]
. The system structure is shown in
Fig. 3
, there are five geographical regions with three major transmission corridors between NETS and NYPS in the system. In this work, the tieline power exchange between NETS and NYPS is 700MW in total. Each machine has been represented by a sixthorder detailed dynamic model. All sixteen machines are equipped with the IEEE type ST1 excitation control. System loads are modeled as constant PQ load models. There are a total of 50 uncertain parameters as input random variables within the test system being investigated (15 PV generators are considered to be completely independent and 35 loads are considered to be correlative).
16machine, 68bus test system
Before assessing the impact of multiple correlated input random variables on probabilistic index of power system, we first verify whether the root mean square correlation error index (∆
ρ_{rms}
in Eq. (15)) of multiple input random variables can be achieved by the proposed method.
The values of ∆
ρ_{rms}
of simple random sampling with random permutation (SRSRP), LHS with random permutation and (LHSRP) LHS with Cholesky decomposition based permutation (LHSCD) of different sampling dimensions are shown in
Tables 5
and
6
. The results in
Table 5
demonstrate the root mean square correlation error index ∆
ρ_{rms}
of 15 generators power input random variables. It is clear that LHSCD has a smaller ∆
ρ_{rms}
value compared with SRSRP and LHSRP when the input variable sampling values are independent. The correlation error index ∆
ρ_{rms}
of 35 system loads random variables are listed in
Table 6
. The results demonstrate that although LHSCD can’t obtain the exact desired correlation coefficients between multiple correlated random variables, the correlation error index ∆
ρ_{rms}
is very small and its maximum value produced by LHSCD is 1.45% (with N=500). Furthermore, the correlation error index ∆
ρ_{rms}
is decreasing with the increase of sampling samples size N. So the obtained samples can reflect the expected correlation between random variables.
Comparisons of ∆ρrmsof different methods in 16machine, 68bus Test System with 15 independent input random variables
Comparisons of ∆ρ_{rms} of different methods in 16machine, 68bus Test System with 15 independent input random variables
Comparisons of ∆ρrmsof different correlation coefficients in 16machine, 68bus test system with 35 correlated input random variables by LHSCD
Comparisons of ∆ρ_{rms} of different correlation coefficients in 16machine, 68bus test system with 35 correlated input random variables by LHSCD
This test system displays four interarea modes, as shown in
Table 7
. The system is unstable when no PSS is equipped at the nominal system condition. They are 6 unstable local oscillation modes (damping ratios are shown in italics in
Table 7
). So the probabilistic analysis of small signal stability without PSS is not considered in this study. Eigenvalue analysis shows that all local modes and interarea modes receive acceptable damping (5%) when some local PSSs are installed in the system
[19]
. However, for some extreme operating conditions, two interarea modes might be poorly damped (M1 and M2 Bold in the table). M1 with frequency at about 0.38 Hz dominates the power oscillation between NYTS, NYPS, and the rest of the system, while M2 oscillating at a frequency of about 0.70 Hz mainly reflects the power oscillation between NYTS and NYPS. Hence, the probabilistic analysis presented here focuses on only the lowest likely poorly damped low frequency interarea oscillation mode M1, which is the critical interarea oscillation mode of the system.
All of the oscillation modes in 16machine, 68bus test system with and without PSS at the nominal system condition
All of the oscillation modes in 16machine, 68bus test system with and without PSS at the nominal system condition
A correlation coefficient matrix is used to represent load correlation. When the correlation coefficient varies from 0 to 1.0 with a step of 0.1, the probabilistic insecurity index (PIseI) of the critical interarea oscillation mode is shown in
Table 8
at different levels of load uncertainty. The mean values and standard deviations of the critical interarea oscillation mode are also given in
Tables 9
and
10
. The results in
Table 10
are also illustrated in
Fig. 4
. LHSMCS method is adopted to obtain all of the test results.
PIseI of the critical oscillation mode M1 (N=500)
PIseI of the critical oscillation mode M1 (N=500)
The mean values of the critical oscillation mode M1 (N=500)
The mean values of the critical oscillation mode M1 (N=500)
The standard deviations of the critical oscillation mode M1 (N=500)
The standard deviations of the critical oscillation mode M1 (N=500)
The standard deviations of the critical mode
From the simulation results listed in
Table 8
to
10
and
Fig. 4
, it is observed that:

1) The system insecurity probability PIseI (ζi< 5% ) is higher with the increasing of uncertainty level and

correlation coefficient.

2) Under a certain load uncertainty level, the PIseI of critical oscillation mode is enlarged with the increase of load correlation.

3) Uncertainty and correlation interaction effect is nonlinear.

4) Load correlation has no significant influence on the mean value of critical oscillation mode, it slightly reduces the damping ratio and frequency of the critical oscillation mode. But the standard deviation of critical oscillation mode obviously increases with enhancement of load correlation. It is noted that the relationship between standard deviation and load correlation coefficient is almost linear as shown inFig. 4. When the loads are highly interrelated, the standard deviations own the amplification properties. That means the instability risk increase when the load correlation is strongly correlated.
The PDF curves of the critical interarea oscillation mode M1’s real part and damping ratio with different load correlation coefficient is shown in
Figs. 5
and
6
.
The PDF of critical mode’s real part
The PDF of critical mode’s damping ratio
From these two PDF curves, it is observed that different strength of load correlation is corresponding to different probability density. The highandnarrow type PDF curve, which is corresponding to low load correlation, has demonstrated the smaller volatility of the critical mode’s real part and damping ratio, while the pyknic type PDF curve has reflected much more fluctuation of the critical mode under strong load correlation. The curves in
Fig. 6
also have proven that the PIseI of the critical mode is increasing with the strengthening of load correlation. It is noted that all other modes are with similar characteristics.
Load correlation reflects the simultaneity of power demand and it also means the reduction of load random characteristics. In fact, when load correlation is more and more strong, the consistency of the load change is more obvious and the probability of power system appearing heavy load increases. When the system becomes heavy load, the PIseI even PIstI corresponding to the critical oscillation mode becomes larger. This study confirmed this viewpoint.
5. Conclusion and Future Work
This paper has proposed a method by using LHS technique as a kind of variance reduction technique to solve the PSSSA problems in the presence of load correlation. The test results have illustrated that LHSMCS can achieve a better sampling efficiency compared with SRSMCS, and provide an alternative PSSSA method with much smaller size. The simulation results also show that, the correlations between loads have a significant effect on PSSSA results. When load correlation increases, the instability probability of critical oscillation mode increases accordingly. The PSSSA study considering load correlation factors is more realistic in power system operation.
In the future research, we will focus on applying this method to a largescale practical power system. In addition, Wind speed correlation
[22
,
25]
and its effect on power system probabilistic small signal stability will be taken into account.
BIO
Jian Zuo He received the B.Eng. degree in electrical engineering from Huazhong University of Science and Technology (HUST), Wuhan, China in 2008, where he is currently pursuing the Ph.D. degree. His research interests include power system dynamic stability analysis and wide area damping control.
Yinhong Li She received the B.Eng. degree and the Ph.D. degree in electrical engineering from Huazhong University of Science and Technology (HUST), Wuhan, China in 1998 and 2004, respectively. Dr. Li is an Associate Professor in the College of Electrical and Electronic Engineering at HUST. Her research interests include relay coordination, voltage stability, and load modeling.
Defu Cai He received the B.Eng. degree and the Ph.D. degree in electrical engineering from Huazhong University of Science and Technology(HUST), Wuhan, China in 2009 and 2014, respectively. He now works at State Grid Hubei Electric Power Research Institute. His research interests include renewable energy, microgrid, electric vehicle and power system analysis.
Dongyuan Shi He received the B. Eng. and the Ph.D. degrees in electrical engineering from Huazhong University of Science and Technology(HUST), Wuhan, China, in 1996 and 2002, respectively. He is currently an Associate Professor in the College of Electrical and Electronic Engineering at HUST. His research interests include power system analysis and application of artificial intelligence in power systems.
Rueda J. L.
,
Guaman W. H.
,
Cepeda J. C.
,
Erlich I.
,
Vargas A.
2013
“Hybrid Approach for Power System Operational Planning With Smart Grid and Small Signal Stability Enhancement Considerations”
IEEE Trans. Smart Grid
4
(1)
530 
539
DOI : 10.1109/TSG.2012.2222678
Burchett R.C.
,
Heydt G.T.
1978
“Probabilistic Methods For Power System Dynamic Stability Studies”
IEEE Trans. Power App. Syst.
PAS97
(3)
695 
702
DOI : 10.1109/TPAS.1978.354540
Yi H.Q.
,
Hou Y.H.
,
Cheng S.J.
,
Zhou H.
,
Chen G.G.
2007
“Power system probabilistic small signal stability analysis using two point estimation method”
Proceeding of 42nd International universities power engineering conference
New York
Dong Z.
,
Mishra Y.
,
Zhang P.
2009
“Power system sensitivity analysis considering induction motor loads”
Proceeding of 8th International Conference on Advances in Power System Control, Operation and Management (APSCOM 2009)
Hong Kong
Wang K.W.
,
Chung C.Y.
,
Tse C.T.
,
Tsang K.M.
2000
“Improved probabilistic method for power system dynamic stability studies”
IEE ProceedingsGeneration, Transmission and Distribution
147
(1)
37 
43
DOI : 10.1049/ipgtd:20000025
Li M. Y.
,
Ma J.
,
Dong Z. Y.
2012
“Uncertainty Analysis of Load Models in Small Signal Stability”
Proceedings of the IEEE Power and Energy Society General Meeting
San Diego, CA
Preece R.
,
Milanovic J. V.
2012
“The Probabilistic Collocation Method for dealing with uncertainties in power system small disturbance studies”
Proceeding of the IEEE Power and Energy Society General Meeting
San Diego, CA
Rueda J. L.
,
Colome D. G.
,
Erlich I.
2009
“Assessment and Enhancement of Small Signal Stability Considering Uncertainties”
IEEE Trans. Power Systems
24
(1)
198 
207
DOI : 10.1109/TPWRS.2008.2009428
Huazhang H.
,
Chung C. Y.
2012
“Coordinated Damping Control Design for DFIGBased Wind Generation Considering Power Output Variation”
IEEE Trans. Power Systems
27
(4)
1916 
1925
DOI : 10.1109/TPWRS.2012.2190110
Xu Z.
,
Dong Z. Y.
,
Zhang P.
2005
“Probabilistic small signal analysis using Monte Carlo simulation”
Proceeding of the IEEE Power and Energy Society General Meeting
New York
Dissanayaka A.
,
Annakkage U. D.
,
Jayasekara B.
,
Bagen B.
2011
“RiskBased Dynamic Security Assessment”
IEEE Trans. Power Systems
26
(3)
1302 
1308
DOI : 10.1109/TPWRS.2010.2089809
Billinton R.
,
Dange H.
2008
“Effects of Load Forecast Uncertainty on Bulk Electric System Reliability Evaluation”
IEEE Trans. Power Systems
23
(2)
418 
425
DOI : 10.1109/TPWRS.2008.920078
Li W.
,
Billinton R.
1991
“Effect of bus load uncertainty and correlation in composite system adequacy evaluation”
IEEE Trans. Power Systems
6
(4)
1522 
1529
DOI : 10.1109/59.116999
Yu H.
,
Chung C. Y.
,
Wong K. P.
,
Lee H. W.
,
Zhang J. H.
2009
“Probabilistic Load Flow Evaluation with Hybrid Latin Hypercube Sampling and Cholesky Decomposition”
IEEE Trans. Power Systems
24
(2)
661 
667
DOI : 10.1109/TPWRS.2009.2016589
Yan C.
,
Jinyu W.
,
Shijie C.
2013
“Probabilistic Load Flow Method Based on Nataf Transformation and Latin Hypercube Sampling”
IEEE Trans. Sustainable Energy
4
(2)
294 
301
DOI : 10.1109/TSTE.2012.2222680
Zhen S.
,
Jirutitijaroen P.
2011
“Latin Hypercube Sampling Techniques for Power Systems Reliability Analysis with Renewable Energy Sources”
IEEE Trans. Power Systems
26
(4)
2066 
2073
Jirutitijaroen P.
,
Singh C.
2008
“Comparison of Simulation Methods for Power System Reliability Indexes and Their Distributions”
IEEE Trans. Power Systems
23
(2)
486 
493
DOI : 10.1109/TPWRS.2008.919425
Jirutitijaroen P.
,
Singh C.
2008
“Reliability Constrained MultiArea Adequacy Planning Using Stochastic Programming With Sample Average Approximations”
IEEE Trans. Power Systems
23
(2)
504 
513
Rogers G.
2000
Power System Oscillations
Springer US
Boston, London
Cai D.
,
Chen J.
,
Shi D.
,
Duan X.
,
Li H.
,
Yao M.
2012
“Enhancements to the Cumulant Method for probabilistic load flow studies”
Proceeding of Power and Energy Society General Meeting
San Diego, CA
McKay M. D.
,
Beckman R. J.
,
Conover W. J.
1979
“Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code”
Technometrics
21
(2)
239 
245
Li H.
,
Cai D.
,
Li Y.
2013
“Probabilistic Assessment of Voltage Stability Margin in Presence of Wind Speed Correlation”
Journal of Electrical Engineering & Technology
8
(4)
719 
728
DOI : 10.5370/JEET.2013.8.4.719
Kundur P.
1994
Power system stability and control
McGrawhill
New York
EIKeib A. A.
2009
“Probabilistic Reliability Evaluation of Power Systems Including Wind Turbine Generators Considering Wind Speed Correlation”
Journal of Electrical Engineering & Technology
4
(4)
485 
491
DOI : 10.5370/JEET.2009.4.4.485