A Probabilistic Assessment Method for Voltage Stability Considering Large Scale Correlated Stochastic Variables

Voltage stability has always been one of the most important concerns. As the increasing integration of large-scale renewable energy sources in power systems, the correlation between load demands and renewable energy systems becomes more and more complex and important for probabilistic voltage stability. There are two significant issues for probabilistic voltage stability assessment: (i) how to choose the reasonable power increment direction which determines the reliability of voltage stability assessment when considering the actual operating characteristics of the power system; and (ii) how to obtain the samples characterized with the specified distribution and the desired correlation. We propose methodologies to define the reasonable power increment direction with theoretical proof. Moreover, power method transformation combined with Latin hypercube sampling and twice-permutation technique is proposed for probabilistic voltage stability assessment. Case studies with two modified IEEE test systems show that the proposed method is accurate and efficient.


I. INTRODUCTION
Voltage stability has always been one of the most important concerns in power system planning and operations [1].Over the past decades, deterministic voltage stability was studied in many previous literatures [2]- [4].Load margin is an efficient voltage stability index which is convenient to be understood and easy to be used for operators [5], [6].For the deterministic voltage stability evaluation, the smallest load margin associated with the closest bifurcation point is always applied [7], [8].However with the integrating of large-scale renewable energy sources in power systems and load variations, there are more and more stochastic factors The associate editor coordinating the review of this manuscript and approving it for publication was Arash Asrari .
affecting power system planning and operations.The effects of these uncertainties need to be thoroughly explored in order to fully ensure system security.Therefore, the probabilistic voltage stability analysis was studied and the results provided rich information to discover more useful insights than the deterministic analysis did [9], [10].For the probabilistic voltage stability assessment, the covariance matrix between stochastic variables can be used to measure the most probable bifurcation point with an iterative method [11].Authors in [12] evaluated voltage stability with global sensitivity analysis, in which the load increment direction is assumed as normal distribution.Authors in [13] calculated the voltage collapse point by increasing only one single bus load and remaining the other constant.Authors in [14] adopted nonconforming load model as the load direction.The load VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see http://creativecommons.org/licenses/by/4.0/model is assumed as normal distribution to estimate voltage stability [15].Authors in [16] estimated voltage stability with a direction by predicting load-path.However, two issues are still not solved well based on the previous literature.Firstly, the correlation between the stochastic variables are ignored, which are more and more complex and important for the probabilistic voltage stability analysis.Secondly, the probabilistic voltage stability is studied without fully explaining the actual power increment direction.The stochastic variables and the correlation between them are not carefully considered when choosing the power increment direction.Under such circumstances, the load margin assessment may not truly reflect the actual stochastic characteristics of the system and therefore cannot be effectively used to evaluate probabilistic voltage stability.
In order to deal with these two problems, based on a reasonable power increment direction for generators and loads considering the stochastic variables and the correlation between them, we propose the power method combined with Latin hypercube sampling and twice-permutation technique (PLT) in this paper.The power increment direction is proven in theory that the actual stochastic characteristics of the power system are maintained.The PLT realizes the nonlinear correlation matrix transformation between the real stochastic variables and the standard normal variables, and obtains the samples of large scale input stochastic variables with desired correlation coefficients.
The rest of this paper is organized as follows.Models of typical stochastic sources and assessment of load margin are introduced in Section II.The method used to generate the samples of the stochastic variables with desired correlation is proposed in Section III.The performances of the proposed method are studied with a modified IEEE 14-bus system and a modified IEEE 118-bus system in Section IV.Finally, Section V concludes the paper.

II. PROBABILISTIC VOLTAGE STABILITY MODEL A. STOCHASTIC SOURCES
There are more and more stochastic sources in power systems than before.Modeling of the stochastic characteristic and the correlation between the sources is essential for probabilistic voltage stability assessment.Generally, the wind generation, the solar generation and the load are described as PQ buses and the power factor is fixed.The active power of the wind generation, P W , is decided by the wind speed modeled as the Weibull distribution, the active power of the photovoltaic generation, P PV , is typically modeled as Beta distribution, and the active power of the load, P L , is commonly modeled as normal distribution [17], [18].The correlation matrix, R, is symmetric.

