A Fast Calculation Method for Economic Dispatch of Electro-Thermal Coupling System Considering the Dynamic Process of Heat Transfer

Abstract

The dynamic spatial and temporal characteristics of heat transfer within heating network pipelines are important factors affecting the accuracy of economic dispatch decision-making results of electro-thermal coupling systems. However, the pipeline heat transmission process is described by partial differential equations, which makes it difficult to solve quickly. Therefore, this study introduces a model for calculating the economic dispatch of the electro-thermal coupling system (EDETCS) that takes into account the pipeline transmission process. Firstly, based on the implicit upwind difference method, a two-port model of branch heat transfer dynamics is established. Secondly, the general term formula of the two-port model is derived. Finally, the established two-port model is applied to the EDETCS. The findings from the example analysis indicate that, in contrast to the conventional calculation method, the proposed model improves the calculation speed while ensuring the accuracy of the solution.

1. Introduction

Compared with thermal plants, combined heat and power (CHP) plants have more advantages in improving energy efficiency, energy saving and emission reduction [1]. However, under the new energy situation, CHP units operated in the traditional “determining electricity by heat” mode can easily reach the minimum rigid technical output of the unit, which seriously restricts the flexible supply of the power grid [2,3]. By considering the pipeline transmission process in the EDETCS, we can effectively manage renewable energy consumption and minimize system operating costs. This is achieved by leveraging the complementary nature of power and heat demand alongside the ability to store the heat of the heating network [4,5]. However, the pipeline heat transmission process is described by partial differential equations, and it is very difficult to solve such optimization problems with partial differential equation constraints [6,7,8,9].

At present, researchers have conducted some research on the pipeline transmission process. In references [10,11], considering the heat loss effect and ignoring the heat delay effect, the linear relationship between the inlet and outlet water temperature of the pipeline is obtained. Reference [12] deduced the partial differential equation of heat transmission in the pipeline through the law of thermodynamics and solved it by the finite difference method. However, its accuracy depends on a smaller time and space step; the calculation amount is large, and the calculation time is long. References [13,14] proposed the node method to analyze and model the heating network. The principle is to track the time when the water body starts from the previous node to the next node and calculate the temperature at the head of the pipe according to the historical temperature of different nodes. On this basis, the temperature at the end of the pipe is calculated. The node method reduces the computing time, but large errors will occur when the selected time resolution is too large [15].

In addition to the research ideas of the above-simplified mechanism model, the existing research work has begun to try the analogy circuit method to establish the equivalent model of heat networks. In [16], an equivalent model of a heat network based on the Fourier transform is proposed, and the partial differential dynamic model of heat transfer is mapped to the frequency domain for solution. Reference [17] proposed a heat network equivalent model based on Laplace transform. This model transforms the partial differential equation into the s-domain to solve and realizes the simplification of the model to the algebraic equation. Based on the thermoelectric analogy method, Reference [18] introduced the power system analysis method to equivalent the loss, energy storage and temperature influence of each temperature micro-element in the heating network pipeline through components such as ‘heat resistance’, ‘heat capacity’ and ‘heat source’. Finally, the functional relationship between the outlet temperature and the initial temperature and time is obtained by mathematical derivation. References [19,20] proposed the unified energy path theory. The physical mechanisms involved in gas and heating networks are similar to those in the power transmission system. In the frequency domain, a unified set of mathematical equations for power, heat, and gas networks has been developed. Finally, the linear two-port lumped parameter model in the frequency domain is obtained. In references [21,22], the distributed parameter circuit model of electric, heat and gas energy flow is established by the Laplace transform, and then the lumped parameter transmission model with the branch as a unit is obtained. Although the above method transforms the time-domain dynamic model of heat transfer described by partial differential equations into a linear two-port model in the frequency domain or s-domain, it reduces the difficulty of solving the dynamic model of the pipeline transmission process. However, it is necessary to carry out a transformation on the actual known initial conditions and boundary conditions before solving, and after the calculation is completed, the calculation results need to be inversely transformed back to the time domain, and the transformed link reduces the solution efficiency.

On the basis of the above research, considering the timeliness of economic dispatch, this paper establishes a two-port model of the pipeline heat transmission process, aiming to reduce the calculation scale of the EDETCS without reducing the accuracy. Firstly, based on the implicit upwind difference structure, a two-port model of heat transfer dynamics in the heating network is established to reduce the calculation scale. Secondly, the general term formula of the two-port model is derived, which makes the model more convenient and more operable for the modeling of complex heat networks. Then, the proposed two-port model is embedded in the EDETCS. Finally, the example analysis demonstrates the effectiveness and advantages of the two-port model.

