Probabilistic Energy Flow Calculation through the Nataf Transformation and Point Estimation

Featured Application: The work presented in this paper provides methods for probabilistic energy ﬂow calculations and will serve as a reference for the operation and maintenance of integrated energy systems. Abstract: With the increasing capacity of renewable energy sources, uncertainties regarding renewable energy and other dynamic loads in integrated energy systems (IESs) are increasing. Thus, it is necessary to study the probabilistic energy ﬂow (PEF) of IESs. However, existing PEF calculation methods such as the point estimate method (PEM) are computationally ine ﬃ cient when there are many random variables and estimated points; moreover, relatively large errors can occur when the estimated points are outside their limits. Hence, this paper presents a calculation method that addresses these problems. Because there are correlations among the variables, the Nataf transformation is employed to control the correlation quickly and e ﬀ ectively. A model for an IES that is interconnected with natural gas and electricity systems and accounts for the uncertainties of wind plants, photovoltaic power plants, and dynamic gas loads is presented. Correlations between wind plants and photovoltaic power plants are handled using the Nataf transformation. Finally, a modiﬁed PEM is developed to solve the PEF. For situations in which the estimated points exceed their boundaries, the power transformation and equal constraint transformation methods are used. The results of time-domain simulations demonstrate the e ﬀ ectiveness of the proposed approach.


Introduction
The increasing demand for energy has resulted in increased concerns regarding energy sources, and the exploration of new energy supply modes has attracted considerable interest. The concept of integrated energy systems (IESs) challenges that of traditional energy systems. IESs, which can fully utilize energy in various forms and optimize resource allocation, have rapidly developed in recent years [1]. Compared with traditional energy systems, these systems have significant advantages [2,3]. For example, all types of energy are closely coupled, which makes them complementary and mutually beneficial, and they can realize cascade utilization and collaborative optimization of energy. Presently, research on IESs is mostly focused on three aspects: energy flow analysis and calculation [4], multi-energy system optimization scheduling [2,3], and multi-energy system planning and multi-market games [5,6].
The ability to integrate various types of energy structures is a compelling advantage of IESs [7,8]. However, owing to the significant amount of renewable energy sources (RESs) involved, crucial Currently, the approximation method is an important approach to solving the PPF problem. This method mainly includes the first-order second-moment method (FOSMM) [24] and point estimate method (PEM) [25]. The FOSMM is fast and precise. Because it adopts the Taylor series expansion method, only the linear term is considered and the higher-order term is ignored; therefore, when the system is highly nonlinear or the random variables are very asymmetrical, a large error will be caused. PEM is now widely used for solving the PPF problem; through Gaussian interpolation, the estimated points can be easily obtained, resulting in fewer power flow calculations. PEM is very fast. However, the original PEM cannot deal with correlated variables [26] and the estimated points might exceed the estimator boundary constraints, which results in a less accurate result [27]. Therefore, PEM is restricted to specific applications.
However, studies examining cases wherein a great number of probability variables with different marginal distributions co-exist in IESs are scarce. In practice, the probability variables, such as those of RESs, are highly correlated with each other; for example, wind speeds and solar radiation at different sites within the same geographical area are typically dependent owing to similar weather conditions. Thus, large errors will occur if we ignore the correlations among different probability variables. The commonly used approaches to deal with the correlation of random variables are transformation methods such as the third-order polynomial normal transformation (TPNT) [28], rotational transformation and orthogonal transformation (RTOT) [29], Rosenblatt transformation [30], unscented transformation [26], and copula function method (CFM) [20].
TPNT transfers variables into the third-order polynomial normal space, which causes mapping errors and low accuracy in some cases. RTOT transfers the correlation variables from the linear correlation coefficient matrix into the generated samples. Unlike the unscented transformation, which does not require a linear transformation and is suitable for any arbitrary nonlinear function, RTOT is essentially a linear transformation method. CFM uses marginal distribution functions to obtain joint distribution functions; this is a convenient method for adding relevance among variables with good precision, so it is now widely adopted. The Nataf transformation [31], which is a type of Gaussian copula, is one of the most important CFMs and famous for its high precision and speed; it can link any distribution space to the normal distribution space, and this is extremely beneficial for the PEM. The estimated points and weights can be easily obtained through the Gauss-Hermite integration. Moreover, because the calculation of PEF is extremely time consuming, the PEM can considerably reduce the computation time. Therefore, the combination of the Nataf transformation with PEM is an excellent choice for PEF calculations. However, the probability variables can be bounded by constraints. For instance, wind farms have cut-in and cut-out wind speed constraints, and when the estimated points are outside the range of these constraints, they lose their physical significance and cause errors. In particular, when dealing with 2-point and 3-point estimation situations, because of the small number of estimated points, the errors can become considerable. Hence, it is necessary to limit the estimated points to remain within the boundaries of the constraints.
The goal of this study is to solve the PEF in an IES. We comprehensively consider the uncertainty regarding wind energy, solar energy in power systems, and gas loads in natural gas systems to establish a model for calculating the PEF. Considering the correlation between different wind-driven generators and solar energy generators in reality, we propose using the Nataf transformation to solve the correlation problem accurately and efficiently. As is well known, solving the PEF in an IES is more time consuming than solving the PPF in a power system. The Monte Carlo method is an inefficient solution to the PEF; hence, the method developed in our study improves the solution efficiency using the multipoint estimation method (MPEM). This paper proposes an extension to the existing MPEM that reduces the computational complexity and cost. For cases where the estimated points are outside of the range, the approach proposed in this paper uses the power transformation method and equality constraint transformation to improve 2-point and 3-point estimations and reduce the calculation error.
The remainder of this paper is organized as follows: A PEF model for IESs is described in Section 2. A PEM based on the Nataf transformation is detailed in Section 3, and 2-point estimation methods considering boundary constraints are discussed in Section 4. Section 5 presents a case discussing PEF calculations, and Section 6 discusses the ideas proposed in this paper and directions of research worth studying in future. The concluding remarks are given in Section 7. Notations and a list of abbreviations are provided in the Appendix A.

