Interval Energy Flow Analysis in Integrated Electrical and Natural-Gas Systems Considering Uncertainties

: As integrated electrical and natural-gas systems (IENGS) are popularized, the uncertainties brought by variation of electrical load, power generation, and gas load should not be ignored. The aim of this paper is to analyze the impact of those uncertain variables on the steady-state operation of the whole systems. In this paper, an interval energy ﬂow model considering uncertainties was built based on the steady-state energy ﬂow. Then, the Krawczyk–Moore interval iterative method was used to solve the proposed model. To obtain precise results of the interval model, interval addition and subtraction operations were performed by a ﬃ ne mathematics. The case study demonstrated the e ﬀ ectiveness of the proposed approach compared with Monte Carlo simulation. Impacts of uncertainties brought by the variation of electrical load, power generation, and gas load were analyzed, and the convergence of energy ﬂow under di ﬀ erent uncertainty levels of electrical load was studied. The results led to the conclusion that each kind of uncertainties would have an impact on the whole system. The proposed method could provide good insights into the operating of IENGS with those uncertainties.


Introduction
With the increasing demand for energy and the growing environmental issues, it is urgent to utilize energy more efficiently. One way to overcome this problem is to make full use of the integrated electrical and natural-gas systems. The integrated electrical and natural-gas systems (IENGS) can realize this goal through the strengthening of the coordination between electrical systems and natural-gas systems [1]. For example, coupling units such as gas-fired generators and P2G (power to gas) in the electrical and natural-gas systems achieves effective allocation of energy by real-time scheduling. IENGS has the advantage of high energy efficiency, less environmental pollution, and sustainable energy supply, which is a trend of energy development in the future.
As one of the most important analytical tools for IENGS, energy flow analysis lays the foundation of the optimal operation [2][3][4][5][6], planning [7][8][9][10], energy market [11], and optimal scheduling [12] of the IENGS. There are some studies focused on energy flow analysis with different energy coupling links. A decomposing energy flow strategy for a multienergy carrier system has been proposed in Ref. [13] and the interdependent equipment, such as combined heating and power (CHP), boiler, and electrical pump, is taken into account. Ref. [14] considers the effect of temperature in the natural-gas system and a distributed slack node technique in the electricity network and proposes an integrated formulation for the steady-state analysis of electricity and natural-gas coupled systems based on the Newton-Raphson formulation. The interactions between power and gas systems by microturbines have been modeled and dynamic impact of interdependency has been studied in Ref. [15]. A unified steady-state energy flow formulation is developed with bidirectional energy conversion: power to

•
Interval energy flow model of IENGS is built to reflect the various uncertainties brought by electrical load, power generation, and gas load, and the Krawczyk-Moore interval iterative method is used to solve interval energy flow equations and verify its correctness with the comparison to Monte Carlo simulation.

•
The influence of each kind of uncertainty brought by electrical load, power generation, and gas load has been investigated, and especially the impact of uncertainty in power system on gas system and the impact of uncertainty in gas system on power system are analyzed.

•
The convergence of energy flow is analyzed when the uncertainty level of electrical load is increased to verify the approach's effectiveness.
The remainder of this paper is organized as follows. The interval energy flow model considering uncertainties is built based on the steady-state energy flow in Section 2. In Section 3, the Krawczyk-Moore interval iterative method is used to calculate energy flow. Case study and simulation results are illustrated in Section 4, and Section 5 draws the conclusion of the paper.

Energy Flow Model in IENGS Considering Uncertainties
An IENGS consists of the power system, the natural-gas system, and coupling units. Electrical power is transmitted from the power supply side to the load side by the power network, and the natural-gas system is applied to transmit gas from gas source to gas load through pipe lines and compressors. Interconversion of electricity and gas is achieved by equipment such as CHP units, electric boilers, gas-fired generators, and P2G device. This paper takes the gas-fired generator as a coupling unit for the interconnected linkage. The IENGS considered in this paper is shown in Figure 1.
Energies 2019, 12, x FOR PEER REVIEW  3 of 19 electric boilers, gas-fired generators, and P2G device. This paper takes the gas-fired generator as a coupling unit for the interconnected linkage. The IENGS considered in this paper is shown in Figure 1.

