Steady-State Power Flow Analysis of Cold-Thermal-Electric Integrated Energy System Based on Uniﬁed Power Flow Model

: The integrated energy system includes various energy forms, complex operation modes and tight coupling links, which bring challenges to its steady-state modeling and steady-state power ﬂow calculation. In order to study the steady-state characteristics of the integrated energy system, the topological structure of the cold-thermal-electric integrated energy system is given ﬁrstly. Then, the steady-state model of the power subsystem, the thermal subsystem, the cold subsystem and the distributed energy station are established, the uniﬁed power ﬂow model is established, and the Newton Raphson algorithm is used to solve the uniﬁed power ﬂow model. Finally, the inﬂuence of the key technical parameters on the steady-state power ﬂow of the integrated energy system is analyzed. Research results show that the photovoltaic power generation plays a supporting role in the voltage of each bus; with the increase of electric load power, the unit value of bus voltage decreases continuously; the water supply temperature of the source node has a greater impact on the steady-state ﬂow in the pipeline and the water supply temperature of each node; the pipeline length of the heat network has a greater impact on the end temperature of the pipeline, the water supply temperature, and the return water temperature of each node. The analysis results can support the planning, design, and optimal operation of the integrated energy system.


Introduction
Currently, a variety of energy production systems, such as the electric power system, water conservancy system, wind power system, thermal power system, cold energy system, etc., are often operated independently and regulated separately. Each system is lacking in necessary coupling connection and coordinated control, so that the energy utilization rate is low, and energy waste phenomenon is serious. The inefficient energy utilization is hardly in line with the goals of social sustainability and reducing carbon emission [1][2][3]. Therefore, the Energy Internet has been proposed to solve existing problems in energy production, transmission, and utilization [4,5]. The integrated energy system (IES) acts as a basic part of the Energy Internet, and can realize efficient utilization of different energy sources within a region (such as a city, community, industrial park, etc.). IES refers to an integrated system formed by the energy generation, transmission and distribution, conversion, storage, consumption, and other links in the process of planning, construction, and operation [6,7].

Model of Cold-Thermal-Electric IES
The topology of the cold-thermal-electric IES is shown in Figure 1, which was adopted as the research object in this paper. The 33 node three-phase balance electric power systems [15] and the modified 9 node thermal power system [16] were adopted. The cooling power system is directly connected with the load in order to reduce energy transmission losses. storage, consumption, and other links in the process of planning, construction, and operation [6,7]. It greatly promotes energy consumption and improves the energy efficiency; thus, it has important development prospects and utilization scenarios in the future. On the basis of system modeling, the simulation analysis of the IES to reveal its operation mechanism and dynamic characteristics is the premise of the application. For the modeling and simulation analysis of the IES, it is necessary to solve the problems of multi-time scale, high dimension, and nonlinearity. A mathematical model of electric power subsystem and natural gas subsystem in IES is established in literature [8]. It presents a global calculation method of multi energy flow of electricity and gas, which can give an accurate solution by using fewer iterations. In reference [9], considering the time and space correlation between load and renewable energy generation, a multi-energy flow calculation method based on nataf transformation and point estimation was proposed. In reference [10], the mathematical model of production, transmission, and consumption links of the IES was established, and the overall mathematical model of the IES was obtained by connecting each link with energy hub. In reference [11], based on the relationship between electricity and heat output of the cogeneration system, different operation modes of the energy hub were analyzed. On this basis, a calculation method of multi-energy flow is proposed, but only the heat load was considered in the paper, and the heat network model was not considered. In reference [12], the mathematical model of the three subsystems of electricity, heat, and gas in the IES was established, the optimal energy flow model was established by selecting appropriate optimization variables, and the intelligent optimization algorithm was used to solve the model. Among the modeling methods of the IES, the concept of the energy hub, proposed by Professor G. Anderson of the Zurich Federal Institute of Technology, has attracted wide attention from academic circles [13,14]. The energy hub highly abstracts and classifies the energy supply and application, which embodies the concept of multi-energy coordination. Although the steady-state modeling of the IES has been carried out in preliminary research at home and abroad, the existing research of the cold-thermal-electric coupling link contains less energy conversion equipment, thus, there is still a lack of practicality in the project. In addition, the accuracy and solving speed of the model still need to be improved.
This article mainly studies the unified power flow of the cold-thermal-electric IES. In Section 2, the electric power subsystem, thermal subsystem, cold subsystem, and the coupling link are modeled to establish a cold-thermal-electric IES model. In Section 3, the unified energy flow model is constructed and Newton-Raphson is used to solve the model. On this basis, Section 5 studies the influence of the key technical parameters on the steady-state power flow of the IES. Finally, the whole paper is summarized in Section 6.

