A Novel Methodology for Electric-Thermal Mixed Power Flow Simulation and Transmission Loss Analysis in Multi-Energy Micro-Grids

For urban electrical and thermal energy supply and consumption terminals, the multi-energy system is a way to realize the energy conversion and consumption efficiently, cleanly and economically. To further study the power flow character of multi-energy systems, additional factors such as electric-thermal coupling equipment, thermal flow rules should be taken into consideration. The traditional Newton-Raphson iteration method which has been commonly used for electrical power flow calculation must be expanded to electric-thermal mixed power flow. This paper realized a new expression of hot water temperature drop with the thermal supply pipe network. This expression modified the Sukhov cooling operator and intuitively revealed the impact factor of transmission loss in hot water pipes. As a result, mathematical descriptions of electricity, heat and hydraulic flow in a multi-energy system were established. With the modified Sukhov cooling operator, the mathematical model of the expanded Newton-Raphson fast iterative solution algorithm, together with its Jacobian matrix elements affected by the electro-thermal coupling relationship was derived. A regional electric-heat supply system was selected as an example to verify the effectiveness of the method. Results showed that the transmission loss in electrical grid is related with the power load while in thermal pipes it was mainly related with the ambient temperature.


PREFACES
In recent years, various forms of energy coupling and complementary multi-energy systems have developed rapidly, and many researchers have begun to research on multi-energy operational and planning strategies and continue to reveal its cooperative coupling mechanisms (Omalley and Kroposki, 2013;Mancarella, 2014). In order to discover the high-efficiency, clean, and economical characteristics of multi-energy system, it is necessary to firstly establish the expression of its internal mechanisms, especially the mathematical model of multi-energy power flow, and simplify it by linearizing the nonlinear factors (Geidl and Andersson, 2007;Loukarakis and Mancarella, 2017). In terms of system modeling, one of the most representatives is the energy hub model proposed by Geidl et al. in 2007(Geidl, 2007, which details the modeling of energy conversion units and energy transmission units in multi-energy flow systems. Reference (Good, et al., 2015) introduces and details a high resolution domestic multi-energy model comprised of physically based space heating, domestic hot water (DHW), cooking and electrical appliance models. The presented model can be utilized by many actors involved in energy systems for various purposes. As a nonlinear, and strongly coupled complex system with multi-variables, the simplification and linearization of the mathematical model of the integrated energy system is essential for its operation simulation and planning configuration (Shao et al., 2017;Huang et al., 2020). In terms of system simulation and calculation of multi-energy power flow, the representative is Dr. Liu from Cardiff University who has made detailed research on simulation method of thermoelectric coupling system (Liu, 2013;Liu et al., 2016). Based on the model of energy hub, the difference between decoupling and coupling solution of thermoelectric system is dentally compared (Moeini-Aghtaie et al., 2014;Shabanpour-Haghighi and Seifi, 2015), and local optimization is strives to obtain the best power flow solution. Considering the thermal network characteristics could imply significant convergence difficulties for optimization solvers (Loukarakis and Mancarella, 2017). In (Lingmin et al., 2020), Particle Swarm Optimization (PSO) algorithm is used to achieve energy flow optimization; while in (Kriechbaum et al., 2020) cumulative exergy minimization together with load flow calculations is used to determine the optimal system design of a multi-cell municipal energy system. To deal with nonlinear heat transfer constraint, Reference (He et al., 2020) proposes a heat current model to comprehensively reflect the heat transfer, migration and storage characteristics. Based on the conclusions obtained from power flow calculation of multi-energy system, control strategy of multi-microgrids (Kargarian and Rahmani, 2015), and even operational guidance can be provided for absorb renewable energy (Salpakari et al., 2016). But fast and efficient solutions to multiple energy flows are the basis and prerequisite for any application.
Therefore, in view of the requirement of fast solution and high accuracy in multi-energy flow calculation, this paper simplifies the cooling operator in the iteration process, establishes the mathematical model of extended Newton-Raphson fast iteration solution algorithm, and deduces the calculation method of Jacobi matrix element value affected by electricthermal coupling relationship. An example of a regional electrical-thermal supply system is used to demonstrate the effectiveness of the method. The first part describes the expression of the improved Sukhof cooling formula and its accuracy verification; the second part describes the mathematical model of the electric power network, thermal network and electric-thermal coupling equipment; the third part describes in detail the flow of electric-thermal hybrid power flow calculation using the extended Newton-Raphson algorithm and the treatment of the coupled balancing nodes; in the fourth part, an example shows the effectiveness of this method and demonstrates the difference and regularity of electrical-heat network losses. The fifth part gives a summary of the whole paper.
(1) Mathematical description of cooling rules in heat network According to Sukhov's cooling formula, in a thermal network shown in Figure 1, the water temperature at the end terminal of thermal pipe with initial temperature T start and length L (Wu, 2006) is: among them, λ is the conductivity coefficient of the pipe, T a is the external environment temperature, C P is the specific heat ratio of hot water, ρ is the density of hot water, and m is the flow rate of hot water in the thermal pipe. This formula accurately describes the temperature drops in the thermal pipe network, but it contains exponent non-linear factors, which will greatly increase the iteration amount in the process of calculation. Meanwhile, the heat energy at the beginning and end of a thermal pipe can be expressed as: where Φ store start and Φ store end are the thermal power contained in the start and end terminal and ρ is the density of hot water. Substituting (Eq. 1) into (Eq. 3) obtains: At a load node in the thermal network, hot water flows back through the water return pipe after passing through the heat load. The difference in thermal power between supply-water and backwater is the available supplied thermal power of this load: where Tr is the temperature of backwater. Substituting (Eq. 4) into (Eq. 5) obtains: where Φ ava start C p ρm(T start − T r ) is the available heat power at the beginning of the thermal pipe. The heat power loss after hot water flows through the pipe with length L can be expressed as: By use of first-order Taylor expansion to Eq. 7 at L 0, the following equation can be obtained: It can be seen from (Eq. 8) that the heat power loss on the pipeline mainly relates to the temperature difference between the heating water temperature and the ambient temperature, the conductivity coefficient of the pipeline network and the length of the pipeline, and is also proportional to the length of the pipeline but independent of the flow rate, thus realizing the decoupling between the available heat power and the flow rate. This method can greatly reduce the iteration amount of power flow calculation in the heat network. The new expression of (Eq. 1) can be obtained as follows: This is the modified Sukhov cooling formula, which linearizes the relationship between temperature and the reciprocal flow rate.
Subtract T a from both sides of the equal sign in Eq. 9 to get: Define the relative temperatureT ' end T end − T a , T ' start T start − T a , that is, use the difference between hot water temperature and ambient temperature at both terminals of the pipeline to represent the temperature at both terminals of the pipeline, and let Ψ 1 − λL C p ρm be the cooling operator, then Eq. 10 can be simplified as: In a heat transfer pipe with a conductivity of 0.2W/(m·K), the comparison between the modified cooling factor and the actual value with different lengths and flow rates is shown in Table 1 and Figure 2. It can be seen that the modified cooling operator proposed here has a high accuracy within hundreds of meters range.

ELECTRICAL-THERMAL NETWORK MODEL
An electrically-thermally hybrid energy supply system, consisting of a electrically power network, a thermal network and electrically-thermally coupled devices, is shown in Figure 3. Thermal power network includes supply network and backwater network, which connect heat source and load and transfer energy with hot water. Electric power network mainly consists of power lines, which transmit active and reactive power in the form of electromagnetic energy. Coupling devices are connected to both the heat network and the electrical power grid to convert other forms of energy into electrical and thermal energy (e.g., gas-fired cogeneration), or to heat with electrical energy (e.g., heat pump, electric boiler, etc.).This section mainly describes their independent mathematical models.

Electrical Power Flow Model
The current vector, network admittance matrix and voltage vector of each node in the electrical power network model satisfy the following equation: among them, I i is the current flow of node i, U k is the voltage vector of node k, Y ik is the admittance matrix, and N is the number of nodes in this grid.
In order to obtain the node power equation, (Eq. 12) is rewritten by the equation between node power and current. The equation between node power and current is as follows: among them, P i and Q i are the active and reactive power of node i, respectively; I i is the conjugate vector of the voltage vector for node i. Substituting (Eq. 13) into (Eq. 12) gives: Rewriting Eq. 14 into the complex domain form and expanding it according to the real and imaginary parts respectively, the active and reactive power of the nodes in vector form can be obtained in Eq. 15:

Thermal Power Flow Model for Heat Network
In a thermal pipe network, according to the principle of flow balance, the net flow in each node's pipeline should be equal to the injected flow at that node, which is: where m q is the net outflow vector of the nodes, A is the nodebranch correlation matrix of the thermal network, and m is the flow variable of each branch. Meanwhile, in the heat pipe network circuit, the sum of the pressure drop of hot water in an enclosed circuit at the same direction is 0, then: BKm|m| 0 among them, B is the loop-branch correlation matrix of the heat pipe network, and K is the resistance coefficient matrix of the heat pipe network.  For heat sources or load nodes in a heating network, the heating power is: where Φ is the heat power of the node; T s is the supplying hot water temperature vector; and T o is the heat load out coming water temperature vector. For heat source nodes and load nodes where the effluent is not mixed with other return water branches, T o is equal to the return water temperature, that is: T o T r For heat loads without return water mixing T r,source Node i is not a mixing node (19) Substituting (Eq. 16) into (Eq. 18) gives.
For the hot water mix nodes in the heating pipe network shown in Figure 4, according to the principle of energy conservation, the thermal energy flowing into the node is equal to the thermal energy flowing out of the node. For node k, the energy conservation equation is shown in (Eq. 21): among them, m i,in and m j,out are respectively the pipeline flow with each incoming branch and the pipeline flow with each outgoing branch. T ' i,in and T ' out are the relative water temperature at the beginning of each pipe flowing into node k, and the relative temperature of hot water at node k.
Applying the temperature relationship in (Eq. 11) to Figure 4 obtains: where Ψ i is the Sokhov cooling factor of each incoming branch. Applying (Eq. 21) and (Eq. 22) to the supply and return water network, respectively, obtains: where T ' s and T ' r are the hot water relative temperature vectors of each node in the heating supply network and return water network, C s and C r are the matrices related to the heating and return pipe network and the flow, respectively. b s and b r are phasors related to heating supply and backwater temperatures, respectively.
For node i, find node j and pipe k where water flows from node j to node i through pipe k with flow rate m k , C s and b s can be obtained by the following expressions: Node i is not a mixing node Node i is not a mixing node 0 OtherWise (27) For the backwater network, C r and b r can be obtained by: C rii m i Node i is a mixing node 1 Node i is not a mixing node (28) Coupling Device Model

Gas Turbine
For gas-turbines, or other gas-fired internal-combustion cogeneration units, the conversion relationship between "fuelpower generation" and "fuel-heating" can be described as: where p out E and Φ out T are the output electric power and output heat power; η E is the electric power generating efficiency of the unit; c m is the unit's thermoelectric ratio; λ fuel and f fuel are the heat density of the fuel and the fuel consumption rate.

Steam Engine
For steam engines, the conversion relationship between "fuelelectric power generation" and "fuel-heating" can be described as: where p out E and Φ out T are the output electric power and output heat power; η E is the efficiency of the unit at full power generation; Z is the unit's thermoelectric ratio; λ fuel and f fuel are the heat density of the fuel and the fuel consumption rate.

Heat Pump
For various types of heat pump devices, the thermoelectric ratio is defined to characterize their hotspot conversion relationships: where Φ out T is the output thermal power of the electric heating device; p in E is the electrical power consumed; c m is the thermoelectric ratio. Since other media(air, soil and water) are used for heat exchange, the thermoelectric ratio is greater than 1.

Electric Boiler
An electric boiler is a device that uses electrical energy to make hot water for supply. The heating efficiency is defined to characterize its hotspot conversion relationship: where Φ out T is the output thermal power of the electric heating device; p in E is the electric power consumed; η ET is the heating efficiency, which is generally less than 1.

Circulating Pump
Circulating pump is a device that consumes electric power to provide water pressure to the heat pipe network to urge it to circulate heat. It is generally distributed at the outlet of each heat source. Its mathematical model is: where P P and η p are the power and efficiency consumed by the circulating pump; m p is the discharge flow of pressurized heat source; H p is the head required for pressure drop and thermal load pressure drop of the hot water circulation network. Circulating pumps generally consume much less power than heat loads and sources, so it is usually ignored in power balance calculation.

Newton-Raphson Algorithm Principle
Newton-Raphson algorithm is a non-linear mathematical method with good convergence and fast iteration speed. It has been widely used in electrical power flow calculation. This paper will extend it to the electric-thermal network mixed power flow calculation. By simultaneously establishing (Eqs. 15,17,20,23,24), the unified power flow deviation equation of the electric-heat hybrid energy supply network is obtained as: among them, x θ, |U|, m, T ' s , T ' r T is the state quantity to be solved from the equation Δf (x) 0, including the voltage phase angle, amplitude, hot water pipe network flow, hot water supply temperature and return water temperature of the thermal load nodes. The first formula of (Eq. 32) does not include the power system balance node, and the third formula does not include the thermal system balance node. P SP , Q SP , and Φ SP represent the active (except the power balance node), reactive and thermal power values (except the thermal balance node) given to the known nodes. A sl is a node-branch correlation matrix with reduced order after the thermal balance node is removed from the thermal pipe network. Let the true solution of Δf (x) 0 be x p . If the state value x k after k iterations is sufficiently close to the true solution x p , that is: Substituting (Eq. 39) into (Eq. 38) and performing Taylor expansion at x k , leaving behind only the first-order Taylor expansion, then: Thereby, Accordingly, △x k can be solved and superimposed on x k to obtain a new state quantity x k+1 . The k+1 iteration is: After n iterations, for a given error allowable value ε (usually 10 -3 ), if: Then the solution x n of Δf (x) 0 is obtained. Δf '(x k ) is a Jacobian matrix whose elements can be updated continuously according to the new value of x k during each iteration, which is also the core of Newton-Raphson algorithm to solve the nonlinear equation smoothly.

Jacobian Matrix Analysis
In (Eq. 38), the first two equations are the active and reactive power equations of the power system, and the last four equations are the equations of thermal system thermal power, thermal loop pressure, hot water temperature, and return water temperature. Let x e θ, |U| T andx h m, T ' s , T ' r T , then the Jacobian matrix can be expressed as: where the diagonal elements J ee and J hh represent the relationship between power flow and the state variable of electrical and thermal system itself, respectively.
among them: The non-diagonal elements are mainly affected by the electro-thermal coupling device acting as the balance nodes. For J eh , it is only necessary to consider the electric power of the power nodes in the power grid that are coupled to the power grid by the thermal balance nodes, while the other nodes in electrical power grid are not affected by thermal balance in heat network. The expression for determining the coupled electric power according to the heat source balance node is: Where A t−sourse,p is the row of the node-branch correlation matrix of the heat pipe network related to the balance node heat source. The subscript "p" represents the balance node. Therefore, define the operator, and find the partial derivative of (Eq. 49), J eh can be expressed as: For J he , if the electrical power system is connected to the large power grid, the change in electric power is balanced by the large power grid, so J he 0.If the power grid is an isolated system, it is balanced by the selected electrical power balance nodes, which will affect the power balance of the heat network. For the balance node of the electrical power grid, the electrical power can be calculated from the state quantity process iteration values, as follows: The subscript "p" denotes the balancing node and k denotes the k-th node associated with the balancing node in the electrical power system. Similarly, with the aid of operator φ p , for the balance node of electrical power system, the heat power coupled to the heat network can be expressed as:  Substitute it into (Eq. 44), J he can be expressed as:: Solving Algorithm Procedure For power system, nodes are divided into three categories: PQ node, PV node and power balance node. For thermal system, nodes are divided into heat source node, heat load node and heat source balance node. The selection principle and known and unknown quantities in calculation are shown in Table 2.
In an electric-thermal hybrid system, if the power balance node is not coupled with the thermal node and the thermal balance node is not coupled with the electric node, the power grid and the heat network can carry out the power flow solution independently. If the balance node of the electrical power system is coupled with the thermal system but the balance node of the thermal system is not coupled with the power system (e.g. gas or coal-fired boiler), then the power system is solved first, and then the coupled thermal power of the balancing node of the power system is substituted for the thermodynamic system as a given value, and finally the thermal system is solved. If the balance node of power system is not coupled to the thermal system, but the thermal system balance node is coupled to the power system (such as the power system is connected to the grid, the large power grid is used as the balance node, and the thermal system selects the gas unit as the balance node), then the thermal system is solved first to obtain the thermal balance node. The electric power coupled to the power system is used as a given value of the power node of the power system, and finally the power system is solved. In view of the most complex situation, in which the balance node of electrical power system couples to thermal system and also the balance node of thermal system couples to electrical power system, the main steps of hybrid power flow calculation of electricthermal network are as follows: (1) Analyze the electric-thermal network system and substitute the known quantity, select the balance node of electrical power system and thermal system, and define the mathematical descriptions of the electric-thermal coupling relationship; (2) Set the initial iteration value of the state quantities to be determined, calculate the return water temperature of the heat source according to (Eq. 21), calculate the power P t−sourse of the heat balance node according to (49), and calculate the power Φ e−sourse,p of the electrical balance node according to (Eq. 52).
(3) Calculate the deviation state function Δf (x k ) according to (Eq. 38), where P k t−sourse and Φ k e−sourse,p calculated in step 2 or step 7 are included in the given values of the power of the corresponding nodes in the grid and heat network in the first and third equations respectively. Take T o T r,sourse as the heat source return water temperature obtained in step 2 or step 8. (4) Calculate each element of the Jacobian matrix; (5) Solve the state quantities equation JΔx k −Δf (x k ) to obtain the required state quantity deviationΔx k ; (6) A new state quantity is obtained through x k+1 x k + Δx k iteration; (7) Calculate the return water temperature of the heat source according to (Eq. 21) after iteration, calculate the power of the heat balance node according to (Eq. 49) update, and recalculate the power of the electrical balance node coupling according to (Eq. 52); (8) Determine whether the latest solution after iteration meets the error range. If not, go back to step 3. (9) According to the obtained state quantities, the fuel consumption rate, the branch current of the power grid and the loss of the thermal pipeline are calculated according to the content in Section Coupling Device Model. (10) Output the calculation result.  The corresponding algorithm flow is shown in Figure 5.

EXAMPLE ANALYSES
In order to verify the accuracy of the electrical-thermal multienergy system modeling and power flow calculation method proposed in this paper, an off-grid electric-thermal energy supply system on Barry Island in (Liu et al., 2013) was selected for research. The system supplies energy to several electric and heating loads in five sub-regions. The voltage level of the main electrical power network is 11 kV. Its topology, branch length and numeral identification are shown in Figure 6. This example includes three thermo-electric coupling nodes, namely gas turbine, steam turbine and gas internal combustion engine. The gas turbine outlet voltage is 33 kV, and it is connected to the main power network through a stepdown transformer with a capacity of 15 MVA. The steam engine is selected as the balance node of the thermal system, and the gas turbine is used as the balance node of the electrical power system. According to the regional division, the size of the electrical and thermal loads are shown in Table 3. The division and known quantities of each source node are shown in Table 4. The power factor of the power system is 1. Outlet water temperature of each heat load is 30°C. The ambient temperature is 10°C and natural gas has a heat density of 9.7 kW·h/m 3 .The parameters of thermal lines can be found in (Xing et al., 2012), and those of electrical power lines can be found in Table 5.

Analysis of Simulation Results
Set the initial value of the power angle for the eight unbalanced nodes be 0, the initial value of the voltage amplitude for the six load nodes be 1.05 p.u., the initial value of the flow for the 32 heating network branches be 1L/s, and the initial values of the heating temperature and return water temperature for the 29 thermal loads be 70°C and 30°C, respectively. The error tolerance is 10 -3 . The results come out after 14 iterations. The simulation is carried out on a computer with Intel (R) core (TM) i5-3210M CPU@2.50 GHz, 4.00G RAM. The solution time of the conventional algorithm and with the modified Sukhov cooling operator is 9.2 and 5.3 s respectively. Table 6 shows the solved voltage-phase angle of each node in the electrical power system. It can be found that the voltage drop of each node is small, and there is no voltage exceeding the limit of 1.06 p.u. Table 7 shows the heating temperature and return water temperature of each load node in the thermal system. It can be found that the heating flow rates of source 1 to source 3 are 6.26, 4.81, and 2.25 L/s. Figure 7 shows a comparison of the electro-thermal power and gas consumption rates of the three coupled sources. It can be found that Source 1 has the largest power output in terms of the electro-thermal supply.

Adaptability Analysis
In order to prove the adaptability of this algorithm, the following scenarios are set for analysis based on the above example. Scenario 2: The electrical load remains the same and each thermal load increases by 0.02 MW; Scenario 3: The thermal load remains the same and each electrical load increases by 0.1 MW; Scenario 4: Each electrical load increases by 0.1 MW and each thermal load increases by   0.02 MW. Table 8 shows the calculation results of the algorithm described in this paper in various scenarios. It can be seen that with the electrical load increases, source 1 (the balancing node of the power system) realizes a new electricalthermal balance through the change of the main active power output; with the thermal load increases, source 2 (the balancing node of the thermal system) contributes the change of the main thermal power output to achieve new electro-thermal balance. When the electric and thermal loads increase, Source 1 and Source 2 respectively lead the power system and thermal system to achieve a new balance. Source 1, as the largest electric and thermal source, undertakes the main power output changes, and a fast solution by this method with modified Sokhov cooling operator can still be achieved.

Loss Analysis
For electricity-heating service providers, how to reduce the losses on the line and pipe network is directly related to the operating efficiency, so it is very important. The loss of the electric heating network in the four scenarios described in Section Adaptability Analysis is shown in Figure 8. It can be found that the transmission loss rate of the heat network is much larger than the transmission loss rate of the power grid. In terms of loss changes, with the increase of the supplied electrical load, the loss of the power network increases significantly, but for the thermal network, the increase of the supplied thermal load does not have a significant effect on the transmission loss of the heat pipe network. Figure 9 shows transmission loss and transmission efficiency of thermal pipeline networks under different external ambient temperature conditions in scenario 1. It