Steady-State Energy Flow Model without Consideration of Uncertainty
To build the steady-state energy flow model of IENGS, the power flow model is given firstly, and then the gas flow model is built with reference to the modeling process of power flow since there are many similarities between them. Lastly, the coupling unit model is given.

Power Flow Model
For each bus i in the power system, the power flow calculation model under polar coordinates is expressed as: Where Pi and Qi are active power and reactive power injected into bus i, respectively. Ui is the voltage amplitude of bus i, and δij is the voltage phase angle difference between bus i and j. Gij and Bij are conductance and susceptance between bus i and j in bus admittance matrix, respectively. nE is the total bus number in the whole power system. Buses in the power system can be classified into PQ buses, PV buses, and slack buses, and each bus satisfies Kirchhoff's first law in the power system, that is, nodal power balance equations of bus i can be expressed as: ,, where Ps,i and Qs,i are active power and reactive power generated by bus i, respectively. Pl,i and Ql,i are active power and reactive power consumed by bus i, respectively.

Gas Flow Model
Natural-gas systems have many similarities with power systems. A structural comparison between gas systems and power systems is illustrated in Table 1.

Steady-State Energy Flow Model without Consideration of Uncertainty
To build the steady-state energy flow model of IENGS, the power flow model is given firstly, and then the gas flow model is built with reference to the modeling process of power flow since there are many similarities between them. Lastly, the coupling unit model is given.

Power Flow Model
For each bus i in the power system, the power flow calculation model under polar coordinates is expressed as: where P i and Q i are active power and reactive power injected into bus i, respectively. U i is the voltage amplitude of bus i, and δ ij is the voltage phase angle difference between bus i and j. G ij and B ij are conductance and susceptance between bus i and j in bus admittance matrix, respectively. n E is the total bus number in the whole power system. Buses in the power system can be classified into PQ buses, PV buses, and slack buses, and each bus satisfies Kirchhoff's first law in the power system, that is, nodal power balance equations of bus i can be expressed as: where P s,i and Q s,i are active power and reactive power generated by bus i, respectively. P l,i and Q l,i are active power and reactive power consumed by bus i, respectively.

Gas Flow Model
Natural-gas systems have many similarities with power systems. A structural comparison between gas systems and power systems is illustrated in Table 1. There exist lots of pipelines to realize long-distance transmission of natural gas similar to the transmission line in power systems. An accurate model for calculating the flow of natural-gas pipelines is given in [14]. Ignoring the influence of temperature and the environment on the gas flow of pipelines, the gas flow model of pipelines in this paper is given as [31]: where f mn is the gas flow passing through nodes from m to n. v m and v n are the gas pressure at nodes m and n, respectively. h mn is a parameter which is related to gas pipelines and fluids and influenced by various factors such as the length and inner diameter of pipelines, the friction coefficient of natural gas flowing in the pipeline, and the natural-gas temperature. sgn(v m -v n ) is a symbolic function used to characterize the flow direction of natural gas in a pipeline, which is expressed as [27]: 2. Gas flow model of compressors Power transformers are typically placed at certain buses in order to reduce power losses caused by long-distance transmission in power systems. Natural-gas compressors are used to compensate the natural-gas pressure drop caused by the pipeline transmission process. Compressors increase the gas pressure level by an electric motor or a gas-fired turbine. The gas required for the gas-fired turbine compressor extracted from the natural-gas network is expressed as [14]: where H c,m is the energy consumption of the compressor. B c,mn and Z c,mn are the compressor parameters, which depend on the temperature, efficiency, and compression factor of the compressor. τ c,m is the gas required for the gas-fired turbine compressor extracted from the natural-gas network, and α, β, and γ are compressors gas consumption coefficients. In this paper, γ c is set 0 as is in Ref. [27], so the equation is converted to a linear function.

Node type
Analogous with power systems, nodes in the natural-gas network can be divided into two types: (1) known-pressure node. Gas supply of this node is sufficient and the pressure is a given reference value to compute all other unknown nodal pressures, which is similar to slack buses in power systems. (2) known-injection node. Gas consumption and injection are known while pressure is unknown for the node, which is analogous to the PQ buses in power systems.

