# 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 MCS-based 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 large-scale 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 LHS-based 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 LHS-based 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

Multi-machine power system dynamic behavior is represented by a set of nonlinear differential-algebraic-equations (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 (x0, y0), the full system linearization algebraic-differentiation- equations (2) can be obtained by Taylor series expansion by ignoring the higher-order terms.

where A = ▽xf, B = ▽yf, C = ▽xg, and D = ▽yg where. Then, the system state Eq. (3) is obtained by eliminating the algebraic variable Δy

where As = A−BD −1C is the system state matrix, thus implicitly assuming that D is nonsingular.

The eigenvalues λ = [λ1, λ2, …, λn] of the system state matrix As are calculated by solving the Eq. (4).

The eigenvalues may be real or complex. If As 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 large-scale 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 Pe and the system PQ load powers PL, QL 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 pi which belongs to generator or load power vector P=[p1, p2,…, pn]T, the PDF f(pi) can be expressed as:

where μi and σi2 are mean and variance of nodal power pi 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 CZ is the correlation coefficient matrix of input random variable vector Z=[z1, z2, … , zn]T. The correlation coefficient ρwij can be calculated by (10).

where σzi and σzj are the standard deviations of zi and zj, Cov(zi, zj) is the covariance of zi and zj, and ρzy is the correlation coefficient of wi and wj.

For a series of ideal independent random variables, the correlation coefficient matrix is a standard unit matrix. But in most situations, CZ 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 CZ is positive definite.

## 2.5. SRS-MCS based PSSSA

Simple random sampling Monte Carlo simulation (SRS-MCS) 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 SRS-MCS can obtain high accuracy results, this method is timeconsuming. The computational burden of PSSSA based on SRS-MCS is extremely heavy, especially when applied in a large-scale 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 LHS-MCS 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 X1, X2, … XM in a probabilistic problem, the probability cumulative function of Xm is expressed as Ym = Fm(Xm). For the sample size N, the range of Ymis 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 Xm are calculated by the inverse function of

where N is the number of maximum sampling.

A row of sampled values Xm1, Xm2, … XmN 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 L0 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 L0.

The original ordering matrix L0 is generated randomly, and X arranged according to this random generated L0 is not completely independent in each row. The associated M × M correlation coefficient matrix of original ordering matrix L0 can be represented by CX , 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 L1 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 L1 instead of randomly generated L0 in random permutation. Since F has an identity correlation matrix, the updated order matrix L1, which indicates the orders of the element in a row of F, has a much smaller correlation than L0. And each row in the sampling matrix X is considered to approximate independence.

Assume CX′ is the correlation coefficient matrix of arranged matrix X. It should be noted that after the permutation process, the correlation coefficient matrix CX′ is not accurately equal to CX. There are small differences between correlation coefficient matrix CX′ and CX, because CX is just the order correlation matrix of X.

Through the above two steps, LHS-based 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 CZ. The basic idea is similar to the inverse process of removing correlation in the permutation process of LHS.

Let Y=[y1, y2,…, ym] be the m independent standard normal variables, a lower triangular matrix E is obtained by Cholesky decomposition of CZ. Then Z′ are the standard normal random variables with correlation matrix CZ.

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 CZ′.

For non-standard 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 CZ′ is very close to the desired correlation coefficient matrix CZ. Assume ∆ρrms denotes the root mean square correlation error index of CZ′and CZ, and its squared value can be expressed in (15) as follows to evaluate the error degree of correlation.

where ρkj and ρkj′ are the off-diagonal elements in correlation coefficient matrix CZ and CZ′, respectively. M is the number of input random variables.

## 3.3 Procedure of LHS-MCS method

The procedure of solving PSSSA problems with LHSMCS is as follows.

1) Read the study power system basic data for PSSSA, the sample size N, the number K of PV generators, the number M of loads, the PDF of input random variables and the desired correlation coefficient matrix CZ. 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 generators K plus loads M. 2) Generate two normal distribution primary sampling matrix X and Y with dimension M × N and K × N by LHS. The same process of X and Y will be omitted for space saving. 3) Randomly generate an M × N original ordering matrix L0 of primary sampling matrix X to indicate rank number for permutation. 4) Calculate the undesired correlation coefficient matrix CX of random sampling matrix X with ordering matrix L0. 5) Decompose CX and CZ by the Cholesky decomposition into CX = G·GT and CZ = E · ET. 6) Compute the new ordering matrix of F = G−1L0. 7) Form the new independent random sampling matrix X' with new ordering matrix F. 8) Process correlation procedures. Z′ = EX′ is the correlated sampling matrix with correlation coefficient matrix CZ. 9) Set n = 1. 10) Perform deterministic SSSA with the ith column of the correlated sampling matrix Z′ to obtain the eigen-values and damping ratios. 11) Set n = n + 1 ; if n ＜ 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 two-area, four-machine system is used to confirm the efficiency of the recommended LHS-based method. The 16-machine, 68-bus 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.93-GHz CPU and 3GB of RAM.

