Probabilistic load flow using improved three point estimate method

Electrical Power and Energy Systems 117 (2020) 105618

Contents lists available at ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

Probabilistic load flow using improved three point estimate method Yulong Che a b

a,b

a,⁎

a

, Xiaoru Wang , Xiaoqin Lv , Yi Hu

T

a

School of Electrical Engineering, Southwest Jiaotong University, Chengdu, Sichuan 610031, People’s Republic of China School of Automation and Electrical Engineering, Lanzhou Jiaotong University, Lanzhou, Gansu 730070, People’s Republic of China

A R T I C LE I N FO

A B S T R A C T

Keywords: Probabilistic load flow Point-estimate method Monte Carlo simulation Decorrelation Non-normal

An improved three point-estimate method (I3PEM) scheme is proposed to estimate the probability moments of probabilistic power flow (PLF) in this paper. I3PEM is obtained by the existing three point-estimate method (3PEM), two point-estimate method (2PEM) and Chebyshev inequality. With only first three order moments (mean, standard deviation and skewness) of input random variables, the additional calculations in I3PEM are performed by a newly added pair of estimate points to improve the accuracy of estimated moments for output random variables. It is avoided the calculation of higher order center moments and the non-real number solutions that may occur with normalized center distances. The case studies of IEEE 14-bus and 118-bus system are employed to verify the performance of I3PEM in the presence of Normal distribution inputs, correlated inputs and non-normal inputs. Under the different coefficient of variation (CV) conditions, the results of MCS are used as benchmark. The accuracy of results obtained by I3PEM have been validated by comparing with those obtained from the basic 3PEM and 2PEM. Further, the advantages and applicability of I3PEM are compared.

1. Introduction With the development of renewable energy power generation and electric vehicles, the uncertainty factors in the power system have greatly increased, and the calculation of deterministic load flow cannot fully reflect the uncertainty operation of the grid [1,2]. The probabilistic load flow (PLF) firstly proposed by Borkowska in 1974 has become one of effective means to analyze the stochastic characteristics of grid [3]. The essence of PLF is to solve high-dimensional nonlinear equations with random variables, whose purpose is to obtain the change information of bus voltage and line current. The PLF calculation can obtain more abundant statistical information, such as expectation, standard deviation, distribution interval, the most likely numerical value, etc. [4], compared with the traditional deterministic load flow. Currently, PLFs have been successfully applied in areas such as renewable energy generation [5], stochastic programming [6] and risk assessment [7] of grid. Therefore, further research and development of PLF calculation methods are of great significance for the application of PLF. The commonly used methods of PLF calculation are mainly divided into three categories [8,9]: Monte Carlo simulation (MCS) method, analytical method and approximation method. MCS [10] transforms the probability problem into a series of deterministic computational problems by sampling the random variables. The MCS method is

⁎

computationally flexible and suitable for complex stochastic models. And the MCS method is often used as a benchmark for evaluating the pros and cons of other probabilistic analysis methods because it is the most accurate algorithm [11]. However, the MCS method are time consuming and computationally inefficient as large amount of random sampling and iterative calculations are required to achieve satisfactory calculation results. In order to improve the computational efficiency of MCS, there are mainly two improved aspects: sampling strategy more suitable for PLF calculation [11,12] and effective robust correlation control technology [13,14]. Analytical methods include convolution method [15], Fourier transform method [16] and cumulant method [17]. Convolution operation and correlation transformation in convolution methods are relatively complex and computationally intensive. The commonly used cumulant method is to linearize the nonlinear equation of load flow at the reference operating point to obtain a linear model of the output and input random variable. The cumulant of output variables are obtained from linear model and the cumulant of input variables. And then the probability distribution functions (PDF) of the output variable are obtained by combining the appropriate series expansion method [18]. The advantage of the cumulant method is that the amount of calculation is small since only one times deterministic power flow calculation is performed. The calculation error of the cumulant method is independent of the grid size, while it is related to the fluctuation range of

Corresponding author. E-mail address: [email protected] (X. Wang).

https://doi.org/10.1016/j.ijepes.2019.105618 Received 22 April 2019; Received in revised form 23 August 2019; Accepted 7 October 2019 0142-0615/ © 2019 Elsevier Ltd. All rights reserved.

Electrical Power and Energy Systems 117 (2020) 105618

Y. Che, et al.

calculations. Under the different coefficient of variation (CV) conditions, the results of MCS are used as benchmark. The case studies of IEEE 14-bus and 118-bus system are employed to verify the performance of I3PEM in the presence of different type of random variables with Normal or Non-normal (bi-modal) distributed, correlated or not. The accuracy of results obtained by I3PEM have been validated by comparing with those obtained from the basic 3PEM and 2PEM. The remainder of this paper is organized as follows. The most basic PEMs such as 2PEM and 3PEM are presented in Section 2, as well as the I3PEM scheme are proposed. Section 3 presents the procedure of PEMs and MCS for the PLF with correlated inputs variables. Case studies and discussion, based on the IEEE 14-bus and 118-bus test systems in the presence of different type of random variables considering Normal or Non-normal distributed, correlated or not, respectively, are presented in Section 4. Finally, Section 5 is the conclusion part.

the input variables. Linearization in cumulant method inevitably introduces errors. When the input variable fluctuates greatly, the error will be large. The main limitation of cumulant method is that the correlation between random variables cannot be handled well. Therefore, variable correlation processing in analytical methods is a hot topic [19,20]. Approximation methods approximate the statistical features of state variables in system by using that of the input random variables. These methods avoid large-scale repeated sampling, so the calculation speed is faster. It is also took more attention because it can better take into account the correlation between input variables. Point estimate methods (PEMs) [21], unscented transformation method (UTM) [22] and Taguchi method [23,24] are the currently common approximation methods [25]. The advantage of UTM is that the correlated random variables can be processed directly. And like PEMs, there is no need to know the probability distribution of the input random variables. However, there is no better algorithm to determine the optimal or better parameter δ for the UTM. The amount of calculation of UTM is equivalent to the Km + 1 scheme of PEMs, but it is difficult to extend to the case of K = 4, 6, …. The Taguchi method was used in Ref. [23] without considering that the input random variables might be correlated. Ref. [24] extended the Taguchi method to unbalanced electrical distribution systems with correlated inputs, while it didn’t account how to obtain the PDF or CDF of output random variables. As with the MCS method, the most popular PEMs use the deterministic calculation to solve the probability problem. Only some statistical moment information of the probability function are used in PEMs, such as mean, variance, skewness and kurtosis. The advantage of PEMs compared to the MCS method is that it requires less data, then it has less computational burden and time consuming. Moreover, improved PEMs can better handle the correlation of input random variables. However, the high-order moment error of the output random variable by PEMs is larger, and the probability distributions can be not until obtained with the help of suitable approximation/expansion method [26,27]. Ref. [21] pointed out that the use of fewer estimate points m sometimes leads to inaccuracies in estimates of mean and standard deviation. In theory, the increased number of estimate points m can improve the accuracy of PEMs. However, as the value of m increases by one, the higher-order moment of the input random variable is required. The normalized center moment calculated from the fourth order or more center moment does not have much physical meaning in practice. Meanwhile, it will bring the following problems [28,29]:

