# Integration and Modeling of Multi-Energy Network Based on Energy Hub

Authors: Min Mou; Yuhao Zhou; Wenguang Zheng; Yurong Xie
Journal: Complexity
 (2022)

---

## Abstract

The energy conversion units and energy storage equipment connected to the multi-energy system are becoming diversified, and the uncertain factors brought by distributed wind power and photovoltaic power generation make the system energy flow structure more complex, which brings great difficulties to the modeling and application of traditional energy hub modeling methods. This study deeply analyzes the multi-energy flow coupling structure and operation mechanism of multi-energy systems, and carries out the power flow calculation and analysis of multi-energy systems based on an energy hub, so as to ensure the safe and stable operation of regional energy. Based on the physical characteristics of energy systems such as power systems, thermal systems, and gas systems, this article studies the comprehensive power flow model including the electric-gas-thermal multi-energy coupling network and proposes the power flow decomposition of the energy supply subsystem and its applicable equation based on Newton–Raphson method. The effectiveness of the proposed method under different operation modes is verified by case studies. The calculation results show that under constant load, the energy hub running in fixing thermal by electricity (FEL) and fixing electricity by thermal (FTL) mode has little influence on the voltage of each node in the power sub-network. Within the constraint range, the natural gas flow obtained from the natural gas subsystem is coupled with the power subsystem to meet the load demand. The influence on the power flow at each node of the heat network is not obvious.


---

## Body

## 1. Introduction

With the popularization and development of multi-energy and low-carbon technologies, such as combined heat and power generation (CHP), heat pump, electric boiler and absorption refrigerator, the interaction between power network, heat pipe network and natural gas network in various links is closer, and the coordination between networks is stronger. Therefore, a new energy system with multiple types of energy such as natural gas, comes into being. Its main feature is multi-energy complementary coordination [1], unified planning, and unified scheduling of various energy sources such as electricity, gas, and heat. At the same time, based on the idea of interconnection, many energy supply systems can realize the dual interaction of energy and information for different individuals and regions like the internet [2]. This coordinated interconnection can not only improve the flexibility, security and economy of the system but also change the energy supply mode. The complementary advantages and synergistic benefits brought by this change can make up for the shortcomings of the separate operation of the traditional sub-energy system, which is of great significance to improve the energy utilization efficiency of the whole system. Therefore, it is highly praised by various countries [3].In the integrated energy system, the energy is no longer a relatively independent relationship, but a coupling relationship. Obviously, the traditional energy network analysis methods can not reflect the real nature of the coupling network. Therefore, it is necessary to study the modeling of multi-energy systems to tap the energy supply advantages and potential of different energy systems [4]. The concept of power flow is not only limited to the AC/DC hybrid power flow algorithm of micro-grid [5], but also refers to a multi-energy flow calculation method that satisfies the analysis of coupled energy demand under the framework of an integrated energy system. It can effectively evaluate the network operation state and is of great significance to realize safe and stable operation.In the power flow analysis and optimization research of electric-gas combined energy systems, a study [3] proposes the optimal power flow calculation method of the electric-gas coupled energy system, in which the dual interior point method is used to calculate the optimal power flow. Although this method has the advantage of fast calculation speed, it is not suitable for solving non-convex optimization problems. In the power flow analysis and optimization research of electric-thermal combined energy systems, a study [6] proposes the power flow decomposition and comprehensive solution method for electric-thermal coupling systems. On this basis, one study [7] establishes the power flow model of the CHP energy supply system and analyzes its economic and social benefits.In the power flow analysis and optimization research of the system considering the interaction and combination of multiple energy sources, the coupling relationship between electricity, gas, and heat is considered in detail in [8], and the corresponding power flow solution method is proposed. The study [9] uses the power flow decoupling method to transform the integrated energy system into an independent solution problem similar to the traditional form and puts forward a calculation method suitable for large-scale multi-energy power flow distribution. The study [10] comprehensively considers the actual operation constraints of each energy sub-network in the multi-energy system and realizes the multi-objective optimization calculation of comprehensive economy and emission. The study [11] coordinates the interests of all parties in the integrated energy system, constructs an optimization model with energy efficiency, economy, and environmental protection as indicators, and sets the Pareto frontier to obtain the optimal solution. Load demand will affect the operation of the integrated energy system with climate change. Therefore, the study [12] establishes the supply and demand power balance model of regionally integrated energy system under different scenarios according to different seasons and climate, selects the economy and environmental protection of the objective function, and considers the multiple operating conditions of the equipment, so as to solve the global optimal energy flow. Overall, there is still a lack of intuitive, easy, and general modeling methods for regional integrated energy networks.This study proposes a multi-energy network integration and modeling method based on energy hub, which is not only the basis for studying the scheduling and planning of an integrated energy system, but also the basis of ensuring the reliability and safe operation of the multi-energy coupling system. Through power flow calculation, the value of each state quantity in the system can be obtained. The operation state of the equipment and various equipment parameters in the system can be adjusted during operation, which provides help for the dispatching organization to formulate a reasonable energy supply scheme and ensure the safe and optimal operation of the system.


## 2. Power Flow Modeling of Multi-Energy System