Nodal gas flow balance equation
Kirchhoff's first law is satisfied by each node in the natural-gas network, namely, the sum of gas flow inlet is equal to the sum of gas flow output for each node, which can be expressed as: where M p and M c are the set of pipeline nodes and compressor nodes adjacent to the node m, respectively, and f l,m and f s,m are gas load and source at the node m, respectively.

Coupling Model
This paper considers the gas-fired generator as a coupling unit for the interconnected linkage. The gas flow required by the coupling is expressed as [14]: where P G,i is the power output of gas-fired generator at bus i in the power system, and F gas,m is the gas flow of the gas-fired generator of node m in the natural-gas system. In Equation (11), the fuel consumption of the gas-fired generator is related to net heating value (lower heating value) and the gross heating value. Here, we regard it as an ideal process and choose the gross heating value as in Ref. [14]. GHV is gross heating value, and α g , β g , and γ g are heat rate coefficients of gas-fired generators.

Steady-State Energy Flow Balance Equation
Considering the coupling unit, Equations (3) and (9) are rewritten as (12) and (13): So, the energy balance equation of the integrated power and natural-gas systems is corrected to Equations (4), (12) and (13).

Uncertainties in Interval Forms
There are a large number of uncertain factors in the operating IENGS. In this paper, three kinds of uncertainties are taken into consideration caused by variations of electrical load, power generation, and gas load. Interval mathematics provides an effective tool to describe and address uncertain variables. One of the advantages of applying interval mathematics is that only the upper and lower bounds of uncertainties are needed, especially when the probability distribution function of uncertainties is hard to obtain accurately.
The interval number and interval function are expressed by "[ ]" to distinguish a certain number and function. The number of interval [x] = [x, x] is the set of all x that satisfies x ≤ x ≤ x, which can be formulated by: where x ∈ R, x ∈ R, and x ≤ x. In addition, x and x are called the lower and upper limits of x, respectively, and R is the set of real numbers. The arithmetic operations are described in [32]. Hence, the three kinds of uncertainties (variations of electrical load, power generation, and gas load) in interval forms are expressed as: The variables in interval forms are formulated by:

Energy Flow Model Using Interval Mathematics
Now that the uncertainties and variables are obtained in interval forms, the energy flow equations considering uncertainties can be formulated by interval mathematics. The interval energy flow model of IENGS can be expressed as: Also, (17) can be simplified as nonlinear interval equations formulated by:

Krawczyk-Moore Iterative Method in Energy Flow Analysis
The interval energy flow model can be solved by the Krawczyk-Moore iterative method. The Krawczyk-Moore interval iterative method is based on Moore's research [33], which has been used in interval power flow analysis [34,35]. To apply this method, an interval operator should be built, which can be expressed as: , which is a certain vector, and Y is a nonsingular matrix, and The advantage of the Krawczyk-Moore iterative method is that the inverse of Jacobian matrix need not be calculated in each iteration by using the interval operator K(x 0 , [x]), since calculation for the inverse of Jacobian matrix needs mass computation and sometimes is quite difficult. The ith iteration process of the Krawczyk-Moore method is as follows: When the Krawczyk-Moore iterative method is used in the interval energy flow model above, [x] is expressed as (21) and (22) in the power system and gas system, respectively: Meanwhile, the corresponding F (x) is expressed as (23) and (24), respectively: where J E and J G are the Jacobian matrices of power flow and gas flow, respectively.

Improvement of the Krawczyk-Moore Iterative Process with Affine Arithmetic
When the Krawczyk-Moore iterative method is used to solve the interval energy flow model, the solution may be too conservative since there exits correlation between interval variables in [J E ] and [J G ] [36]. In this situation, affine mathematics can describe the correlation and reduce the conservativeness effectively.