2. Two-Port Model of Heat Transfer in Pipes

2.1. Dynamic Model of Heat Transfer Based on Implicit Up-Wind Difference Scheme

China’s heating network is generally composed of a two-level network. The two exchange heat through the heat exchange station. Because the secondary pipe is very short, we do not consider the secondary heating network; only the primary heating network is modeled [23].

The pipeline transmission process in the heat pipe network can be described by Equation (1) [24]:

$${\displaystyle \frac{\partial T_{\mathrm{w}}}{\partial t}}+{\displaystyle \frac{m}{A\rho }}{\displaystyle \frac{\partial T_{\mathrm{w}}}{\partial x}}+{\displaystyle \frac{\lambda }{A\rho c}}(T_{\mathrm{w}}-T_{\mathrm{a}})={\displaystyle \frac{\lambda }{\rho c}}{\displaystyle \frac{{\partial }^2{T}_{\mathrm{w}}}{\partial x^2}}$$

(1)

where m is the mass flow, and when the system adopts mass regulation mode, m is constant. $\rho$ is the water density. $A$ is the pipe’s sectional area, $c$ is the water-specific heat capacity, λ is the thermal conductivity, T_w is water temperature, and T_a is ambient temperature.

Compared with thermal convection, the heat conduction effect between pipe fluids is very weak. Therefore, the heat conduction term in Equation (1), that is, the second-order partial derivative term, can be ignored without obvious error. The pipeline transmission process in the heating network can be written as [17]:

$${\displaystyle \frac{\partial T}{\partial t}}+{\displaystyle \frac{m}{A\rho }}{\displaystyle \frac{\partial T}{\partial x}}+{\displaystyle \frac{\lambda }{A\rho c}}T=0$$

(2)

where T is the difference between T_w and T_a.

The heating system in China mainly adopts the mass regulation mode. Therefore, this paper will derive the model of heat transfer under the premise of known mass flow.

Firstly, the first-order implicit upwind difference quotient [12] is used to approximate the first-order partial derivative and the difference grid is shown in Figure 1:

$${\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{T_i^n-T_i^{n-1}}{\tau }}$$

(3)

$${\displaystyle \frac{\partial T}{\partial x}}={\displaystyle \frac{T_i^n-T_{i-1}^n}{h}}$$

(4)

$$I={\displaystyle \frac{L}{h}},N={\displaystyle \frac{P}{\tau }}$$

(5)

where τ is the time step, h is the spatial step, i is the number of spatial nodes, and n is the number of time nodes. $T_i^n$ are the water temperature at node (n, i).

Combining Equation (3), Equation (4) and heat transfer Equation (2), the difference scheme of pipeline temperature is obtained:

$$T_i^n=K_1{T}_i^{n-1}+K_2{T}_{i-1}^n$$

(6)

where K₁ and K₂ are two parameters defined to simplify the representation:

$$K_1={\displaystyle \frac{A\rho ch}{A\rho ch+h\tau \lambda +cm\tau }}$$

(7)

$$K_2={\displaystyle \frac{cm\tau }{A\rho ch+h\tau \lambda +cm\tau }}$$

(8)

Then, the boundary conditions and initial conditions are determined:

Usually, the initial temperature and the inlet water temperature are known, and this paper assumes that the temperature of each node at the time t = 0 is equal, so the boundary conditions and initial conditions can be determined as follows:

$$T_x^0=T_0^0 , T_0^t=T^{\mathrm{in}}$$

(9)

where: $T_x^0$ is the initial temperature of each node, $T_0^0$ is the initial temperature at the head of the pipe, $T_0^t$ is the inlet water temperature at the head of the pipe, and $T^{\mathrm{in}}$ is the temperature of water at the inlet of the pipe changing with time.

2.2. Two-Port Model of Branch Heat Transfer Dynamics

According to the initial conditions, boundary conditions Equation (9) and the derived linear model Equation (6), the temperature curve at the end of the pipe can be iteratively solved. However, this finite difference method will calculate the temperature of each node for iteration during the calculation process. When selecting a smaller time and space step, the computation involves a substantial quantity of data and takes a significant amount of time. To enhance computational efficiency, this section will establish a linear correlation between the temperatures of various nodes at different timestamps within the difference grid and the boundary conditions and initial conditions according to the difference equation of temperature and finally construct the two-port model of branch heat transfer dynamics to reduce the scale of model calculation.

