A Novel Probabilistic Power Flow Algorithm Based on Principal Component Analysis and High-Dimensional Model Representation Techniques

Because the penetration level of renewable energy sources has increased rapidly in recent years, uncertainty in power system operation is gradually increasing. As an efficient tool for power system analysis under uncertainty, probabilistic power flow (PPF) is becoming increasingly important. The point-estimate method (PEM) is a well-known PPF algorithm. However, two significant defects limit the practical use of this method. One is that the PEM struggles to estimate high-order moments accurately; this defect makes it difficult for the PEM to describe the distribution of non-Gaussian output random variables (ORVs). The other is that the calculation burden is strongly related to the scale of input random variables (IRVs), which makes the PEM difficult to use in large-scale power systems. A novel approach based on principal component analysis (PCA) and high-dimensional model representation (HDMR) is proposed here to overcome the defects of the traditional PEM. PCA is applied to decrease the dimension scale of IRVs and eliminate correlations. HDMR is applied to estimate the moments of ORVs. Because HDMR considers the cooperative effects of IRVs, it has a significantly smaller estimation error for high-order moments in particular. Case studies show that the proposed method can achieve a better performance in terms of accuracy and efficiency than traditional PEM.


Introduction
In recent decades, the penetration level of renewable energy sources (RESs), such as wind and solar power, has increased rapidly. RESs have undispatchable and limited prediction characteristics, which brings challenges for power system operation. For a power system with large-scale RES integration, power generation becomes uncertain, which causes system state variables, such as line flows and bus voltage, to vary. Clearly, in this situation, the state variables should not be seen as having deterministic values, and the traditional deterministic power flow (DPF) is incapable of analyzing system states.
Instead of DPF, the concept of probabilistic power flow (PPF) emerges [1]. In PPF, some power injections, including RES generation and loads, are seen as input random variables (IRVs), while state variables are seen as output random variables (ORVs). Given the probability information of IRVs, the moments or probability distribution of state variables can be obtained by a PPF algorithm. PPF is an effective tool that can provide valuable information for applications, such as operational security assessment [2], planning and economic dispatch [3]. Therefore, the study of the PPF algorithm has seen increased attention.
The remainder of this paper is organized as follows: Section 2 presents the methodology of the proposed method, followed by the procedure in Section 3. In Section 4, case studies are examined. The proposed method is compared with MCS and Zhao's PEM. The influence of PCA and the adaptability of the proposed method are also analyzed in Section 4. In Section 5, a discussion of the proposed method is presented. Finally, conclusions are drawn in Section 6.

Problem Formation
The load flow equation shows the relationship between IRVs and ORVs, and the Alternating Current (AC) load flow equations are as follows where P i and Q i are the active and reactive power injections at bus i, respectively, V i and V j are the voltage magnitudes of bus i and j, respectively, θ ij is the voltage angle difference between bus i and j, and G ij and B ij are parameters related to the line that connects bus i and j. Clearly, the above AC load flow equations can be expressed with a qualitative equation where X is the vector of IRVs, including active and reactive power injections. Y is the vector of ORVs, including bus voltage magnitude V and angle θ and line active flows P and reactive flows Q, and G(·) represents the nonlinear AC power flow equations. Given the probability information of the IRVs, we want to know the moments of the ORVs. The lth central moment of the ORVs is as follows where E(·) is the expectation operator, l denotes the l th central moment, and f X (X) represents the joint probability density function (PDF) of X.
Here are some challenges for directly solving (3). The first is that the integration expressed in (3) is multi-dimensional; secondly, if X follows non-Gaussian probability distributions, the expression of f X (X) is complicated; moreover, the difficulty of integration is related to the non-linearity of G l (X). Proper transformation and dimensionality reduction are needed to perform this integration. In this paper, PCA and HDMR are combined to perform the integration.

