Integrated Heat and Electricity Dispatch for District Heating Networks With Constant Mass Flow: A Generalized Phasor Method

Using thermal inertia in district heating systems (DHSs) to improve the dispatch flexibility and economy of integrated heat and electricity systems (IHESs) is a research hotspot and difficulty. In most existing studies, the partial differential equations (PDEs) of thermal inertia are approximated by discrete-time models, making it difficult to accurately describe the continuous dynamic processes. In this paper, we propose a novel generalized phasor method (GPM) for thermal inertia in DHSs with constant mass flow. Based on the analytical solution of the PDEs and the Fourier transform, the intractable PDEs are transformed into a series of complex algebraic equations represented by phasors. The GPM has higher accuracy compared to traditional discrete models because it is essentially a continuous model in the time domain. Then, we present a different representation of an integrated heat and electricity dispatch (IHED) model combining a DHS model in phasor form and a traditional electrical power system model. The IHED model is a convex programming problem and can be easily solved. The effectiveness of the proposed GPM and dispatch model is verified in three test systems. Compared with traditional methods for modeling the thermal inertia, the proposed GPM is more accurate.

Abstract-Using thermal inertia in district heating systems (DHSs) to improve the dispatch flexibility and economy of integrated heat and electricity systems (IHESs) is a research hotspot and difficulty. In most existing studies, the partial differential equations (PDEs) of thermal inertia are approximated by discrete-time models, making it difficult to accurately describe the continuous dynamic processes. In this paper, we propose a novel generalized phasor method (GPM) for thermal inertia in DHSs with constant mass flow. Based on the analytical solution of the PDEs and the Fourier transform, the intractable PDEs are transformed into a series of complex algebraic equations represented by phasors. The GPM has higher accuracy compared to traditional discrete models because it is essentially a continuous model in the time domain. Then, we present a different representation of an integrated heat and electricity dispatch (IHED) model combining a DHS model in phasor form and a traditional electrical power system model. The IHED model is a convex programming problem and can be easily solved. The effectiveness of the proposed GPM and dispatch model is verified in three test systems. Compared with traditional methods for modeling the thermal inertia, the proposed GPM is more accurate.
Index Terms-Energy management, generalized phasor method, integrated heat and electricity dispatch, thermal inertia. Index of dispatch periods.

Parameters and Constants
Cross-sectional area of pipeline b. B n Shunt susceptance at EPS bus n. C w Specific heat of water.
Upper and lower bounds Upper and lower bounds of sin Δθ b .
Active power flow at the start/end of EPS line b.
Reactive power flow at the start/end of EPS line b.
Temperature at the start/end of pipeline b.