When t = 1, the difference grid of the pipe temperature is shown in Figure 2. According to Equation (6), the water temperature at node i $(1\le i\le I)$ is:

$$T_i^1=K_1{T}_i^0+K_2{T}_{i-1}^1$$

(10)

Substituting Equation (9) into Equation (10):

$$T_i^1=K_1{T}_i^0+K_2{T}_{i-1}^1$$

(11)

The temperature at node (1, i) is obtained:

$$\begin{array}{ll}{T}_i^1& =T_0^0{\displaystyle \sum _{b=1}^i{K}_1{(K_2)}^{b-1}}+T_0^1{(K_2)}^i\\ & =a_{i,0}^1{T}_0^0+a_{i,1}^1{T}_0^1+\cdots +a_{i,N}^1{T}_0^N\end{array}$$

(12)

$$\Rightarrow T_i^1=\left[a_{i,0}^1,a_{i,1}^1,\cdots a_{i,N}^1\right]{\mathit{T}}^{\mathrm{in}}$$

(13)

where:

$$a_{i,0}^1={\displaystyle \sum _{b=1}^i{K}_1{(K_2)}^{b-1}}$$

(14)

$$a_{i,1}^1={(K_2)}^i$$

(15)

$$a_{i,2}^1=\cdots =a_{i,N}^1=0$$

(16)

$${\mathit{T}}^{\mathrm{in}}={\left[T_0^0,T_0^1,\cdots ,T_0^N\right]}^T$$

(17)

When t = 2, the difference grid of the pipe temperature is shown in Figure 3. According to Equation (6), the water temperature at node I $(1\le i\le I)$ is:

$$T_i^2=K_1{T}_i^1+K_2{T}_{i-1}^2$$

(18)

Substituting Equation (9) into Equation (18), the temperature of each node is calculated in turn, and the temperature at node (2, i) is obtained:

$$\begin{array}{lll}{T}_i^2& =& T_0^0{(K_1)}^2{\displaystyle \sum _{b=1}^{i}b{(K_2)}^{b-1}}+\\ & & T_0^{1}i{K}_1{(K_2)}^i+T_0^2{(K_2)}^i\\ & =& a_{i,0}^2{T}_0^0+a_{i,1}^2{T}_0^1+\cdots +a_{i,N}^2{T}_0^N\end{array}$$

(19)

$$\Rightarrow T_i^2=\left[a_{i,0}^2,a_{i,1}^2,\cdots a_{i,N}^2\right]{\mathit{T}}^{\mathrm{in}}$$

(20)

where:

$$a_{i,0}^2={(K_1)}^2{\displaystyle \sum _{b=1}^{i}b{(K_2)}^{b-1}}$$

(21)

$$a_{i,1}^2=i{K}_1{(K_2)}^i$$

(22)

$$a_{i,2}^2={(K_2)}^i$$

(23)

$$a_{i,3}^2=a_{i,4}^2=\cdots =a_{i,N}^2=0$$

(24)

Through the above method, it can be obtained that when t = 3, the temperature at node i $(1\le i\le I)$ is:

$$\begin{array}{lll}{T}_i^3& =& T_0^0{(K_1)}^3{\displaystyle \sum _{b=1}^i\frac{b(b+1)}{2}{(K_2)}^{b-1}}+\\ & & T_0^1{\displaystyle \frac{i(i+1)}{2}}{(K_1)}^2{(K_2)}^i+\\ & & T_0^{2}i{K}_1{(K_2)}^i+T_0^3{(K_2)}^i\\ & =& a_{i,0}^3{T}_0^0+a_{i,1}^3{T}_0^1+\cdots +a_{i,N}^3{T}_0^N\end{array}$$

(25)

$$\Rightarrow T_i^3=\left[a_{i,0}^3,a_{i,1}^3,\cdots a_{i,N}^3\right]{\mathit{T}}^{\mathrm{in}}$$

(26)

where:

$$a_{i,0}^3={(K_1)}^3{\displaystyle \sum _{b=1}^i\frac{b(b+1)}{2}{(K_2)}^{b-1}}$$

(27)

$$a_{i,1}^3={\displaystyle \frac{i(i+1)}{2}}{(K_1)}^2{(K_2)}^i$$