Affine Arithmetic
Affine mathematics is an improved form of interval mathematics, which can describe the interconnection between uncertain variables, thus making the calculation results of interval energy flow more precise [37]. A quantityx is represented as an affine form (25): where x 0 and x i are the center value and partial deviation ofx, respectively, and ε i is the noisy symbol that lies in the interval [−1,1]. The interval form and the affine form of a quantity can be converted to each other. If [x] is converted tox, then x 0 = (x + x)/2 and x 1 = (x − x)/2 in two order form, and

Improvement in Iterative Process with Affine Arithmetic
Affine mathematics can optimize the process of calculation considering the interdependency between variables. Therefore, K(x) in Equation (19) is rewritten as: where x 1 is the partial deviation of [x], denoted as Considering there is correlation in [J E ] and [J G ], affine mathematics can describe the correlation in energy flow analysis quantitatively.

Improvement of [J E ]
Considering the correlation in [J E ] of the power system, there exits correlation of sin δ and cos δ in submatrices of H, N, J, and L. Take the submatrix H as the instance, the matrix elements in submatrix H can be expressed as [34]: where µ H,ij = arctan −B ij /G ij , which is a fixed value. Other submatrices of N, J, and L of Jacobian matrix in the power system can be calculated in the same way.

Improvement of [J G ]
Considering the correlation in [J G ] of the natural-gas system, v 2 m − v 2 n can be formulated by square difference formula and considering the correlation between v m and v n , [ f mn ] can be expressed as: In this way, by using square difference formula, only a calculation of interval multiplication is needed to be performed when [J G ] is calculated by [ f mn ], while addition and subtraction is performed by affine mathematics. Hence, the results will be more accurate by using (30) and (31):

Interval Energy Flow Analysis Process
The analysis process of interval energy flow in this paper is shown in Figure 2, which can be summarized in eight steps.  Step 1: Initialize the upper and lower bounds of variables. In this step, the maximum limit and minimum limit of variables in energy flow analysis are chosen as the initial values of the upper and lower bounds via (32) Step 2: Convert the interval variables to the affine form via (34) and (35).
Step 3: which is expressed as: Step 1: Initialize the upper and lower bounds of variables. In this step, the maximum limit and minimum limit of variables in energy flow analysis are chosen as the initial values of the upper and lower bounds via (32) and (33).
Step 2: Convert the interval variables to the affine form via (34) and (35).
Step 3: , which is expressed as: if i ≥ 1, Y i is expressed as: Step  (20).
Step 5: Calculate parameters of coupling unit. Firstly, calculate the power output of gas turbines, then calculate the gas flow consumed by gas-fired generators via (10) and (11).
Step 6: Calculate K G x i G in the gas system, which is similar to Step 4.
Step 7: Calculate the next iteration value of energy flow, which can be expressed as: Step 8: Determine convergence of energy flow results via (41) and (42) where ε err,E and ε err,G are the iterative errors in power flow and gas flow, respectively. If the results converge, then output the results, otherwise jump to Step 3 to continue the iteration.

Thirty-Nine-Bus Electrical System Coupled with 20-Node Gas System
The IENGS consisting of the IEEE 39-bus system and the Belgian natural-gas system shown in Figure 3 are presented to analyze the suitability of the proposed approach. The scales of the two subsystems match each other. The Belgian natural-gas system, whose detailed data can be found in Ref. [14], consists of 20 gas nodes and 24 gas pipelines, in which there are 7 gas source nodes, 8 gas load nodes, and 2 pressure nodes driven by external energy sources. There are 3 gas-fired generators at buses 30, 31, and 34 in the power system, which are supplied from nodes 4, 12, and 16 in the gas system, respectively. Bus 31 in the power system and node 1 in the gas system are selected as slack nodes.
Gas source and gas load information of the gas network is shown in Table 2, and detailed information of pipelines can be found in [38]. In addition, parameters of gas-fired generators are shown in Table 3. It should be noted that although a model with the possibility of using a quadratic one is set up, this model is transformed to a linear model because γ = 0. The value of GHV in (11) is 1015 BTU/SCF. In this paper, the British notation is used instead of the SI notation, and their conversion relations are 1 BTU/SCF = 37.25 × 10 3 J/m 3 , and 1 MSCF = 28.32 m 3 . Gas source and gas load information of the gas network is shown in Table 2, and detailed information of pipelines can be found in [38]. In addition, parameters of gas-fired generators are shown in Table 3. It should be noted that although a model with the possibility of using a quadratic one is set up, this model is transformed to a linear model because γ = 0. The value of GHV in (11) is 1015 BTU/SCF. In this paper, the British notation is used instead of the SI notation, and their conversion relations are 1 BTU/SCF = 37.25  10 3 J/m 3 , and 1 MSCF = 28.32 m 3 .   In this section, five scenarios are designed to analyze the energy flow in the IENGS. Scenario 1: without uncertainties. Scenario 2: with a ± 5% variation on electrical load. Scenario 3: with a ± 5% variation on the power generator of bus 36. Scenario 4: with a ± 5% variation on gas load. Scenario 5: considering uncertainties including Scenario 2, 3, and 4.