Various forms of energy can be divided into electricity, gas, heat, etc. they are connected to the multi-energy system as energy supply subsystems. In addition, there are various energy hub models and coupling equipment to realize the interaction and transformation between a variety of energy sources. Therefore, when solving the power flow state of the multi-energy system, it is necessary to consider that this interaction and transformation will make the sub-power flow of a single system affected by other systems. It can also be understood that the power flow equation of an energy subsystem contains the state variables of other subsystems. This study mainly studies the multi-energy system for electricity-gas-heat interaction [13]. Therefore, the comprehensive power flow model of the regional comprehensive energy system can be described as(1)Fexe,xh,xg,xeh=0Fhxe,xh,xg,xeh=0Fgxe,xh,xg,xeh=0Fehxe,xh,xg,xeh=0,where Fe, Fh, Fg, and Feh, respectively, represent the power equations of the electric, thermal, and gas energy supply subsystems and their energy hubs in the multi-energy system; xe represents the variables related to the operation of the power subsystem [14]; xh represents the variables related to the operation of thermal subsystem; xg represents variables related to the operation of natural gas subsystem; xeh represents the variables related to the operation of the coupling elements in the energy hub.Due to the complex coupling relationship between each part of the energy network, the complexity of the coupled system is higher than that of the single energy system, and the solution method becomes more complex. For the coupling relationship, the electric-thermal-gas interaction part can be decoupled, and then the iterative solution results of each energy network can be coupled through the energy hub. When the iteration termination conditions are met, the power flow calculation results of each energy network can be output to realize the power flow solution of the multi-energy system. For this decomposition solution method, there is a more mature method, which is called the alternating decomposition method [15].


## 3. Power Flow Solution of Energy Supply Subsystem

When solving the multi-energy power flow distribution, the extended Newton–Raphson method is usually used as the iterative method for power flow solution [16]. The deviation of the power flow equation in the power flow iterative process of multi-energy system can be expressed by (2), and the iterative correction equation is shown in (3):(2)ΔPΔQΔΦΔhΔTsΔTrΔfT,(3)Δxk=−Jk−1Fxkxk+1=xk+Δxkk=0,1,2,⋯n,where ΔP is the active power deviation in the iterative process; ΔQ is reactive power deviation; ΔΦ is the thermal power deviation of the thermal system; Δh is the head pressure deviation of the thermal system; ΔTs, ΔTr respectively represents the deviation of heating water temperature and return water temperature of thermal system; Δf is the natural gas flow deviation vector of the natural gas system. xk is the vector composed of state variables at the kth iteration; Δxk is the correction vector for the state variable at the kth iteration; Jk is the Jacobian correction matrix at the kth iteration; Fx is the power flow equation.


### 3.1. Power Subsystem

During the operation of multi-energy system, each node of the power subsystem can be numbered, as shown in Table1.Table 1
Node number of power subsystem.
Node typePQ nodePV nodeBalance nodeNode number1, 2, …,mm + 1, m + 2, …, n − 1nNumber of nodesmn − m − 11The state variables in the iterative correction equation can be expressed by (4)–(6):(4)xEk=UkθkT,(5)Uk=ΔU1kU1kΔU2kU2k⋯ΔUmkUmkT,(6)θk=θ1kθ2k⋯θn−1kT,where Uk and θk respectively represents the substitution vector of the voltage amplitude and the voltage phase at the kth iteration.The Jacobian matrix and its elements in the corresponding iterative correction process are shown in (7)–(11):(7)JEk=HkNkMkLk,(8)Hij=∂ΔPi∂θj=−UiUjGijsinθij−Bijcosθij,i≠jUi2Bii+Qi,i=j,(9)Nij=∂ΔPi∂UjUj=−UiUjGijcosθij+Bijsinθij,i≠j−Ui2Gii−Pi,i=j,(10)Mij=∂ΔQi∂θj=UiUjGijcosθij+Bijsinθij,i≠jUi2Gii−Pi,i=j,(11)Hij=∂ΔQi∂UjUj=−UiUjGijsinθij−Bijcosθij,i≠jUi2Bii−Qi,i=j.


### 3.2. Natural Gas Subsystem

When calculating the power flow of natural gas subsystem, the natural gas system can be compared with the power system, and the nodes can be numbered, as shown in Table2.Table 2
Node number of natural gas subsystem.
Node typeUnbalanced nodeBalance nodeNode number1, 2, …,mm+ 1, m+ 2, …, nNumber of nodesmn −mWhen solving the power flow of natural gas subsystem, Newton–Raphson completes the power flow iteration process. According to the network modeling of natural gas subsystem, the node branch incidence matrix meets the node(12)Ag,ij=1flows into nodeifrom pipelinej−1flows out of nodeifrom pipelinej0no natural gas connection between pipelinejand nodei.The Weymouth equation can be replaced by (13), where πi=pi2, πj=pj2. And the branch pressure drop can be expressed by the node branch matrix as (14).(13)fij=Kijsijπi−πj=KijsijΔπij,(14)Δπ=−AgTπ.Combining (12)–(14), the unbalance deviation when using the Newton–Raphson method to solve the natural gas network power flow can be expressed as (15), the correction equation can be expressed as (16), and the elements of the Jacobian matrix are shown in (17).(15)Mπ=σ1π⋯σnπ=AgKs−AgTπ−f,(16)φπi=−Ji−1Mπiπi+1=πi+φπi,(17)Jii=12∑j∈ikijsijΔπij−12,i=mJij=−12kijsijΔπij−12Nodeiis connected to nodejJij=0,Nodeiand nodejare not connected.To sum up, after analogy with the power flow of the power subsystem, the power flow calculation process of the natural gas subsystem shall be solved uniformly and iteratively according to the following steps:(1)
Input the original data of the natural gas system, such as network topology and node parameters; classify and number the nodes;(2)
The initial pressure value of each non-equilibrium node is given, and the number of iterations is zero;(3)
After each iteration, calculate the deviation and check each component of the unbalance; if the maximum value of each component after iteration is less than the given accuracy, stop the iteration; the pressure of the non-equilibrium node, the natural gas flow of the equilibrium node and the flow of each branch pipeline are calculated.(4)
Calculate the Jacobian matrix and modify the state variables;(5)
Add one to the number of iterations and return to step (3).