Model of Cold-Thermal-Electric IES
The topology of the cold-thermal-electric IES is shown in Figure 1, which was adopted as the research object in this paper. The 33 node three-phase balance electric power systems [15] and the modified 9 node thermal power system [16] were adopted. The cooling power system is directly connected with the load in order to reduce energy transmission losses.     Thus, the electric power subsystem, thermal subsystem, cooling power subsystem, and coupling link component were modeled in this chapter.
All energy conversion components are modeled in the coupling link, including a combined heat and power (CHP), heat pump (HP), electric boiler (EB) and other electric-thermal coupling components, absorption chiller (AC) and other electric-cooling coupling components, electric chiller (EC) and other electric-cooling coupling components, and an electric power and cooling joint supply (PC) device and other thermal-electric-cooling coupling components.

Model of Electric Power Subsystem
The traditional model of the power distribution system is adopted for the electric power subsystem. The influence of the three-phase imbalance on power flow calculation is ignored [17]. If the power distribution network has n e nodes, node1 is the balancing node, node2 to node1 + n pv are PV node, and the rest of the nodes belong to the PQ node. The model of the power distribution network is the node power equation reflecting the relationship among node power, node voltage and phase angle:

Model of Thermal Power Subsystem
Energy is transmitted with water as media in the thermal subsystem, thus, they are modeled by the hydraulic model and thermal power force model.

Hydraulic Model
Based on the graph theory, the characteristics of the pipeline in the thermal subsystem are described, and the flow law in the thermal subsystem is modeled by Kirchhoff's law for reference. The continuity equation of energy flow can be described by Equation (5): Head loss is the pressure change per unit length caused by pipeline friction. In a closed circuit, the sum of head losses is zero: where B refers to the loop-pipeline incidence matrix of the thermal subsystem and h f refers to head loss (m), which can be expressed with Equation (7): where K refers to the resistance coefficient of the pipeline.

Thermal Power Model
The solution of the thermal power system is mainly related to the following three temperatures: supply water temperature T s (the temperature from heat supply network to all thermal load node), outlet water temperature T o (outlet water temperature of the thermal load node), and return water The thermal of all nodes can be expressed as follows: Considering the heat loss of the pipeline, the temperature drop of the water flow at the first and last nodes of the pipeline during its transmission is as follows: The temperature at the joint of multiple pipes can be calculated according to Equation (10):

Model of Cooling Power Subsystem
The model of the cooling power subsystem is similar to that of thermal power subsystem, which also can be described by Equations (5)- (10). No further description is given here.

Model of Coupling Link
The coupled integrated energy station is adopted for the electricity-thermal-cooling coupling link, the energy conversion and flow relationship is shown in Figure 2. The thermal of all nodes can be expressed as follows: Considering the heat loss of the pipeline, the temperature drop of the water flow at the first and last nodes of the pipeline during its transmission is as follows: The temperature at the joint of multiple pipes can be calculated according to Equation (10):

Model of Cooling Power Subsystem
The model of the cooling power subsystem is similar to that of thermal power subsystem, which also can be described by Equations (5)- (10). No further description is given here.

Model of Coupling Link
The coupled integrated energy station is adopted for the electricity-thermal-cooling coupling link, the energy conversion and flow relationship is shown in Figure 2. Energy concentrator theory is adopted for obtaining the following mathematical model of the energy station: The thermal-to-electric ratio and thermal-to-cooling ratio of the energy station are obtained as follows, under the operation mode of determining electricity/cold through heat:  Energy concentrator theory is adopted for obtaining the following mathematical model of the energy station: Energies 2019, 12, 4455