SSEF Calculation in an IES
The SSEF calculation is the basic step in PEF calculations; similar to the power flow model of power systems, the SSEF model can be expressed using the following equation: Here, P s i and Q s i are respectively the active and reactive power that the source injects into node i, P l i and Q l i are respectively the active and reactive power demanded by the ith node load, P c i and Q c i are respectively the active and reactive power consumed by the compressor, and P g i and Q g i are respectively the active and reactive power that the gas generator injects into node i. In addition, n represents the number of nodes directly connected to the ith node, Q i is the injection of natural gas at the ith node, and f im and f in represent the injections at the downstream and upstream nodes, respectively. The gas consumption of the compressor is denoted by F j : when the compressor takes air from node i, the correlation coefficient G ij is 1; otherwise, G ij is 0. If the energy required by the compressor comes from the power system, then G ij is also 0. Finally, p, q and f NG represent the active, reactive, and gas flow unbalances in a node.
The specific calculation of the SSEF and parameters for a natural gas network are discussed elsewhere [10]. This study uses separate Newton iterations to solve the SSEF problem. That is, the SSEFs of the power and natural gas systems are each iterated independently using Newton's method until convergence is achieved. Then, the unbalanced quantity between these energy systems is iterated through the energy coupling system unit until overall convergence is achieved, as shown in Figure 1.
The main steps of SSEF computation are as follows: Step 1: Initialize the variables with the information of the natural gas network and construct its correlation matrix.
Step 2: Using the initial pressure of the natural gas network, solve each variable of the natural gas network (such as the pressure of each node and the gas flow of each pipeline) using Newton's method.
Step 3: Compute the value of the energy interaction of the natural gas network with the power system in the coupling system, which consists of the gas turbine and compressor.
Step 4: Construct the nodal admittance matrix of the power system and compute its initial values from the energy interaction obtained in Step 3. Solve the power flow of the power system using Newton's method and obtain its information, which includes node voltage, active power, reactive power, and voltage phase angle.
Appl. Sci. 2019, 9, 3291 5 of 23 Step 5: Compute the energy interaction value again from the results of Step 4 and then use it to compute the natural gas network information using Newton's method. Repeat Steps 2-5 until the relative error of the current and previous energy interaction values meets the accuracy requirement (e.g., less than a given tolerance).
Appl. Sci. 2019, 9, 3291 5 of 25 Step 5: Compute the energy interaction value again from the results of Step 4 and then use it to compute the natural gas network information using Newton's method. Repeat Steps 2-5 until the relative error of the current and previous energy interaction values meets the accuracy requirement (e.g., less than a given tolerance).  Figure 1. Steady-state energy flow calculation in an integrated energy systems.