Principal Component Analysis
PCA [23,24] is an efficient tool for data analysis, especially for high-dimensional data with potential correlations. The mechanism of PCA is orthogonal transformation. Through orthogonal transformation, a high-dimensional dataset is transformed to a low-dimensional dataset while maintaining its main characteristics. In this section, we introduce the procedure of PCA, as it can eliminates the correlations among IRVs, and improve the overall efficiency of proposed algorithm.
In this paper, the power injections of each bus are considered IRVs, specifically, the power injections of RES and load. The distribution model of IRVs can significantly influence the results of PPF, especially the distribution of RESs. However, distribution modelling itself is not the focus of this paper; moreover, the proposed method has good adaptability to any distribution model, thus, samples of IRVs can be derived from commonly used distributions or from historical data.
Suppose m power injection samples are available for each IRV and the dimension of the IRV vector is n; then, an n×m IRV sample matrix X can be constructed: The first step is centralization. By subtracting the mean value of each IRV µ X , the centred sample matrix X C can be obtained Suppose that the IRVs have linear correlations, and the covariance matrix Σ corresponding to X C is used to describe the correlation where σ 2 i is the variance of random variable i and σ 2 i,j is the covariance between random variables i and j. σ 2 i,j is calculated by where x k i and x k j are the kth samples of IRVs i and j, respectively, and µ i and µ j are the means of IRVs i and j, respectively.
Clearly, the covariance matrix is an n×n symmetric positive semi-definite matrix, and its eigenvalues can be calculated by det(Σ − λ i I) = 0 i= 1,2,· · · ,n (8) where det(.) is the determinant, λ i is an eigenvalue, and I is the identity matrix. Sorting the eigenvalues in descending order, we can calculate the corresponding eigenvectors e i by Thus, we can construct the projection matrix E E= [e 1 e 2 · · · e n ] Transforming the dataset X C via E, the PC matrix Z is obtained Z is an n×m matrix that has the same dimension as X or X C . The ith row of the matrix Z is the ith PC, and we use z i to express it. z i is a linear combination of IRVs that has no actual physical meaning.
Here are some characteristics of PCs: (1) Under the assumption that only linear correlation exists, the covariance between any two PCs is 0; therefore, correlation is eliminated; (2) The variance of z i is λ i . Clearly, the variance in elements in Z has a descending order. According to PCA, the larger the variance is, the more information is contained in a PC. Therefore, in Z, the importance of z i follows a descending order.
Clearly, matrix Z and matrix X have the same dimension: if all of the components in Z are retained, the scale of the PPF problem is not reduced. Thus, we must determine how many PCs should be retained, which can be achieved by calculating the PCs' cumulative contributions.
The contribution α i of ith PC is the proportion of its corresponding eigenvalue compared to all eigenvalues. The greater the proportion, the stronger the ability of the PC to integrate the information of the original IRVs. α i is calculated as follows The cumulative contribution of the first k components is calculated as The greater the value of M k , the more fully the first k PCs represent the information of the original IRVs. When choosing the retained PCs, M k should be compared to a threshold, such as M k ≥ 95%. Suppose the first k components of the projection matrix are chosen as E k = [e 1 , e 2 , · · · e k ], then the retained PC matrix is Z k = E T k X c . Z k is a k×m matrix, and it retains the most information of the original IRVs.
For a PEM-like algorithm, the calculation burden is closely related to the dimensions of the IRVs. By performing PCA on the IRV vector X, we can obtain an uncorrelated k×m matrix Z k , where k<<n. Clearly, the calculation burden is not directly related to the dimension of IRVs but relates to the dimension of the retained PC matrix. Using PCA before the moment estimation procedure, the total efficiency of the algorithm is significantly improved; moreover, since there is no correlation between PCs, the impact of the correlations between IRVs is handled well.
A PC is a linear combination of IRVs, and its distribution is arbitrary. Sometimes, the abscissas cannot be directly generated from an arbitrary probability distribution, as it may be outside the region in which the random variable is defined [16]. Here, we perform a Nataf transformation [25] on PCs. This transformation will connect the principal component space to a standard Gaussian space, which will facilitate the subsequent implementation of HDMR and ensure the validity of chosen abscissas.
Assuming that the marginal cumulative distribution function (CDF) where Φ(·) and Φ −1 (·) are the CDF and inverse CDF of the standard normal variable u i . Therefore, through PCA and marginal transformation, Y can be expressed as a function of the independent standard Gaussian vector U:

High-Dimensional Model Representation
HDMR [26][27][28] is an analysis tool for capturing the high-dimensional relationships between IRVs and ORVs. The core idea of HDMR is using several low-dimensional functions to approximate the original high-dimensional functions. Applying HDMR to (14), the function H(U) is split into a summation of several sub-function terms H 0 is a constant term that represents the mean response on Y. H i (u i ) is a first-order term that represents the effects of a single IRV u i on Y; H i 1 ,i 2 (u i 1 , u i 2 ) is a second-order term that represents the cooperative effects of the IRVs u i 1 , u i 2 on Y. Clearly, higher-order terms represent the cooperative effects of multiple IRVs on Y, and H i 1 i 2 ···i n (u i 1 , u i 2 , · · · u i n represents the cooperative effects of all IRVs The detailed expressions of the components in each term are as follows where U c is the reference vector; as we have transformed the original IRV space to the standard Gaussian space by PCA and marginal transformation, U c is [0, · · · , 0] here. U i,c is formulated after eliminating the sub-vectors of U c corresponding to u i ; therefore, H(u i , U i,c = H(0, · · · , 0, u i , 0, · · · , 0). The expression of the higher-order component can be obtained recursively. The expression and calculation process become complex after including the terms in (15). A practical method for doing this is to use a truncated form of (16); given a highest order for term s, only cooperative effects with no more than s variables are involved.
Considering terms up to first order, and combining (16) and (17), H(U) is expressed as Similarly, considering terms up to second order, H(U) is expressed as where R 2 and R 3 are the residue errors. Ignoring R 2 and R 3 , we obtain the first-order and second-order approximations of H(U). Clearly, the first-order approximation only considers a single variable's individual effect on Y, and all cooperative effects of IRVs are ignored, while the second-order approximation considers the cooperative effects of any two IRVs. It is worth mentioning that the first-order approximation is exactly the same as that used in Zhao's PEM scheme [15]; therefore, Zhao's PEM scheme could be considered as a special case of HDMR.
Since load flow equations are non-linear equations, ignoring the cooperative effects of IRVs will inevitably cause estimation errors; therefore, it is necessary to use a form of approximation that considers cooperative effects. According to our experience, ignoring the cooperative effects of IRVs will cause an unacceptable level of estimation error for high-order moments, which will seriously affect the accurate description of the ORVs' distribution, especially when the IRVs follow a non-Gaussian distribution: considering the cooperative effects of two IRVs can significantly eliminate this error and is sufficient in practical use.
Using the second-order approximation, the lth central moment of Y can be calculated by In (19), the expectation of the univariate and bivariate functions can be calculated by Gauss-Hermite quadrature where r is the number of abscissas, w GH,m is the weight and α GH,m is the abscissa. Table 1 shows the typical abscissas and their corresponding weights.

Solution Procedure
In this section, the solution procedure of the proposed method is summarized step by step as follows: 1. Prepare the necessary data, including samples of uncertainty sources (RES and load), which could be generated from distribution models or collected from historical data and the parameters of the power system; 2. Obtain the centralized sample matrix Xc and its corresponding covariance matrix Σ; 3. Obtain the projection matrix according to (8), (9) and (10); 4. Obtain the retained PC matrix Z k according to (11), (12), and (13) and the predetermined threshold value; 5. According to the dimension of Z k , choose abscissas in standard Gaussian space, then use (11) and (14) to project these abscissas into the original correlated IRV space; 6. Use the projected abscissas as IRVs and calculate the corresponding DPFs; 7. Calculate the moments of the ORVs from (20) and (21).

Case Study
The proposed method is tested on two modified test systems. The DLF program is based on MATPOWER 6.0 [29] and is executed on a PC with a 2.8-GHz CPU and 16 GB RAM.
The thirdorder moment and fourth-order moment, namely skewness and kurtosis, are important parameters to describe the shape of a PDF, we calculate the ORV's moments up to the fourth-order moment. The performance of the proposed method is compared with that of HDMR alone, Zhao's PEM proposed in [21], and 50000 iterations of MCS. Additionally, three contribution thresholds of PCA (M k ≥ 99%, M k ≥ 95% and M k ≥ 90% respectively) are set subjectively to measure the influence of the contribution threshold on the accuracy and efficiency of the proposed method. The number of abscissas of HDMR is chosen as 5, and, for consistency, that of Zhao's PEM is adjusted to be the same.
The average relative error index (AREI) [11] is adopted to indicate the accuracy of the moment of the ORVs estimated using each method. The average relative error index is defined as follows where i corresponds to the ith central moment; N is the number of types of ORV; j is the corresponding series number; and γ MC i,j and γ EV i,j refer to the results calculated by MCS and the evaluated method, respectively.

IEEE-30 Test System
A modified IEEE-30 [30] test system is used here to test the accuracy performance of the proposed method. This system contains six generators, 41 lines, and 20 loads.
Four wind farms are integrated in the test system. The probability distribution of wind farm output is assumed to follow the beta distribution, as this distribution is a proper choice to describe the short-term uncertainties of wind power [31], and the power factor of wind farms is assumed to be 0.95. The penetration level of wind power output is approximately 30%. The located bus numbers and uncertainty descriptions of wind farms are shown in Table 2. The uncertainty of the load at each bus is assumed to follow a Gaussian distribution, with the means of active and reactive power injections equal to those of the base case data and standard deviations equal to 5% of the means. Accordingly, the total number of IRVs is 48 (including the active and reactive injections of four wind farms and 20 loads). All IRVs are assumed to be independent of each other in this case.
The performance of the evaluated methods is shown in Table 3. ε 3 and ε 4 , which represent the estimation error of high-order moments of evaluated methods, are shown in Figure 1.    The data presented in Table 3 and Figure 1 show the following results: (1) HDMR has the minimal estimation error of the first four moments compared to other candidates, but at the cost of a much longer calculation time (more than 50 s); (2) The estimation accuracy of Zhao's PEM scheme for the first two orders of moments is satisfactory, but the estimation error is extremely high when estimating third-and fourth-order moments; (3) Combining PCA with HDMR can be very effective in improving computational efficiency. When the contribution threshold is set at 90%, the calculation time of PCA+HDMR is even shorter than that of the PEM; (4) Reducing the contribution threshold results in a decrease in computational accuracy. This is because some of the information useful for moment estimation is discarded by PCA. Therefore, in practice, the contribution threshold of PCA should be set carefully. When setting the contribution threshold, we need to balance the conflicting goals of precision and efficiency. The result shows that 99% is a good choice in this case.
Note that here, the estimation error of the third-order moment is relatively high for all evaluated methods. This is because the absolute value of the third-order moment is usually very close to 0; therefore, a small variation will result in a large relative error; however, this kind of large relative error has a minor influence on describing the PDF of a random variable.

Sensitivity to Correlations
In this section, the outputs of wind farms are assumed to be correlated. We will test the adaptability of the proposed method to different correlation scenarios. Three correlation scenarios are assumed in this section: low correlation, medium correlation and high correlation. The correlation matrix is as follows: As seen in Table 3, we found that a contribution threshold setting of 99% is good in terms of both accuracy and efficiency; therefore, the contribution threshold of these three scenarios is set at 99%.
The AREIs of the proposed method in the three correlation scenarios are shown in Table 4. ε 3 and ε 4 of these three scenarios are shown in Figure 2.
As the data show in Table 4 and Figure 2, in the three correlation scenarios, the estimated error of the proposed algorithm does not fluctuate greatly, so it can be concluded that the proposed method has good adaptability to different correlation scenarios.   As the data show in Table 4 and Figure 2, in the three correlation scenarios, the estimated error of the proposed algorithm does not fluctuate greatly, so it can be concluded that the proposed method has good adaptability to different correlation scenarios.

IEEE-118 Test System
In this section, the adaptability of the proposed method to large power systems is tested. The IEEE 118 system contains 54 generators, 186 lines, and 100 loads. Eight wind farms are integrated in this test system. Accordingly, the total number of IRVs is 216. The located bus numbers and uncertainty descriptions of the wind farms are shown in Table 5. The contribution threshold of PCA is set at 99%, and other assumptions and settings are the same as those of the IEEE-30 test system. The performance of the evaluated methods is shown in Table 6. The data presented in Table 6 show the following results: (1) Similar to the results of the IEEE-30 test system, HDMR has the best accuracy performance, followed by PCA+HDMR, and the PEM has the worst performance. Although HDMR has the highest calculation accuracy, its computational burden becomes unacceptable in large systems (up to 406 s in the IEEE-118 test system); (2) Zhao's PEM still cannot effectively estimate third-and fourth-order moments. Clearly, ignoring the cooperative effects of the IRVs significantly increases the estimation error of high-order moments; (3) PCA+HDMR performs high-order moment estimation with a small loss of accuracy. The calculation time in this test system is only 6.52 s. With a comprehensive consideration of calculation accuracy and calculation efficiency, PCA+HDMR is seen to be more suitable for practical use than the other methods.

Non-Linear Dependency
In this paper, dependency among IRVs are assumed to be linear. For some cases in which IRVs have non-linear dependencies, using PCA may cause additional errors, since non-linear dependencies cannot be fully eliminated by PCA. In this situation, kernel PCA (kPCA) [32] may be used instead of PCA to increase accuracy. The basic idea of kPCA is to map the data in the input space to a feature space via a non-linear map and then apply PCA in the feature space.

Probability Distribution Approximation
For practical engineering applications, knowing only the moments of the ORVs may not be sufficient: their PDF or CDF may also be of interest. Gram-Chalier series [10] and Cornish-Fisher series [11] are commonly used in the literature to generate PDF or CDF curves according to moment information. However, there are some defects in such series-based methods. One is that the terms of the series are hard to determine in advance, and the other is that these two series may generate illegitimate negative probability distribution functions [13]. Versatile distributions [33], the Johnson system [8] and Gaussian mixture models [3,34] are good candidates that can provide better performance than series-based methods.

Conclusions
To overcome the defects of the traditional PEM, a new PEM approach based on PCA and HDMR is proposed in this study. The proposed method uses PCA to decrease the dimension scale of IRVs and eliminate correlations among IRVs, which makes HDMR adaptable to power systems that have a large number of IRVs. HDMR is used to estimate ORV moments; in this manner, a defect of the traditional PEM, that it has difficulty estimating high-order moments, is eliminated. Case studies show that the proposed method can achieve better performance in terms of accuracy and efficiency than the PEM and MCS when considering a large number of correlated, non-Gaussian IRVs.
In the future, other methods, such as kPCA and parallel computing, can be studied and combined with the proposed method to further enhance its performance. The proposed method may also prove useful for other PPF-based applications.
Author Contributions: H.L. conceived and designed the study; this work was performed under the guidance of Z.Z. and X.Y. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.