I. INTRODUCTION
T ODAY, the relationships between different energy systems, e.g., electrical power systems (EPSs) and district heating systems (DHSs), are getting closer. For example, in north China, most of the heat loads are supplied by combined heat and power (CHP) units. In microgrids, industrial, commercial and residential customers may need various forms of energy at the same time, such as electricity and heat. Under this trend of increasingly close coupling of various types of energy, multienergy systems have become a research hotspot [1]- [3].
Integrated heat and electricity systems (IHESs) are some of the most common multi-energy systems. Many studies have verified that dispatching an EPS and DHS together can improve the economic efficiency and promote renewable energy consumption [4]- [9]. Reference [4] presented a steady-state model of a DHS and studied integrated heat and electricity dispatch (IHED). In [5], IHED with heat storage was studied to consume more renewable energy. Reference [6] studied IHED with a heat pump to increase the schedulability of the isolated community microgrid. The authors of [7] studied IHED and demonstrated that the electric boilers in DHSs can help EPSs accommodate more wind energy. References [8], [9] used thermal energy storage to make an IHES more flexible and reduce operating costs. However, these studies were based on steady-state DHS models and rarely considered the thermal inertia characteristics. Because the time delay between sources and loads of a DHS is much larger than that of an EPS, a steady-state DHS model will result in an unreliable result.
Thermal inertia was considered by many researchers to make IHED more accurate and further explore the dispatch flexibility of IHESs. An accurate thermal inertia model for a district heating network (DHN) was described by partial differential equations (PDEs). In mathematics, the finite difference method [10] is a common method for approximating and solving partial differential equations. It replaces derivatives and partial derivatives with temporal and spatial differences. In [11], [12], a set of difference schemes was proposed to solve PDEs for a thermal dynamic process along a pipeline. However, the finite difference method has rarely been used for IHED because its computational burden is too large.
To overcome the weakness of finite difference methods, researchers used discrete models in the time domain to approximate PDEs without spatial differences. The node method (NM) [13], [14] is the most widely used method for considering the thermal inertia of pipelines. The key idea of the NM is to calculate the lossless temperature at the end of a pipeline by averaging the temperature of the water outflowing from the pipeline in each dispatch period. In [15], [16], Li et al. proposed iterative solution strategies with the NM to solve the IHED. Reference [17] proposed a modified NM that can better consider the variable mass flow and studied the IHED. References [18]- [20] also studied IHED and adopted the NM to take into account the thermal inertia. In [21], [22], a simplified NM was used to study IHED in which the time delay of a pipeline was approximate to an integer multiple of the dispatch interval. Reference [23] studied the transfer function of a DHN in the Laplace domain and used the inverse Laplace transform to establish a discrete time-domain relationship between the heat source and heat load. This relationship is equivalent to the NM and will produce the same error as the NM.
However, the actual thermal inertia process is continuous in time. Therefore, using a discrete-time model to describe the thermal inertia inevitably leads to errors. When the time delay of a pipeline is not an integer multiple of the dispatch interval, the process of averaging in the abovementioned NM and modified NM will lead to errors.
Quality regulation, which means the control method of fixing the mass flow rates and varying the supply temperatures to satisfy the heat demand, is widely used in China, Russia and some Nordic countries [26]- [27]. Under this background, quality regulation was adopted by many researchers [18]- [23]. In this paper, the mass flow in the DHN is also considered constant.
As with the temperature in DHSs, the instantaneous voltage and current in AC circuits are also related to time and space. The phasor method for AC circuits was proposed by Steinmetz C P in 1893 [24], [25]. It transforms the differential and integral equations into complex algebraic equations, thereby simplifying the analysis and calculation of the circuit. Inspired by the traditional phasor method in the field of EPSs, we propose a novel generalized phasor method (GPM) for DHSs with constant mass flow by transforming the PDEs into complex algebraic equations. Based on the new modelling method, a new IHED model considering the thermal inertia of the DHS is developed.
The main contributions of this paper are summarized below. 1) A novel GPM for DHNs with constant mass flow is proposed. The intractable time domain PDE is converted into a series of linear complex algebraic equations denoted by the phasors in the frequency domain. This GPM not only can avoid spatial differences but also has a higher accuracy compared to the NM method. 2) Based on the Fourier transform and the GPM, a novel IHED model is proposed by combining a DHS model in the frequency domain and an EPS model in the time domain. All constraints in the IHED model are convex, and the numbers of variables and constraints in the GPM have similar orders of magnitude to those in the traditional NM, so the IHED model can be solved efficiently. The remainder of this paper is organized as follows. In Section II, based on the analytical solutions of the PDEs, we derive the GPM for a DHN considering the thermal inertia. In Section III, we propose an IHED model combined with a frequency-domain DHS model and a time-domain EPS model. In Section IV, simulation experiments are conducted on three test systems to demonstrate the accuracy and effectiveness of the GPM and the IHED model. Conclusions are drawn in Section V.

II. GPM FOR A DHN
In this section, we propose the GPM for a DHN with constant mass flow. First, we define the phasor variables used in the GPM. Then, the PDEs and the algebraic equations in the time domain are transformed into a series of complex algebraic equations represented by phasors. Finally, a simulation procedure based on the GPM is proposed.

A. Variables in the GPM
In the first step of the derivation of the GPM, we rewrite some variables and parameters of the DHN into Fourier series. They are the thermal energy generation (1), heat load demand (2), temperature at the node (3), temperature at the start of the pipeline (4), and temperature at the end of the pipeline (5).
All Fourier series are represented by phasor variables to facilitate the derivation. In (1)-(5), the function Re(·) denotes taking the real part of a complex number. Taking (1) as an example, q HS i (τ ) is approximated by the superposition of trigonometric functions.Q HS i,kΩ is the phasor variable associated with kΩ frequency, and Ω is the fundamental frequency whose value is discussed in Section III.
A very important point is that all phasors should be defined in rectangular coordinates. In this case, the proposed GPM is linear with the real and imaginary parts of the phasors. If the phasors are defined in polar coordinates, the nonlinear constraints are difficult to apply in the dispatch model.
Next, we introduce the DHN constraints one by one in phasor form. All complex algebraic equations mentioned below are equivalent to a real-part equation and an imaginary-part equation that are both linear.