Comparison with Monte Carlo Simulation
To validate the effectiveness of the proposed approach, certain energy flow results of Scenario 1, interval energy flow solutions of Scenario 2 obtained by the proposed approach, and Monte Carlo (MC) stochastic simulation are calculated. MC stochastic simulation randomly samples load variation for 10 4 trials to get maximum/minimum voltage magnitude, voltage angle of the power system, and gas pressure of the gas system. It is assumed that solutions of MC simulation are regarded as the real interval results if sufficient numbers of samples are applied. Partial results are listed in Tables 4 and 5 for the sake of brevity and readability.  The results of Scenario 1 are lying in the solutions obtained by the proposed approach from the above comparison. Furthermore, the results calculated by the proposed approach is close to MC simulation results and former results include all the latter results, which can illustrate the completeness of the proposed approach. The detailed results are illustrated in Figure 4, which reflects the comparison between the proposed approach and MC simulation intuitively. Figure 4a shows the variation of voltage magnitude from bus 1 to bus 29 in the power system. It should be noted that upper bound and lower bound of voltage magnitudes from bus 30 to bus 39 obtained by the proposed approach and MC simulation have the same values, which are also equal to the results calculated by certain energy flow in Scenario 1. That is because bus 31 is the slack bus and the others are PV buses, whose voltage magnitude cannot change in the energy flow calculation process. Figure 4b shows the variation of voltage angle from bus 1 to bus 39 in the power system. Figure 4c shows the variation of gas pressure from node 1 to node 20 in the gas system. Voltage angle of bus 31 in Figure 4b and gas pressure of node 1 in Figure 4c will not change under three situations, as bus 31 in the power system and node 1 in the gas system are slack nodes. certain energy flow in Scenario 1. That is because bus 31 is the slack bus and the others are PV buses, whose voltage magnitude cannot change in the energy flow calculation process. Figure 4b shows the variation of voltage angle from bus 1 to bus 39 in the power system. Figure 4c shows the variation of gas pressure from node 1 to node 20 in the gas system. Voltage angle of bus 31 in Figure 4b and gas pressure of node 1 in Figure 4c will not change under three situations, as bus 31 in the power system and node 1 in the gas system are slack nodes. The comparison in Figure 4 demonstrates that the proposed approach can estimate the upper and lower bounds of energy flow in IENGS considering uncertainties. The results obtained by MC simulation are included by the proposed approach and the two kinds of results are very close. However, it takes less time to get the results by the proposed approach compared with MC simulation, hence the proposed approach has sound applicability for energy flow analysis considering uncertainty. The comparison in Figure 4 demonstrates that the proposed approach can estimate the upper and lower bounds of energy flow in IENGS considering uncertainties. The results obtained by MC simulation are included by the proposed approach and the two kinds of results are very close. However, it takes less time to get the results by the proposed approach compared with MC simulation, hence the proposed approach has sound applicability for energy flow analysis considering uncertainty.

Uncertain Energy Flow Analysis
To analyze the impact of uncertainties brought by variation of electrical load, power generation, and gas load, the results of Scenarios 2, 3, and 4 are calculated and those solutions are shown in Figure 5 together.

