Constant Jacobian Matrix-Based Stochastic Galerkin Method for Probabilistic Load Flow

: An intrusive spectral method of probabilistic load ﬂow (PLF) is proposed in the paper, which can handle the uncertainties arising from renewable energy integration.


Introduction
Load flow studies are essential for power system operators and planners, and accurate real-time load flow analysis is the basis of advanced application of modern EMS, such as optimal power flow, voltage control and so on.However, with the rapid increase in renewable energy (primarily wind and solar energy) which has strong volatility and intermittency, deterministic load flow (DLF) analysis cannot reflect the uncertainties of a power system properly.Therefore, the uncertainty of modern power grid needs to be qualified to analyze its influence on the power system.
Uncertainty quantification (UQ) [1] is the science of quantitative characterization and reduction of uncertainties in applications.It tries to determine how likely certain outcomes are if some aspects of the system are not exactly known.Probabilistic load flow (PLF) evaluation is an efficacious UQ tool to tackle the uncertainties in power systems and obtain the statistical quantities of random output variables.Traditional PLF methods can be classified into four categories: simulation methods, approximate techniques, analytical methods and heuristic procedures.
The Monte Carlo method (MCM) [2] is the most straightforward among simulation methods, but it requires DLF calculations at numerous sampling nodes even though some improvements have been proposed (such as Latin Hypercube [3] and Quasi-Monte Carlo [4] samplings).Approximate techniques provide approximate description of the statistical quantities of random output variables, and they include the point estimate method (PEM) [5], first-order second-moment method (FOSM) [6], non-parametric density estimate method [7], state variable method (such as polynomial transformation [3] and unscented transformation [8]).However, the accuracy of an approximate technique worsens as the order of estimated statistical quantities becomes higher.Analytical methods Energies 2016, 9, 153 2 of 18 include the FFT method [9], the cumulants method (based on Gram-Charlier [10,11] Cornish Fisher [12] and Legendre [13] expansion) etc.; they have a high calculation speed but require some algorithm simplifications which might introduce large errors.Heuristic methods mainly refer to the fuzzy logic method [14]; it can handle the uncertain variables without statistical information (such as probability distribution function (PDF) and cumulative distribution function (CDF)).
A spectral method based on generalized polynomial chaos (gPC) expansion [15][16][17][18][19][20][21] has become one of the most widely adopted UQ methods nowadays, which is efficient at uncertainty reduction and quantitative characterization for complex stochastic systems.The gPC-based methods represent system uncertain variables by truncated gPC series expansions, and gPC coefficients are computed by a stochastic collocation (SC) [18][19][20][21] or stochastic Galerkin (SG) approach [18,20].The SC method is non-intrusive, and it solves a set of decoupled equations by repeatedly calling an existing deterministic solver at sampling nodes and obtaining gPC coefficients by post processing.It is easy to apply but has inferior accuracy due to "aliasing errors" [16].The intrusive SG method forms coupled deterministic equations by Galerkin projection and computes gPC coefficients directly.By ensuring the residue of the Galerkin system equations to be orthogonal to the linear space spanned by gPC basis functions, SG method offers the optimal accuracy with the least number of equations among all the gPC-based methods [16].However, the derived Galerkin system equations are usually coupled and high dimensioned, which makes them cumbersome to solve or even unsolvable.
In this paper, a spectral method based on the high-precision SG method was presented to solve PLF, and the cumbersome SG method is modified by utilizing the P-Q decoupling properties of load flow.Thus the proposed method can achieve high calculation speed and accuracy.
This paper is organized as follows: the basic ideas of gPC expansion and the PLF Galerkin system equations constructed by Galerkin projection are introduced in Section 2. Section 3 presents the modified SG method to solve the coupled equations of the PLF Galerkin system, while detailed derivations are given in Appendix.Section 4 studies the performances of the proposed method.Finally, Section 5 concludes this paper.