Uncertainty Analysis and Uncertainty Modeling
With the rapid development of new energy sources in recent years, the capacity for new energy and power generation in IESs is increasing quickly, which causes considerable uncertainty in these systems. At the same time, the demand for gas is also volatile and uncertain in some places; therefore, the uncertainty associated with IESs mainly comprises uncertainties pertaining to wind power, solar power, and gas load.

Probability Model of Wind Power Plants
The probability model of wind speed satisfies the two-parameter Weibull distribution. Wind speed probability v can hence be computed as where p k is the shape parameter of the Weibull distribution and p c is the scale parameter. The power generated by wind turbines is proportional to the wind speed, as follows: where P and  are the relevant scaling coefficients, () px is the power of the wind turbines, ci x is the cut-in wind speed, and co x is the cut-out wind speed.

Uncertainty Analysis and Uncertainty Modeling
With the rapid development of new energy sources in recent years, the capacity for new energy and power generation in IESs is increasing quickly, which causes considerable uncertainty in these systems. At the same time, the demand for gas is also volatile and uncertain in some places; therefore, the uncertainty associated with IESs mainly comprises uncertainties pertaining to wind power, solar power, and gas load.

Probability Model of Wind Power Plants
The probability model of wind speed satisfies the two-parameter Weibull distribution. Wind speed probability v can hence be computed as where k p is the shape parameter of the Weibull distribution and c p is the scale parameter. The power generated by wind turbines is proportional to the wind speed, as follows: where P and δ are the relevant scaling coefficients, p(x) is the power of the wind turbines, x ci is the cut-in wind speed, and x co is the cut-out wind speed. The probability model of light intensity can be approximated as a beta distribution for a given period of time. Its probability density function is The output characteristic of a photovoltaic power plant is as follows: In the above formula, G stc and p stc j are the illumination intensity and relevant scaling coefficient, respectively, and G stc (t) and G stcmax are respectively the actual light intensity and maximum light intensity over a certain period. Parameters α and β are the shape parameters of the beta distribution.

Probability Model of Gas Load
The uncertainty of the gas network load can be fit using the following normal distribution function: where µ is the mean of the load, δ is its standard deviation, and x is the gas consumed under the gas load.

Nataf Transformation Theory
The Nataf transformation can transform random variables from any distribution to a normal distribution, which is convenient for correlation control, and has high precision. This method was proposed by Liu in [30], who provided the correlation coefficient transformation table for the Nataf transformation, which proved extremely beneficial for researchers in this area. For ease of understanding, the main steps of the Nataf transformation are briefly described below.
First, we set the input random variable vector Here, Y = y 1 , y 2 . . . y n is the standard normal random vector, Φ(·) is the standard normal cumulative distribution function, and Φ −1 (·) is the inverse cumulative distribution function.
The joint probability density function of random vectors can be derived using the derivative rule of implicit functions as follows: Appl. Sci. 2019, 9, 3291 7 of 23 Assuming the correlation coefficient matrix of the input random vector X is ρ, the expressions for each component of the correlation coefficient matrix can be derived as follows: where ρ 0ij is a component of the correlation coefficient of standard normal random vector Y. When we obtain all the values of ρ 0ij , the Cholesky decomposition of ρ 0 is ρ 0 = L 0 L T 0 . The relevant normal random vectors Y can be transformed into independent normal random vectors U and U = L −1 0 Y. Moreover, the following inverse transformation can be performed:

Combining of Nataf Transformation with PEM
The m-PEM are defined according to the number of estimated points. Generally speaking, the bigger of m, the more accurate results will be got, but the computational burden will become larger accordingly. The PEM is a type of Gauss-Hermite integration. For integrals with a normal distribution, the nodes and weights of the Gauss-Hermite integration can be easily obtained from tables.
In this method, the equal probability rule is used to transform the space of the normal distribution into an arbitrary distribution Assuming that the response function of pointin X in any systemis h is h = h(X 1 , X 2 · · · , X n ), we discuss two cases, considering the number of variables X in h.

Univariate response function
Using Formulas (7)-(12), the following relations can be obtained: where X = (x 1 , x 2 · · · , x k ), and k is the number of a set of samples. Let h = h(X) = h(x 1 , x 2 · · · , x k ); then the probability statistic of h is Substituting Equation (13) into Equations (14)-(15) yields the following formula: Because y is a variable of normal distribution, φ(y) is the probability density function of normal distribution, and the expansion of Equations (16) and (17) using the Gauss-Hermite integration is Equation (19) shows that the cost of computing a univariate response function for an m-PEM is proportional to m, where m is the number of estimation nodes. In addition, y i is the point estimated by m-PEM and p i is the weight coefficient of y i . Hence, m-PEM is expressed as follows: where ξ i is a position parameter of y i , which is the point estimated by m-PEM. If y satisfies the standard normal distribution, µ y = 0. Parameters ξ i and p i can be easily obtained from Table 1, which lists their values for a Gauss-Hermite integration.

Multivariable Response Function
Formulas (13)- (20) are used in a single-variable point estimation algorithm. In practice, the input of the response function usually has multiple variables. Assuming there are n variables in h with G = [G 1 (X 1 ), G 2 (X 2 ), · · · G n (X n )], h = h(X 1 , X 2 · · · , X n ), then the following equations can be used: where X i = (x i1 , x i2 · · · , x ik ) , i = 1 · · · n , and k is the number samples in the set.
Then the variance and mean of h can be obtained from Substituting Equation (21) into Equations (22) and (23) yields the following formulas: Appl. Sci. 2019, 9, 3291 9 of 23 Because y is a variable with a normal distribution, φ(y) is the probability density function of normal distribution. Hence, the expansion of Formula (24)-(25) by the Gauss-Hermite integration is Equations (26) and (27) show that the computation cost of an n-response function for m-PEM is proportional to m × n, where, m is the number of estimation nodes. In addition, y j denotes the points estimated by m-PEM and p j is the weight coefficient of y j . Hence, m-PEM is expressed as follows: where ξ i is the position parameter of y i . If y satisfies the standard normal distribution, then µ y = 0. Again, ξ i and p i can be easily obtained from Table 1.

MPEM with Less Computational Cost
Because the cost of calculating PEF is proportional to m × n using above method, this paper proposes the following approximation method that can be used to calculate m-PEM conveniently and the reduce number of PEF calculations to (m − 1) × n + 1: where, For 5-PEM and 7-PEM, when k = 1, we have

Algorithm for PEM Combined with the Nataf Transformation
Step 1: Obtain i initial sets of independent normal distribution samples Y 0 .
Step 2: Initialize correlation matrix ρ in the original probability space and determine relevance matrix ρ 0 in the normal space using Equation (10).
Step 3: Obtain the Cholesky decomposition of ρ 0 :ρ 0 = L 0 L T 0 and obtain new normal distribution samples with relevance Y:Y = L 0 Y 0 .
Step 4: Select the estimated numbers of each variable for the PEM corresponding to n-PEM.
Step 5: According to the Gauss-Hermite integration values in Table 1, obtain the estimation points and weight coefficients; obtain k points in the standard normal space Y.
Step 6: Determine all values of x i in the original probability space X using the inverse Nataf transformation.
Step 7: Obtain the response function of point x in any system h.
Step 8: Determine if h is a univariate or multivariable response function to obtain the estimated mean and density functions using Equations (18)- (20) or Equations (29)-(31).

