Generalized Energy Flow Analysis Considering Electricity Gas and Heat Subsystems in Local-Area Energy Systems Integration

To alleviate environmental pollution and improve the efficient use of energy, energy systems integration (ESI)—covering electric power systems, heat systems and natural gas systems—has become an important trend in energy utilization. The traditional power flow calculation method, with the object as the power system, will prove difficult in meeting the requirements of the coupled energy flow analysis. This paper proposes a generalized energy flow (GEF) analysis method which is suitable for an ESI containing electricity, heat and gas subsystems. First, the models of electricity, heat, and natural gas networks in the ESI are established. In view of the complexity of the conventional method to solve the gas network including the compressor, an improved practical equivalent method was adopted based on different control modes. On this basis, a hybrid method combining homotopy and the Newton-Raphson algorithm was executed to compute the nonlinear equations of GEF, and the Jacobi matrix reflecting the coupling relationship of multi-energy was derived considering the grid connected mode and island modes of the power system in the ESI. Finally, the validity of the proposed method in multi-energy flow calculation and the analysis of interacting characteristics was verified using practical cases.


Introduction
The energy industry has gone through great change due to the rapid development of technology and the economy in the 21st century.In recent years, fossil energy exhaustion, environmental pollution, large-scale renewable energy application, etc. have proven to be key drivers in the development of a modern energy system which needs to meet the requirements of low-carbon emissions and efficient utilization of friendly sustainable energy sources.However, different energy systems are designed, planning and operating independently of the traditional style, which separates the coupling correlations between different types of energy and largely limits the flexibility of the energy system.In this case, ESI is proposed as it is an important type of new energy system, and covers electric power, heat and natural gas systems with the purpose of coordinating the whole system in the processes of generation, transmission and consumption, which are considered to be the core technology of the "third industrial revolution" [1].Therefore, the development of an integrated modeling and analysis framework for ESI represents an essential need for future research.

Heat source
Heat load Gas source Numerous functions in ESI will be largely different compared to conventional energy systems, especially the energy flow analysis function.Over the past few decades, power flow analysis, which is essential to analyze the steady-state as well as dynamic performance in power systems, has provided valuable references for its planning, operation, maintenance and protection.Similarly, in other energy systems, the quality of energy flow analysis contributes greatly to the efficiency of the operation and planning of any energy utility since many potentially costly decisions are based on such calculation results.With the advent of ESI, dramatic changes have occurred in energy utilization.Generalized energy flow (GEF), which analyzes the flow of various energy systems jointly as the cornerstone theory, has a tremendous guiding significance for energy suppliers, aggregators, energy marketers, independent system operators, etc.Hence, the steady state analysis of the multiple energy flow will play a key role in ESI.
Many studies have been done on energy flow and much work has been done concerning the interaction of various energy systems.Several conceptual approaches for modeling the integration of energy systems have been published, including energy hubs [4], multi-energy systems and distributed multi-generation [5], community energy [6], smart energy systems [7], and integrated energy systems [8].Reference [9] proposed a steady energy flow study to quantify the interdependency of energy infrastructures, which computed an equilibrium point of electricity and natural gas coupled systems by using an energy flow algorithm.In Reference [10], a combined analysis was developed to investigate the performance of electricity and heat networks as an integrated whole.Reference [11] proposed a new sequential calculation method to solve natural gas and power flow.In Reference [12] Figure 1.The relationship between wide-area energy systems integration (ESI) and local-area ESI.
Numerous functions in ESI will be largely different compared to conventional energy systems, especially the energy flow analysis function.Over the past few decades, power flow analysis, which is essential to analyze the steady-state as well as dynamic performance in power systems, has provided valuable references for its planning, operation, maintenance and protection.Similarly, in other energy systems, the quality of energy flow analysis contributes greatly to the efficiency of the operation and planning of any energy utility since many potentially costly decisions are based on such calculation results.With the advent of ESI, dramatic changes have occurred in energy utilization.Generalized energy flow (GEF), which analyzes the flow of various energy systems jointly as the cornerstone theory, has a tremendous guiding significance for energy suppliers, aggregators, energy marketers, independent system operators, etc.Hence, the steady state analysis of the multiple energy flow will play a key role in ESI.
Many studies have been done on energy flow and much work has been done concerning the interaction of various energy systems.Several conceptual approaches for modeling the integration of energy systems have been published, including energy hubs [4], multi-energy systems and distributed multi-generation [5], community energy [6], smart energy systems [7], and integrated energy systems [8].Reference [9] proposed a steady energy flow study to quantify the interdependency of energy infrastructures, which computed an equilibrium point of electricity and natural gas coupled systems by using an energy flow algorithm.In Reference [10], a combined analysis was developed to investigate the performance of electricity and heat networks as an integrated whole.Reference [11] proposed a new sequential calculation method to solve natural gas and power flow.In Reference [12] the concept of a probabilistic load flow, which has been widely applied in electric power systems, was extended to probabilistic energy flow analysis including electric power and natural gas flows.In addition to steady energy flow, plenty of studies have focused on optimal energy flow studies.In References [4,13], the optimal energy flow was solved by the centralized optimization algorithm and the distributed optimization algorithm, respectively, based on the energy hub concept.Reference [14] suggested a multi-objective algorithm to optimize its energy flow by considering the dynamics of the system and the equipment's thermal inertias.With respect to nonlinear algebraic equation solutions, one standard approach of solving the algorithm is the Newton-Raphson method [15].However, traditional power flow relying on Newton-Raphson is not robust if no starting point sufficiently close to a solution is provided.To address this issue, several methods were developed to reliably compute nonlinear equations, such as the Gauss Seidel iterations method [16], and the homotopy-enhanced method [17]; however, their convergence rate was relatively limited.Therefore, the improvement of existing approaches is still an urgent problem to be solved regardless of the mathematical model or solution algorithm.
By taking account of current scientific research and the engineering applications of ESI, the contributions of this paper can be summarized in three aspects as follows: (1) To the best of our knowledge, most previous research has concentrated either on the interaction of electricity-gas or electricity-heat.In this paper, the GEF of a local-area ESI is proposed, which designs the whole energy flow framework of gas, heat, and electricity jointly.(2) Complicated nonlinear mathematical models of some devices in energy systems create dilemmas in the implementation of energy flow calculations, especially regarding the compressor facility that belongs to gas energy system.By considering that issue, we effectively simplified the model as an equivalent transformation which reduces computational complexity.(3) By considering the advantage of the algorithm and the complexity of the ESI model, a hybrid algorithm, which combines the homotopy method and Newton's method, was adopted to solve nonlinear equations of GEF.
This paper is organized as follows.Section 2 presents the fundamental concept and mathematic model of each energy subsystem.Section 3 describes the GEF calculation method based on homotopy and the Newton-Raphson algorithm, and the Jacobi matrix, which maps multi-energy coupling relationships, is proposed.Section 4 shows the experimental results used to verify the effectiveness of the proposed algorithm in different scenarios.Finally, our conclusions and future research avenues are given in Section 5.

