An Integrated LHS–CD Approach for Power System Security Risk Assessment with Consideration of Source–Network and Load Uncertainties

Abstract

Large-scale wind power integrated into power grids brings serious uncertainties and risks for power system safe operation, and it is imperative to evaluate power system security risk pertinent to high-level of uncertainties. In this paper, a comprehensive source–network–load probabilistic model, representing the typical uncertainties penetrated in power generation transmission consumption portion, is firstly set for power system operation. Afterwards an integrated LHS–CD approach based on the Latin hypercube sampling (LHS) and Cholesky decomposition (CD) is tailored to effectively conduct the security risk assessment, in which the LHS is utilized to stratified sample the uncertainties of wind power and thermal power, transmission line outage, and load demands, while the CD part is adopted to address the correlations of uncertainties by rearranging the sampled matrix generated by LHS. Moreover, static voltage risk and transmission line overloaded risk index are properly defined for quantitatively evaluating power system operational security risk. Simulation results of a modified New England 39-bus system confirm that the proposed integrated LHS–CD approach is effective and efficient for power system security risk assessment with consideration of source–network–load demand uncertainties.

1. Introduction

Power system operation usually involves various uncertainties stemming from, for instance, stochastic load demands, power source output fluctuations, and transmission line failure, etc. With the rapid increase of wind power generation worldwide, the variability and uncertainty of wind power bring heavy power fluctuations into power networks and increase power system operation risk, such as the voltage fluctuations and transmission line overloaded caused by the high level of wind power penetration [1,2]. To effectively evaluate power system security risk with consideration of various uncertainties is significant for power grid safe operation.

The objective of power system security risk assessment is to simultaneously consider the probability and severity of certain accidents from a risk perspective [1]. Although security risks cannot be completely eliminated due to always-existing uncertain perturbances in power systems, the risks can be quantized and managed at an acceptable level during the planning, design, and operation phases. In [3], a risk-based security assessment model was presented for the operational security risk assessment of a power system with the network topology, loads, and renewable generation uncertainties and solved by copula function-based Monte Carlo (MC) methodology. A bi-level optimization model was proposed in [4] to conduct the risk assessment of transmission systems, in which the lower level is to dispatch power generation for minimizing the total load shedding while the upper-level is to maximize the severity risk of the worst N-k contingency. In recent decades, the application of risk assessment technology for wind power penetrated power systems has attracted high academic interest [5]. In [6], an uncertain model for short-term wind power generation was established, and the conditional risk value-based safe distance (S-D) was intended to reveal the tail risk of system operation. On the basis of S-D, four new risk indicators were defined to reflect the risks of short-term wind power output changes in the near future. Prof. Karki et al. proposed a short-term wind power forecast method and quantified the wind power generation commitment risk in [7,8]. A risk-based admissibility assessment method was proposed in [9] to quantitatively evaluate how much wind power generation can be accommodated by a large-capacity power system under a given unit commitment strategy. In [10], Nataf transformation was used to establish a wind speed correlation model, and the effect of wind speed correlation on voltage stability margin was studied. Prof. Mehdizadeh introduced the correlation autoregressive moving average time series and copula function to predefine the correlations among wind farms, and afterward evaluated the power system static security risk with consideration of wind power correlations [11]. However, there are few reports simultaneously considering correlated wind power and load demands, as well as the uncertain network for power system risk assessment.

So far, power system security risk assessment methods can be classified into the following two branches: The analytical method and the sampling simulation method. For power systems with high complexity, the analytic method is usually very complicated and with low accuracy due to the linearization process. As a representative sampling simulation method, straightforward MC has been widely used in power system risk assessment and some improved variations of MC have also been developed [12,13,14]. The auxiliary importance sampling density function was employed in [15] to reduce the computational effort of MC for power system risk assessment. However, with the rapid development of renewable generation, many more uncertainties have penetrated power system operation, and MC has encountered the disadvantages of a huge sample size and low efficiency [16]. In recent years, the Latin hypercube sampling (LHS) method has been introduced for power system risk assessment [17], which is a hybrid of stratified strategy and random sampling [18]. In [19], an improved sequential simulation method combining the discrete LHS and importance sampling was presented for composite system reliability assessment [20].

In this paper, an integrated LHS–Cholesky decomposition (LHS–CD) based approach is proposed to quantitatively evaluate power system operational security risk with consideration of various representative uncertainties. Firstly, a source–network–load probability model considering uncertain network outages, probabilistic wind power, thermal power, and load demand is first established. And afterward, the LHS method is used to sample the uncertainties of wind power and thermal generators, probabilistic transmission line outages, and load demands; meanwhile, the Cholesky decomposition is adopted in order to rearrange the sampling matrix for dealing with wind power correlations and load demand correlations. Finally, voltage risk and transmission line overloaded risk index are defined to evaluate the power system security risk. The contributions of this manuscript are threefold.

