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 electric-thermal 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:

$T_{\text{end}}=\left(T_{\text{start}}-T_{\text{a}}\right)e^{-\frac{\lambda L}{C_{\text{p}}\rho m}}+T_{\text{a}}(1)$

among them, $\lambda$ is the conductivity coefficient of the pipe, T_a is the external environment temperature, C_P is the specific heat ratio of hot water, $\rho$ 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.

FIGURE 1

FIGURE 1. Typical heat network pipe.

Meanwhile, the heat energy at the beginning and end of a thermal pipe can be expressed as:

${\mathit{\Phi }}_{\text{start}}^{\text{store}}=C_{\text{p}}\rho m{T}_{\text{start}}(2)$

${\mathit{\Phi }}_{\text{end}}^{\text{stroe}}=C_{\text{p}}\rho m{T}_{\text{end}}(3)$

where ${\mathit{\Phi }}_{\text{start}}^{\text{store}}$ and ${\mathit{\Phi }}_{\text{end}}^{\text{store}}$ are the thermal power contained in the start and end terminal and $\rho$ is the density of hot water. Substituting (Eq. 1) into (Eq. 3) obtains:

${\mathit{\Phi }}_{\text{end}}^{\text{store}}={\mathit{\Phi }}_{\text{start}}^{\text{store}}{e}^{-\frac{\lambda L}{C_{\text{p}}\rho m}}+\left(1-e^{-\frac{\lambda L}{C_{\text{p}}\rho m}}\right)C_{\text{p}}\rho m{T}_{\text{a}}(4)$

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:

${\mathit{\Phi }}_{\text{end}}^{\text{ava}}=C_{\text{p}}\rho m\left(T_{\text{end}}-T_{\text{r}}\right)(5)$

where Tr is the temperature of backwater. Substituting (Eq. 4) into (Eq. 5) obtains:

${\mathit{\Phi }}_{\text{end}}^{\text{ava}}={\mathit{\Phi }}_{\text{start}}^{\text{ava}}{e}^{-\frac{\lambda L}{C_{\text{p}}\rho m}}+\left(1-e^{-\frac{\lambda L}{C_{\text{p}}\rho m}}\right)C_{\text{p}}\rho m\left(T_{\text{a}}-T_{\text{r}}\right)(6)$

where ${\mathit{\Phi }}_{\text{start}}^{\text{ava}}=C_{\text{p}}\rho m\left(T_{\text{start}}-T_{\text{r}}\right)$ 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:

$\Delta \mathit{\Phi }={\mathit{\Phi }}_{\text{start}}^{\text{ava}}-{\mathit{\Phi }}_{\text{end}}^{\text{ava}}=\frac{T_{\text{start}}-T_{\text{a}}}{T_{\text{start}}-T_{\text{r}}}\left(1-e^{-\frac{\lambda L\left(T_{\text{start}}-T_{\text{r}}\right)}{{\Phi }_{\text{start}}^{\text{ava}}}}\right){\mathit{\Phi }}_{\text{start}}^{\text{ava}}(7)$

By use of first-order Taylor expansion to Eq. 7 at L = 0, the following equation can be obtained:

$\Delta \mathit{\Phi }=\lambda \left(T_{\text{start}}-T_{\text{a}}\right)L(8)$

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:

$T_{\text{end}}=\frac{{\mathit{\Phi }}_{\text{start}}^{\text{ava}}-\Delta {\mathit{\Phi }}^{\text{ava}}}{C_{\text{p}}\rho m}=T_{\text{start}}-\frac{\lambda \left(T_{\text{start}}-T_{\text{a}}\right)L}{C_{\text{p}}\rho m}(9)$

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:

$\frac{T_{\text{end}}-T_{\text{a}}}{T_{\text{start}}-T_{\text{a}}}=1-\frac{\lambda L}{C_{\text{p}}\rho m}(10)$

