Multi-Energy Flow Calculation Considering the Convexification Network Constraints for the Integrated Energy System

To alleviate environmental pollution and improve the energy efficiency of end-user utilization, the integrated energy systems (IESs) have become an important direction of energy structure adjustment over the world. The widespread application of the coupling units, such as gas-fired generators, gas-fired boilers, and combined heat and power (CHP), increases the connection among electrical, natural gas, and heating systems in IESs. This study proposes a mixed-integer nonlinear programming (MINLP) model combining electrical, natural gas, and heating systems, as well as the coupling components, such as CHP and gas-fired generators. The proposed model is applicable for either the radial multi-energy network or the meshed multi-energy network. Since the proposed MINLP model is difficult to be solved, the second-order cone and linearized techniques are used to transform the non-convex fundamental matrix formulation of multi-energy network equations to a mixed-integer convex multi-energy flow model, which can improve the computational efficiency significantly. Moreover, the potential convergence problem of the original model can also be avoided. A simulation of IEEE 14-node electrical system, 6-node natural gas system, and 23-node heating system are studied to verify the accuracy and computational rapidity of the proposed method.


INTRODUCTION
With the rapid development of the economy, energy and environmental problems have become increasingly prominent. How to achieve clean and efficient use of energy has become the focus of research in recent years. The integrated energy system (IES) OMalley and Kroposki, 2013) incorporates the production, transmission, distribution, conversion, storage, and consumption of many kinds of energy; can realize the comprehensive management and economic dispatching of electricity, heat, gas, etc.; and provides an essential solution for learning the total utilization of energy. The energy efficiency of natural gas-based combined heat and power (CHP) units (Yang et al., 2010) is more than 80%. It is an efficient and environmentally friendly energy supply mode and has become an integral coupling unit among electric, gas, and heating networks. Under the IES, all kinds of energy conversion equipment, such as CHP, gas turbine, and gas boiler, make electricity, heat, and nature closely coupled and realize the interaction and conversion of multi-energy. The integrated energy system recognizes the exchange and transformation of thermal/electric/gas energy, but the coupling of the three energy sources has dramatically changed the system's trend. How to effectively calculate the distribution of multi-energy flow (multi-energy flow, MEF) (Pan et al., 2016) is of great significance to guide the investment planning and operation decision of IESs.
At present, aiming at the problems related to MEF, the research at home and abroad is mainly focused on the joint analysis of electricity/gas or electricity/heat energy networks. In the study by Zhang (2005), a sequential solution of hybrid power flow is proposed by combining the existing natural gas hydraulic calculation method with the power system power flow calculation method, and the energy concentrator model is established in the study by Geidl and Andersson (2007), Arnold et al. (2008), Geidl (2007). The centralized optimization algorithm and distributed optimization algorithm are used to solve the electric/gas hybrid optimal power flow, respectively.
From the point of view of the reliability of energy supply, the transmission delay and compressibility of natural gas are considered, and the optimal short-term operation of the electricity/gas coupling system is studied in reference Correa-Posada and Sanchez-Martin (2015). In reference Gu et al. (2015), an optimization model of the electro-thermal energy integrated system considering the constraints of power network and thermal network is established, and the benefit of wind power heating is studied. The research significance, application prospect, and critical technologies of the electric heating combined system with large capacity heat storage are reviewed in reference Xu et al. (2014). Considering the internal coupling of MEF, the decomposition algorithm of electric/thermal/gas hybrid optimal power flow based on energy hub is proposed in reference Moeini-Aghtaie et al. (2014), Shabanpour-Haghighi andSeifi (2015), but the method cannot guarantee the optimal global solution. In reference Xu et al. (2015), a hierarchical energy management model of the regionally integrated energy system is established by considering timescale and network constraints, but the consideration of thermal part is limited to adjustable heat load.
The above research shows that MEF computing under IESs has been widely concerned, but there are still the following problems: 1) Most of the research objects are electricity/gas systems or electricity/thermal systems, and there is a lack of research on electricity/gas/thermal interconnection systems. 2 In the steady natural flow model, the velocity at the inlet and outlet of the gas pipeline is the same, and there is a quadratic Weymouth function relationship between the velocity of the natural gas pipeline and the pressure difference of the gas pipeline. However, the Weymouth equation is non-convex and nonlinear, which brings difficulties to the MEF calculation.
3) The original heating network power flow equation is a nonlinear equation. The coupling relationship between temperature and flow is strong and contains an exponential equation, making the computational complexity high so that the numerical stability is difficult to be guaranteed. 4) At present, the research on IESs is mainly focused on the distribution of energy in the system, and there is a lack of research on the coupling interaction of electricity/heat/gas systems.
In order to solve the above problems, the MEF calculation method of IESs with electricity, heat, and gas is studied in this article. First, the modeling of many kinds of electro-thermal coupling units such as CHP and the gas turbine is studied, and the mathematical models of subsystems and coupling links in IESs are established, and the Weymouth equation is linearized reasonably by making use of the characteristic of short pipeline in the natural gas network. The traditional method of solving power flow in a power system is improved to establish a model which is easier to solve. As for the heating system, a method based on Taylor's second-order expansion is implemented to avoid the nonlinear equation in this article. Considering different operation modes of CHP units, two models of cogeneration are established, including the backpressure model and the pumping model. On this basis, the multi-energy flow solution model of the joint electric/ heating/gas network is established, and the practicability and rapidity of the method proposed in this article are proven by practical examples.