Generalized Energy Flow Model
ESI is composed of a power system, thermal system, natural gas system and multi-energy conversion equipment, which are distributed in local-area ESI (Figure 2).
Electrical energy flow transfers from electrical sources to the electrical load in power system networks.The district heating system, which is not suitable for long distance transmission, consists of a heat source, thermal load, supply and return networks, and delivers heat energy (in the form of hot water or steam) from the point of heat generation to the end consumers.The natural gas system is comprised of a gas source, gas supply pipe, a compressor and a gas load.Compressions are installed in gas pipelines to compensate for the loss of pressure due to both the friction of the pipes, as well as to provide the pressure needed to transport gas from one location to another [18].A combined heat and power (CHP) unit, electric boiler, gas boiler, and other elements make up the energy conversion facility to meet the demand of different energy sources.Additionally, the system is equipped with electrical energy storage, thermal energy storage and a gas storage station to improve economic and operational flexibility.The storage elements are modeled as sources or as loads depending on their operating condition.

Model of Steady State Power Flow in Power Subsystems
A classic AC power flow calculation method is employed in the steady-state electrical power flow model of a power subsystem in ESI.The steady state operation of a power system is formulated by stipulating that at each system's node, the power injected by generators, the power demanded by loads, and the power exchanged through the transmission elements connected to the node must add up to zero.The polar formulations of power flow can be denoted as per Reference [15] (1) where i and j are the index of electrical system nodes; Pgen.i and Qgen.i represent the generation of active and reactive power at the ith node, respectively; Pload.i and Qload.irepresent the demand of active and reactive power at the ith node; G is the conductance of the nodal admittance matrix; B is the susceptance of the nodal admittance matrix; Vi is voltage magnitude at ith node; and θi is the voltage angle at the ith node.

Model of Steady State Energy Flow in District Heat Subsystems
The model of the steady state heat flow consists of the hydraulic and thermal models.Usually, hydraulic analysis is conducted before thermal analysis [19].

Hydraulic Model
Similar to the node power balance equation in the power system, the mass flow that enters into a node is equal to the mass flow that leaves the node plus the flow consumption at the node., , q in q out q m m m  

Model of Steady State Power Flow in Power Subsystems
A classic AC power flow calculation method is employed in the steady-state electrical power flow model of a power subsystem in ESI.The steady state operation of a power system is formulated by stipulating that at each system's node, the power injected by generators, the power demanded by loads, and the power exchanged through the transmission elements connected to the node must add up to zero.The polar formulations of power flow can be denoted as per Reference [15]: where i and j are the index of electrical system nodes; P gen.i and Q gen.i represent the generation of active and reactive power at the ith node, respectively; P load.i and Q load.irepresent the demand of active and reactive power at the ith node; G is the conductance of the nodal admittance matrix; B is the susceptance of the nodal admittance matrix; V i is voltage magnitude at ith node; and θ i is the voltage angle at the ith node.

