Research on Optimal Energy Flow Calculation of an Electro-Thermal Coupled System via an Interior-Point Approach

With the contradiction between energy consumption and ecological protection becoming increasingly prominent, it has become a serious issue to study an energy system with high efficiency, low pollution, and strong controllability. The application of the backpressure unit and the peak-shaving electric boiler makes the power grid and the heat grid connected and becomes an electro-thermal coupling system, which can be studied as a representative of the integrated energy system. In this study, to make full use of energy, a power system model considering a large number of distributed energy injections and a multi-source thermal system model are constructed. After the energy flow is solved, the optimization method of the electrothermal coupled system is studied to minimize the fuel consumption and obtain the minimum operating cost. In addition, a more comprehensive algorithm is proposed in this study, and the actual model is used in the fourth section of the article to verify that the algorithm can be practically applied.

1 Introduction

The development of the integrated energy system (IES) can make the allocation of resources more rational and optimal, which is a good way to deal with the huge energy consumption and resource depletion in the current era, and it is also a foundation for the complete transformation of the country’s energy structure in the future (Moeini-aghtaie et al., 2014). In the traditional energy systems such as power grids and thermal networks, each energy system has an independent operating system, the coupling between them is not tight, and the energy efficiency is low. The coupling between integrated energy systems is tight, the energy utilization rate is high, and the ability to absorb renewable energy is greatly improved (Liu et al., 2016; Sun et al., 2020; Zhang et al., 2022). When the power grid is actually running to calculate the power flow, the forward-backward method and the Newton–Raphson method are mostly used (Hu et al., 2020). These two methods have a strong theoretical system in the power grid, and their ideas are also relatively mature. However, when it is applied to an electro-thermal combined system, the solution program becomes complicated, the number of iterations is large, and it is difficult to realize the combination between the various systems. At present, many scholars at home and abroad have proposed energy flow solving methods for integrated energy systems (Jia et al., 2014; Xia et al., 2017; Yu et al., 2021). Liu et al. (2021) have refined the heat network and considered the method of solving the energy flow in the case of heat loss. Guo et al. (2019) constructed a unified energy flow solution method based on energy hubs. Li et al. (2015) have considered the effect of adding energy storage elements to the whole system.

The modeling of each part of the network and the energy flow solution of the electro-thermal coupled system are the basis of the optimal energy flow research. As the most important part of this study, the control variables are not given as the general energy flow calculation but need to be determined through optimization, thus forming a nonlinear programming problem with multiple constraints, and the solution methods are also diverse to simplify the nonlinear programming method represented by the gradient method, and the interior point method is used more. The research on the optimal power flow has been very extensive. Zhang et al. (2019) considered the construction of the heat network model more comprehensively and built a refined energy flow solution model on the basis of the heating system. Jiang and Liang. (2015) considered the characteristics of poor performance in traditional AC-DC hybrid power flow calculation and uses the Gaussian method for an iterative solution to improve the applicability of the entire system. To realize the optimal allocation of energy and reduce energy consumption, this study constructs a simplified integrated energy system model and uses a unified solution method to achieve a fast and accurate solution, simple solution process, high calculation accuracy, and simple programming and does not require too many iterations (Krause et al., 2011; Malley and Kroposki, 2013; Chen et al., 2020). When a large number of PV nodes are injected into the system, the PV nodes can be equivalent to improve the traditional forward and backward generation method. In addition, the interior point method chosen in this study in the process of solving the optimal energy flow does not need to find the initial value in the feasible region, and it will adapt quickly when the system scale increases (Wang et al., 2019; Wang et al., 2020; Wang et al., 2021a; Wang et al., 2021b). In the fourth section of the study, an optimal energy flow calculation model including a distributed power grid and a multi-heat source radiant heat grid is established, and two examples are used to analyze the characteristics of the two solutions proposed.

2 Research on the Energy Flow Calculation Model and Algorithm of the Combined Electric-Heat System

2.1 Power System Model

The running state of the power system is not constant, and a power flow equation can be listed in each state. If the column is written as the node voltage, it can be expressed as follows:

${\dot{U}}_i=U_i\angle {\theta }_i=U_i\left(\cos \theta +j\sin \theta \right).(1)$

The node power equation of power flow calculation can be expressed as follows:

$\{\begin{array}{c}\mathrm{\Delta }P=P_{G,i}-P_{Load,i},\\ \mathrm{\Delta }Q=Q_{G,i}-Q_{Load,i},\end{array}(2)$

where $P_{G,i}$ and $Q_{G,i}$ are the active and reactive power delivered by the generator into the network; $P_{Load,i}$ and $Q_{Load,i}$ are the load active and reactive power corresponding to the node.

2.2 Thermal System Model

The thermodynamic system model is composed of a hydraulic model and a thermodynamic model. According to the hydraulic model of the heat network, the amount of water flowing through the pipeline and the injected water flow at the heat load node can be obtained. The following flow continuity equation is used to express as follows:

$A_{l}m=m_c,(3)$

where $A_l$ is the network correlation matrix of the heat load node relative to each pipeline, and $m$ is the water flow vector through the pipe.

A thermal model can be used to determine the temperature of each node of the thermal network, including the following equations. Eq. 4 is the equation of node injection heat power, mass flow, and node temperature. Eq. 5 is the relationship between the temperature of the head and tail of the pipeline. Eq. 6 is the equation for the outflow of hot water from the heat source through different heat networks and finally mixed at heat loads.

$\varphi ={\varphi }_i-{\varphi }_{Ei}=c_{p}m\left(T_{il}-T_{ol}\right),(4)$

$T_{end}=\left(T_{start}-T_e\right)e^{-\frac{\lambda d}{c_{p}m}}+T_e+\mu ,(5)$

$\left({\displaystyle \sum _{a=1}^n{m}_{out,a}}\right)T_{out}={\displaystyle \sum _{b=1}^m{m}_{in,b}{T}_{in,b}},(6)$

where $\varphi$ is the node injection heat power vector, ${\varphi }_i$ is the node heat load power vector, ${\varphi }_{Ei}$ is the node electric boiler injection heat power vector, $T_{il}$ is the temperature of hot water provided by each node, $T_{ol}$ is the temperature of hot water flowing back, $T_{end}$ is the temperature of the hot water as it leaves the pipe, $T_{start}$ is the temperature at which the hot water enters the pipe, $T_e$ is the ambient temperature, $\lambda$ is the coefficient of heat conduction per unit length of the pipeline, $d$ is the transmission distance of the pipeline, $\mu$ is the temperature compensation parameter, $m_{out,a}$ is the branch flow of the branch out of the pipeline for the branch $a$, and $T_{out}$ is the temperature of hot water flowing back to the pipeline at different locations.

There is often more than one heat source in a thermal system. Figure 1A is the schematic diagram of the radiant heat system model with multiple heat sources. In this model, two heat sources are used to heat the heat load at the same time. Then, the two heat sources are transformed into two equivalent heat sources shown in Figure 1B to heat the heat load. Therefore, the multi-heat source radiation type heat network is converted into multiple single heat source radiation-type heat networks. The transformation is as follows:

$m_{k_{ij}}=\frac{m_{H_1}}{m_{H_1}+\dots +m_{H_i}}{m}_{k_i}.(7)$

FIGURE 1

FIGURE 1. (A) Multi-source radiant thermal system model. (B) Single-heat source radiant heat network model.

This is rewritten in the matrix form as

$\left[\begin{array}{c}{m}_{k_{1\theta }}\\ \cdots \\ m_{k_{N\theta }}\end{array}\right]=\frac{m_{H_1}}{m_{H_1}+\cdots +m_{H_i}}\left[\begin{array}{c}{m}_{k_1}\\ \cdots \\ m_{k_i}\end{array}\right],(8)$

where, $m_{k_{ij}}$ (i = 1, … , n, n represents the number of pipes branching from the node, j = 1, … , m, m represents the number of heat sources in the heat network) is the corresponding single heat source radiation heat network after conversion. $m_{H_i}$ (i = 1, ..., n) is the pipe flow input by the corresponding heat source to the node, and $m_{k_i}$ (i = 1, ..., n) is the branch flow before conversion, that is, the multi-heat source radiation type heat network corresponds to the pipeline flow.

${\displaystyle \sum _{i=1}^m{\varphi }_{l_i}={\varphi }_l},(9)$

where, ${\varphi }_{l_i}$ (i = 1, ..., m, m is the number of heat sources in the heat network) is the equivalent heat load power after the multi-heat source radiation heat network is equivalent to multiple single heat source radiation heat networks, and ${\varphi }_l$ is the equivalent heat before the node thermal load power.

2.3 Thermal-Electro Coupled Model

Coupling elements of power systems and thermal systems such as thermoelectric units and heat pumps have the property of generating heat while supplying power. A backpressure unit is a widely used thermoelectric unit, which does not contain a condenser. The exhaust heat in the back pressure unit will be fully utilized, and the thermal efficiency is high. The disadvantage is that the generating capacity of the backpressure unit should be based on heat production, and the regulation sensitivity is not strong. Therefore, the backpressure unit is generally selected as the heat source to play a central role in the thermal system, and the electric boiler is selected to meet the demand for peak shaving.

The thermal power and electrical power emitted by the backpressure unit satisfy the following equation:

${\zeta }^{-1}\cdot {\varphi }_{BY}=P_{BY},(10)$

where $\zeta$ is the ratio of heat production to electricity production, which is generally a constant value, ${\varphi }_{BY}$ is the heat power produced by the backpressure unit, and $P_{BY}$ is the electricity generated by the backpressure unit.

The equation of thermal power and electrical power emitted by a peak-shaving electric boiler is as follows:

${\delta }^{-1}\cdot {\varphi }_{EB}=P_{EB},(11)$

where $\delta$ is the heat-to-electricity ratio, which is generally a constant value, ${\varphi }_{EB}$ is the thermal power emitted by the peak-shaving electric boiler, and $P_{EB}$ is the electric power generated by the peak-shaving electric boiler.

2.4 Calculation Method of the Energy Flow in the Electro-Thermal Coupled System

At present, Newton–Raphson’s energy flow solution method is the most widely used. However, the Newton–Raphson method has the shortcomings of strict requirements on the initial value and more iterations. Therefore, this study proposes an improved forward-backward generation method, which improves the drawbacks that the traditional forward-backward generation method cannot deal with PV nodes, and is extended to the radiant thermal system with multiple heat sources, which reduces the number of iterations and the requirement for initial values. The improvement of the thermal system simplifies the model and improves the calculation speed.

When distributed energy is added to the grid, PV nodes will become more. But the traditional push-back method cannot handle this. Therefore, it can be improved to allow power flow calculations to continue. Figure 2 shows the flow solution diagram of the improved forward-backward substitution method.

FIGURE 2

FIGURE 2. Flow chart of the power flow solution.

The power loss formula is as follows:

$\mathrm{\Delta }S=\mathrm{\Delta }P+j\mathrm{\Delta }Q=I^2\left(R_i+j{X}_i\right)=\frac{P_i^2+Q_i^2}{U_N^2}\left(R_i+j{X}_i\right)(12)$

where $\Delta S$ is power loss, $I$ is the line current, $R_i$, $X_i$ are line parameters, and $U_N$ is the rated voltage of the line.

Next, the voltage drop formula of each line is calculated as follows,

$\mathrm{\Delta }U=Z_i\cdot \mathrm{\Delta }{S}^{\ast }=\left(R_i+j{X}_i\right)\left(\mathrm{\Delta }P-j\mathrm{\Delta }Q\right)=R_i\mathrm{\Delta }P+X_i\mathrm{\Delta }Q,(13)$

where $\Delta U$ is the longitudinal component of the voltage drop, $Z_i$, $R_i$, $X_i$ are line parameters, $\Delta P$ is the active loss, and $\Delta Q$ is the reactive loss.

For the $PV$ nodes, if $\Delta P=0$, the following equation is satisfied:

$\mathrm{\Delta }Q=X_i^{-1}\cdot \mathrm{\Delta }U(14)$

where $\Delta Q$ is the reactive power loss, and $X_i$ is the line parameter.

Since the injection of distributed generation into the distribution network will increase the number of $PV$ nodes in the power system, a large number of $PV$ nodes cannot be calculated by using the forward-backward substitution method. If the initial reactive power value at the power source is set to 0, the $PV$ node can be equivalent to a $PQ$ node first to facilitate subsequent solutions, and the reactive power at the power source is corrected at each iteration to ensure the real-time data.

The correction equation is as follows:

$\mathrm{\Delta }Q=\gamma \left(Q_i-{\displaystyle \sum Q_{ij}}\right)\left(1+e^{-t}\right)=\gamma \left(\frac{-P_{i}R-\mathrm{\Delta }{U}_b{U}_i}{X_a}-Q_{ij}\right)\left(1+e^{-t}\right),(15)$

where $\gamma \in \left(-1,1\right)$ is the calculation step length. Generally $\gamma =0.1$, $Q_i$ is the reactive power calculated by the $PV$ node, $Q_{ij}$ is the reactive power of all branches and nodes connected to the $PV$ node, $P_i$ is the power of the $PV$ node, and $X_a$ are line parameters.

The fact that whether the voltage of each node satisfies the convergence condition of the following equation is judged. If not, the voltage value is updated before the power flow calculation, and the equation is corrected until the voltage satisfies the convergence condition.

The convergence condition is as follows:

$\max \left\{\left|\mathrm{\Delta }{U}_i^{\left(k\right)}\cdot {\epsilon }^{-1}\right|\right\}=\max \left\{\left|\left(U_i^{\left(k\right)}-U_0\right)\cdot {\epsilon }^{-1}\right|\right\}<1,(16)$

where $\epsilon$ is the convergence accuracy, $U_i^{\left(k\right)}$ is the voltage value obtained by the kth iteration, and $U_0$ is the initial voltage value.

The aforementioned idea can be extended to the thermal system to be equivalent to the power system. The temperature of the heat network is analogous to the voltage of the power grid, the heat flow is analogous to the current, the heat load node is equivalent to the load, the heat source node is equivalent to the power supply, and the node with known heating temperature is an analogy to the balance node in the power system. Then, according to the coupled part modeling mentioned in 2.3, the sum of thermal power of all equilibrium nodes in the thermal system is transformed into the power of electric load to realize the transformation process of the power flow solution.

3 Research on the Optimal Energy Flow Calculation Model and Algorithm of the Electric-Heat Combined System

3.1 Problem Description

One of the main research problems in this study is the construction of the optimal energy flow model for the entire system and the solution method. The purpose of solving the optimal energy flow is to make the entire system run efficiently and to obtain security guarantees at the same time. Adjustment through certain algorithms of the output of each energy unit can achieve the optimal distribution of energy among various systems. When the model is built, the grid is used as the basis. The thermal grid and the grid are connected together through coupling elements at a certain node, and the thermal power of the thermal grid is converted into electrical power and introduced into the grid. To calculate the optimal energy flow, the objective function must be selected first. In this study, the goal is to minimize the operating cost of the entire system. This is a nonlinear programming problem that considers multiple constraints of independent networks and coupled parts. The following assumptions should be made before modeling:

1) The operating conditions of each power source and a heat source are known (if there is a parallel connection between power sources or heat sources, it will not be considered);

