Modeling of a District Heating System and Optimal Heat-Power Flow

With ever-growing interconnections of various kinds of energy sources, the coupling between a power distribution system (PDS) and a district heating system (DHS) has been progressively intensified. Thus, it is becoming more and more important to take the PDS and the DHS as a whole in energy flow analysis. Given this background, a steady state model of DHS is first presented with hydraulic and thermal sub-models included. Structurally, the presented DHS model is composed of three major parts, i.e., the straight pipe, four kinds of local pipes, and the radiator. The impacts of pipeline parameters and the environment temperature on heat losses and pressure losses are then examined. The term “heat-power flow” is next defined, and the optimal heat-power flow (OHPF) model formulated as a quadratic planning problem, in which the objective is to minimize energy losses, including the heat losses and active power losses, and both the operational constraints of PDS and DHS are respected. The developed OHPF model is solved by the well-established IPOPT (Interior Point OPTimizer) commercial solver, which is based on the YALMIP/MATLAB toolbox. Finally, two sample systems are served for demonstrating the characteristics of the proposed models.


Introduction
With the gradual depletion of fossil energy and the degradation of the environment, renewable energy generation resources, mainly solar power and wind power, have received extensive attention [1][2][3].On the other hand, combined heat and power (CHP) plants and electric boilers (EBs) are widely employed [4], which intensifies the coupling between a power distribution system (PDS) and a district heating system (DHS) concerned.Thus, it is becoming necessary to take the PDS and DHS as a whole in both planning [5][6][7] and operation [8,9] aspects.To this end, it is very demanding to develop mathematical models and analytical methods for energy flows analysis.
The integration of PDS and DHS is discussed in some recent publications.The coupling devices, i.e., the CHP plant and EBs, and heat balance constraints are included in the presented rolling scheduling strategy for cogeneration units in [10].In [11], CHP plants, heat storage tanks, and EBs are co-optimized to enhance the accommodating capability of wind power generation in the PDS dispatch.A comprehensive model for small-scale integrated heating and cooling energy systems is developed in [12] in order to minimize the system operation cost.However, due to the absence of energy flows analysis, the heat losses and pressure losses in [10][11][12] are assumed as constants or are just overlooked.However, these losses are considerable in an actual DHS and have significant impacts on the distribution of heat flow [13,14].
To address the above mentioned problems, the well-established Newton-Raphson method has been used in some publications to analyze the heat flow and power flow simultaneously.A Newton-Raphson based combined analysis method is developed to investigate the performance of integrated power and heating systems (IPHS) in [15].The Newton-Raphson method is employed to simulate the dynamic process of pipe water flow rates, water temperatures, and room air temperatures in thermal and hydronic systems [16].The DHS model is described by hydraulic and thermal sub-models, and the Newton-Raphson method used to solve in [17].The models and that are methods presented in [15][16][17] are feasible for analyzing a DHS in both hydraulic and thermal fields, but the computational burden is heavy.Moreover, the impacts of pipeline parameters on the losses are not addressed in [15][16][17], but represent a very important issue in the planning and operation of DHS.Pipeline parameters are the most important optimization variables in DHS planning studies [18][19][20][21], and the existing Newton-Raphson based methods do not address the impacts of pipeline parameters on expansion and operation costs.
In general, a heating pipe (HP) in a DHS consists of a straight pipe (SP), a local pipe (LP) [22], and a radiator [23].It is pointed out in [24] that the LP and radiator play important roles in an actual DHS, and have impacts on the performance of IPHS as well as the SP.However, to the best of our knowledge, the LP and radiator have not yet been addressed in the coordinated planning and operation of IPHS.This is mainly due to the very complicated procedure of solving transient state models.Especially, if the coordinated planning problem is formulated as a mixed-integer nonlinear programming problem, then the computational complexity will be very high.To solve the problem efficient, a steady state model with more components of the HP and access to the existing PDS models is demanding.
Given the above background, an equivalent steady state model is presented in this work with more components of the HP as well as the combination of the hydraulic sub-model and the thermal one.In addition, accurate formulas of heat losses and pressure losses are attained, with the influences of pipeline parameters and the environment temperature being considered.Then, the optimal heat-power flow (OHPF) model is developed to analyze the energy flows of IPHS.Two sample systems are next employed to demonstrate the features of the proposed models.It should be mentioned that the OHPF model could be used in more complicated systems, such as integrated heat-gas-power energy systems [25], by including the corresponding operational constraints [26].
The rest of this paper is organized as follows: the steady state model of DHS is presented and simplified in Section 2. An OHPF model that is based on the steady state model is developed in Section 3. Case studies and simulation results are presented in Section 4. Finally, conclusions are given in Section 5.

