Fast Decomposed Method for Dynamic Energy Flow Calculation in Integrated Electricity and Heat System

The development of integrated electricity and heat system (IEHS) promotes the renewable energy penetration and the complementary coupling between the electric power system (EPS) and the district heat system (DHS). The efficient energy flow calculation (EFC) in the IEHS has become a hotspot as the system scale and complexity increase. This paper firstly constructed the dynamic DHS model using Euler difference scheme. On this basis, the fast decomposed method is proposed to decompose the system with various scales and topologies, which improves the computation efficiency. Moreover, the modification mode is then developed to modify the EFC results in the decomposed systems and ensure the accuracy. Numerical simulation verifies the accuracy of the proposed method and its superiority in computation efficiency.


INDEX TERMS
The growing energy crisis and environmental problems have promoted a new requirement to build and efficient and clean energy structures [1][2]. As a typical form of energy evolution, the integrated electricity and heat system (IEHS) combines the electric power system (EPS) and the district heating system (DHS) through energy cascade utilization and complementary optimization, which promotes the renewable energy absorption and operational flexibility [3][4]. Due to the superiority in utilization efficiency and emission reduction, the IEHS has been promoted broadly. Nevertheless, the electricity and thermal transmission follow different energy laws and mathematic representations, which hinders the comprehensive and efficient analysis of IEHS. The energy flow calculation (EFC) provides the data for operational optimization, state VOLUME XX, 2017 1 estimation, and security control, which is importance for the harmonious operation of IEHS [5][6]. Massive studies have been performed for EFC in IEHS using different models and methods, as summarized in Table I. Reference [7] constructed the energy hub model for the optimal operation, where the uncertainty of parameters and thermal storage were both analyzed. On this basis, [8] took the network constraints into consideration when solving the unit commitment problem, where the operation flexibility of IEHS can be improved. Reference [9] further modeled the hydrogen storage to formulate the stochastic scheduling problem of the integrated energy system, and the operational risks were quantitatively evaluated. The mentioned literatures built the energy flow model for different applications, but the nonlinearity in the independent energy system was neglected. Therefore, [10] developed a unified model to solve the model nonlinearity based on Newton method, and [11] further applied the unified model into an integrated electricity, gas, and heat system. The unified model overcomes the nonlinearity using iteration, but its computation burden enlarges sharply when the system scale increase [12]. Moreover, the Newton-based unified model relies on the suitable initial values to converge, which is difficult to be obtained in engineering [10][11][12]. Although a decoupling method was proposed in [13][14] to solve the convergence problem, the computation complexity still exists especially for the large-scale system. Besides, the aforementioned literature all adopted the static model to describe the energy flow distribution in IEHS, while the dynamics in the IEHS were neglected.
To solve the mentioned gaps, [15] developed the node method to describe the thermal inertia, and [16] further classified the interaction between the EPS and the DHS into four stages using the node method. Despite that the node method employs the finite element method to model the dynamics in IEHS, the time-varying thermal states and the bidirectional hydraulic states in DHS were not involved. Thus, [17] presented the partial differential equations (PDEs) to construct the dynamic DHS model, where the dynamic pipe temperature transmission process was explicitly studied. Reference [18] further optimized the differential step to ensure higher accuracy and efficiency, and the model was then used for real-time optimization. On this basis, the uncertainty of renewable energy resources was included in [19], where an iterative method was proposed based on the dynamic IEHS model. Literatures above study the EFC for IEHS from different aspects, however, several problems still exist: 1) The computation amounts greatly increase along with the increment of system scale, which is not suitable for simultaneous analysis; 2) The analysis for the meshed network is still blocked, where the bidirectional hydraulic states in IEHS are complicated; 3) The improvement of the efficiency in EFC mainly based on the static model, while the researches using dynamic model are few.
In this paper, the dynamic IEHS model is established based on the finite difference method. After that, a fast decomposed method is proposed to decompose the system into several subsystems. Then, the modification mode between the decomposed systems is comprehensively investigated, which relates the state variables at different time-step. To validate the performance of the proposed methods, numerical studies are further conducted. The contributions of this paper can be summarized as below.
1) This paper proposed the dynamic model for IEHS using finite difference method with Euler difference scheme.
2) A fast decomposed method is developed to reduce the studied system scale and computational burden.
3) The modification mode between the EFC results in the decomposed systems is established to ensure the accuracy The remainder of this paper is organized as follows. Section II describes the modeling of dynamic IEHS. Section III introduces the fast decomposed method and the modification mode. Section IV and Section V are the case study and the conclusion, respectively.