(1)

A completed source–network–load probabilistic model with consideration of wind power and load demand correlations is proposed, and the proposed model could well reflect representative uncertainties of power system generation–transmission–consumption portion.

(2)

With the help of introducing LHS to sample the source–network–load demand uncertainties and adopting CD to handle wind power and load demand correlations by rearranging the sampling matrix, an integrated LHS–CD approach is tailored to handle power system generation–transmission–consumption uncertainties for security risk assessment.

(3)

Simulation results of the modified New England 39-bus system demonstrate that the proposed integrated LHS–CD approach is effective and efficient for power system security risk assessment with comprehensive consideration of source–network–load uncertainties.

The paper is organized as follows: Section 2 establishes the source–network–load probabilistic model for power system security risk assessment, and Section 3 presents the LHS and Cholesky decomposition strategy for correlated uncertainties processing. Section 4 proposes the nodal voltage and transmission line overloaded index to quantify power system security risks. In Section 5, simulations of New England 39-bus system validate the efficiency and effectiveness of the proposed LHS–CD approach for risk assessments, and conclusions are drawn in the last section.

2. Source–Network–Load Probabilistic Model for Power System Security Risk Assessment

This section presents the source–network–load probabilistic model for describing the uncertainties of transmission line outages, load demands, wind power generation, and thermal power generation for power system security risk assessment.

2.1. Discrete Probabilistic Model of Transmission Line Outage

In general, the operation state of power network transmission line can be modeled by a two-state running–outage model for power system risk assessment. For a transmission line r, if its outage probability is P_r, the state S_r of the line r can be expressed in discrete form as Formula (1):