The Steady State Model of DHS
The steady state model to be developed for DHS consists of three main parts: the SP model, the LP model, and the radiator model.

The SP Model
In [27], a distributed parameter model (DPM) for the SP is presented, with the heat energy of mass flow and the pressure of pipe node being regarded as the state variables of DHS.The DPM can be depicted as Figure 1 and is formulated as Equations ( 1) and (2).
where Z L and Y L are the equivalent pressure resistance and thermal conductance of SP, respectively; Φ 1 and p 1 are the inlet heat flow and pressure of SP; Φ 2 and p 2 are the outlet heat flow and pressure of SP; D, d L , λ and L denote the mass flow, inner diameter, friction factor, and length of SP, respectively; ρ and c denote the density and specific heat of hot water, respectively; k is the heat transferring coefficient of the pipe wall; T 0 is the environment temperature; and, z 0 and y 0 are the unit pressure resistance and unit thermal conductance of SP, respectively.In this work, the DPM of SP is employed as an important part of the HP in a DHS.It is pointed out in [24] that the LP and radiator play important roles in an actual DHS, and have impacts on the performances of the IPHS and SP.Thus, the modeling of the LP and radiator will be carried out next to enhance DPM with reference to the existing differential analysis methods for DHS [14,22,23].

The LP Model and Radiator Model
In the DHS modeling and controlling presented in [24], the detailed models of LPs and radiator are addressed.Regretfully, these models are mainly transient state models and are hard to compute.In order to attain a simple steady state model, the typical dual-port network and three-port network are employed, as shown in Figure 2. In this work, four kinds of LP are included, i.e., the tee pipe, reducer union, elbow pipe, and valve.Based on a typical DHS, the formulas of the involved pressure resistances and thermal conductance in Figure 2 can be obtained.
The detailed derivation procedures of some complicated equations are moved to Appendix A, so as to significantly simplify the presentation in the main body of the paper.

The Tee Pipe Model
In Figure 2, a typical three-port network is employed to describe the tee pipe.In [24], a detailed lumped parameter model for the tee pipe is presented, and the pressure resistances of the tee pipe can be attained by comparing the models in Figure 2 and [24].
where Z T1 and Z T2 are the pressure resistances of the two outlet pipes in the tee pipe, respectively; K f is the local loss coefficient; d T and D 1 are the inner diameter and inlet mass flow of the tee pipe, respectively; and, Φ 2,1 and Φ 2,2 respectively denote the heat flows at the two outlet pipes in the tee pipe.

The Reducer Union Model and Elbow Pipe Model
In Figure 2, a typical dual-port network is employed to depict the reducer union and the elbow pipe, respectively.By comparing the reducer union model and the elbow pipe model presented in [24], the pressure resistances of the reducer union Z R and elbow pipe Z E can be attained as ) where D 2,1 and D 7 are, respectively, the inlet mass flows of the reducer union and elbow pipe; 2 denotes the local loss coefficient of the reducer union; d R1 and d E are the inlet inner diameters of the reducer union and elbow pipe, respectively; and, Φ 3 and Φ 8 are, respectively, the outlet heat flows of the reducer union and elbow pipe.

The Valve Model
In Figure 2, a variable resistor model is used to represent the valve.With reference to Equations ( 3)-( 5) and some publications [28,29], the variable local loss coefficient K f is presented here to formulate the variable pressure resistance of the valve Z V , as shown in Equations ( 6) and (7).
where D 5 and K m ∈ [0, 1] are, respectively, the mass flow and the opening coefficient of the valve [29]; d V is the inlet inner diameter of the valve; and, α and β are both fitting coefficients.Φ 6 is the outlet heat flow of the valve.