gPC Expansion of Stochastic Variables
Define tφ i pξqu as the univariate gPC polynomial of random variable ξ P R, index i as the polynomial order.Taking as the set of D-variate gPC basis functions of muti-dimensional random variable ξ " rξ 1 , ..., ξ D s P R D , the order of ψÑ i pξq equals ˇˇˇÑ i ˇˇˇ" i 1 `. . .`iD , ψÑ i pξq can be defined as follows: The inner product of a random function is defined as: where x¨, ¨y denotes the inner product, and ρpξq denotes the PDF of variable ξ.Set where δ ij is the Kronecker delta function.The orthogonal property of ψÑ i pξq is as follows: Energies 2016, 9, 153 3 of 18 According to UQ theory [15][16][17], any R-valued random variable defined on a probability space pΘ, Σ, Pq can be well approximated by a truncated gPC series.Thus Spξqcan be approximated by an order-P gPC series.In Equation (7), which contains K gPC basis functions, ŝ| Replacing the index vector Ñ i by a scalar kpk P r1, Ksq: Equation ( 7) can be rewritten as: With ψ 1 pξq " 1, the mean and variance of Spξq can be approximated by: Similarly, other statistical quantities of Spξq can also be approximated by applying their definitions directly to the gPC approximation.
Generally in the modern power grid, the uncertain factors (such as fluctuations of wind speeds and loads, etc.) are mutually dependent, and the impactions of their correlations are significant and cannot be ignored.Setting the PDF of D-dimensioned random variable ξ " rξ 1 , . . ., ξ D s as ρpξq, and constructing the gPC basis function tψ k pξqu according to ρpξq can achieve the optimal approximation convergence rate; the detailed procedure can be seen in [16,17].The typical correspondence between the type of gPC basis function and their underlying random variables ξ is shown in Table 1.There are two ways to choose the gPC model.First, if all the random variables are the same type, we can directly choose the gPC basis, and Table 1 shows some common conditions.Second, if the types of random variables are different or not common, we can directly construct the gPC basis.The two methods have been studied deeply in [16,17] and will not be elaborated on here.In this paper, in order to facilitate the following study and description, inverse transformation method and Cholescky decomposition [22,23] are used to transform various correlated random inputs into independent Gaussian random variables.

Construct Galerkin System of PLF Based on gPC
The original non-linear load flow equations for a power system can be expressed as:

#
f P i pθ, V, P SP i q " 0 @i P S PV Y S PQ f Q i pθ, V, Q SP i q " 0 @i P S PQ (11) The uncertainties of the power system can be depicted by random variables ξ " rξ 1 , ..., ξ D s; random inputs such as Gaussian distributed power loads and Weibull distributed wind speeds are expressed as P l,i pξq, Q l,i pξq and vpξq.Wind power generation can be derived by applying the wind speed in the power curve P W pvq of a wind turbine generator (see Equation (31)).Therefore, the P-order gPC series of P SP i and Q SP i can be pre-calculated, and their gPC coefficients Pi,k , Qi,k are treated as known constants: The unknown random outputs can also be expressed as: Therefore, the following main task is to calculate the gPC coefficients θi,k , Vi,k by SG method.Apply Equations ( 12) and ( 13) in load flow Equation (11) to obtain the PLF model, which is also called the residual PLF equations in the SG method.Express the residual equations in a vector-matrix form: where N 1 " N PQ `NPV , N 2 " 2N PQ `NPV , vector F P " r f P 1 , . . ., f P N1 s T contains the active power balance equations of PQ and PV buses, vector F Q " r f Q 1 , . . ., f Q N PQ s T contains the reactive power balance equations of PQ buses.Define the matrices of coefficients θ, V as follows: Galerkin projection ensures residual Equation ( 14) to be orthogonal to the linear space spanned by gPC basis functions tψ k pξqu , forming a N 2 K dimensional PLF Galerkin system with G P P R N 1 ˆK, G Q P R N PV ˆK as follows: The inner product in the right hand of Equation ( 16) is actually multivariate stochastic integration.Numerical quadrature such as Gaussian quadrature can be used to effectively evaluate the integral by a weighted sum of function values at specified grid points, that is: where Θ D " " ξ 1 , . . ., ξ N ı and W " rW 1 , . . ., W N s are the quadrature nodes and weights that constructed by the cubature rule.
Among all the cubature rules, the tensor product rule is used to construct a full-grid point Θ D and its corresponding weights W [24], and the number of quadrature nodes equals N " pP `1q D which increases exponentially with the dimensionality D of the integral.Thus, improved methods are proposed, such as a nested grid (with N " ř P i"1 2 i pD ´1 `iq!{pD ´1q!i!) [17], testing nodes (with N " K) [20] and sparse grid [16,25].The details of these cubature rules for constructing quadrature nodes are out of the scope of this paper.This paper applies the widely used sparse-grid technique, and software for generating sparse grids is available online [26].The sparse-grid technique uses a subset of the points from the tensor product rule and rescales the weight values appropriately; it can significantly decrease the computation expense especially when the dimensionality D is much larger than 1: Set ψ n k " ψ k pξ n q, we can derive a constant matrix ψ by calculating the value of gPC basis function tψ k pξqu K k"1 at N sets of quadrature nodes: Set are the ith elements in the vectors F n P and F n Q respectively.Approximating Equation ( 16) by Gaussian quadrature, the PLF Galerkin system can be obtained: Reshape matrices G P , G Q , θ and V into vectors: where reshpA, m, nq reshape the m-by-n matrix whose elements are taken column-wise from matrix A.