### 3.3. Thermal Network Subsystem

Similar to the natural gas system, the thermal subsystem can also be classified and numbered by a method similar to the power subsystem, as shown in Table3.Table 3
Node number of the thermal subsystem.
Node typeLoad nodeHeat source nodeBalance nodeNode number1, 2, …,mm + 1, m + 2, …, n − 1nNumber of nodesmn − m − 11The hydraulic model and thermal model can be obtained from the thermal system model. Combined with these two models, they can be transformed into the general power flow model of the thermal subsystem [17], as shown in (18):(18)F1ml,Tiin,Tiout=ClAsmlTiin−Tiout−qiF2ml=Bhsml2F3ml,Tiin=CsTiin−asF4ml,Tiout=CrTiout−ar.When the Newton–Raphson method is used for unified iterative calculation of power flow in a thermal system, the selected state variable is (19), the correction equation is (20), and the Jacobian matrix and its elements are (21)–(24). The subscript of each matrix is the number of rows and columns of the matrix, mmix is the number of nodes with pipeline crossing in the load node, mnon−mix is the node without pipeline crossing in the heat supply network load node.(19)x=mlnpipeTiinn−m−1Tioutn−m−1,(20)Δxi=−Ji−1Fxixi+1=xi+Δxi,(21)J=J11J12J13J21J22J23J31J32J33,(22)J11=∂F1;F2∂ml=2Bhsmlnloop×npipeClAsTiin−Toin−1×npipeJ12=∂F1;F2∂Tiin=0nloop×ndiagClAsmln−1×n−1J13=∂F1;F2∂Tiout=0nloop×ndiag−ClAsmln−1×n−10mEn−m−1,(23)J21=∂F3∂ml=0m×mnon−mix0mnon−mix×n−mnon−mix−Ti,jinψk+Ti,iin&0mmix×n−mnon−mixJ22=∂F3∂Tiin=Csm×mJ23=∂F3∂Tiout=0m×m,(24)J31=∂F4∂ml=Ti,iout−Toi&−Ti,ioutψk+Toi&00mnon−mix×m0m×n−mJ32=∂F4∂Tiin=0m×mJ33=∂F4∂Tiout=Crm×m.Similar to the power flow calculation of the natural gas subsystem, the specific steps in the power flow calculation of thermal subsystem can be summarized as follows:(1)
Input the original data of the thermal system, including network topology and node parameters; classify and number the nodes, and select the balance node of the heating network [18, 19];(2)
The initial values of heating temperature and return water temperature of each pipeline flow and load node are given, and the number of iterations is zero;(3)
Calculate the unbalance after each iteration and check each component of the unbalance; if the maximum value of each component after iteration is less than the given accuracy, stop the iteration; calculate the balance node power and return water temperature;(4)
Calculate the Jacobian matrix and modify the state variables, and calculate the return water temperature of each node;(5)
Add one to the number of iterations and return to step (3).


## 3.1. Power Subsystem

During the operation of multi-energy system, each node of the power subsystem can be numbered, as shown in Table1.Table 1
Node number of power subsystem.
Node typePQ nodePV nodeBalance nodeNode number1, 2, …,mm + 1, m + 2, …, n − 1nNumber of nodesmn − m − 11The state variables in the iterative correction equation can be expressed by (4)–(6):(4)xEk=UkθkT,(5)Uk=ΔU1kU1kΔU2kU2k⋯ΔUmkUmkT,(6)θk=θ1kθ2k⋯θn−1kT,where Uk and θk respectively represents the substitution vector of the voltage amplitude and the voltage phase at the kth iteration.The Jacobian matrix and its elements in the corresponding iterative correction process are shown in (7)–(11):(7)JEk=HkNkMkLk,(8)Hij=∂ΔPi∂θj=−UiUjGijsinθij−Bijcosθij,i≠jUi2Bii+Qi,i=j,(9)Nij=∂ΔPi∂UjUj=−UiUjGijcosθij+Bijsinθij,i≠j−Ui2Gii−Pi,i=j,(10)Mij=∂ΔQi∂θj=UiUjGijcosθij+Bijsinθij,i≠jUi2Gii−Pi,i=j,(11)Hij=∂ΔQi∂UjUj=−UiUjGijsinθij−Bijcosθij,i≠jUi2Bii−Qi,i=j.


## 3.2. Natural Gas Subsystem

When calculating the power flow of natural gas subsystem, the natural gas system can be compared with the power system, and the nodes can be numbered, as shown in Table2.Table 2
Node number of natural gas subsystem.
Node typeUnbalanced nodeBalance nodeNode number1, 2, …,mm+ 1, m+ 2, …, nNumber of nodesmn −mWhen solving the power flow of natural gas subsystem, Newton–Raphson completes the power flow iteration process. According to the network modeling of natural gas subsystem, the node branch incidence matrix meets the node(12)Ag,ij=1flows into nodeifrom pipelinej−1flows out of nodeifrom pipelinej0no natural gas connection between pipelinejand nodei.The Weymouth equation can be replaced by (13), where πi=pi2, πj=pj2. And the branch pressure drop can be expressed by the node branch matrix as (14).(13)fij=Kijsijπi−πj=KijsijΔπij,(14)Δπ=−AgTπ.Combining (12)–(14), the unbalance deviation when using the Newton–Raphson method to solve the natural gas network power flow can be expressed as (15), the correction equation can be expressed as (16), and the elements of the Jacobian matrix are shown in (17).(15)Mπ=σ1π⋯σnπ=AgKs−AgTπ−f,(16)φπi=−Ji−1Mπiπi+1=πi+φπi,(17)Jii=12∑j∈ikijsijΔπij−12,i=mJij=−12kijsijΔπij−12Nodeiis connected to nodejJij=0,Nodeiand nodejare not connected.To sum up, after analogy with the power flow of the power subsystem, the power flow calculation process of the natural gas subsystem shall be solved uniformly and iteratively according to the following steps:(1)
Input the original data of the natural gas system, such as network topology and node parameters; classify and number the nodes;(2)
The initial pressure value of each non-equilibrium node is given, and the number of iterations is zero;(3)
After each iteration, calculate the deviation and check each component of the unbalance; if the maximum value of each component after iteration is less than the given accuracy, stop the iteration; the pressure of the non-equilibrium node, the natural gas flow of the equilibrium node and the flow of each branch pipeline are calculated.(4)
Calculate the Jacobian matrix and modify the state variables;(5)
Add one to the number of iterations and return to step (3).