2. Point estimate methods Let Z = H (X) = h (X1 , X2 , …, Xn ) denotes a nonlinear function between the input random variable X and the output random variable Z , where X represents the set of random variables Xk, k = 1, 2, …, n . The core idea of PEMs is to obtain the moments of the output random variable Z according to the probability distribution of the input random variable X [33]. According to the statistical moment of X , m × n times calculations are performed at m points to obtain the moments of Z. Let the location xk, i (k = 1, 2,…, n, i = 1, 2,…, m) of m estimate points be extracted from the input random variables Xk. And wk, i is the weight (or concentrations) corresponding to each estimate point located at(μ1 , μ 2 , …, xk, i , …, μn ) . The location xk, i is the ith value in the input random variable Xk of the function H to be estimated, which can be represented by the mean and standard deviation of the input data, where

xk, i = μk + ξk, i σk

(1)

where ξk, i is the standard location. μk and σk are the mean and standard deviation of the input random variable Xk. The weight wk, i is a weighting factor that takes into account the relative importance of this estimate point in the output random variable [31]. Since the sum of wk, i is one, we have n

m

∑ ∑ wk , i = 1

(2)

k=1 i=1

(1) It is difficult to obtain high-order statistics for complex probability distribution functions; (2) When m > 3, the analytical solutions of the standard locations and weights for input random variable cannot be obtained; (3) The used higher order moment may probably cause the non-real or negative value of the locations, as well as that of the corresponding weights.

The standard location ξk, i and weight wk, i can be obtained by solving the following equations [23]

⎧ ⎪

m

∑ wk, i = 1/ n, k = 1, 2, …, n i=1

⎨m j ⎪ ∑ wk, i (ξk, i ) = λk, j , j = 1, …, 2m − 1 ⎩i=1

Two point-estimate method (2PEM) with m = 2 and three pointestimate method (3PEM) with m = 3 are generally selected in practical problems [21,30,31]. 2PEM and 3PEM are also called 2m and 2m + 1 concentration scheme, respectively. 2PEM only uses the first threeorder moments information of input random variable and requires the deterministic calculations of 2m times [30]. 3PEM utilizes the first fourorder moments information of input random variable, and the total number of calculations is only 2m + 1 times. This is because only one calculation is equivalent to the effect of m calculations when the estimate points of all random variables take the mean value [31,32]. To address the above issues, an improved three point-estimate method (I3PEM) scheme is proposed for PLF calculation to obtain higher accuracy in this paper. Only considering the first four-order moment information of input random variables, the newly added pair of estimate points are used to obtain the higher accuracy by additional

(3)

where λk, j is the ratio of the jth order central moments Mj (xk ) to the jth power of standard deviation (σk ) j , that is

⎧

λk , j =

Mj (xk ) (σ k ) j

⎨ M (x ) = ∫+∞ (x − μ ) jf (x )dx k k −∞ ⎩ j k

(4)

From the above (4), λk,1 = 0 , λk,2 = 1. λk,3 and λk,4 are the skewness and kurtosis of random variables xk , respectively. The skewness and kurtosis of the random variable are 0 and 3 for the normal probability distribution, respectively. Once the standard location ξk, i and the weight wk, i of each estimate point location are calculated, the function Z = H (X) can be estimated at the points (μ1 , μ 2 , …, xk, i , …, μn ) . And then the estimated values of the moments of the output random variable Z can be obtained 2

Electrical Power and Energy Systems 117 (2020) 105618

Y. Che, et al. n

E (Z l ) ≈

m

∑ ∑ wk,i·[H (μ1, μ2, …, xk,i, …, μn )]l k=1 i=1

accurate when the more values of points m in theory. But for every increment of points m, it means a higher order moment need to be calculated. In case of five point-estimate method (or 4m + 1 concentration scheme), for example, until the eighth order central moment of the input random variables are necessary to known, so a greater computational burden is needed [32]. At the same time, the calculated λk, j of the higher order center moment above the fourth order does not have much physical meaning in practice. Moreover, the calculated ξk, i and wk, i probably are non-real solutions [28,29]. The I3PEM scheme is proposed in this paper based on the existing 3PEM. Only a new pair of estimate points are added for additional calculations to improve the calculation accuracy of 3PEM, while it is no need to calculate the higher order moments above the fourth order of the input random variables. There are two improved thinking. One is that the newly added estimation points only consider the first two order moments (mean and variance), and the other is that the newly added estimation points consider the first three order moments (mean, variance and skewness). Their implementation and differences are introduced as follows. Firstly, the one implementation process of I3PEM is that a new pair of estimate points with uniform probability are added for additional calculations based on the original 3PEM in Section 2.1. By (2), let the weight corresponding to this pair of each estimate point take the value be wk′, i = 1/2n . According to the definition of the location in (1), and combining with Chebyshev's inequality, it is known that the location xk, i is less likely to appear outside the several times standard deviation of the mean μk [34]. Therefore, only the first two order moments of this newly added pair of estimate points are considered, that is the mean and standard deviation, while the moments of the third order and above ′ and ξk,2 ′ for this new are ignored. The values of standard locations ξk,1 pair of estimate points are taken as follows:

(5)

where l = 1, 2, 3 …. when l = 1, E (Z ) is the mean of Z. When l = 2 , the standard deviation of Z is

E (Z 2) − E (Z )2

σZ =

(6)

With the above calculated moments as in (5), the PDF and CDF of output random variable Z are computed using the approximation/expansion methods, such as Gram-Charlier expansion and Cornish-Fisher expansion. 2.1. Two point-estimate method (2PEM) When m = 2, that is, only two locations are extracted for each input random variable. This PEM is called two point-estimate method, 2PEM. Since 2n estimate points are used for calculation, it is also called 2nconcentration scheme. By combing (2) and (3), the standard location ξk, i ( k = 1, 2, …, n, i = 1, 2) and the corresponding weight wk, i of each location xk, i can be obtained.

ξk, i =

λk2,3 λk,3 + (−1)3 − i n + , i = 1, 2 2 4

wk , i =

i (−1)iξk,3 − i 1 (−1) ξk,3 − i 1 = , i = 1, 2 n ξk,1 − ξk,2 n 2 n + (λk,3/2)2

(7)

(8)

The significant advantage of 2PEM is that it is simple and less computationally. However, the standard locations ξk,1 and ξk,2 depend on the number n of input random variables from (7). The locations xk,1 and xk,2 are farther from the mean by the ratio n when n increases. Then the location xk, i , (i = 1, 2) may fall outside of its domain.

ξk′,1 = −ξk′,2 =

2.2. Three point-estimate method (3PEM)

k, i

⎨ ⎩

=

λk ,3 2

+ (−1)3 − i λk,4 −

3λk2,3 4

xk′, i = μk + (−1) k n σk ·k = 1, 2, …, n, i = 1, 2

2

n

E (Z l ) ≈

1⎧ ∑ ∑ wk,i [H (μ1, μ2, …, xk,i, …, μn )]l + w0 [H (Xμ)]l 2 ⎨ k=1 i=1 ⎩ +

1 2n

2

n

