Optimal Operation for Integrated Electricity–Heat System with Improved Heat Pump and Storage Model to Enhance Local Energy Utilization

As the need for clean energy increases, massive distributed energy resources are deployed, strengthening the interdependence of multi-carrier energy systems. This has raised concerns on the electricity-heat system’s co-operation for lower operation costs, higher energy efficiency, and higher flexibility. This paper discusses the co-operation of integrated electricity–heat system. In the proposed model, network constraints in both systems are considered to guarantee system operations’ security: the branch flow model is utilized to describe the electricity network, while a convexified model considering variable mass flow and temperature dynamics is adopted to describe the heat network. Additionally, novel models for heat pumps and the stratified water tank are proposed to represent the physical system more accurately. Finally, to preserve the information privacy of separate systems, a distributed algorithm is proposed based on the alternating direction method of multipliers (ADMM). Numerical studies show that the co-operation could provide a more economical and reliable solution than the decoupled operation of the heat network and electricity network. Moreover, the ADMM-based algorithm could derive solutions very close to the optimum provided by centralized optimization.


Introduction
It is a non-stoppable trend that civilization would gradually transit to a greener and more sustainable society. This increasing need has driven the rapid installation of renewable energy integrated in power systems. Taking distributed photovoltaic power generation in China as an example, its newly installed capacity has reached 12.2 GW, a year-on-year increase of 41.3% [1], in 2019. To enhance renewable energy integration, the power system faces multiple challenges [2]. Moreover, the mismatching of demands and supplies in local energy systems is getting worse, which might harm the security of electricity system operation in the future. People have taken measures to enhance system flexibility, such as deploying battery storage systems and demand management. As a result, the increase in the expense of electricity system operation and planning is inevitable. For instance, even battery cost has rapidly decreased for decades, the unsubsidized levelized cost of battery storage in commercial and industrial is still up to $432-590 per MWh [3].
Apart from planning and managing the electricity system solely, integrating multi-energy systems provides a new solution to harnessing heterogeneous energy resources efficiently. To our knowledge, storage of thermal energy is more applicable while the delivery of electricity is of higher efficiency. For northern China, a large amount of thermal generation units produce electricity as well as heat in cold winters. Furthermore, as electricity substitution takes place in many places in China, electricity-driven

•
The nonlinear characteristics of the energy conversion units are addressed. For modeling heat pumps, the relationship of heat production and power consumption is expressed by a quadratic function and transformed to conic constraints by the second-order conic relaxation. Additionally, a new expression of stratified heat storage water tank is proposed to be incorporated in the proposed network model. • A network constrained optimal operation model for the integrated power and heat system is then proposed based on mixed-integer second-order conic programming (MISOCP). For the operation of PDN, ACOPF is modeled and described by the branch flow model. For the operation of DHN, a convex description of the thermal network based on conic relaxation is applied, in which water flow could be adjusted in different time slots while considering temperature dynamics. • For a decentralized operation between thermal network and electricity network without a centralized operator, an ADMM based approach is proposed to reconstruct the centralized co-operation model into separated network models run by independent heat system operator (HSO) and distribution system operator (DSO). In the proposed structure for decentralized optimization, HSO and DSO can independently operate DHN and PDN with limited information exchange. As shown by the numerical study, a near-optimal solution could be found through iterations, which is very close to the centralized optimization result.
The rest of this paper is organized as follows. Preliminaries and notations of the integrated energy system are introduced in Section 2. The proposed co-operation model is described in detail in Section 3. In Section 4, the convexification method of DHN is introduced. Additionally, the distributed algorithm for the co-operation model based on ADMM is proposed. Case studies are presented in Section 5. Section 6 goes into a more in-depth discussion about the rationality and effectiveness of the proposed model. Finally, the conclusion is given in Section 7 by highlighting the essential findings of our work.

Preliminaries and Notations
In this work, we concentrate on the optimal operation of integrated electricity and heat systems considering network constraints.

Integrated Electricity and Heat System
The details of the integrated system would be discussed in this subsection. As is mentioned in this paper, the integrated system consists of two sub-network system: the PDN and DHN. The paradigm of the system is illustrated in Figure 1.

Preliminaries and Notations
In this work, we concentrate on the optimal operation of integrated electricity and heat systems considering network constraints.