B. Derivation of the GPM 1) Thermal Inertia Equation of the Pipeline:
A widely used PDE (6) is applied to model the thermal inertia of hot water flows [28]. In (6), τ denotes time, and l denotes length. The second term denotes heat convection, and the third term denotes heat transfer between the environment and water. The conductive heat transfer within the water is not considered in (6) because it is relatively weak compared with the other heat transfer terms. m b is considered a constant.
By applying the superposition theorem, the PDE (6) is decomposed, and the PDE for each phasor is obtained, as shown in (7) and (8).
Using the method of characteristics [29], we can obtain the analytical solution of (7) and (8), as shown in (9) and (10), respectively. Here, By rewriting equations (9) and (10) into phasor form, we obtain the phasor form of the relationship between the temperatures at the start and end of the pipeline as follows: Equation (11) describes the thermal inertia of the pipelines in phasor form. The intractable PDEs (6) are transformed into complex linear algebraic equations (11). When k ≥ 1, equation (11) consists of a real-part equation (12) and an imaginary-part equation (13) that are both linear.

2) Temperature Mixing Equation:
Equation (14) is the temperature mixing equation in the time domain and denotes the energy balance at each node.

3) Temperature Continuity Equations:
The temperature at the start of the pipeline is equal to the temperature of the start node of the pipeline, as shown in (16).
By substituting equations (3) and (4) into (16), we obtain the phasor form of (16) as follows: Remark 1: If the mass flow rate m b is not a constant, the superposition theorem is still valid, and solutions (9) and (10) can be obtained. However, the coefficients A b,kΩ and B b are related to m b (τ ) by integral equations. Therefore, the current GPM is difficult to use in the IHED with a variable mass flow.

C. Application in Simulation
So far, we have established a linear DHN model in phasor form, as shown in (11), (15) and (17). The GPM can be used to simulate the thermal inertia of the DHN. The greatest advantage of the proposed phasor method over a finite difference method is that it requires fewer calculations. The following three steps can be used to calculate the time-continuous expression of the DHN variables: Step 1: Transform the input data into phasors using the Fourier transform.
Step 3: Restore the phasor results at different frequencies to the time domain, and superimpose them. The error of the GPM mainly comes from the truncation error and spectrum aliasing of the Fourier series (1)- (5). Other parts of the GPM are strictly equivalent, which convert equations (6), (14) and (16) to equations (11), (15) and (17). In Section IV, we further verify the accuracy of the model and compare it with the traditional NM.
Compared with a steady-state model in the time domain, the GPM requires listing equations for each frequency. However, the phasors are not related to time, so the GPM does not need to list equations at each dispatch point. In the next section, we apply this GPM to formulate the IHED model.

III. IHED MODEL
In this section, to apply the GPM to the IHED, we propose a preprocessing method for DHS data and constraints (21) to determine the initial state of the DHS. Then a IHED model is proposed, consisting of a quadratic objective function, a linear DHS model based on the GPM and a second-order cone (SOC) EPS model based on an existing SOC relaxation method.

A. Preprocessing of DHS Data
In addition to the forecast data for the dispatch day, we also need to consider historical data because of the time delay characteristics of the DHS.
In this paper, we consider historical data for the day before the dispatch day. The time lag of a large-scale DHS is approximately a few hours, so considering one day's historical data is sufficient.
Assume that there are K dispatch periods per day and that the dispatch interval Δτ satisfies KΔτ = 24 hours. The dispatch point is denoted by τ w , and the dispatch points of the dispatch day are τ w = wΔτ , w = 1, 2, …, K. The dispatch points of the previous day are τ w = wΔτ, w = −K+1, −K + 2, . . . , 0.
There are 2K scheduling points in two days. According to the sampling theorem, we can fit a K-order Fourier series. This matches the symbol K used in Section II.
3) Initial Value Constraints of the Node Temperatures: The initial boundary values are necessary for solving the PDEs, as well as solving the GPM. In this paper, the historical node temperatures are considered to formulate the initial boundary value constraints of the PDEs (6). We only need to consider the initial values of the temperatures at the heat source nodes (21), because the state of the mass-constant DHS is uniquely determined by the heat load power and the temperatures at the heat source nodes.
Remark 2: According to the properties of RDFT, the imaginary parts ofṪ N n,0 andṪ N n,KΩ are equal to 0. That is, the imaginary parts ofṪ N n,0 andṪ N n,KΩ are not considered as decision variables in the IHED.