## 4.1 The two-area, Four-machine system case

The two-area, four-machine 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 sixth-order detailed dynamic model. All four machines are equipped with the IEEE type ST-1 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].

**Fig. 1.**Two-area, Four-machine test system.

SRS-MCS with 20,000 samples is used as a reference to evaluate the accuracy of the proposed method. The system has three electromechanical modes, one inter-area mode and two local modes. The inter-area mode is unstable when no PSS is equipped at the nominal system condition. The mean values and standard deviations based on SRS-MC without and with PSS are shown in Tables 1 and 2.

**Table 1.**Electromechanical oscillation modes with SRSMCS accurate calculation (without PSS)

**Table 2.**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 SRS-MCS, while μLand σ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.

**Table 3.**REIs of there electromechanical oscillation modes with LHS-MCS 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 inter-area 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.

**Fig. 2.**Comparison between the PDF of inter-area 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 LHS-MCS 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 LHS-MCS and SRS-MCS 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 LHS-MCS method, are much smaller than those obtained by SRS-MCS, especially in small sample sizes. For instance, the mean value and standard deviation of obtained by LHS-MCS are 6.35×10−3 and 4.11×10−3 , while those obtained by SRS-MCS 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.

**Table 4.**REIs comparisons of two methods with 50 times

Moreover, the computing time for the LHS-MCS method with sample size 500 is 78.2s while SRS-MCS 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 LHS-MCS method, the SRS-MCS takes much longer time. That means with the same small sample size, LHS-MCS is a more effective method compared with LHS-MCS.

From all of the simulation results, we can conclude that LHS-MCS is very efficient in PSSSA problems under small sample size.

## 4.2 16-machine, 68-bus test system case

The effectiveness of the LHS-MCS 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 16-machine, 68-bus test system which contains multiple correlated loads. This system is the reduced-order 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 tie-line power exchange between NETS and NYPS is 700MW in total. Each machine has been represented by a sixth-order detailed dynamic model. All sixteen machines are equipped with the IEEE type ST-1 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).

**Fig. 3.**16-machine, 68-bus 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 (SRS-RP), LHS with random permutation and (LHS-RP) LHS with Cholesky decomposition based permutation (LHS-CD) 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 LHS-CD has a smaller ∆ρrms value compared with SRS-RP and LHS-RP 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 LHS-CD 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.

**Table 5.**Comparisons of ∆ρrms of different methods in 16-machine, 68-bus Test System with 15 independent input random variables

**Table 6.**Comparisons of ∆ρrms of different correlation coefficients in 16-machine, 68-bus test system with 35 correlated input random variables by LHS-CD

This test system displays four inter-area 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 inter-area 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 inter-area oscillation mode M1, which is the critical inter-area oscillation mode of the system.

**Table 7**All of the oscillation modes in 16-machine, 68-bus 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 inter-area oscillation mode is shown in Table 8 at different levels of load uncertainty. The mean values and standard deviations of the critical inter-area oscillation mode are also given in Tables 9 and 10. The results in Table 10 are also illustrated in Fig. 4. LHS-MCS method is adopted to obtain all of the test results.

**Table 8.**PIseI of the critical oscillation mode M1 (N=500)

**Table 9.**The mean values of the critical oscillation mode M1 (N=500)

**Table 10.**The standard deviations of the critical oscillation mode M1 (N=500)

**Fig. 4.**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 in Fig. 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 inter-area oscillation mode M1’s real part and damping ratio with different load correlation coefficient is shown in Figs. 5 and 6.

**Fig. 5**The PDF of critical mode’s real part

**Fig. 6.**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 high-and-narrow 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 LHS-MCS can achieve a better sampling efficiency compared with SRS-MCS, 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 large-scale 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.