Integrated Electricity and Heat System
The details of the integrated system would be discussed in this subsection. As is mentioned in this paper, the integrated system consists of two sub-network system: the PDN and DHN. The paradigm of the system is illustrated in Figure 1.
In the PDN, the electricity is supplied from the main grid and distributed generation, such as fuel-fired CHPs and renewable power plants. In the DHN, heat demand is supplied by heat pumps and CHPs. Storage systems and flexible demands are considered as flexible resources for the integrated system.  CHP units act as the source of both electricity and heat system. There are mainly two types of CHP, i.e., back-pressure CHP units and condensing CHP units [4]. The back-pressure CHP units are most widely used to deliver power and heat, in which all input steam is utilized to supply heat after exhaustion. Besides, back-pressure CHPs usually have a smaller capacity and are less flexible, since the heat-electricity ratio is fixed. In our paper, the subject of our research is the local integrated electricity-heat system, where back-pressure CHP units are commonly used.
Heat pumps play an increasing role in the heating sector [6]. A heat pump acts as electricity demand in the PDN system, while a heat source in the DHN. Electricity is utilized for heat pumps to transfer the heat from a colder space to warmer one by a vapor compression cycle. The heat could be taken from sources like air, ground heat, and water, etc.
Renewable power plants, such as PVs and wind turbines, could generate electricity power at a very low marginal price. In the PDN, distributed PV stations and rooftop PV units are commonly used in the distribution power system communities. Thus, later in our model, only PV units are considered.
A stratified thermal storage water tank is typical to heat storage [9], which is widely used in Europe practically. In this type of water tank, hot water is withdrawn (or injected) from the top of the tank, while the cold water would be injected (or withdrawn) into the bottom of the tank to keep the volume constant. Physically, the hotter water is with lower density and thus naturally stays in the upper part of the tank, and the colder water would stay in the lower part of the tank. Thus, if well controlled, the water tank would be stratified and able to store thermal energy. However, due to the In the PDN, the electricity is supplied from the main grid and distributed generation, such as fuel-fired CHPs and renewable power plants. In the DHN, heat demand is supplied by heat pumps and CHPs. Storage systems and flexible demands are considered as flexible resources for the integrated system. CHP units act as the source of both electricity and heat system. There are mainly two types of CHP, i.e., back-pressure CHP units and condensing CHP units [4]. The back-pressure CHP units are most widely used to deliver power and heat, in which all input steam is utilized to supply heat after exhaustion. Besides, back-pressure CHPs usually have a smaller capacity and are less flexible, since the heat-electricity ratio is fixed. In our paper, the subject of our research is the local integrated electricity-heat system, where back-pressure CHP units are commonly used.
Heat pumps play an increasing role in the heating sector [6]. A heat pump acts as electricity demand in the PDN system, while a heat source in the DHN. Electricity is utilized for heat pumps to transfer the heat from a colder space to warmer one by a vapor compression cycle. The heat could be taken from sources like air, ground heat, and water, etc.
Renewable power plants, such as PVs and wind turbines, could generate electricity power at a very low marginal price. In the PDN, distributed PV stations and rooftop PV units are commonly used in the distribution power system communities. Thus, later in our model, only PV units are considered.
A stratified thermal storage water tank is typical to heat storage [9], which is widely used in Europe practically. In this type of water tank, hot water is withdrawn (or injected) from the top of the tank, while the cold water would be injected (or withdrawn) into the bottom of the tank to keep the volume constant. Physically, the hotter water is with lower density and thus naturally stays in the upper part of the tank, and the colder water would stay in the lower part of the tank. Thus, if well controlled, the water tank would be stratified and able to store thermal energy. However, due to the heat conduction on the water interface with different temperatures, there would be a third layer, i.e., the thermocline layer, between the hot-water and the cold-water layers. In the thermocline layer, the temperature changes quickly. As time goes, this layer would gradually become thicker, degrading the tank's valid heat capacity. Figure 2 shows the structure of this type of water tank and its connection with DHN.
Energies 2020, 13, x FOR PEER REVIEW 5 of 23 heat conduction on the water interface with different temperatures, there would be a third layer, i.e., the thermocline layer, between the hot-water and the cold-water layers. In the thermocline layer, the temperature changes quickly. As time goes, this layer would gradually become thicker, degrading the tank's valid heat capacity. Figure 2 shows the structure of this type of water tank and its connection with DHN. The modeling of DHN and PDN involves plenty of variables and parameters. For clarity, the following subsection would introduce the basic background knowledge for these networks and the important notations utilized in the problem formulation.

Notations of District Heating Network
In the subsection, the main notations of the adopted DHN model would be introduced. The modeling of DHN is based on the description in [20]. The pipeline structure of the DHN is presented in Figure 3. The DHN consists of two sub-systems: the supply sub-system and the return sub-system. Warm water flows through the supply sub-system from source nodes to demand nodes for delivering heat to end-users. In contrast, the used cool water returns from demand nodes and flows back to source nodes through the return sub-system. , Figure 3. Structure of the district heating network. The modeling of DHN and PDN involves plenty of variables and parameters. For clarity, the following subsection would introduce the basic background knowledge for these networks and the important notations utilized in the problem formulation.

Notations of District Heating Network
In the subsection, the main notations of the adopted DHN model would be introduced. The modeling of DHN is based on the description in [20]. The pipeline structure of the DHN is presented in Figure 3. The DHN consists of two sub-systems: the supply sub-system and the return sub-system. Warm water flows through the supply sub-system from source nodes to demand nodes for delivering heat to end-users. In contrast, the used cool water returns from demand nodes and flows back to source nodes through the return sub-system. heat conduction on the water interface with different temperatures, there would be a third layer, i.e., the thermocline layer, between the hot-water and the cold-water layers. In the thermocline layer, the temperature changes quickly. As time goes, this layer would gradually become thicker, degrading the tank's valid heat capacity. Figure 2 shows the structure of this type of water tank and its connection with DHN. The modeling of DHN and PDN involves plenty of variables and parameters. For clarity, the following subsection would introduce the basic background knowledge for these networks and the important notations utilized in the problem formulation.