4) Boundary Constraints of the DHN:
The node temperatures cannot exceed the limit. The bound constraints are expressed in phasor form as follows: The operating constraints of gas boilers are as follows: Constraints (24) and (25) are the output range constraints and the ramping constraints, respectively.
Remark 3: Constraint (24) implies that the gas boilers are always on, but in some cases, the on-off states of gas boilers can be adjusted, and the IHED should consider mixed-integer linear constraints such as . The GPM is linear, so it can be easily combined with the on-off model of gas boilers. After considering the on-off states of gas boilers, the DHS model becomes a mixed-integer linear model. If other parts of the IHED are linear or mixed-integer linear, the IHED is a mixed-integer linear programming (MILP) problem that can be solved efficiently by commercial solvers such as Gurobi [31]. In addition, if a SOC EPS model and a quadratic objective function are considered, the IHED becomes a mixedinteger SOC programming (MISOCP) problem. Although the MISOCP problem is more complex, it can still be solved by Gurobi.
The operating constraints of CHP units are as follows:   Fig. 1 shows the feasible operating region of CHP unit i, which is a convex polygon. Constraints (26) and (27) use the linear combination of vertexes (H κ,i , P κ,i ) to express the electrical power and heat power of CHP units. Constraints (28) and (29) bound the combination factor ζ τ w κ,i to ensure that p ES,τ w i and h HS,τ w i are within the feasible operating region. Constraints (30) and (31) are the ramping constraints. For a more detailed model of CHP units, please refer to [33]. In [33], a concave piecewise linear model is presented to consider the partial-load performance of CHP units. The concave piecewise linear model can be converted into a mix-integer linear model, and then it can be used with the proposed GPM.

C. EPS Constraints
The EPS constraints are expressed in (32)-(56). For the computation complexity challenge, a SOC relaxation of the AC power flow is adopted to model the EPS. This SOC relaxation (32)-(51) was proposed in [34], and it is valid for both mesh and radial power networks.
∀n ∈ N E , w = 1, . . . , K. (33) In the above formulations, the set I ES is the union of the sets I T U , I CHP and I W ind . Constraints (32)-(33) represent the active and reactive power balance. Constraints (34)-(37) define the active power loss and reactive power loss of the lines. The drop of the squared voltage magnitude along line b is restricted by (38). Constraints (39)-(49) are the convex approximation of the AC power flow constraint v τ w

D. Objective Function
The objective is to minimize the total operating cost of the IHES, and this cost includes the costs of CHP units, non-CHP thermal units and gas boilers. The total operating cost is expressed as Generally, equation (58) is convex [34], [35]. Thus, the objective function (57) is a convex quadratic function.

E. Matrix Form of the IHED Model
The proposed IHED model includes frequency-domain and time-domain variables. In the following matrix form, the vector of all phasor variables (frequency-domain variables) is denoted by x fd , and the vector of all time-domain variables is denoted by x td . The matrix form of the IHED model is expressed as In the above model, (62) and (63) are the matrix forms of time-domain constraints (26)-(56). (64) is the matrix form of frequency-domain constraints (11), (15), (17), (21) and (22). (65) is the matrix form of inverse Fourier transform (20).
The model (61)-(65) is a second-order cone programming (SOCP) problem and can be easily solved.