B. VOLTAGE STABILITY MODEL
The index of load margin is one kind of overall voltage stability indices (VSIs), and in general the accuracy of the overall voltage stability indices (VSIs) is better than that of the line and bus VSIs [19].The max load margin is always adopted to evaluate power system voltage stability as it is an efficient index, which is convenient to be understood and easy to be used for operators [20].Optimization based methods are popular to calculate the max load margin [21].In this paper, the model based on optimal power flow is adopted to calculate the max load margin.The optimization model is as follows, which is a nonlinear optimal problem and the interior point method is used for solving it.max λ (1.1) 3) ) (1.12) In the optimization model above, the objective is to maximize the load margin λ while satisfying power flow constraints (1.2)-(1.3)and other operating constraints (1.4)-(1.12).P Gi and Q Gi are the active power and reactive power of traditional generator at bus i respectively.P Li and Q Li are the active power and reactive power of the load at bus i respectively.P Ri and Q Ri are the active power and reactive power of the renewable energy source at bus i respectively.P 0 Gi , P 0 Li and Q 0 Li are the base power of the traditional generator and load at bus i respectively.V i is the voltage of bus i. θ ij is the phase angle difference between bus i and bus j, G ij and B ij are the conductance and susceptance of the line connecting the bus i and bus j respectively.P ij and Q ij are the active and reactive power flow from bus i to bus j. t ij is the transformation ratio of line from bus i to bus j.In this optimal model, the renewable energy sources keep without power increment due to their intermittent outputs.While the powers of traditional generators and loads increase with load margin λ and the power increment direction of P d Gi , As illustrated in Fig. 1, the load margin is sensitive to the power increment direction.Different power increment directions of generator and load power lead to different load margins.
The actual operating characteristics of the power system need to be considered to determine the reasonable power increment direction used.The correlation between the stochastic variables reflects the operating characteristics of the power system.In order to maintain the correlation of the input stochastic variables, the direction of base generator and load power is adopted.
Theorem 1: the correlation coefficient of the stochastic loads maintains invariable when loads increase with the direction of base load power.
Proof: the correlation coefficient between any load i and load j with the increment direction of base load power is Theorem 2: the correlation coefficient between the stochastic loads and the renewable energy sources maintains invariable when loads increase with the direction of base load power.
Proof: in this situation, only the load power increases with the direction of base load power.The correlation coefficient between any load i and renewable energy source j is

III. MATH GENERATION CORRELATED STOCHASTIC SAMPLES FOR PROBABILISTIC VOLTAGE STABILITY ASSESSMENT
Power system probabilistic analysis is a powerful tool to discover the power system stochastic characteristics [22], [23].Monte Carlo simulation (MCS) is a traditional and reliable mathematical technique for probabilistic analysis and MCS with simple random sampling (SRS) has been used in power system probabilistic analysis [17], [24].Though MCS is always computational expensive, the result obtained by MCS is often used as a benchmark compared with other methods.Obtaining accurate samples of stochastic variables with desired correlation is most important for MCS, which is focused in this paper.When the improvement of the computational cost of MCS is needed, there are methods to be selected [24].
With the rapidly increasing on the number of stochastic variables in the power system, the correlation between the stochastic variables cannot be ignored for probabilistic analysis.Cholesky decomposition is a feasible technique and was widely used to deal with the correlation between variables in condition of positive definitive correlation matrix [16], [17] [24]- [26].But the correlation matrix is not always positive definitive in general especially when the number of the correlated variables is large [27], which makes the method infeasible.In order to deal with this problem, the polynomial transformation was used [28].However, the polynomial expansion terms are difficult to be solved especially for large scale variables.An optimal permutation technique based on Evolutionary algorithm is proposed in [27].However, the evolutionary algorithm is time consumed.We proposed a twice-permutation technique based on singular value decomposition (SVD) in [29].This method is valid and efficient even for non-positive definitive correlation matrix.
In this paper, the power method combined with Latin hypercube sampling and twice-permutation technique is proposed for probabilistic voltage stability assessment.In this proposed method, power method based transformation realizes the nonlinear correlation matrix transformation between the real stochastic variables and the standard normal variables.Sampling efficiency is improved by Latin hypercube sampling (LHS).The desired correlation coefficients are maintained by the twice-permutation technique with high efficiency and reliability.