Model of Steady State Energy Flow in District Heat Subsystems
The model of the steady state heat flow consists of the hydraulic and thermal models.Usually, hydraulic analysis is conducted before thermal analysis [19].

Hydraulic Model
Similar to the node power balance equation in the power system, the mass flow that enters into a node is equal to the mass flow that leaves the node plus the flow consumption at the node.
∑ m q,in − ∑ m q,out = m q (3) where m q,in , m q,out are the mass flow within each pipe; m q is the vector of mass flow through each node injected from a source or discharged to a load; A h is the network incidence matrix, whose individual elements describes the directions of flow in a pipe; and m i is the corresponding mass flow.
The relation between the flow m ij and the head losses h ij along in pipe ij is: where κ ij is the vector of the resistance coefficients of each pipe and generally depends largely on the diameter of a pipe.The loop pressure equation states that the sum of head losses H around a closed loop must equal zero.
where B is the loop incidence matrix that relates the loops to the pipes; and h is the corresponding loop pressure.

Thermal Model
The thermal model is used to determine the general temperatures at each node.Heat power is calculated using Equation (7), where Φ ij is the vector of heat power consumed or supplied at each node; C p is the specific heat of water; and m ij is the vector of the mass flow rate through each node injected from a supply or discharged to a load; T s,i indicates the temperature before the hot stream is injected into the load; and T r,i indicates the temperature after the hot stream flow out of the load.
The temperature at the outlet of a pipe T end is calculated using the temperature drop Equation ( 8).
where T start and T end are the temperatures at the start node and the end node of a pipe; T a is the ambient temperature; λ is the overall heat transfer coefficient of each pipe per unit length; and L is the length of each pipe.The temperature of the water leaving a node with more than one incoming pipe is calculated as the mixture temperature of the incoming flows using Equation (9).
where T out is the mixture temperature of a node; m out is the mass flow rate within a pipe leaving the node; T in is the temperature of flow at the end of an incoming pipe; and m in is the mass flow rate within a pipe coming into the node.

Model of Steady State Energy Flow in Gas Subsystems
The steady-state modeling of a natural gas system is formulated by the equations related to the gas flow and thermal behavior of the gas in the elements comprising the system.The development of the gas flow equation in hydraulically smooth pipes is commonly found in numerous publications.The main difference between them being how the friction coefficient term and natural gas characteristics are considered in the formulation.In this context, the Weymouth equation is deployed to simulate compressible gas flow in pipelines which has been most widely used in the design of distribution gas networks no matter in commercial or academic field [20][21][22].
The conventional analysis approach of the natural gas steady state usually adopts the method of solving the nonlinear equations group jointly, which contains both the pipeline and compressor models [11].Due to the complexity of the compressor model, the calculation process always expends a large number of computing resources.Therefore, we isolated the pipeline branch with the compressor from the network, which separately models the network with or without the compressor.Specifically, the gas flow considering the pipeline with the compressor can be calculated directly by the method we employed which regards the two terminal nodes of the branch with the compressor as the equivalent load.This approach results in the mathematical model of the compressor no longer being involved in the network calculation, which greatly simplifies the calculation process.

Natural Gas Network Model without Compressor
The model without compressor is usually used to describe the relationship between the nodal pressure and the pipeline flow in the natural gas network.The steady state flow f r of the natural gas pipeline r can be defined as: where K r is the pipeline constant; ∆p r 2 = p i 2 − p j 2 represents the pressure drop of pipeline r; s ij is used to characterize the flow direction of natural gas, when p i > p j , s ij is +1, otherwise −1.D ij is diameter of the pipeline between node i and j. l ij is the length of the pipeline between i and j.Z a is the average compressibility factor; T a is average temperature of natural gas; and G is gas relative density.
The flow continuity equation of the natural gas network is expressed as: where A g is the node-branch incidence matrix of the natural gas network without the pipeline containing the compressor; f indicates the pipeline natural gas flow; L indicates the flow of each node Π i = p i 2 , ∆Π r = ∆p r 2 , then the pressure drop vector of the natural gas pipeline are noted as:

Natural Gas Flow Calculation Method with the Compressor
The pipeline considering the compressor driven by gas turbine is shown in Figure 3, where f com is the flow through the compressor; f cp is flow consumed by gas turbine itself; and f mi is the inlet flow of the compressor and f on is outlet flow of compressor.
Energies 2017, 10, 514 6 of 17 Specifically, the gas flow considering the pipeline with the compressor can be calculated directly by the method we employed which regards the two terminal nodes of the branch with the compressor as the equivalent load.This approach results in the mathematical model of the compressor no longer being involved in the network calculation, which greatly simplifies the calculation process.