F. Discussion of the Number of Variables and Constraints
Compared with solving the optimization problem, the preprocessing (18) requires very little calculation time and can be ignored. Moreover, the EPS model described in Section III.C is a traditional discrete time-domain model. Therefore, we only need to discuss the number of variables and constraints in the DHS model. We assume that there are N nodes, N S heat source nodes, M branches and S heat sources in the DHS. The results of a comparison of the number of variables and constraints in the GPM and a steady-state model [4] are shown in Table I.
For a fair comparison, we consider a complex variable as two real variables and a complex constraint as two real constraints. In terms of the variables and constraints of the DHN, the GPM has approximately twice the number of variables and constraints as the steady-state model. These additional constraints and variables are mainly due to the fact that two days of data are considered in the GPM, while only one day of data is considered in the steady-state model. Based on the above discussion, we can conclude that the numbers of variables and constraints in the GPM are not high and have similar orders of magnitude to those in the steady-state model.

IV. CASE STUDIES
In this section, the accuracy of the proposed GPM is tested and compared with the analytical solution and the traditional NM. Experiments for the IHED model (61)-(65) are conducted on a small-scale IHES and a large-scale IHES. Moreover, the IHED model based on the GPM is compared with the IHED model based on the NM. All tests are performed on a laptop with a six-core processor running at 2.20 GHz with 16 GB of memory. Programs are coded using MATLAB R2018a. The IHED model is solved by Gurobi 9.0.

A. Simulation of a Single Pipeline
The simulation is conducted on a single heating pipeline to show the accuracy of the proposed GPM. The parameters of the pipeline are in [36]. The input data shown in Fig. 2 are available from a real DHS in Jilin, China. The red line is the data of the temperature at the start of the pipeline for two days.
Based on the sampling points of the input curve, we use the GPM and the traditional NM to calculate the temperature at the end of the pipeline between 1:00 and 24:00. Assuming time delays of 1 hours and 1.25 hours and sampling intervals of 1.5 hours, 1 hour and 0.5 hour, the simulation results are shown in Fig. 3. Δτ denotes the sampling interval, and τ delay denotes the time delay of water flowing from the start to the end of the pipeline. The simulation process of the GPM is described in Section II.C. According to the sampling theorem, the order of the phasors is half of the total number of sampling points.  The NM can obtain only a few discrete points, while the GPM can obtain a continuous curve. To be fair, we compare only the solutions at each discrete point. Table II shows the numerical results of the GPM and the NM. When the time delay is 1 hour and the sampling interval is 1 hour or 0.5 hour, the solution of the NM is equal to the analytical solution and the solution of the GPM. When the time delay is 1.25 hours, the smaller the sampling interval, the closer is the solution of GPM and NM to the analytical solution. Furthermore, with the same sampling interval, the solution of the GPM is closer to the analytical solution compared to the solution of the NM.
As shown in Fig. 3 and Table II, the proposed GPM is more accurate than the traditional NM when the time delay is not an integer multiple of the sampling interval. The reasons behind this finding are explained below. The NM is a discrete model in the time domain, and the temperature at the end of the pipeline is calculated by a weighted summation of the temperature at the start of the pipeline (66). When the time delay of the pipeline is not an integer multiple of the time interval, this process of weighted summation will cause errors. The phasor method that we proposed is essentially continuous in the time domain and does not include the process of weighted summation, so this method is more accurate than the NM. The accuracy of the phasor method mainly depends on the fitting accuracy of the Fourier series. Generally, the temperature in the DHS changes continuously and slowly, so the Fourier series often has good accuracy.
B. Small-Scale IHES Fig. 4 shows the configuration of the test IHES. This smallscale IHES consists of a six-bus EPS and a six-node DHS. The time delay of each pipeline is marked in Fig. 4. The time needed for hot water to flow from the heat source to each load is 2 hours. The electrical load demand, heat load demand and available wind power are shown in Fig. 5. The ambient temperature is set to 10°C. For more detailed parameters, please refer to [37]. For convenience, we assume that both the predicted data and dispatch plans vary in the diurnal period. The dispatch interval is 1 hour. Experiments are conducted for the following two cases.    Fig. 6, along with the total heat load curve and the wind power dispatch curve. Fig. 6(a) shows that when the GPM is used, the heat source is decoupled from the heat load. The optimal cost of IHED based on the GPM is $247,981.90. If IHED is modeled using the steadystate model, the optimal objective value is $260,655.94. The proposed IHED model can obtain more economic benefits and utilize 212.49 MWh more wind energy than the steady-state model, as shown in Fig. 6(b).
The dispatch plan of the steady-state IHED is not only costly but may also cause problems such as insufficient heating temperatures. Based on the initial state and the dispatch plan of CHP units solved by the steady-state IHED, we use analytical methods to calculate the temperature at node 2, as shown in Fig. 7. The red curve solved by the steady-state model is the optimization result of the steady-state IHED, and the blue curve is the simulation  result based on the dynamic model (6). The red curve satisfies the minimum temperature constraint. However, the blue curve is lower than the minimum temperature of 65°C for more than half a day, which means that the steady-state model will result in insufficient heating temperatures for all loads.
After solving the IHED model based on the GPM, the operator can also obtain the continuous temperature distribution of the whole DHN in time and space. The pipeline between node 1 and node 2 is adopted as an example to visualize the dynamic characteristics of DHNs in continuous time and space. Fig. 8 illustrates the temperatures along the pipeline in continuous time and space. Both the temperature decay along the pipeline and the temperature change over time can be seen in Fig. 8.