The Radiator Model
Due to the time coupling between temperature and heat transferring, the present radiator model is a differential equation [24,30].This makes thermal analysis in a DHS very complicated.In order to rebuild a simpler model for the radiator, a two-terminal network is employed, as shown in Figure 2.
where Φ 6 , Φ 7 , Φ L , and Y M denote the inlet heat, outlet heat, heat load, and thermal conductance of the radiator, respectively; p 7 is the pressure of the radiator.
In Equation ( 8), the thermal conductance of the radiator is expressed by the discrete-time decoupled heat load Φ L , and thus makes the heat flow analysis much easier than the existing differential model for the radiator.
Based on the above stated three models, an equivalent simplified τ-type model can be attained, as detailed below.

Simplifications of the DHS Model
Despite the formulas that have been derived, the steady state model in Figure 2 is still computation-intensive, especially for the heat flow analysis.Therefore, some simplification is demanded to make the proposed model more tractable.Without loss of generality, a representative HP is employed to demonstrate the simplification process.
First, a representative HP is transformed into an equivalent τ-type steady state model, as shown in Figure 3.It is assumed that the composition of the representative HP is uncertain, because the composition is variable with the structure of HP.Then, the universal formulas of pressure conductance Y Pij and thermal conductance Y Hij can be attained as where ε Rij , ε Vij , ε Eij , and ε Mij , respectively, denote the binary variables (1 or 0) that are indicating whether or not HP contains the reducer union, valve, elbow pipe, and radiator.For any specific HP, the binary variables and other parameters in the above equations will become constants as long as its composition is specified.Therefore, the constants C 1ij and C 2ij could be used to simplify Equations ( 9) and (10), respectively.In next, Z Lij and Y Lij in Equations ( 9) and (10) will be further expanded and simplified.More specifically, "sinh(•)" and "cosh(•)" in Equation (1) are expanded in Taylor's series with the high order terms omitted.Thus, Z Lij and Y Lij can be formulated as where C 3ij denotes the constant of the representative HP.
Then, Y Pij and Y Hij can be formulated as where C ij denotes the constant of the representative HP.
Up to now, the simplified τ-type model of the representative HP is attained, as shown in Figure 3.The heat losses and pressure losses of HP can be formulated as Equations ( 15) and ( 16), respectively: ∆Φ ij and ∆p ij can be calculated by the constants of the representative HP, i.e., C 2ij , C ij , Φ L ij , c, and k, and the environment temperature T 0 .In other words, different pipeline parameters and temperatures will lead to different optimization strategies for a DHS regarding the coordinated planning and operation.Note that the heat flow and pressure distribution of a DHS decline during the process of hot water transferring, which, respectively, due to the heat transfer with the environment and the frictions of HPs.For the convenience of presentation, ∆Φ ij > 0 is used to describe the heat transferring process from hot water to environment, and ∆p ij > 0 represents pressure reduction.Note that the temperature of hot water is limited to be higher than T 0 and the hot water is transferred without additional pressure, then ∆Φ ij and ∆p ij are both assumed to be non-negative numbers in this work.

The Optimal Heat-Power Flow Model
Based on the simplified τ-type model, a unified analysis method of heat and power flow in the IPHS will be addressed in this section.

The Connection Matrix of HP
Traditionally, the DHS adopts a loop structure and consists of two symmetric networks: the supply pipe network and the return pipe one [12].Note that hot water is pressured by the circulating pump (CP) to the loop in a certain direction.To supply the heat demand, the hot water that is drained out of the CHP plant is transferred through the supply pipes network.After the heating process, hot water returns back to the CHP plant through the return pipes network.Hence, the mass flow of the HP can be assumed to be a non-negative number.The symmetric network means that more than one pipe connecting to the end of any HP in a DHS.Thus, the connection matrix of the HP becomes more crucial to simulate the distribution of heat flow and formulate the heat flow balance.In order to attain the connection matrix, an illustration, as shown in Figure 4, is employed: With respect to Figure 4, the connection matrix of the HP can be formulated as where ξ i,l donates an element of the HP connection matrix; Ω H is the set of nodes of the DHS; and, Φ i , Φ j and Φ k , respectively, donate the inlet heat flow of HP i, HP j, and HP k.
It should be mentioned that the heat flow balance is established in the view of the HP rather than based on a DHS node.While hot water supplies heat demand by transferring through the radiator that is connected with a HP, this is significantly different from the situation of the power load that is connected at a specified bus point.To balance the heat flow, the HP is regarded as a generalized node, as shown in Figure 4.As to the pressure distribution, the analysis method is similar to the one that is used in voltage calculation for many years, i.e., the Distflow [31,32].Referring to [31,32], it is known that the incidence matrix of the DHS is very vital to the analysis of the pressure distribution.The incidence matrix can be attained using the method that is presented in [33], and the details of the method will not be presented here due to the limited space.

