Probabilistic load flow with detailed wind generator models considering correlated wind generation and correlated loads

Renewable Energy 94 (2016) 96e105

Contents lists available at ScienceDirect

Renewable Energy journal homepage: www.elsevier.com/locate/renene

Probabilistic load flow with detailed wind generator models considering correlated wind generation and correlated loads Neeraj Gupta Department of Electrical Engineering, National Institute of Technology, Hamirpur, Himachal Pradesh, India

a r t i c l e i n f o

a b s t r a c t

Article history: Received 6 July 2015 Received in revised form 19 January 2016 Accepted 4 March 2016

The enhancement in the penetration of intermittent generation necessitates the need to include uncertain behaviour in the conventional power flow programs. In this paper, four different wind generation models have been incorporated in probabilistic load flow for calculating the probability distribution of the reactive power consumed by the wind generators for three different scenarios; i) uncorrelated wind and uncorrelated loads ii) uncorrelated wind and correlated loads and iii) correlated wind and correlated loads The above mentioned scenarios have been implemented in probabilistic load flow using point estimate method in the IEEE-118 bus test system and accuracy of the results have been validated by comparing these results with those obtained by Monte Carlo simulation studies. © 2016 Elsevier Ltd. All rights reserved.

Keywords: Probabilistic load flow Point estimate method Wind generators

1. Introduction Presently, to address the twin concerns of global warming and depleting fossil fuels, a lot of importance is being given to generation of power from renewable energy sources (RES). Among various RES, wind generating system (WGS) has already reached a state of mature technology and as a result, significant number of WGS has already been installed around the world. Further, in almost every country of the world, efforts are underway to exploit the full power generation potential of WGS. Now because of random and wide variation of wind velocity, the power output from a WGS is intermittent and of fluctuating nature. When this fluctuating power is injected into the grid, it causes variations in bus voltages and line power flows of transmission system. These variations are going to be quite significant in the future (if not already) because of significant and increasing penetration of WGS in the grid. Therefore, for successful integration of WGS in the grid, these possible variations need to be properly analysed, estimated and quantified. This can be achieved through load flow analysis of the grid in the presence of uncertain power generation from WGS and towards this goal, various probabilistic load flow (PLF) methods have already been suggested in the literature [1e6]. Out of the above works, Probability Density Function (PDF) and Cumulative Distribution Function (CDF) of the quantities of interest have been calculated in Refs. [1e5], while in Ref. [6] only the mean

E-mail addresses: [email protected], [email protected]. http://dx.doi.org/10.1016/j.renene.2016.03.030 0960-1481/© 2016 Elsevier Ltd. All rights reserved.

and the standard deviations (of the quantities of interest) have been calculated using an extended Point estimate method (PEM). Further, in all these works the obtained results from PLF have also been compared with the results obtained by Monte-Carlo simulation (MCS) studies. In Refs. [5], Fourier Transform based convolution using DC power flow method has been used to find out the PDF and CDF of the variables of interest. On the other hand, in Refs. [1e4] initially the moments and cummulants of the variables of interest have been calculated using a suitable technique (such as PEM [1,2] and weighted sum of cummulants of input variables [3,4]) and subsequently the PDF and CDF have been computed using an appropriate series expansion (such as Cornish-Fisher [1,2] and Gram-Charlier [3,4]). In Refs. [5e7], different PEM based methods have been proposed for probabilistic analysis of a power system in the presence of WTGs. In Refs. [5], PEM along with Nataf transformation has been used while in Refs. [6] and [7], discrete PEM and extended PEM have been employed respectively. Although, all the above works have considered load uncertainties, correlation in loads has been considered in Refs. [2] and [6] only. For considering the uncertainties in WGS, basically two approached have been used in the literature. In the first approach, a WGS has been modeled as an uncertain real power injection using Beta distribution [1e3]. In the second approach the uncertainty in the wind speed has been considered using Weibull distribution [4e6] and Rayleigh distribution [8] and subsequently the corresponding injected real power by the WGS has been calculated using the speed-power relationship of the wind generator. Further, in Refs. [1,2,8,7] the correlation among the wind generator has also

N. Gupta / Renewable Energy 94 (2016) 96e105

been considered. Also, in some of the works such as [1,5] the uncertainty of the conventional generator has also been taken into account. From the above discussion it can be seen that in the literature, almost all types of uncertainties have been considered in PLF with WGS. However, in all of the above works, a WGS has been represented as an uncertain injected real power only and no detailed model of the WGS has been considered in these works. As a result it is not possible to calculate the PDF of the reactive power consumed by the WGS by using the approaches described above. Now, it is well known that the reactive power consumed by an individual WGS is quite significant. Therefore, with increasing penetration of WGS in the grid, the aggregate amount of reactive power consumed by all the wind generators is also going to be quite substantial. Thus, sooner or later it will be quite necessary to determine the PDF of this reactive power consumption such that adequate reactive power planning strategy can be adopted for successful integration of WGS into the grid. To address the above issue, in this paper detailed wind generator models [9] are considered to carry out PLF with wind generators. The basic objective is to calculate the PDF/CDF of the reactive power consumed by WGS in the presence of uncertain loads and uncertain wind power generation. To compute the moments of variables of interest, PEM [10,11] has been used in this work and subsequently the CDF has been determined by using the CornisheFisher expansion [1]. Further, the correlation among the loads, as well as among the wind generators has also been considered in this paper. However, it has been assumed in this work that there is no correlation between the loads and the wind generators. This paper is organized as follows: In Section 2, the basic procedure of PEM based PLF is described. In this section, both three point and five point estimate methods are discussed. In Section 3, the different wind generator models considered in this work are described in detail. In Section 4, the procedure for considering correlation in PLF is discussed. Lastly, in Sections 5 and 6, main results and conclusions of this work are presented, respectively.