Analysis of electrical load variation
In Scenario 2, electrical load variation has an influence on the IENGS. That is to say, the uncertainty produced by the power system is delivered to the gas system, making gas pressure of its nodes fluctuate. Power generation in bus 31, which is supplied by the gas-fired generator from coupling node 12 in the gas system, varies with electrical load. The power generation of the gas-fired generator changes with its gas consumption. Then, the variation of the coupling node 12 causes the variation of gas load accordingly, which will make gas pressure of all gas nodes fluctuate. Hence the energy flow will change over electrical load variation.

Uncertain Energy Flow Analysis
To analyze the impact of uncertainties brought by variation of electrical load, power generation, and gas load, the results of Scenarios 2, 3, and 4 are calculated and those solutions are shown in Figure  5 together.

Analysis of electrical load variation
In Scenario 2, electrical load variation has an influence on the IENGS. That is to say, the uncertainty produced by the power system is delivered to the gas system, making gas pressure of its nodes fluctuate. Power generation in bus 31, which is supplied by the gas-fired generator from coupling node 12 in the gas system, varies with electrical load. The power generation of the gas-fired generator changes with its gas consumption. Then, the variation of the coupling node 12 causes the variation of gas load accordingly, which will make gas pressure of all gas nodes fluctuate. Hence the energy flow will change over electrical load variation.

Analysis of power generator variation
In Scenario 3, the uncertainty is reflected in the power system and gas system, similar to the situation in Scenario 2. When the output of power generator fluctuates by uncertainty, all the results of electrical nodes will fluctuate, as well as the variation of all gas nodes. That is because, like Scenario 2, output of variation of coupling bus 31 will make the gas load fluctuate in gas node 12, thus gas pressure of all gas nodes fluctuates in Scenario 3. Besides, the impact on the power system by the variation of a single generator in Scenario 3 is less than the variation of multiple loads in Scenario 2.

Analysis of power generator variation
In Scenario 3, the uncertainty is reflected in the power system and gas system, similar to the situation in Scenario 2. When the output of power generator fluctuates by uncertainty, all the results of electrical nodes will fluctuate, as well as the variation of all gas nodes. That is because, like Scenario 2, output of variation of coupling bus 31 will make the gas load fluctuate in gas node 12, thus gas pressure of all gas nodes fluctuates in Scenario 3. Besides, the impact on the power system by the variation of a single generator in Scenario 3 is less than the variation of multiple loads in Scenario 2.

Analysis of gas load variation
In Scenario 4, gas load variation will make IENGS fluctuate. Furthermore, it should be noted that the variation of the power system brought by the uncertainty in the gas system is less than the variation brought by the uncertainty of electrical load in Scenario 2. That is because in Scenario 4, when gas loads fluctuate, gas consumption of gas-fired generators will vary at gas nodes 4, 12, and 16, which will change the power generation in bus 30, 31, and 34 of the power system. Therefore, the uncertainty is delivered to the whole power system. Now that there are three uncertain power sources in Scenario 4, the variation in power system is more than variation brought by a single uncertain power source in Scenario 3.

Analysis of all uncertainties
When all three uncertainties are taken into consideration, the operating parameters U, δ, and v will fluctuate apparently to an extent. In addition, the variation of energy flow is much larger than Scenario 2-4 as is shown in Figure 5. To compare the uncertainty level of operating parameters in Scenario 5 with Scenario 2-4, the interval breadth ratios λ is proposed, which can be calculated by 3,4), where x is denoted as U, δ, or v. Solutions of Scenario 5 are set as reference values and λ x,i > 0. When λ x,i < 1, it means that fluctuation of the solution x in Scenario i is smaller than that in Scenario 5. The smaller λ x,i is, the lighter the variation of parameter x is in Scenario i compared with Scenario 5. When λ x,i > 1, it means that variation of the solution x in Scenario i is larger than that in Scenario 5.
The solutions of U, δ, or v in Scenarios 2, 3, and 4 are shown in Figure 6. Most values of λ U,i , λ δ,i and λ v,i (i = 2, 3, 4) are less than 1, from which it can be seen that the variation in Scenario 5 is larger than in Scenarios 2, 3, and 4. That is because the number of the uncertainty sources is larger than the number in Scenarios 2, 3, and 4, which can illustrate that accumulation of uncertainty sources will make the IENGS much more uncertain.
power source in Scenario 3.