## 3.3. Thermal Network Subsystem

Similar to the natural gas system, the thermal subsystem can also be classified and numbered by a method similar to the power subsystem, as shown in Table3.Table 3
Node number of the thermal subsystem.
Node typeLoad nodeHeat source nodeBalance nodeNode number1, 2, …,mm + 1, m + 2, …, n − 1nNumber of nodesmn − m − 11The hydraulic model and thermal model can be obtained from the thermal system model. Combined with these two models, they can be transformed into the general power flow model of the thermal subsystem [17], as shown in (18):(18)F1ml,Tiin,Tiout=ClAsmlTiin−Tiout−qiF2ml=Bhsml2F3ml,Tiin=CsTiin−asF4ml,Tiout=CrTiout−ar.When the Newton–Raphson method is used for unified iterative calculation of power flow in a thermal system, the selected state variable is (19), the correction equation is (20), and the Jacobian matrix and its elements are (21)–(24). The subscript of each matrix is the number of rows and columns of the matrix, mmix is the number of nodes with pipeline crossing in the load node, mnon−mix is the node without pipeline crossing in the heat supply network load node.(19)x=mlnpipeTiinn−m−1Tioutn−m−1,(20)Δxi=−Ji−1Fxixi+1=xi+Δxi,(21)J=J11J12J13J21J22J23J31J32J33,(22)J11=∂F1;F2∂ml=2Bhsmlnloop×npipeClAsTiin−Toin−1×npipeJ12=∂F1;F2∂Tiin=0nloop×ndiagClAsmln−1×n−1J13=∂F1;F2∂Tiout=0nloop×ndiag−ClAsmln−1×n−10mEn−m−1,(23)J21=∂F3∂ml=0m×mnon−mix0mnon−mix×n−mnon−mix−Ti,jinψk+Ti,iin&0mmix×n−mnon−mixJ22=∂F3∂Tiin=Csm×mJ23=∂F3∂Tiout=0m×m,(24)J31=∂F4∂ml=Ti,iout−Toi&−Ti,ioutψk+Toi&00mnon−mix×m0m×n−mJ32=∂F4∂Tiin=0m×mJ33=∂F4∂Tiout=Crm×m.Similar to the power flow calculation of the natural gas subsystem, the specific steps in the power flow calculation of thermal subsystem can be summarized as follows:(1)
Input the original data of the thermal system, including network topology and node parameters; classify and number the nodes, and select the balance node of the heating network [18, 19];(2)
The initial values of heating temperature and return water temperature of each pipeline flow and load node are given, and the number of iterations is zero;(3)
Calculate the unbalance after each iteration and check each component of the unbalance; if the maximum value of each component after iteration is less than the given accuracy, stop the iteration; calculate the balance node power and return water temperature;(4)
Calculate the Jacobian matrix and modify the state variables, and calculate the return water temperature of each node;(5)
Add one to the number of iterations and return to step (3).


## 4. Multi-Energy Network Power Flow Calculation

### 4.1. Power Calculation of Coupling Part Energy Hub

Multi-energy systems generally work in the grid connection mode, that is, they are directly connected with the power grid and other energy equipment through the tie line and obtain the energy transmitted by them [20]. Under the grid connection mode, the active and reactive powers of each equipment have been given, and the frequency of the power subsystem is controlled by the large power grid [21]. In the process of power flow iteration, the calculation methods after coupling are different due to different operation modes. When the electric heating load is known, the interaction value between the energy hub and each network can be obtained. When the PQ node of the power system is taken as the balance node, its voltage and parameter values can be obtained, as shown in (25) and (26); Thus, the pressure and flow value of the equipment in the natural gas network coupled with the electric and thermal network can be obtained, as shown in (27).(25)Pe,iMT=vMTηge,iMTPg,iPg,iMT=vMTPg,i=Pe,iMT/ηge,iMT,(26)Pe,iAC=vACPe,iPh,iMT=ηgh,iMT∗Pe,iMT/ηge,iMTPh,iAC=vACPe,iηiAC=Pe,iACPg,iGB=1−vMTηGB=Lh,i−Ph,iMT−Ph,iAC/ηGB,(27)fMT=Pg,iMT/LHVfGB=Pg,iGB/LHV,where LHV is the low calorific value of fuel, also known as net thermal efficiency. Similarly, when the natural gas balance node is selected to be solved first, the electrical network parameters of the coupling equipment can be obtained by using (25) to (27) through the known quantity of the coupling matrix of the joint energy hub.


### 4.2. Power Flow Decomposition Algorithm for Multi-Energy Network