97

their corresponding weights are described below.

2.1. Three point estimate Method (3PEM) 1. Find the coefficients of skewness and kurtosis of xl using eq. (2) and (3) respectively [11];

ll;3 ¼

ll;4 ¼

i h E ðxl  ml Þ3

(2)

s3l i h E ðxl  ml Þ4

(3)

s4l

N P where E½ðxl  ml Þp  ¼ ðxl ðtÞ  ml Þp  Pr ðxl ðtÞÞ; p ¼ 3,4; N is the t¼1 number of observations for xl; xl(t) is the tth observation of xl and Pr(xl(t)) is the probability of xl(t).

2. Calculate xl;1 and xl;2 by using eq. (4). Also set xl;3 ¼ 0

xl;k

ll;3 þ ð1Þ3k  ¼ 2

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 ll;4  l2l;3 ; k ¼ 1; 2: 4

(4)

3. Obtain the three point estimates of PDF (denoted as xl,1, xl,2 and xl,3 respectively) from eq. (1). Further obtain the corresponding weighting factors wl,1, wl,2 and wl,3 from eqns. (5) and (6) below.

wl:k ¼

ð1Þ3k  ; k ¼ 1; 2 xl;k xl;1  xl;2

(5)

wl:3 ¼

1 1  n ll;4  l2

(6)

l;3

2. Point estimate based PLF The point estimate method can be used to calculate the statistical moments of a random quantity which, in turn, is a function of one or several random variables [10]. In this method, ‘h’ points on the PDF are first estimated and subsequently, from these estimated points the complete PDF is constructed. The general theory of ‘h’ point estimation method is given in Ref. [12]. However, in this presented work 3 and 5 point estimate methods have been used and the detailed procedures of PLF using these three methods are given below. For this purpose, it has been assumed that in a power system there are total ‘n’ number of random input variables. For instance, if there are ‘L’ load buses in a power system, each having both real and reactive power loads which are randomly fluctuating, then n ¼ 2L. The objective of PLF is to calculate the PDFs of bus voltage magnitudes and angles from the PDFs of these ‘n’ variables. Let the lth random variable xl (l ¼ 1,2,…n) having PDF fl be considered [11]. The PEM uses two, three or ‘h’ estimated points of xl i.e. xl,1, xl,2 or xl,h as defined in Eq. (1) to replace fl by matching the first hþ1 moments of fl.

xl;k ¼ ml þ xl;k sl

for k ¼ 1; 2:::h

(1)

In eq. (1), ml and sl are mean and standard deviation of xl respectively and xl;k can be obtained as explained in the following two sub-sections for three, five and seven point estimate methods. The procedure for estimating the points of each variable xl with

2.2. Five point estimate Method 1. Find the standard central moments as [10]:

i h E ðxl  ml Þi ll;i ¼ ; i ¼ 3; …; 2m: sil

(7)

where, m ¼ 4, for five point estimate method (5PEM). 2. Find the standard locations xl;q , where q ¼ 1,…m, by obtaining the roots of the polynomial given in eq. (8).

pðxÞ ¼ C0 þ

m X

Cj x j

(8)

j¼1

In eq. (8), Cm ¼ 1 and the coefficients C0,C1,…,Cm1 are the solutions of the system of equations shown below:

2

0 6 1 6 4 « ll;m

1 ll;3 «

ll;mþ1

ll;3 ll;4 « …

… … 1 …

2 3 3 ll;mþ1 C0 6 6 7 7 7 ll;mþ1 76 C1 7 6 ll;mþ1 7 54 « 5 ¼ 4 « 5 « ll;2m1 ll;2m Cm1 ll;m

32

(9)

98

N. Gupta / Renewable Energy 94 (2016) 96e105

3. After the calculation of standard locations xl;q , obtain xl,q from eq. (1), where, q ¼ 1,…,m with m ¼ 4 for 5PEM. The weighting factors wl,q and wm are determined by the following equations:

3 2 x l;1 wl;1 6 wl;2 7 6 x2l;1 6 7 6 6 7¼6 « « 6 7 6 m1 4 wl;m1 5 6 4 xl;1 wl;m xm 2

l;1

wm ¼ 1 

n X m X

xl;2 x2l;2 « xm1 l;2 xm l;2

… … 1 … …

3 xl;m 1 2 0 3 7 x2l;m 7 7 6 6 1 7 7 « 7 l 7 6 l;3 7 7 6 4 5 xm1 « 5 l;m m ll;m xl;m

wl;k

(10)

(11)