PEM Considering Boundary Constraints
In this study, the boundary constraints indicate if the estimated wind speed is lower than the cut-in wind speed or exceed the cut-out wind speed. When estimated wind speed is outside of these boundary constraints, the estimated points cannot represent the whole sample well. Especially for 2-PEM and 3-PEM, the estimated points obtained from Table 1 will cause a large error. Therefore, the accuracy of PEM in normal space will degrade. In [32], Hong proposed 2-PEM and 3-PEM for an arbitrary probability space; this method, unlike the PEM proposed by Zhao and Ono [33], which should be transformed to a normal distribution space, makes it easier to constrain the estimation points in the sample space. Using Hong's method, this paper proposes power transformation and equality constraint transformation methods for 2-PEM and 3-PEM when the estimated points are outside of the range.

Hong's 2-PEM and 3-PEM
Hong's 2-PEM and 3-PEM are classic methods that are used widely [34,35]. To facilitate the understanding of the proposed method, these two methods are described briefly below.
Assuming x is the sample in an arbitrary probability space, µ x is the mean, and σ x is the standard deviation of x, λ x,i is the ratio of the ith center distance to the ith power of the standard deviation.
In addition, f (x) is the probability of x. The value of λ x,i is found using Reference [32] has shown the entire derivation process, and the results of the 2-PEM and 3-PEM are as follows: Here, ξ i and p i represent the position parameter and weight coefficients and are written as follows: Here, n is the number of random variables. The 2-PEM and 3-PEM equations for PEF calculations are hence as follows: The variances of the estimated probability density function can be derived as follows: If the estimated points are outside of the range when using Zhao's or Hong's PEM, this problem can be solved by processing the sample space before using PEM in the sample space, which is described in the next section.

Power Transformation Method
The power transformation method can transform an asymmetric probability distribution into a symmetric or approximately symmetric distribution. A symmetrical distribution function can reduce the possibility of estimated points falling outside the boundary constraints. The definition of a power transformation is The skewness of y is The wind speed samples x are discrete, so the values of y are discrete: It can be inferred that when λ y,3 is equal to 0, the power-transformed function is a symmetric function. Then, the estimated points can be directly obtained as follows: Substituting y i into function f with the range constraints to obtain f (x i ) yields However, if x i remains outside of the range, the power transformation method will be invalid.
Usually, it is difficult to find α such that λ y,3 is 0, so the following iterative method is used to find an α such thatis within the limits: Step 1. Calculate the variable's original skewness λ y,3 . Let a j = ka j−1 , j ≥ 0, and a 0 = 1; if λ y,3 is greater than 0, set the coefficient value to k = 1+ λ x,3 ) −1/9 . If the skewness coefficient is less than 0, set it to k= (1 − λ x,3 ) 1/9 .
Step 2. Set the maximum number of steps j to 100; calculate each step a j to obtain y i,j (i = 1, 2) and λ j,x,3 . Then, calculate y i,j and determine whether x i,j remains continuously within the limits. If it meets this criterion, the iteration stops.
Step 3. If α cannot be determined using the above steps, find the smallest λ l,x,3 . Let the number of iterations used to obtain λ l,x,3 be l; if λ l,x,3 is smaller than 0.01, a suitable value for α cannot be found such that is x i is within the limits.
Step 4. If λ l,x,3 is greater than 0.01, set n = 1 and change Calculate y i,j and then determine whether x i,j is within the allowed range from iteration l − 1 to iteration l+3 n . If x i, j is outside the limits, set n = n + 1 and repeat the above steps until x i,j within the range; then, α is determined. However, if λ l,x,3 is less than 0.01, then α cannot be found.

Equal Constraint Transformation
The equal constraint transformation method either transforms the original probability space into a space that is equivalent to the constrained probability space or it transforms the original space into the constrained probability space directly. The advantage of this method is that the estimation points are sure to be within the constraint, so it is a simple and more general method.
Suppose F(x) is a constraint function of x, which has the following format: Transform x such that it is an equal constraint to F(x) as follows: Then, 2-PEM and 3-PEM are used to calculate x to obtain the estimated points x 1 . . . , x m ; next, an inverse transformation is performed to obtain the estimated points x 1 · · · , x m in the original space as follows: It is clear that when k = 1, the transformed constraint space is the original constraint space.