of 16
The thermal-to-electric ratio and thermal-to-cooling ratio of the energy station are obtained as follows, under the operation mode of determining electricity/cold through heat: The conversion efficiency of all energy conversion components in the energy station and distribution coefficient of the energy flow are shown in Table 1. The thermal-to-electric ratio c m,e of the energy station is 0.54 through calculation, and the thermal-to-cooling ratio c m,c is 3.87.

Unified Power Flow Model and Calculation Algorithm
Modeling of each decoupling subsystem is the basis of power flow analysis of the IES. It is necessary to set up the system power flow solution model and select the appropriate power flow solution algorithm in order to realize the calculation of power flow of cooling-thermal-electric IES.

Unified Power Flow Model
The steady-state mixed power flow model of the IES with power subsystem, thermal power subsystem and cooling power subsystem can be expressed through the following form: The coupling relationship between above systems can be described by the coupling judgment matrix CO: Each element in the equation represents the mutual influence of the energy system. When each element of CO matrix is not equal to zero, the above three energy subsystems are coupled completely; they are dependent on mutual operation state. Coupling mode analysis of the IES is the basis of energy flow analysis. Different mixed power flow solutions should be considered for different coupling modes.

Solution Algorithm
The theoretical basis for solving the mixed power flow problem lies in that different energy networks all satisfy energy conservation law and material conservation law. However, a large amount of high-dimensional nonlinear equation and decision-making variable are added due to the simultaneous equation of different energy system, which makes the solution of the equations complicated. At present, the Newton-Raphson method is used as a more efficient algorithm for solving the unified power flow model of the IES [16,18].
The essence of power flow calculation is the process of obtaining the system operation point under a given series of conditions for the IES. Since the cooling power subsystem does not have a common node network, modeling of the IES is not conducted for cooling power subsystem. Thus, the unified power flow model of electric-thermal coupling system can be described as follows: where line 1 and 2, respectively, express the active deviation and reactive deviation of the electric power system, line 3-6, respectively, express node thermal power deviation of the thermal power subsystem, loop pressure drop deviation of the heat supply network, heat supply temperature deviation, and back heating temperature deviation, P, Q, and Φ are, respectively, the active power, reactive power, and thermal power of the system, A S refers to the reduced order incidence matrix of the heat supply network, C S and C r are the matrix, respectively, of the heat supply network and heat supply network structure and flux, B s and b r are column vectors related to heat supply temperature and output temperature, respectively, and x = [θ,U,m,T s,load ,T r,load ] T is the system state quantity.
The unified power flow model is solved through Newton-Raphson Method, and the iteration form is shown as follows: where ∆F refers to the deviation value of the power flow equation set, J refers to Jacobi matrix, x(k) and ∆x(k), respectively, refer to the state variable in the k iteration and deviation value of the state variable, and the superscript k represents iterations. Jacobi matrix J can be expressed as follows: where diagonal block J ee and J hh , respectively, refer to the relationship between own power flow of separate electric and thermal systems and own state variable and non-diagonal block J eh and J he , respectively, refer to the coupling relationship among the electric and thermal subsystem. The power flow equation set was solved simultaneously to get the current operation point of the system according to the initial conditions.

Example Calculation
The topology shown in Figure 1 and the modeling method in Section 2 were used to establish the unified power flow model for modeling and simulation, and the parameter list of electric power load and thermal power load is given as shown in A1 [19,20] and A2 [17] of Appendix A.
The pipeline flux, node supply water temperature and node return water temperature of the thermal power subsystem are obtained, as shown in Figures 3-5, respectively, after alternating iteration operation, and the electric power system bus voltage is shown in Figure 6.                  coupling link because the thermal-to-electric ratio of the cooling-thermal-electric coupling link is higher. The coupling link also can provide more electric power at the same time under the condition of meeting the thermal power requirements. Therefore, it has large support role to the electric power system bus voltage.
The load power supplied to the cooling power system is 103.39 kW by calculation through the thermal-to-cooling ratio.

Steady State Characteristic Analysis
Due to the close coupling relationship of the integrated energy system, the system model presents the characteristics of multi-input, multi-output, and dynamic change, which causes a change of a certain parameter in the system to directly affect the output characteristics and output state. Thus, it is necessary to carry out a sensitivity analysis and comparison of some parameters of the system. However, it is difficult to analyze the sensitivity change of the whole system, and it is simple and meaningful to analyze the sensitive parameters of a single system. Therefore, this paper will analyze some sensitive parameters of a single system.