A. POWER METHOD BASED TRANSFORMATION
Let X i be the stochastic variable of any distribution, it can be expressed as standard form below by transforming it's expected value to zero.
where E(•) is the expected value function and Var(•) is the variance value function.Y i is the standard form of X i , which has expected value of zero and variance value of one.
The standard stochastic variable can be approximated with polynomials of mth-order as below [30].
where Z i is a standard normal distribution.It is the third-order polynomial transformation when m = 3 [31], and is fifthorder polynomial transformation when m = 5 [32].c ir is the VOLUME 8, 2020 polynomial coefficient obtained by solving equations with m associated moments of Y i as shown below.
where E[Y m i ] is the mth-order moment of the stochastic variable Y i , which can be calculated when the marginal distribution or the samples of the stochastic variable are given.In this paper m is set as 5 and the simplifying equations are referenced in [32].
When samples of the correlated non-normal stochastic variable X i with correlation matrix R are desired, the intermediate correlation matrix R * of the standard normal variables Z i with the form in ( 5) must be calculated first.The equations about each nondiagonal correlation coefficient r ij in R and its corresponding intermediate correlation coefficient r * ij in R * are built in [32].

B. SAMPLING AND PERMUTATION
In order to improve the probabilistic voltage stability assessment efficiency, the samples of the n stochastic variables with a size of k shown in ( 6) can be obtained by LHS method.
The quality and accuracy of the probability analysis solutions can be affected by the correlation coefficients between the stochastic variables.Therefore, the samples for probabilistic voltage stability assessment must have the desired correlation coefficients.Cholesky decomposition based permutation method is always adopted due to its less calculating cost.However, there must be two assumptions to be satisfied.The first assumption is that samples S is independent, and the second assumption is that the correlation matrix of stochastic variables is positive definite [33].In fact, the independence assumption of initial samples S generated by LHS is always violated [26], and the desired correlation matrix is not positive definite sometimes especially when there are large scale stochastic variables in power systems [27].Therefore, the twice-permutation technique in [29] is used in this paper.The twice-permutation technique can decorrelate initial samples generated by LHS and handle with the desired correlation matrix, which even is non-positive definite.
The twice-permutation technique is described in (7).
where P, E and R * satisfy SVD equation (8).Q, D and R S satisfy SVD equation ( 9).R S is the correlation matrix of samples S.
Permutating each row of S according to matrix S * , a sample matrix Z with the same rank correlation matrix to that of matrix S * is generated [34].Then the correlation matrix of Z is close to R * .
Finally the stochastic samples X with approximate desired correlation matrix R is calculated by using ( 2) and ( 3) with substituting Z.Each column in the sample matrix X forms one set of samples to be used as inputs to the deterministic max load margin (1).

C. PROBABILISTIC VOLTAGE STABILITY ASSESSMENT METHOD
Once the stochastic samples with specified distribution and desired correlation are obtained, the probabilistic voltage stability assessment can be carried out by repeatedly deterministic calculating.Then the probability characteristics of max load margin can be obtained by statistics.The flowchart of the probabilistic voltage stability assessment method is shown in Fig. 2. The corresponding calculating steps are below.
Step 1: Input cumulative distribution function of each stochastic variable and the desired correlation matrix R.
Step 2: Generate n×k sample matrix S by LHS.
Step 4: Obtain Z by permutating S according to S * .
Step 6: Solve deterministic max load margin (1) repeatedly with all stochastic samples X.
Step 7: Calculate the probability characteristics of max load margin results by statistics.
The following mean square error index of the correlation matrix is used to measure the statistics accuracy of the samples.
where r ij is the ith row and the jth column element of the desired correlation matrix, r x ij is the ith row and the jth column element of the correlation matrix of the samples, and n is the total number of the stochastic variables.
The relative errors of means and standard deviations are adopted in (11) and (12).They are used to measure the statistics accuracy of the max load margin.where the µ a and σ a are the accurate mean and standard deviation of the max load margin respectively, which are obtained from MCS with SRS of a large enough sample size.µ s and σ s are the mean and standard deviation of the max load margin obtained by the method shown in Fig. 2. Due to the random sampling process, error indices of multiple tests are adopted to evaluate the steadiness of the proposed method as below.where m is the total number of tests.µ ρ corr , µ ε µ and µ ε σ are the means of the average errors of m tests, respectively.σ ρ corr , σ ε µ and σ ε σ are the standard deviations of average errors of m tests, respectively.

IV. SIMULATION ANALYSIS
In order to investigate the performance of the proposed probabilistic voltage stability assessment method, studies on two modified IEEE test systems are explored.Matlab is used and the codes are performed with 1.60GHz CPU.