(28)

$$a_{i,2}^3=i{K}_1{(K_2)}^i$$

(29)

$$a_{i,3}^3={(K_2)}^i$$

(30)

$$a_{i,4}^3=a_{i,5}^3=\cdots =a_{i,N}^3=0$$

(31)

Through the temperature equation at each moment, the matrix form of the temperature at node (n, i) can be obtained:

$$T_i^n={\mathit{A}}_i^n{\mathit{T}}^{\mathrm{in}}$$

(32)

where ${\mathit{A}}_i^n=\left[a_{i,0}^n,a_{i,1}^n,\cdots a_{i,N}^n\right]$.

By merging the temperature expressions at each time point, the temperature at each time point at node i is obtained, that is, the two-port model of branch heat transfer:

$${\mathit{T}}_i={\mathit{A}}_i{\mathit{T}}^{\mathrm{in}}$$

(33)

$${\mathit{T}}_i={\left[T_i^1,T_i^2,\cdots ,T_i^N\right]}^T$$

(34)

$${\mathit{A}}_i=\left[\begin{array}{cccc}{a}_{i,0}^1& a_{i,1}^1& \cdots & a_{i,N}^1\\ a_{i,0}^2& a_{i,1}^2& \cdots & a_{i,N}^1\\ ⋮& ⋮& \ddots & ⋮\\ a_{i,0}^N& a_{i,1}^N& \cdots & a_{i,N}^1\end{array}\right]$$

(35)

2.3. The General Term Formula of Coefficient Matrix of Two-Port Model

In Section 2.2, a two-port model of branch heat transfer dynamics is established. However, various pipelines exhibit differences in parameters like diameter and length. Therefore, the coefficient matrix of the model needs to be re-derived, especially in large-scale heating systems, which will greatly increase the calculation workload. For this reason, this section will derive the general term formula of the coefficient matrix to reduce the calculation workload.

Summarize the coefficients of the differential form of water temperature at each time. We find that the elements of the coefficient matrix of the temperature expression at node (n, i) are the linear combination of the Pascal matrix and the parameters K₁, K₂.

Let M = max{I,N}, let M-order diagonal matrix Q(i) = diag [1,1,…,0], and the first i diagonal elements are 1. Let R(n) = [0,…,0,1,0,…,0] be an M-order row vector, and its nth element is 1. Let M-order column vector K = [1,K₂,(K₂)²,…,(K₂)^M−1]^T. Let P(M) be an M-order Pascal matrix. Then P(M)Q(i) means that the first i column elements of P(M) are retained, and other elements are set to zero, R(n)P(M)Q(i) represents the nth row element of the matrix and [P(M)Q(i)]; R(n)P(M)Q(i)R(i) represents the nth row i column element of the matrix [P(M)Q(i)].

Then, The formula for the general term of the coefficient matrix ${\mathit{A}}_i^n$ can be articulated as follows:

$$a_{i,b}^n=\left\{\begin{array}{l}{(K_1)}^n\mathit{R}(n)\mathit{P}(M)\mathit{Q}(i)\mathit{K}\hspace{1em}\hspace{1em}\hspace{1em} b=0\\ {(K_1)}^{n-b}{(K_2)}^i\mathit{R}(n-b+1)\mathit{P}(M)Q(i)R(i)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} 0<b\le N\end{array}\right.$$

(36)

By merging the coefficient matrices at each moment, the general term formula of the coefficient matrix A_i of the model can be obtained.

The establishment of the general term formula makes it unnecessary to derive the temperature expression of each moment and each node in the initial modeling and only needs to assign the coefficient matrix A_i through the loop statement. Moreover, for different pipelines, the elements of the model coefficient matrix A_i are only different in parameters K₁ and K₂. It is only necessary to change the parameters and assign them again without re-derivation, which greatly reduces the calculation scale and makes the model establishment of large-scale complex heat networks easier.

2.4. Two-Port Model of Heat Transfer Dynamics in Heating Network Pipeline

According to the previous derivation, when the initial conditions and boundary conditions are the initial temperature of each node of the pipeline and the inlet water temperature, respectively, the outlet temperature of the pipeline is:

$${\mathit{T}}^{\mathrm{out}}={\mathit{A}}_I{\mathit{T}}^{\mathrm{in}}$$

(37)

$${\mathit{T}}^{\mathrm{out}}={\left[T_I^1,T_I^2,\cdots ,T_I^N\right]}^T$$