∑ ∑ [H (μ1, μ2, …, x ′k,i, …, μn )]l ⎫ k=1 i=1

⎬ ⎭

(13)

It is shown by (13) that the final expected values of the output random variables are the arithmetic mean after linear superposition of all function values. Secondly, the other implementation process of I3PEM is that also based on the original 3PEM. However, a new pair of estimate points without uniform probability are added for additional calculations. Unlike the first implementation process of I3PEM, the weight wk′, i, (i = 1, 2) and the standard location ξk′, i of the newly added pair of locations xk′, i = μk + ξk′, i σk consider the influence of the skewness. And the weight wk′, i and standard location ξk′, i directly take the values of (7) and (8) of 2PEM, respectively. Thus, (5) becomes

(9)

3−i

(−1) ⎧ wk , i = , i = 1, 2 (ξk , i (ξk ,1 − ξk ,2)) ⎪ n n 1 ⎨w = ∑ w = 1 − ∑ 0 k ,3 λk ,4 − λk2,3 ⎪ k=1 k=1 ⎩

(12)

Thus, the estimated values of the moments for the output random variable Z in (5) becomes

, i = 1, 2

ξk,3 = 0

(11)

′ and ξk,2 ′ of this pair of It can be seen that the standard locations ξk,1 newly added estimate point are variants of the standard location ξk, i of 2PEM in (7). Then, the additional pair of locations is

When m = 3 and the standard location ξk,3 = 0 , this PEM is called three point-estimate method, 3PEM. ξk,3 = 0 means xk,3 = μk that can be seen from (1), thus n of 3n locations are the same point (μ1 , μ 2 , …, μk , …, μn ) . Therefore, by updating the corresponding weight n of this location xk,3 = μk to w0 = ∑k = 1 wk,3 , the calculation effect of only the one time at this location is equivalent to that of the same n times. The number of calculations of 3PEM is only once more than 2PEM, so it is also called 2n + 1-concentration scheme. The analytical expression for the standard location ξk, i, ( k = 1, 2, …, n, i = 1, 2, 3) and the corresponding weight wk, i of each location xk, i are:

⎧ξ

n

(10)

The standard locations ξk, i of 3PEM depend only on the skewness λk,3 and kurtosis λk,4 as in (9), while they have nothing to do with the number n of input random variables. Although the number of calculations of 3PEM is only once more than 2PEM, 3PEM is more accurate than 2PEM because it considers the kurtosis λk,4 of the input random variables. It can also be seen from (9) that when λk,4 − 3λk2,3/4 < 0 , 3PEM will obtain the non-real number solution of the standard locations ξk, i .

2

n

E (Z l ) ≈

1⎧ ∑ ∑ wk,i [H (μ1, μ2, …, xk,i, …, μn )]l + w0 [H (Xμ )]l 2 ⎨ k=1 i=1 ⎩ n

+

2

∑ ∑ w′k,i [H (μ1, μ2, …, x ′k,i, …, μn )]l ⎫ k=1 i=1

2.3. Improved three point-estimate method (I3PEM)

⎬ ⎭

(14)

Obviously, (14) is better than (13) because the newly added estimate points as in (14) consider more information of input random

The estimate value of output random variable Z can be more 3

Electrical Power and Energy Systems 117 (2020) 105618

Y. Che, et al.

converts the set of input correlated random variables into an uncorrelated set of random variables. Once the uncorrelated variables are obtained, the PEMs can be applied. Let a vector of n correlated random variables x = (x1, x2, …, x n )T with a mean vector μ x = (μ1 , μ 2 , …, μn )T , and a correlation matrix as

variables, such as the skewness. Compared with the corresponding standard location and weight of each location as in (13), it is not uniform probability for the corresponding standard location ξk′, i and weight wk′, i of each estimate point xk′, i = μk + ξk′, i σk as in (14). ξk′, i and wk′, i as in (14) take the value as in (7) and (8), respectively, and they both consider the role of the skewness of input random variables. I3PEM can be seen as a combination of 2PEM and 3PEM. The computational complexity of the I3PEM is the same as that of the 4m + 1 concentration scheme. However, the 4m + 1 scheme which needs to calculate to the seventh-order moment and there is not analytical formula to compute the concentration values [31]. Because I3PEM avoids the calculation of higher order moments greater than fourth order moment and the possible non-real number solutions of λk′, j , I3PEM is more effective than the 4m + 1 scheme.

⎡ 1 ρ12 1 ⎢ρ Px = ⎢ 21 ⋮ ⋮ ⎢ ⎣ ρn1 ρn2

2 ρ12 σx1 σx2 ⎡ σx 1 ⎢ρ σ σ σx22 C x = ⎢ 21 x2 x1 ⎢ ⋮ ⋮ ⎢ ρn1 σxn σx1 ρn2 σxn σx2 ⎣

3.1. PEM for the PLF with correlated input variables The main idea of PLF calculation for power system using PEM is to select m points for each input random variable according to the probability distribution characteristics. And then the deterministic load flow calculations combing with the selected points and mean of other input random variables are performed. The deterministic load flow equation of the power system is

(19)

1. Carry out the Cholesky decomposition of C x as C x = LLT , and let B=L−1 . 2. Transform the first four central moments of the correlated input random variables x to the new space as

n

(15)

μq = L−1μ x

where Pi and Qi are the injected active and reactive power of bus i respectively; Ui and δi are the voltage amplitude and phase angle of bus i respectively; δij = δi − δj ; Gij and Bij are the real and imaginary parts of the bus admittance matrix respectively; n is the number of system nodes. Therefore, the output variable Z, such as the bus voltage and the branch transmission power to be solved, can be expressed as a nonlinear function of power (including active power P and reactive power Q), which can be expressed as Z = H (P, Q) .When the PEMs are used to calculate the PLF, one deterministic load flow calculation is performed at each point (μ1 , μ 2 , ...,xk, i , ...,μn ) . Therefore, the solution of load flow equation is

Cq = L−1C x (L−1)T = I n