l¼1 k¼1

2.3. PLF using PEM Once the points and weights corresponding to 3PEM and 5PEM are estimated, the following procedure is adopted for PLF: 1. Form the input matrices X1,X2,…,Xk as [13];

2

x1;k 6 mx1 6 Xk ¼ 4 « mx1

mx2 x2;k « mx2

… … 1 …

3 mxn mxn 7 7 « 5 xn;k

(12)

where, k ¼ 1,…,m, m ¼ 2,4 for 3PEM and 5PEM respectively. 2. For each row of Xk, a deterministic load flow is carried out to determine the corresponding values of the output variables of interest. 3. Repeat step 2 for all the rows of the matrices X1,X2,…,Xk. As a result, a total of ‘nm’ load flow computations would be carried out. 4. For each output variable of interest yi,lk set the jth moment as,

  P j   n P m j E yi;lk ¼ wl;k E yi;lk ; j ¼ 1; ::; no: of moments : l¼1 k¼1

(13) In this work the first eight moments have been used. It is to be noted that in eq. (13), yi,lk denotes the value of ith variable of interest corresponding to (lk)th load flow where l ¼ 1,…n and k ¼ 1,…m. 5. Lastly, a deterministic load flow is carried out with a vector X mean ¼ ½mx1 ; mx2 ; …; mxl ; ; ::; mxn , and let yi,m denotes the value of ith output variable of interest corresponding to this load flow and j the moment Eðyi;lk Þ is updated as:

    j  E yji;lk ¼ E yji;lk þ wm yi;m

(14) Fig. 1. Flowchart of PEM based PLF.

6. With the above calculated moments, the CDF and PDF of yi is computed using CornisheFisher series [1]. The flowchart of PEM based PLF is shown in Fig. 1 3. Wind generation models In this work, four different models of wind generators have been considered. For each of these models, the active power can be obtained from the power curve provided by the manufacturer for a given wind speed Vw. A typical turbine curve with quadratic approximation is shown in Fig. 2. The corresponding power output

expressions are given in eq. (15) [8].where, PR is the rated power of turbine, Vc,VR, and VF are cut in, rated and cut out speeds of turbine respectively.

3.1. Simple PQ model of induction generator In this model, the reactive power Q is obtained by the following expression [14].

N. Gupta / Renewable Energy 94 (2016) 96e105

Power delivered (kW)

PR

where R1 is the stator resistance, R2 is the rotor resistance, Xl1 is the stator leakage reactance, Xl2 is the rotor leakage reactance and Xm is the magnetizing reactance of induction machine. From eq. (18), the slip is calculated as,

Shedding the wind

(Rated power)

Cut in wind speed

Rated wind speed

VR

VC

Wind speed (m/s)

99

Furling or cut out wind speed

VF

ffi b±pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b2  4ac s ¼ min 2a

(19)

Knowing s, Q is computed as [9];

Xm Xl2 s2 ðXm þ Xl2 Þ þ Xl1 s2 ðXm þ Xl2 Þ2 þ R22 ðXm þ Xl1 Þ 2 Q¼ h  i2 Vt 2 R2 R1 þ s Xm  ðXm þ Xl2 ÞðXm þ Xl1 Þ

(20)

þ½R2 ðXm þ Xl1 Þ þ sR1 ðXm þ Xl2 Þ2

Fig. 2. Idealized power curve.

8 0 > > > > 2 2 > > > < PR Vw  Vc VR2  Vc2 P¼ > > > > PR > > > : 0

Q zVt2

9 > > > > > > = Vc < Vw  VR > Vw  V

> > > VR < Vw  VF > > > > ; VF < Vw

Xc  Xm X þ 2 P2 Xc Xm Vt

3.4. Semi variable speed induction generator

(15)

(16)

where, X is the sum of stator and rotor leakage reactances, Xc is the reactance of capacitor bank, Xm is magnetizing reactance, and Vt is the terminal voltage magnitude. However, it is to be noted that, the stator and rotor resistances have been neglected in eq. (16).

aR2eq þ bReq þ c ¼ 0

(21)

where. a ¼ PðR21 þ ðXl1 þ Xm Þ2 Þ  Vt2 R21 2  X2 V 2 b ¼ 2R1 PXm m t and.

3.2. Synchronous generator model In this case, either Q is specified or power factor cosf is specified [9]. If Q is specified, then Q ¼ Qsp (specified value), else Q is calculated as.

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  cos2 f Q ¼P cosf

This type of induction generator consists of pitch controlled wind turbine and wound rotor induction generator, the rotor circuit of which is connected to variable resistance whose value is varied by power electronic devices [9]. Thus, the rotor resistance is unknown and determined by the controller. To overcome this problem, the quadratic equation given in eq. (18) is recasted in terms of R2/s (denoted as Req). Hence, even when the R2 and s are unknown the quantity R2/s can be computed by solving the quadratic equation in Req. The quadratic equation for Req can be written as [9];

(17)

2 ð c ¼ PR21 ðXl2 þ Xm Þ2 þ PðXm Xm þ Xl2 ÞðXm þ Xl1 ÞÞ2  R1 ðXm þ Xl2 Þ2 Vt2