Analysis of all uncertainties
When all three uncertainties are taken into consideration, the operating parameters U, δ, and v will fluctuate apparently to an extent. In addition, the variation of energy flow is much larger than Scenario 2-4 as is shown in Figure 5. To compare the uncertainty level of operating parameters in Scenario i is larger than that in Scenario 5. The solutions of U, δ, or v in Scenarios 2, 3, and 4 are shown in Figure 6. Most values of , Ui  , ,i   and , vi  (i = 2, 3, 4) are less than 1, from which it can be seen that the variation in Scenario 5 is larger than in Scenarios 2, 3, and 4. That is because the number of the uncertainty sources is larger than the number in Scenarios 2, 3, and 4, which can illustrate that accumulation of uncertainty sources will make the IENGS much more uncertain. In addition, from the comparison between Scenario 2 and Scenario 4, it can be concluded that the impact on gas system by the variation of electrical buses is much less than that of gas nodes, since only a gas node changes in Scenario 2. Similarly, the impact on gas system by the variation of single power generator is much less than that of gas nodes by the comparison between Scenario 3 and Scenario 4. In the same view, it can be seen that the impact on power system by the variation of gas nodes is much less than that of electrical buses by the comparison of Scenario 2 and Scenario 4.

Convergence Analysis of Uncertain Energy Flow
This part is aimed at the stability of the whole system when the uncertainty of a certain variable increases. Take Scenario 2, for example: the uncertainty level of electrical load is  5%. Here,  W% is expressed as the uncertainty level of electrical load, and the value of W is gradually increased to analyze the convergence of the proposed approach. As the variation of electrical load expands, the In addition, from the comparison between Scenario 2 and Scenario 4, it can be concluded that the impact on gas system by the variation of electrical buses is much less than that of gas nodes, since only a gas node changes in Scenario 2. Similarly, the impact on gas system by the variation of single power generator is much less than that of gas nodes by the comparison between Scenario 3 and Scenario 4. In the same view, it can be seen that the impact on power system by the variation of gas nodes is much less than that of electrical buses by the comparison of Scenario 2 and Scenario 4.

Convergence Analysis of Uncertain Energy Flow
This part is aimed at the stability of the whole system when the uncertainty of a certain variable increases. Take Scenario 2, for example: the uncertainty level of electrical load is ± 5%. Here, ± W% is expressed as the uncertainty level of electrical load, and the value of W is gradually increased to analyze the convergence of the proposed approach. As the variation of electrical load expands, the interval breadth of results becomes larger. Figure 7 shows [U], [δ], and [v] of uncertain energy flow when W is increased. It can be seen that as the W% increases, the convergence of the solution changes. When W is just over 40, the solution of uncertain energy flow becomes divergent. As bus 31 is the slack bus, and bus 30, 32-39 are PV buses, the voltage magnitudes of these nodes have the same value when W increases, which is no longer shown in Figure 7a. The voltage angle of bus 31 in Figure 7b and the gas pressure of node 1 in Figure 7c will remain unchanged when W changes, as bus 31 in the power system and node 1 in the gas system are slack nodes.

Conclusions
In this paper, an interval energy flow model considering uncertainties is built based on the steady-state energy flow model. The uncertainties include variation of electrical load, power generation, and gas load that may occur in the IENGS. To solve this model, the Krawczyk-Moore interval iterative method is used to calculate energy flow, and affine mathematics is applied in the calculation process to get less conservative solutions. The case study shows the effectiveness of the proposed approach by comparing it with MC simulation. By calculating energy flow in different cases, several conclusions can be made as follows: 1. Uncertainty caused by the power system (variation of electrical load or power generation) will be delivered to the gas system by coupling nodes. Similarly, uncertainty caused by the gas As the variation of electrical load becomes larger, the lower bound of power output in bus 31 will be 0. So, the gas-fired turbine in bus 31 will not deliver power, which means the turbine will stop operating. The gas consumption will not change at node 12 in the gas system. Therefore, the upper bound of gas pressure in the gas system will not change, while its lower bound will still change in Figure 7c.