The OHPF Model
Regarding the HP connection matrix, incidence matrixes, and generalized node, the OHPF model can be formulated mathematically.To coordinate different energy flows, the term "heat-power flow" is defined and employed in the OHPF model with the heat flow and power flow being uniformly modelled.
The presented OHPF model is a quadratic constrained optimization problem, as detailed below.

The Objective Function
The objective function of the OHPF model is defined as

The Constraints of DHS
The operation constraints of a DHS can be formulated as

The Constraints of PDS
The operation constraints of a PDS can be formulated as

The Constraints of Coupling Devices
The operation constraints of three kinds of coupling devices can be formulated as where ∆ OHPF means the optimization objective of the OHPF model; Γ is the run time of an IHPS; Ω L and Ω P are, respectively, the set of HPs in the DHS and the set of feeders in the PDS; ∆Φ l and ∆p l , respectively, denote the heat losses and the pressure losses of HP l; ∆P i , and ∆Q i , respectively, denote the active and reactive power losses of feeder i; ϕ H and ϕ P , respectively, denote the unit price of heat power and that of electrical power; Φ l is the inlet heat flow of HP l; p i is the pressure at node i; Φ L l and D l , respectively, denote the heat load and mass flow of HP l; C l and C 2l , respectively, denote the constants of HP l; Φ CHP l , Φ EB l and p CP l respectively denote the thermal power output of the CHP plant, the thermal power output of the EB, and the pressure compensation of the CP in HP l; A H i,l and A P i,j are, respectively, the elements of the incidence matrixes of the DHS and PDS; Φ l and Φ l are, respectively, the upper and lower heat flow limits of HP l; p i and p i denote the upper and lower limits of the pressure at node i; P i and Q i are, respectively, the active power and reactive power in feeder i; P S i and Q S i, respectively denote the injection active power and reactive power at the power-outflow terminal of feeder i; P L i and Q L i, are, respectively, the active and reactive loads at the power-outflow terminal of feeder i; P CHP i , P CP i and P EB i, respectively denote the power output of the CHP plant, the power demand of the CP, and the power demand of the EB at the power-outflow terminal of feeder i; C CHP is the heat-power ratio of the CHP plants; µ CP is the efficiency factor of the CPs; C EB is the heat-power ratio of the EBs; P i and Q i , respectively, denote the upper limits of active and reactive power in feeder i; R i , X i and ∆V i , respectively, denote the resistance, reactance, and voltage losses of feeder i; V B is the reference voltage of the PDS; V i is the voltage at the power-outflow terminal of feeder i; P S i and P S i respectively denote the upper and lower limits of injection active power at the power-outflow terminal of feeder i; Q S i and Q S i denote the upper and lower limits of injection reactive power at the power-outflow terminal of feeder i; Φ CHP l and Φ CHP l respectively denote the upper and lower limits of the thermal power output of the CHP plant in HP l; p CP l and p CP l respectively denote the upper and lower limits of the pressure compensation of the CP in HP l; P EB i and P EB i denote the upper and lower limits of the power demand of the EB at the power-outflow terminal of feeder i; Λ EB , Λ CHP and Λ CP respectively denote the location transformation matrixes of the EBs, the CHP plants, and the CPs.
In the OHPF model, Equation ( 18) stands for the optimization objective function with heat losses and active power losses included; Equation (19) represents the heat losses and active power losses; Equation (20) represents active power balance and reactive power balance; Equation (21) represents pressure losses; Equation (22) denotes the mass flow balance; Equations ( 23) and ( 24), respectively, denote the upper and lower output limits of D l , Φ l and p i ; Equation (25) represents the active power balance and reactive power balance; Equation ( 26) represents the constraints of P i , Q i , and V i ; Equation (27) represents voltage losses; Equation (28) represents the constraints of injection active and reactive power; Equation ( 29) represents the models of CHP plants [34], CPs [27], and EBs [35], sequentially; Equation (30) includes the upper and lower output limits of CHP plants, CPs and EBs.
The well-established IPOPT (Interior Point OPTimizer) commercial solver based on YALMIP/MATLAB toolbox is employed to solve the developed quadratic OHPF model.