A. ELECTRIC POWER SYSTEM
Due to different operational timescales, the EPS can rapidly become static in milliseconds, while the thermal states in the DHS are still in the dynamic process. Thus, the EPS is commonly considered static in IEHS analysis, and the alternating current model is adopted to describe the voltage and active/reactive power distribution, which can be expressed as [13]: ( cos sin )

B. DISTRICT HEATING SYSTEM
The DHS realizes the thermal power exchange using the water flow through the supply/return network with the same structure. Thus, the DHS model usually involves the hydraulic and thermal parts. This paper focuses the DHS operating under the quality regulation mode, where the hydraulic states in the DHS are fixed. Due to the lower control difficulty and better operation stability, the quality regulation is widely applied in the practical engineering, particularly in North Europe, North China, and Russia [11][12] [19]. The structure of DHS is shown in Fig. 1

1) HYDRAULIC MODEL
The hydraulic model describes the distribution of mass flow and node pressure analogy to the Kirchhoff's current and voltage law [13]. Firstly, the mass flow that comes into a node equals to the sum of mass flow that leaves and consumes at the node, which is expressed as: where A is the reduced order node-branch incidence matrix, aij=1/-1 if pipe j starts/end at node i, else aij=0. Secondly, the pressure drop around a closed loop equals to zero, which is expressed as: where B is the loop-branch incidence matrix, bij=1/-1 if the direction of pipe j is consistent/inconsistent with loop i, else bij=0. The second equation in (4) means that the pressure drop along the pipe mainly depends on its mass flow and inherent properties [13].

2) DYNAMIC THERMAL MODEL
The dynamic process of temperature transfer along the pipe is usually formulated as the partial difference equation (PDE) with the static heat conduction in a fluid ignored, which is expressed as [18,19]: In this paper, the finite difference method using the Euler difference scheme is proposed to solve (5), and the spatialtemporal sphere to be studied is discretized into several difference points, which is expressed as: where L is the pipe length, Γ is the temporal interval. On this basis, the partial differential item in (5) can be rewritten as (6) using the forward difference quotient, where O represents the High order remainder term.
Differentiate (5) at point (xi, tk), and then substitute (7) into (5), we can transform the partial difference equation into a linear function as (8): where λ1, λ2, and λ3 are constant coefficients for simplification. Equation (9) indicates that the current temperature not only depends on the current states, but the previous states are also should be considered. Thus, the results of (9) is often calculated with a series of initial and boundary conditions. Moreover, the temperature mixture equations are employed to calculate the node temperature, which is expressed as [13]: (10) where E in n is the set of pipes that ends at node n, V is the set of nodes in DHS, E out n is the set of pipes that starts at node n, T out i is the temperature at the end of pipe i. Besides, the nodes in DHS are modeled as the heat exchangers, which can be expressed as [19]: ( ) Finally, the node temperature after mixing equals to the temperature at the start of certain pipes, which is expressed as:

C. COUPLING UNIT
According to the coupling mechanism, the EPS and the DHS are mainly connected with two types of equipment: the cogeneration equipment like combined heat and power (CHP) units and the energy conversion equipment like electric boiler (EB) and heat pump (HP). In this paper, the electric boiler and heat pump are considered, which both consume electric power to generate thermal power. Thus, the EB and HP are treated as the source in the DHS and load in the EPS, and the influence of DHS on the EPS is reflected through the coupling units. The model of EB and HP are given as (13) with the thermalelectric ratio.
We assume that the IEHS operate at the heat-load-following mode [10][11][12][13], in which condition the DHS is an equivalent load of EPS, and the electric power is determined according to the states in DHS. The coupling relationships between the DHS and the EPS are as shown in  To analyze the EFC in IEHS more clearly, the nodes in DHS are classified according to its known states, which is similar with that in the EPS. As shown in Table 1, the nodes in DHS are sorted into four categories. The DHS has two source types, the thermoregulation node is worked for temperature regulation, whose function is the same as the PV bus in the EPS; the slack node is worked for thermal power regulation, whose function is the same as the slack bus in the EPS.

III. FAST DECOMPOSED METHOD FOR DYNAMIC ENERGY FLOW CALCULATION
The EFC for the IEHS includes 3 parts: the EFC for the This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2021.3116810, IEEE Access VOLUME XX, 2017 1 DHS, the coupling units, and the EPS, which are performed consecutively at the heat-load-following mode [10][11][12][13]. In addition, the researches about the EFC in the EPS containing the distributed methods, the decentralized methods, and the parallel methods are very mature. Therefore, the improvement in this paper mainly focuses on the EFC for the DHS [12][13], [19].