Combined with the modeling and power flow calculation methods in the above sections, this section introduces the power flow decomposition algorithm of multi-energy system. The multi-energy system power flow decomposition algorithm takes into account the interactive coupling of electricity, heat, gas, and other energy sources, realizes power flow calculation by decoupling them, and can adjust the correlation value between the hub and the energy network, so as to obtain a more accurate multi-energy flow distribution based on the interaction between electricity, heat, and gas [22]. Figure 1 shows the method of solving multi-energy system power flow and the specific steps can be expressed as follows:(1)
Input the network parameters of each energy sub-network, node type, node injection energy, and the electric heating load demand of the system;(2)
Select the energy hub model and its operating parameters, generate the conversion matrix, and convert it into the corresponding network node parameters;(3)
According to the known quantity of electric heating load, the electrical demand is converted by using the conversion matrix;(4)
The power flow distribution is solved independently by using the energy supply subsystem model and its balance node parameters [23]; if the result converges, go to step (5); if not, correct the power exchange between the energy hub and the network and return to step (3);(5)
According to the power flow calculation results of each subsystem, the joint coupling equipment model pushes back and forth to calculate the parameters of the coupling network; if the constraints are met, go to step (6); if the results do not meet the constraints, adjust the parameter setting of the network coupling node and return to step (2);(6)
Output the result and end the operation.Figure 1
Flow chart of power flow decomposition algorithm.


## 4.1. Power Calculation of Coupling Part Energy Hub

Multi-energy systems generally work in the grid connection mode, that is, they are directly connected with the power grid and other energy equipment through the tie line and obtain the energy transmitted by them [20]. Under the grid connection mode, the active and reactive powers of each equipment have been given, and the frequency of the power subsystem is controlled by the large power grid [21]. In the process of power flow iteration, the calculation methods after coupling are different due to different operation modes. When the electric heating load is known, the interaction value between the energy hub and each network can be obtained. When the PQ node of the power system is taken as the balance node, its voltage and parameter values can be obtained, as shown in (25) and (26); Thus, the pressure and flow value of the equipment in the natural gas network coupled with the electric and thermal network can be obtained, as shown in (27).(25)Pe,iMT=vMTηge,iMTPg,iPg,iMT=vMTPg,i=Pe,iMT/ηge,iMT,(26)Pe,iAC=vACPe,iPh,iMT=ηgh,iMT∗Pe,iMT/ηge,iMTPh,iAC=vACPe,iηiAC=Pe,iACPg,iGB=1−vMTηGB=Lh,i−Ph,iMT−Ph,iAC/ηGB,(27)fMT=Pg,iMT/LHVfGB=Pg,iGB/LHV,where LHV is the low calorific value of fuel, also known as net thermal efficiency. Similarly, when the natural gas balance node is selected to be solved first, the electrical network parameters of the coupling equipment can be obtained by using (25) to (27) through the known quantity of the coupling matrix of the joint energy hub.


## 4.2. Power Flow Decomposition Algorithm for Multi-Energy Network

Combined with the modeling and power flow calculation methods in the above sections, this section introduces the power flow decomposition algorithm of multi-energy system. The multi-energy system power flow decomposition algorithm takes into account the interactive coupling of electricity, heat, gas, and other energy sources, realizes power flow calculation by decoupling them, and can adjust the correlation value between the hub and the energy network, so as to obtain a more accurate multi-energy flow distribution based on the interaction between electricity, heat, and gas [22]. Figure 1 shows the method of solving multi-energy system power flow and the specific steps can be expressed as follows:(1)
Input the network parameters of each energy sub-network, node type, node injection energy, and the electric heating load demand of the system;(2)
Select the energy hub model and its operating parameters, generate the conversion matrix, and convert it into the corresponding network node parameters;(3)
According to the known quantity of electric heating load, the electrical demand is converted by using the conversion matrix;(4)
The power flow distribution is solved independently by using the energy supply subsystem model and its balance node parameters [23]; if the result converges, go to step (5); if not, correct the power exchange between the energy hub and the network and return to step (3);(5)
According to the power flow calculation results of each subsystem, the joint coupling equipment model pushes back and forth to calculate the parameters of the coupling network; if the constraints are met, go to step (6); if the results do not meet the constraints, adjust the parameter setting of the network coupling node and return to step (2);(6)
Output the result and end the operation.Figure 1
Flow chart of power flow decomposition algorithm.


## 5. Case Study

### 5.1. Basic Data