Solve the PLF Galerkin System in a Decoupled Manner
Starting from an initial guess of θ0 and V0 , θ, V are computed by using the Newton iteration method to solve the Galerkin system Equation ( 20): V ˜s, and Jacobian matrix J has the following structure: Set θ n j " θ j pξ n q, V n j " V j pξ n q, and if i, j P r1, N 1 s, the pk 1 , k 2 qth element in the submatrix J i,j P R KˆK can be calculated as follows, and more detailed derivation of J is given in the Appendix.
Energies 2016, 9, 153 6 of 18 At each iteration, in J need to be evaluated at N quadrature nodes, and their sum needs to be calculated for each element in J. Therefore, the Newton's iteration of SG is very expensive.Even if the power system is small and J is sparse, it can still be unsolvable.For example, Figure 1 shows the structure of J in the IEEE 14-bus test system; if D " 4, P " 3, K " 35, there will be 124 sub-blocks and 65,920 non-zeros in J, and p65920 ˆN ˆLq calculation times are needed to renew J if L iterations are used.
,  J R i j can be calculated as follows, and more detailed derivation of J is given in the Appendix.
f in J need to be evaluated at N quadrature nodes, and their sum needs to be calculated for each element in J .Therefore, the Newton's iteration of SG is very expensive.Even if the power system is small and J is sparse, it can still be unsolvable.For example, Figure 1 shows the structure of J in the IEEE 14-bus test system; if  , there will be 124 sub-blocks and 65,920 non-zeros in , J and (65920 N L)   calculation times are needed to renew J if L iterations are used.Utilizing the P-Q decoupled characteristic in power flow calculation, this paper forms constant Jacobian matrices ' J and '' J to avoid the regeneration of Jacobian, the calculation speed will be significantly improved.
' B , '' B are the constant Jacobian matrices in the traditional P-Q fast decoupled load flow method [27], which only contains network admittances: ; 0; ; 0 Applying Equation ( 24) in the derivation of J above, J changed into two constant matrices which only need to be calculated once (see the Appendix).Figure 2 shows the structure of matrices ' J and '' J after the matrix J in Figure 1 is decoupled; there will be 6615 and 2565 nonzeros in ' J and '' J , and (9180 N)  calculation times are needed, which is far less than before.Utilizing the P-Q decoupled characteristic in power flow calculation, this paper forms constant Jacobian matrices J 1 and J 2 to avoid the regeneration of Jacobian, the calculation speed will be significantly improved.B 1 , B 2 are the constant Jacobian matrices in the traditional P-Q fast decoupled load flow method [27], which only contains network admittances: Applying Equation (24) in the derivation of J above, J changed into two constant matrices J 1 P R N 2 ˆN2 , J 2 P R N 1 ˆN1 which only need to be calculated once (see the Appendix).Figure 2 shows the structure of matrices J 1 and J 2 after the matrix J in Figure 1 is decoupled; there will be 6615 and 2565 nonzeros in J 1 and J 2 , and p9180 ˆNq calculation times are needed, which is far less than before.For ease of expression, N constant matrices Thus ' J and '' J can be expressed as: For ease of expression, N constant matrices rΦ 1 , . . ., Φ N s are constructed at N quadrature nodes: Φ n " W n ψpn, :q T ψpn, :q P R KˆK , n P r1, Ns Energies 2016, 9, 153 7 of 18 Thus J 1 and J 2 can be expressed as: V ˜can be computed directly by Newton iteration method: Solve : Update : The proposed PLF method (referred to as SG-PLF method) only needs to compute deterministic PLF Galerkin system Equation ( 20) and construct constant Jacobian matrices Equation ( 26) at limited quadrature nodes, then solve Equation ( 20) by the Newton iteration method only once to reach high accuracy.The flow chart of the SG-PLF method is depicted by Figure 3, which consists of two parts.