−1 3 ) λ xr ,3 σx3r λ qk,3 = ∑ (Lkr r=1 n

−1 4 ) λ xr ,4 σx4r λ qk,4 = ∑ (Lkr r=1

(20)

3. Compute the standard location ξk, i and the corresponding weight wk, i as explained in Section 2. 4. From the 2PEM, 3PEM and I3PEM vectors in the newly transformed space q , the transformed points are constructed in the form (μq1 , μq2 , …, qk, i , …, μqn ) . 5. Transform the points defined in step 4 back into the original space by using the inverse transformation x=B - 1q . 6. Solve one deterministic power flow problem described in (16) for each one of the points resulting from step 5. This step generates the solution vectors Z(k , i) . 7. The whole result for moment E (Z j ) can be obtained in (17). The means and standard deviations of output random variables also be computed. And the approximation/expansion method desired over these moments can be applied to approximate the PDF and CDF of output random variables, Z .

(16)

where Z (k , i) is the output random variable associated with the location of i-th point of the random variable xk ; H is the nonlinear relationship between the input and output random variables in the load flow equation. The total calculation number of deterministic power flow depends on the number of estimate points of PEM. Combined with the 2PEM, 3PEM and the proposed I3PEM in Section 2, the solutions of the PLF calculation for the power system can be obtained, as follows.

E (Z ) ≅ E (Z ) + wk, i Z (k , i) E (Z j ) ≅ E (Z j ) + wk, i Z (k , i) j

⋯ ρ1n σx1 σxn ⎤ ⋯ ρ2n σx2 σxn ⎥ ⎥ ⎥ ⋱ ⋮ ⎥ ⋯ σx2n ⎦

where σx1, σx2, ⋯, σxn are the standard deviations of the random variables x1, x2 , …, x n respectively. The procedure using PEMs to solve the PLF with correlated input random variables is as follows [36,37]:

⎧ Pi = Ui ∑ Uj (Gij cos δij + Bij sin δij ), i = 1, 2, …, n ⎪ j=1

Z (k , i) = H (μ1 , μ 2 , ...,xk, i , ...,μn )

(18)

where T denotes the transpose of a matrix, ρij is the correlation between ith and jth random variable. Then the variance-covariance matrix can be obtained form the correlation matrix and known standard deviations.

3. Probabilistic load flow

n ⎨ ⎪Qi = Ui ∑ Uj (Gij sin δij − Bij cos δij ), i = 1, 2, …, n j=1 ⎩

⋯ ρ1n ⎤ ⋯ ρ2n ⎥ ⋱ ⋮ ⎥ ⎥ ⋯ 1 ⎦

3.2. MCS for the PLF with correlated input variables (17) Given that the MCS is used as a benchmark methodology, the required N samples must be created for each input random variables in the load flow. Since the PEMs consider the correlation among input random variables, the samples of input random variables generated for the MCS must also be correlated. Based on the Normal to anything (NORTA) process [38], the procedure applied to generate correlated samples for MCS is described as follows [36,37]:

By using (17), the estimated origin moments of Z are used to calculate the wanted statistical information of the output random variables. In this paper, Cornish-Fisher series expansion [26,27] is used to calculate the probability distribution function of the output variables based on the mean, standard deviation, and the first fourth order origin moments. The previously described is the step process of uncorrelated input random variables. In the presence of correlation among input random variables, such as loads and wind generation, the well-known rotational transformation method is required [14,35]. This precautionary step

1. Let a vector of n correlated and arbitrarily distributed random variables x = (x1, x2, …, x n )T with a correlation matrix as in (18). 4

Electrical Power and Energy Systems 117 (2020) 105618

Y. Che, et al.

Since the NORTA process is based on normal distribution to generate the random correlation values, the matrix Px must be adjusted to Py in normal distribution space. A multiplication factor G is possible established to relate each correlation coefficient ρij′ in Py to its corresponding counterpart ρij in Px .

ρij′ = G (ρij ) ρij

Table 1 Mean and standard deviation results for IEEE 14-bus system (CV = 0.15, Selected Values). δ10 [p.u.]

P1 [p.u.]

Q3 [p.u.]

P2,3 [p.u.]

Q2,3 [p.u.]

2PEM

1.0509 0.0038 1.0509 0.0038 1.0509 0.0038 1.0509 0.0039

−0.2636 0.0168 −0.2636 0.0168 −0.2636 0.0168 −0.2636 0.0167

2.3255 0.1946 2.3255 0.1945 2.3254 0.1944 2.3250 0.1928

0.2523 0.0886 0.2523 0.0885 0.2523 0.0885 0.2521 0.0875

0.7328 0.0845 0.7328 0.0845 0.7328 0.0845 0.7326 0.0836

0.0363 0.0085 0.0363 0.0083 0.0363 0.0084 0.0363 0.0083

3PEM I3PEM

2. Cholesky decomposition is applied to Py in order to obtain a lowertriangular matrix, so that LLT = Py . And let B=L−1 , such that w = By , with w and y being a vector of uncorrelated and correlated standard normal variables respectively. 3. Generate the sample ws = (w1, s, w2, s, …, wn, s )T of independent standard normal variables. 4. By reversing the orthogonal transformation, each vector ys with correlated standard normal samples is obtained from the sample ws in step 3 above, that is, ys = B−1ws . 5. By reverting the normal transformation, the sample x s of the original vector x with correlated input random variables can be obtained, that is xs = Fx−s1 (Φ(ys )) . 6. Solve one deterministic load flow calculation for the vector x s obtained in Step 5. 7. Steps 3–6 are repeated and performed until a sufficient number N of simulation. 8. The means and standard deviations of output random variables are computed and other statistical information interested also be obtained.

MCS

4.1. IEEE 14-bus system 4.1.1. Case study with normal distribution inputs in IEEE 14-bus system It is assumed that the active power in IEEE 14-bus system obey the Normal distribution, in which the load adopts constant power factor model. The results of the MCS simulation N = 10,000 times are used as the reference. When the CV is selected as 0.15 for IEEE 14-bus system, the results of mean μ and standard deviation σ for some selected output variables are shown in Table 1. The selected output results of load flow are the amplitude and angle of voltage at bus no.10, the injection power at bus no.1, the injection reactive power at bus no.3, and the active power and reactive power in the line no.2–3. The mean and standard deviation of these PLF output results obtained by 2PEM, 3PEM, I3PEM are very close compared to those obtained by MCS. Table 2 shows the average error index of IEEE 14-bus system for three different PEM schemes. The average error is the mean value ε¯μ in mean error and the mean ε¯σ in standard deviation error. The smaller ε¯μ and ε¯σ index, the better the performance of the method. The output variables as shown in Table 2 are amplitude V and angle δ of voltage, injected active power P and reactive power Q, line active power Pij and reactive power Qij. It can be seen that the overall performance of 3PEM is better than that of 2PEM. The ε¯μ and ε¯σ index obtained by I3PEM are obviously smaller than 3PEM. Therefore, the accuracy of PLF results performed by I3PEM on IEEE 14-bus systems are better than those of 3PEM and 2PEM. It should be noted that the advantages of I3PEM are not obviously seen from Table 1, however it can be clearly from Table 2 that the accuracy of PLF results by I3PEM are better than that of the other two (2PEM and 3PEM). Figs. 1 and 2 show that the mean error ε¯μ in mean and the mean error ε¯σ in standard deviation of three kinds of PEMs increase with CV growth in the case of five CVs (0.05, 0.10, 0.15, 0.20

1. Compute the eigenvalues d1, d2, ⋯, dn and the corresponding eigenvectors v1, v2, ⋯, vn of the matrix cov (U) . 2. Generate the n-vector Γ of independent Gaussian random variables represented by zero mean values and standard variance equal to d 1 , d2 , ⋯, dn . 3. Generate the vector U as:

U = μ (U) + (v1, v2, ⋯, vn )Γ

(22)

3.3. Measures to assess the performances of methods To compare the performance of different PEMs scheme, the following error indexes for each output random variable are defined N

∑i r |μ MCS − μPEM |/|μ MCS | Nr

× 100 [%]

(23)

N

∑i r |σMCS − σPEM |/|σMCS | Nr

× 100 [%]

μ σ μ σ μ σ μ σ

variables considering Normal or Non-normal distributed, correlated or not. Define the coefficient of variation (CV) as the ratio between the standard deviation and the mean. The mean of the input random variables associated with the location of estimate points are the load data given by [40]. The PLF calculation is done on a computer with Intel(R) Core(TM) i3 2.53 GHz and RAM of 4 GB, using MATPOWER [40] in MATLAB environment.

In order to compare with the adjustment procedure for taking into account the correlation, such as the above Cholesky decomposition, the following procedure is applied to generate the corresponding vector U of inputs random variables presented by the mean value μ (U) and the covariance matrix cov (U) [32].

ε¯σ =

V10 [p.u.]

(21)

The value of G for each conversion ρij → ρij′ can be obtained from the expression provided in [39]. If x i and x j are normally distributed, the value of G is equal to 1.

ε¯μ =

IEEE 14-bus system μ and σ Results

(24)

Table 2 Average error of IEEE 14-bus system (CV = 0.15).

where μ MCS and σMCS are used as reference values, which are respectively the mean and standard deviation calculated by the MCS; similarly, μPEM and σPEM are the mean and standard deviation calculated by the given PEM, respectively; Nr is the number of random variables. 4. Case studies and discussion In this section, the IEEE 14-bus and 118-bus system are employed to verify the proposed PEMs in the case of different type of random 5

IEEE 14-bus system

V

δ

P

Q

Pij

Qij

2PEM

ε¯μ [%]

0.0016

0.0127

0.0223

0.0427

0.0839

0.1128

3PEM

ε¯σ [%] ε¯μ [%]

0.4215 0.0016

0.3673 0.0127

0.9109 0.0223

0.8949 0.0426

0.6569 0.0839

0.6918 0.1126

I3PEM

ε¯σ [%] ε¯μ [%]

0.5075 0.0015

0.3319 0.0114

0.8596 0.0193

0.8020 0.0352

0.6505 0.0815

0.7332 0.1025

ε¯σ [%]

0.4299

0.3369

0.8489

0.8085

0.6301

0.6785

Electrical Power and Energy Systems 117 (2020) 105618

Y. Che, et al.

14-bus system obey the Normal distribution. And the load adopts constant power factor model. Case studies with correlated inputs in IEEE 14-bus system have two scenarios as follow: Scenario A: consistent correlation. The consistent correlation coefficient of all correlated loads is ρ = 0.75. Scenario B: inconsistent correlation. The whole system of IEEE14bus system is divided into 2 areas. Area 1 includes node no.1–6. Area 2 includes node no.7–14. Loads in the same area are correlated, and loads that are not in the unified area are independent of each other. Arranged by the size of the load nodes, the correlation coefficient matrix is

Fig. 1. Mean error in mean Vs coefficient of variation of IEEE 14-bus system.

⎡1.0 ⎢ 0.5 P1 = ⎢ 0.8 ⎢ 0.7 ⎢ ⎣ 0.6

0.5 1.0 0.7 0.6 0.8

0.8 0.7 1.0 0.5 0.6

0.7 0.6 0.5 1.0 0.9

0.6 ⎤ 0.8 ⎥ 0.6 ⎥ 0.9 ⎥ 1.0 ⎥ ⎦

⎡1.0 ⎢ 0.5 0.8 P2 = ⎢ ⎢ 0.7 ⎢ 0.6 ⎢ 0.4 ⎣

0.5 1.0 0.7 0.6 0.8 0.5

0.8 0.7 1.0 0.5 0.6 0.5

0.7 0.6 0.5 1.0 0.9 0.6

0.6 0.8 0.6 0.9 1.0 0.6

(25)

0.4 ⎤ 0.5 ⎥ 0.5 ⎥ 0.6 ⎥ 0.6 ⎥ 1.0 ⎥ ⎦

(26)

Table 3 shows the average error indexes in the presence of light, strong and inconsistent correlation among the input random variables when the CV is equal to 0.15. It is clear that the presence of correlation increases the average error of all PEMs schemes. In particular, the average error of the standard deviation is assumed to be non-negligible values. Therefore, the suitable adjustment procedure is necessary to take into account to the correlation. The procedure using PEMs with rotational transformation to solve the PLF with correlated input random variables can be seen in Section 3.1. The procedure applied to generate correlated samples for MCS using NORTA method can be seen in Section 3.2. With same as Scenario A and Scenario B in Table 3 when the CV is equal to 0.15, Table 4 shows the average error indexes for correlated input random variables applying the procedure able to considering correlation. Compared the results in Table 4 with those in Table 3, it is evident that the procedure considering correlation improves the accuracy of the results. Under Scenario A and Scenario B, Figs. 3 and 5 show that the mean error ε¯μ in mean of three kinds of PEMs increased with CV growth respectively, and Figs. 4 and 6 show that the mean error ε¯σ in standard deviation of three kinds of PEMs increased with CV growth respectively. It can be seen that ε¯μ and ε¯σ of I3PEM are all smaller than the other two. Hence, the accuracy of PLF results performed by I3PEM on IEEE 14-bus systems for corelated input variables applying the procedure able to consider correlation are better than those of 3PEM and 2PEM.

Fig. 2. Mean error in standard deviation Vs coefficient of variation of IEEE 14bus system.

and 0.25), respectively. It can be seen that ε¯μ and ε¯σ of I3PEM are all smaller than the other two. 4.1.2. Case study with correlated inputs in IEEE 14-bus system Based on the above case, it is assumed that the active power in IEEE Table 3 Average error of IEEE 14-bus system for corelated input variables (CV = 0.15). IEEE 14-bus system

V

δ

P

Q

Pij

Qij

Scenario A 2PEM

ε¯μ [%]

0.0040

0.0477

0.0661

0.7298

0.0581

0.8994

3PEM

ε¯σ [%] ε¯μ [%]

42.1831 0.0040

56.1187 0.0478

49.5702 0.0661

41.5098 0.7302

37.1878 0.0581

39.0798 0.8999

I3PEM

ε¯σ [%] ε¯μ [%]

42.2291 0.0040

56.1361 0.0478

49.5957 0.0661

41.5970 0.7300

37.1920 0.0581

39.2507 0.8997

ε¯σ [%]

42.2061

56.1274

49.5829

41.5534

37.1899

39.1650

Scenario B 2PEM

ε¯μ [%]

0.0010

0.0169

0.0034

0.2595

0.0521

0.3408

3PEM

ε¯σ [%] ε¯μ [%]

30.7784 0.0010

36.9337 0.0169

32.6279 0.0033

26.9497 0.2599

30.5386 0.0521

25.2971 0.3411

I3PEM

ε¯σ [%] ε¯μ [%]

30.8376 0.0010

36.9587 0.0169

32.6620 0.0033

27.0618 0.2597

30.5441 0.0521

25.5205 0.3410

ε¯σ [%]

30.8080

36.9462

32.6449

27.0057

30.5413

25.4084

6

Electrical Power and Energy Systems 117 (2020) 105618

Y. Che, et al.

Table 4 Average error of IEEE 14-bus system for correlated input variables applying the procedure able to consider correlation. (CV = 0.15). IEEE 14-bus system

V

δ

P

Q

Pij

Qij

Scenario A 2PEM ε¯μ [%]

0.0051

0.1556

0.1685

0.3268

0.1514

0.4374

3PEM

ε¯σ [%] ε¯μ [%]

0.1381 0.0051

0.1231 0.1559

0.1175 0.1687

0.2463 0.3274

0.1954 0.1515

1.0890 0.4394

I3PEM

ε¯σ [%] ε¯μ [%]

0.1968 0.0045

0.2288 0.1393

0.1875 0.1532

0.1826 0.2846

0.2009 0.1450

0.5019 0.3921

ε¯σ [%]

0.1476

0.2076

0.1797

0.1556

0.2203

0.5268

Scenario B 2PEM

ε¯μ [%]

0.0066

0.1865

0.1848

0.3843

0.1687

0.5406

3PEM

ε¯σ [%] ε¯μ [%]

0.1522 0.0066

0.3214 0.1866

0.2229 0.1848

0.2483 0.3845

0.2528 0.1687

0.3617 0.5413

I3PEM

ε¯σ [%] ε¯μ [%]

0.1735 0.0071

0.2806 0.1797

0.1913 0.1818

0.1504 0.3197

0.2488 0.2083

0.5536 0.5113

ε¯σ [%]

0.1300

0.2796

0.1810

0.1905

0.2744

0.4001

Fig. 5. Mean error in mean Vs coefficient of variation of IEEE 14-bus system in Scenario B.

Fig. 3. Mean error in mean Vs coefficient of variation of IEEE 14-bus system in Scenario A.

Fig. 6. Mean error in standard deviation Vs coefficient of variation of IEEE 14bus system in Scenario B. Table 5 Average error of IEEE 14-bus system for bimodal PDF’s (CV = 0.15). IEEE 14-bus system

V

δ

P

Q

Pij

Qij

2PEM

ε¯μ [%]

0.0007

0.0387

0.0520

0.0766

0.1479

0.1719

3PEM

ε¯σ [%] ε¯μ [%]

6.4984 0.0007

7.1785 0.0387

1.5188 0.0519

4.3807 0.0769

9.5385 0.1479

6.0236 0.1721

I3PEM

ε¯σ [%] ε¯μ [%]

6.5820 0.0007

7.2135 0.0387

1.5675 0.0519

4.4861 0.0767

9.5419 0.1479

6.0481 0.1720

ε¯σ [%]

6.5402

7.1960

1.5432

4.4334

9.5402

6.0358

bi-modal PDF’s that denote the active and reactive load powers are located at bus no.9 and 13. Table 5 shows the average error indexes in the presence of the bi-modal PDF’s. With the same to [33], the bi-modal PDF’s are obtained by distributing the active and reactive load power values around two points, and the standard deviations are equal to 15%. Fig. 7 shows the active load powers at bus nos 9 and 13 respectively. Figs. 8–10 shows the PDFs of voltage at bus no. 9, 10 and 14 respectively. Figs. 11 and 12 shows the PDFs of active and reactive load flow in line 9–14 respectively. Fig. 13 show the PDF of power loss. Also in this case study, the results obtained with I3PEM present a better approximation than the 2PEM and 3PEM.

Fig. 4. Mean error in standard deviation Vs coefficient of variation of IEEE 14bus system in Scenario A.

4.1.3. Case study with non-normal distribution inputs in IEEE 14-bus system In order to verify the performance of the improved PEMs in which the PDF’s of the real and reactive load powers are not all Gaussian, the 7

Electrical Power and Energy Systems 117 (2020) 105618

Y. Che, et al.

Fig. 7. PDF of active load powers at bus nos 9 and 13. Fig. 10. PDF of voltage at no.14.

Fig. 8. PDF of voltage at no.9. Fig. 11. PDF of active load flow in line no.9–14.

Fig. 9. PDF of voltage at no.10. Fig. 12. PDF of reactive load flow in line no.9–14.

4.2. IEEE 118-bus system 4.2.1. Case study with normal distribution inputs in IEEE 118-bus system Same as the case study with normal distributed inputs in the IEEE 14 bus system, Table 6 is the average error indexes of three kinds of PEMs for IEEE 118-bus system. The CV of IEEE118-bus system is selected as

0.15, and the results of the MCS simulation 10,000 times are used as the reference. Although the ε¯σ index obtained by I3PEM are larger than those of 3PEM, the ε¯μ index obtained by I3PEM are slightly smaller than those of 3PEM. In the case of five CVs (0.05, 0.10, 0.15, 0.20 and 0.25), 8

Electrical Power and Energy Systems 117 (2020) 105618

Y. Che, et al.

Fig. 15. Mean error in standard deviation Vs coefficient of variation of IEEE 118-bus system.

Fig. 13. PDF of power loss. Table 6 Average error of IEEE 118-bus system (CV = 0.15). IEEE 118

V

δ

P

Q

Pij

Qij

2PEM

ε¯μ [%]

0.0011

0.0573

0.1139

1.0440

0.4332

0.2248

3PEM

ε¯σ [%] ε¯μ [%]

2.1646 0.0011

0.6042 0.0581

1.0358 0.1161

3.6389 1.0386

0.6148 0.4352

2.1605 0.2247

I3PEM

ε¯σ [%] ε¯μ [%]

0.9852 0.0011

0.4241 0.0580

0.7161 0.1155

0.7535 1.0390

0.5472 0.4343

0.9618 0.2244

ε¯σ [%]

1.1956

0.4377

0.8772

1.9676

0.5503

1.2245

Fig. 16. Mean error in mean Vs correlation coefficient of IEEE 118-bus system with correlated inputs.

Fig. 14. Mean error in mean Vs coefficient of variation of IEEE 118-bus system.

Figs. 14 and 15 show that the mean error ε¯μ in mean and the mean error ε¯σ in standard deviation of IEEE 118-bus system with three kinds of PEMs increase with CV growth respectively. As the number of input random variables increases, the results of ε¯σ (the mean value in standard deviation error) with I3PEM is not as good as than those of 3PEM, while the results of ε¯μ (the mean value in mean error) with I3PEM is slightly better than 3PEM and 2PEM. Fig. 17. Mean error in standard deviation Vs correlation coefficient of IEEE 118-bus system with correlated inputs.

4.2.2. Case study with correlated inputs in IEEE 118-bus system Similar to the Scenario A of case study with correlated inputs in IEEE 14-bus system, this section shows the case study with correlated inputs in IEEE 118-bus system and correlation effect on the three kinds of PEMs scheme. In the case of different correlation coefficient 9

Electrical Power and Energy Systems 117 (2020) 105618

Y. Che, et al.

(between 0.1 and 0.8) and the CV of IEEE118-bus system is selected as 5%, Figs. 16 and 17 show the mean error ε¯μ in mean and the mean error ε¯σ in standard deviation in IEEE-118 bus system with correlated inputs respectively. Compared with the larger value of the mean error ε¯σ in standard deviation, the mean error ε¯μ in mean with three kinds of PEMs scheme are small and those difference are also small. As the correlation coefficient increases, the mean error ε¯μ in mean by I2PEM is gradually smaller than that of 3PEM and 2PEM, and the mean error ε¯σ in standard deviation by I2PEM gradually increases between 3PEM and 2PEM. 4.3. Discussion Compared with the size of IEEE 118-bus system, the size of IEEE 14bus system is relatively small. It is known from the case studies in Section 4.1 that the performance of I3PEM is better in the results of PLF calculation with Normal distribution inputs, correlated inputs and bimodal distribution inputs. The value of input random variable number determines the size of one system. In (14) the standard locations of the new pair of estimate points for I3PEM are related to n , and the new locations xk′, i (k = 1, 2, …n, i = 1, 2) are away from the mean μk by the scale relationship n when n increases. The Chebyshev inequality gives the divergence degree of the deviation of the random variable from the mean [35]. From the standardization technique with z-score of Normal distribution [41], the z-score equation is:

z=

xk − μk σk

Fig. 19. Comparison of mean error in standard deviation of IEEE 14-bus and 118-bus system by I3PEM.

2PEM. Figs. 18 and 19 shows the ε¯μ and ε¯σ results obtained by I3PEM under different CV values for IEEE 14-bus and 118-bus system respectively. From the above analysis, the ε¯μ and ε¯σ results of IEEE 118-bus system obtained by I3PEM are greatly increased as the CV increases. However, those mean error results of IEEE 114-bus system are not changed much with the increasing of CVs. These results show that the deterministic load flow calculated by the Newton-Raphson method may not converge and it takes more time, when the locations extracted by the PEMs locate outside the multiple standard deviation and exceed the domain of the random variables. This also reveals the shortcomings of the 2m concentration scheme. Since the newly pair of estimate points of I3PEM are added based on 3PEM, the calculation time of I3PEM will be more than that of 2PEM and 3PEM. Table 7 shows the CPU time for PLF with Normal distribution inputs of IEEE 14-bus and 118-bus systems using three kinds of PEMs and MCS simulation for 10,000 times, respectively. The CPU time is taken as the average computational time value of PLF calculations at the condition of 5 CVs (0.05, 0.10, 0.15, 0.20 and 0.25, respectively). Although 3PEM has one more calculation for deterministic power flow than 2PEM, 3PEM is faster than 2PEM. The calculation time of I3PEM is close to two times of that of 3PEM. However, the calculation time of I3PEM is obviously faster compared with the calculation time of MCS. Some studies need an equivalent system or do not require a large system, such as the impact of the traction power supply system on the grid is local [44–46]. Therefore, I3PEM is a better choice for the small-scale systems.

(27)

Comparing (27) with (1), it is found that z corresponds to the standard location. From the 3σ criterion (or 68–95-99.7 rule) [41], 68%, 95%, 99.73%, and 99.9937% of observations fall within 1, 2, 3, and 4 standard deviations of the mean in Normal distribution. Combining (7), (11) and (14), let λk,3 be equal to zero, then ξk max = n max = 4 , and thus, n max = 16. In other words, I3PEM is more suitable for the system with the number of input random variables less than 16, like the IEEE 14-node system. This result is also confirmed by other literature that many reported results indicate unacceptable accuracy of the PEMs in the case of large-scale and complex systems [29,42,43]. The number of input random variables (PQ nodes) are 99 and 11 for the IEEE 118-bus system and 14-bus system respectively. For the case studies of Normal distribution inputs and correlated inputs of IEEE 118bus system in Section 4.2, the mean error ε¯σ in standard deviation by I2PEM gradually increases between 3PEM and 2PEM, while the mean error ε¯μ in mean with I3PEM is slightly better than that of 3PEM and

5. Conclusion The improved three point-estimate method (I3PEM) scheme is proposed to estimate the probability moments of PLF in this paper. On the basis of the original 3PEM and combining with the Chebyshev inequality, only a new pair of estimate points considering the skewness of input random variables is added for additional calculations to improve the calculation accuracy of 3PEM. The case studies of IEEE 14-bus and 118-bus system are employed to verify the performance of I3PEM in the presence of different type of random variables with Normal distribution Table 7 CPU time with normal distribution inputs (in Seconds).

Fig. 18. Comparison of mean error in mean of IEEE 14-bus and 118-bus system by I3PEM. 10

CPU Time

2PEM

3PEM

I3PEM

MCS

IEEE14 IEEE118

0.1870 2.7864

0.1840 2.7714

0.3650 5.5076

82.5094 131.9674

Electrical Power and Energy Systems 117 (2020) 105618

Y. Che, et al.

inputs, correlated inputs and non-normal inputs. Considering five kinds of different CVs, the results of the MCS method are taken as benchmark. And the results of I3PEM are compared with those obtained by the existing 2PEM and 3PEM. The following conclusions were obtained:

[14] [15] [16]

1. The additional calculations by the newly added pair of estimate points of I3PEM is equivalent to the consequent of the higher the number of points used by a PEM, which yields that the more accurate estimation of the statistical moments of the output random variables is achieved. However, I3PEM is more concise and efficient because it avoids the calculation of higher order moments greater than fourth order and the possible non-real number solutions, or even an infinite set of solutions. 2. The performance of I3PEM is better than that of 3PEM and 2PEM in the results of PLF calculation with Normal distribution inputs, correlated inputs and bi-modal distribution inputs, like the IEEE 14-bus system. However, with increasing the scale size of system (or the number of input random variables), the mean error in standard deviation by I3PEM gradually increases between 3PEM and 2PEM, and the mean error in mean with I3PEM is slightly better than that of 3PEM and 2PEM, like the IEEE 118-bus system. 3. The standard locations of 2PEM (or 2m concentration scheme) are related to the square root of the number of input variable. The deterministic load flow calculated by the Newton-Raphson method may not converge and it takes more time, when the locations extracted by the PEMs locate outside the multiple standard deviation and exceed the domain of the random variables.

[17]

[18]

[19]

[20]

[21] [22]

[23]

[24]

[25]

[26]

Declaration of Competing Interest

[27]

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

[28]

[29]

Acknowledgement [30]

This work was supported by Innovative Ability Enhancement Project of Gansu Provincial Higher Education (No. 2019B-049).

[31] [32]

References

[33] [34]

[1] Aien Morteza, Rashidinejad M, Firuz-Abad MF. Probabilistic optimal power flow in correlated hybrid wind-PV power systems: a review and a new approach. Renew Sustain Energy Rev 2015;41(1):1437–46. [2] Ruiz-Rodriguez FJ, Hernández JC, Jurado F. Voltage behaviour in radial distribution systems under the uncertainties of photovoltaic systems and electric vehicle charging loads. Int Trans Electr Energy Syst 2018;28(2):e2490. [3] Borkowska B. Probabilistic load flow. Power Appar Syst IEEE Trans 1974;93(3):752–9. [4] Aien Morteza, Rashidinejad M, Fotuhi-Firuzabad M. On possibilistic and probabilistic uncertainty assessment of power flow problem: a review and a new approach. Renew Sustain Energy Rev 2014;37(3):883–95. [5] Ruiz-Rodriguez FJ, Hernández JC, Jurado F. Voltage unbalance assessment in secondary radial distribution networks with single-phase photovoltaic systems. Int J Electr Power Energy Syst 2015;64:646–54. [6] Nikmehr N, Najafi Ravadanegh S. Heuristic probabilistic power flow algorithm for microgrids operation and planning. Gener Transmiss Distrib IET 2015;9(11):985–95. [7] Prusty B Rajanarayan, Jena D. An over-limit risk assessment of PV integrated power system using probabilistic load flow based on multi-time instant uncertainty modeling. Renew Energy 2018;116:367–83. [8] Tung Y-K, Yen B-C. Hydrosystems engineering uncertainty analysis. New York: McGraw-Hill; 2005. [9] Liu Y, Gao S, Yang S, et al. Review on algorithms for probabilistic load flow in power system. Autom Electr Power Syst 2014;38(23):127–35. [10] Liu Jun S. Monte Carlo strategies in scientific computing. Springer; 2001. [11] Kroese DP, Taimre T, Botev ZI. Handbook of Monte Carlo methods. John Wiley & Sons; 2013. [12] Huang H, Chung CY, Chan KW, et al. Quasi-Monte Carlo based probabilistic small signal stability analysis for power systems with plug-in electric vehicle and wind power integration. IEEE Trans Power Syst 2013;28(3):3335–43. [13] Chen Yan, Wen J, Cheng S. Probabilistic load flow method based on Nataf

[35] [36] [37]

[38]

[39] [40]

[41] [42] [43]

[44] [45] [46]

11

transformation and Latin hypercube sampling. IEEE Trans Sustain Energy 2013;4(2):294–301. Xu Xiaoyuan, Yan Z. Probabilistic load flow evaluation considering correlated input random variables. Int Trans Electr Energy Syst 2016;26(3):555–72. Allan RN, Al-Shakarchi MRG. Linear dependence between nodal powers in probabilistic AC load flow. Proc Inst Electr Eng 1977;124.6:529–34. Allan RN, Leite dSAM, Burchett RC. Evaluation methods and accuracy in probabilistic load flow solutions. Power Apparat Syst IEEE Trans PAS 1981;100(5):2539–46. Zhang Pei, Lee ST. Probabilistic load flow computation using the method of combined cumulants and Gram-Charlier expansion. IEEE Trans Power Syst 2004;19(1):676–82. Ruiz-Rodriguez FJ, Hernández JC, Jurado F. Probabilistic load-flow analysis of biomass-fuelled gas engines with electrical vehicles in distribution systems. Energies 2017;10(10):1–23. Hernández JC, Ruiz-Rodriguez FJ, Jurado F. Modelling and assessment of the combined technical impact of electric vehicles and photovoltaic generation in radial distribution systems. Energy 2017;141(1):316–32. Tao Wang, Chun Wang, Chenghao Li. Probabilistic load flow calculation based on Copula function and Rosenblatt transformation considering correlation among input variables. Power Syst Protect Control 2018;46(21):18–24. Hong HP. An efficient point estimate method for probabilistic analysis. Reliab Eng Syst Saf 1998;59:261–7. Aien Morteza, Fotuhi-Firuzabad M, Aminifar F. Probabilistic load flow in correlated uncertain environment using unscented transformation. IEEE Trans Power Syst 2012;27(4):2233–41. Hong YY, Lin FJ, Yu TH. Taguchi method-based probabilistic load flow studies considering uncertain renewables and loads. IET Renew Power Gener 2016;10(2):221–7. Carpinelli G, Rizzo R, Caramia P, et al. Taguchi’s method for probabilistic threephase power flow of unbalanced distribution systems with correlated wind and photovoltaic generation systems. Renew Energy 2018;117:227–41. Prusty BR, Jena D. A critical review on probabilistic load flow studies in uncertainty constrained power systems with photovoltaic generation and a new approach. Renew Sustain Energy Rev 2017;69:1286–302. Usaola J. Probabilistic load flow with wind production uncertainty using cumulants and Cornish-Fisher expansion. Int J Electr Power Energy Syst 2009;31(9):474–81. Ruiz-Rodriguez FJ, Hernandez JC, Jurado F. Probabilistic load flow for photovoltaic distributed generation using the Cornish-Fisher expansion. Electr Pow Syst Res. 2012;89:129–38. Zhang Z, Li GY, Wei JQ. Probabilistic evaluation of voltage quality in distribution networks considering the stochastic characteristic of distributed generators. Proc CSEE 2013;33.13:150–6. Mohammadi M, Shayegani A, Adaminejad H. A new approach of point estimate method for probabilistic load flow. Int J Electr Power Energy Syst 2013;51(10):54–60. Su CL. Probabilistic load-flow computation using point estimate method. IEEE Trans Power Syst 2005;20(4):1843–51. Morales Juan M, Perez-Ruiz J. Point estimate schemes to solve the probabilistic power flow. IEEE Trans Power Syst 2007;22(4):1594–601. Caramia P, Carpinelli G, Varilone P. Point estimate schemes for probabilistic threephase load flow. Electr Power Syst Res 2010;80(2):168–75. Rosenblueth E. Point estimates for probability moments. PNAS 1975;72.10:3812–4. Rice John A. Mathematical statistics and data analysis. 3rd ed. Brooks/Cole Thomson Learning; 2006. Edgar NM, Nicolaos C, Poul S. High dimensional dependence in power systems: a review. Renew Sustain Energy Rev 2018;94:197–213. Morales JM, Baringo L, Conejo AJ, et al. Probabilistic power flow with correlated wind sources. IET Gener Transm Distrib 2010;4(5):641–51. Delgado C, Domínguez-Navarro JA. Point estimate method for probabilistic load flow of an unbalanced power distribution system with correlated wind and solar sources. Int J Electr Power Energy Syst 2014;61:267–78. Cario MC, Nelson BL. Modeling and generating random vectors with arbitrary marginal distributions and correlation matrix. Evanston, Illinois: Department of Industrial Engineering and Management Sciences, Northwestern University; 1997. Tech. Rep. Der Kiureghian A, Asce M, Liu PL. Structural reliability under incomplete probability information. J Eng Mech 1986;112(1):85–104. Zimmerman RD, Murillo-Sanchez C, Matpower Deqing G. a matlab power system simulation package. 3.0.0ed. NY, USA: Cornell University Power System Engineering Research Center (PSERC); 2005. Diez David M, Barr Christopher D, Çetinkaya-Rundel Mine. OpenIntro statistics. 3rd ed. CreateSpace Independent Publishing Platform; 2015. Mohammadi M, Basirat H, Kargarian A. Nonparametric probabilistic load flow with saddle point approximation. IEEE Trans Smart Grid 2018;9(5):4796–808. Vahid-Pakdel MJ, Mohammadi-Ivatloo B. Probabilistic assessment of wind turbine impact on distribution networks using linearized power flow formulation. Electr Power Syst Res 2018;162:109–17. Yulong Che X, Wang X Lv, et al. Impact analysis of traction loads on power grid based on probabilistic three-phases load flow. J Eng 2019;16:2306–11. Haji Moosa Moghimi, Xu W. Online determination of external network models using synchronized phasor data. IEEE Trans Smart Grid 2018;9(2):635–43. Mohammad Mohammadi. Probabilistic harmonic load flow using fast point estimate method. IET Gener Transm Distrib 2015;9(13):1790–9.