SYSTEM MODELING
The integrated energy system with electricity, gas, and heat is composed of a power system, thermal system, natural gas systems, and the coupling units such as CHP, gas turbine, and gas boilers.
The power system mainly includes generator, electric load, and transmission line; the thermal system mainly comprises a heat source, heat load, supply, and reflux pipeline; and the natural gas system includes explicitly gas source, gas load, and gas transmission pipeline.
The classification and variables of each system node are shown in Table 1.
For electrical networks, P, Q, V, and θ are the active power, reactive power, voltage amplitude, and phase angle of the power system nodes, respectively. In the natural gas network, f and π are the natural gas injection flow and pressure of the nodes, respectively. T s , T r , ϕ, and m are the supply temperature,

System type
Node type Known Unknown f π reflux temperature, heat power, and water flow of the heat source or heat load of the thermal system, respectively.

Heating Network
The steady-state power flow calculation model of the heating network is divided into two parts: hydraulic model and thermodynamic model (Liu, 2013).

Hydraulic Model
The flow of hot water in the network should meet the fundamental law of the network: the flow of each pipeline should satisfy the flow continuity equation at each node, that is, the injection flow at the node is equal to the outflow; in a closed-loop composed of pipes, the sum of the head loss of water flowing in each pipeline is 0, that is, where A s is the node-branch correlation matrix of the heating network, m is the flow of each pipeline, m q is the flow out of each node, B h is the loop-branch correlation matrix of the heating network, h f is the head loss vector, and its calculation method is as follows: where K is the resistance coefficient matrix of the pipeline.

Thermodynamic Model
For each heat load node, the heating temperature T s represents the temperature before the hot water is injected into the load node, the output temperature T o represents the temperature when the hot water flows out of the load node, and the heat recovery temperature T r represents the temperature after the hot water from the load node is mixed with the water of other pipes at the pipe node. The thermodynamic model is as follows: Equation 3 is the expression of the node thermal power Φ of the thermal network, and C p is the specific heat capacity of water; Eq. 4 represents the relationship between the temperature at the end of the pipe T end and the temperature T start at the beginning, T a is the ambient temperature, λ is the heat conduction coefficient of the pipe, and L is the length of the pipe; Eq. 5 represents the temperature relationship before and after the hot water is mixed at the node, and m out , T out , m in , T in represent the flow rate and temperature of the water flowing out and into the pipe, respectively.

Natural Gas Network
A natural gas system can be specified by several equations related to various elements of this system, including pipelines, compressors, sources, and loads. The input-output flow balance of each node should be considered for a feasible operational condition. The amount of gas flow through a pipeline connected between nodes i and j in standard cubic meter per hour (SCM/h) can be expressed by (De Wolf and Smeers, 2000) the following: where π i and π j are pressures at the nodes i and j, respectively. sign(π i , π j ) is a sign function of pressures where its value is +1 if π i > π j and −1 otherwise. The pipeline constant can be computed by Eq. 7, where z is the gas compressibility factor, T is the gas temperature, and δ is the density of gas relative to air. D ij and L ij are the interior diameter and length of the pipeline, respectively. λ ij is the friction factor of the gas pipeline. It is computed by Eq. 8, where ε g is absolute rigidity of the pipe. The node flow balance equation of the natural gas network is as follows: where f l is the natural gas flow consumed by the natural gas load, including gas turbine, gas boiler, CHP, and conventional natural gas load. For the compressor inlet node, it should also include the compressor outlet flow, that is, the compressor outlet flow as the load of the inlet node, f s is the injection flow of the gas source, and f ij is the flow of all non-compressor branches connected with the node.