Notations of District Heating Network
In the subsection, the main notations of the adopted DHN model would be introduced. The modeling of DHN is based on the description in [20]. The pipeline structure of the DHN is presented in Figure 3. The DHN consists of two sub-systems: the supply sub-system and the return sub-system. Warm water flows through the supply sub-system from source nodes to demand nodes for delivering heat to end-users. In contrast, the used cool water returns from demand nodes and flows back to source nodes through the return sub-system.  Mass flow rate is used to describe the mass of water passing through a pipe per unit of time. At each node, the injected mass flow rate m j,t is defined. Note that the nodal mass flow rate m j,t Energies 2020, 13, 6729 6 of 23 is positive for source nodes, while negative for demand nodes. T S j,t and T R j,t represent the supply and return temperature at node j.
The network connects all nodes through pipes. In these pipes, water is pushed forward because of the pressure difference between source and demand nodes. φ j,t is defined as the nodal pressure. As for the pipe from i to j, the mass flow rate in this pipe is denoted by m ij,t . Note that the direction of water flow in the pipe is not fixed due to the variation of demand and source profiles. Without loss of generality, y ij,t is introduced to define the direction of water flow. y ij,t = 1 means the water flow from node i to node j (i.e., the predefined direction), while y ij,t = 0 denotes the opposite situation. In the supply sub-system, temperatures at the inlet and outlet ends of the pipe i j are denoted by T S,i ij,t and T S,o ij,t , respectively, while T R,i ij,t and T R,o ij,t are counterparts for the corresponding pipe in the return sub-system. Finally, water flows from several pipes would mix into one node, as shown by node j in Figure 3. Thus, T Sm j,t and T Rm j,t are introduced to represent the mixing temperatures at node j for supply and return sub-systems, respectively.

Notations of Power Distribution Network
In our study, PDN operations are modeled as an AC-OPF problem by branch flow model [16]. P mn,t and Q mn,t denote active and reactive power flows along the distribution line mn with the resistance r mn and the reactance x mn . The squared current magnitude of this line mn at time t is denoted by mn,t . The squared voltage magnitude of bus m at time t is denoted by v m,t . p net m,t and q net m,t represent the net active and reactive power injected into bus m at time t.

Problem Formulation
In this section, the centralized co-operation model of the DHN-PDN system is proposed. This scheme assumes an independent operator monitoring and managing the integrated system. All information needed would be gathered and sent to the ISO. The objective of this model would benefit the whole system considering both sides of the costs and revenues, while it meets the operational constraints of both PDN and DHN.

Objective
The goal of the optimal operation model of this integrated system is to minimize the total operation cost in both power and heat systems. In the PDN, the operation cost consists of power generation costs from renewable power plants and CHPs. While in the DHN, the operation cost consists of heat generation cost from CHPs and thermal storage devices. Finally, the revenue from flexible heat demands is also considered in this objective. The objective function is presented in Equation (1).
where λ Grid denotes the grid supply electricity price. C E,PV r is the unit cost of PV unit, which could be rather a small value or 0. C E,CHP r and C H,CHP r are unit cost of CHP unit for generating electricity and heat. C HS r denotes the operation cost of thermal storage tank.

Component Constraints
•

Combined Heat and Power
As in [4], the power production of back-pressure CHP depends on its heat production. Thus, this dependence is modeled by linear constraints (2). The complete constraints of CHP are shown in Equations (2)-(4).

Heat pumps
To quantify the utilization efficiency of the heat pump, COP is defined as the ratio of heat production to electricity consumption. Previous models treat the COP of heat pumps as a constant value. For example, COP = 4 means that 1 unit of electricity could generate 4 units of heat. In fact, as the heat demand increases, the COP would decline. This means the economic efficiency of heat pumps is decreasing as their demand grows.
However, some experimental results indicate that the COP is not constant, thus requiring more precise models for heat pumps. In fact, the COP would decrease as heat production increases [5]. This relationship could be modeled by the quadratic function, as shown in Equation (5).
The above Equation can be convexified by relaxing the equality into second-order conic constraints Equation (6). Then, boundary constraints for heat pumps are given in Equations (7) and (8).
, and d r = −C r + 1+b r 2 2 . According to Equation (5), COP could be expressed as the Equation (9), where A, B, and C are unknown parameters. By proper fitting techniques for the parameters in Equation (6), we could acquire a more accurate expression of heat pumps. For instance, the experimental data from [6] are then used to fit these unknown parameters of heat pumps, as shown in Figure 4, which would definitely significant difference in optimal operation results. For example, in the constant COP model, a heat pump with COP = 4 needs 150 kWh of power consumption to produce 600 kWh heat, while in the variable COP model, 252.9 kWh of electricity is acquired at COP = 2.37.
Energies 2020, 13, x FOR PEER REVIEW 8 of 23  The constraint on PV's output is simple, as shown in (10). At each time slot, the capability of PV's output could be denoted by p HP r,t . However, due to the limitation of network constraints and the balancing of the demand and supply, the operator sometimes needs to curtail the output of PV to p HP r,t , in order of reliable and economical operation for the PDN.
• Heat Storage As mentioned above, hot and cool water could be stored simultaneously in a stratified thermal storage water tank. However, the heat loss between two layers of water could not be neglected. Equations (11)- (14) are constraints for heat storage operation, where η HS,l denotes the loss in the thermocline layer, while η HS denotes the efficiency of charging and discharging in the water tank.
To incorporate the storage model into the network model, z HS r indicating the state of charging and discharging is introduced.