Comparison with other PLF Methods
This section compares the calculation procedures of MCM, SC method (SCM) and the proposed modified SG method (SGM); see Figure 4.
MCM and SCM are both non-intrusive methods, they both start from the original load flow Equation ( 11) without using gPC approximation a priori, and compute the deterministic solutions at Firstly, based on the system basic data and PDF of input random variables, we form the PLF model Equation ( 14) based on gPC expansion of random variables and transform it into Equation ( 16) by Galerkin projection, then further approximate Equation ( 16) by Gaussian quadrature, and finally obtain the PLF Galerkin system Equation (20).Secondly, we solve the Galerkin system for coefficients θ, V directly by a Newton iteration method, and the gPC expansion of random output variables can be derived easily.

Comparison with other PLF Methods
This section compares the calculation procedures of MCM, SC method (SCM) and the proposed modified SG method (SGM); see Figure 4.
MCM and SCM are both non-intrusive methods, they both start from the original load flow Equation ( 11) without using gPC approximation a priori, and compute the deterministic solutions at a set of nodes.The main difference of MCM and SCM lies in how to select the nodes.MCM draws sampling nodes randomly according to the PDF of ξ and is uncertainty-dimension independent; the MCM may be the only choice when the number of random variables to be dealt with is too high (in ultra-high-dimensional problem).However, SCM constructs collocation nodes with a nodes generator which uses cubature rules such as tensor product, nested grid or sparse grid techniques.After repeatedly simulating deterministic load flow equations, MCM provides the statistical information by gathering and handling the resulting outputs, whereas SCM reconstructs the gPC coefficients in a post-processing step such as numerical integration.For example, SCM constructs and weights rW 1 , . . ., W N 2 s of random output variable Spξq (in Equation ( 8)), and then the gPC coefficients can be computed by: Energies 2016, 9, 153 9 of 19 (in Equation ( 8)), and then the gPC coefficients can be computed by: The SGM is an intrusive method as it directly computes the gPC coefficients by simulating a larger-size coupled PLF only once.With gPC approximation, it starts from the residual function Equation (14).The proposed SG-PLF method constructs nodes by sparse grid technique and modified the SGM by utilizing the P-Q decoupling properties of power system, which can avoid the heavy computation burden of forming Jacobian matrix in each iteration.Figure 5 shows a detailed illustration of the modified calculation model of the proposed SG-PLF method.ψ ', '' J J The SGM is an intrusive method as it directly computes the gPC coefficients by simulating a larger-size coupled PLF only once.With gPC approximation, it starts from the residual function Equation (14).The proposed SG-PLF method constructs nodes by sparse grid technique and modified the SGM by utilizing the P-Q decoupling properties of power system, which can avoid the heavy computation burden of forming Jacobian matrix in each iteration.Figure 5 shows a detailed illustration of the modified calculation model of the proposed SG-PLF method.The SGM is an intrusive method as it directly computes the gPC coefficients by simulating a larger-size coupled PLF only once.With gPC approximation, it starts from the residual function Equation (14).The proposed SG-PLF method nodes by sparse grid technique and modified the SGM by utilizing the P-Q decoupling properties of power system, which can avoid the heavy computation burden of forming Jacobian matrix in each iteration.Figure 5 shows a detailed illustration of the modified calculation model of the proposed SG-PLF method.ψ ', '' J J  or po t t   , so there are no significant difference between the total computational time of SCM and SGM , this can be verified in the following case study.If the original load flow Equation ( 11) are solved by P-Q fast decouple method, and the time cost of solving is t or .then the total time cost by MCM is pt or ˆN3 q, where N 3 is much larger than N 1 and N 2 .If the node generators in SCM and SGM are both using the sparse grid technique, where N 1 " N 2 " N, set the time for constructing nodes is t no and the time cost of post-processing is t po .Thus, the total time cost of SCM is pt no `tor ˆN `tpo q.In SGM, the times needed for constructing constant Jacobian matrices are t cj and for solving the modified model are t mo , so the total time required for SGM is pt no `tmo `tcj q.
The dimension of modified PLF model in SGM is K ˆN times the size of the original PLF model, and the equations need to be solved in SCM is N times the size of the original PLF model.t cj can be ignored, and t mo is slightly larger than t or ˆN `tpo , so there are no significant difference between the total computational time of SCM and SGM , this can be verified in the following case study.
In conclusion, on the one hand, since N 3 is much larger than N 1 and N 2 , MCM needs numerous deterministic load flow calculations to obtain accurate results, the SG-PLF method has a distinct speed advantage over MCM.On the other hand, SCM only requires no errors at each collocation node and will introduce aliasing errors caused by the introduction of the nodal sets, whereas the SG-PLF method ensures that the residue of the PLF Equation ( 14) is orthogonal to the linear space spanned by the gPC basis functions, as in Equation ( 16).In this sense, the accuracy of SGM is optimal.
Therefore, the proposed SG-PLF method has accuracy advantage over the non-intrusive SCM and speed advantage over the MCM.

Verification Index
In order to evaluate the computational accuracy and efficiency of SG-PLF method, a series of PLF studies are carried out on the modified IEEE14-bus and IEEE 118-bus test systems.The base case data of these two test systems are available in [28].Two types of input random variables are considered in this case study, load and wind power included.Relative error ε (RE) [5] and average root mean square (ARMS) [12] error τ are computed using the MCM results of 30,000 trials as reference.RE is defined as: The ARMS is defined as: where x refers to the type of random output variables (V, θ and P line etc.), s may refer to the mean µ and standard deviation σ. s x M denotes the result calculated by MCM and is taken as the reference value, s x G is the result obtained by SG-PLF method.ε x s,mean and ε x s,max are the average and maximum values of ε x s .Choosing T sampling nodes of variable ξ arbitrarily: rξ 1 , . . ., ξ T s, X Gt and X Mt are the outputs of the SG-PLF method and MCM at the tth node ξ t pt P r1, Tsq.

IEEE 14-Bus Test System
The case considered in this section is the IEEE 14-bus test system which is modified to include wind generation.The uncertainties of forecasted loads are modeled as normal distribution with σ 5% Energies 2016, 9, 153 10 of 18 to µ, whose means equal the base case data.The load correlation matrix between bus 4 and bus 5 is ρ l .Two wind farms are located at bus 13 and bus 14, and their wind speed correlation matrix is ρ v : The wind power is mainly affected by wind speed distribution, power curve and control strategy of wind turbine generator.The wind speed is assumed to follow the Weibull distribution with scale parameter as 10.7 and shape parameter as 3.97.The wind speed is transformed into wind power according to the power curve of wind turbine generator P W pvq [3]: the cut-in, rated and cut-off speeds are v ci " 4 m{s, v r " 15 m{s, v co " 25 m{s respectively, and the rated power P Wr is 0.7 p.u.The wind turbine generator is controlled by constant power factor strategy.

Impact of the gPC Expansion Order
Figure 6 shows the CDFs obtained by SG-PLF method under different gPC expansion order P (defined in Equation ( 7)).The results shows that the order-3 SG-PLF method provides good fitting to the CDF obtained by MCM, gives an intuitive confirmation that the order-3 SG-PLF method has the ability to well handle over the PLF problem as almost the same accuracy as the MCM.

Impact of the gPC Expansion Order
Figure 6 shows the CDFs obtained by SG-PLF method under different gPC expansion order P (defined in Equation ( 7)).The results shows that the order-3 SG-PLF method provides good fitting to the CDF obtained by MCM, gives an intuitive confirmation that the order-3 SG-PLF method has the ability to well handle over the PLF problem as almost the same accuracy as the MCM.2, which verified the good computing performance of the order-3 SG-PLF method.
Therefore, the gPC expansion order P can affect the calculation performance of the proposed method: when P is too low, the gPC expansion approximation Equation ( 8) will induce large truncation error; when P is too high, the number of quadrature nodes N will increase significantly, so it will lead to increasing of calculation but no progress in computation accuracy, because all the high order gPC coefficients in gPC expansion approximation Equation (8) will be zeros.In this paper, the expansion order is chosen by numerical experiments, but direct methods to determine the order are to be further studied in our future work.Furthermore, the resulted ε γ s,mean and ε γ s,max under different order P are shown in Table 2, which verified the good computing performance of the order-3 SG-PLF method.Therefore, the gPC expansion order P can affect the calculation performance of the proposed method: when P is too low, the gPC expansion approximation Equation (8) will induce large truncation error; when P is too high, the number of quadrature nodes N will increase significantly, so it will lead to increasing of calculation but no progress in computation accuracy, because all the high order gPC coefficients in gPC expansion approximation Equation (8) will be zeros.In this paper, the expansion order is chosen by numerical experiments, but direct methods to determine the order are to be further studied in our future work.
To further investigate the calculation efficiency of using the SGM to solve PLF, the calculation time and nodes number N of SGM, SCM and MCM are compared in Table 3, in which, the gPC-based SCM and SGM both construct the nodes according to the sparse grid technique, and the time cost with different P are compared.Figure 7 shows the PDFs and CDFs of the buses whose power injection have large fluctuations, and the results corresponding to the order-3 SG-PLF method are almost not visible because they coincide with those of MCM.time cost with different P are compared.Figure 7 shows the PDFs and CDFs of the buses whose power injection have large fluctuations, and the results corresponding to the order-3 SG-PLF method are almost not visible because they coincide with those of MCM.
(a) (b)  From the results above, the order-3 SG-PLF method can maintain high accuracy with feasible computation expense.This conclusion is further confirmed by comparing the ARMS results of the order-3 SG-PLF method in Table 4. Thus in the following study, the order-3 SG-PLF method is used for PLF investigation and verification.Take the frequency table of 4 V as an example, set the results of MCM as the reference, and further compare the calculation accuracy of order-3 SCM and order-3 SGM, shown in Table 5: From the results above, the order-3 SG-PLF method can maintain high accuracy with feasible computation expense.This conclusion is further confirmed by comparing the ARMS results of the order-3 SG-PLF method in Table 4. Thus in the following study, the order-3 SG-PLF method is used for PLF investigation and verification.Take the frequency table of V 4 as an example, set the results of MCM as the reference, and further compare the calculation accuracy of order-3 SCM and order-3 SGM, shown in Table 5: In this section, the results above show that: (1) the MCM can achieve the highest accuracy, but it needs the longest calculation time; and (2) although the computational speed of SGM is slightly slower than SCM, the intrusive SGM has obvious precision advantage than SCM (which is also analyzed in Section 4).
Therefore, through comprehensive consideration of calculation speed and the proposed SG-PLF method is the optimal among three methods.

Impact of the Fluctuation of Power Injection
In order to appraise how the variation level of uncertain variable impacts on the calculation performance of the proposed order-3 SG-PLF method, Table 6 shows the ε γ s,mean and ε γ s,max of proposed method under different load fluctuations at bus 4 (set the ratio of σ to µ is from 10% to 40%).From the results, it can be inferred that the variance level of power injection has no discernible effect on the calculation accuracy of the proposed method.

Impact of Wind Speed Correlation
Choose the statistical moments of active line power (mean µ-P line and standard deviation σ-P line ) as the research objects in this section.The uncertain wind powers of two wind farms are directly injected into both the end of lines 13-14, while lines 1-2 have no direct connection to wind power.
Under different wind speed correlations (WSC) between two wind farms, Figures 8 and 9 present the calculated mean and standard deviation (SD) values of the REs of all buses and illustrate how the calculation accuracy of the proposed 3-order SG-PLF method will change with the WSC growing from 0.1 to 0.9.The results show that WSC has a minor effect on the calculation accuracy of the proposed method.Figures 10 and 11 present the values of µ-P line and σ-P line with WSC growing from 0.1 to 0.9.Results show that, at different sites, the WSC has different impact on system condition.µ-P 13´14 remains basically unchanged but σ-P 13´14 is proportional to WSC, while µ-P 1´2 , σ-P 1´2 are all proportional to WSC.

IEEE 118-Bus Test System
In this section, PLF problem is solved for the IEEE 118-bus test system with the purpose of assessing how an increase in the uncertainty level (dimension and variance of random variables)

IEEE 118-Bus Test System
In this section, PLF problem is solved for the IEEE 118-bus test system with the purpose of assessing how an increase in the uncertainty level (dimension and variance of random variables)

IEEE 118-Bus Test System
In this section, PLF problem is solved for the IEEE 118-bus test system with the purpose of assessing how an increase in the uncertainty level (dimension and variance of random variables)

IEEE 118-Bus Test System
In this section, PLF problem is solved for the IEEE 118-bus test system with the purpose of assessing how an increase in the uncertainty level (dimension and variance of random variables)

IEEE 118-Bus Test System
In this section, PLF problem is solved for the IEEE 118-bus test system with the purpose of assessing how an increase in the uncertainty level (dimension and variance of random variables) would affect the performance of proposed method.From bus 24 to bus 27, each bus is connected to a wind farm, the rated power of each wind turbine generator P Wr is 0.2 p.u.The random active loads at bus 60, bus 78, 79 and bus 82 are assumed to follow Gaussian distribution.The variance of the random loads are larger than the previous case, with σ is 10% to µ.The correlation matrices of wind speeds and loads are in Equation (33).Other parameters of wind farms are consistent with the previous case.
To solve this case, the calculation time of MCM is 16.436 s while the order-3 SG-PLF method only takes 3.306 s.The RE and ARMS results of the proposed method are shown in Table 7, where the maximum RE values of V and θ appeared at bus 26 and bus 49. Figure 12 shows the resulting CDFs and PDFs of V 26 and θ 49 , which proved the calculation accuracy of the proposed order-3 SG-PLF method again.Therefore, when the variance and dimension of uncertain variables rise, the proposed method can maintain high accuracy whilst providing significant computational savings.
To solve this case, the calculation time of MCM is 16.436 s while the order-3 SG-PLF method only takes 3.306 s.The RE and ARMS results of the proposed method are shown in Table 7, where the maximum RE values of V and  appeared at bus 26 and bus 49. Figure 12 shows the resulting CDFs and PDFs of 26 V and 49 ,  which proved the calculation accuracy of the proposed order-3 SG-PLF method again.Therefore, when the variance and dimension of uncertain variables rise, the proposed method can maintain high accuracy whilst providing significant computational savings.Based on all of the previous simulation results, the details of some influencing factors on the performance of the proposed SG-PLF method are summarized in Table 8, where  denotes increase,  means decrease,  presents minor impact.Based on all of the previous simulation results, the details of some influencing factors on the performance of the proposed SG-PLF method are summarized in Table 8, where Ò denotes increase, Ó means decrease, Ñ presents minor impact.

Conclusions
An intrusive spectral method was proposed in this paper to solve the PLF problem, and simulations in IEEE 14-bus and test systems demonstrate its accuracy, efficiency and practicability.The main characteristics of the proposed method are as follows: 1.
High accuracy.The accuracy of proposed method mainly depends on the gPC expansion approximation Equation (7); by properly choosing the type of gPC basis function and the expansion order P, the proposed methods can offer the optimal accuracy with the least number of equations among all the gPC-based methods.

2.
Low computation burden.By fully utilizing the P-Q decoupling properties of load flow, the proposed method formed two constant Jacobian matrices in Newton's iterations, which reduces the computation cost significantly.

3.
Comprehensibility the results.The gPC expansion can efficiently handle various random variables and form a functional relationship between the system random inputs and outputs, which is also known as the uncertainty propagation in UQ theory.Thus it can efficiently illustrate the effect of the random inputs on the outputs in PLF problem.
Since the future operating conditions of power system involve more uncertainty due to the growing penetration of distributed generation and renewable energy sources, the proposed method may be increasingly desirable.We believe that it would encourage more interest in UQ research for power systems, with more sophisticated uncertainties and less expenditures of calculation time and resources.

Figure 1 .
Figure 1.Structure of Jacobian matrix for an IEEE 14-bus test system with D = 4 before decoupling.

Figure 1 .
Figure 1.Structure of Jacobian matrix for an IEEE 14-bus test system with D = 4 before decoupling.

Figure 2 .
Figure 2. Structure of Jacobian matrix for IEEE 14-bus test system with D = 4 after decoupling.

Figure 2 .
Figure 2. Structure of Jacobian matrix for IEEE 14-bus test system with D = 4 after decoupling.

Figure 3 .
Figure 3. Flow chart of the proposed SG-PLF method.

Figure 3 .
Figure 3. Flow chart of the proposed SG-PLF method.

Figure 4 .
Figure 4. Calculation procedure of PLF by SGM, SCM and MCM.

Figure 4 .
Figure 4. Calculation procedure of PLF by SGM, SCM and MCM.
order SG-PLF method under different WSC.
order SG-PLF under different WSC.
order SG-PLF under different WSC.
order SG-PLF method under different WSC.
order SG-PLF under different WSC.
order SG-PLF under different WSC.
order SG-PLF method under different WSC.
order SG-PLF under different WSC.
order SG-PLF under different WSC.
order SG-PLF method under different WSC.
order SG-PLF under different WSC.
order SG-PLF under different WSC.

Table 1 .
Correspondence of gPC and their underlying random variables.
The dimension of modified PLF model in SGM is K N  times the size of the original PLF model, and the equations need to be solved in SCM is N times the size of the original PLF model.
cj t can be ignored, and mo t is slightly larger than N

Table 2 .
Mean and maximum of ER of SG-PLF under 1-3 order.

Table 2 .
Mean and maximum of ER of SG-PLF under 1-3 order.

Table 4 .
Mean and maximum of ARMS of order-3 SG-PLF method.

Table 4 .
Mean and maximum of ARMS of order-3 SG-PLF method.

Table 5 .
The frequency table of different methods.

Table 6 .
Mean and maximum of RE of order-3 SG-PLF with different load variance.

Table 7 .
Mean and maximum ERs and ARMSs of 3-order SG-PLF.affect the performance of proposed method.From bus 24 to bus 27, each bus is connected to a wind farm, the rated power of each wind turbine generator Wr P is 0.2 p.u.The random active loads at bus 60, bus 78, bus 79 and bus 82 are assumed to follow Gaussian distribution.The variance of the random loads are larger than the previous case, with  is 10% to  .The correlation matrices of wind speeds and loads are in Equation (33).Other parameters of wind farms are consistent with the previous case. would

Table 7 .
Mean and maximum ERs and ARMSs of 3-order SG-PLF.

Table 8 .
Summary on the factors influencing the performance of the proposed SG-PLF method.