Natural Gas Network Model without Compressor
The model without compressor is usually used to describe the relationship between the nodal pressure and the pipeline flow in the natural gas network.The steady state flow fr of the natural gas pipeline r can be defined as: where Kr is the pipeline constant; Δpr 2 = pi 2 − pj 2 represents the pressure drop of pipeline r; sij is used to characterize the flow direction of natural gas, when pi > pj, sij is +1, otherwise −1.Dij is diameter of the pipeline between node i and j. lij is the length of the pipeline between i and j.Za is the average compressibility factor; Ta is average temperature of natural gas; and G is gas relative density.
The flow continuity equation of the natural gas network is expressed as: where Ag is the node-branch incidence matrix of the natural gas network without the pipeline containing the compressor; f indicates the pipeline natural gas flow; L indicates the flow of each node Πi = pi 2 , ΔΠr = Δpr 2 , then the pressure drop vector of the natural gas pipeline are noted as:

Natural Gas Flow Calculation Method with the Compressor
The pipeline considering the compressor driven by gas turbine is shown in Figure 3, where fcom is the flow through the compressor; fcp is flow consumed by gas turbine itself; and fmi is the inlet flow of the compressor and fon is outlet flow of compressor.The node o represents the compressor outlet node, and the node i is the entrance node.Then, the mathematical model of the pipeline with the compressor can be expressed as: The node o represents the compressor outlet node, and the node i is the entrance node.Then, the mathematical model of the pipeline with the compressor can be expressed as: where k cp is the compression ratio; K mi , K on are constants of the inlet and outlet pipes; p m , p n , p i , p o are the correspondent nodal pressure; q gas is the natural gas heating value; T gas is the natural gas temperature; and a is the polytropic exponent.
In this paper, we propose a practical strategy to calculate pipeline flow when considering the compressor.The compressor control pattern can be divided into four categories, the details are shown in Figure 4: where kcp is the compression ratio; Kmi, Kon are constants of the inlet and outlet pipes; pm, pn, pi, po are the correspondent nodal pressure; qgas is the natural gas heating value; Tgas is the natural gas temperature; and a is the polytropic exponent.
In this paper, we propose a practical strategy to calculate pipeline flow when considering the compressor.The compressor control pattern can be divided into four categories, the details are shown in Figure 4: Pattern 1: When po is known, taking po, pn as input to calculate fon using the first line of Equation ( 14), which can be regarded as the initial value of fmi.By combining with pm to calculate pi and kcp, and using the second line of Equation ( 14) to calculate fcp, the new value of fmi can be obtained.Calculating iteratively until the variation of |fmi| is less than the mismatch tolerance ε.
Pattern 2: When pi is known, use the fourth line of Equation ( 14) to calculate fmi as the initial value of fon, before calculating po and kcp, using the second line of Equation ( 14) to calculate fcp.Update fon = fmi − fcp until the variation of |fon| is less than the mismatch tolerance ε.
Pattern 3: When fcom or fon are known, po can be calculated by the first line of Equation ( 14), which then transforms into Pattern 1.
Pattern 4: When kcp is known, the initial po is obtained easily, so pi = po/kcp.fmi and fon can be calculated by using the first and the fourth lines of Equation ( 14), respectively, so fcp = fmi − fon.The second line of Equation ( 14) is used to calculate fcom, namely fon, and then the first line of Equation ( 14) is used to calculate the new value of po.Calculate iteratively until the variation of |po| is less than ε.
The calculation process flow chart is shown in Figure 3. fmi and −fon can be obtained by the abovementioned calculation method as the equivalent load feq of node m and node n, respectively.

Analogue Analysis of GEF and Its Relevant Rules
The energy flow formulation of a district heat network or a gas network is similar to that of an electrical network.The AC electrical power flow model for electrical networks was well established in Reference [15], and the flow in other subsystems also possess the corresponding rules and laws according to its own characteristics, in order to further compare their differences and similarities of GEF.The rules in GEF can be summarized as follows in Table 1.Pattern 1: When p o is known, taking p o , p n as input to calculate f on using the first line of Equation ( 14), which can be regarded as the initial value of f mi .By combining with p m to calculate p i and k cp , and using the second line of Equation ( 14) to calculate f cp , the new value of f mi can be obtained.Calculating iteratively until the variation of |f mi | is less than the mismatch tolerance ε.
Pattern 2: When p i is known, use the fourth line of Equation ( 14) to calculate f mi as the initial value of f on , before calculating p o and k cp , using the second line of Equation ( 14) to calculate f cp .Update f on = f mi − f cp until the variation of |f on | is less than the mismatch tolerance ε.
Pattern 3: When f com or f on are known, p o can be calculated by the first line of Equation ( 14), which then transforms into Pattern 1.
Pattern 4: When k cp is known, the initial p o is obtained easily, so p i = p o /k cp .f mi and f on can be calculated by using the first and the fourth lines of Equation ( 14), respectively, so f cp = f mi − f on .The second line of Equation ( 14) is used to calculate f com , namely f on , and then the first line of Equation ( 14) is used to calculate the new value of p o .Calculate iteratively until the variation of |p o | is less than ε.
The calculation process flow chart is shown in Figure 3. f mi and −f on can be obtained by the above-mentioned calculation method as the equivalent load f eq of node m and node n, respectively.