Operation Constraints of DHN Model
Constraints Equations (15)-(37) describe the operation of DHN. Equations (15) and (16) illustrate the relationship of the net heat injection at each node with the nodal status (i.e., the temperature difference and the mass flow rate) in DHN. c w is the constant of heat storage in the water (c w = 4.2 kJ/(kg·K)), and ∆T j,t is the temperature difference between supply end and return end at node j at time t.
for storage nodes.
The dynamic of DHN operations can be modeled in two processes: the hydraulic process and the thermal dynamic process. First, the water flow in heat networks should satisfy the following hydraulic constraints. Here, Equation (17) shows the conservation of water flow in the network. At each node, all inflow mass flow rates should always equal to all outflow rates. m ij,t is the mass flow rate along the pipe i j, while m j,t represents the injected mass flow rate at node j. Equations (18) and (19) limit the pressure loss in pipes, where φ j,t is the water pressure at node j and κ ij,t is the pressure loss along pipe i j. Equations (20) and (21) are bounds for the hydraulic variables.
Energies 2020, 13, 6729 Temperature dynamics constraints of DHN operations are presented in Equations (23)-(37). Equations (23) and (24) describe the mixing temperature at node j. Here, we assume that all flows are thoroughly mixed at one node before flowing out. Thus, the mixing temperature T Sm j,t at supply sub-system in node j is dominated by all inflows weighted by their mass flow rates. T S,o ij,t is the temperature at the outlet of the pipe i j in the supply network, T S j,t is the nodal temperature at node j in the supply network.
Note that, as for source nodes, the flow direction is from the return pipe network to the supply pipe network, while load nodes are in the opposite situation. Therefore, the mass flow of the source node is treated as the inflow in the supply network, while the mass flow of the load node is treated as the inflow in the return network, which is expressed in Equations (25) and (26). m S j,t and m R j,t are the nodal mass flow rates injected into the demand node and source node, respectively.
Then, the nodal mixing temperature will dominate the temperatures of all outflows from this node. Whether in either supply or return network, the temperature of the pipeline inflow is dominated by the mixing temperature of its starting node. This principle is expressed in Equations (27) and (28), where Equation (28) denotes the reverse flow situation. When it comes to nodal flows, for source nodes, the nodal return flow is an outflow, while for demand nodes, the nodal supply flow is an outflow. Therefore, as in Equations (29) and (30), respectively, the mixed temperature in the return network dominates nodal return temperature at source nodes, while the mixed temperature in the supply network dominates nodal supply temperature at demand nodes. Similarly, Equation (31) denotes the temperature domination rules in nodes with storage. For nodes with heat storage, the node becomes a source node when heat storage is in discharging mode, while a demand node when in charging mode.
Last but not least, the temperature loss in the pipelines depends on the ambient temperature T a t , the pipe mass flow rate m ij,t , the heat conductance ij and the pipe length L ij , which is modeled in Equations (32) and (33).
Finally, constraints Equations (34)-(37) set upper and lower bounds of temperatures from all nodes, pipes, and mixing points. Thus, DHN can operate in a feasible region.

Operation Constraints of PDN Model
The PDN is a radial network that could be precisely modeled by the branch flow model [16]. All operation constraints for PDN are presented in Equations (38)-(43). Equations (38) and (39) are nodal active and reactive power balance constrains, respectively.
where other than the root bus, the net active power is given as p net n,t = l∈P CHP v n,t = v m,t − 2(r mn P mn,t + x mn Q mn,t ) + r 2 mn + x 2 mn mn,t , ∀n ∈ P (40) mn ≤ mn,t ≤ mn , ∀n ∈ P, ∀m ∈ A n (43) Energies 2020, 13, 6729 11 of 23

Convexification of Our Proposed Operation Model
The proposed model in Section 2 is highly nonlinear if mass flow rates and temperature profiles are both treated as variables. In [22], a convexification method is proposed to address this challenge. In this section, this method would be introduced to help convexify the nonlinear operation model into a mixed-integer second-order conic programming (MISOCP).

•
Heat and its Balance Heat nodal balance constraints Equation (15) are bilinear if both the temperature difference and the nodal mass flow rate are seen as variables. In practice, the temperature difference ∆T j,t could be fixed as a constant, such as 40 • C at demand nodes and 41-43 • C for source nodes, depending on the network insulation level. As a result, Equation (15) is linearized for all nodes in DHN.

Pressure Loss
The pressure loss along pipes satisfies the constraint Equation (22), which could be easily transferred to the SOC constraint by relaxing the equality and transfer the expression to Equation (45). Then, it could be written in the standard form of SOC constraint in Equation (46). •

Heat Loss
Theoretically, in Equations (47) and (48), the heat loss along a pipe could be expressed as the product of the heat capacity of water c w , the absolute value of mass flow rate M ij,t , and the absolute temperature loss ∆T S ij,t along the pipe.
Due to the reverse water flow in the heat network, a set of restrictions are introduced in Equations (49)-(54) of limiting ancillary variables M ij,t , ∆T S ij,t , and ∆T R ij,t to represent the absolute value.
Note that the temperature loss could be calculated from Equations (32) and (33). The nonlinear exponential function in Equations (32) and (33) can be linearized as Equations (55) and (56). Finally, the heat loss could be simplified as a constant if the inlet pipe temperature is substituted by a reference temperature. Usually, the inlet temperature would vary in a small range. If this reference temperature is carefully chosen, the heat loss can be simplified without losing much accuracy. Later in Section 4, we would go into more detail to verify the conclusion made here by numerical studies.
Thus, the heat loss Equation could be simplified into the following bilinear constraints Equations (57) and (58). Further, the second-order constraints after relaxing these Equations are shown in Equations (59) and (60). •