Case Studies
Two sample IPHS systems are served for demonstrating the feasibility of the proposed steady state model as well as the effects of the OHPF model in coordinated operation.

Sample System 1
The network of the first sample system (sample system 1) is shown in Figure 5, in which the CHP plants and CPs are employed as coupling devices between PDS and DHS.

Accuracy of the OHPF Model
The network structure and parameters of sample system 1 are all taken from [24].Specified parameters are listed in Table 1.In order to compare the results between heat-power flow and optimal heat-power flow, the upper and lower output limits of the coupling devices are first assumed to be the same and take the values specified in [24].Since EBs are not addressed in [24], it is assumed that "P EB I = P EB I = 0" and "Φ EB l = Φ EB l = 0".This means that Φ CHP l and p CP l are constants and that they cannot be optimized.Thus, the heat flow and pressure distribution results that were attained by the OHPF model could be used to compare with those that are presented in [24], as shown in Figure 6.From Figure 6, it can be found that the optimization results that were attained by the OHPF model are almost the same as those that are presented in [24].This means that the accuracy of the OHPF model is acceptable, although some simplifications have been made to make the steady state model more tractable, as detailed in Section 2.

Temperature Sensitivity Analysis based on the OHPF Model
To investigate the effect of T 0 on OHPF, three scenarios with different values of T 0 , i.e., −15 • C, −5 • C, and 5 • C, are examined.Equation ( 18) is formulated to minimize the heat losses and active power losses, with T 0 remaining unchanged as long as the OHPF model can be successfully solved.That is because the outputs of the coupling devices could be adjusted in feasible regions to obtain the same minimum ∆Φ and ∆P, and hence the impacts of T 0 can be waived.In order to analyze temperature sensitivity, it is specified that "ϕ H = ϕ P = 0" so as to prevent ∆Φ and ∆P from decreasing during the process of seeking the optimal solution of Equation (18).The optimal values of ∆Φ, ∆p, D, ∆P, and voltage under different scenarios are attained, as shown in Table 2 and Figure 7.In Table 2, ∆Φ decreases significantly with the increase of T 0 , and ∆p maintains in step with ∆Φ.This is due to more drastic thermal interactions happening at low T 0 , which increases ∆Φ when hot water is transferred through a HP.Different levels of thermal interactions have impacts on ∆p by decreasing the mass flow in the corresponding HP, but are not as significant as those on ∆Φ.When a lower T 0 comes along with a higher output of the CHP plant, both in heat and power, different levels of ∆P and voltage of PDS will occur, as shown in Figure 7.In summary, the state variables of the OHPF model, i.e., ∆Φ, ∆p, D, ∆P and voltage, are all sensitive to T 0 .

Sample System 2
Based on the mentioned methods, a software package named "Multi-Flow" is developed to online attain the coordinated optimal operation strategy for any IPHS, as shown in Figure 8.In "Multi-Flow", the OHPF model is coded in MATLAB, and its reliability and efficiency are tested by numerous sample systems.Once an IPHS is built up by component modules, optimal heat-power flow can then be carried out.The involved parameters are included in the menu-bar for setting different simulation scenarios.MATLAB and YALMIP toolboxes are employed to provide computing services for "Multi-Flow".
In sample system 2, a larger IPHS is built based on the software "Multi-Flow" that was developed by us.Three DHSs are integrated into a PDS with three photovoltaic (PV) generating units being included.The impacts of electric heating on the accommodation capability for PV generation are examined for three different load-level scenarios.The daily power and heat loads of sample system 2 are given in Figures A1 and A2 of Appendix B.
To demonstrate the impacts of different load levels on OHPF, Φ and voltage at different time points, i.e., 5:00, 13:00, and 21:00, are calculated and are shown in Figure 9.The values of ϕ H and ϕ P are specified to be 0, as what has been done in sample system 1.Simulation results show that Φ of DHS and the voltage of PDS are both different under various load levels.Specifically, the higher levels of heat and power loads at 13:00 lead to lower heat flow and smaller voltage fluctuation.This is because a higher heat demand will increase the heat outputs of CHP plants.As a result, the power outputs from CHP plants increase significantly to support the power demand and improve the voltage distribution.As shown in Figure 8, EBs can be employed to enhance the accommodating capability for PV power generation.The hourly PV power outputs in a day are shown in Figure A3 of Appendix B. To attain the optimal operation strategy for sample system 2, Equation ( 18) is employed to minimize the operational costs.Moreover, in order to maintain the optimization objective of Equation (18) as unchanged, electricity purchasing costs are not considered, and the punishment costs for PV generation curtailment are not included as well.Then, the following three scenarios with different ratios of the electric load over the heating load are examined to find out the impacts of EBs on the optimal operation strategy: (1) S1: No electric heating load; (2) S2: 50% of the heating load is supplied by EBs; and, (3) S3: The heating load is fully supplied by EBs.
The power outputs of PV units and CHP plants under these three scenarios are attained by the presented OHPF model, as shown in Figure 10.From Figure 10, it is obvious that the power outputs from the PV units increase with the raise of the electric-heating ratio.The power outputs of the PV units under these three scenarios are shown in Figure A3 of Appendix B. The reduced demand of hot water led to a declined thermal power output of the CHP plant.As a result, the heating flow through the DHS reduced, and the heating losses also reduced accordingly.In a nutshell, electric heating makes contributions to the economic operation of a DHS and the accommodation capability enhancement of PV power generation.
Regarding a PDS, heating with EBs could increase the power load demand, and will have impacts on the PDS planning process.In analyzing the economics of EBs, the costs for expansion planning and operation should be taken into account.In-depth study will be carried out in the future to figure out the economics of employing EBs and the coordinated planning for an IHPS.