2) Case 2 (With Fixed Return Temperatures):
After fixing the temperature of the return water at the heat source node, we solve the IHED model again. To better demonstrate the advantages of the GPM, the heat loss is ignored in this case. The results are shown in Fig. 9. The red dashed curve is obtained by shifting the total heat load curve to the left by 2 hours. The heat output curve and the heat load curve have a strict translation relationship, and the time difference is 2 hours, which is equal to the time needed for the hot water to flow from the heat source to the heat load. This finding is consistent with intuitive physical insights when the return temperature and the mass flow rate are fixed.
In this case, we also use the NM to model the thermal inertia of pipelines and compare the results with those of the GPM. The solution of the IHED model based on the NM is shown in Fig. 10. The heat output curve is not a translation of the heat   TABLE III  CALCULATION TIME AND ERROR BETWEEN THE TOTAL HEAT OUTPUT AND SHIFTED TOTAL HEAT LOAD load curve because the time delay of some pipelines is not an integer multiple of the dispatch interval. If the scale of the DHS is larger, the error of the NM will increase. As shown in Table III, this case verifies that the IHED model based on the GPM has a higher accuracy than the IHED model based on the traditional NM.

C. Large-Scale IHES
Based on a real system in Jilin Province of Northeast China, we design a large-scale IHES that consists of a 319-bus EPS and a 55-node DHS. The parameters of this test system can be found in [38]. IHED is modeled based on the GPM. In this case, the SOC relaxation of the AC power flow requires a lot of calculation time (more than 500 s), so we use the DC power flow (67)-(68) to model the EPS and verify the computational efficiency of the GPM. The results of the IHED model are shown in Fig. 11. In this large-scale test system, the heat source is also decoupled from the heat load, and more dispatch flexibility is utilized.
∀b ∈ B E , w = 1, . . . , K. (68) Moreover, we generate more scenarios to test the calculation time of the proposed IHED by changing the values of the predicted heat load demand and the available wind power. The scaling factors u and o are defined to scale the available wind power and the predicted heat load demand, respectively. The calculation time in each scenario is shown in Table IV. It can be concluded from Table IV that the computational efficiency of the GPM is comparable to that of the NM. When the EPS is modelled by the DC power flow model, the IHED model based on our GPM requires very little calculation time. This is because the proposed GPM and the DC power flow model are both linear and concise, with a small number of constraints and decision variables.

V. CONCLUSION
This paper presents a GPM for DHNs with a constant mass flow to consider thermal inertia. All constraints on DHNs are formulated using phasor variables. The main feature of the proposed GPM is that it can describe the continuous changes of temperature in time and space, while traditional methods are discrete in the time domain. The GPM is derived based on the analytical solutions of PDEs, and its accuracy is satisfactory. Moreover, the GPM is linear and has only a small number of constraints and decision variables, so this method can be easily applied in IHED. By applying the GPM and the Fourier transform, we propose a novel IHED model combining a frequency-domain DHS model and a time-domain EPS model. This model can accurately consider the dynamic processes of DHSs. In addition to traditional dispatch plans of units, the IHED model can be solved to obtain the continuous temperature distribution of a whole DHN in time and space to help operators make decisions. Three test systems are used to demonstrate that the GPM is more accurate compared to traditional methods and has a similar computational efficiency to those of traditional methods. Moreover, the IHED model can further utilize the flexibility of DHNs to improve the economics and reduce wind curtailment.