Case Study
In this study, we simulated and analyzed the IES interconnected with an IEEE 30-bus power system and a 13-bus natural gas system. The overall diagram of the system is shown in Figure 2. There are three wind power plants and two solar power plants in the IES. The natural gas load at number 3 is a fluctuating load that satisfies the normal distribution. Specific data values for natural gas systems can be found in [1]. A gas turbine is equivalent to a power source in the power system and a load in the natural gas network. The conversion between input gas flow and output electric power is as follows: H g,j = a g,i + b g,i P G,i + c g,i P 2 where H g,j is the input heat value of gas turbine node j, P G,i is the output power of the gas turbine, F m,i gas is the equivalent gas load of gas node m in natural gas, GHV is a fixed high calorific value (GHV= 0.2275), and a g,i , b g,i , and c g,i are determined by the heat consumption curve of the gas turbine. In this case study, we simplify the heat consumption curve of the gas turbine and regard it as a linear relationship, that is, a g,i = 0, b g,i = 7, and c g,i = 0.
Hence, the relationship between the gas consumed by steam turbines and the electricity supplied to the power grid is as follows: The compressor is equivalent to load for the power system, and its power consumption is as follows: where P i c is the power consumed by compressors, p j and p m are the pressures of nodes j and m, respectively, Q c is the gas flow rate from node m to node j, B is a coefficient (B = 306.2746), k is the compression ratio power, and k = 2. Hence,

Case 1: Estimation Points within the Limit Range
In this case, the mean of wind speed is 8.5 m/s, and each wind power plant follows the distribution given in Equation (9) . Moreover, each photovoltaic power plant follows the distribution given in Equation

