1 Introduction

Load flow calculation is one of the most important tools to analyze static characteristics of electric power systems [1]. Traditional deterministic load flow method calculates system states and load flows in terms of specific values of bus injections for a determinate network configuration. As a result, uncertainties in power systems such as generator and load fluctuations, component faults [2] cannot be considered. With increasing concerns on emission reduction and sustainable utilization of renewable energy, distributed generation including wind turbine (WT) and photovoltaic panel developed very fast in recent years all over the world, challenging the conventional power system operation and planning [3]. Different from traditional generations, outputs of wind power generation (WPG) are usually not controllable due to the wind stochasticity and intermittency in natural environment. Therefore, the integration of renewable generation increases the difficulty for grid operators to estimate system reserve, dispatch and manage system load flow [4].

Probabilistic load flow (PLF) method was first proposed in 1970s to take into account the facts of unscheduled outages, load forecast inaccuracies and measurement errors [5]. PLF can yield the probabilistic density functions (PDFs) and cumulative distribution functions (CDFs) of bus voltages and branch flows. PLF results are able to provide grid operators and planners with an explicit perspective on present or future system conditions, helping them to reduce arbitrariness on decision making. The most used technique in PLF calculation is Monte Carlo simulation (MCS) [6, 7]. It actually involves repeated load flow processes with various values of input variables yielded by their PDFs. Based on the MCS, a Latin supercube sampling method have been applied to improve the efficiency of MCS [8]. However, MCS is still remarkably time consuming because a huge amount of simulations are required for an accurate enough outcome. Therefore, MCS is generally not adequate in practical applications especially for complex and bulk power system. However, due to its simple processes and high accuracy, it is usually used as a baseline for methods comparison.

Different from MCS, the PDF convolution technique is an analytical PLF calculation method based on the assumption that all input variables are independent with each other. By applying DC load flow model [5] or linearized AC load flow model [9], bus voltages and branch flows can be represented as a linear combination of input variables, which makes convolution procedure feasible. However, this technique is still inefficient for real-life power systems due to its inherent complexity [10]. Fast Fourier transform based convolution technique [11] was then proposed to reduce the calculation amount, but it is still unable to effectively handle the problem.

The combined cumulant (semi-variant) algorithm and Gram–Charlier (GC) expansion theory can be used to convert the complicated convolution calculation into a simple arithmetic process according to the superior properties of cumulants [12]. This method greatly enhances the calculation speed and is able to accurately approximate the PDF or CDF of output variables [13,14,15,16,17]. Based on the GC theory, several improved approaches were employed to address PLF analysis, such as the stochastic collocation interpolation [18], and polynomial chaos expansion [19, 20]. However, the GC expansion based PLF (GC-PLF) calculation method may mistakenly obtain negative values in tail regions of PDF, which limits system planners to determine line security and maximum power rating [21, 22]. Reference [21] adopted an improved C-type GC expansion to avoid gaining negative values results. Reference [22] introduced the principle of maximum entropy (ME) from information theory to reconstruct PDF from an entirely different perspective according to the moments of output variables, guaranteeing positive PDF results. Work in [22] showed the superiority of ME to GC in both mean results and operational limit predictions. However, the algorithm did not consider the integration of distributed generation.

It should be noticed that the stochastic output of WPG is not a normal distribution and has a complex relationship with wind speed which usually follows Weibull distribution. Therefore, WPG is completely different from traditional generation units, and may cause severe fluctuation and skew on PDF curves of bus voltages and branch flows. So far, most researches obtained the PDF of WPG based on an assumption that the relationship between wind speed and WPG is linear [14,15,16]. This can introduce unacceptable errors for the cumulant calculation of WPG. In addition, GC PDF reconstruction method is applied directly to obtain the PDF of bus voltages and branch flows, but the results are inaccurate in some cases. Lastly, dependence between active and reactive power injections of WPG caused by constant power factor control strategy is still not well considered, which also introduces extra errors to PDF calculations.

In this paper, a ME based PLF (ME-PLF) calculation algorithm for power systems integrated with WPG is proposed to improve the PDF accuracy of bus voltages and branch flows. The proposed approach can also deal with linear dependence between active and reactive power injections of WPG by the modified cumulant calculation framework.

The rest of the paper is organized as follows. In Sect. 2, the ME reconstruction method along with probability theories are reviewed. In Sect. 3, the mathematical models of the basic loads and WPG, the modified cumulant calculation based on linearized load flow formulation and the procedure of ME-PLF calculation algorithm are detailed. Case studies are presented in Sect. 4. Finally, Sect. 5 gives the conclusions.