Sensitivity Analysis of Power Subsystem
The sensitive parameters of the electric power subsystem are considered in this section, including the output of photovoltaic power generation. The output of the photovoltaic power generation system is affected by solar radiation intensity, photovoltaic array area, inverter efficiency, and other factors. The maximum output power is shown as follows: where ε refers to the photovoltaic conversion efficiency (%) of the photovoltaic array, A PV refers to the area (m 2 ) of the photovoltaic array, G T refers to solar radiation intensity (kW/m 2 ), and T C refers to the solar cell temperature ( • C). Therefore, the output of the photovoltaic power generation is mainly affected by solar radiation intensity G T . The change group of the solar radiation intensity can be estimated through statistics of previous solar radiation intensity data in the area, thereby obtaining the interval of the photovoltaic power generation output. A group of photovoltaic power generation device is added in bus 7 of the power subsystem as shown in Figure 7.
Energies 2019, 12, x FOR PEER REVIEW 8 of 17 Figure 6 shows that the voltage around power subsystem bus 19 is greatly increased compared, with 33 node power distribution network connected with thermal/cooling network not through the coupling link because the thermal-to-electric ratio of the cooling-thermal-electric coupling link is higher. The coupling link also can provide more electric power at the same time under the condition of meeting the thermal power requirements. Therefore, it has large support role to the electric power system bus voltage.
The load power supplied to the cooling power system is 103.39 kW by calculation through the thermal-to-cooling ratio.

Steady State Characteristic Analysis
Due to the close coupling relationship of the integrated energy system, the system model presents the characteristics of multi-input, multi-output, and dynamic change, which causes a change of a certain parameter in the system to directly affect the output characteristics and output state. Thus, it is necessary to carry out a sensitivity analysis and comparison of some parameters of the system. However, it is difficult to analyze the sensitivity change of the whole system, and it is simple and meaningful to analyze the sensitive parameters of a single system. Therefore, this paper will analyze some sensitive parameters of a single system.

Sensitivity Analysis of Power Subsystem
The sensitive parameters of the electric power subsystem are considered in this section, including the output of photovoltaic power generation. The output of the photovoltaic power generation system is affected by solar radiation intensity, photovoltaic array area, inverter efficiency, and other factors. The maximum output power is shown as follows: where ε refers to the photovoltaic conversion efficiency (%) of the photovoltaic array, APV refers to the area (m 2 ) of the photovoltaic array, GT refers to solar radiation intensity (kW/m 2 ), and TC refers to the solar cell temperature (°C). Therefore, the output of the photovoltaic power generation is mainly affected by solar radiation intensity GT. The change group of the solar radiation intensity can be estimated through statistics of previous solar radiation intensity data in the area, thereby obtaining the interval of the photovoltaic power generation output. A group of photovoltaic power generation device is added in bus 7 of the power subsystem as shown in Figure 7.    Table 2. Comparing the data in Table 2, the photovoltaic power generation is regarded as a distributed power supply and connected into the power subsystem. It plays a supporting role to the voltage of all buses, the output is lager, and the supporting role to the bus voltage is larger. Since the capacity of photovoltaic power generation set is smaller, the total power generation capacity can only meet the load demand of bus 7. However, a small part of electric load in bus 7 should be supplied by the power grid. The demand for electric power load of bus 7 is decreased from the perspective of power grid, thus, the voltage of all buses is increased to a small degree, and the voltage at bus 7 drops slowly.

Sensitivity Analysis of Thermal Power Subsystem
In this section, by adjusting the key technical parameters of the thermal subsystem and solving the steady-state power flow, the influence characteristics of different parameters on the system are analyzed. steady-state power flow under different source-node supply water temperatures was calculated. The pipeline steady-state flux, node supply water temperature and node return water temperature were obtained, as shown in the Figures 8-10.
In this section, by adjusting the key technical parameters of the thermal subsystem and solving the steady-state power flow, the influence characteristics of different parameters on the system are analyzed.