$$S_r=\{\begin{array}{lll}0& \left(\mathrm{Outage}\right)& 0\le U_r<P_r\\ 1& \left(\mathrm{Run}\right)& U_r>P_r\end{array}$$

(1)

where U_r represents a random number uniformly distributed in [0,1] interval. According to the two-state model in Formula (1), the samples for the transmission line operation state are generated based on the idea of Figure 1, where X_r represents the sampled transmission line state, and the value in [0,a) indicates the outage state (S_r = 0) and in [a,b) indicates the run state (S_r = 1). Y_r represents the state probability in [0,1].

2.2. Power Source Probabilistic Model

2.2.1. Probabilistic Model of Correlated Wind Power Output

The two-parameter Weibull distribution is widely used to model wind speed uncertainties as [21]:

$$f\left(v\right)=\left(\frac{k}{c}\right){\left(\frac{v}{c}\right)}^{k-1}\exp \left[-{\left(\frac{v}{c}\right)}^k\right]$$

(2)

where v represents wind speed; k and c are the shape and scale parameters. According to Formula (2), the cumulative probability distribution function of wind speed is:

$$F\left(v\right)=1-\exp \left[-{\left(\frac{v}{c}\right)}^k\right]$$

(3)

Assume U_w is a uniform distribution in the [0,1] interval and equivalent to the cumulative probability distribution function F(v) of wind speed (U_w = F(v)), the wind speed v can be derived from Formula (4) by using the inverse transformation operator.

$$v=c{\left(-\ln U\right)}^{\frac{1}{k}}$$

(4)

where U = 1 − U_w also satisfies the uniform distribution in [0,1] interval.

It is well known that the rotor of a wind turbine starts to rotate only if the wind speed striking the turbine blades reaches the cut-in speed V_in, then the wind turbine will generate non-zero wind power output with increased wind speed. When the wind speed reaches the so-called rated wind speed V_rate, the wind turbine automatically adjusts the blade angle to hold the power at a constant level, and the corresponding power output is around the rated power P_rate. When the wind speed further increases to the cut-out speed V_cut, the wind turbine stops rotating in case of any damages to wind blades [22,23]. The typical wind speed–power curve of a wind turbine is shown in Figure 2 and can be described as Formula (5).

$$P\left(v\right)=\{\begin{array}{ll}0& 0\le v<v_{in}\\ \left(A+Bv+C{v}^2\right)P_{rate}& v_{in}\le v<v_{rate}\\ P_{rate}& v_{rate}\le v<v_{cut}\\ 0& v\ge v_{cut}\end{array}$$

(5)

where A, B, and C are constant parameters, which can be calculated from Formula (6) [24].

$$\begin{array}{l}A=\frac{1}{{\left(v_{in}-v_{rate}\right)}^2}\left[v_{in}\left(v_{in}+v_{rate}\right)-4v_{in}{v}_{rate}{\left(\frac{v_{in}+v_{rate}}{2v_{rate}}\right)}^3\right]\\ B=\frac{1}{{\left(v_{in}-v_{rate}\right)}^2}\left[4\left(v_{in}+v_{rate}\right){\left(\frac{v_{in}+v_{rate}}{2v_{rate}}\right)}^3-\left(3v_{in}+v_{rate}\right)\right]\\ C=\frac{1}{{\left(v_{in}-v_{rate}\right)}^2}\left[2-4{\left(\frac{v_{in}+v_{rate}}{2v_{rate}}\right)}^3\right]\end{array}$$

(6)

When a power network contains multiple wind farms, especially located in the same wind belt, the correlation of wind farms is significant and needs to be considered for wind power generation analysis. Assuming that a system includes p wind farms, their correlation coefficient could be described by a matrix ${\mathit{\eta }}_{WF}$ with the element η_ij representing the correlations between the ith and jth wind farm. Considering that the self-correlation coefficient for a wind farm is always 1, i.e., η_ii = 1, and the mutual correlation coefficient of any two wind farms is equal, i.e., η_ij = η_ji, the correlation matrix is described as Formula (7).

$${\mathit{\eta }}_{WF}={\left[\begin{array}{cccc}1& {\eta }_{12}& \cdots & {\eta }_{1p}\\ {\eta }_{12}& 1& \cdots & {\eta }_{2p}\\ ⋮& ⋮& & ⋮\\ {\eta }_{1p}& {\eta }_{2p}& \cdots & 1\end{array}\right]}_{p\times p}$$

(7)

It is clear that for a system containing p wind farms, the correlation matrix in Formula (7) has p(p − 1)/2 different coefficients, which can be obtained from the wind power historical data.

2.2.2. Probabilistic Model of Thermal Generator Operation State

The operation–failure two-state model is also used for the conventional thermal power generators running and outage state, i.e., S_r = 0, for outage state when $0\le U_r<P_r$ and S_r = 1 for running state when $U_r>P_r$, where U_r represents a random number uniformly distributed in [0,1] and P_r is the threshold of outage probability.

2.3. Probabilistic Model of Correlated Load Demands

The load demand in power systems is usually uncertain with time-varying fluctuations, which are generally described by a normal distribution function. In this manuscript, the probability density functions for the uncertain active load and reactive load are modeled as Formulas (8) and (9).

$$\begin{array}{l}f\left(P_i\right)=\frac{1}{\sqrt{2\pi }{\sigma }_{P_i}}\exp \left[-\frac{{\left(P_i-{\mu }_{P_i}\right)}^2}{2{\sigma }_{P_i}^2}\right]\\ f\left(Q_i\right)=\frac{1}{\sqrt{2\pi }{\sigma }_{Q_i}}\exp \left[-\frac{{\left(Q_i-{\mu }_{Q_i}\right)}^2}{2{\sigma }_{Q_i}^2}\right]\end{array}$$

(8)

where P_i is the active power of load demand I, μ_Pi and σ_Pi are the mean value and standard deviation of active power, respectively; Q_i is the reactive power of load demand i, μ_Qi and σ_Qi are the mean value and standard deviation of reactive power, respectively.

Influenced by factors such as weather conditions and electricity usage habits, load demands at nearby buses of power network usually have correlations. Assuming a system has m load buses subjected to the normal distribution in Formula (8) and any two load buses have a correlation, then the load covariance matrix C_load of the whole system can be formed with m × m elements and the load vector X of the whole network obeys the m-dimensional correlated normal distribution, which can be generated from the standard normal distribution random vector Z as Formula (9).

$$\begin{array}{l}{\mathit{C}}_{load}=\mathit{L}{\mathit{L}}^T\\ \mathit{X}=\mathit{L}\mathit{Z}+\mathit{\mu }\end{array}$$

(9)

where L is the lower triangular matrix decomposed from the covariance matrix C_load and μ is the mean vector.

3. Integrated LHS–CD Approach for Power System Security Risk Assessment

Latin hypercube sampling (LHS) is a layered sampling method, which can well cover the distribution space of random variables with small-scale sampling points. Based on effective combinations of the LHS and Cholesky decomposition (CD) strategy, an integrated LHS–CD approach is proposed for power system security risk assessment with the advantages of high efficiency. The proposed integrated LHS–CD approach mainly includes two steps to address the correlated uncertainties: (1) Sampling and (2) permutation. The objective of sampling is to ensure that the sampled points could cover the random distribution features of uncertain input variables, and the purpose of permutation is to control the correlations among these sampled points.

3.1. Sampling Based on LHS

Assuming there are K random input variables X = [X₁, X₂, ..., X_k]^T, and their probabilistic cumulative functions are Y_k = F_k(X_k), (k = 1,2, ..., K). The schematic diagram of the LHS is shown in Figure 3. For an assumed sample size N, the range of the probability of Y_k from 0 to 1 is equally divided into N non-overlapping intervals with a length of 1/N. Any point in the interval can be selected as the sample value of Y_k. In this paper, the midpoint is selected and the corresponding X_k at this point is calculated using the inverse function as Formula (10).

$$\begin{array}{cc}{X}_{kn}=F_k^{-1}\left(\frac{n-0.5}{N}\right)& n=1,2\cdots N\end{array}$$

(10)

N points are sampled for the uncertain variable X_k (k = 1, 2, ..., K) and randomly arranged in a row, then a K × N primary sampled matrix can be obtained for K input random variables as described by Formula (11).

$$X_{KN}=\left[\begin{array}{cccc}{X}_{11}& X_{12}& \cdots & X_{1N}\\ X_{21}& X_{22}& \cdots & X_{2N}\\ ⋮& ⋮& & ⋮\\ X_{K1}& X_{K2}& \cdots & X_{KN}\end{array}\right]$$

(11)

3.2. Permutation Based on CD

It is evident that (1) if the input uncertain variables are independent of each other in theory, the correlation of the sampled matrices should be correspondingly minimized; (2) when there are assumed correlations among random input variables, their sampled matrices should satisfy the very similar correlations. However, since the sampled matrix X_KN is randomly arranged and the correlation is uncontrollable, the primary sampled matrix in Formula (11) needs to be properly rearranged to satisfy the assumed correlations. CD strategy has the merits of light computation burdens and high precisions, and thus in this manuscript, the CD is used to rearrange Formula (11) for obtaining satisfied correlations according to the following procedures.

3.2.1. CD-Based Permutation for Independent Random Input Variables

The CD-based permutation procedure for independent random input variables is as follows:

(1)

Generate a K × N random ordering matrix L_r, where the element in each row is an integer number 1 to N, which represents the rank position according to the element numerical value in a descending order in that row.

(2)

Calculate the correlation coefficient matrix C_Lr of the random ordering matrix L_r. Since C_Lr is a positive definite symmetric matrix, it could be decomposed by CD to obtain the lower triangular matrix Q.

$$C_{Lr}=Q{Q}^T$$

(12)

(3)

Then a matrix G can be calculated based on the lower triangular matrix Q.

$$G=Q^{-1}{L}_r$$

(13)

It is not difficult to prove that the correlation coefficient matrix of matrix G is a unit matrix:

$$C_G=\mathrm{cov}\left(G,G^T\right)=\mathrm{cov}\left(Q^{-1}{L}_r,L_r{}^T{\left(Q^{-1}\right)}^T\right)=Q^{-1}\cdot \mathrm{cov}\left(L_r,L_r{}^T\right)\cdot {\left(Q^{-1}\right)}^T=Q^{-1}{C}_{Lr}{\left(Q^{-1}\right)}^T=Q^{-1}Q{Q}^T{\left(Q^{-1}\right)}^T=I$$

(14)

(4)

With the help of above transformation, there is no correlation between the row and column elements of matrix G, so the sample matrix of independent variables can be arranged according to G. However, since the elements in G may not be positive integers, it cannot be arranged directly according to G. The common solution is to obtain the ordering matrix L₀ of G, and then the primary sample matrix X_KN can be arranged according to L₀ for obtaining a new sample matrix X₀, which meets the characteristic of independent input variables.

3.2.2. CD-Based Permutation for Correlated Random Input Variables

The core idea of permuting correlated random input variables X = [X₁, X₂, ... X_K]^T is to create a standard normal distribution sampled matrix Z_c,KN, which has the same correlation as that of original input variables X. Since the cumulative probability function F(·) of the input random variable X and the cumulative probability function Φ(·) of the standard normal distribution are both monotones increasing functions, the ordering matrix of the standard normal distribution sampled Z_c,KN and the input random variables X is the same. Therefore, as long as the ordering matrix L_Z of Z_c,KN is obtained, the primary sample matrix X_KN in Formula (11) can be rearranged according to L_Z to obtain satisfactory correlations.

Generating a non-independent random vector Z_c = [Z₁, Z₂, ... Z_k ... Z_K]^T subject to the standard normal distribution to satisfy Formula (15).

$$\begin{array}{cc}{Z}_k={\Phi }_k^{-1}\left(F_k\left(X_k\right)\right)& k=1,2\cdots K\end{array}$$

(15)

where Φ_k(·) is the cumulative probability distribution function of the standard normal distribution, and F_k(·) is the cumulative probability distribution function of the input random variable X_k. Assuming that the correlation coefficient matrix of input random variables X = [X₁, X₂, ... X_K]^T is C_X and the correlation coefficient matrix of non-independent random vector Z_c = [Z₁, Z₂, ... Z_K ]^T is C_Z:

$${\mathit{C}}_X={\left[\begin{array}{cccc}1& {\rho }_{12}& \cdots & {\rho }_{1K}\\ {\rho }_{21}& 1& \cdots & {\rho }_{2K}\\ ⋮& ⋮& & ⋮\\ {\rho }_{K1}& {\rho }_{K2}& \cdots & 1\end{array}\right]}_{K\times K}$$

(16)

$${\mathit{C}}_Z={\left[\begin{array}{cccc}1& {\rho }_{12}^{\prime }& \cdots & {\rho }_{1K}^{\prime }\\ {\rho }_{21}^{\prime }& 1& \cdots & {\rho }_{2K}^{\prime }\\ ⋮& ⋮& & ⋮\\ {\rho }_{K1}^{\prime }& {\rho }_{K2}^{\prime }& \cdots & 1\end{array}\right]}_{K\times K}$$

(17)

And ρ_ij is the correlation coefficient for input random variables X calculated by Formula (18).

$${\rho }_{ij}=\frac{\mathrm{cov}\left(X_i,X_j\right)}{{\sigma }_i{\sigma }_j}$$

(18)

where σ_i and σ_j are the standard deviations of random variables X_i and X_j, respectively. And the off-diagonal elements of C_X and C_Z satisfy the following relation:

$${\rho }_{ij}^{\prime }=T\left({\rho }_{ij}\right){\rho }_{ij}$$

(19)

where T(ρ_ij) is related to the distribution of X_i and X_j, if X_i and X_j are the normal distributions, T(ρ_ij) = 1; if X_i and X_j are Weibull distributions, T(ρ_ij) is approximately equal to Formula (20) [25].

$$T\left({\rho }_{ij}\right)=1.063-0.004{\rho }_{ij}-0.2\left(\frac{{\sigma }_i}{{\mu }_i}+\frac{{\sigma }_j}{{\mu }_j}\right)-0.001{\rho }_{ij}{}^2+0.337\left(\frac{{\sigma }_i^2}{{\mu }_i^2}+\frac{{\sigma }_j^2}{{\mu }_j^2}\right)+0.007{\rho }_{ij}\left(\frac{{\sigma }_i}{{\mu }_i}+\frac{{\sigma }_j}{{\mu }_j}\right)-0.007\frac{{\sigma }_i{\sigma }_j}{{\mu }_i{\mu }_j}$$

(20)

where μ_i and μ_j are the mean value of random variables X_i and X_j, respectively.

The correlation coefficient matrix C_Z of the non-independent random vector Z_c can be easily obtained by CD to get its lower triangular matrix B, i.e., C_Z = BB^T. Then, Z_c can be obtained from a vector of independent standard normal variables w = [w₁, w₂, ..., w_K]^T as Formula (21).

$${\mathit{Z}}_c=\left[\begin{array}{c}{Z}_1\\ Z_2\\ ⋮\\ Z_K\end{array}\right]=\mathit{B}\mathit{w}=\left[\begin{array}{cccc}{b}_{11}& & & \\ b_{21}& b_{22}& & \\ ⋮& ⋮& \ddots & \\ b_{K1}& b_{K2}& \cdots & b_{KK}\end{array}\right]\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_K\end{array}\right]$$

(21)

Then, the independent standard normal distribution random variables w_i (i = 1,2…K) are sampled N times to obtain the sample matrix W_KN, and the normal distribution sample matrix Z_c,KN with the correlation coefficient matrix C_Z can be obtained by Formula (22).

$${\mathit{Z}}_{c,KN}={\left[\begin{array}{cccc}{Z}_{11}& Z_{12}& \cdots & Z_{1N}\\ Z_{21}& Z_{22}& \cdots & Z_{2N}\\ ⋮& ⋮& & ⋮\\ Z_{K1}& Z_{K2}& \cdots & Z_{KN}\end{array}\right]}_{K\times N}=\mathit{B}{\mathit{W}}_{KN}=\left[\begin{array}{cccc}{b}_{11}& & & \\ b_{21}& b_{22}& & \\ ⋮& ⋮& \ddots & \\ b_{K1}& b_{K2}& \cdots & b_{KK}\end{array}\right]\left[\begin{array}{cccc}{W}_{11}& W_{12}& \cdots & W_{1N}\\ W_{21}& W_{22}& \cdots & W_{2N}\\ ⋮& ⋮& & ⋮\\ W_{K1}& W_{K2}& \cdots & W_{KN}\end{array}\right]$$

(22)

The ordering matrix L_z of Z_c,_KN is a K × N matrix, and the elements in each row are a non-repeated number of integers 1 to N. Afterward, by arranging the sampled matrix X_KN according to L_z, the rearranged matrix will have the required correlations of random input variables X = [X₁, X₂, ... X_K]^T.

The main steps of the permutation method for handling the correlated input random variables are as follows:

(1)

Calculate the correlation coefficient matrix C_Z from C_X according to Formulas (16)–(20);

(2)

Use CD to obtain the lower triangular matrix B of C_Z, i.e., C_Z = BB^T;

(3)

Sample K independent random variables of standard normal distribution N times to obtain W_KN;

(4)

Calculate the normal distribution sample matrix Z_c,KN = BW_KN in Formula (22), which has the same correlations of the original input random variables X;

(5)

Obtain the ordering matrix L_z of Z_c,KN;

(6)

Rearrange the sample matrix X_KN generated in Formula (11) according to L_z, and the perfect sampled matrix is finally obtained.

3.3. LHS–CD for Generating the Correlated Sampling Matrix for Power System Source–Network–Load Uncertainties

In order to conduct the risk analysis, the uncertain input variables shall be sampled properly in terms of a sampling matrix, which should reflect the completed power system source–network–load uncertainties. Based on the probabilistic model of transmission network outages, wind power, and thermal power output, as well as the load demands presented in Section 2, the sampling matrices of thermal generators and transmission lines are generated as independent input variables according to the procedure in Section 3.2.1, and the sampling matrices of correlated wind speeds (which was later converted into wind power samples based on the wind speed–power curve according to Formulas (5) and (6)) and correlated load demands are generated according to the steps in Section 3.2.2. The final input variable sampling matrix X_F is generally formed as:

$${\mathit{X}}_F={\left[X_{W1},\cdots X_{Wp},X_{G1},\cdots X_{Gq},X_{L1},\cdots X_{Ls},X_{D1},\cdots X_{Dm}\right]}^T$$

(23)

where X_W1, ..., X_Wp are for wind power samples; X_G1, ..., X_Gq are thermal generators running state samples; X_L1, ... X_Ls are transmission line operation–outage state; X_D1, ... X_Dm are load demand samples; p, q, s, and m are the numbers of wind power generation, thermal generators, transmission lines, and load demands, respectively, and the total number of input random variables is K, i.e., K = p + q + s + m. Since each row of X_F represents one sampling vector for power system source–network and load uncertainties, it could be used in the power flow procedures to conduct power system security risk analysis in the next section.

4. Index for Power System Security Risk Assessment

4.1. Static Voltage Risk Index

The risk index of bus voltage can quantify the severity of overvoltage or under voltage, which is defined as the product of the severity SE_Uik of bus i voltage deviation and the corresponding probability P_Uik:

$$R_{Ui}={\displaystyle \sum _{k=1}^N{P}_{Uik}S{E}_{Uik}}={\displaystyle \sum _{k=1}^N{P}_{Uik}}\left|\frac{U_{ik}-U_{i\_rate}}{U_{i\_rate}}\right|$$

(24)

where N is the number of sample points, and R_Ui is the overridden voltage risk index of the bus i; U_ik is the voltage calculated from power flow for the kth sampling point; U_{i_rate} is the rated voltage of the bus i; P_Uik is the probability of sample k. For an LHS–CD approach with N sampling points, P_Uik is calculated as $P_{Uik}=\frac{1}{N}$.

4.2. Transmission Line Overloaded Risk Index

The line overloaded risk index R_Lj of line j can be defined as the product of the transmission line overloaded severity SE_Ljk and the corresponding probability P_Ljk, and Formula (25) can be used to quantify the transmission line overloaded risk for overall N sampling points of power system uncertainties.

$$R_{Lj}={\displaystyle \sum _{k=1}^N{P}_{Ljk}S{E}_{Ljk}}={\displaystyle \sum _{k=1}^N{P}_{Ljk}}\times \frac{\max (S_{jk}-S_{j\_rate},0)}{S_{j\_rate}}$$

(25)

where S_jk is the apparent power flow of transmission line j for the kth sampling point; S_{j_rate} is the rated apparent power. Similarly, for an LHS–CD approach with N sampling points, P_Ljk is calculated as P_Ljk = 1/N.

4.3. Main Steps of Proposed LHS–CD Approach for Power System Static Risk Assessment with Consideration of the Source–Network–Load Uncertainties

The main steps of the proposed LHS–CD approach for power system risk assessment are shown in Figure 4.

5. Simulation Results

The effectiveness of the proposed method is verified on the New England 10-generator 39-bus system [26], and the system base power is 100 MVA and the base voltage is 345 kV. As indicated in Figure 5, there are four wind farms connected to buses 11, 13, 18, and 22, respectively. Wind speed of these four wind farms is assumed as the Weibull distribution with parameters detailed in

Table 1. Wind speed parameters of four wind farms.

Wind Farm Connected Bus	Shape Parameter k	Scale Parameter c	Cut-In Speed Vin (m/s)	Cut-Out Speed Vcut (m/s)	Rated Speed Vrate (m/s)	Rated Power Prate (MVA)
11	2	12	3	25	12	200
13	2	12	3	25	12	200
18	2	10	3	25	12	200
22	2	10	3	25	12	200

, and the correlation coefficient of four wind farms is ${\mathit{\eta }}_{WF}=\left[\begin{array}{cccc}1& 0.1& 0.3& 0.2\\ 0.1& 1& 0.2& 0.1\\ 0.3& 0.2& 1& 0.3\\ 0.2& 0.1& 0.3& 1\end{array}\right]$. The outage probability P_r of thermal generators and transmission lines are 0.001 and 0.02, respectively. The probabilistic loads are modeled by a multi-dimensional normal distribution, the baseload given in [26] is assumed as the mean value with the standard deviation fixed as 10% of the mean value, and the correlation coefficients of the 21 correlated loads are assumed as 0.5. In order to verify the efficiency and accuracy of LHS–CD, the standard MC is adopted as a benchmark for risk assessment, and both the LHS and MC methods are run on a PC with a 2.8 GHz Intel Core i5-8400 CPU and 8 GB RAM. The power flow calculation of LHS–CD and MC risk assessment are both conducted by the modified Matpower and Psat in MATLAB R2016a.

5.1. Comparison of MC and LHS–CD for Wind Speed Sampling

In order to verify the effectiveness of the LHS–CD approach, wind speed is sampled by LHS–CD and MC, respectively, for the same scale of 3000 sampled points. The actual wind speed and the sampled wind speed probability curves of Wind Farm 1 by LHS–CD and MC are, respectively, shown in Figure 6.

From the comparisons in Figure 6, it is clear that with the same number of sampling points, wind speed generated by LHS is much more close to the actual wind speed, while the wind speed obtained by MC has large deviations. This is because LHS–CD is a stratified sampling method, which can fully reflect the distribution features of random variables. However, MC is just a random sampling method, and its sample points cannot completely cover the density distribution function when the sampling scale is not large enough. Compared with MC, LHS–CD results are smoother and more accurate with the small-scale sampled points.

5.2. LHS–CD Convergence for Power Security Risk Assessment

In order to clearly compare the performance of MC and LHS–CD for system security risk assessment, we investigated the convergence of risk index under different sampling scales. Taking the voltage risk index of Bus 37 as an example, the convergence curves of MC and LHS–CD with increasing sampling scale are shown in Figure 7. It is clear that LHS–CD has a faster convergence than MC, and the former could reach a stable value with several hundreds of sampling points. Since LHS–CD is based on a stratified sampling strategy, only a few sampling points are required to well reflect the entire probability distribution characteristics of input variables. It is clear that the LHS–CD approach has the good convergence for sampling the system uncertainties and evaluating system security risk.

In order to verify the advantage of the LHS–CD-based risk assessment method, MC is deployed as a benchmark for power system security analysis considering source–network–load system uncertainties. According to convergence theory in [15], the required number of samples for MC is estimated to be 20,000 for satisfying an accuracy of 0.01, while the sampling scale of LHS–CD is determined by the mean and standard deviation errors ${\epsilon }_{\mu }^x=\left(\left|{\mu }_{MC}-{\mu }_{LHS}\right|\right)/{\mu }_{MC}$ and ${\epsilon }_{\sigma }^x=\left(\left|{\sigma }_{MC}-{\sigma }_{LHS}\right|\right)/{\sigma }_{MC}$ benchmarked with MC, where ${\mu }_{MC}$, ${\mu }_{LHS}$ and ${\sigma }_{MC}$, ${\sigma }_{LHS}$ are the mean and standard deviation of random input variables sampled by MC and LHS, respectively. After testing, LHS–CD can meet the same accuracy when there are only 500 sampling points, and this also indicates that the LHS–CD approach has good accuracy.

5.3. Influence of Wind Power Correlation and Load Correlation

To clearly reflect the influence of correlated wind speeds on total wind power, the probability of total wind power of four wind farms is estimated by LHS–CD under two conditions: Regardless of wind speed correlations and with consideration of wind speed correlations. As shown in Figure 8, there are noticeable differences in the probability curve of total wind power with and without consideration of the wind speed correlations. When the correlations are taken into account, the probability of the total wind power is lower in the middle part and higher on the two sides. This is because when wind speed of one wind farm is larger (or smaller), the wind speeds in other wind farms are correspondingly larger (or smaller) due to the positive correlations, so the probability of the total wind power is raised for the larger (or smaller) part, and the probability density curve considering correlations is smoother. Based on these observations, the wind speed correlations for different wind farms in nearby areas could affect the total wind power output, and they should be considered in the security risk assessment of power systems penetrated with multiple wind farms.

Owing to the influence of climate, geographical factor, and electricity usage habits, etc., there are obvious correlations among multiple load demands at different locations. The total load demand is sampled by LHS–CD, and the total load demand probability with and without consideration of load correlations are shown in Figure 9. It is quite clear that there are distinguished differences in the total load demand whether to consider the load correlations.

In order to directly verify the influence of correlation on power system security risk, all transmission lines overloaded risk index in New England 39-bus system is evaluated with and without considering the wind power and load demand correlations. As demonstrated by the compared results in Figure 10, the transmission line overloaded risk index for the condition ignoring the wind power and load demand correlations is quite different from the results of considering these correlations. If correlation is not taken into account, critical errors will be caused, making the evaluation results inaccurate and inconsistent with the actual situation. Therefore, the wind power correlation and load correlation should be considered in the static security risk assessment of power system with source–network–load uncertainties.

5.4. Accuracy of Security Risk Assessment

In order to verify the accuracy and efficiency of the proposed LHS–CD method, the standard MC is also used to conduct the security risk assessment under the same conditions. According to the sampling scale analysis obtained by the previous simulation in Section 5.2, 500 samples are performed for LHS–CD method and 20,000 samples for the MC method to obtain the wind speed curve of Wind Farm 1, as shown in Figure 11. It is clear that the two methods can roughly accurately simulate the actual wind speed.

In addition, when comparing Figure 11 with Figure 6, it can be seen that increasing the sampling scale of MC can effectively improve its sampling accuracy, and the wind speed simulation curve is closer to the actual wind speed. While the number of LHS–CD samples is reduced from 3000 to 500, it can still accurately simulate the wind speed. Furthermore, LHS–CD and MC are further used to sample the uncertainties of source–network–load mentioned in Section 2, and then the defined risk indexes, including the voltage risk index of 39 buses and the overloaded risk index of 46 transmission lines, are calculated according to Formulas (24) and (25) and shown in Figure 12 and Figure 13.

From voltage risk indexes of 39 buses in Figure 12 and the overloaded risk indexes of 46 transmission lines in Figure 13, it is obvious that the risks calculated by LHS and MC are consistent with each other. These comparisons demonstrate that the LHS–CD-based risk assessment method is quite accurate. However, in terms of time-consuming for security risk assessment, the LHS–CD-based security risk method costs 246.10 s, of which 28.91 s is spent on sampling and permutation and 217.19 s is spent on solving power flow. Meanwhile, the MC-based security risk method spends 9245.29 s, of which 156.27 s is spent on sampling and permutation and 9089.02 s is spent on solving power flow. These investigations indicate that the LHS can accurately reflect the uncertain features of random input variables with only a small number of samples, and therefore, the proposed LHS–CD-based risk assessment method is highly accurate and efficient.

6. Conclusions

In this paper, a completed power system source–network–load demand probabilistic model representing the uncertainties of wind power and thermal power output, transmission line outage, and load demand is firstly established, and afterward, an integrated LHS–CD approach combining Latin hypercube sampling with Cholesky decomposition is proposed to evaluate power system probabilistic security by the defined nodal voltage risk and transmission line overloaded risk. Simulation results of New England 39-bus system have indicated that (1) the wind speed correlations and load correlations have a great impact on the total wind power output and total load demand, and these correlations should be considered in power system security risk assessment; (2) LHS can well reflect the overall probabilistic distribution of random variables via small-scale sampling points, and based on the defined voltage risk index and line overloaded risk index, the proposed integrated LHS–CD approach is effective and efficient for evaluating power system operation risk with consideration of source–network–load uncertainties.