From eq. (21), Req is computed as;

ffi b±pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b2  4ac Req ¼ max 2a

(22)

Knowing Req, Q can be computed as 3.3. Pitch regulated induction generator In this case, according to the wind speed variations the pitch angle controller regulates the wind turbine blade angle [9]. In this model, for the calculated value of P, a quadratic equation involving induction generator slip (s) is solved to compute slip. The quadratic equation is given below [9];

as2 þ bs þ c ¼ 0

(18)

R2 ðX þ Xl1 Þ  ðXm þ Xl2 Þ  eq m  2 Xm  ðXm þ Xl2 ÞðXm þ Xl1 Þ 2 Q¼ h  i2 Vt 2 Req R1 þ Xm  ðXm þ Xl2 ÞðXm þ Xl1 Þ

2 þ Req ðXm þ Xl1 Þ þ R1 ðXm þ Xl2 Þ

(23)

The block diagram of wind generation system used is shown in Fig. 3, in which f(Pi,g) and f(Qi,g) are the PDFs of active and reactive powers produced and consumed, respectively, by the ith wind generation system. f(Vi,g) is the voltage PDF at ith bus.

where. 4. Inclusion of correlation in loads and wind generation a ¼ PR21 ðXl2 þ Xm Þ2 þ PðXm Xl2 þ Xl1 ðXl2 þ XmÞÞ2  Vt2 R1 ðXl2 þ Xm Þ2 2  V 2R X2 b ¼ 2PR1 R2 Xm t 2 m and c ¼ PR22 ðXl1 þ XmÞ2 þ PðR1 R2 Þ2  Vt2 R1 R22

In the operation of power systems, RES are essentially intermittent, random [1] and also spatially correlated in a significant manner [2] within a given geographical area, as they are influenced

100

N. Gupta / Renewable Energy 94 (2016) 96e105

Fig. 3. Block diagram of wind generation system.

by same physical phenomena [15]. The inclusion of correlation in wind generation and the loads for point estimate method is explained as follows. Let a vector of z correlated random variables x ¼ (x1,x2…xz)T be b x as considered with the corresponding mean vector m m b x ¼ ðm1 ; m2 …mz ÞT , where the subscript 'T' denotes the transpose of a vector. Let the correlation matrix be denoted as

2 x1 x2 x3 x4 3 1 r12 / r1z x1 x 6 r21 1 / r2z 7 7 Px ¼ 2 6 « 4 « « 1 « 5 xz rz1 rz2 / 1

(24)

where rij is the correlation between ith and jth random variable. From correlation matrix and known standard deviations the variance-covariance matrix can be obtained as [16].

2

s2x1 6r s s 21 x2 x1 Cx ¼ 6 4 « rz1 sxz sx1

r12 sx1 sx2 s2x2 « rz2 sxz sx2

/ / 1 /

3 r1z sx1 sxz r2z sx2 sxz 7 7 5 « 2 sxz

(25)

where sx1,sx2…,sxz are the standard deviations of the random variables x1,x2…,xz respectively. For including the correlation in the input variables using point estimate method, following procedure is adopted [7]:

(26)

2. Compute the matrix B as

B ¼ L1

(27)

3. Using the above matrix B, the original vector x of correlated random variables is transformed into a vector of uncorrelated variables y as

y ¼ Bx 4. Transform the central moments of x to the new space as

r¼1

lyl;4 ¼

n X

(29)

B4lr lxr;4 s4xr

r¼1

where, Blr represents the element located at lth row and rth column of the matrix B, lxr,3 and lxr,4 are the skewness and kurtosis of the rth element of the original correlated vector x respectively, sxr is the standard deviation of the rth element of the original correlated vector x, while lyl,3 and lyl,4 are the skewness and kurtosis of the lth b y is vector of element of the transformed vector y respectively, and m mean values of y. 5. In the transformed space, compute the points (yl,k) and corresponding weights (wl,k) using the point estimate method as explained in Section-II. 6. Form the (2m þ 1) and (4m þ 1) vectors (corresponding to 3PEM and 5PEM respectively) in the transformed space y, following step 1 of Section II-C, as

i h Y l ¼ my1 ; my2 ; …; myl1 ; yl;k ; mylþ1 ; …; myz

(30)

where, k ¼ 1,2 for 3PEM, k ¼ 1,…4 for 5PEM, l ¼ 1,…z 7. Transform the vector Yl into the original space by using the transformation

1. Carry out the Cholesky decomposition of Cx as

C x ¼ LLT

b y ¼ L1 m b x; m n X lyl;3 ¼ B3lr lxr;3 s3xr

(28)

X l ¼ B1 Y l :

(31)

8. For remaining nz uncorrelated variables the vector Xl is generated as given in eq. (12). 9. Follow steps 2, 3 and 4 of Section II-C for obtaining the CDFs of the output variables.

5. Results For calculating the CDFs of the reactive power absorbed by the wind generators, all the four wind generator models described in Section III have been incorporated in PLF. The PLF has been solved by using both 3PEM and 5PEM methods and the results obtained by

N. Gupta / Renewable Energy 94 (2016) 96e105