A. SYSTEM DECOMPOSITION
The increment of variables due to the large and complex system scale challenges the efficiency of EFC for the DHS. In this part, a fast decomposed method (FDM) is proposed to reduce the computation burden using topology decomposition. As shown in Fig. 3, the original DHS is decomposed into several subsystems at the intermediate nodes. Among the subsystems, one is called the master system (with the original slack node), and the others are called the slave systems. The decomposed nodes are seen as the equivalent loads for the master system, which transfer the energy flow to their connected slave systems. In the slave systems, the decomposed nodes are seen as the equivalent sources for the slave system, which obtain the energy flow from the master systems to provide the thermal power for the inner loads. Fig.  3 gives the examples of the decomposition of radial and meshed DHSs for better understanding.

1) RADIAL SYSTEM
As shown in Figure 3(a), the radial DHS is decomposed at node 3, which is an intermediate node in the original DHS. For the master DHS, the decomposed node is an equivalent load node (node 3); for the slave DHS, the decomposed node is an equivalent slack node (node 5). The overall energy flow distribution is the combination of the results in the master DHS and the slave DHS.

2) MESHED SYSTEM
As shown in Figure 3(b), at least two intermediate nodes are required to decompose the meshed DHS. In the master DHS, the load nodes (node 9 and node 7) are virtual loads to reflect the load consumption. In the slave DHS, the two virtual source nodes (node 2 and node 3) are the slack node and the thermoregulation node, respectively. The energy flow distribution in the original meshed DHS is also the combination of the results in the decomposed systems.