Electric Network
The power system model adopts the AC system model, and the electric load of the heating system and the natural gas system is taken into account. The balance equation of active power and reactive power is as follows: where P G and Q G are the active power and reactive power generated by the generator, respectively; P l and Q l are the active power and reactive power consumed by the load in the electric network, respectively.

Coupling Units
The coupling units of electricity, gas, and heat integrated energy systems mainly include CHP, gas boiler, gas turbine, electric compressor, and water pump. CHP is the main coupling unit in electric, gas, and thermal systems. CHP is a kind of unit which not only generates electricity by a steam turbine but also supplies heat to thermal users by steam after power generation. The ratio of heat generation power to electricity generation power of CHP with gas turbine and reciprocating internal combustion engine in prime mover can be regarded as a constant, and the ratio of heat to electricity can be expressed as follows: where P CHP and ϕ CHP are the power generated and thermal power of CHP, respectively, and c m is the thermoelectric ratio of CHP. According to the working mode, CHP can be divided into back pressure type and pumping type. The heat-power ratio of the back pressure CHP unit is constant, while the extraction-condensation CHP unit changes the heat-power ratio of the CHP unit by adjusting the amount of steam extracted, and its output is mainly related to the natural gas flow consumed, that is, where f CHP is the natural gas flow consumed by CHP and q gas is the calorific value of natural gas. η e and η h represent the power generation efficiency and heat generation efficiency of CHP, respectively. The gas turbine is the equipment that consumes natural gas to generate electric energy. The relationship models between the consumed natural gas flow and the generated power are as follows: where f G is the natural gas flow consumed by the gas turbine, P G is the electric power generated by the gas turbine, and a G , b G , and c G are consumption coefficients of the gas turbine.

MEF Model
Based on the above, an IES-oriented MEF model is constructed as follows: In Eq. 15, the first row and the second row represent the active power deviation and reactive power deviation of the power system, respectively. The third to the sixth lines represent the node thermal power deviation of the thermal system, the pressure drop deviation of the heating network loop, the heating temperature deviation, and the regenerative temperature deviation, respectively. The seventh line represents the node flow deviation of the natural gas system. P G , Q G , Φ, and f l are the active power, reactive power, thermal power, and natural gas load given by the system, respectively. P L and Q L are the total active load and reactive load, respectively. A s is the incidence matrixes of the heating network. C s and C r are matrices related to the structure and flow of heating network and regenerative network, respectively, while b s and b r are column vectors related to heating temperature and output temperature, respectively.
It is worth noting that the power balance equation of power network includes trigonometric function, the flow balance equation of gas network involves non-convex nonlinear pipeline flow equation, and the heat network temperature balance equation includes exponential equation so that the IES-oriented MEF model is a non-convex nonlinear problem, which makes it difficult for us to solve MEF directly with traditional methods, and the calculation accuracy is difficult to be guaranteed.
Therefore, in this article, the convex optimization of the nonconvex nonlinear equations will be carried out below, which can not only ensure accuracy but also greatly simplify the calculation process and shorten the calculation time.

Rectangular Formulation
Let Y denote the nodal admittance matrix, which has components Y ij G ij + iB ij for each line. The complex voltage (also called voltage phasor) V i at bus i can be expressed in the rectangular form as V i e i + if i .
With the above notation, the power flow conservation at each bus is given in the so-called rectangular formulation as follows: Generation and voltage bounds at each bus are as follows: Here, P l i , Q l i is the summary of all kinds of active and reactive load in the electric system.

Second-Order Cone Program Relaxation of Alternative Formulation
Note that the rectangular formulation of AC power flow is a nonconvex quadratic optimization problem. However, quite importantly, we can observe that all the nonlinearity and nonconvexity come from one of the following three forms: (1) To capture this nonlinearity, we define new variables C, S to get rid of the variable e and variable f (Kocuk et al., 2016). The new variables are defined as follows: Frontiers in Energy Research | www.frontiersin.org July 2021 | Volume 9 | Article 718151 C ij e i e j + f i f j , S ij e i f j − e j f i .
With a change of variables, we can introduce an alternative formulation of the power flow problem in the electric system as follows: Through the above transformation, we convex the original Cartesian coordinate Eqs 16-20 into Eqs 24-29, which is more convenient to solve.