(38)

$${\mathit{A}}_I=\left[\begin{array}{cccc}{a}_{I,0}^1& a_{I,1}^1& \cdots & a_{I,N}^1\\ a_{I,0}^2& a_{I,1}^2& \cdots & a_{I,N}^1\\ ⋮& ⋮& \ddots & ⋮\\ a_{I,0}^N& a_{I,1}^N& \cdots & a_{I,N}^1\end{array}\right]$$

(39)

The above two-port model eliminates the intermediate nodes between the source and the load in the calculation process and further realizes the dimension reduction modeling and operation of the model. The coefficient matrix can be calculated offline without online iteration.

3. Formulation of EDETCS

3.1. Objective Function

The objective of the dispatch model is to reduce the system’s operating cost to a minimum:

$$F_{\min }=C_{\mathrm{G}}+C_{\mathrm{chp}}+C_{\mathrm{w}}$$

(40)

where $C_{\mathrm{G}}$ is the operating cost of thermal units, $C_{\mathrm{chp}}$ is the operation cost of CHP units, and $C_{\mathrm{w}}$ is the wind curtailment cost.

The three costs are:

$$C_{\mathrm{G}}={\displaystyle \sum _{t=1}^N{\displaystyle \sum _{j=1}^{N^{\mathrm{G}}}(a_j^{\mathrm{G}}{P}_{j,t}^{\mathrm{G}}+b_j^{\mathrm{G}})}}$$

(41)

$$\begin{array}{ll}{C}_{\mathrm{chp}}=& {\displaystyle \sum _{t=1}^N{\displaystyle \sum _{j=1}^{N^{\mathrm{chp}}}[}}{a}_j^{\mathrm{chp}}{(P_{j,t}^{\mathrm{chp}})}^2+b_j^{\mathrm{chp}}{P}_{j,t}^{\mathrm{chp}}+\\ & c_j^{\mathrm{chp}}+d_j^{\mathrm{chp}}{(H_{j,t}^{\mathrm{chp}})}^2+e_j^{\mathrm{chp}}{H}_{j,t}^{\mathrm{chp}}+\\ & f_j^{\mathrm{chp}}{H}_{j,t}^{\mathrm{chp}}{P}_{j,t}^{\mathrm{chp}}]\end{array}$$

(42)

$$C_{\mathrm{w}}={\displaystyle \sum _{t=1}^N{\displaystyle \sum _{j=1}^{N^{\mathrm{W}}}{\lambda }_j^s(P_{j,t}^{\mathrm{w},\mathrm{pre}}-P_{j,t}^{\mathrm{w}})}}$$

(43)

where $N^{\mathrm{G}}$, $N^{\mathrm{chp}}$, $N^{\mathrm{W}}$ are the set of units’ number. $a_j^{\mathrm{G}}$, $b_j^{\mathrm{G}}$ are the cost coefficients of thermal units. $P_{j,t}^{\mathrm{chp}}$, $H_{j,t}^{\mathrm{chp}}$ are the CHP’s power and heat output. $a_j^{\mathrm{chp}}$, $b_j^{\mathrm{chp}}$, $c_j^{\mathrm{chp}}$, $d_j^{\mathrm{chp}}$, $e_j^{\mathrm{chp}}$, $f_j^{\mathrm{chp}}$ are the correlation coefficients of the operating costs of CHP units. ${\lambda }_j^s$ is the wind curtailment cost coefficient. $P_{j,t}^{\mathrm{w},\mathrm{pre}}$ is the predicted maximum output of wind power. $P_{j,t}^{\mathrm{w}}$ is the actual output of wind power.

3.2. Constraints

(1) Heat transfer dynamics constraints:

The end temperature of the pipe k satisfies:

$${\mathit{T}}_k^{\mathrm{out}}={\mathit{A}}_{I(k)}{\mathit{T}}_k^{\mathrm{in}}$$

(44)

where ${\mathit{T}}_k^{\mathrm{out}}$ is the outlet temperature matrix. k represents the pipe number.

(2) Temperature mixing constraints:

$${\displaystyle \sum _{k\in S_n^{\mathrm{pipe}-}}{m}_k{T}_{k,t}^{\mathrm{out}}=T_{l,t}^{\mathrm{in}}}{\displaystyle \sum _{l\in S_n^{\mathrm{pipe}+}}{m}_l}$$

(45)