B. MODIFICATION MODE
With the system scale reduced, the state variables at the decomposed node should be consistent in each decomposed system. According to topological properties, the modification mode for radial and meshed DHSs are presented below, respectively. Since the mass flow and node pressure are fixed under the quality regulation model, only the node temperature and the thermal power need to be modified.
The decomposed DHSs in the radial DHS are only connect through one intermediate node, which is seen as the load node in the master system and the slack node in the slave system. Take the system in figure 3(a) as an example, the procedure of the decomposition for the radial DHS is as follows. Besides, during the decomposition process, the initial conditions of the temperature are maintained (a) Give the assumed supply temperature (boundary condition) at time k of node 5 (T k1,j n5,s ) in the slave DHS at j th iteration.
(b) Perform the energy flow calculation in the slave DHS. Obtained the result of the thermal power of node 5 at time k at j th iteration (ϕ k1,j 5 ). (c) Treat the provided thermal power of node 5 in the slave DHS as the consumed thermal power of node 3 at jth iteration (ϕ k2,j 3 =−ϕ k1,j 5 ). (d) Perform the energy flow calculation in the master DHS. Obtain the result of supply temperature of node 3 (T k2,j n3,s ) in the master DHS (e) Modify the supply temperature of node 5 in the slave DHS at time k at (j+1) th iteration using (14), where δT is the temperature convergence limit and equal to 10-4 in this paper 2 (f) Determine whether the convergence condition in (15) satisfied. If so, the procedure is over, else return to step (a).
The decomposed DHSs in the meshed DHS are connect through at least two intermediate nodes, among which one is chosen as the slack node and the others are chosen as the thermoregulation nodes. Take the system in figure 3(b) as an example, where node 2 is chosen as the slack node and node 3 is chosen as the thermoregulation node. The procedure of the decomposition for the meshed DHS is as follows. Besides, the treatment of the initial conditions is the same as those in the radial DHS.
(a) Give the assumed supply temperature of node 2 (T k1,j n2,s ) and the provided thermal power at time k of node 3 (ϕ k1,j 3 ) at time k in the slave DHS at j th iteration.
(b) Perform the energy flow calculation in the slave DHSs. Obtained the results of the provided thermal power of node 2 (ϕ k1,j 2 ) and the supply temperature of node 3 (T k1,j n3,s ) at time k at j th iteration.
(c) Treat the provided thermal power of node 2 and node 3 in the slave DHS as the consumed thermal power of node 7 and 9 in the master DHS at j th iteration (ϕ k2,j 7 =−ϕ k1,j 2 , ϕ k2,j 9 =−ϕ k1,j 3 ). Obtain the result of supply temperature of node 7 (T k2,j n7,s ) and node 9 (T k2,j n9,s ) in the master DHS. (e) Modify the supply temperature of node 2 and the provided thermal power of node 3 in the slave DHS at time k at (j+1) th iteration according to (16 (16) where T k2,j n9 ,r is the return temperature of node 9 at time k at (j+1) th in the master DHS.
(f) Determine whether the convergence condition in (17) satisfied. If so, the procedure is over. Else, return to step (a).

Ⅳ. CASE STUDY
To validate the superiority of the proposed method, numerical simulations are performed in the radial and meshed DHS, respectively. The program is coded by Matlab R2018b on a PC with Intel i7 4710 CPU and 4GB RAM, and the proposed FDM is compared with the static model in [8] and the tradition dynamic EFC method (TM) in [19]. The structure of the studied IEHS is shown in figure 4, which consists of the well-know the DHS in Barry Island (node 31 is selected as the slack node) [13] and the IEEE 33-bus EPS [19]. The DHS and the EPS are connected through one heat pump and two electric boilers. The meshed DHS is firstly decomposed into two radial DHSs at node 5 and node 25, which aims to validate the effectiveness of the decomposition for meshed DHS. On this basis, the decomposed DHSs are further decomposed into the smaller radial DHSs, which aims to validate the effectiveness of the decomposition for radial DHS. Finally, the dynamic EFC for the IEHS is performed, and the results is compared with those of TM. The process of two-step decomposition for the DHS is given in figure 5. The ambient temperature is 0℃ and the whole analysis period is 3h. The temporal step and the spatial step are 1min and 50m, respectively. It should be noted that the results by the method in [19] is only used for comparison due to the lack of measured data, and the comparison with real measurements will be performed in our future work to further validate the correctness of the proposed method.  As shown in Fig. 4, four small radial DHSs exist after the two-step decomposition with one master DHS and three slave DHSs. For the master DHS, the slack node is the same as that in the original DHS; for the slave DHSs, the decomposed nodes are chosen as the slack nodes, respectively. Thus, the dynamic EFC in this condition contains two steps: 1) Firstly, the EFC is performed until the modification between the master DHS and the slave DHS1, the slave DHS2 and the slave DHS3 finished. On this basis, the EFC is then performed until the modification between the master DHS(b) and the slave DHS(a) finished.
The results of dynamic EFC are shown in Fig. 6, where the results by the static model in [8]. According to the Fig 6. (a) and (b), the results by the proposed FDM and TM in [19] are almost the same. The mean relative error (MRE) of the FDM is 0.263%, indicating the accuracy of the proposed method. The reason for the accuracy of FDM mainly contains two aspects: 1) There is no simplification during the EFC process and the proposed method only reduce the computation scale. 2) The global convergence condition in FDM is the same as that in the traditional method once the modification between the decomposed systems finished. As shown in figure 6(a), the supply temperature of node 2 locates at the lower right at that of node 1, indicating the thermal loss and transmission delay, respectively. The transmission delay is about 7min. The transmission delay in the return network is similar as that in the supply network, as VOLUME XX, 2017 1 shown in figure 6(c). However, node 1 is the equivalent load node in the return network. Therefore, its temperature is lower than that of node 2 in the return network. The comparison between results by the dynamic model and the static model is shown in Fig. 6(c) and (d). Since the temperature distribution in static model only depends on the source output, the obtained temperature is similar with that of the given boundary conditions. In contrast, the initial conditions influence the temperature distribution in dynamic model significantly. Therefore, the calculated results of the proposed FDM and the static model are obviously different, with the MRE of 2.65%. However, the evolution trends of the supply temperature and the supplied thermal power at the slack node are similar, as shown in Fig. 6(d).
The modification processes in the dynamic EFC at 1.5h in the DHS are shown in Fig. 7, which explicitly demonstrates the mismatch of supply temperature at node 5 and node 25 during each iteration.
The average computation time over the whole period by TM is 4.7108s, while that by FDM is 2.7681s. The improvement is over 40% because the computation burden during each iteration is greatly reduced due to the smaller scale of the decomposed DHSs. Moreover, the decomposed DHSs are all radial, thus the judgement of the node location and bidirectional mass flow is no longer needed, which is vital for the EFC in the meshed systems.
The calculated voltage amplitude and the active power loss at time 1.5h is shown in figure 8, which indicates that the results are almost the same in the EPS with the error of 6.16-7%. The differences of the results by the FDM and the TM for the EPS is mainly resulted from the calculated electric output of the coupling units, which is determined by the energy flow distribution in DHS. Since the error in the DHS is relatively small, the electric power flow is slightly influenced.

Ⅴ. CONCLUSION
To provide the accurate and efficient information for simultaneous analysis in the large-scale IEHS, a dynamic model is constructed based on the finite difference method using the forward Euler scheme. Then, a fast decomposed method is proposed to reduce the system scale. In the proposed method, the DHS with various scales and topologies are decomposed into several radial systems at certain nodes, which the temperatures at different time-steps are related. On this basis, the modification mode is developed to ensure the consistentence of the states at the decomposed node, and thus the accuracy can be guaranteed. Case study validates the accuracy of the modification mode in both radial and meshed systems, with the improvement over 40% in computation efficiency.
In our future work, we will extend the proposed method into the DHS under quantity regulation mode. Besides, more attention will be paid into the accurate modeling of equipment considering the dynamic characteristics during the energy flow calculation.