IEEE14 node, 9-node natural gas system, and 12 node thermal system are used as examples. In the coupling equipment part, CHP unit, gas boiler GB and central air conditioning AC are selected as the energy hub for coupling each energy system, and the grid-connected working mode is used for the system. The schematic diagram of the calculation example is shown in Figure2. In the figure, EB, GB, and Hb are used to represent the electrical system bus, natural gas system bus and thermal bus, and the numbers represent the number of nodes. The following describes the configuration of network parameters and node configuration.(1)
In IEEE 14 node power system, node parameter setting is shown in Table4 and node classification is shown in Table 5.Among them, node 1 is the public connection point between the multi-energy system and the large power grid; Other PV nodes 2, 6, and 8 can be used as access points for new energy power generation equipment or coal-fired units and other power generation units; Node 3 is the connection point between the power system and the gas turbine to realize the coupling interaction with node 8 of the natural gas subsystem; Node 14 serves as the power load node connected to the air conditioning AC to realize the coupling with the heating network and the gas network; The remaining nodes except node 7 are regarded as load nodes.(2)
In the 9-node natural gas system, the node classification is shown in Table6, and the network node parameter setting of the natural gas system is shown in Table 7.Among them, GB1 is the constant pressure node of the system, GB8 and GB9 are constant pressure and constant current nodes, and other gas distribution nodes can be simplified to be expressed only by the flow. At the same time, as an output node of the natural gas network, GB8 is connected with the gas turbine (MT) to generate electric energy, which realizes the coupling and interaction between the two electrical energy sources. The gas turbine generates heat, which realizes the coupling and interaction with the thermal system; as the access point of the gas-fired boiler (GB), GB9 realizes the coupling of the natural gas system and thermal network.(3)
In the 12-node thermal system, the load of each node is shown in Table8, and the pipeline parameters and heating distance are shown in Table 9.Figure 2
Schematic diagram of calculation example of the multi-energy system.Table 4
Power system node parameters.
Node numberGenerating power (MW)Load power (MW)Initial value of voltage (V)Voltage setting value at PV node (V)160011.06265 + 42.4i21.7 + 12.7i11.0453Pmt94.2 + 19i11.014047.8 − 3.9i10507.6 − 1.6i10685 + 12.24i11.2 + 7.5i11.0770010817.36i011.099029.5 + 16.6i101009 + 5.8i101103.5 + 1.8i101206.1 + 1.6i1013013.5 + 5.8i10140Pac10Table 5
Power system node classification.
Node classificationNode numberBalance node1PV node2, 3, 6, 8PQ node4, 5, 9, 10, 11, 12, 13, 14Table 6
Node classification of natural gas subsystem.
Node numberNode classification1Gas source node3, 4, 5, 8, 9Gas distribution nodeTable 7
Node parameter values of the natural gas subsystem.
Node numberSet value1π1=5MPa2f2=03f3=200m3/h4f4=150m3/h5f5=250m3/h8π8=5.5MPa,f8=fMT9π9=6MPa,f9=fGBTable 8
Load of each node of the thermal subsystem.
Node numberHeat load powerOutput temperatureNode numberHeat load powerOutput temperature10.154070.14020.154080.14030.154090.254040.1540100.154050.1540110.24060.154012——Table 9
Thermal pipe network parameters.
Pipeline directionLengthDiameter12-15002001-24002002-36002004-34002001-52002001-61502002–71802002–81502003–91002003–101102004–1190200.In the thermal system, only node 12 is the heat source node, and the other nodes are the heating load nodes. The heat source nodes are connected with the energy hub, and are associated and coupled with the power grid and gas network through the air conditioning AC, gas turbine MT and gas boiler GB in the energy hub.


### 5.2. Energy Hub Parameters

The load parameters in the energy hub are taken respectively, and the electrical loadLe is 180 kW, the thermal load Lh is 300 kW and other loads Lo are 100 kW; the performance and parameters of the energy hub of the coupling part, such as gas turbine, transformer, distributed energy, and central air conditioning, are shown in Table 10.Table 10
Equipment parameters in heat supply network.
Parameter variableTypeValueMicro gas turbine (MT)Upper power limit150Generating power0.35Thermal power0.43Gas-fired boiler (GB)Upper power limit500Efficiency0.81Central air-conditioning (AC)Upper power limit200Efficiency0.62Transformer (T)Efficiency Nt0.95Efficiency Noe0.1Distributed energyEfficiency Nog0.1Efficiency No0.98The energy distribution coefficient shown in (28) can be coupled with other energy conversion parameters of the distributed device, and the energy distribution coefficient can be reduced.(28)180300100=0.95×1−vAC−voe0.35×vMT0.10.62×vAC0.43×vMT+0.81×1−vMT−vog0.1voevog0.98PePgPo,where the distribution coefficient voe and vog are almost zero. When the power flow is calculated first, the fMT and fGB can be solved by PeMT and PeAC, so as to solve the coupled natural gas network power flow; when the natural gas power flow is calculated first, the initial value can be set first, and then the corresponding electrical interaction value can be calculated, so as to solve the coupled power network power flow. After each interaction, the value of the coupling network is obtained, and the constraints of the integrated power flow network of the multi-energy system can be verified through back calculation, so as to adjust the distribution coefficient and stabilize the power flow distribution.


### 5.3. Result Analysis

The impact of different operation modes on the system power flow is analyzed to determine the power flow distribution of multi-energy systems. For more intuitive analysis, it is assumed that in the whole multi-energy system, only the input and output of the coupling node and balance node change, and the electric, gas, and heat energy input and load demand of other nodes are constant. The two operation modes of the energy hub are fixing thermal by electricity (FEL) and fixing electricity by thermal (FTL). In the electric heating mode, the natural gas can be obtained from the electrical quantity, and in the thermal heating mode, the electrical quantity can be obtained from the natural gas. By comparing the two modes, the initial solution of the power flow calculation of the system can be obtained, as shown in Table11.Table 11
Initial solution of power flow in multi-energy systems.
Operation modePower subsystemNatural gas subsystemPmtPacfmtfgbFEL9050⟶27.6518.51FTL97.6543.37←3020


#### 5.3.1. Power Subsystem

From the above data, we can get the voltage value of each node in the power subsystem under FEL and FTL modes. The results are shown in Figure3.Figure 3
Voltage of each node in different modes.Through the voltage results of each node shown in Figure3, it can be found that the voltage amplitude of each node of the power subsystem in FEL mode is slightly higher than that in FTL mode, but the difference is extremely small and remains within the stable range. This shows that under the condition of constant load, the operation of the energy hub in FEL and FTL modes has little effect on the voltage of each node of power sub-network.


#### 5.3.2. Natural Gas Subsystem

The pressure change at each node of the natural gas system under the two different operation modes is shown in Figure4.Figure 4
Pressure of each node in different modes.From the node pressure results of the natural gas subsystem under different operation modes shown in Figure4, it is not difficult to see that under the assumption that the load of the branch node does not change, the air pressure value of each node is stable within the rated working pressure. Here, the energy hub plays a regulating role, coupling and interconnecting the natural gas flow obtained in the natural gas subsystem with the power subsystem within the constraint range, and jointly output to meet the load demand.