Temperature of Source-Node Supply Water
Four groups of different source-node supply water temperatures are set under the condition that the parameters of other technologies are consistent, including: 100 °C, 90 °C, 80 °C, and 70 °C. The steady-state power flow under different source-node supply water temperatures was calculated. The pipeline steady-state flux, node supply water temperature and node return water temperature were obtained, as shown in the Figures 8-10.  In this section, by adjusting the key technical parameters of the thermal subsystem and solving the steady-state power flow, the influence characteristics of different parameters on the system are analyzed.

Temperature of Source-Node Supply Water
Four groups of different source-node supply water temperatures are set under the condition that the parameters of other technologies are consistent, including: 100 °C, 90 °C, 80 °C, and 70 °C. The steady-state power flow under different source-node supply water temperatures was calculated. The pipeline steady-state flux, node supply water temperature and node return water temperature were obtained, as shown in the Figures 8-10.   Figure 10 shows that the source-node supply water temperature is lower and the steady-state flux in the pipeline is higher. The same conclusion also can be obtained by the thermal equation:   Figure 10 shows that the source-node supply water temperature is lower and the steady-state flux in the pipeline is higher. The same conclusion also can be obtained by the thermal equation: When the load thermal power Φ is constant and the node outlet water temperature T 0 is the same, the node supply water temperature T s is lower, the flux m q for the load is larger, and the pipeline flux in the heat supply network system is also larger. Figure 8 shows that the supply water temperature of all nodes is also reduced when the source-node supply water temperature is reduced. However, the relative relationship among all node supply water temperatures remains unchanged. Figure 9 shows that the return water temperature of all load nodes is still the same, with the set outlet water temperature under different source-node supply water temperatures, namely 50 • C. The return water temperature of the non-load node is slightly increased with the increase of the source-node supply water temperature, and the change trend is small.
Parameter changes of the thermal power subsystem have influence on the electric power subsystem through the coupling link. Since the thermal power subsystem scale is much smaller than the power subsystem, the source node supply water temperature changes in the thermal subsystem may generate almost no influence on the electric power subsystem. For example, it can be seen that when the water supply temperature of the source-node of the thermal subsystem changes, the steady flow of the pipeline and the return water temperature of the node have a great influence, but have little effect on the power subsystem. The output of the power subsystem and the cold subsystem cannot be changed under the operation mode of determining electricity with heat or determining cold with heat.

Pipeline Length of Heat Supply Network
All technical parameters were consistent, taking the parameters of the example in Section 4 as the standard, and with four groups of different pipe lengths of the heat supply network: 0.5 times standard length, 0.75 times standard length, standard length, and 1.25 times standard length. The steady-state power flow under different pipeline lengths of the heat supply network was calculated. Steady-state pipeline flux, node supply water temperature and node return water temperature were obtained, as can be seen in Figures 11-13.   It can be seen from Figures 11-13 that when the pipe length of the heat supply network increases, the steady-state pipe flow increases by a small margin, and the supply water temperature and the return water temperature of each node decrease by a small margin. This is because with the increase of the pipe length of the heating network, the loss of heat in the transmission process with water as    It can be seen from Figures 11-13 that when the pipe length of the heat supply network increases, the steady-state pipe flow increases by a small margin, and the supply water temperature and the return water temperature of each node decrease by a small margin. This is because with the increase of the pipe length of the heating network, the loss of heat in the transmission process with water as It can be seen from Figures 11-13 that when the pipe length of the heat supply network increases, the steady-state pipe flow increases by a small margin, and the supply water temperature and the return water temperature of each node decrease by a small margin. This is because with the increase of the pipe length of the heating network, the loss of heat in the transmission process with water as the medium increases, which is shown as the increase of temperature loss, as shown in the drop formula of water flow temperature at the first and last nodes of the pipe: The pipeline length L is larger, the temperature of the pipeline terminal is lower. Therefore, both supply water temperature and return water temperature of all nodes are decreased, and the supply water temperature is decreased more obviously. Therefore, the temperature difference of the thermal power supplied by the thermal-supply pipe network is decreased, and the pipeline flux is increased accordingly, thereby ensuring adequate supply of thermal load.
The changes of the heat supply network pipeline length in the thermal power subsystem may produce certain influence on steady-state pipeline flux in the thermal power subsystem as well as node supply water temperature and return water temperature. There is almost no influence on the power subsystem and cooling power subsystem.