1. 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 challenges in terms of energy uncertainties arise. The uncertainties come from the intermittency of RESs; additionally, the volatility of the loads in different energy structures aggravate the energy uncertainties in IESs [
9], which creates significant challenges for the operation and maintenance of IESs. We can determine the probabilistic distribution of each node state in the IESs by calculating the probabilistic energy flow, which provides guidance to address possible risks. Therefore, the probabilistic energy flow (PEF) in IESs must be studied.
The basis of PEF is the calculation of steady-state energy flow (SSEF), which has been studied extensively. Because of the similarity between the SSEF model of the power system and the gas network, Newton’s method can be used to solve the SSEF of the gas network. When the initial point of Newton’s method is near the equilibrium point, the method has good and fast convergence. Moreover, the IES mostly works near the equilibrium point. Therefore, Newton’s method is a commonly used method for solving the energy flow problem. Moreover, the solutions obtained using this method for multiple energy systems can be divided into two forms according to the Jacobian matrices used in the solution, i.e., separate Newton iterations for each system [
10] or a unified Newton iteration for all systems [
11]. In [
10], each energy system (e.g., power system or natural gas system) is represented by its own Jacobian matrix and independent iterations are performed using Newton’s method. The unbalanced quantity is iterated using the energy coupling unit until the whole energy system achieves convergence. This method is relatively simple, but the computation speed is somewhat slow. In [
11], the power and natural gas networks are represented by a single Jacobian matrix for the unified Newton iterations. Then, the integrated energy flow of the entire network can be obtained. This method is more complicated in cases where there is a compressor in the gas grid, but the iterations are faster.
Presently, research on the impact of uncertainties in IESs is limited. However, the calculations performed for these systems are similar to the probabilistic power flow (PPF) calculations in traditional power systems. They include the Monte Carlo method [
12], analytical method [
13], stochastic spectral analysis method [
14], and approximation method [
15].
Generally, a Monte Carlo simulation (MCS) requires several samples to obtain good precision according to the law of large numbers; this simulation requires considerable time to compute. However, it is possible to achieve the most accurate results with an MCS; therefore, it is used as reference to measure the precision of other methods. The studies on MCSs mainly focus on improving sampling efficiency to reduce the number of sampling representatives, for instance stratified sampling [
16,
17] and importance sampling [
18]. Others focus on improving the calculation efficiency using methods such as parallel computing [
19]. The sampling sequence can also be changed to reduce the number of samples, such as in the quasi-Monte Carlo method [
20]. This method obtains random numbers using the Sobol sequence, which generates more uniformity than random numbers generated by the pseudo random number technique.
The analytical method is used to solve the PPF problem with a linearized alternating current model and convolution. The number of calculations in this method increases exponentially when the system is relatively large and several components need to be considered. Therefore, this method is not as widely used as the other three methods.
Stochastic spectral analysis includes the Karhunen–Loeve expansion method (KLEM) [
21], polynomial chaos expansion [
22], and surrogate modelling [
23]. KLEM expands the system response into a deterministic part and a stochastic part. The stochastic part consists of a series of unrelated random variables and deterministic coefficients (eigenvalues and eigenfunctions). The disadvantage of KLEM is that it is not applicable to the covariance function of unknown stochastic processes. Polynomial chaos expansion expands the structural probability using the orthogonal polynomial of the basic variable and uses the chaotic polynomial as the basis of the expansion. Therefore, it can be applied to both Gaussian and non-Gaussian outputs. However, when the number of uncertain parameters increases, the number of polynomial coefficients also increase, which greatly increases the cost of computation. Hence, stochastic spectral analysis method is a new method developed recently and remains an on-going research.
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.
2. PEF Model for IESs
2.1. 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, and are respectively the active and reactive power that the source injects into node , and are respectively the active and reactive power demanded by the node load, and are respectively the active and reactive power consumed by the compressor, and and are respectively the active and reactive power that the gas generator injects into node . In addition, represents the number of nodes directly connected to the node, is the injection of natural gas at the node, and and represent the injections at the downstream and upstream nodes, respectively. The gas consumption of the compressor is denoted by : when the compressor takes air from node , the correlation coefficient is 1; otherwise, is 0. If the energy required by the compressor comes from the power system, then is also 0. Finally, , and 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.
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).
2.2. 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.
2.2.1. Probability Model of Wind Power Plants
The probability model of wind speed satisfies the two-parameter Weibull distribution. Wind speed probability
can hence be computed as
where
is the shape parameter of the Weibull distribution and
is the scale parameter. The power generated by wind turbines is proportional to the wind speed, as follows:
where
and
are the relevant scaling coefficients,
is the power of the wind turbines,
is the cut-in wind speed, and
is the cut-out wind speed.
2.2.2. Probability Model of Photovoltaic Power Plants
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, and are the illumination intensity and relevant scaling coefficient, respectively, and and are respectively the actual light intensity and maximum light intensity over a certain period. Parameters and are the shape parameters of the beta distribution.
2.2.3. 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
is the gas consumed under the gas load.
4. 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.
4.1. 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
is the sample in an arbitrary probability space,
is the mean, and
is the standard deviation of
,
is the ratio of the
center distance to the
power of the standard deviation. In addition,
is the probability of
. The value of
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,
and
represent the position parameter and weight coefficients and are written as follows:
Here,
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.
4.2. 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 wind speed samples
are discrete, so the values of
are discrete:
It can be inferred that when
is equal to 0, the power-transformed function is a symmetric function. Then, the estimated points can be directly obtained as follows:
Substituting
into function
with the range constraints to obtain
yields
However, if remains outside of the range, the power transformation method will be invalid.
Usually, it is difficult to find such that 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 . Let , , and ; if is greater than 0, set the coefficient value to . If the skewness coefficient is less than 0, set it to .
Step 2. Set the maximum number of steps to 100; calculate each step to obtain and . Then, calculate and determine whether 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 . Let the number of iterations used to obtain be ; if is smaller than 0.01, a suitable value for cannot be found such that is is within the limits.
Step 4. If is greater than 0.01, set and change to in iteration . Calculate and then determine whether is within the allowed range from iteration to iteration . If is outside the limits, set and repeat the above steps until within the range; then, is determined. However, if is less than 0.01, then cannot be found.
4.3. 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
is a constraint function of
, which has the following format:
Transform
such that it is an equal constraint to
as follows:
Then, 2-PEM and 3-PEM are used to calculate
to obtain the estimated points
; next, an inverse transformation is performed to obtain the estimated points
in the original space as follows:
It is clear that when k = 1, the transformed constraint space is the original constraint space.
5. 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:
where
is the input heat value of gas turbine node
,
is the output power of the gas turbine,
is the equivalent gas load of gas node
in natural gas,
is a fixed high calorific value (
), and
,
, and
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,
,
, and
.
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
is the power consumed by compressors,
and
are the pressures of nodes
and
, respectively,
is the gas flow rate from node
to node
,
is a coefficient (
),
is the compression ratio power, and
. Hence,
5.1. 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
,
,
,
,
. Moreover, each photovoltaic power plant follows the distribution given in Equation (10), where
is the standard light intensity value,
,
,
,
, and
. The mean of the gas load 3 is
, its standard deviation is
, 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.
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,
Table 2 and
Table 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.
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.
5.2. 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
,
,
,
,
, 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
Table 4 and
Table 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.
For 2-PEM, 3-PEM, 5-PEM, and 7-PEM, the computation quantities are 12, 18, 25, and 37, respectively.
To limit points that are outside of the boundary constraints, the power transformation and equality constraints transformation methods are used to recalculate Case 2.
Table 6,
Table 7 and
Table 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.