2 ME reconstruction method

In this section, theoretical background including moment and cumulant is first introduced. Then, the principle of ME and the model of ME reconstruction method are presented to find the most unbiased probability distribution.

2.1 Moment and cumulant

Defining p(x) as the PDF of a continuous random variable x. For a positive integer n, the function xn is integrable with respect to p(x) over (−∞, +∞), the integral shown in (1) is called the nth raw moment of p(x) [12]. If the p(x) of the random variable x is unknown, the raw moment αn can be gained from its historical data by statistic way, as shown in (2).

$$\alpha_{n} = \int_{ - \infty }^{ + \infty } {x^{n} p(x){\text{d}}x}$$
(1)
$$\alpha_{n} = E\left( {x^{n} } \right)$$
(2)

where E denotes the calculation of expectation.

Defining characteristic function ϕ(t) as the Fourier transform of p(x):

$$\phi (t) = \int_{ - \infty }^{ + \infty } {\text{e}^{{\text{j}}tx}p(x){\text{d}}x}$$
(3)

where \(\text{j} = \sqrt { - 1}\).

If the Nth raw moment of p(x) exists, the characteristic function ϕ(t) can be developed in MacLaurin’s series for small value of t, as follows:

$$\phi (t) = 1 + \sum\limits_{n = 1}^{N} {\frac{{\alpha_{n} }}{n!}} (\text{j}t)^{n} + O(t^{N} )$$
(4)

The nth order cumulants γn are then defined by:

$$\ln \phi (t) = \sum\limits_{n = 1}^{N} {\frac{{\gamma_{n} }}{n!}} (\text{j}t)^{n} + O(t^{N} )$$
(5)

The cumulants and raw moments of a random variable can be converted to each other [23] by:

$$\left\{ \begin{array}{ll} \gamma_{n} = \alpha_{n} & n = 1 \hfill \\ \gamma_{n + 1} = \alpha_{n + 1} - \sum\limits_{k = 1}^{n} {{\text{C}}_{n}^{k} \alpha_{k} \gamma_{n - k + 1} } & n \ge 2 \hfill \\ \end{array} \right.$$
(6)

If a random variable Y has a linear relationship with other m independent random variables, e.g., Y = a0 + a1X1 + ··· + amXm, then the cumulants of Y can be obtained by:

$$\gamma_{Y,n} = \left\{ \begin{aligned} &a_{0} + a_{1} \gamma_{{X_{1} ,n}} + \cdots + a_{m} \gamma_{{X_{m} ,n}} & n = 1 \hfill \\ &a_{1}^{n} \gamma_{{X_{1} ,n}} + a_{2}^{n} \gamma_{{X_{2} ,n}} \cdots + a_{m}^{n} \gamma_{{X_{m} ,n}} & n \ge 2 \hfill \\ \end{aligned} \right.$$
(7)

where γY,n denotes the nth order cumulant of Y; \({{\gamma}_{X_{1}, n}}\), \({{\gamma}_{X_{2}, n}}\), …, \({{\gamma}_{X_{m}, n}}\) denote the nth order cumulants of X1, X2, …, Xm, respectively; a n m denotes am to the nth power. This calculation property makes cumulants much simpler than raw moments in theoretical deduction.

2.2 ME reconstruction method

The principle of ME [24] states that subject to the given information, the probability distribution which best represents the current state of knowledge is the one with maximum Shannon entropy. It conveys an idea that the ME estimate is maximally noncommittal with respect to unknown information, any other estimations will lead to biased results.

In detail, ME is a nonlinear probability density function reconstruction method that maximizes Shannon entropy H of the probability distribution [25] shown in (8) while at same time satisfying the constraints shown in (9).

$$\hbox{max} \, H = - \int {p(x)\ln p(x){\text{d}}x}$$
(8)
$$\int {p(x)} \varphi_{n} (x){\text{d}}x = \alpha_{n} \quad n = 0, \, 1, \, \ldots , \, N$$
(9)

where α0 = 1; α1, α2, …, αN are known raw moments of p(x); φ0(x) = 1; φn(x) = xn, n = 1, 2, …, N.

The entropy distributions will be of the following form through Lagrange multiplier method:

$$p(x) = \text{exp}\left( { - \sum\limits_{n = 0}^{N} {\lambda_{n} \varphi_{n} (x)} } \right)$$
(10)

where \(\lambda_{n}\) is Lagrange multiplier, which can be solved by (11), and is obtained by substituting (10) into (9).

$$\int {\varphi_{n} (x)\exp \left( { - \sum\limits_{n = 0}^{N} {\lambda_{n} \varphi_{n} (x)} } \right){\text{d}}x = \alpha_{n} } \quad n = 0,1, \ldots ,N$$
(11)

In the numerical calculation for solving (11), Newton iteration method is used [26], and the upper and lower limits of integral are described by interval (μ−σ, μ+σ), where μ and σ are the mean and standard deviation of random variable x; is an arbitrary real number which determines the length of integral interval. The values of μ and σ can be obtained in cumulant calculation process, which will be discussed in the next section. In this paper, the value of is selected with 6 considering the efficiency and accuracy of PLF calculation.

3 ME-PLF calculation algorithm

In this section, the probabilistic model of basic loads and WPG injections are first introduced. Then, the linearized load flow formulation and modified cumulant calculation process considering linear dependence between active and reactive power injections of WPG are elaborated. Lastly, the calculation procedure of ME-PLF calculation algorithm for power system integrated with WPG is presented.

3.1 Probabilistic model of power injections

In this paper, two significant uncertainties in power systems are considered, i.e. load fluctuations and WPG injections.

For active load PD and reactive load QD, normal distribution is commonly used to describe their probabilistic characteristics [5]. The PDFs of PD and QD can be formulated as:

$$\left\{ \begin{aligned} f(P_{D} ) = \frac{1}{{\sqrt {2\uppi} \sigma_{{P_{D} }} }}\exp \left[ { - \frac{{\left( {P_{D} - \mu_{{P_{D} }} } \right)^{2} }}{{2\sigma_{{P_{D} }}^{2} }}} \right] \hfill \\ f(Q_{D} ) = \frac{1}{{\sqrt {2\uppi} \sigma_{{Q_{D} }} }}\exp \left[ { - \frac{{\left( {Q_{D} - \mu_{{Q_{D} }} } \right)^{2} }}{{2\sigma_{{Q_{D} }}^{2} }}} \right] \hfill \\ \end{aligned} \right.$$
(12)

where \({\mu}_{P_{D}}\) and \( \sigma_{{P_{D} }}^{2} \) are the expectation and variance of PD; \({\mu}_{Q_{D}}\) and \( \sigma_{{Q_{D}}}^{2} \) are the expectation and variance of QD. So, their raw moments can be immediately obtained by (1).

For WPG, uncertainties mainly result from the intermittent characteristic of wind speed, which can be described by Weibull distribution [27]. In this paper, the PDF for wind of speed is assumed to be a Weibull distribution [28] as follows:

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

where \(\upsilon\) is the wind speed; k is the shape parameter; c is a scale parameter which indicates the average wind speed. Generally, k and c are set to be 2.0 and 8.5 [21, 29], or other values [14, 16]. These will be further compared in case studies.

The WPG can be calculated from the wind speed by the wind power curve [29], which is formulated as follows:

$$P_{w} (\upsilon ) = \left\{ \begin{array}{ll} 0 &\quad \upsilon < \upsilon_{ci} {\text{ or }}\upsilon > \upsilon_{co} \, \hfill \\ \frac{{P_{R} }}{{\upsilon_{R}^{3} - \upsilon_{ci}^{3} }}(\upsilon^{3} - \upsilon_{ci}^{3} ) &\quad \upsilon_{ci} \le \upsilon < \upsilon_{R} \hfill \\ P_{R} &\quad \upsilon_{R} \le \upsilon \le \upsilon_{co} \hfill \\ \end{array} \right.$$
(14)

where Pw is the active power output of WT; PR is the rated output power; \(\upsilon_{ci}\) is the cut-in wind speed; \(\upsilon_{R}\) is the rated wind speed; \(\upsilon_{co}\) is the cut-out wind speed. The parameters can be selected in terms of [29].

In this paper, the active power generation moments of WT can be calculated from the historical data record generated by (13) and (14). It is simpler than the mathematical approach of gaining WPG distribution and computing its moments.

Asynchronous generator model is common applied in most WTs. WT usually absorbs a large amount of reactive power to provide the excitation current for establishing the magnetic field while generating active power. So WT can be simplified as PQ bus and have a constant power factor through the automatic capacitor switching [30]. The reactive power of WT can be formulated as:

$$Q_{w} = P_{w} \tan (\delta )$$
(15)

where tan(δ) is the proportional coefficient between active and reactive power. Obviously, there is a linear dependence relationship between the active and reactive power injections of WPG.

3.2 Modified cumulant calculation

The original AC load flow can be formulated as:

$$\left\{ \begin{aligned} &{\varvec{W}} = \varvec{F}({\varvec{V}}) \hfill \\ &{\varvec{Z}} = \varvec{G}({\varvec{V}}) \hfill \\ \end{aligned} \right.$$
(16)

where W is the bus power injection vector; V is the bus voltage vector; Z is the branch power flow vector; F and G denote the power injection function and line flow function, respectively.

Expand (16) at the operating point in terms of Taylor’s series and neglect the high order items, as follows:

$$\left\{ \begin{aligned} &\Delta {\varvec{W}} = {\varvec{J}}_{0} \Delta {\varvec{V}} \hfill \\ &\Delta {\varvec{Z}} = {\varvec{G}}_{0} \Delta {\varvec{V}} \hfill \\ \end{aligned} \right.$$
(17)

where J0 = ∂F(V0)/∂V; G0 = ∂G(V0)/∂V.

Equation (17) can also be expressed as:

$$\left\{ \begin{aligned}& \Delta {\varvec{X}} = {\varvec{J}}_{0}^{ - 1} \Delta {\varvec{V}} = {\varvec{S}}_{0} \Delta {\varvec{V}} \hfill \\ &\Delta {\varvec{Z}} = {\varvec{G}}_{0} {\varvec{S}}_{0} \Delta {\varvec{V}} = {\varvec{T}}_{0} \Delta {\varvec{V}} \hfill \\ \end{aligned} \right.$$
(18)

where S0 is the inverse matrix of J0, which is called the sensitivity matrix; T0 = G0S0 is another sensitivity matrix [11].

So the calculation of (V, Z) can be solved in three steps, first is to calculate the operating point (V0, Z0) through deterministic load flow calculation, then is to calculate the indeterminate item (∆V, ∆Z) by PLF calculation method, and last is to add the two parts together, as follows:

$$\left\{ \begin{aligned} &{\varvec{V}} = {\varvec{V}}_{0} + \Delta {\varvec{V}} = {\varvec{V}}_{0} + {\varvec{S}}_{0} \Delta {\varvec{W}} \hfill \\ &{\varvec{Z}} = {\varvec{Z}}_{0} + \Delta {\varvec{Z}} = {\varvec{Z}}_{0} + {\varvec{T}}_{0} \Delta {\varvec{W}} \hfill \\ \end{aligned} \right.$$
(19)

From (19), if the power system bus number is R, the state variables of bus voltage and branch flow (i.e. vi and zi) can be formulated as:

$$\left\{ \begin{aligned} v_{i} = v_{i0} + \sum\limits_{r = 1}^{R} {s_{ir0} \Delta w_{r} } \hfill \\ z_{i} = z_{i0} + \sum\limits_{r = 1}^{R} {t_{ir0} \Delta w_{r} } \hfill \\ \end{aligned} \right.$$
(20)

According to (20), the output variables are actual the linear sums of the power injections of which the weight coefficients are corresponded to the sensitivity elements assuming all the injections are independent from each other. Therefore the cumulant method can be adopted to obtain the numerical characteristics of the output variables.

For bus i integrated with WT, some power injections become dependent according to (15), so that (20) should be modified to consider the linear dependence between WPG’s active power injection Δpwi and reactive power injection Δqwi.

The power injection of bus i is corresponded to two items ∆wk and ∆wm in (20), as follows:

$$\left\{ \begin{aligned} \Delta w_{k} = \Delta p_{Di} * \Delta p_{wi} \hfill \\ \Delta w_{m} = \Delta q_{Di} * \Delta q_{wi} \hfill \\ \end{aligned} \right.$$
(21)

where * represents convolution calculation; ∆pDi and ∆qDi are the probabilistic part of basic active and reactive loads.

According to (7), (21) can be converted to a much simpler calculation formulation for cumulants, as follows:

$$\left\{ \begin{aligned} \Delta w_{k}^{(n)} = \Delta p_{Di}^{(n)} + \Delta p_{wi}^{(n)} \hfill \\ \Delta w_{m}^{(n)} = \Delta q_{Di}^{(n)} + \Delta q_{wi}^{(n)} \hfill \\ \end{aligned} \right.$$
(22)

where \(\Delta p_{Di}^{(n)}\) and \(\Delta q_{Di}^{(n)}\) are the nth cumulant of basic active and reactive loads, respectively; \(\Delta p_{wi}^{(n)}\) and \(\Delta q_{wi}^{(n)}\) are the nth cumulant of WPG’s active and reactive power injections, respectively; \(\Delta w_{k}^{(n)}\) and \(\Delta w_{m}^{(n)}\) are the nth cumulant of ∆wk and ∆wm.

Considering the linear dependence relationship as shown in (23), Δpwi and Δqwi should be put together into a power injection group (1 + tan(δi))Δpwi, which is considered to be independent of all other independent power injections [31]. So (22) can be modified into two forms considering the sensitivity coefficients, as shown in (24) and (25).

$$\Delta q_{wi} = \tan (\delta_{i} )\Delta p_{wi}$$
(23)
$$\left\{ \begin{array}{ll} \Delta w_{vk}^{(n)} = \Delta p_{Di}^{(n)} + \left( {1 + \frac{{s_{im0} \tan (\delta_{i} )}}{{s_{ik0} }}} \right)\Delta p_{wi}^{(n)} \hfill \\ \Delta w_{vm}^{(n)} = \Delta q_{Di}^{(n)} \hfill \\ \end{array} \right.$$
(24)
$$\left\{ \begin{array}{ll} \Delta w_{zk}^{(n)} = \Delta p_{Di}^{(n)} + \left( {1 + \frac{{t_{im0} \tan (\delta_{i} )}}{{t_{ik0} }}} \right)\Delta p_{wi}^{(n)} \hfill \\ \Delta w_{zm}^{(n)} = \Delta w_{vm}^{(n)} = \Delta q_{Di}^{(n)} \hfill \\ \end{array} \right.$$
(25)

Equations (24) and (25) represent the modification for the power injections in case of computing the bus voltage vi and branch flow zi, respectively. Thus the reactive power injection of WPG is eliminated and the power injection of WPG can be only corresponded to active power injection. So the probabilistic power injection matrix ΔW can be modified into two forms in terms of (24) and (25), as follows:

$$\Delta {\varvec{W}}_{v} = \left[ {\Delta w_{1} \, \cdots \, \Delta w_{vk} \, \cdots \, \Delta w_{vm} \, \cdots \, \Delta w_{R} } \right]^{\text{T}}$$
(26)
$$\Delta {\varvec{W}}_{z} = \left[ {\Delta w_{1} \, \cdots \, \Delta w_{zk} \, \cdots \, \Delta w_{zm} \, \cdots \, \Delta w_{R} } \right]^{\text{T}}$$
(27)

For other buses integrated with WTs, the relevant elements of ΔW can be modified according to (21)-(27) similarly. The cumulants of output variables can be gained by (7), as follows:

$$\left\{ \begin{aligned} &\Delta {\varvec{V}}^{(n)} = {\varvec{S}}_{0}^{n} \Delta {\varvec{W}}_{v}^{(n)} \hfill \\ &\Delta {\varvec{Z}}^{(n)} = {\varvec{T}}_{0}^{n} \Delta {\varvec{W}}_{z}^{(n)} \hfill \\ \end{aligned} \right.$$
(28)

where ∆V(n) and ∆Z(n) are the nth order cumulants of ∆V and ∆Z, respectively; \(\Delta{\varvec{W}}_{v}^{(n)}\) and \(\Delta {\varvec{W}}_{z}^{(n)}\) are the nth order cumulants of ∆Wv and ΔWz, respectively; S n0 and T n0 denote elements of S0 and T0 raised to the power n.

3.3 Procedure of ME-PLF calculation algorithm

According to (28), the cumulants of ∆V and ∆Z from the 1st to Nth order can be calculated, these cumulants can then be converted to raw moments in terms of (6). The raw moments are regarded as the constraints for ME reconstruction model (8) and (9), so the PDFs of ∆V and ∆Z can be obtained by solving (11). So far, we get the PDFs of the probabilistic part of V and Z. Finally, add them to V0 and Z0 (i.e. the expectations of V and Z) so as to get the PDFs of V and Z.

The flow chart of ME-PLF calculation algorithm is shown in Fig. 1.

Fig. 1
figure 1

Flow chart of proposed ME-PLF calculation algorithm

Full size image

4 Case studies

In this section, IEEE 118-bus test system integrated with WPG is first studied to demonstrate the accuracy and efficiency of the proposed method. Essential differences between ME and GC reconstructions are analyzed in detail. Properties of the accuracy of ME-PLF calculation algorithm with different orders are also studied with different WT locations. In addition, the impact of power factors on PLF results is elaborated. The 1999 to 2000 winter peak Polish 2383-bus system [32] is then introduced to confirm the effectiveness of the proposed method in a large scale power system. All case studies are performed with MATLAB platform. The programs for solving ME model are produced referring to [26], and power flow simulation are completed by MATPOWER [33] (Intel i5-4460 3.20 GHz, 8 GB, Windows 7).

Assume that the variances of all bus active and reactive basic loads are 10% of their expectations. The capacity of each WT is set to be 1.5 MW, with the cut-in wind speed \(\upsilon\)ci = 5 m/s, the rated wind speed \(\upsilon\)R = 15 m/s and the cut-out wind speed \(\upsilon\)co = 25 m/s. The proportional coefficient tan(δ), shape parameter k and scale parameter c of Weibull distribution for wind speed are set in three scenarios, as shown in Table 1. As the comparative baseline, maximum sampling number of the MCS is set to be 50000.

4.1 Test case of IEEE 118-bus system

4.1.1 Accuracy and efficiency of ME-PLF calculation algorithm

In this case, bus 101 is integrated with 20 × 1.5 MW WTs based on IEEE 118-bus system. Wind speed parameters k and c in Scenarios 1 and 2 are set to be different to study the effects of wind power characteristics on accuracies of ME-PLF calculation algorithm and GC-PLF calculation method. Previous works on PLF calculation integrated with WTs presented GC effectiveness in terms of Scenario 1 [14, 16]. However, [27] and [28] note that the shape parameter k should be a value between 1.8 ~ 2.3. Generally, it is set to be 2.0 [21, 29]. So wind speed parameters in Scenario 2 is more practical than those in Scenario 1. The discrete PDFs for WPG of a WT in Scenarios 1 and 2 are shown in Fig. 2.

Fig. 2
figure 2

Discrete PDF for WPG of a WT

Full size image

The voltage magnitudes V101 and V102, the branch active load flows P100-101 and P101-102 are calculated by GC-PLF calculation method with the constraint order of 8th (GC-8th) and ME-PLF calculation algorithm with the constraint orders from 4th to 8th (ME-4th to ME-8th). In order to compare the results, the average root mean square (ARMS) [12] is calculated using MCS as a reference. ARMS is used for comparing the accuracies of different calculation methods. The value of ARMS is defined as:

$$ARMS = \frac{{\sqrt {\sum\limits_{i = 1}^{M} {\left( {PLF_{i} - MCS_{i} } \right)^{2} } } }}{M}$$
(29)

where PLFi is the ith point’s value on the cumulative distribution curve calculated using the ME-PLF calculation algorithm or GC-PLF calculation method; MCSi is the ith point’s value on the CDF calculated using the MCS; M denotes the number of points. M = 1000 is chosen to ensure legible plots and with the same M used for all comparative calculations.

Table 2 and Table 3 show the ARMS results of V101, V102, P100-101 and P101-102 in Scenario 1 and Scenario 2, respectively. Figure 3 compares ME-PLF calculation algorithm and GC-PLF calculation method for V101 and P101-102 particularly. From Fig. 3, Tables 2 and 3, it is found that GC-PLF calculation method and ME-PLF calculation algorithm get approximate ARMS for voltage magnitudes. GC-8th has a larger ARMS than ME-4th for branch flows calculation in Scenario 2, while in Scenario 1 the differences are slight. ME-4th can give satisfactory results, and an apparent decline in ARMS emerges with ME-6th. It is worth noting that the decline of ARMS is non-monotonic with the increase of order of ME-PLF calculation algorithm, because the cumulant calculation errors exist resulted from the linearized AC load flow model.

Fig. 3
figure 3

Comparison of ARMS for V101 and P101–102

Full size image

Figures 4 and 5 show the PDF curves of V101 and P101–102 in Scenario 1 and Scenario 2, respectively. GC-8th gains negative PDF with the tail regions in branch flow curves, and has a considerable distortion with the MCS baseline in Scenario 2. ME-6th gives a reliable result and ME-8th obtains a more accurate PDF curves refer to MCS. Despite of the change of scenarios, ME-PLF calculation algorithm always gets accurate PDF curves.

Fig. 4
figure 4

Probability density curves comparison in Scenario 1

Full size image
Fig. 5
figure 5

Probability density curves comparison in Scenario 2

Full size image

Obviously, the accuracy of GC-PLF calculation method can be influenced by the change of scenarios, because the two scenarios adopt different sets of wind speed parameters, which makes a direct impact on wind power injection characteristics shown in Fig. 2. GC-PLF calculation method can give reliable calculation results in Scenario 1, but in most cases, the wind speed is with the parameters of Scenario 2, which will lead to significant errors by GC-PLF calculation method. ME-PLF calculation algorithm gives accurate PDF curves in both scenarios.

Table 4 shows the 9th to 12th cumulants relative errors of P100-101 and P101-102 PDF curves obtained by ME-PLF calculation algorithm and GC-PLF calculation method with MCS baseline based on the same 1st to 8th order cumulants constraints. The result illustrates that the PDF results gained by ME-PLF calculation algorithm can give a better estimation of the unknown higher cumulants information than those of GC-PLF calculation method, indicating that ME-PLF calculation algorithm is able to get more accurate PDF results from another point of view.

From the study, accuracy of GC-PLF calculation method is much more susceptible than ME-PLF calculation algorithm with different wind speed parameters, because the error of GC reconstruction mainly caused by truncation, which is substantially different from ME in dealing with the unknown information. The truncation with GC expansion means a biased estimation for the true distribution by partly discarding the uncertainties of unknown information. Whereas ME aims at finding an estimation that maximizes PDF entropy which is maximally noncommittal with regard to unknown information so as to avoid an arbitrary reconstruction.

For time consumptions, 50000 times MCS is 242.6 s. Both GC-8th and ME-8th methods are more efficient, with 0.96 s and 1.04 s, respectively.

4.1.2 Impacts of WT location

So far, ME-4th can obtain preferable PDF curves than GC-8th, and ME-6th is adequate to achieve accurate PDF curves in most situations. In the following study, buses 101, 102, and 106 are integrated with 20 × 1.5 WM, 20 × 1.5 WM and 40 × 1.5 WM WTs, respectively. The WT parameters are formed with Scenario 2.

Table 5 shows the ARMS of GC-PLF calculation method and ME-PLF calculation algorithm, Fig. 6 shows the P101-102 and P100-106 PDF curves of the two methods. It is noted that ME-4th gets a worse curve than GC-8th for P101-102, and ME-6th obtains an approximate PDF curve with GC-8th, both of which are inadequate to reach accurate enough results comparing with MCS.

Table 1 Wind parameter scenarios
Full size table
Table 2 ARMS results comparison in Scenario 1
Full size table
Table 3 ARMS results comparison in Scenario 2
Full size table
Table 4 Comparison of cumulant estimation error results
Full size table
Table 5 Maximum ARMS results comparison of multi-bus WTs case
Full size table
Fig. 6
figure 6

Probability density curves comparison of branch flows

Full size image

The reason for this phenomenon is that both sides of the branch from bus 101 to bus 102 are integrated with wind power injections, which can form a symmetry effect on the branch flow distribution. This feature makes the skew and higher odd orders cumulants of the branch flow distribution sharply decreased, meanwhile rises the higher even orders cumulants. Therefore, the PDF of branch flow is tend to be an approximate normal distribution, which makes the GC-PLF calculation method effectual in this case. From Fig. 7, ME-8th obtains the most accurate PDF for P101-102 comparing with others.

Fig. 7
figure 7

Probability density curves comparison of P101–102

Full size image

4.1.3 Dependence between wind active and reactive power associated with wind power factor

The constant wind power factor control leads to a linear dependence relationship between active and reactive power injections of WPG, and the factor is usually negative. Some studies consider the wind reactive power injection absorbed by compensators so that wind power integration buses only correspond with the active power injection, which can avoid discussing the dependence issue. This treatment is not practical, because it cannot consider wind power factor fluctuation and introduces extra errors to PDF calculation results.

This case will discuss the effects of wind power factor change on PDFs of voltages and branch flows, so it is necessary to consider the reactive power injection. If we regard the active and reactive injections as dependent random variables, the PDF curves should be different from the true ones.

Figure 8 shows the PDF curves of V101 and P101-102 in case of Scenario 2 (20 × 1.5 WM WTs are integrated into bus 101) considering two different assumptions to the correlation of wind active and reactive power injections, i.e., independent and linear dependent. The results show that the correlation of power injections mainly influences the shapes of PDF curves of voltage magnitude. Figure 8 also demonstrates that ME-8th can consider the linear dependence between active and reactive power injections of WPG comparing with MCS baseline. So when the wind power factor alters, the PDFs change of voltages and branch flows can be researched through the proposed ME-PLF calculation algorithm.

Fig. 8
figure 8

Probability density curves comparison considering dependence

Full size image

Figure 9 shows the PDF curves of V101 and P101-102 in Scenario 3 (20 × 1.5 WM WTs are integrated into bus 101). The proportional coefficient tan(δ) alters from − 0.3 to − 0.98. It can be seen that PDF shapes of voltage magnitudes have an evident change comparing with Scenario 2, while little on branch flows. ME-8th can give reasonable PDF results in this case while GC-8th is still inadequate to calculate the PDFs accurately.

Fig. 9
figure 9

Probability density curves comparison in Scenario 3

Full size image

4.2 Test case of Polish power system

A 2383-bus test case from MATPOWER platform [29] is utilized to confirm the effectiveness of ME-PLF calculation algorithm in a large scale power system. This system represents the 1999 to 2000 winter peak Polish power system. It has 2383 buses, 327 generators and 2898 branches [28]. In this case, buses 2015 and 1684 are both integrated with 25 × 1.5 MW WTs. The WT parameters are formed with Scenario 2. PDFs of V1640, V1684, V2015, P2015-1640, P2015-1684 and P1684-1683 are calculated by GC-8th and ME-4th to ME-8th.

Table 6 shows the ARMS for PDFs of voltage magnitudes and branch flows. Figure 10 shows the PDF curves of V2015, P2015-1640 and P2015-1684. From ARMS and PDF, it can be seen that both ME-PLF calculation algorithm and GC-PLF calculation method can keep a constant error level with the rapid increase of power system scale. The characteristics of PDFs for bus voltages and branch flows are similar with the study of IEEE 118-bus case. ME-PLF calculation algorithm still gives better results than GC-PLF calculation method.

Table 6 Comparison of cumulant estimation error results in Polish system
Full size table
Fig. 10
figure 10

Probability density curves comparison in Polish system

Full size image

For time consumptions, MCS takes up a tremendous time, with 4909 s, GC-8th and ME-8th are both efficient, with 27.07 s and 24.12 s, respectively. In conclusion, time consumption of ME-PLF calculation algorithm increases slowly with the enlargement of system scale without sacrificing the superiority of calculation accuracy against GC-PLF calculation method, which demonstrates the practicality of ME-PLF calculation algorithm facing PLF calculation of large scale power system.

5 Conclusion

This paper proposes a ME-PLF calculation algorithm for power system integrated with WPGs. The proposed algorithm gives more reliable and accurate PDF results of bus voltages and branch flows comparing with GC-PLF calculation method in scenarios with different WT parameters. It can also avoid gaining negative values in tail regions of PDFs, whereas GC-PLF calculation method sometimes mistakenly gets negative values. PLF can be analyzed along with the changes of wind power factor by the modified cumulant calculation framework, which can take into account linear dependence between active and reactive power injections of WPGs. The conclusions of case studies can be summarized as follows:

  1. 1)

    With the same input information, ME-PLF calculation algorithm can obtain preferable PDFs in scenarios with all WT parameters against GC-PLF calculation method. It is because ME-PLF calculation algorithm not only satisfies the given constraints, but also gives a maximally noncommittal estimation of the unknown information through maximizing the Shannon entropy of PDF. So the proposed method can output better PDFs than GC-PLF calculation method with limited information. Besides, it can also deal with the limitation of GC-PLF calculation method that mistakenly gaining negative values in tail regions of PDFs.

  2. 2)

    ME-6th is always more accurate than GC-8th, and is adequate to achieve accurate PDFs in most situations. For some special circumstances, ME-8th is requisite for a reliable outcome. Further increase of order of ME-PLF calculation algorithm will no longer improve the accuracy due to cumulant calculation errors introduced by linearized model of AC load flow.

  3. 3)

    Linear dependence between active and reactive power injections of WPGs greatly affects PDF curves of voltage magnitudes. When the proportional coefficient tan(δ) alters from -0.3 to -0.98, PDF curves of voltage magnitudes get wider. In another word, the voltage fluctuation becomes larger.

  4. 4)

    For large scale power systems, ME-PLF calculation algorithm keeps the advantage of gaining correct and accurate PDF against the GC-PLF calculation method. In addition, both ME-PLF calculation algorithm and GC-PLF calculation method consume much less time than MCS without sacrificing accuracy. Therefore, the proposed method is able to deal with PLF calculation for practical systems reliably and effectively.