1. Sample the input random variables i.e. active and reactive loads, from the given probability distributions. 2. Obtain the output vector yi for each input vector xi by solving Newton Raphson Load Flow (NRLF). 3. Repeat steps 1 and 2 a sufficient number of times to obtain a large sample of output vector yi. 4. From the output vector yi calculate the frequency distribution of each output variable yi,k and from the frequency distribution, find the probability distribution by converting the frequency values to probability values. Here, yi,k denotes the kth component of the vector yi with k ¼ 1,2,…,NT, NT being the total number of output variables of interest. Due to the high accuracy of MCS, it is used as a benchmark for comparing the accuracy of other methods. For investigating the feasibility of the proposed methods, three different cases have been considered in this work: (a) uncorrelated loads and uncorrelated wind (b) correlated loads and uncorrelated wind (b) correlated loads and correlated wind. It is to be noted that any correlation between wind and load has been neglected in this paper. All these three scenarios have been investigated on IEEE-118 bus system. In this system it has been assumed that a total of 50 wind generators (WG) have been connected equally at bus numbers 94 to 98 (i.e. at each bus, 10 WGs have been assumed to be connected). Further, at each bus the following types of generators have been assumed: PQ type WGS at bus no. 94, synchronous generator type WGS at bus no. 95, pitch regulated WGS at bus nos. 96 and 98 and semi-variable speed type WGS at bus no. 97. It is to be noted that this distribution of different type of WGS among different buses has been assumed totally arbitrarily, any other distribution could have been assumed as well. The generator parameters and the correlation matrices between the loads and between the WGS are given in Appendix A and B respectively. Further, in this work the wind speed is assumed to have a Rayleigh distribution (a Weibull distribution with shape parameter of 2) [18]. For each of the buses from 94 to 98, the wind turbine power curve has been assumed to have the same shape as that shown in Fig. 1. The parameters of these power curves for these five buses are given in Appendix C. For modeling the uncertainties in the loads, the active power and the reactive power consumed by any load have been assumed to be normally distributed with the mean values equal to the load data given in Ref. [19] and the standard deviation equal to 10% of the corresponding mean value of the load. All the simulation studies have been carried out on an Intel Dual-core 2.67 GHz, 2 GB RAM machine, using the MATLAB environment. For implementing the proposed technique in Matlab following procedure is adopted. 1. Standard central moments of PDFs of active and reactive power loads (with known mean and standard deviation), points and weights are calculated, as explained in Section 2.1 And 2.2.

2. For wind, samples are generated using weibull distribution function and then output power is obtained by using wind turbine curve. The output power of WGS is taken as negative input power (power injection model) to the system and the weights and points are calculated. However, the value of reactive power and voltage at bus, where WGS is connected, is initially taken as zero and unity respectively. 3. Form the input matrix of the points (as given by eq. (12)). 4. For each row of input matrix, a Newton-Ralphson load flow is calculated, which is repeated for all the rows of the matrices and output variables of interest (voltage and line flows) are obtained. 5. For each iteration of NRLF, reactive power consumed by WGS is calculated and voltage at that bus (where WGS is connected) is updated. 6. For each output variable of interest, the moments are calculated and scaled by the weights. Finally, these moments are added with the moments obtained by running NRLF for mean value and scaled with 1-(sum of weights) as explained in Section 2.3.

5.1. Uncorrelated loads and uncorrelated wind generation In this case, no correlation has been assumed either among the wind generators or among the loads. The PLF has been computed using both 3PEM and 5PEM methods and the CDFs of reactive power consumption at bus no. 94 for both these methods are shown in Fig. 2. Further for the purpose of comparison, CDF has also been computed with MCS study. For performing the MCS study, deterministic load flow computations have been carried out 1,00,000 times with random variations of the wind speeds and the loads within their respective ranges. The CDF obtained with MCS studies is also shown in Fig. 4. Form Fig. 4, it is observed that the CDFs obtained by these three methods are almost overlapping each other. For qualitative assessment of the accuracies of 3PEM and 5PEM methods vise
a -vis MCS study, maximum error has been computed following the procedure described below. Initially, each of the three CDFs obtained by 3PEM, 5PEM and MCS study is discretized at points xi with corresponding value of probability (pi), where, i ¼ 1,…,N, N being the total number of points used for discretization. Let the values of (pi) obtained from the curves of MCS, 3PEM and 5PEM methods be denoted as pMCS , p3PEM , i i and p5PEM respectively. Lastly, the maximum errors in the estimated i values of probability (pi) for 3PEM and 5PEM methods are calculated as:

1 0.8

probability

these two methods have also been compared with those obtained by the Monte Carlo simulation (MCS) studies. In MCS method, sampling is carried out to obtain the state of each component in the system, which includes various system equipments such as loads, generators, transmission lines, transformers, etc. Sampling can be carried out by either using inversetransform method or acceptance rejection method [17]. Let the vector of input active and reactive power variables for ith sample be represented as xi and corresponding output voltages, angles and power flows be represented as yi. The main steps involved in MCS are:

101

MCS 5PEM 3PEM

0.6 0.4 0.2 0

5.8

6

6.2

6.4

6.6

6.8

(p.u.)

7

7.2