#### 5.3.3. Thermal Subsystem

The supply and return water temperature of each node of the thermal system under the above two operation modes are shown in Figures5 and 6.Figure 5
Water supply temperature at each node.Figure 6
Return water temperature of each node.It can be seen from Figures5 and 6 that when the energy hub operates in two different modes, the variation range of heating water temperature and return water temperature at each node of the heat supply network is small, i.e., the impact of power flow at each node of the heat supply network is not obvious, which is related to the selection of the scale of the heat supply network system, and also because the delay change of the heat supply network is not reflected in time. Node 12 is a heat source node, and its heating temperature is the highest. Because it provides heat to other nodes, its temperature change is the largest.


## 5.1. Basic Data

IEEE14 node, 9-node natural gas system, and 12 node thermal system are used as examples. In the coupling equipment part, CHP unit, gas boiler GB and central air conditioning AC are selected as the energy hub for coupling each energy system, and the grid-connected working mode is used for the system. The schematic diagram of the calculation example is shown in Figure2. In the figure, EB, GB, and Hb are used to represent the electrical system bus, natural gas system bus and thermal bus, and the numbers represent the number of nodes. The following describes the configuration of network parameters and node configuration.(1)
In IEEE 14 node power system, node parameter setting is shown in Table4 and node classification is shown in Table 5.Among them, node 1 is the public connection point between the multi-energy system and the large power grid; Other PV nodes 2, 6, and 8 can be used as access points for new energy power generation equipment or coal-fired units and other power generation units; Node 3 is the connection point between the power system and the gas turbine to realize the coupling interaction with node 8 of the natural gas subsystem; Node 14 serves as the power load node connected to the air conditioning AC to realize the coupling with the heating network and the gas network; The remaining nodes except node 7 are regarded as load nodes.(2)
In the 9-node natural gas system, the node classification is shown in Table6, and the network node parameter setting of the natural gas system is shown in Table 7.Among them, GB1 is the constant pressure node of the system, GB8 and GB9 are constant pressure and constant current nodes, and other gas distribution nodes can be simplified to be expressed only by the flow. At the same time, as an output node of the natural gas network, GB8 is connected with the gas turbine (MT) to generate electric energy, which realizes the coupling and interaction between the two electrical energy sources. The gas turbine generates heat, which realizes the coupling and interaction with the thermal system; as the access point of the gas-fired boiler (GB), GB9 realizes the coupling of the natural gas system and thermal network.(3)
In the 12-node thermal system, the load of each node is shown in Table8, and the pipeline parameters and heating distance are shown in Table 9.Figure 2
Schematic diagram of calculation example of the multi-energy system.Table 4
Power system node parameters.
Node numberGenerating power (MW)Load power (MW)Initial value of voltage (V)Voltage setting value at PV node (V)160011.06265 + 42.4i21.7 + 12.7i11.0453Pmt94.2 + 19i11.014047.8 − 3.9i10507.6 − 1.6i10685 + 12.24i11.2 + 7.5i11.0770010817.36i011.099029.5 + 16.6i101009 + 5.8i101103.5 + 1.8i101206.1 + 1.6i1013013.5 + 5.8i10140Pac10Table 5
Power system node classification.
Node classificationNode numberBalance node1PV node2, 3, 6, 8PQ node4, 5, 9, 10, 11, 12, 13, 14Table 6
Node classification of natural gas subsystem.
Node numberNode classification1Gas source node3, 4, 5, 8, 9Gas distribution nodeTable 7
Node parameter values of the natural gas subsystem.
Node numberSet value1π1=5MPa2f2=03f3=200m3/h4f4=150m3/h5f5=250m3/h8π8=5.5MPa,f8=fMT9π9=6MPa,f9=fGBTable 8
Load of each node of the thermal subsystem.
Node numberHeat load powerOutput temperatureNode numberHeat load powerOutput temperature10.154070.14020.154080.14030.154090.254040.1540100.154050.1540110.24060.154012——Table 9
Thermal pipe network parameters.
Pipeline directionLengthDiameter12-15002001-24002002-36002004-34002001-52002001-61502002–71802002–81502003–91002003–101102004–1190200.In the thermal system, only node 12 is the heat source node, and the other nodes are the heating load nodes. The heat source nodes are connected with the energy hub, and are associated and coupled with the power grid and gas network through the air conditioning AC, gas turbine MT and gas boiler GB in the energy hub.


## 5.2. Energy Hub Parameters

The load parameters in the energy hub are taken respectively, and the electrical loadLe is 180 kW, the thermal load Lh is 300 kW and other loads Lo are 100 kW; the performance and parameters of the energy hub of the coupling part, such as gas turbine, transformer, distributed energy, and central air conditioning, are shown in Table 10.Table 10
Equipment parameters in heat supply network.
Parameter variableTypeValueMicro gas turbine (MT)Upper power limit150Generating power0.35Thermal power0.43Gas-fired boiler (GB)Upper power limit500Efficiency0.81Central air-conditioning (AC)Upper power limit200Efficiency0.62Transformer (T)Efficiency Nt0.95Efficiency Noe0.1Distributed energyEfficiency Nog0.1Efficiency No0.98The energy distribution coefficient shown in (28) can be coupled with other energy conversion parameters of the distributed device, and the energy distribution coefficient can be reduced.(28)180300100=0.95×1−vAC−voe0.35×vMT0.10.62×vAC0.43×vMT+0.81×1−vMT−vog0.1voevog0.98PePgPo,where the distribution coefficient voe and vog are almost zero. When the power flow is calculated first, the fMT and fGB can be solved by PeMT and PeAC, so as to solve the coupled natural gas network power flow; when the natural gas power flow is calculated first, the initial value can be set first, and then the corresponding electrical interaction value can be calculated, so as to solve the coupled power network power flow. After each interaction, the value of the coupling network is obtained, and the constraints of the integrated power flow network of the multi-energy system can be verified through back calculation, so as to adjust the distribution coefficient and stabilize the power flow distribution.