Case 1: Estimation Points within the Limit Range
In this case, the mean of wind speed is 8.5 m/s, and each wind power plant follows the distribution given in Equation (9), where k p = 2, c p = 10, x ci = 3.5 m/s , x co = 18m/s, P = 0.315MW. Moreover, each photovoltaic power plant follows the distribution given in Equation (10), where x is the standard light intensity value, x = G stc (t)/G stcmax , α = 3, β = 2, Γ(α+β) Γ(α)+Γ(β) = 20, and p stc = 0.5MW. The mean of the gas load 3 is µ = 204m 3 /h, its standard deviation is δ = 2, We then have Figure 3 shows the probability density of the sampling points and probability distribution of the power plants and photovoltaic power plants after the Nataf transformation. The appearance of this figure noisy is because the pseudo random variables produced by the computer do not strictly satisfy the probability density function of samples. Figure 3c,d show that the results agree well the theoretical sampling points, indicating that the Nataf transformation has a very good precision.   Figure 4 shows the correlation between each wind power plant and photovoltaic power plant after the Nataf transformation. The distribution of the sampling points is all concentrated in an ellipse. A flatter ellipse indicates a stronger correlation. It can be seen that the Nataf transformation controls the correlation well.    Figure 4 shows the correlation between each wind power plant and photovoltaic power plant after the Nataf transformation. The distribution of the sampling points is all concentrated in an ellipse. A flatter ellipse indicates a stronger correlation. It can be seen that the Nataf transformation controls the correlation well. Next, these points are used to calculate the PEF. Figure 5 shows the results of m-PEM combined with the Nataf transformation. To show specific errors, Tables 2 and 3 show the errors in the variance and mean. The calculation results show that this method has a very good accuracy when calculating the mean value, especially when calculating the mean value of the grid voltage. When calculating variance, it performs worse, which may be related to the number of Monte Carlo iterations. Because the PEF solution to an IES is very time-consuming, the number of Monte Carlo iterations must be as low as possible (3000). Next, these points are used to calculate the PEF. Figure 5 shows the results of m-PEM combined with the Nataf transformation. To show specific errors, Tables 2 and 3 show the errors in the variance and mean. The calculation results show that this method has a very good accuracy when calculating the mean value, especially when calculating the mean value of the grid voltage. When calculating variance, it performs worse, which may be related to the number of Monte Carlo iterations. Because the PEF solution to an IES is very time-consuming, the number of Monte Carlo iterations must be as low as possible (3000). For 2-PEM, 3-PEM, 5-PEM, and 7-PEM, the computational quantities are 12, 18, 25, and 37, respectively. and mean. The calculation results show that this method has a very good accuracy when calculating the mean value, especially when calculating the mean value of the grid voltage. When calculating variance, it performs worse, which may be related to the number of Monte Carlo iterations. Because the PEF solution to an IES is very time-consuming, the number of Monte Carlo iterations must be as low as possible (3000). For 2-PEM, 3-PEM, 5-PEM, and 7-PEM, the computational quantities are 12, 18, 25, and 37, respectively. To compare the results of PEF with and without the Nataf transformation, we used a common pseudo-random sample method for wind and photovoltaic farms. Figure 6 shows the correlation of common pseudo-random sampling points. It can be seen that the correlation of samples is very small: the correlation graph of the sample is approximately distributed in a circle in the plane. If the random error of sampling is ignored, the samples are independent of each other. 10 15 20 25 30 ind speed 2(m/s) 10 15 20 25 Wind speed 3(m/s)   Table 3. Average relative error of the power network. To compare the results of PEF with and without the Nataf transformation, we used a common pseudo-random sample method for wind and photovoltaic farms. Figure 6 shows the correlation of common pseudo-random sampling points. It can be seen that the correlation of samples is very small: the correlation graph of the sample is approximately distributed in a circle in the plane. If the random error of sampling is ignored, the samples are independent of each other. To compare the results of PEF with and without the Nataf transformation, we used a common pseudo-random sample method for wind and photovoltaic farms. Figure 6 shows the correlation of common pseudo-random sampling points. It can be seen that the correlation of samples is very small: the correlation graph of the sample is approximately distributed in a circle in the plane. If the random error of sampling is ignored, the samples are independent of each other. Next, the results of the two sampling methods were compared using MCS. Figure 7 shows the comparison results of MCS with and without the Nataf transformation. The effect of sampling without the Nataf transformation is even greater for the power system, because the correlation of the variables in this example is in the power network. The cumulative probability distributions of Figure  7d,e strongly support this point. Next, the results of the two sampling methods were compared using MCS. Figure 7 shows the comparison results of MCS with and without the Nataf transformation. The effect of sampling without the Nataf transformation is even greater for the power system, because the correlation of the variables in this example is in the power network. The cumulative probability distributions of Figure 7d,e strongly support this point. Next, the results of the two sampling methods were compared using MCS. Figure 7 shows the comparison results of MCS with and without the Nataf transformation. The effect of sampling without the Nataf transformation is even greater for the power system, because the correlation of the variables in this example is in the power network. The cumulative probability distributions of Figure  7d,e strongly support this point.

Case 2: Estimation Points Outside the Limits
In this case study, there is no wind in the environment, so the wind speed is lower than normal speed (mean = 7.14 m/s). Hence, the wind speed meets following conditions: Just as in Case 1, the relevance between each power plant and photovoltaic power plant is added using the Nataf transformation. Figure 8 shows the results of points that are outside of the boundary constraints, and Tables 4 and 5 show the errors in the variance and mean. It can be seen that the outside-of-range point has great impact on 2-PEM and 3-PEM. In these cases, the average error of the variance is more than 20%. The impacts on 5-PEM and 7-PEM are smaller.

Case 2: Estimation Points Outside the Limits
In this case study, there is no wind in the environment, so the wind speed is lower than normal speed (mean = 7.14 m/s). Hence, the wind speed meets following conditions: where k p = 2, c p = 65, x ci = 3.5m/s, x co = 18m/s, P = 0.315MW, and the photovoltaic data and gas load 3 are the same as those in Case 1. Just as in Case 1, the relevance between each power plant and photovoltaic power plant is added using the Nataf transformation. Figure 8 shows the results of points that are outside of the boundary constraints, and Tables 4 and 5 show the errors in the variance and mean. It can be seen that the outside-of-range point has great impact on 2-PEM and 3-PEM. In these cases, the average error of the variance is more than 20%. The impacts on 5-PEM and 7-PEM are smaller.     To limit points that are outside of the boundary constraints, the power transformation and equality constraints transformation methods are used to recalculate Case 2. Tables 6-8 show the results of limiting the wind speed estimation points. It can be seen when the power transformation method works, the two methods have similar accuracy, but the equality constraint transformation method has wider applicability. In general, both methods effectively improve the accuracy of calculation.  To limit points that are outside of the boundary constraints, the power transformation and equality constraints transformation methods are used to recalculate Case 2. Tables 6-8 show the results of limiting the wind speed estimation points. It can be seen when the power transformation method works, the two methods have similar accuracy, but the equality constraint transformation method has wider applicability. In general, both methods effectively improve the accuracy of calculation.