7.4

7.6

Fig. 4. CDF of reactive power consumption at bus no. 94 for case A.

7.8 -3 x 10

102

N. Gupta / Renewable Energy 94 (2016) 96e105

1

Em5PEM

¼ max pMCS  p5PEM ; ci ¼ 1; 2…N i i

0.8

In this study the value of N has been taken as 1000. The values of Em3PEM and Em5PEM have been found to be 1.9981  104 and 1.9784  104 respectively. Following the same procedure the maximum error for reactive power absorbtion at bus nos. 95 to 98 have also been calculated and the results are shown in Table 1. The CDFs of the bus voltage of bus no. 108 obtained by the three methods are shown in Fig. 5 for this case. Again as observed from Fig. 5, the three CDFs are almost overlapping each other. Further, following the procedure for the calculation of maximum errors described above, the values of Em3PEM and Em5PEM for the voltages of all the buses have been calculated and the results are shown in Fig. 6(a) and (b) respectively. In these two figures, it is to be noted that, these two quantities i.e. Em3PEM and Em5PEM are denoted as err_3PEM and err_5PEM respectively. From these figures it can be seen that the maximum errors are quite less (of the order of 1  104), which proves that the 3PEM and 5PEM methods are accurate enough for calculating the CDFs of the desired variables. Now, in Fig. 6(a) and (b) there are 118 values of Em3PEM and Em5PEM each. To quantify the overall maximum voltage error, the maximum of these 118 values of Em3PEM has been calculated and is 3PEM . Similarly, the maximum of 118 values of Em denoted as Verr 5PEM 5PEM . It is to be noted has also been calculated and is denoted as Verr that in the calculation of the CDFs of the bus voltages the maximum error between the MCS and 3PEM(5PEM) would be always less than 3PEM ðV 5PEM Þ. The values of these two quantities have or equal to Verr err been found to be 2.1027  104 and 2.1026  104 respectively. The CDFs of the real power flow over the line between bus no. 69 and 116 are shown in Fig. 7. Again for this case also, an excellent agreement among the CDFs can be observed. Further, the values of Em3PEM and Em5PEM for the real power flow over all the lines have been calculated and the results are shown in Fig. 8(a) and (b). From Fig. 8(a) and (b), it can be seen that the maximum errors are quite less in this case also. Now in the 118 bus system there are a total of 186 lines. Therefore, in Fig. 8(a) and (b), there are all together 186 values of Em3PEM and Em5PEM each. Following the same procedure as described earlier, the maximum of these 186 values of 3PEM ) and similarly the Em3PEM has been calculated (denoted as Perr 5PEM value of Perr has also been calculated. The values of these two quantities have been found to be 5.523  103 and 5.521  103 respectively. As these two values are quite small, the CDFs of line power flows can be estimated with sufficient accuracy through 3PEM and 5PEM methods.

Table 1 Maximum error (104) for the reactive power absorption. Bus no.

94 95 96 97 98

Case A

Case B

Case C

Em3PEM

Em5PEM

Em3PEM

Em5PEM

Em3PEM

Em5PEM

1.9981 1.9880 1.9734 1.9923 1.9832

1.9784 1.9543 1.9641 1.9872 1.9771

2.0180 2.0392 1.9985 2.1024 2.0009

1.9985 2.0010 1.9921 2.0012 1.9981

2.0318 2.1724 2.1293 2.1343 2.1002

2.0184 2.1658 2.1100 2.1298 2.0931

probability

Em3PEM

¼ max pMCS  p3PEM ; ci ¼ 1; 2…N i i

MCS 5PEM 3PEM

0.6 0.4 0.2 0

0.95

0.955

0.96

0.965

(p.u.)

0.97

0.975

0.98

0.985

Fig. 5. CDF of voltage at bus no. 108 for case A.

5.2. Correlated loads and uncorrelated wind generation In this case it has been assumed that the active and reactive power loads at bus nos. 81 to 84 are correlated with each other. The correlation among the loads has been incorporated following the procedure given in Section IV. CDFs of reactive power consumption at bus no. 95 for the three methods (3PEM, 5PEM and MCS) are shown in Fig. 9. The maximum errors for reactive power absorptions at bus nos. 94 to 98 have been calculated and are shown in Table 1. The CDFs of the bus voltage of bus no. 71 obtained by the three methods are shown in Fig. 10. It can be observed that in this case also the three CDFs are almost overlapping each other. The value of 3PEM ðV 5PEM Þ is 2.6334  104 (2.6322  104). Verr err Again the small values of these maximum errors establish the suitability of 3PEM and 5PEM methods even in the presence of correlation among the loads. The CDFs of the reactive power flow in the line between bus no. 5 and 6 are shown in Fig. 11. In this case also, an excellent agreement among the CDFs can be seen. Following the same procedure as described in the previous subsection, for each of the 186 numbers of line reactive power flows, the errors Em3PEM(Em5PEM) have been calculated. The maximum of these 186 numbers of 3PEM ðQ 5PEM Þ. The values of these Em3PEM(Em5PEM) is denoted as Qerr err two quantities have been found to be 5.7249  103 and 5.7182  103 respectively. As these values are found to be quite low, 3PEM and 5PEM methods are also suitable for obtaining the CDFs of line reactive power flows.

5.3. Correlated loads and correlated wind generation In this case, apart from the correlation among the active and reactive power loads at bus nos. 81 to 84, correlation among the WGS located at buses 94e98 has also been considered. To incorporate these correlations, the procedure described in Section IV has again been followed. The CDFs of reactive power consumption at bus no. 96 calculated using the three methods (3PEM, 5PEM and MCS) are shown in Fig. 12. The calculated maximum errors for reactive power absorptions at bus nos. 94 to 98 and are given in Table 1. The low values of these maximum errors once again establish the suitability of 3PEM and 5PEM methods even in the presence of correlation among the loads as well as among the WGS. The CDFs of the bus voltage of bus no. 114 obtained by the three methods are depicted in Fig. 13. It can be observed that in this case also the three CDFs are in excellent agreement with each other. The 3PEM ðV 5PEM Þ are 1.9769  104 (1.9732  104) for this values of Verr err

N. Gupta / Renewable Energy 94 (2016) 96e105

103

-4

err_5PEM

4

x 10

5PEM

2 0 0

20

40

60

Bus number

err_3PEM

100

80

100

(a)

-4

4

80

x 10

3PEM

2 0 0

20

40

60

Bus number (b)

Fig. 6. Calculated absolute errors for bus voltages.

0.6 0.4

0.6 0.4 0.2

0.2

1.6

1.7

1.8

(p.u)

1.9

2

0 0

2.1

10

x 10

0.005

0.01

(p.u.)

Fig. 9. CDF of reactive power consumption at bus no. 95 for case B.

Fig. 7. CDF of active power flow in line no. 69e116 for case A.

err_5PEM

0 1.5

MCS 5PEM 3PEM

0.8

-3

5PEM

5 0 0 10

err_3PEM

probability

0.8

1

MCS 5PEM 3PEM

probability

1

x 10

20

40

60

-3

80

100

120

140

160

180

100

120

140

160

180

Line number (a)

3PEM

5 0 0

20

40

60

80

Line number (b)

Fig. 8. Calculated absolute errors for line active power flows.

0.015

104

N. Gupta / Renewable Energy 94 (2016) 96e105

1

0.8

probability

0.8

probability

1

MCS 5PEM 3PEM

0.6 0.4 0.2

0.6 0.4 0.2

0 0.983

0.984

0.985

0.986

0.987

0.988

(p.u)

0.989

0.99

0.991

Fig. 10. CDF of voltage at bus no. 71 for case B.

1

0 0.945

0.95

0.955

0.96

(p.u.)

0.965

0.97

0.975

Fig. 13. CDF of voltage at bus no. 114 for case C.

MCS 5PEM 3PEM

0.8

probability

MCS 5PEM 3PEM

0.6 0.4 0.2 0 -0.4

-0.35

-0.3

-0.25

-0.2

-0.15

(p.u.)

-0.1

-0.05

0

0.05

0.1

Fig. 11. CDF of reactive power flow in line no. 5e6 for case B.

1

probability

0.8

Fig. 14. CDF of active power flow in line no. 4e6 for case C.

MCS 5PEM 3PEM

6. Conclusion In this work, 3PEM and 5PEM methods have been applied for computing the CDFs of reactive power absorbed by the wind generators in the presence of load and wind generation uncertainties. For calculating the reactive power absorption, four WGS models (two simplified and two detailed models) have been used in this paper. Based on the simulation studies carried out in this work, following conclusions can be drawn:

0.6 0.4 0.2 0 3

4

5

6

7

(p.u.)

8

9

10

11

12

x 10

-3

Fig. 12. CDF of reactive power consumption at bus no. 96 for case C.

case. The CDFs of the active power flow in the line between bus no. 3PEM and P 5PEM are 4 and 6 are shown in Fig. 14. The values of Perr err 8.6243  103 and 8.6210  103 respectively. In all the above three cases, it has been found that the 3PEM and 5PEM methods are capable of computing the CDFs of the any desired quantity quite accurately. A comparison of the computa
-vis the MCS tional times required by these two method vis-a studies is given in Table 2 for all the above three cases. It can be seen from this table that it takes approximately 6 h to complete the MCS study with 1,00,000 deterministic load flow solutions, while the time taken by 3PEM and 5PEM methods is even less than 3 min. However, the time taken by 5PEM method is marginally more as compared to that taken by 3PEM method.

1. Both 3PEM and 5PEM methods are capable of computing the CDFs of any desired quantity (bus voltage, line real and reactive power flow and reactive power absorbed by WGS) accurately in the presence of load and wind generation uncertainties with and without any correlation among them. 2. The computation burden required by these two methods is significantly less as compared to that required by the MCS study.

Table 2 Time required for simulation studies (in seconds). Scenario/scheme

3PEM

Uncorrelated loads and uncorrelated wind generation Correlated loads and uncorrelated wind generation Correlated loads and correlated wind generation

120.327 147.179 21597

5PEM

MCS

125.174 153.542 23466 137.573 159.184 24752

N. Gupta / Renewable Energy 94 (2016) 96e105

References

Appendix A. induction generator circuit parameters

Type of WGS

Pitch regulated fixed speed [9]

Semi-variable speed [20]

Rating (MW) Rated (kV) R1 Xl1 R2 Xl2 Xm Xc cosf ¼ 0.8

0.5 0.69 0.005986 0.08212 0.01690 0.107225 2.5561 2.5561

1 0.69 0.005671 0.15250 0.00462 0.096618 2.8985 2.8985

Appendix B. Correlation matrix for load and wind generators Correlation matrix between the loads at bus no. 81 to 84.

2

P load

1:0000 6 0:4529 6 ¼4 0:9058 0:0234

0:4529 1:0000 0:4978 0:1036

0:9058 0:4978 1:0000 0:0068

3 0:0234 0:1036 7 7 0:0068 5 1:0000

Correlation matrix between the wind generators at bus no. 94 to 98.

2

P wind

1:0000 6 0:4529 6 ¼6 6 0:9058 4 0:0234 0:4664

0:4529 1:0000 0:4978 0:1036 0:6099

0:9058 0:4978 1:0000 0:0068 0:4561

0:0234 0:1036 0:0068 1:0000 0:0368

105

3 0:4664 0:6099 7 7 0:4561 7 7 0:0368 5 1:0000

Appendix C. Power curve parameters

Bus no.

Vcutin (m/s)

Vrated (m/s)

Vcutout (m/s)

Prated (MW)

94 95 96 97 98

4.5 4.5 4.2 4.1 4.3

15 16 14 15 16

50 55 57 45 47

0.50 0.55 0.60 0.40 0.45

[1] J. Usaola, Probabilistic load flow in systems with wind generation, IET Gener. Transm. Distrib. 3 (12) (2009) 1031e1041. [2] J. Usaola, Probabilistic load flow with correlated wind power injections, Electr. Power Syst. Res. 80 (5) (2010) 528e536. [3] J. Usaola, Probabilistic load flow with wind production uncertainty using cumulants and cornish- fisher expansion, Int. J. Electr. Power Energy Syst. 31 (9) (2009) 474e481. [4] Y. Yuan, J. Zhou, P. Ju, J. Feuchtwang, Probabilistic load flow computation of a power system containing wind farms using the method of combined cumulants and gram-charlier expansion, IET Renew. Power Gener. 5 (6) (2011) 448e454. [5] C. Chen, W. Wu, B. Zhang, H. Sun, Correlated probabilistic load flow using a point estimate method with nataf transformation, Int. J. Electr. Power Energy Syst. 65 (2015) 325e333. [6] X. Ai, J. Wen, T. Wu, W.-J. Lee, A discrete point estimate method for probabilistic load flow based on the measured data of wind power, in: IEEE Industry Applications Society Annual Meeting, 2012, pp. 1e7. [7] J. Morales, L. Baringo, A. Conejo, R. Minguez, Probabilistic power flow with correlated wind sources, IET Generation, Transm. Distrib. 4 (5) (2010) 641e651. [8] D. Villanueva, J. Pazos, A. Feijoo, Probabilistic load flow including wind power generation, IEEE Trans. Power Syst. 26 (3) (2011) 1659e1667. [9] K. Divya, P.N. Rao, Models for wind turbine generating systems and their application in load flow studies, Electr. Power Syst. Res. 76 (9e10) (2006) 844e856. [10] J. Morales, J. Perez-Ruiz, Point estimate schemes to solve the probabilistic power flow, IEEE Trans. Power Syst. 22 (4) (2007) 1594e1601. [11] C.-L. Su, Probabilistic load-flow computation using point estimate method, IEEE Trans. Power Syst. 20 (4) (2005) 1843e1851. [12] A.C. Miller, T.R. Rice, Discrete approximations of probability distributions, Manag. Sci. 29 (3) (1983) 352e362. [13] N. Gupta, V. Pant, B. Das, Probabilistic load flow incorporating generator reactive power limit violations with spline based reconstruction method, Electr. Power Syst. Res. 106 (2014) 203e213. [14] A. Feijoo, J. Cidras, Modeling of wind farms in the load flow analysis, IEEE Trans. Power Syst. 15 (1) (2000) 110e115. [15] J. Morales, R. Minguez, A. Conejo, A methodology to generate statistically dependent wind speed scenarios, Appl. Energy 87 (3) (2010) 843e855. [16] Z. Qin, W. Li, X. Xiong, Generation system reliability evaluation incorporating correlations of wind speeds with different distributions, IEEE Trans. Power Syst. 28 (1) (2013) 551e558. [17] W.L. Martinez, A.R. Martinez, Computational statistics handbook with matlab, Chapman Hall/CRC (2002), 3rd edition (December 22, 2015). [18] J. Carta, P. Ramirez, S. Velazquez, A review of wind speed probability distributions used in wind energy analysis: Case studies in the canary islands,, Renew. Sustain. Energy Rev. 13 (5) (2009) 933e955. [19] Power systems test case archive (online), available: http://www.ee. washington.edu/research/pstca. [20] Z. Lubosny, Wind Turbine Operation in Electric Power Systems Advanced Modeling, Springer Verlag, 2003.