2) The output of each power source and heat source are known;

3) The topologies of the heat grid and the grid are known;

4) The node load size is known;

5) The parameters of the overall network are known.

The mathematical model of the optimal energy flow can be expressed as

$\{\begin{array}{c}\min f\left(u,x\right),\\ s.t.\begin{array}{c}g\left(u,x\right)=0,\\ h\left(u,x\right)\le 0,\end{array}\end{array}(17)$

where $u$ is generally an independent variable that can adjust the active and reactive power output of each power supply and the temperature and flow of hot water flowing out of the heat source, and $x$ is generally a variable that follows $u$, such as the voltage of each node in the power grid and the branch power, the heat network. The temperature of the water flow at each node and the flow of the pipeline at other places are determined; $\min f\left(u,x\right)$ is the minimum system operating cost in this study, $g\left(u,x\right)=0$ is the power balance equation of the power grid and the heat network, and $h\left(u,x\right)\le 0$ is the inequality constraint.

3.2 Objective Function

If the objective function is to minimize the operating cost, the corresponding coal consumption should be minimized. The heat source selected in this study is a backpressure unit, and a peak-shaving electric boiler is used to assist it, so it can be considered that the objective function established in this study is based on heat production and production. Electricity is a quadratic fitting function of the variable, and the objective function is as follows:

$\min {\displaystyle \sum _{i=1}^N\left(c_{0\left(i\right)}+c_{1\left(i\right)}{P}_{BY\left(i\right)}+c_{2\left(i\right)}{P}_{BY\left(i\right)}^2+c_{3\left(i\right)}{D}_{BY\left(i\right)}+c_{3\left(i\right)}{D}_{BY\left(i\right)}^2+c_{5\left(i\right)}{D}_{BY\left(i\right)}{P}_{BY\left(i\right)}\right)},(18)$

where $i$ is the backpressure unit from 1 to n, $c_{0\left(i\right)}\sim c_{5\left(i\right)}$ is the fitting constant, $P_{BY\left(i\right)}$ is the electrical power generated by the backpressure unit, and $D_{BY\left(i\right)}$ is the thermal power generated by the backpressure unit.

3.3 Restrictions

3.3.1 Grid Active and Reactive Power Balance Constraints

The equation constraints of the power grid can be considered as the active and reactive power equations of the power flow calculation, and the balance condition is that the inflow power of each node in the power grid is equal to the sum of the outflow and losses:

$\{\begin{array}{c}{\displaystyle \sum _{a=1}^N{P}_{BY\left(a\right)}-P_a-U_a}{\displaystyle \sum _{b=1}^n{U}_b}\left(G_{ab}\cos {\theta }_{ab}+B_{ij}\sin {\theta }_{ab}\right)=0\\ {\displaystyle \sum _{a=1}^N{Q}_{BY\left(a\right)}-Q_a-U_a}{\displaystyle \sum _{b=1}^n{U}_b}\left(B_{ij}\sin {\theta }_{ab}+G_{ab}\cos {\theta }_{ab}\right)=0\end{array}(19)$

where ${\displaystyle \sum _{a=1}^N{P}_{BY\left(a\right)}}$ and ${\displaystyle \sum _{a=1}^N{Q}_{BY\left(a\right)}}$ are the sum of the active power and reactive power of the backpressure units connected to 1 to N unit, and $P_a$ and $Q_a$ are the active and reactive power of the node.

3.3.2 Power Balance Constraints of the Heating Network

The power balance constraints of the heating network include the electric power balance of backpressure unit nodes, the thermal power of heat source nodes, the thermal power used by users and the thermal power supplied by electric boilers, and the temperature balance of hot water flowing into pipes.

${\varphi }_i-{\varphi }_{Ei}-c_{p}m\left(T_{il}-T_{ol}\right)=0,(20)$

$\left(T_{start}-T_e\right)e^{-\frac{\lambda d}{c_{p}m}}+T_e+\mu -T_{end}=0,(21)$

${\varphi }_i+c_{p}m\left(T_{il}-T_{ol}\right)-{\varphi }_{EB}-D_{BY\left(i\right)}=0,(22)$

where ${\varphi }_i$ is the node heat load power vector, ${\varphi }_{Ei}$ is the node electric boiler injection heat power vector, $T_{il}$ is the temperature of hot water provided by each node, $T_{ol}$ is the temperature of hot water flowing back, $T_{end}$ is the temperature of the hot water as it leaves the pipe, $T_{start}$ is the temperature at which the hot water enters the pipe, $T_e$ is the ambient temperature, $\lambda$ is the coefficient of heat conduction per unit length of the pipeline, $d$ is the transmission distance of the pipeline, and $\mu$ is the temperature compensation parameter.

3.3.3 Constraint Conditions of Power Grid Inequality

For all backpressure units, the upper and lower power limits of active and reactive power are met:

$\{\begin{array}{c}{P}_{BY\left(i\right)}^{\min }\le P_{BY\left(i\right)}\le P_{BY\left(i\right)}^{\max },\\ Q_{BY\left(i\right)}^{\min }\le Q_{BY\left(i\right)}\le Q_{BY\left(i\right)}^{\max },\end{array}(23)$

where $P_{BY\left(i\right)}^{\min }$ and $P_{BY\left(i\right)}^{\max }$ are the minimum and maximum active power output of the back pressure unit, and $Q_{BY\left(i\right)}^{\min }$ and $Q_{BY\left(i\right)}^{\max }$ are the limit values of the reactive power output of the back pressure unit.

Similarly, all peak-shaving electric boilers also meet the upper and lower power limits of active and reactive power:

$\{\begin{array}{c}{P}_{BE\left(i\right)}^{\min }\le P_{BE\left(i\right)}\le P_{BE\left(i\right)}^{\max },\\ Q_{BE\left(i\right)}^{\min }\le Q_{BE\left(i\right)}\le Q_{BE\left(i\right)}^{\max },\end{array}(24)$

where $P_{BE\left(i\right)}^{\min }$ and $P_{BE\left(i\right)}^{\max }$ are the minimum and maximum active power output of the peak-shaving electric boiler, and $Q_{BE\left(i\right)}^{\min }$ and $Q_{BE\left(i\right)}^{\max }$ are the limit values of the reactive power output of the peak-shaving electric boiler.

To meet the stability of power grid operation, it is necessary to restrict the amplitude of the voltage and the phase angle difference of the voltage:

$\{\begin{array}{c}{U}_{a\left(i\right)}^{\min }\le U_{a\left(i\right)}\le U_{a\left(i\right)}^{\max }\\ {\left|{\theta }_{a\left(i\right)}^{\min }\le U_{a\left(i\right)}\right|}_{\min }\le \left|{\theta }_{a\left(i\right)}\le {\theta }_{a\left(j\right)}\right|\le \left|{\theta }_{a\left(i\right)}\le {\theta }_{a{\left(j\right)}_{\max }}\right|\end{array}(25)$

where $U_{a\left(i\right)}^{\min }$ and $U_{a\left(i\right)}^{\max }$ are the upper and lower limits of the voltage amplitude of the $i$ node, and ${\left|{\theta }_{a\left(i\right)}-{\theta }_{a\left(j\right)}\right|}_{\min }$ and ${\left|{\theta }_{a\left(i\right)}-{\theta }_{a\left(j\right)}\right|}_{\max }$ are the limits of the voltage phase angle difference between the two nodes.

3.3.4 Inequality Constraint Conditions of the Heating Network

If the variables $m$, $T_{il}$ ,$T_{ol}$, $T_{end}$, $T_{start}$, $T_e$, and $D_{BY\left(i\right)}$ of the heating network exceed a certain range, the operation of the whole combined system will become unstable, so the upper and lower limits of these variables should be constrained:

$\{\begin{array}{c}{m}_{a\left(i\right)}^{\min }\le m_{a\left(i\right)}\le m_{a\left(i\right)}^{\max },\\ T_{il\left(i\right)}^{\min }\le T_{il\left(i\right)}\le T_{il\left(i\right)}^{\max },\\ T_{ol\left(i\right)}^{\min }\le T_{ol\left(i\right)}\le T_{ol\left(i\right)}^{\max },\\ T_{end\left(i\right)}^{\min }\le T_{end\left(i\right)}\le T_{end\left(i\right)}^{\max },\\ T_{start\left(i\right)}^{\min }\le T_{start\left(i\right)}\le T_{start\left(i\right)}^{\max }.\\ T_{e\left(i\right)}^{\min }\le T_{e\left(i\right)}\le T_{e\left(i\right)}^{\max },\\ D_{BY\left(i\right)}^{\min }\le D_{BY\left(i\right)}\le D_{BY\left(i\right)}^{\max },\end{array}(26)$

3.3.5 Power Balance Constraint of the Electrothermal Coupling Node

When the backpressure unit is used as the hub of the power grid and the heating network and when the peak-shaving electric boiler assists in regulation, there is the following power balance equation at the junction of the two networks:

$\{\begin{array}{c}{\delta }^{-1}\cdot {\varphi }_{EB}-P_{EB}=0,\\ {\zeta }^{-1}\cdot {\varphi }_{BY}-P_{BY}=0.\end{array}(27)$

The aforementioned formula is the conversion relationship between thermal power and electric power of the backpressure unit and electric boiler.

3.4 Calculation Method of the Optimal Energy Flow in the Electrothermal Coupling System.

Among many methods for solving optimal energy flow, the nonlinear programming method and interior point method are widely used. Compared with the former, the interior point method has the advantages of fast calculation speed and better convergence. When the system is huge, the calculation amount is not as complex as other methods, and the initial value is not particularly strict. Therefore, the interior point method is selected to solve the optimal energy flow of the whole joint system in this study. The traditional interior point method is always in the feasible region in the iterative process. Therefore, it is often difficult to choose the initial value. However, with the emergence of the central trajectory method, it is not necessary to find the initial value in the feasible region of the established model, and the range has changed greatly than the previous one. It only needs to meet the conditions of relaxation variables and disturbance factors, and the requirements for the initial value are no longer so strict.

At first, the interior point method appeared to calculate the optimal power flow in the power grid. When constructing the objective function, it is not necessary to list the state variables and control variables separately like other power flow models. In the interior point method, these two variables are not distinguished, and they are uniformly represented by $x$. In addition to the active power, reactive power, voltage amplitude, and phase angle of each node, $x$ in this study also includes the thermal power generated by the heat source in the heating network, the temperature and flow rate of water flowing out to the pipeline, the temperature information of the outlet and inlet of each pipeline, and the water flow rate in the pipeline.

After the aforementioned analysis, the optimal power flow model is simplified as the following formula:

$\{\begin{array}{c}\min f\left(x\right),\\ s.t.\begin{array}{l}p\left(x\right)=0,\hfill \\ q_{\min }\le q\left(x\right)\le q_{\max },\hfill \end{array}\end{array}(28)$

where $f\left(x\right)$ is the minimum operating cost function described in 3.2, $p\left(x\right)=0$ is the equality constraint, $q_{\min }\le q\left(x\right)\le q_{\max }$ is the inequality constraint, and $x$ is the set composed of variables $P_{BY}$, $Q_{BY}$, $P_{BE}$, $Q_{BE}$, $U_a$, ${\theta }_a$, $D_{BY}$, $D_{BE}$, $m$, $T_{il}$, $T_{ol}$, $T_{end}$, and $T_{start}$.

For Equation 28, it is difficult to solve the model. Two relaxation variables $m$ and $n$ can be introduced to transform it into the same equality constraint as $p\left(x\right)=0$:

$\{\begin{array}{c}\min f\left(x\right),\\ s.t.\begin{array}{l}p\left(x\right)=0,\hfill \\ q\left(x\right)+m=q_{\max },\hfill \\ q\left(x\right)-n=q_{\min },\hfill \\ mn>0,\hfill \end{array}\end{array}(29)$

where $m$ and $n$ are relaxation variables, both of which are positive numbers.

If the solution directly according to the aforementioned formula may make the result beyond the feasible range, Eq. 29 can be combined with the objective function to limit the solution result within the determined range, and once it approaches the boundary, the result increases rapidly and returns to the region again.

$\{\begin{array}{c}\min f\left(x\right)-\upsilon {\displaystyle \sum _{a=1}^n\log \frac{m_a}{n_a},}\\ s.t.\begin{array}{l}p\left(x\right)=0,\hfill \\ q\left(x\right)+m=q_{\max },\hfill \\ q\left(x\right)-n=q_{\min },\hfill \end{array}\end{array}(30)$

where $\upsilon$ is the perturbation factor; the perturbation factor is used to eliminate inequality constraints and transform them into equality constraints for subsequent solutions.

For equality constraints, Lagrange functions can be constructed to find the extremum of each variable:

$L=f\left(x\right)-\upsilon {\displaystyle \sum _{a=1}^n\log }\frac{m_a}{n_a}-r^{T}p\left(x\right)-s^T\left[q\left(x\right)+m-q_{\max }\right]-t^T\left[q\left(x\right)-n-q_{\min }\right],(31)$

where $r^T={\left[r_1,r_2,\cdots ,r_n\right]}^T$, $s^T={\left[s_1,s_2,\cdots ,s_n\right]}^T$, and $t^T={\left[t_1,t_2,\cdots ,t_n\right]}^T$ are multipliers introduced in the constructor.

The partial derivatives of $x$, $m$, $n$, $r$, $s$, and $t$ are obtained by using the Lagrange function, and the nonlinear equations are constructed when the derivative value is zero:

$\{\begin{array}{l}\frac{\partial L}{\partial x}={\nabla }_{x}f\left(x\right)-{\nabla }_{x}p\left(x\right)r-{\nabla }_{x}q\left(x\right)\left(s+t\right)=0,\hfill \\ \frac{\partial L}{\partial r}=p\left(x\right)=0,\hfill \\ \frac{\partial L}{\partial s}=q\left(x\right)+m-q_{\max }=0,\hfill \\ \frac{\partial L}{\partial t}=q\left(x\right)-n-q_{\min }=0,\hfill \\ \frac{\partial L}{\partial n}=NTE-\upsilon E=0,\hfill \\ \frac{\partial L}{\partial M}=MSE+\upsilon E=0,\hfill \end{array}(32)$

where $M$, $N$, $S$, and $T$ are diagonal matrices with $a$ elements, and $E$ is an identity matrix with $a$ elements.

Newton's method is used to solve the aforementioned differential equation and simplify it, and the simplified modified equation matrix form is as follows:

$\left[\begin{array}{cc}A& \nabla p\left(x\right)\\ {\nabla }_x^{T}p\left(x\right)& 0\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }x\\ \mathrm{\Delta }r\end{array}\right]=\left[\begin{array}{c}B\\ C\end{array}\right],(33)$

where the matrices $A$, $B$, and $C$ are expressed by the following equations:

$A={\nabla }_x^{2}f\left(x\right)+{\nabla }_x^{2}p\left(x\right)+{\nabla }_x^{2}h\left(x\right)\left(s+t\right)+\nabla h\left(x\right)\left[M^{-1}S-N^{-1}Z\right]{\nabla }_x^{T}h\left(x\right),(34)$

$B={\nabla }_{x}f\left(x\right)-{\nabla }_{x}p\left(x\right)r-{\nabla }_{x}q\left(x\right)\left(s+t\right)+{\nabla }_{x}h\left(x\right)\left[N^{-1}Z-U^{-1}W\right]{\nabla }_x^{T}h\left(x\right),(35)$

$C=-p\left(x\right).(36)$

According to Eq. 33, a group of solutions closest to the optimal solution can be solved.

4 Example Verification

4.1 Example Verification of the Energy Flow Calculation of the Electrothermal Coupling System

To prove that the aforementioned method can be practically applied, an actual model is used in Figure 3 as an example to conduct experiments.

FIGURE 3

FIGURE 3. Integrated energy testing system.

In this example, the heat grid and the power grid are coupled through the cogeneration unit. The heat network part includes two heat source nodes source 1 and source 2 using a backpressure unit, a peak-shaving electric boiler located at node 5, and 14 heat load nodes. The peak modulation ratio is 0.45, the thermoelectric ratio is 1.4, and the heat load is 0.5 MW. The length of each pipeline is set to 1.58 km, $\lambda$ is set to 0.289, and the diameter of the pipeline is set to 100 mm. According to the transformation method of the multi-source radiant heat network in Section 2.2, the multi-source radiant heat network is transformed into several single-source radiant heat networks and then the energy flow solution is performed. Table 1 is the heating system pipeline flow, and Table 2 is the temperature at each node in the heat network.

TABLE 1

TABLE 1. Heat system pipeline flow.

TABLE 2

TABLE 2. Temperature at each node in the heat network.

The thermal system and power system are coupled through node 6. The thermal power of node 6 is 1.7347 MW, which is converted to 2.2551 MW. The power grid part selects 14 node distribution networks, node 1 is set to 1.00 $pu$, node 2 and node 6 are $PV$ nodes, the voltage amplitudes are 0.993 $pu$ and 0.964 $pu$, and its node is $PQ$ node. According to the flow chart in Figure 2, the 1st and 2nd are set to $PQ$ nodes for the first iterative, and then, the program is returned to calculate the reactive power correction. The reactive power correction is substituted for the next iterative calculation. After three iterations, the convergence criteria are met, and the result is output. The power flow calculation results obtained by this method are compared with those obtained by the traditional Newton–Raphson method of power flow calculation, as shown in Figure 4 which verifies the accuracy of this method.

FIGURE 4

FIGURE 4. Comparison of power flow calculation results.

4.2 Example Verification of the Optimal Energy Flow Calculation of the Electrothermal Coupling System

In this chapter, the model of the 14-node grid and 14-node heating network shown in Figure 5 is selected for example verification. Two heat sources $S1$ and $S2$ in the heating network are connected with 6 and 11 nodes in the grid, respectively. The parameters of the two heat source nodes of the back pressure unit and heat load nodes in the combined system are set as shown in Table 3.

FIGURE 5

FIGURE 5. Test model of the electrothermal combined system.

TABLE 3

TABLE 3. Setting of optimal power flow operating conditions.

When setting the objective function, the minimum operating cost is selected as the objective function, and the fitting constant in the objective function is selected in Table 4. After setting the initial values and upper and lower limits of each parameter, the simulation is carried out according to the optimal energy flow calculation method described in Section 3.4, and the running cost of the system is minimized by adjusting the thermal power and electric power of the power supply and heat source. The comparison diagram of pipeline flow calculation results of optimal power flow and general power flow is shown in Figure 6 and Table 5. The comparison diagram of the supply and return temperature of each node is shown in Figure 7 and Table 6.

TABLE 4

TABLE 4. Fitting constant with operating cost as the objective function.

FIGURE 6

FIGURE 6. Comparison chart of pipeline flow calculation results.

TABLE 5

TABLE 5. Water flow of the heat network pipeline (kg/s).

FIGURE 7

FIGURE 7. Comparison diagram of the supply and return temperature.

TABLE 6

TABLE 6. Temperature comparison of heating network nodes.

The node temperature obtained from the optimal power flow calculation is determined by optimizing the control variables, while the control variables in the ordinary power flow are initially given. Therefore, it is more economical to introduce the results of the optimal power flow into the actual operation, and the coal consumption will be lower, thus achieving the purpose of the optimal energy flow calculation.

5 Conclusion

In this study, an improved forward push-back substitution method is proposed for the combined electric and heat system. First, models of individual systems and coupled parts are built, and the calculation model of the electro-thermal combined system is constructed. The improved forward-backward substitution method is used for energy flow calculation to realize the analysis and calculation of the energy flow of the whole combined system. Finally, through the calculation and analysis of the example of the comprehensive energy test system, it is proved that the idea proposed in this study can be practically used. In the part of the study of the optimal energy flow calculation method, the interior point method is used to find the parameters that meet the optimal operation conditions of the system. By comparing the simulation results of energy flow calculation and optimal energy flow calculation, it can be obtained that the interior point method is feasible and economical to calculate the optimal energy flow of the electrothermal combined system.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author Contributions

ZL was responsible for the main idea, simulation part, and writing work of this study. YZ and JX were responsible for the mathematical derivation. QS participated in the coordination of the study and reviewed the manuscript.

Funding

This work was supported by the Key R&D Program of Liaoning Province (2020JH2/10300101), the Liaoning Revitalization Talents Program (XLYC1907138), the Technology Innovation Talent Fund of Shenyang (RC200252), and the Key R&D Program of Shenyang (GG200252).

Conflict of Interest

QS was employed by the State Grid Liaoning Electric Power Co., Ltd.

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors acknowledge the funding of the Key R&D Program of Liaoning Province (2020JH2/10300101), the Liaoning Revitalization Talents Program (XLYC1907138), the Technology Innovation Talent Fund of Shenyang (RC200252), and the Key R&D Program of Shenyang (GG200252).

References

Chen, B. B., Sun, H. B., Yin, G. X., Wu, W. C., and Guo, Q. L. (2020). Energy Circuit Theory of Integrated Energy System Analysis (II): Hydraulic Circuit and Thermal Circuit. Proc. CSEE 40, 2133–2142. doi:10.13334/j.0258-8013.pcsee.200098

CrossRef Full Text | Google Scholar

Guo, Z. G., Lei, J. Y., Ma, X. Y., Yuan, Z. Y., Dong, B., and Yu, H. (2019). Modeling and Calculation Methods for Multi-Energy Flows in Large-Scale Integrated Energy System Containing Electricity, Gas, and Heat. Proc. CSU-EPSA. 31, 97–102. doi:10.19635/j.cnki.csu-epsa.000283

CrossRef Full Text | Google Scholar

Hu, X., Shang, C., Cheng, H. Z., Wang, L., Chen, D. W., and Li, Y. (2020). Review and Prospect of Calculation Method for Energy Flow in Integrated Energy System. Automation Electr. Power Syst. 44, 179–191. doi:10.7500/AEPS20191113001

CrossRef Full Text | Google Scholar

Jia, H. J., Wang, D., Xu, X. D., and Yu, X. D. (2014). Research on Some Key Problems Related to Integrated Energy Systems. Automation Electr. Power Syst. 39, 199–207. doi:10.7500/AEPS20141009011

CrossRef Full Text | Google Scholar

Jiang, P., and Liang, L. (2015). Reactive Power Optimization of Hybrid AC/HVDC Power System Combining Interior Point Algorithm and Genetic Algorithm. High. Volt. Eng. 41, 724–729. doi:10.13336/j.1003-6520.hve.2015.03.004

CrossRef Full Text | Google Scholar

Krause, T., Andersson, G., Fröhlich, K., and Vaccaro, A. (2011). Multiple-energy Carriers: Modeling of Production, Delivery, and Consumption. Proc. IEEE 99, 15–27. doi:10.1109/JPROC.2010.2083610

CrossRef Full Text | Google Scholar

Li, J. L., Tian, L., and Lai, X. K. (2015). Outlook of Electrical Energy Storage Technologies under Energy Internet Background. Automation Electr. Power Syst. 39, 15–25. doi:10.7500/AEPS20150906004

CrossRef Full Text | Google Scholar

Liu, H., Zhao, C. X., Ge, S. Y., Li, J. F., and Liu, J. Y. (2021). Sequential Power Flow Calculation of Power-Heat Integrated Energy System. Automation Electr. Power Syst. 45, 63–72. doi:10.7500/AEPS20200312009

CrossRef Full Text | 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

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 [Guest Editorial]. IEEE Power Energy Mag. 11, 18–23. doi:10.1109/MPE.2013.2266594

CrossRef Full Text | Google Scholar

Rui, W., Qiuye, S., Pinjia, Z., Yonghao, G., Dehao, Q., and Peng, W. (2020). Reduced-order Transfer Function Model of the Droop-Controlled Inverter via Jordan Continued-Fraction Expansion. IEEE Trans. Energy Convers. 35, 1585–1595. doi:10.1109/TEC.2020.2980033

CrossRef Full Text | Google Scholar

Sun, G. Q., Wang, W. X., Wu, Y., Hu, W., and Jing, J. P. (2020). Fast Power Flow Calculation Method for Radiant Electric-Thermal Interconnected Integrated Energy System. Proc. CSEE 40, 4131–4141. doi:10.13334/j.0258-8013.pcsee.190583

CrossRef Full Text | Google Scholar

Wang, R., Sun, Q., Hu, W., Li, Y., Ma, D., and Wang, P. (2021a). SoC-based Droop Coefficients Stability Region Analysis of the Battery for Stand-Alone Supply Systems with Constant Power Loads. IEEE Trans. Power Electron. 36, 7866–7879. doi:10.1109/TPEL.2021.3049241

CrossRef Full Text | Google Scholar

Wang, R., Sun, Q., Ma, D., and Liu, Z. (2019). The Small-Signal Stability Analysis of the Droop-Controlled Converter in Electromagnetic Timescale. IEEE Trans. Sustain. Energy 10, 1459–1469. doi:10.1109/TSTE.2019.2894633

CrossRef Full Text | Google Scholar

Wang, R., Sun, Q., Tu, P., Xiao, J., Gui, Y., and Wang, P. (2021b). Reduced-order Aggregate Model for Large-Scale Converters with Inhomogeneous Initial Conditions in DC Microgrids. IEEE Trans. Energy Convers. 36, 2473–2484. doi:10.1109/TEC.2021.3050434

CrossRef Full Text | Google Scholar

Xia, T., Chen, Y. W., Guo, Q. L., Sun, H. B., Xiong, W., and Liu, Y. Q. (2017). Power Flow Calculation Based on PSASP for Combined Thermal and Electric Energy System. Electr. Power Autom. Equip. 37, 55–61. doi:10.16081/j.issn.1006-6047.2017.06.008

CrossRef Full Text | Google Scholar

Yu, J., Yang, D., and Chen, Z. (2021). Multi-energy Flow Calculation Based on Energy Cell and Parallel Distributed Computation. Int. J. Electr. Power & Energy Syst. 131, 107147. doi:10.1016/j.ijepes.2021.107147

CrossRef Full Text | Google Scholar

Zhang, T., Li, Z., Wu, Q., Pan, S., and Wu, Q. H. (2022). Dynamic Energy Flow Analysis of Integrated Gas and Electricity Systems Using the Holomorphic Embedding Method. Appl. Energy 309, 118345. doi:10.1016/j.apenergy.2021.118345

CrossRef Full Text | Google Scholar

Zhang, Y. Z., Wang, X. J., He, J. H., and Gao, W. Z. (2019). Optimal Energy Flow Calculation Method of Integrated Energy System Considering Thermal System Modeling. Trans. China Electrotech. Soc. 34, 563–570. doi:10.19595/j.cnki.1000-6753.tces.L80501

CrossRef Full Text | Google Scholar