Analogue Analysis of GEF and Its Relevant Rules
The energy flow formulation of a district heat network or a gas network is similar to that of an electrical network.The AC electrical power flow model for electrical networks was well established in Reference [15], and the flow in other subsystems also possess the corresponding rules and laws according to its own characteristics, in order to further compare their differences and similarities of GEF.The rules in GEF can be summarized as follows in Table 1.

GEF Calculation Method
Due to the simultaneous equations of different energy systems, the number of nonlinear equations and the dimension of decision variables are greatly increased, therefore, an inappropriate method used to solve the nonlinear equations will lead to bad convergence or even ill-conditioned flow.In this paper, we used a hybrid method to compute the equilibrium point of ESI, which is comprised of the homotopy method and the Newton-Raphson method.First, the homotopy algorithm was used to calculate the initial value of equations by its global convergence superiority, before the Newton-Raphson algorithm was executed to quickly converge to a stable solution, which met both requirements of computational stability and convergence speed.

GEF Mode
The equations of the GEF model are formulated by the equations of the electrical power flow, the heat flow, the gas flow and the coupling components model.
where the equations of the electrical power system represent active deviation and reactive power deviation; and the equations of the heat system represent the nodal heat power deviation, loop pressure drop deviation, temperature deviation in the heat supply system and the temperature deviation in the regenerative system.The equations of the gas system represent nodal flow deviation.P SP , Q SP , Φ SP and L SP are the given active power and reactive power load, heat load and natural gas load.A sl and A gl are the heat supply network correlation matrix and natural gas network correlation matrix, which has removed the gas compressor.C s and C r are matrix related to heat supply network and regenerative network, respectively.b s and b r are the column vector associated with the heat supply temperature and output temperature, respectively.Details are reported in Reference [13].X = [θ; V; m; T s,load ; T r,load ; Π] are the state variables of the GEF calculation.

Iterative Method
The homotopy analysis method is a numerical solution used to solve nonlinear equations, and can cooperate with many numerical methods/programs to achieve robust computation.It has the advantage of large or global convergence regions and is well suited for finding multiple solutions.Whether the original iteration point is sufficiently close to a provided solution has no influence on convergence.
The main idea is to embed a continuation parameter λ into F(x) of Equation ( 15) to form a new function S(x,λ) which satisfies the property that S(x,0) = G(x), S(x,1) = F(x).The function S(x,λ) is called a homotopy, where the system of equations S(x,λ) = G(x) at λ = 0 is easy to solve or has known solution(s) x 0 , and the system at λ = 1 is identical to the "difficult" problem that we aim to solve.The homotopy equation is defined by S(x,λ) = (1 − λ)G(x) + λF(x).
The Newton homotopy method is used in this case.When G(x) = F(x) − F(x 0 ), S(x,λ) = F(x) + (λ − 1)F(x 0 ), in order to obtain the homotopy path of the function x = x(t).The derivation of S(x,λ) can be noted as: The nonlinear equations are transformed into solving problems of differential equations.
The equation has a unique solution x = x(t), which meets the condition that x(0) = x 0 , x(1) = x* is the solution to Equation (15).Details of the homotopy algorithm are reported in Reference [17].
After the initial value was calculated by the homotopy algorithm, the Newton-Raphson method was used to solve the GEF problem, and the iterative equations can be denoted as follows: where ∆F e = [∆P, ∆Q], ∆F h = [(∆Φ, ∆p), (∆T s , ∆T r )], ∆F g = ∆f represent deviations of electricity, heat and gas, respectively.x e = [θ, V], x h = [m, (T s,load , T r,load )], x g = Π represent state variables of electricity, heat and gas, respectively.The Jacobian matrix of GEF can be noted as: J ee J eh J eg J he J hh J hg J ge J gh J gg where diagonal block matrix J ee , J hh , J gg represent the relationship between the energy flow and its own state variable in electric, thermal, gas system, respectively, and the mathematical expression is the same as the calculation for traditional power, thermal, and natural gas flows.The derivation method of the diagonal block expression is given in References [10,11,15].
In a natural gas network, the slack node is always connected with the gas source, and as the state of the natural gas system changes, the supply and demand fluctuations will be borne by the change of the slack node supply, which does not affect the state of the thermal and power systems, so J eg and J hg in Equation ( 19) are zero term.
In a heat network, the thermal power in slack node is partially supplied by the CHP unit which works in following the thermal load (FTL) model; as the state of the thermal system changes, fluctuations in the thermal power slack node can also change the consumption of natural gas and electricity simultaneously, so J eh , J gh in Equation ( 19) are non-zero term.The output of the CHP unit and gas consumption is denoted by Equation ( 20): where A source,i is the row which is related to the heat source in the node correlation matrix of the thermal system.In electric power grid, operation mode can be divided into two types: (1) Grid Connected Mode When the electrical subsystem in the ESI is connected to the bulk power grid, the fluctuation of the electric power in the ESI is balanced by the bulk power grid.
(2) Island Mode When the electrical subsystem in the ESI operates in island mode, as the state of thermal system changes, the fluctuation is balanced by the CHP unit which works in following the electric load (FEL) model in electric power slack node can also change the consumption of natural gas and heat power simultaneously, so J he , J ge in Equation ( 19) are non-zero term.
The CHP output and the gas consumption of the slack node in the island mode power system are expressed as: As per the analysis above, Equations ( 20) and ( 21) are substituted into the equations of the GEF.After computing the first derivation of GEF, the non-diagonal block of the Jacobi matrix under the condition of grid connected mode and island mode is shown in Table 2.
Table 2. Expressions for non-diagonal blocks of the Jacobi matrix.