Conclusions
A steady state model for an IPHS is presented in this work, with the impacts of pipeline parameters and environment temperature being taken into account.In the presented model, various components of HP are modeled, and the hydraulic sub-model is combined with the thermal one.Then, an OHPF model is proposed to analyze the performance of an IPHS.Two sample systems are served for demonstrating the proposed method.It is shown by the simulation results that the accuracy of the presented OHPF model is acceptable and the optimization results of the OHPF model are sensitive to environment temperature.The optimal operation strategies that are attained by the OHPF model are demonstrated with three scenarios under different load levels.The electric heating load makes contributions to the economic operation of a DHS and the accommodation of PV power generation.
The following issues will be carried out in the future: (1) The OHPF model will be employed to optimize the coordinated planning and the operation strategies for integrated heat-gas-power energy systems by including the operation constraints of the nature gas distribution system.(2) Modeling of heat storage equipment will be addressed to expand the steady state model that is presented in this work, and to include more feasible operation modes of IPHSs.(3) Intelligent traffic networks with various kinds of electric vehicles will be included in the expansion planning of multi-energy systems and charging infrastructures.
Author Contributions: Wentao Yang proposed the methodological framework and mathematical model, performed simulations and drafted the manuscript; Fushuan Wen organized the research team, reviewed and improved the methodological framework and implementation algorithm; Ke Wang and Yuchun Huang reviewed the manuscript and provided suggestions; Md.Abdus Salam reviewed and polished the manuscript.All authors discussed the simulation results and agreed for submission.

Conflicts of Interest:
The authors declare no conflict of interest.

Figure 2 .
Figure 2. The steady state model of a typical district heating system (DHS).

Figure 3 .
Figure 3.The simplifications of the steady state model of a representative heating pipe (HP).

Figure 4 .
Figure 4. Attainment of the HP connection matrix.

Figure 5 .
Figure 5.The network of sample system 1.

Figure 6 .
Figure 6.The heat flow and pressure distribution of DHS in sample system 1.

Figure 7 .
Figure 7. ∆P and voltage of PDS in sample system 1.

Figure 8 .
Figure 8.The network of sample system 2 attained by the software "Multi-Flow".

Figure 9 .
Figure 9. Heat flows and pressure distributions at different time points.

Figure 10 .
Figure 10.The power outputs of photovoltaic (PV) units and combined heat and power (CHP) plants under three scenarios.

Table 1 .
Specified parameters of sample system 1.

Table 2 .
∆Φ, ∆p, and D of DHS under different scenarios in sample system 1.
Rij Binary variables that indicating whether or not HP contains the reducer union ε Vij Binary variables that indicating whether or not HP contains the valve ε Eij Binary variables that indicating whether or not HP contains the elbow pipe ε Mij Binary variables that indicating whether or not HP contains the radiator C 1ij , C 2ij , C 3ij , C ij Constants of the representative HP ξ i,l An element of the HP connection matrix Ω HSet of nodes of the DHS