where $S_n^{\mathrm{pipe}+}$, $S_n^{\mathrm{pipe}-}$ are the set of pipes starting at node n and ending at node n.

(3) Heat exchange station constraints:

$$H_{j,t}^{\mathrm{LN}}=c{m}_j^{\mathrm{LN}}(T_{j,t}^{\mathrm{LN},\mathrm{s}}-T_{j,t}^{\mathrm{LN},\mathrm{r}})$$

(46)

$${\underset{¯}{T}}_j^{\mathrm{LN},\mathrm{r}}\le T_{j,t}^{\mathrm{LN},\mathrm{r}}\le {\overline{T}}_j^{\mathrm{LN},\mathrm{r}}$$

(47)

$${\underset{¯}{T}}_j^{\mathrm{LN},\mathrm{s}}\le T_{j,t}^{\mathrm{LN},\mathrm{s}}\le {\overline{T}}_{j,t}^{\mathrm{LN},\mathrm{s}}$$

(48)

where $m_j^{\mathrm{LN}}$ is the mass flow at the entrance of heat substations, $H_{j,t}^{\mathrm{LN}}$ is the heat consumed by the heat substations. $T_{j,t}^{\mathrm{LN},\mathrm{s}}$ is the supply and $T_{j,t}^{\mathrm{LN},\mathrm{r}}$ is the return temperature from heat substations. ${\overline{T}}_j^{\mathrm{LN},\mathrm{r}}$ is the upper and ${\underset{¯}{T}}_j^{\mathrm{LN},\mathrm{r}}$ is the lower limits of the return temperature, ${\overline{T}}_{j,t}^{\mathrm{LN},\mathrm{s}}$ is the upper and ${\underset{¯}{T}}_j^{\mathrm{LN},\mathrm{s}}$ is the lower limits of the supply temperature.

(4) CHP unit constraints:

$$H_{j,t}^{\mathrm{chp}}=c{m}_j^{\mathrm{chp}}(T_{j,t}^{\mathrm{chp},\mathrm{s}}-T_{j,t}^{\mathrm{chp},\mathrm{r}})$$

(49)

where $H_{j,t}^{\mathrm{chp}}$ is the CHP heat output. $m_j^{\mathrm{chp}}$ is the mass flow at the entrance of CHP units, $T_{j,t}^{\mathrm{chp},\mathrm{r}}$ is the return and $T_{j,t}^{\mathrm{chp},\mathrm{s}}$ is the supply temperature from the heat source.

(5) power balance constraints:

$${\displaystyle \sum _{j=1}^{N^{\mathrm{G}}}{P}_{j,t}^{\mathrm{G}}}+{\displaystyle \sum _{j=1}^{N^{\mathrm{chp}}}{P}_{j,t}^{\mathrm{chp}}}+{\displaystyle \sum _{j=1}^{N^{\mathrm{W}}}{P}_{j,t}^{\mathrm{w}}}=P_t^{\mathrm{LN}}$$

(50)

where $P_t^{\mathrm{LN}}$ is the active load.

(6) CHP unit constraints:

$${\underset{\_}{P}}_j^{\mathrm{chp}}-c_j^{\mathrm{v}}{H}_{j,t}^{\mathrm{chp}}\le P_{j,t}^{\mathrm{chp}}\le {\overline{P}}_j^{\mathrm{chp}}-c_j^{\mathrm{c}}{H}_{j,t}^{\mathrm{chp}}$$

(51)

$$P_j^{\mathrm{co}}+c_j^{\mathrm{m}}{H}_{j,t}^{\mathrm{chp}}\le P_{j,t}^{\mathrm{chp}}\le {\overline{P}}_j^{\mathrm{chp}}-c_j^{\mathrm{c}}{H}_{j,t}^{\mathrm{chp}}$$

(52)

$${\underset{¯}{H}}_j^{\mathrm{chp}}\le H_{j,t}^{\mathrm{chp}}\le {\overline{H}}_j^{\mathrm{chp}}$$

(53)

where ${\overline{P}}_j^{\mathrm{chp}}$ is the maximum and ${\underset{\_}{P}}_j^{\mathrm{chp}}$ is the minimum power output of CHP units. $c_j^{\mathrm{v}}$, $c_j^{\mathrm{c}}$ $c_j^{\mathrm{m}}$ $P_j^{\mathrm{co}}$ are the feasible region boundary coefficient of CHP units. ${\overline{H}}_j^{\mathrm{chp}}$ is the upper and ${\underset{¯}{H}}_j^{\mathrm{chp}}$ is the lower limits of the heat power of CHP units.