Define the relative temperature$T_{\text{end}}^{\text{'}}=T_{\text{end}}-T_{\text{a}}$, $T_{\text{start}}^{\text{'}}=T_{\text{start}}-T_{\text{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 $\Psi =1-\frac{\lambda L}{C_{\text{p}}\rho m}$ be the cooling operator, then Eq. 10 can be simplified as:

$T_{\text{end}}^{\text{'}}=T_{\text{start}}^{\text{'}}\Psi (11)$

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.

TABLE 1

TABLE 1. Accuracy comparison of modified cooling factor.

FIGURE 2

FIGURE 2. Comparison of modified and actual cooling factor.

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.

FIGURE 3

FIGURE 3. Electrical-thermal hybrid power supply system.

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:

${\stackrel{·}{I}}_i={\displaystyle \sum _{k=1}^N{Y}_{ik}\stackrel{·}{U_k}}(12)$

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:

${\mathbf{S}}_i={\mathbf{P}}_i+j{\mathbf{Q}}_i={\stackrel{·}{\mathbf{U}}}_i{\widehat{\mathbf{I}}}_i\begin{array}{cc}& \left(i=\mathrm{1,2},\dots ,N\right)\end{array}(13)$

among them, P_i and Q_i are the active and reactive power of node i, respectively; ${\widehat{\mathbf{I}}}_i$ is the conjugate vector of the voltage vector for node i.

Substituting (Eq. 13) into (Eq. 12) gives:

${\mathbf{P}}_i+j{\mathbf{Q}}_i=\stackrel{·}{{\mathbf{U}}_i}{\displaystyle \sum _{k=1}^N{\widehat{Y}}_{jk}{\widehat{U}}_k\left(i=\mathrm{1,2},\dots ,N\right)}(14)$

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:

$\{\begin{array}{c}\mathbf{P}=\mathrm{Re}\left\{\stackrel{·}{\mathbf{U}}\left(\widehat{Y}\widehat{U}\right)\right\}\\ \mathbf{Q}=\mathrm{Im}\left\{\stackrel{·}{\mathbf{U}}\left(\widehat{Y}\widehat{U}\right)\right\}\end{array}(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:

$\mathit{Am}={\mathit{m}}_{\mathit{q}}(16)$

where m_q is the net outflow vector of the nodes, A is the node-branch 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:

$\mathit{BKm}\left|\mathit{m}\right|=0(17)$

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:

$\mathit{\Phi }=C_p\rho {\mathit{m}}_{\mathit{q}}\left({\mathit{T}}_{\text{s}}-{\mathit{T}}_{\text{o}}\right)(18)$

where $\mathit{\Phi }$ 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=\{\begin{array}{l}{T}_r\mathrm{ For heat loads without return water mixing}\\ T_{r,source}\mathrm{ Node }\mathit{i}\mathrm{ is not a mixing node}\end{array}(19)$

Substituting (Eq. 16) into (Eq. 18) gives.

$\mathit{\Phi }=C_p\rho \mathbf{A}\mathit{m}\left({\mathit{T}}_{\text{s}}-{\mathit{T}}_{\text{o}}\right)(20)$

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):

${\displaystyle \sum _{i=1}^A\left(m_{i,\text{in}}{T}_{i,\text{in}}^{\text{'}}\right)}=\left({\displaystyle \sum _{j=1}^B{m}_{j,\text{out}}}\right)T_{\text{out}}^{\text{'}}(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,\text{in}}^{\text{'}}$ and $T_{\text{out}}^{\text{'}}$ 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.

FIGURE 4

FIGURE 4. Mixing node in thermal network.

Applying the temperature relationship in (Eq. 11) to Figure 4 obtains:

$T_{i,\text{in}}^{\text{'}}{\Psi }_i=T_{\text{out}}^{\text{'}}\text{i }==1\text{,}2\text{,\&,A}(22)$

where ${\Psi }_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:

${\mathbf{C}}_{\text{s}}{\mathbf{T}}_{\text{s}}^{\text{'}}={\mathbf{b}}_{\text{s}}(23)$

${\mathbf{C}}_{\text{r}}{\mathbf{T}}_{\text{r}}^{\text{'}}={\mathbf{b}}_{\text{r}}(24)$

where ${\mathbf{T}}_{\text{s}}^{\text{'}}$ and ${\mathbf{T}}_{\text{r}}^{\text{'}}$ are the hot water relative temperature vectors of each node in the heating supply network and return water network, ${\mathbf{C}}_{\text{s}}$ and ${\mathbf{C}}_{\text{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:

$C=\{{\displaystyle \begin{array}{l}\sum _k{m}_k\mathrm{ Node }\mathit{i}\mathrm{ is a mixing node}\\ \mathrm{1 Node }\mathit{i}\mathrm{ is not a mixing node}\end{array}}(25)$

$C_{sij}=\{\begin{array}{l}-m_k{\psi }_{k\mathrm{ }\mathrm{A}\mathrm{n}\mathrm{d}\mathrm{ }j\mathrm{ }\mathrm{i}\mathrm{s}\mathrm{ }\mathrm{a}\mathrm{ }\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d}\mathrm{ }\mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}}^{\mathrm{N}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{ }i\mathrm{ }\mathrm{i}\mathrm{s}\mathrm{ }\mathrm{a}\mathrm{ }\mathrm{m}\mathrm{i}\mathrm{x}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{n}\mathrm{o}\mathrm{d}\mathrm{e}}\\ 0\mathrm{ OtherWise}\end{array}(26)$

$b_{si}=\{\begin{array}{l}{\mathrm{T}}_{\mathit{sj}}^{\mathit{\text{'}}}{\Psi }_{\mathit{k}}\mathrm{ Node }\mathit{i}\mathrm{ is not a mixing node}\\ {\displaystyle \sum _{j,k}{m}_k{\mathrm{T}}_{\mathit{sj}}^{\mathit{\text{'}}}{\Psi }_{\mathit{k}}{\mathrm{ }}_{\mathrm{Node }\mathit{i}\mathrm{ is not a mixing node}}^{\mathrm{Node }\mathit{i}\mathrm{ is not a mixing node}}}\\ 0\mathrm{ OtherWise}\end{array}(27)$

For the backwater network, C_r and b_r can be obtained by:

$C_{rii}=\{\begin{array}{l}{m}_i\mathrm{ Node }\mathit{i}\mathrm{ is a mixing node}\\ 1\mathrm{ Node }\mathit{i}\mathrm{ is not a mixing node}\end{array}(28)$

$C_{rij}=\{\begin{array}{l}-m_k{\psi }_k\mathrm{ Node }\mathit{i}\mathrm{ is a mixing node}\\ 0\mathrm{ Node }\mathit{i}\mathrm{ is not a mixing node}\end{array}(29)$

$b_{ri}=\{\begin{array}{l}{m}_{qi}{T}_{oi}^{\text{'}}\mathrm{ Node }\mathit{i}\mathrm{ is a mixing node}\\ T_{oi}^{\text{'}}\mathrm{ Node }\mathit{i}\mathrm{ is not a mixing node}\end{array}(30)$

Coupling Device Model

Gas Turbine

For gas-turbines, or other gas-fired internal-combustion cogeneration units, the conversion relationship between “fuel-power generation” and “fuel-heating” can be described as:

$p_{\text{E}}^{\text{out}}={\eta }_{\text{E}}{\lambda }_{\text{fuel}}{f}_{\text{fuel}}(31)$

${\mathbf{\Phi }}_{\text{T}}^{\text{out}}=c_{\text{m}}{p}_{\text{E}}^{\text{out}}(32)$

where $p_{\text{E}}^{\text{out}}$ and ${\mathbf{\Phi }}_{\text{T}}^{\text{out}}$ are the output electric power and output heat power; ${\eta }_{\text{E}}$ is the electric power generating efficiency of the unit; $c_{\text{m}}$ is the unit's thermoelectric ratio; ${\lambda }_{\text{fuel}}$ and $f_{\text{fuel}}$ are the heat density of the fuel and the fuel consumption rate.

Steam Engine

For steam engines, the conversion relationship between “fuel-electric power generation” and “fuel-heating” can be described as:

$Z=\frac{{\mathbf{\Phi }}_{\text{T}}^{\text{out}}}{P_{\text{con}}-p_{\text{E}}^{\text{out}}}(33)$

$P_{\text{con}}={\eta }_{\text{E}}{\lambda }_{\text{fuel}}{f}_{\text{fuel}}(34)$

where $p_{\text{E}}^{\text{out}}$ and ${\mathbf{\Phi }}_{\text{T}}^{\text{out}}$ are the output electric power and output heat power; ${\eta }_{\text{E}}$ is the efficiency of the unit at full power generation; Z is the unit's thermoelectric ratio; ${\lambda }_{\text{fuel}}$ and $f_{\text{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:

${\mathbf{\Phi }}_{\text{T}}^{\text{out}}=c_{\text{m}}{p}_{\text{E}}^{\text{in}}(35)$

where ${\mathbf{\Phi }}_{\text{T}}^{\text{out}}$ is the output thermal power of the electric heating device; $p_{\text{E}}^{\text{in}}$ is the electrical power consumed; $c_{\text{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:

${\mathbf{\Phi }}_{\text{T}}^{\text{out}}={\eta }_{\text{ET}}{p}_{\text{E}}^{\text{in}}(36)$

where ${\mathbf{\Phi }}_{\text{T}}^{\text{out}}$ is the output thermal power of the electric heating device; $p_{\text{E}}^{\text{in}}$ is the electric power consumed; ${\eta }_{\text{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:

$P_{\text{P}}=\frac{m_{\text{p}}g{H}_{\text{p}}}{{10}^6{\eta }_{\text{p}}}(37)$

where $P_{\text{P}}$ and ${\eta }_{\text{p}}$ are the power and efficiency consumed by the circulating pump; $m_{\text{p}}$ is the discharge flow of pressurized heat source; $H_{\text{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.

Electric-Thermal Network Hybrid Power Flow Calculation Method

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:

$\Delta f\left(\mathbf{x}\right)=\{\begin{array}{c}\begin{array}{c}\Delta \mathbf{P}\\ \Delta \mathbf{Q}\end{array}\\ \begin{array}{c}\Delta \mathit{m}\\ \Delta \mathit{p}\end{array}\\ \begin{array}{c}\Delta {\mathbf{T}}_{\text{s}}^{\text{'}}\\ \Delta {\mathbf{T}}_{\text{r}}^{\text{'}}\end{array}\end{array}=\{\begin{array}{l}{\mathbf{P}}^{\text{SP}}-\mathrm{Re}\left\{\stackrel{·}{\mathbf{U}}\left(\widehat{Y}\widehat{U}\right)\right\}\\ {\mathbf{Q}}^{\text{SP}}-\mathrm{Im}\left\{\stackrel{·}{\mathbf{U}}\left(\widehat{Y}\widehat{U}\right)\right\}\\ C_{\text{p}}\rho {\mathbf{A}}_{\text{sl}}\mathit{m}\left({\mathit{T}}_{\text{s}}-{\mathit{T}}_{\text{o}}\right)-{\mathit{\Phi }}^{\text{SP}}\\ \mathit{BKm}\left|\mathit{m}\right|\\ {\mathbf{C}}_{\text{s}}{\mathbf{T}}_{\text{s}}^{\text{'}}-{\mathbf{b}}_{\text{s}}\\ {\mathbf{C}}_{\text{r}}{\mathbf{T}}_{\text{r}}^{\text{'}}-{\mathbf{b}}_{\text{r}}\end{array}(38)$

among them, $\mathbf{x}={\lfloor \mathbf{\theta },\left|\mathbf{U}\right|,\mathbf{m},{\mathbf{T}}_{\text{s}}^{\text{'}},{\mathbf{T}}_{\text{r}}^{\text{'}}\rfloor }^{\text{T}}$ is the state quantity to be solved from the equation $\Delta f\left(\mathbf{x}\right)=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. ${\mathbf{P}}^{\text{SP}}$, ${\mathbf{Q}}^{\text{SP}}$, and ${\mathit{\Phi }}^{\text{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 $\Delta f\left(\mathbf{x}\right)=0$ be ${\mathbf{x}}^{\ast }$. If the state value ${\mathbf{x}}^k$ after k iterations is sufficiently close to the true solution ${\mathbf{x}}^{\ast }$, that is:

${\mathbf{x}}^{\ast }={\mathbf{x}}^k+△{\mathbf{x}}^k(39)$

Substituting (Eq. 39) into (Eq. 38) and performing Taylor expansion at ${\mathbf{x}}^k$, leaving behind only the first-order Taylor expansion, then:

$\Delta f\left({\mathbf{x}}^{\ast }\right)=\Delta f\left({\mathbf{x}}^k+△{\mathbf{x}}^k\right)=\Delta f\left({\mathbf{x}}^k\right)+\Delta f\text{'}\left({\mathbf{x}}^k\right)△{\mathbf{x}}^k=0(40)$

Thereby,

$\Delta f\text{'}\left({\mathbf{x}}^k\right)△{\mathbf{x}}^k=-\Delta f\left({\mathbf{x}}^k\right)(41)$

Accordingly, $△{\mathbf{x}}^k$ can be solved and superimposed on ${\mathbf{x}}^k$ to obtain a new state quantity ${\mathbf{x}}^{k+1}$. The k+1 iteration is:

${\mathbf{x}}^{k+1}={\mathbf{x}}^k+△{\mathbf{x}}^k(42)$

After n iterations, for a given error allowable value $\epsilon$ (usually 10^–3), if:

$\left|\Delta f\left({\mathbf{x}}^n\right)\right|\le \epsilon (43)$

Then the solution ${\mathbf{x}}^n$ of $\Delta f\left(\mathbf{x}\right)=0$ is obtained. $\Delta f\text{'}\left({\mathbf{x}}^k\right)$ is a Jacobian matrix whose elements can be updated continuously according to the new value of ${\mathbf{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 ${\mathbf{x}}_e={\lfloor \mathbf{\theta },\left|\mathbf{U}\right|\rfloor }^{\text{T}}$ and${\mathbf{x}}_h={\lfloor \mathbf{m},{\mathbf{T}}_{\text{s}}^{\text{'}},{\mathbf{T}}_{\text{r}}^{\text{'}}\rfloor }^{\text{T}}$, then the Jacobian matrix can be expressed as:

$\mathbf{J}=[\begin{array}{l}{\mathbf{J}}_{\text{ee}}\\ {\mathbf{J}}_{\text{he}}\end{array}\begin{array}{l}{\mathbf{J}}_{\text{eh}}\\ {\mathbf{J}}_{\text{hh}}\end{array}]=[\begin{array}{l}\frac{\partial {\left[\Delta \mathbf{P},\Delta \mathbf{Q}\right]}^{\text{T}}}{\partial {\mathbf{x}}_{\text{e}}^{\text{T}}}\\ \frac{\partial {\left[\Delta \mathit{\Phi },\Delta \mathbf{p},\Delta {\mathbf{T}}_{\text{s}}^{\text{'}},\Delta {\mathbf{T}}_{\text{r}}^{\text{'}}\right]}^{\text{T}}}{\partial {\mathbf{x}}_{\text{e}}^{\text{T}}}\end{array}\begin{array}{l}\frac{\partial {\left[\Delta \mathbf{P},\Delta \mathbf{Q}\right]}^{\text{T}}}{\partial {\mathbf{x}}_{\text{h}}^{\text{T}}}\\ \frac{\partial {\left[\Delta \mathit{\Phi },\Delta \mathbf{p},\Delta {\mathbf{T}}_{\text{s}}^{\text{'}},\Delta {\mathbf{T}}_{\text{r}}^{\text{'}}\right]}^{\text{T}}}{\partial {\mathbf{x}}_{\text{h}}^{\text{T}}}\end{array}](44)$

where the diagonal elements ${\mathbf{J}}_{\text{ee}}$ and ${\mathbf{J}}_{\text{hh}}$ represent the relationship between power flow and the state variable of electrical and thermal system itself, respectively.

${\mathbf{J}}_{\text{ee}}=\left[\begin{array}{cc}\begin{array}{c}\mathrm{Re}\left({\mathbf{J}}_{{\mathbf{S}}_{\mathbf{\theta }}}\right)\\ \mathrm{Im}\left({\mathbf{J}}_{{\mathbf{S}}_{\mathbf{\theta }}}\right)\end{array}& \begin{array}{c}\mathrm{Re}\left({\mathbf{J}}_{{\mathbf{S}}_{\mathbf{V}}}\right)\\ \mathrm{Im}\left({\mathbf{J}}_{{\mathbf{S}}_{\mathbf{V}}}\right)\end{array}\end{array}\right](45)$

among them:

${\mathbf{J}}_{{\mathbf{S}}_{\mathbf{\theta }}}=\frac{\partial \Delta {\mathbf{S}}_i}{\partial {\mathbf{\theta }}_k}=\{\begin{array}{c}\text{j}\stackrel{·}{U_i}{\widehat{Y}}_{ik}{\widehat{U}}_k\\ \text{j}\stackrel{·}{U_i}{\widehat{Y}}_i{\widehat{U}}_i-\text{j}\stackrel{·}{U_i}{\displaystyle \sum _{n=1}^N{\widehat{Y}}_{in}{\widehat{U}}_n}\end{array}\begin{array}{c}k\ne i\\ k=i\end{array}(46)$

${\mathbf{J}}_{{\mathbf{S}}_{\mathbf{V}}}=\frac{\partial \Delta {\mathbf{S}}_i}{\partial {\left|\mathbf{U}\right|}_k}=\{\begin{array}{c}-\stackrel{·}{U_i}{\widehat{Y}}_{ik}{e}^{-j{\theta }_k}\\ -\stackrel{·}{U_i}{\widehat{Y}}_{ii}{e}^{-j{\theta }_i}-\frac{\text{j}\stackrel{·}{U_i}{\displaystyle \sum _{n=1}^N{\widehat{Y}}_{in}{\widehat{U}}_n}}{{\left|U\right|}_i}\end{array}\begin{array}{c}k\ne i\\ k=i\end{array}(47)$

$\begin{array}{l}{\mathbf{J}}_{\text{hh}}=\left[\begin{array}{ccc}\begin{array}{c}\frac{\partial \Delta \mathit{\Phi }}{\partial {\mathbf{m}}^{\text{T}}}\\ \frac{\partial \Delta \mathbf{p}}{\partial {\mathbf{m}}^{\text{T}}}\\ \frac{\partial \Delta {\mathbf{T}}_{\text{s}}^{\text{'}}}{\partial {\mathbf{m}}^{\text{T}}}\\ \frac{\partial \Delta {\mathbf{T}}_{\text{r}}^{\text{'}}}{\partial {\mathbf{m}}^{\text{T}}}\end{array}& \begin{array}{c}\frac{\partial \Delta \mathit{\Phi }}{\partial {\mathbf{T}}_{\text{s}}^{\text{' T}}}\\ \frac{\partial \Delta \mathbf{p}}{\partial {\mathbf{T}}_{\text{s}}^{\text{' T}}}\\ \frac{\partial \Delta {\mathbf{T}}_{\text{s}}^{\text{'}}}{\partial {\mathbf{T}}_{\text{s}}^{\text{' T}}}\\ \frac{\partial \Delta {\mathbf{T}}_{\text{r}}^{\text{'}}}{\partial {\mathbf{T}}_{\text{s}}^{\text{' T}}}\end{array}& \begin{array}{c}\frac{\partial \Delta \mathit{\Phi }}{\partial {\mathbf{T}}_{\text{r}}^{\text{' T}}}\\ \frac{\partial \Delta \mathbf{p}}{\partial {\mathbf{T}}_{\text{r}}^{\text{' T}}}\\ \frac{\partial \Delta {\mathbf{T}}_{\text{s}}^{\text{'}}}{\partial {\mathbf{T}}_{\text{r}}^{\text{' T}}}\\ \frac{\partial \Delta {\mathbf{T}}_{\text{r}}^{\text{'}}}{\partial {\mathbf{T}}_{\text{r}}^{\text{' T}}}\end{array}\end{array}\right]\\ =\left[\begin{array}{ccc}\begin{array}{c}{C}_{\text{p}}\rho \text{diag}\left\{\left({\mathit{T}}_{\text{s}}-{\mathit{T}}_{\text{o}}\right)\right\}{\mathbf{A}}_{\text{sl}}\\ \mathit{2BK}\left|\mathit{m}\right|\\ \frac{\partial \Delta {\mathbf{b}}_{\text{s}}}{\partial {\mathbf{m}}^{\text{T}}}\\ \frac{\partial \Delta {\mathbf{b}}_{\text{r}}}{\partial {\mathbf{m}}^{\text{T}}}\end{array}& \begin{array}{c}{C}_{\text{p}}\rho \text{diag}\left\{{\mathbf{A}}_{\text{sl}}\mathit{m}\right\}\\ 0\\ {\mathbf{C}}_{\text{s}}\\ 0\end{array}& \begin{array}{c}0\\ 0\\ 0\\ {\mathbf{C}}_{\text{r}}\end{array}\end{array}\right]\end{array}(48)$

The non-diagonal elements are mainly affected by the electro-thermal coupling device acting as the balance nodes. For ${\mathbf{J}}_{\text{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:

Among them,

${\mathit{\Phi }}_{\text{t}-\text{sourse,p}}=C_{\text{p}}\rho {\mathbf{A}}_{\text{t}-\text{sourse,p}}\mathbf{m}\left({\mathbf{T}}_{\text{s,source}}-{\mathbf{T}}_{\text{r,source}}\right)(49)$

Where ${\mathbf{A}}_{\text{t}-\text{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), ${\mathbf{J}}_{\text{eh}}$ can be expressed as:

$\begin{array}{l}{\mathbf{J}}_{\text{eh}}=\frac{\partial \Delta {\left[\mathbf{P},\mathbf{Q}\right]}^{\text{T}}}{\partial {\mathbf{x}}_{\text{h}}^{\text{T}}}\\ =\left[\begin{array}{cc}\begin{array}{c}{\phi }_{\text{p}}{C}_{\text{p}}\rho \text{diag}\left\{\left({\mathit{T}}_{\text{s,source}}-{\mathit{T}}_{\text{r,source}}\right)\right\}{\mathbf{A}}_{\text{t}-\text{sourse,p}}\\ \\ 0\end{array}& \begin{array}{cc}\begin{array}{c}0\\ \\ 0\end{array}& \begin{array}{c}0\\ \\ 0\end{array}\end{array}\end{array}\right]\end{array}(50)$

For ${\mathbf{J}}_{\text{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 ${\mathbf{J}}_{\text{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:

$P_{\text{e}-\text{sourse,p}}=\mathrm{Re}\left\{\text{j}{U}_{\text{p}}\left({\widehat{Y}}_{\text{p,}k}{\widehat{U}}_k\right)\right\}(51)$

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 ${\phi }_{\text{p}}$, for the balance node of electrical power system, the heat power coupled to the heat network can be expressed as:

${\varphi }_{e-source,p}\mathrm{ }=\{\begin{array}{l}{C}_m{P}_{e-source,p}\mathrm{ For gas engines}\\ \mathit{Z}\left({\eta }_{\mathit{E}}{\lambda }_{fuel}{f}_{fuel}-P_{e-source,\mathrm{p}}\right)\mathrm{ For Steam engines}\\ -{\mathit{C}}_{\mathit{m}}{P}_{e-source,\mathrm{p}}\mathrm{ For Heat pumps}\\ -{\eta }_{\mathit{ET}}{P}_{e-source,\mathrm{p}}\mathrm{ For electrical boilers}\end{array}(52)$

Substitute it into (Eq. 44), ${\mathbf{J}}_{\text{he}}$ can be expressed as::

$\begin{array}{l}{\mathbf{J}}_{\text{he}}=\frac{\partial \Delta {\left[\mathit{\Phi },\mathbf{p},{\mathbf{T}}_{\text{s}}^{\text{'}},{\mathbf{T}}_{\text{r}}^{\text{'}}\right]}^{\text{T}}}{\partial {\mathbf{x}}_{\text{e}}^{\text{T}}}\\ =\left[\begin{array}{c}\frac{1}{{\phi }_{\text{p}}}\mathrm{Re}\left\{\text{j}{U}_{\text{p}}\left({\widehat{Y}}_{\text{p,}k}{\widehat{U}}_k\right)\right\}\begin{array}{cc}& \frac{1}{{\phi }_{\text{p}}}\mathrm{Re}\left\{-{\widehat{U}}_{\text{p}}{\widehat{Y}}_{\text{p,}k}{e}^{-j{\theta }_{\text{p,}k}}\right\}\end{array}\\ 0\begin{array}{ccc}& & 0\end{array}\\ \begin{array}{c}0\\ 0\end{array}\begin{array}{ccc}& & \begin{array}{c}0\\ 0\end{array}\end{array}\end{array}\right]\end{array}(53)$

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.

TABLE 2

TABLE 2. Determination and classification of electric/thermal hybrid system.

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 electric-thermal 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_{\text{t}-\text{sourse}}$ of the heat balance node according to (49), and calculate the power ${\mathbf{\Phi }}_{\text{e}-\text{sourse,p}}$ of the electrical balance node according to (Eq. 52).

(3) Calculate the deviation state function $\Delta f\left({\mathbf{x}}^k\right)$ according to (Eq. 38), where $P_{\text{t}-\text{sourse}}^k$ and ${\mathbf{\Phi }}_{\text{e}-\text{sourse,p}}^k$ 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 ${\mathit{T}}_{\text{o}}={\mathit{T}}_{\text{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 $\mathbf{J}\Delta {\mathbf{x}}^k=-\Delta f\left({\mathbf{x}}^k\right)$ to obtain the required state quantity deviation$\Delta {\mathbf{x}}^k$;

(6) A new state quantity is obtained through ${\mathbf{x}}^{k+1}={\mathbf{x}}^k+\Delta {\mathbf{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.

FIGURE 5

FIGURE 5. Flowchart of hybrid power flow calculation.

Example Analyses

In order to verify the accuracy of the electrical-thermal multi-energy 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 step-down 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³.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.

FIGURE 6

FIGURE 6. Topology of the electrical & thermal system.

TABLE 3

TABLE 3. Load parameters of electrical & thermal system.

TABLE 4

TABLE 4. Source parameters of electrical & thermal system.

TABLE 5

TABLE 5. Parameters of the power system.

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.

TABLE 6

TABLE 6. Solutions of the electrical system.

TABLE 7

TABLE 7. Solutions of the thermal system.

FIGURE 7

FIGURE 7. Source output power and gas consumption rate.

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 electrical-thermal 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.

TABLE 8

TABLE 8. Solutions of more scenes.

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 8

FIGURE 8. Electrical and thermal power loss analysis.

Figure 9 shows transmission loss and transmission efficiency of thermal pipeline networks under different external ambient temperature conditions in scenario 1. It can be found that higher ambient temperature leads to smaller loss in the process of hot water transmission. The lower the ambient temperature, the greater the temperature difference between the hot water and the outside world, the greater the loss of the pipe network and the lower the efficiency of the pipe network, which is consistent with the description of (Eq. 8) in this paper. Therefore, for the heating service, the heating temperature can be appropriately adjusted according to the ambient temperature change to achieve more economical operation.

FIGURE 9

FIGURE 9. Relationship between thermal power loss and ambient temperature.

Conclusion

In this paper, the model of electric-thermal mixed power flow and the model of key equipment of electric-thermal coupling node are constructed. On the basis of deriving the improved Suhof cooling operator for thermal pipeline network, an extended Newton-Raphson algorithm is proposed, and the method of obtaining asymmetric elements in Jacobi matrix is described. The main conclusions are as follows:

(1) The use of the modified Sukhov cooling operator proposed in this paper eliminates the need for exponential operations and helps to achieve fast algorithm convergence during the iteration process.

(2) The electric-thermal coupling equipment of the balance node has a great influence on the power flow calculation and solution of the electric-thermal hybrid system.

(3) The expanded Newton-Raphson power flow algorithm proposed in this paper can better adapt to the solution of strongly coupled multi-quantities.

(4) The electric grid loss is related to the transmission load, and the heat network loss is more directly related to the environmental temperature difference. Optimizing the difference between hot water temperature and ambient temperature is an effective way to achieve the economic operation of the heat network.

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.

Author Contributions

HL proposed the theory. YL deduced the theory. XC verified the theory by developing the example calculation. KJ was responsible for the paper writing. While ZZ provided the simulation. All authors contributed to the article and approved the submitted version.

Conflict of Interest

Author ZZ was employed by Chongqing Shibei Power Supply Company, State Grid Corporation of China.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by grant 2019YFE0111500 from the National key R & D plan of China.

References

Geidl, M., and Andersson, G. (2007). Optimal power flow of multiple energy carriers. IEEE Trans. Power Syst. 22, 145–155. doi:10.1109/TPWRS.2006.888988

CrossRef Full Text | Google Scholar

Geidl, M. (2007). Integrated modeling and optimization of multi-carrier energy systems. Doctoral Thesis. Swiss Federal Institute of Technology. doi:10.3929/ethz-a-005377890

CrossRef Full Text | Google Scholar

Good, N., Zhang, L., Navarro-Espinosa, A., and Mancarella, P. (2015). High resolution modeling of multi-energy domestic demand profiles. Appl. Energy 137, 193–210. doi:10.1016/j.apenergy.2014.10.028

CrossRef Full Text | Google Scholar

He, K. L., Chen, Q., Ma, H., Zhao, T., and Hao, J. H. (2020). An isomorphic multi-energy flow modeling for integrated power and thermal system considering nonlinear heat transfer constraint. Energy 211, 119003. doi:10.1016/j.energy.2020.119003

CrossRef Full Text | Google Scholar

Huang, W., Zhang, N., Cheng, Y., Yang, J., Wang, Y., and Kang, C. (2020). Multienergy networks analytics: standardized modeling, optimization, and low carbon analysis. Proc. IEEE 108, 1411–1436. doi:10.1109/JPROC.2020.2993787

CrossRef Full Text | Google Scholar

Kargarian, A., and Rahmani, M. (2015). Multi-microgrid energy systems operation incorporating distribution-interline power flow controller. Elec. Power Syst. Res. 129, 208–216. doi:10.1016/j.epsr.2015.08.015

CrossRef Full Text | Google Scholar

Kriechbaum, L., Gradl, P., Reichenhauser, R., and Kienberger, T. (2020). Modelling grid constraints in a multi-energy municipal energy system using cumulative exergy consumption minimisation. Energies 13, 3900. doi:10.3390/en13153900

CrossRef Full Text | Google Scholar

Lingmin, C., Jiekang, W., Fan, W., Huiling, T., and Yan, X. (2020). Energy flow optimization method for multi-energy system oriented to combined cooling, heating and power. Energy 211, 118536. doi:10.1016/j.energy.2020.118536

CrossRef Full Text | Google Scholar

Liu, X. (2013). Combined analysis of electricity and heat networks. Ph. D. thesis. Cardiff, Wales: Cardiff University.

Google Scholar

Liu, X., Wu, J., Jenkins, N., and Bagdanavicius, A. (2016). Combined analysis of electricity and heat networks. Appl. Energy 162, 1238–1250. doi:10.1016/j.apenergy.2015.01.102

CrossRef Full Text | Google Scholar

Loukarakis, E., and Mancarella, P. (2017). A sequential programming method for multi-energy districts optimal power flow. IEEE Manchester PowerTech 17044946. doi:10.1109/PTC.2017.7980850

CrossRef Full Text | Google Scholar

Mancarella, P. (2014). MES (multi-energy systems): an overview of concepts and evaluation models. Energy 65, 1–17. doi:10.1016/j.energy.2013.10.041

CrossRef Full Text | Google Scholar

Moeini-Aghtaie, M., Abbaspour, A., Fotuhi-Firuzabad, M., and Hajipour, E. (2014). A decomposed solution to multiple-energy carriers optimal power flow. IEEE Trans. Power Syst. 29, 707–716. doi:10.1109/TPWRS.2013.2283259

CrossRef Full Text | Google Scholar

Omalley, M., and Kroposki, B. (2013). Energy comes together: the integration of all systems. IEEE Power Energy Mag. 11, 18–23. doi:10.1109/MPE.2013.2266594

CrossRef Full Text | Google Scholar

Salpakari, J., Mikkola, J., and Lund, P. D. (2016). Improved flexibility with large-scale variable renewable power in cities through optimal demand side management and power-to-heat conversion. Energy Convers. Manag. 126, 649–661. doi:10.1016/j.enconman.2016.08.041

CrossRef Full Text | Google Scholar

Shabanpour-Haghighi, A., and Seifi, A. R. (2015). Simultaneous integrated optimal energy flow of electricity, gas, and heat. Energy Convers. Manag. 101, 579–591. doi:10.1016/j.enconman.2015.06.002

CrossRef Full Text | Google Scholar

Shao, C., Wang, X., Shahidehpour, M., Wang, X., and Wang, B. (2017). An MILP-based optimal power flow in multicarrier energy systems. IEEE Transactions on Sustainable Energy 8, 239–248. doi:10.1109/TSTE.2016.2595486

CrossRef Full Text | Google Scholar

Wu, F. (2006). Analyses and improves on design method of centralized heating system networks. Harbin, China: Harbin Institute of Technology.

Google Scholar

Xing, Y., Bagdanavicius, A., Lannon, S. C., Pirouti, M., and Bassett, T. (2012). “Low temperature district heating network planning with the focus on distribution energy losses.” in International Conference on applied energy ICAE2012-A10103, July 5–8, 2012, Suzhou, China

Google Scholar