Case Studies
The suitability of the proposed approach for conducting steady-state GEF studies was tested on a system network located in the Daxing District, Beijing, China as the ESI demonstration project of Goldwind Techology Co. Ltd.The different energy systems were coupled via two CHP units (in Figure 5).EBi, HBi and GBi served as the nodes of the power, thermal, and gas systems, respectively.A gas compressor was installed in the gas subsystem due to the heavy gas load in node GB6 to ensure pressure and gas supply.For line constants of the ESI, refer to Appendix A.

Simulation of Natural Gas Subsystems
To further verify the efficiency of the proposed gas compressor equivalent approach: First, the natural gas system compressor was regarded as the research object, and it was assumed that the load capacity of GB5 and GB6 were 1800 m 3 /h and 2000 m 3 /h, respectively.We compared the simulation results of our proposed method in Section 2.3, and the method in Reference [19] under the condition that the compressor outlet pressure was known, and the simulation times were 0.0030 s and 0.0055 s, respectively.The simulation result of the flow through the compressor was 3258.8 m 3 /h.In addition, the case in Reference [19] included 22 nodes in the natural gas system with compressor to test, and the simulation times were 0.0062 s and 0.0157 s, respectively.The method proposed in this paper is obviously accurate and fast.The compressor of the natural gas subsystem in Figure 5 was simulated in four modes.The results are shown in Table 3.

Calculation Results Analysis of ESI in Two Scenarios
 Scenario 1 When the ESI power system operates in grid connected mode, node EB13 is connected to the bulk power grid as the slack node.EB11 and EB2 are the PV nodes, the CHP1 unit works in FEL mode, the CHP2 unit works in FTL mode, and HB13 is the thermal system slack node.


Scenario 2 When the ESI power system operates in island mode, the EB13 node is connected to a wind generator with a known active power and voltage amplitude.The working mode of CHP1 and CHP2 is the same as in Scenario 1, and EB11 is the slack node of the power system.
The results obtained for all cases are reported in Tables 4-7.The electrical, thermal power and natural gas consumption of CHP units are demonstrated in Figure 6.Conclusions can be drawn from the simulation results from Scenarios 1 and 2. In the power system, the amplitude of node voltage varied from 0.94-1.06p.u.The equilibrium point obtained for each case was located within the safety

Simulation of Natural Gas Subsystems
To further verify the efficiency of the proposed gas compressor equivalent approach: First, the natural gas system compressor was regarded as the research object, and it was assumed that the load capacity of GB5 and GB6 were 1800 m 3 /h and 2000 m 3 /h, respectively.We compared the simulation results of our proposed method in Section 2.3, and the method in Reference [19] under the condition that the compressor outlet pressure was known, and the simulation times were 0.0030 s and 0.0055 s, respectively.The simulation result of the flow through the compressor was 3258.8 m 3 /h.In addition, the case in Reference [19] included 22 nodes in the natural gas system with compressor to test, and the simulation times were 0.0062 s and 0.0157 s, respectively.The method proposed in this paper is obviously accurate and fast.The compressor of the natural gas subsystem in Figure 5 was simulated in four modes.The results are shown in Table 3.

Calculation Results Analysis of ESI in Two Scenarios
• Scenario 1 When the ESI power system operates in grid connected mode, node EB13 is connected to the bulk power grid as the slack node.EB11 and EB2 are the PV nodes, the CHP1 unit works in FEL mode, the CHP2 unit works in FTL mode, and HB13 is the thermal system slack node.