Temperature Mixing
The most complicated constraints in DHN operations are those related to temperature mixing, where plenty of bilinear terms are introduced, which is of great difficulty to model. To simplify this complicated constraint, we assume that the nodal mixing temperature would be dominated by the inlet temperature with the largest mass flow rate [22]. In the supply network, the mixing temperature might be dominated by the outlet temperature of a certain pipe carrying the water flow into this node. Otherwise, the mixing temperature could be dominated by the nodal temperature for some source nodes. Similarly, in the return network, the dominating flow might flow from pipes or the nodal flows at load nodes. Here, z jt , z + ij,t , z − jk,t are binary variables indicating three situations mentioned above for node j in the supply network, respectively. z jt = 1 means the nodal flow is dominating. As for the pipe flows, z + ij,t = 1 means the flow in pipe i j is largest inflow of node j, while z − jk,t = 1 indicates the largest inflow of node j is the flow in pipe jk. Note that according to the predefined direction, the flow in pipe i j would leave node j unless the flow is reverse (y jk,t = 0). In the return network, z jt , z + ij,t , z − jk,t are also defined. Additionally, the maximal mass flow rates dominating the node in supply and return networks are defined as m S j,t and m R j,t .
Mixing temperature constraints could be simplified as the above constraints Equations (61)-(69). Equations (61)-(63) ensure that the maximal mass flow rates are not always smaller than all inflows. Equations (64)-(66) guarantee that the selected flow is the largest inflow. Once the dominating flow is selected, if it is a pipe flow, Equations (67) and (68) ensure the mixed temperature equals its outlet temperature. If it is a nodal flow, Equations (69) ensures the nodal supply temperature coincides with the mixing nodal temperature. All identical principles are also applied to the return network. For the sake of simplicity, mixing temperature constraints in the return network are not presented here.

ADMM-Based Decentralized Operation
In practice, HSO and DSO usually do not belong to the same utility. The co-operation of DHN and PDN might be impossible to implement due to the requirement of preserving the privacy information of both HSO and DSO. With the development of distributed computing techniques, the above co-operation problem can be solved separately by DSO and HSO. More importantly, only very limited and non-sensitive information is shared between DHN and PDN. ADMM might be the most popular one of distributed optimization techniques, which is thus also applied in our work. Our proposed co-operation model is summarized as follows.  − r∈D U D r h FD r,t . To apply ADMM, the above co-operation problem is first separated into two sub-problems locally optimized by DSO and HSO. Note that (2) and (6) are coupling constraints to combine these two sub-problems. Nevertheless, the standard ADMM requires the coupling constraints to be affine. While in our model, energy conversion constraints (6) for heat pumps is conic, making the standard ADMM non-applicable. To address this issue, ancillary variables p HP r,t are introduced as local copies for electricity consumption of heat pumps as in Figure 5.  ADMM non-applicable. To address this issue, ancillary variables , are introduced as local copies for electricity consumption of heat pumps as in Figure 5. The ADMM-based decentralized algorithm to solve our proposed model is presented as follows. After writing the augmented Lagrangian function, the overall optimization problem could be divided into three sub-problems, which can be solved separately and iteratively.  The ADMM-based decentralized algorithm to solve our proposed model is presented as follows. After writing the augmented Lagrangian function, the overall optimization problem could be divided into three sub-problems, which can be solved separately and iteratively. First, we define local variable sets X E = p CHP , p HP , p PV , p Grid , Ψ PDN for PDS and X H = h CHP , h HP , h HS+ , h HS− , E HS , p HP , Ψ DHN for DHS. Ψ PDN and Ψ DHN are variables related to PDN and DHN, respectively. Note that all local variables are initialized with a value at the first iteration. For the iteration k + 1, three following sub-problems would be solved sequentially.
Sub-Q1 is the PDN operation sub-problem solved by the DSO. It would update primal variables of PDN given the latest coupling variables from HSO.
Finally, after rounds of iteration, the solution would coverage to a consensus point, where both systems find the optimal solution for the integrated system.

System Configuration
To verify the contributions of our work, an integrated electricity and heat system is studied. This test system includes an IEEE standard PDN test system of 33 buses and a district heating system of 32 nodes and 32 pipes. The topology of this test system is illustrated in Figure 6. There are multiple energy sources in this integrated system, including two PV units, three back-pressure CHPs, three heat pumps, and three district stratified water tanks. All characteristics of the integrated system could be found in the Supplementary Materials. The test system is in a domestic community. The electricity demand has two peaks per day, while the heat demand has only one peak at night. Note that the power generation profiles of some intermittent resources like PVs are inconsistent with the electricity demand profile in this community. This would lead to significant needs for flexible resources to facilitate local energy consumption and improve utility efficiency. In this section, we carry out several tests to illustrate our proposed co-operation model's validity and effectiveness.
Energies 2020, 13, x FOR PEER REVIEW 15 of 23 Finally, after rounds of iteration, the solution would coverage to a consensus point, where both systems find the optimal solution for the integrated system.