Natural Gas Network
The original Weymouth equation, which represents the relationship between the natural gas flow with the pressure at the inlet and the outlet of a natural gas pipeline, is a non-convex and nonlinear problem and hard to solve directly. Two approximation methods are used to linearize the Weymouth equation.
Method ①: One-dimensional approximation. From Eq. 6, it can be found that the right side of Eq. 6 is a function of π 2 i and π 2 j . By introducing variables φ i π 2 i and φ j π 2 j , Eq. 6 can be replaced with (Zhou, 2020) the following: In addition, the node pressure constraint can also be replaced with the following: It can be seen that constraint Eq. 30 is a one-dimensional nonlinear equation, which greatly simplifies the linearization process. The upper and lower limits of φ ij can be determined by the following formula: The range [φ min ij , φ max ij ] φ ij can be divided into N ij segments, and the nonlinear Eq. 30 can be transformed into Eqs. 36-39: (36) where p is the one-dimensional piecewise linear segment number index, z ij,p is the binary variable used in the one-dimensional linear approximation method, and K ij p is the constant of the onedimensional linear equation.
Method ②: Taylor expansion approximation. We linearize Eq. 6 by Taylor expansion (Manshadi and Khodayar, 2015) as follows: where π 0 i is the initial pressure value of the node i at the beginning of the optimization period and is a known quantity.
The linearization using the Taylor series is valid only if the difference in natural gas pressure between the inlet and outlet of the pipeline is assumed to be limited, that is, there is no significant pressure drop in the pipeline. This is a reasonable assumption for the short pipelines used in microgrids. The limitation on the node pressure on the gas pipeline network guarantees the accuracy of the approximation. The Weymouth equation is linearized around the initial point procured by solving optimal energy flow within the microgrid considering no disruptions.
By Method ①, we replace the non-convex and nonlinear pipeline Eqs 6-9 with the linear model consisting of Eqs 36-39. By Method ②, we transform the original Weymouth Eq. 6 into the linear model of Eq. 40. Both of them are mature methods and help to solve the problem quickly.

Heating Network
The heating network model studied in this article is the radiant heat network model. The original heating network power flow in Eqs 1-5 is a nonlinear equation. The coupling relationship between temperature and flow is strong and contains an exponential equation, making the computational complexity high so that the numerical stability is difficult to be guaranteed.
Therefore, this study adopts the method (Sun et al., 2020) based on Taylor's second-order expansion; the specific contents are as follows: To get the flow m i of the pipe i in Figure 1, we define n 1 m 1 /m i , n 2 m 2 /m i ···n i−1 m i−1 /m i , n i 1; thus, When calculating n 1 n 2 ,···,n i−1 , the heat loss of the pipe network is ignored, and the heat energy flowing through the pipe i is set to φ i . While to figure the pipe flow m i into the load node, we set n 1 ≈ as the temperature change of the water supply network node is small. k is Frontiers in Energy Research | www.frontiersin.org July 2021 | Volume 9 | Article 718151 the number of the pipe through which the flow m i flows and φ k is the heat energy flowing through the pipe k. Thus, we can define Then the Eq. 42 can be simplified as follows: The value of C is generally less than 10 −4 orders of magnitude, then formula (42) can be approximated to In this article, we define Eqs 42, Eqs 43 as model 2 and model 1, respectively.
The flow rate of each pipeline flowing into the load node can be calculated directly by Eq. 42 or Eq. 43, and the supply temperature, return temperature, and pipe flow rate of each node of the heating network can be obtained according to Eqs 1, 3-5.
Model 1 and model 2 decouple the water supply network and the backwater network, decouple the temperature and flow, the model is simple, the amount of calculation is small, and the derivation is based on Taylor's second-order expansion, and the solution accuracy is higher.
Through the above transformations and simplifications, the original MINLP (mixed-integer nonlinear programming) is reduced to a MICP (mixed-integer convex programming) problem, which can be effectively solved by commercial solvers.