(7) Ramping constraints:

$${\underset{¯}{D}}_j^{\mathrm{G}}\le P_{j,t}^{\mathrm{G}}-P_{j,t-1}^{\mathrm{G}}\le {\overline{D}}_j^{\mathrm{G}}$$

(54)

where ${\overline{D}}_j^{\mathrm{G}}$ is upward and ${\underset{¯}{D}}_j^{\mathrm{G}}$ is the downward ramping capability of thermal units.

(8) Thermal unit constraint:

$${\underset{¯}{P}}_j^{\mathrm{G}}\le P_{j,t}^{\mathrm{G}}\le {\overline{P}}_j^{\mathrm{G}}$$

(55)

where ${\overline{P}}_j^{\mathrm{G}}$ is the upper and ${\underset{¯}{P}}_j^{\mathrm{G}}$ is the lower limits of the thermal unit output.

(9) Wind power output constraints:

$$0\le P_{j,t}^{\mathrm{w}}\le P_{j,t}^{\mathrm{w},\mathrm{pre}}$$

(56)

4. Case Study

4.1. Double Pipe Heating System

A heating system consisting of a supply pipe, a return pipe, a heat source, and a heat substation has been chosen to validate the proposed two-port model. Figure 4 shows the system topology diagram. The accuracy and superiority of the two-port model are verified by comparing the calculation results with the implicit upwind difference method [12].

Table 1. Parameters of Double pipe heating pipeline.