System Configuration
To verify the contributions of our work, an integrated electricity and heat system is studied. This test system includes an IEEE standard PDN test system of 33 buses and a district heating system of 32 nodes and 32 pipes. The topology of this test system is illustrated in Figure 6. There are multiple energy sources in this integrated system, including two PV units, three back-pressure CHPs, three heat pumps, and three district stratified water tanks. All characteristics of the integrated system could be found in the Supplementary Materials. The test system is in a domestic community. The electricity demand has two peaks per day, while the heat demand has only one peak at night. Note that the power generation profiles of some intermittent resources like PVs are inconsistent with the electricity demand profile in this community. This would lead to significant needs for flexible resources to facilitate local energy consumption and improve utility efficiency. In this section, we carry out several tests to illustrate our proposed co-operation model's validity and effectiveness.  Figure 6. Topology of the test system.

Optimal Results of DHN-PDN Co-Operation
The proposed model's optimal operation results are presented in Figures 7 and 8, with the daily operation cost for the integrated energy system is RMB 94826. The operation results in the PDN are illustrated in Figure 7. As we can see, at peak hours of PV production, the electricity demand is in the valley, causing a significant imbalance of local generation and consumption. Compared to PVs, CHPs with higher operation cost are ramped down to leave room for the cheaper renewable generation. Additionally, heat pumps increase the power consumption to help the local consumption, transmitting this clean electricity power into the heating system. Note that the increase of heat pump consumption would cause a declination of the COP. Thus, the co-operation makes a tradeoff between higher COP and larger renewable energy consumption. Note that the hyphens in front of the numbers on the axes in the following figures represent minus signs, which are created by OriginLab 2019b.
The operation result of heat production and storage is illustrated in Figure 8, where "H_TSc" and "H_TSd" denote the sum-up heat storage. As we know, back-pressure CHPs in our test system have a linear relationship between power production and heat production. Thus, in peak hours of renewable energy generation, the ramping down of CHPs would cause the reduction of heat

Optimal Results of DHN-PDN Co-Operation
The proposed model's optimal operation results are presented in Figures 7 and 8, with the daily operation cost for the integrated energy system is RMB 94826. The operation results in the PDN are illustrated in Figure 7. As we can see, at peak hours of PV production, the electricity demand is in the valley, causing a significant imbalance of local generation and consumption. Compared to PVs, CHPs with higher operation cost are ramped down to leave room for the cheaper renewable generation. Additionally, heat pumps increase the power consumption to help the local consumption, transmitting this clean electricity power into the heating system. Note that the increase of heat pump consumption would cause a declination of the COP. Thus, the co-operation makes a tradeoff between higher COP and larger renewable energy consumption. Note that the hyphens in front of the numbers on the axes in the following figures represent minus signs, which are created by OriginLab 2019b. production at the same time, while heat pumps consume more power, causing an increase in heat production. Then, the heat surplus or deficit in the whole day is balanced by heat storage systems. In a word, the optimization model would figure out an optimal portfolio of flexible resources, where the heating system consumes a portion of the renewable output through the coupling.  It is known that COP drops as the power consumption deviates from the nominal value. However, at peak hours of PV generation, renewable plants usually generate electricity at a low marginal cost. As a result, operation optimization should make a tradeoff between the lower efficiency and the maximal consumption of cheaper energy. Figure 9 shows the optimal operation results for heat pump "HP3" and changes of COP during the operation day. The operation of HP3 tries to consume more power in the daytime, while at another time when PVs are off, heat pumps would return to the optimal operating point with relatively high efficiency. Energies 2020, 13, x FOR PEER REVIEW 16 of 23 production at the same time, while heat pumps consume more power, causing an increase in heat production. Then, the heat surplus or deficit in the whole day is balanced by heat storage systems. In a word, the optimization model would figure out an optimal portfolio of flexible resources, where the heating system consumes a portion of the renewable output through the coupling.  It is known that COP drops as the power consumption deviates from the nominal value. However, at peak hours of PV generation, renewable plants usually generate electricity at a low marginal cost. As a result, operation optimization should make a tradeoff between the lower efficiency and the maximal consumption of cheaper energy. Figure 9 shows the optimal operation results for heat pump "HP3" and changes of COP during the operation day. The operation of HP3 tries to consume more power in the daytime, while at another time when PVs are off, heat pumps would return to the optimal operating point with relatively high efficiency. The operation result of heat production and storage is illustrated in Figure 8, where "H_TSc" and "H_TSd" denote the sum-up heat storage. As we know, back-pressure CHPs in our test system have a linear relationship between power production and heat production. Thus, in peak hours of renewable energy generation, the ramping down of CHPs would cause the reduction of heat production at the same time, while heat pumps consume more power, causing an increase in heat production. Then, the heat surplus or deficit in the whole day is balanced by heat storage systems. In a word, the optimization model would figure out an optimal portfolio of flexible resources, where the heating system consumes a portion of the renewable output through the coupling.
It is known that COP drops as the power consumption deviates from the nominal value. However, at peak hours of PV generation, renewable plants usually generate electricity at a low marginal cost. As a result, operation optimization should make a tradeoff between the lower efficiency and the maximal consumption of cheaper energy. Figure 9 shows the optimal operation results for heat pump "HP3" and changes of COP during the operation day. The operation of HP3 tries to consume more power in the daytime, while at another time when PVs are off, heat pumps would return to the optimal operating point with relatively high efficiency.