Discussion
Liu presented a correlation coefficient transformation table for the Nataf transformation [30], which was convenient for Zhao and Ono's MPEM. This study uses the work of these researchers to solve PEF. On the one hand, for 2-PEM and 3-PEM based on Zhao and Ono's MPEM, estimated points outside of the range will cause a relatively large error. This study utilized the space conversion concept to overcome this limitation. The results show that spatial transformation is an important and effective method for solving PEF. One study [26] provided a good example of this transformation. Its aim was to link the probability space of an elliptic space with the probability space of the sample; the PPF was then solved. On the other hand, for 5-PEM and 7-PEM based on Zhao and Ono's MPEM, direct calculations will lead to an exponential growth in the computation cost; therefore, we developed a modified MPEM.
Future research on probabilistic methods will focus more on accuracy and less on the cost of computation. As for the work discussed in this paper, the question as to whether 5-PEM and 7-PEM can allow estimation points outside the limit range without affecting accuracy is worth further investigation. In addition, combining PEF problems with economic dispatch, market operations, and optimization will also be directions for future research.

Conclusions
A PEM based on the Nataf transformation for computing PEF in an IES was proposed in this paper. This method is appropriate for large-scale electricity, gas, and thermal power interconnected systems because the energy flow calculations for such systems are complicated and time consuming; therefore, it is essential to reduce the cost of computation. From the results of the case studies, the following conclusions can be drawn:

1.
The proposed modified MPEM is effective in situations where the numbers of estimation points and variables are relatively large. The combination of the Nataf transformation with MPEM can obtain good accuracy. The proposed method makes the computational cost of m-PEM equal to that of (m − 1)-PEM, thus improving the computational efficiency.

2.
The power transformation and equality constraint transformation methods are used to reduce the large errors when the estimation points are outside of the boundary constraints. Simulation results show the effectiveness of the proposed methods.
Note that these methods are not only suitable for IESs but also for other scientific applications, such as the reliability analysis of structures and geotechnical reliability analysis.

Conflicts of Interest:
The authors declare no conflict of interest.  Reactive power injected from the gas generator into node i Q i Injection of natural gas at the ith node f im Injection to the downstream node f in Injection to the upstream node F j Gas consumption of the compressor G ij Correlation coefficient p Active unbalance q

Appendix A
Reactive unbalance f NG Gas flow unbalance P and δ Relevant scaling coefficients p(x) Power of the wind turbines x ci Cut-in wind speed x co Cut-out wind speed Standard normal cumulative distribution function Φ −1 (.) Inverse cumulative distribution function ρ 0ij Component of the correlation coefficient of the standard normal distribution vector G = [G 1 (X 1 ), G 2 (X 2 ), · · · , G n (X n )] Arbitrary probability distribution space h = h(X 1 , X 2 · · · , X n ) Response function φ(y) Probability density function of normal distribution m Point number of PEM y i Point estimated by m-PEM p i Weight coefficient of y i k Number of samples in a set ξ i Position parameter of y i p j Weight coefficient of y j λ y, 3 Skewness of y l Number of iterations to obtain λ y,3 H g, j Input heat value of gas turbine node j P G,i Output power of gas turbine F m,i gas Equivalent gas load of gas node m in natural gas GHV Fixed high calorific value a g,i , b g,i , c g,i Coefficient determined by the heat consumption curve of a gas turbine P i c Power consumed by compressors p j Pressure of node j p m Pressure of node m Q c Gas flow rate from node m to node j B Coefficient of compressors k Compression ratio power