Case Introduction
In the example of the integrated energy system used in this article, the power system is an IEEE standard 14-bus system (Kersting, 1991), as shown in Figure 2, in which the generator node is replaced by CHP. In order to match the capacity and power of CHP, the load is increased to 3.6 times, and the active load with a certain capacity is arranged at the balance node. The natural gas system is a 6-node system, and its line parameters refer to the study by Cong Liu et al. (2009).
The thermal system adopts the 23-node system in the reference Sun et al. (2020). The topological structure of each system is shown in Figures 2-4. In this article, the FMINCON solver is used to calculate the multi-energy flow in the MATLAB platform.
The coupling units between the power grid and the heating network mainly include the CHP units of each coupling node and the coupling unit of the power grid, the gas network includes the CHP of each coupling node and the gas turbine, and the coupling unit between the heat network and the gas network mainly includes the CHP of each coupling node. The connecting    Table 2.
In the model of gas turbine, the model parameters of a G c G 0, b G 167.1m 3 /MWh. In the model of CHP, this article considers two models: Model I is the pumping type, and Model II is the backpressure type. The parameters of CHP units are shown in Table 3.

Results and Analysis
In this study, we use two models to calculate the mass flow of each pipeline, water supply temperature, and return temperature of each node in the heating network. As is mentioned above, we define Eq. 43 as Model 1 and Eq. 42 as Model 2. Besides, we define Model 3, the control group, as the original model without linearization. The results of the heating network are shown in Tables 4-6.
In the process of solving, we find that the solving time of mode 1 and mode 2 is millisecond, while that of mode 3 is 21 s. In terms of solving time, the two models have the advantages of small amount of calculations and fast calculation speed.
From the above data, we can see that the maximum error of mass flow in Models 1 and 2 is 0.0058%, the maximum error of water supply temperature is 0.0075%, and return water is the same as that of Model 3. Both Models 1 and 2 have a high precision of solution. While it comes to each pipe's mass flow, Model 2 is more precise because the constant C has influence. Therefore, we use the results of Model 2 for the following calculation.     This study considers two models of CHP units and two linearized methods of the gas network. We define the onedimensional approximation method in the gas system as Method ① and the Taylor expansion approximation method as Method ②. So there are four running combinations in the following calculation, just as Table 7 shows. The gas network and electric system results are shown in Tables 8-12, including pressure, mass flow, voltage amplitude, and voltage phase angle.
By comparing the results of combinations 1 and 3, and combinations 2 and 4, we can see that no matter which mode the CHP unit runs in, there is a great difference between the results of Methods ① and ②, although the solving speed of both methods is in millisecond. The decisive factors of the accuracy of the two methods are different. Method ① depends on the number of segments and the value of constant K, while Method ② is only related to the initial value.
For Method ①, if we want to pursue accuracy, we need to divide more segments to achieve a higher degree of approximation. Still, at the same time, it also increases more unknowns and the amount of calculation, and the operation time becomes longer. Although Method ② is only suitable for networks with short pipeline length, it has a small amount of calculation. It involves fewer unknowns, while the natural gas network pipeline used in this study is shorter, so Method ② is more suitable for this study.
By comparing the results of combinations 1 and 2 and combinations 3 and 4, it can be seen that whichever method solving the gas system is adopted, power generated by the CHP unit is different between the two models. Model II generates more power than Model I, and at the same time, gas turbines generate less power in Model I than in Model II. Besides, when the operation mode of the CHP unit changes, the electric power generated by the CHP unit and the natural gas consumed change, which leads to the change of the power flow distribution of the power grid and the gas network, while in Model I, the electric power generated by the CHP unit depends on the thermal power generated by the CHP unit in the heating network. If the thermal power changes, the power flow distribution of the gas network, and the power grid will also be affected by it.
It can be seen that in the electric/thermal/gas integrated energy system, each network is connected into an inseparable whole through coupling units such as CHP and gas turbine, and the change of network state will affect the trend of other networks through the coupling unit.

CONCLUSION
To calculate the integrated energy system combined with electricity, gas, and heat more quickly and accurately, a comprehensive multi-energy flow calculation model based on convexification is proposed in this study. The model is established     according to the different characteristics of electricity, heat, and gas networks. The natural gas network pipeline model is linearized reasonably, which greatly reduces the complexity of the model. In the heat network, two kinds of radiant heat network models that can be solved quickly are established, and the original model is transformed into the problem of solving univariate first-order equation and univariate quadratic equation, respectively; the amount of calculation of which is small, and the calculation speed is fast without any convergence problem. Considering the backpressure type and the pumping type of CHP, four running combinations combined with two models of CHP and two methods of the gas system have been established. Finally, the simulation results show that the algorithm can complete the convergence quickly, proving the algorithm's rapidity and practicability. In addition, the algorithm used in this study takes into account the interaction between different networks, which can quickly get the distribution of the comprehensive power flow of the system.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.