Comparison of Different Operation Models
To investigate the benefits of co-operation for DHN-PDN, a group of comparison is carried out. In the decoupled operation (DO) case, the power system and heat system are run by their operators separately, where the power consumption price of heat pumps are fixed and predefined (set to 0.5 RMB/kWh). Additionally, four scenarios of different PV output profiles are applied to analyze the benefit of co-operating DHN and PDN.  Table 1 presents the comparison of operation costs and renewable energy consumption in two different cases. In this table, we define the cost reduction (CR) ratio as the relative difference in the total operation costs in two cases. It is evident that the total operation cost is lower in Case I compared to Case II by 2.3% to 4.4% according to different scenarios. The reason is: when the systems are operated separately in Case II, the operation cost of DHN decreases, since the flexibility in the DHN is not utilized to facilitate renewable energy utilization in the PDN. However, as a result, due to the lack of flexibility in the PDN, the decoupled operation case is unable to make the best use of local renewable resources, causing the expense of PDN to rise, which ends up a higher total operation cost

Comparison of Different Operation Models
To investigate the benefits of co-operation for DHN-PDN, a group of comparison is carried out. In the decoupled operation (DO) case, the power system and heat system are run by their operators separately, where the power consumption price of heat pumps are fixed and predefined (set to 0.5 RMB/kWh). Additionally, four scenarios of different PV output profiles are applied to analyze the benefit of co-operating DHN and PDN.

Comparison of Different Operation Models
To investigate the benefits of co-operation for DHN-PDN, a group of comparison is carried out. In the decoupled operation (DO) case, the power system and heat system are run by their operators separately, where the power consumption price of heat pumps are fixed and predefined (set to 0.5 RMB/kWh). Additionally, four scenarios of different PV output profiles are applied to analyze the benefit of co-operating DHN and PDN.  Table 1 presents the comparison of operation costs and renewable energy consumption in two different cases. In this table, we define the cost reduction (CR) ratio as the relative difference in the total operation costs in two cases. It is evident that the total operation cost is lower in Case I compared to Case II by 2.3% to 4.4% according to different scenarios. The reason is: when the systems are operated separately in Case II, the operation cost of DHN decreases, since the flexibility in the DHN is not utilized to facilitate renewable energy utilization in the PDN. However, as a result, due to the lack of flexibility in the PDN, the decoupled operation case is unable to make the best use of local renewable resources, causing the expense of PDN to rise, which ends up a higher total operation cost  Table 1 presents the comparison of operation costs and renewable energy consumption in two different cases. In this table, we define the cost reduction (CR) ratio as the relative difference in the total operation costs in two cases. It is evident that the total operation cost is lower in Case I compared to Case II by 2.3% to 4.4% according to different scenarios. The reason is: when the systems are operated separately in Case II, the operation cost of DHN decreases, since the flexibility in the DHN is not utilized to facilitate renewable energy utilization in the PDN. However, as a result, due to the lack of flexibility in the PDN, the decoupled operation case is unable to make the best use of local renewable resources, causing the expense of PDN to rise, which ends up a higher total operation cost for the integrated system. Nevertheless, co-operation as shown in Case I would utilize the flexibility of both systems to consume as much renewable energy output as possible, because these resources are of the lowest marginal costs compared to other resources like gas turbines, which results in a reduction in total operation cost. The decentralized co-optimization through ADMM is studied in this part. In this case, the SOs separately operate two energy systems but negotiate and update the scheduling iteratively. In Figure 11, we could see that the optimization gap drops quickly in the first several iterations. There would be a minimal gap (≤0.1%) after 50 iterations, and the obtained solution would be very similar to the centralized one. As the algorithm iterates, both systems reach a consensus, where the solution converges.
Energies 2020, 13, x FOR PEER REVIEW 18 of 23 for the integrated system. Nevertheless, co-operation as shown in Case I would utilize the flexibility of both systems to consume as much renewable energy output as possible, because these resources are of the lowest marginal costs compared to other resources like gas turbines, which results in a reduction in total operation cost. The decentralized co-optimization through ADMM is studied in this part. In this case, the SOs separately operate two energy systems but negotiate and update the scheduling iteratively. In Figure  11, we could see that the optimization gap drops quickly in the first several iterations. There would be a minimal gap (≤0.1%) after 50 iterations, and the obtained solution would be very similar to the centralized one. As the algorithm iterates, both systems reach a consensus, where the solution converges. Decentralized operation case would highly preserve the privacy information for local operators. Taking our proposed model as an example, coupling variables in the decentralized optimization case are , , , , ℎ , and , , which are only public information shared by each other. For DSO and HSO, they do not need to share network parameters of two systems during the optimization process. While in the case of centralized optimization, all information should be accessible to two independent system operators of the integrated system.

Discussion
From the case study, it is evident that the proposed co-operation model and the ADMM-based decentralized algorithm could facilitate local energy accommodation and increase the total revenue of the local energy system. However, to efficiently solve the network constrained operation problem, several simplifications have been made. In this section, we would make a deeper discussion on the rationality of the assumptions and figure out the effectiveness of the proposed model.

Rationality of Assumptions in the Network Model
As is shown in Section 3, some assumptions are made to reformulate the operation model as a MISOCP form. To verify the rationality of these assumptions made in our proposed operation model, we compare this model with the CFVT model proposed in [20]. The CFVT model is widely used in Decentralized operation case would highly preserve the privacy information for local operators. Taking our proposed model as an example, coupling variables in the decentralized optimization case are p CHP r,t , p HP r,t , h CHP r,t and p HP r,t , which are only public information shared by each other. For DSO and HSO, they do not need to share network parameters of two systems during the optimization process. While in the case of centralized optimization, all information should be accessible to two independent system operators of the integrated system.