•
Scenario 2 When the ESI power system operates in island mode, the EB13 node is connected to a wind generator with a known active power and voltage amplitude.The working mode of CHP1 and CHP2 is the same as in Scenario 1, and EB11 is the slack node of the power system.
The results obtained for all cases are reported in Tables 4-7.The electrical, thermal power and natural gas consumption of CHP units are demonstrated in Figure 6.Conclusions can be drawn from the simulation results from Scenarios 1 and 2. In the power system, the amplitude of node voltage varied from 0.94-1.06p.u.The equilibrium point obtained for each case was located within the safety zone of operation, whose voltage limits were not violated.Node EB10, which has a heavy load, is far away from the power source node, so the node voltage is relatively lower than the others.In the heat system, due to the presence of two heat sources, the minimum temperature of the node occurs at the middle position of the heat system, which is in node HB7.In the gas system, the compressor enables each node to reach the supply and demand balance based on pressure guarantee.In these two scenarios, the compression ratio of gas compressor were 1.2938 and 1.294, respectively, which can be calculated by Equation (12).Although the electrical, thermal and gas load of the two scenarios had no change (on account of the electrical power of the slack nodes in the power system affecting the heat source power in the thermal system in Scenario 2), the two sets of CHP electrical, thermal power and natural gas consumption were changed at the same time.

Convergence Analysis
For the sake of testing and verifying the convergence of the GEF calculation approach, the mean absolute value of Δx was selected as the convergence index.We compared the hybrid method with the Newton-Raphson method to analyze the GEF under two scenarios.
The GEF calculated by the Newton-Raphson, which was initiated randomly, needs 12 times and 64 times to converge to the iterative solution in these two scenarios.However, the proposed method whose state variables were initialized by the homotopy algorithm as per the guidelines given in Section 3.2 only took 5 times and 30 times to converge to the equilibrium point which meets the mismatch tolerance, respectively.The convergence curve of each algorithm in the two scenarios is presented in Figure 7.

Convergence Analysis
For the sake of testing and verifying the convergence of the GEF calculation approach, the mean absolute value of ∆x was selected as the convergence index.We compared the hybrid method with the Newton-Raphson method to analyze the GEF under two scenarios.
The GEF calculated by the Newton-Raphson, which was initiated randomly, needs 12 times and 64 times to converge to the iterative solution in these two scenarios.However, the proposed method whose state variables were initialized by the homotopy algorithm as per the guidelines given in Section 3.2 only took 5 times and 30 times to converge to the equilibrium point which meets the mismatch tolerance, respectively.The convergence curve of each algorithm in the two scenarios is presented in Figure 7.
In contrast, Figure 7 shows that the iterative number calculated by either algorithm to achieve convergence under grid connected mode was far less than that of island mode.When the power system of the ESI operates in island mode, the thermal power of the slack node in the heat system affects the related electrical nodal injection power and the related gas nodal load.Similarly, the electrical power of the slack node in the power system influences its relevant nodal injection power in the heat system and the nodal load in the gas system.Therefore, the state values of the relevant nodes in the GEF model change after each iteration in island mode, which makes the convergence rate slow and the number of iteration increases.In contrast, Figure 7 shows that the iterative number calculated by either algorithm to achieve convergence under grid connected mode was far less than that of island mode.When the power system of the ESI operates in island mode, the thermal power of the slack node in the heat system affects the related electrical nodal injection power and the related gas nodal load.Similarly, the electrical power of the slack node in the power system influences its relevant nodal injection power in the heat system and the nodal load in the gas system.Therefore, the state values of the relevant nodes in the GEF model change after each iteration in island mode, which makes the convergence rate slow and the number of iteration increases.

Interaction Analysis of ESI
To analyze the interactive characteristics of ESI, in Scenario 2, the heat load at node HB1 and node HB2 was increased to 0.6 MW and 0.4 MW, respectively.The proposed method of GEF was applied to analyze the change of ESI state, which had bidirectional interaction between the electrical and thermal systems in this scenario.Figure 8a,b demonstrate the change of supply temperature and bus voltage at each node under the condition of load variation in the heat system.Due to the energy mutual interaction between the power and thermal systems in ESI, it is difficult to obtain the interdependency state of ESI without the proposed approach to analyze the GEF.In Figure 8a,b, we can see that as the heat load increased, the state of the power and heat systems changed simultaneously.In the power system, the nodal voltage, which was adjacent to EB12, decreased slightly.In the heat system, due to the change of hot water flow distribution, the heat temperature of node HB2, HB7 and HB8 rose, with the rising value reaching approximately 0.5 °C.