Parameter	A/m2	λ/(W/(m∙°C)	m/(kg/s)	L/km	Ta/°C
Value	0.05	0.25	36	5	0

lists the pipe parameters. The system works in the quality regulation mode, and the heat load is constant power operation. Set the time step τ = 10 s, space step h = 10 m; heat load is 5 MW.

Figure 5 illustrates the computed water temperature profiles. The water temperature of the water supply pipeline tends to be stable and close to the heating temperature curve after 115 min of cooling. This is consistent with the time that the hot water flows through 5000 m at a mass flow rate of 36 kg/s, reflecting the “heat loss” and “delay” phenomena in the pipeline heat transmission process. This shows that the two-port model of heat transfer dynamics can truly reflect the pipeline transmission process in the heat network pipeline. Moreover, compared with the calculation results of the difference method, the temperature curves obtained by the two calculation methods are almost completely coincident, and the mean relative error of water supply pipe outlet temperature and return pipe outlet temperature are 0.15% and 0.23%, respectively, which further verifies the accuracy of describing the pipeline transmission process.

The duration of the calculations is presented in

Table 2. Time of calculation.

Model	Time/s
The two-port model	0.012
The implicit upwind difference method	0.22

. The solution time of the two-port model is only 0.012 s, while the calculation time of the difference method is 0.22 s, and the speed is increased by more than 18 times, which verifies the superiority of the two-port model in calculation speed.

4.2. District Heating Network (DHN)

The DHN shown in Figure 6 is used to illustrate the two-port model’s effectiveness, and the calculation time is compared with the different method used in Reference [12]. The heating network comprises three water supply lines, three return lines, a heat source, and three heat substations. The system works in the quality regulation mode, and the heat load is constant power operation.

Table 3. Pipe parameters.

Pipe	A/m2	λ/(W/(m∙°C)	m/(kg/s)	L/km	Ta/°C
1, 6	0.05	0.25	36	5	0
2, 5	0.05	0.25	14	3.5	0
3, 4	0.05	0.25	7	2.5	0

lists the pipe parameters. The heat load is set to 2.5 MW, 1 MW and 1 MW, respectively.

Figure 7 and Figure 8 show the calculated water temperature curves of each node under the two methods. The results show that the temperature curves of the two models are almost the same, and the maximum relative errors of water supply pipe end temperature and return pipe end temperature are 0.23% and 0.41%, respectively, which indicates that the precision of both methods in calculations is comparable.

Table 4. Time of calculation.

Model	Time/s
The two-port model	0.057
The implicit upwind difference method	2.68

shows the time it takes to calculate the two models. We can find that the calculation time of the two-port model is much lower than that of the difference method. Compared with the calculation time in 4.1, the calculation time of the difference method is increased by more than 10 times, while the calculation time of the two-port model is only increased by less than 5 times. This is because, in systems with more heating network pipelines, the finite difference method will introduce state variables that are proportional to the number of spatial nodes. The two-port model only needs to introduce variables that are proportional to the number of pipelines, which reduces the total number of variables and reduces the calculation scale.

The above analysis shows that the two-port model can greatly reduce the solution time, and it has more significant advantages in the calculation speed when the heating system structure is more complex.

4.3. EDETCS

In this section, two models are applied to the EDETCS, considering the pipeline transmission process.

Model 1: A steady-state model of heating network used in Reference [10];

Model 2: The two-port model of heat transfer dynamics.

Figure 9 depicts the layout of the testing system, which includes two thermal units, G1 and G2, alongside a wind farm designated as W. The CHP unit serves as the system’s heat source. Relevant system details can be found in

Table 5. Thermal power unit parameters.

Unit	Pmin/MW	Pmax/MW	Cost Coefficients
			aG/Yuan	bG/(Yuan/MW)
G1	100	220	13.5	177
G2	10	100	40	130

and

Table 6. CHP unit parameters.

Output of CHP	Pmin/ MW	Pmax/ MW	Cost Coefficients
			achp/(Yuan/MW2), bchp/(Yuan/MW), cchp/Yuan, dchp/(Yuan/MW2), echp/(Yuan/MW), fchp/(Yuan/MW2),
electric power	45	125	0.0032, 17.7, 181, 0.00085, 4.2, 0.00125

, while

Table 7. Pipe parameters.

Node	Node	A/m2	λ/(W/(m∙°C)	m/(kg/s)	L/km	Ta/°C
1	2	0.5	0.25	360	3.5	0
2	3	0.5	0.25	150	1.75	0
3	4	0.5	0.25	100	1.75	0
2	5	0.5	0.25	210	1.75	0
3	6	0.5	0.25	50	0.75	0

provides information on pipeline specifications. Figure 10 displays the electric load, heat load, and projected peak output from wind power.

The operating costs and wind curtailment costs of the two models are shown in

Table 8. Operation cost.

Model	Operation Cost/Yuan	Wind Curtailment Cost/Yuan
Model 1	101,596	10,303
Model 2	97,408	4379

. The operating costs of model 1 and model 2 are 97,408 yuan and 10,1596 yuan, respectively. Compared with model 1, the system applying model 2 can reduce the cost by 4188 yuan.

Figure 11 shows the wind power scheduling decisions under the two models. The results show that the wind curtailment phenomenon occurs in model 1 at 1:00–3:00 and 5:00–6:00. The reason is that the CHP heat power corresponds to the heat load and cannot be flexibly adjusted. Compared with model 1, model 2 consumes more wind power. The reason is that model 2 does not need to follow the heat load and can reduce the CHP power of the unit at 1:00–3:00 and 5:00–6:00 to absorb wind power.

Figure 12 and Figure 13 show the system’s thermal scheduling and power scheduling results using the two models. In model 1, the results show that the CHP heat output changes with the heat load. As a result, the CHP power dispatch is also limited by the heat output. However, in model 2, CHP heat output and heat load do not correspond, and the two are decoupled so that the system can make full use of CHP units for economic optimization and increase wind power consumption.

Since the generation cost of G2 is the highest, the G2 scheduling strategies of the two models are maintained at the minimum output of 10 MW. Figure 14 shows the scheduling decision of thermal power unit G1. The results show that compared with model 1, the thermal power unit G1 of model 2 increases its output at 6:00–14:00 and 21:00–24:00 because G1 has the lowest power generation cost. This shows that the model considering the pipeline transmission process can enhance the operational flexibility and economy of the system operation.

5. Conclusions

To enhance the efficiency of addressing the power optimization challenges involving partial differential equation constraints—like the EDETCS, which factors in the dynamic aspects of heat transfer—this paper presents a two-port heat transfer model utilizing implicit upwind differences. The key findings are as follows:

(1) Compared with the implicit upwind difference method, the heat transfer two-port model can greatly reduce the computational model complexity and improve the computational efficiency while ensuring accuracy.

(2) Implementing this two-port model within the EDETCS allows for precise representation of the pipeline heat transmission process, effectively leverages the heat storage capabilities of the network, boosts wind energy utilization, and lowers the system’s operational costs.

In the future, we will try to apply the two-port model to the simulation analysis and real-time control of the actual electro-heat integrated energy system.