Discussion
From the case study, it is evident that the proposed co-operation model and the ADMM-based decentralized algorithm could facilitate local energy accommodation and increase the total revenue of the local energy system. However, to efficiently solve the network constrained operation problem, several simplifications have been made. In this section, we would make a deeper discussion on the rationality of the assumptions and figure out the effectiveness of the proposed model.

Rationality of Assumptions in the Network Model
As is shown in Section 3, some assumptions are made to reformulate the operation model as a MISOCP form. To verify the rationality of these assumptions made in our proposed operation model, we compare this model with the CFVT model proposed in [20]. The CFVT model is widely used in the industry and can be taken as the benchmark. Figure 12 illustrates the relative gap between the result from the CFVT model and the proposed model at each pipe at each time slot. We can see that the gap is no greater than 2.5% in any situation. The result shows that the difference in DHN dynamics between the two models is very tiny. the industry and can be taken as the benchmark. Figure 12 illustrates the relative gap between the result from the CFVT model and the proposed model at each pipe at each time slot. We can see that the gap is no greater than 2.5% in any situation. The result shows that the difference in DHN dynamics between the two models is very tiny.
In our model, more accurate expressions are proposed for modeling both DHN and PDN, where both of them involve conic constraints. Additionally, the energy conversion characteristic of heat pumps is modeled by quadratic functions and then convexified through the SOC relaxation. One important problem is the exactness of these conic constraints. The exactness of conic constraints in the AC power flow model can be ensured in the PDN [16]. However, in the DHN, the exactness of conic constraints is not guaranteed, which might cause the obtained solution physically infeasible. To address this challenge, penalty terms to guarantee the exactness are added to the objective. Figure 12. Comparison of heat loss along pipes between two models. Table 2 shows the impacts of penalty terms on the exactness of all types of conic constraints. The exactness of our proposed heat pump model is always satisfied no matter there are penalty terms or not. This is because minimizing the power consumption expense has been incorporated in the objective of the proposed model. However, the exactness of pressure loss constraints and heat loss constraints could only be achieved by introducing a small penalty coefficient.

Conclusions
This paper proposes a co-operation model for integrated power and heat systems considering network constraints and the distributed algorithm based on ADMM for the decentralized optimization for PDN and DHN. In the co-operation model, a convex approximation for the PDN as well as DHN is applied. Additionally, a novel quadratic heat pump model and the heat storage model are proposed to illustrate the physical properties more accurately. Last but not least, a decentralized optimization algorithm based on ADMM is proposed. In our model, more accurate expressions are proposed for modeling both DHN and PDN, where both of them involve conic constraints. Additionally, the energy conversion characteristic of heat pumps is modeled by quadratic functions and then convexified through the SOC relaxation. One important problem is the exactness of these conic constraints. The exactness of conic constraints in the AC power flow model can be ensured in the PDN [16]. However, in the DHN, the exactness of conic constraints is not guaranteed, which might cause the obtained solution physically infeasible. To address this challenge, penalty terms to guarantee the exactness are added to the objective. Table 2 shows the impacts of penalty terms on the exactness of all types of conic constraints. The exactness of our proposed heat pump model is always satisfied no matter there are penalty terms or not. This is because minimizing the power consumption expense has been incorporated in the objective of the proposed model. However, the exactness of pressure loss constraints and heat loss constraints could only be achieved by introducing a small penalty coefficient.

Conclusions
This paper proposes a co-operation model for integrated power and heat systems considering network constraints and the distributed algorithm based on ADMM for the decentralized optimization for PDN and DHN. In the co-operation model, a convex approximation for the PDN as well as DHN is applied. Additionally, a novel quadratic heat pump model and the heat storage model are proposed to illustrate the physical properties more accurately. Last but not least, a decentralized optimization algorithm based on ADMM is proposed. From the numerical results based on a PDN33-DHN32 system, the following conclusions can be summarized:

•
The proposed model could generate solutions with high accuracy and optimality. On one hand, precise models for both electricity and thermal networks are adopted. Particularly, the VFVT DHN model provides more flexibility than models with the mass flows fixed without loss of accuracy. On the other hand, improved descriptions on heat pumps with variable COP as well as the stratified water tank for heat storage are included and reformulated to integrate with the network model, making the proposed model closer to the actual situation.

•
Flexibility of the electricity distribution system could be enhanced through the co-operation of integrated energy system. As a result, the co-operation model of an integrated system could facilitate local energy accommodation and reduce the total operation cost. Compared to the DO, CO's total operation cost is lower by 2.3% to 4.4% under different scenarios. Additionally, results show that the grid supply electricity is about 149.1MWh per day in DO, 23% larger than 121.4 MWh in CO.

•
Decentralized optimization based on ADMM could provide a near-optimal solution but with little sharing of private information. With the distributed techniques, the co-operation could be achieved in a decentralized way with little loss of efficiency. Due to the nonconvexity of MISOCP, there would be a rather small gap (≤0.1%) between the decentralized case with the centralized one. Nevertheless, the operation costs in both cases are significantly decreased compared to the decoupled operation of DHN and PDN.
For future work, market mechanisms for the integrated electricity-heat system have drawn great attention. We are interested in exploring efficient market mechanisms and market participants' behaviors in the local energy system considering PDN-DHN constraints.