Interaction Analysis of ESI
To analyze the interactive characteristics of ESI, in Scenario 2, the heat load at node HB1 and node HB2 was increased to 0.6 MW and 0.4 MW, respectively.The proposed method of GEF was applied to analyze the change of ESI state, which had bidirectional interaction between the electrical and thermal systems in this scenario.Figure 8a,b demonstrate the change of supply temperature and bus voltage at each node under the condition of load variation in the heat system.In contrast, Figure 7 shows that the iterative number calculated by either algorithm to achieve convergence under grid connected mode was far less than that of island mode.When the power system of the ESI operates in island mode, the thermal power of the slack node in the heat system affects the related electrical nodal injection power and the related gas nodal load.Similarly, the electrical power of the slack node in the power system influences its relevant nodal injection power in the heat system and the nodal load in the gas system.Therefore, the state values of the relevant nodes in the GEF model change after each iteration in island mode, which makes the convergence rate slow and the number of iteration increases.

Interaction Analysis of ESI
To analyze the interactive characteristics of ESI, in Scenario 2, the heat load at node HB1 and node HB2 was increased to 0.6 MW and 0.4 MW, respectively.The proposed method of GEF was applied to analyze the change of ESI state, which had bidirectional interaction between the electrical and thermal systems in this scenario.Figure 8a,b demonstrate the change of supply temperature and bus voltage at each node under the condition of load variation in the heat system.Due to the energy mutual interaction between the power and thermal systems in ESI, it is difficult to obtain the interdependency state of ESI without the proposed approach to analyze the GEF.In Figure 8a,b, we can see that as the heat load increased, the state of the power and heat systems changed simultaneously.In the power system, the nodal voltage, which was adjacent to EB12, decreased slightly.In the heat system, due to the change of hot water flow distribution, the heat temperature of node HB2, HB7 and HB8 rose, with the rising value reaching approximately 0.5 °C.Due to the energy mutual interaction between the power and thermal systems in ESI, it is difficult to obtain the interdependency state of ESI without the proposed approach to analyze the GEF.In Figure 8a,b, we can see that as the heat load increased, the state of the power and heat systems changed simultaneously.In the power system, the nodal voltage, which was adjacent to EB12, decreased slightly.In the heat system, due to the change of hot water flow distribution, the heat temperature of node HB2, HB7 and HB8 rose, with the rising value reaching approximately 0.5 • C.

Conclusions
This paper proposed a generalized energy flow analysis method in energy systems integration.The method of modeling the pipeline with a compressor as the equivalent load was used to reduce the computational complexity.In addition, the set of nonlinear equations of GEF were obtained based on rules of each energy flow, with the hybrid method based on the homotopy method and the Newton-Raphson algorithm, executed to solve the nonlinear equations, which harness respective strengths.Finally, the simulation results showed that the proposed algorithm could converge to numerical solutions at a reasonable speed while the power system of ESI runs in grid connected or island modes.Due to the energy interaction between the power and thermal systems in ESI in island mode, more iterative times were needed for calculation.In contrast, GEF could effectively analyze the whole flow distribution variation when one of the energy systems changed its state.

Figure 1 .
Figure 1.The relationship between wide-area energy systems integration (ESI) and local-area ESI.

Figure 3 .
Figure 3. Pipeline with compressor driven by gas turbine.

Figure 3 .
Figure 3. Pipeline with compressor driven by gas turbine.

Figure 4 .
Figure 4. Mass rate calculating flow chart of pipeline with compressor.

Figure 4 .
Figure 4. Mass rate calculating flow chart of pipeline with compressor.

Figure 5 .
Figure 5. Case of energy system integration.

Figure 5 .
Figure 5. Case of energy system integration.

Figure 6 .
Figure 6.Electricity power, heat power and natural gas consumption of CHP.

Figure 6 .
Figure 6.Electricity power, heat power and natural gas consumption of CHP.

Figure 7 .
Figure 7. (a) Convergence of the calculation method in island mode; (b) Convergence of the calculation method in grid connected mode.

Figure 8 .
Figure 8.(a) Voltage change of load nodes in the power system; (b) Supply temperature change of load nodes in the heat system.

Figure 7 .
Figure 7. (a) Convergence of the calculation method in island mode; (b) Convergence of the calculation method in grid connected mode.

Figure 7 .
Figure 7. (a) Convergence of the calculation method in island mode; (b) Convergence of the calculation method in grid connected mode.

Figure 8 .
Figure 8.(a) Voltage change of load nodes in the power system; (b) Supply temperature change of load nodes in the heat system.

Figure 8 .
Figure 8.(a) Voltage change of load nodes in the power system; (b) Supply temperature change of load nodes in the heat system.

Table 1 .
Analogue analysis of generalized energy flow (GEF) rules and principles.

Table 3 .
Simulation results of the compressor.

Table 3 .
Simulation results of the compressor.

Table 4 .
Voltage magnitude of nodes in the power system.

Table 5 .
Node pressure of the gas network.

Table 6 .
Pipeline mass rate of the gas network.

Table 7 .
Supply temperature and return temperature of nodes in the heat system.

Table 7 .
Supply temperature and return temperature of nodes in the heat system.

Table A2 .
Coefficients of the power system.

Table A3 .
Load data of the heat system.