## 5.3. Result Analysis

The impact of different operation modes on the system power flow is analyzed to determine the power flow distribution of multi-energy systems. For more intuitive analysis, it is assumed that in the whole multi-energy system, only the input and output of the coupling node and balance node change, and the electric, gas, and heat energy input and load demand of other nodes are constant. The two operation modes of the energy hub are fixing thermal by electricity (FEL) and fixing electricity by thermal (FTL). In the electric heating mode, the natural gas can be obtained from the electrical quantity, and in the thermal heating mode, the electrical quantity can be obtained from the natural gas. By comparing the two modes, the initial solution of the power flow calculation of the system can be obtained, as shown in Table11.Table 11
Initial solution of power flow in multi-energy systems.
Operation modePower subsystemNatural gas subsystemPmtPacfmtfgbFEL9050⟶27.6518.51FTL97.6543.37←3020


### 5.3.1. Power Subsystem

From the above data, we can get the voltage value of each node in the power subsystem under FEL and FTL modes. The results are shown in Figure3.Figure 3
Voltage of each node in different modes.Through the voltage results of each node shown in Figure3, it can be found that the voltage amplitude of each node of the power subsystem in FEL mode is slightly higher than that in FTL mode, but the difference is extremely small and remains within the stable range. This shows that under the condition of constant load, the operation of the energy hub in FEL and FTL modes has little effect on the voltage of each node of power sub-network.


### 5.3.2. Natural Gas Subsystem

The pressure change at each node of the natural gas system under the two different operation modes is shown in Figure4.Figure 4
Pressure of each node in different modes.From the node pressure results of the natural gas subsystem under different operation modes shown in Figure4, it is not difficult to see that under the assumption that the load of the branch node does not change, the air pressure value of each node is stable within the rated working pressure. Here, the energy hub plays a regulating role, coupling and interconnecting the natural gas flow obtained in the natural gas subsystem with the power subsystem within the constraint range, and jointly output to meet the load demand.


### 5.3.3. Thermal Subsystem

The supply and return water temperature of each node of the thermal system under the above two operation modes are shown in Figures5 and 6.Figure 5
Water supply temperature at each node.Figure 6
Return water temperature of each node.It can be seen from Figures5 and 6 that when the energy hub operates in two different modes, the variation range of heating water temperature and return water temperature at each node of the heat supply network is small, i.e., the impact of power flow at each node of the heat supply network is not obvious, which is related to the selection of the scale of the heat supply network system, and also because the delay change of the heat supply network is not reflected in time. Node 12 is a heat source node, and its heating temperature is the highest. Because it provides heat to other nodes, its temperature change is the largest.


## 5.3.1. Power Subsystem

From the above data, we can get the voltage value of each node in the power subsystem under FEL and FTL modes. The results are shown in Figure3.Figure 3
Voltage of each node in different modes.Through the voltage results of each node shown in Figure3, it can be found that the voltage amplitude of each node of the power subsystem in FEL mode is slightly higher than that in FTL mode, but the difference is extremely small and remains within the stable range. This shows that under the condition of constant load, the operation of the energy hub in FEL and FTL modes has little effect on the voltage of each node of power sub-network.


## 5.3.2. Natural Gas Subsystem

The pressure change at each node of the natural gas system under the two different operation modes is shown in Figure4.Figure 4
Pressure of each node in different modes.From the node pressure results of the natural gas subsystem under different operation modes shown in Figure4, it is not difficult to see that under the assumption that the load of the branch node does not change, the air pressure value of each node is stable within the rated working pressure. Here, the energy hub plays a regulating role, coupling and interconnecting the natural gas flow obtained in the natural gas subsystem with the power subsystem within the constraint range, and jointly output to meet the load demand.


## 5.3.3. Thermal Subsystem

The supply and return water temperature of each node of the thermal system under the above two operation modes are shown in Figures5 and 6.Figure 5
Water supply temperature at each node.Figure 6
Return water temperature of each node.It can be seen from Figures5 and 6 that when the energy hub operates in two different modes, the variation range of heating water temperature and return water temperature at each node of the heat supply network is small, i.e., the impact of power flow at each node of the heat supply network is not obvious, which is related to the selection of the scale of the heat supply network system, and also because the delay change of the heat supply network is not reflected in time. Node 12 is a heat source node, and its heating temperature is the highest. Because it provides heat to other nodes, its temperature change is the largest.


## 6. Conclusion

This study models and analyzes the integrated network power flow when the multi-energy system is operating and analyzes the energy conversion relationship of the energy hub model. The Newton–Raphson method is used to calculate the power flow of each energy supply subsystem. At the same time, for the electrical-thermal interaction coupling unit, the electric heating load is decoupled in the multi-energy network. When the power flow is solvable, the cross-back derivation algorithm is used to adjust the interaction value between the hub and the network, so as to obtain an accurate multi-energy power flow distribution. Finally, according to the comprehensive power flow calculation method of multi-energy systems, the simulation examples and analysis of the corresponding electrical-thermal interaction system are given, and the power flow calculation results under different operating modes are compared and analyzed.



---

*Source: 2698226-2022-09-05.xml*