A. THE MODIFIED IEEE 14-BUS SYSTEM
The standard IEEE 14-bus system is show in [35].For the modified test system, two wind farms are connected to bus 4 and bus 5 and two photovoltaic plants are connected to bus 9 and bus 10 respectively.The power factor of each wind farm is set as 0.95 and that of each photovoltaic plant is set as 1.0.The thermal limit of lines is set to be 100 MVA.Lower and upper bounds of bus voltage are set as 0.8 p.u. and 1.2 p.u., respectively.
The detail parameters of the wind farms and the photovoltaic plants are given in Table 1.
The Details of the load parameters are shown in Table 2 The reactive power of load is determined by the fixed power factor.
There are 15 stochastic variables in total in this case.And the correlation matrix of the stochastic variables is as the same in [29] and the corresponding intermediate correlation matrix with power method transformation is non-positive definite.The stochastic variables on the buses connected directly are assumed dependent.The correlation coefficients between different types of variables are all considered.The probabilistic results of 10 000 times of MCS with SRS are assumed to be accurate and used to calculate the error of solutions obtained by PLT, power method transformation combined with SRS and twice-permutation (PST), and Nataf transformation.Different sample sizes are studied.
The 100 tests mean square error index of the correlation matrix is shown in Fig. 3.The results without considering correlation matrix R S of samples S are also shown in Fig. 3.
The error indices of 100 tests for different types of output stochastic variables are given in Fig. 4 and Fig. 5.Moreover, the error indices with the sample size of 1500 are given in Table 3.
With sample size of 1500, the probability density and cumulative distribution of load margin are shown in Fig. 6.The results of PLT is closer to those of MCS.

B. THE MODIFIED IEEE 118-BUS SYSTEM
The standard IEEE 118-bus system is introduced in The mean of load active power is set as the deterministic load active power of the standard test system and the standard variation is set as 5%.As in the 14-bus study case above, the load reactive power is determined by fixed power factor.There are 131 stochastic variables in total in this case.The correlation matrix of the stochastic variables is the same as that in [29] and the corresponding intermediate correlation matrix with power method transformation is non-positive definite.The correlation coefficients between different variables are all considered.
The error indices of PLT and PST with sample size 1500 are given in Table 3.The same conclusions as those in the modified IEEE 14-bus system can be drawn from the simulation results.

C. RESULTS DISCUSSION
It can be seen from Fig. 3 that twice-permutation technique is valid even when the correlation matrix is non-positive definite.PLT has smaller errors compared with PST and Nataf transformation.More error produces when the stochastic samples obtained by LHS is assumed independent.The  small error of the correlation matrix shows that the adopted power increment direction reflects the correlation between the stochastic variables and shows the actual operating characteristics of the power system.
In Fig. 4, µ ε µ and µ ε σ of the max load margin with PLT are all smallest, which shows that the results with PLT is more accurate than those with PST and Nataf when the same sample size is adopted.In Fig. 5, σ ε µ and σ ε σ of the max load margin with PLT are also all smallest, which shows that the results with PLT is more stable than those with PST and Nataf.The same conclusions can also be drawn from the error indices with sample size 1500 given in Table 3.
The simulation results show that power method transformation combined with Latin hypercube sampling and twice-permutation technique can address the correlations between the input stochastic variables and provide more accurate results stably, even when the correlation matrix is not positive definite.In addition, the consuming time of the twice-permutation for 14-bus system and 118-bus system are 0.035s and 0.167s respectively, and consuming time of the once-permutation for 14-bus system and 118-bus system are 0.031s and 0.161s respectively.The twice-permutation technique is efficient.

V. CONCLUSION
In this paper, a reasonable power increment direction of generator and load power is proposed for probabilistic voltage stability assessment.With the proposed direction the actual operating characteristics of the power system are maintained from the viewpoint of the correlation between stochastic variables.In order to obtain the samples characterized with the specified distribution and the desired correlation between large scale stochastic variables for probabilistic voltage stability assessment, a method based on power method transformation combined with Latin hypercube sampling and twice-permutation is proposed.The power method transformation realizes the nonlinear correlation matrix transformation between the real stochastic variables and the standard normal variables.The twice-permutation technique decorrelates initial samples and handles with the desired correlation matrix, which even is non-positive definite.The method is tested with modified IEEE 14-bus and modified IEEE 118bus systems.The simulation results show that the proposed method is valid and feasible.

FIGURE 1 .
FIGURE 1. Illustration of voltage stability load margin for different direction of generator and load power.

FIGURE 2 .
FIGURE 2. Flowchart of the probabilistic voltage stability assessment method.

FIGURE 3 .
FIGURE 3. Error indices of 100 tests for the correlation matrix.

FIGURE 6 .
FIGURE 6. Cumulative distribution and probability density of max load margins.

TABLE 1 .
Parameters of wind farms and photovoltaic plants.

TABLE 2 .
Details of load parameters.