Optimal placement and sizing of heat pumps and heat only boilers in a coupled electricity and heating networks

Multi-energy systems are reported to have a better environmental and economic performance relative to the conventional, single-carrier, energy systems. Electrification of district heating networks ...


Introduction
Multi-energy system (MES) increases reliability and efficiency by exploiting the synergy of various energy carriers, such as electricity, heat and gas, interacting at a building, district and/or national level [1]. Coupling of electricity and thermal networks, for example, enables to balance the fluctuating renewable electrical sources using low cost thermal storage [2].
Moreover, electrification of district heating network through integration of renewables, heat pumps and combined heat and power (CHP) plants is given more emphasis in order to improve the efficiency of the heating sector in Europe [3,4].
A heat pump (HP) uses electricity to transport heat from lower temperature region to higher temperature region while a CHP produces both heat and electricity from fuel sources, such as waste, biomass and natural gas. Incorporating more of these technologies, in other words, means strengthening the coupling between the electricity and heating networks. Lund et al. [5] also identified that the capability of thermal networks to interact with smart electricity grids is one of the characteristics of the future 4 th generation district heating networks.
Although it is believed that MESs have better environmental and economic performance relative to the "classical" singlecarrier energy systems, it is still challenging to get the details of their operational parameters as the conventional modelling and simulation tools are not designed for such hybrid grids [6]. A heat pump (HP) owned by district heating network operator, for example, is assumed as a fixed demand from the electricity grid point of view. A district heating operator, on the other hand, considers all the other heat sources and decides on the most economical operational strategy of the HP to meet the demand. As the two distribution system operators act independently, any possible difference between the assumed and actual operational strategy of the HP may compromise the efficiency of the overall system [7]. Moreover, the HP acts as a distributed generation for district heating network (DHN) while it is a distributed load from the electricity network point of view. Optimal placement of a distributed generation, usually near to the load center, helps to decrease losses in the distribution network [8]. Thus, the HP will be economical if it is installed near to the electricity generation and near to the thermal load at the same time. In the case of the electricity generation being located far from the thermal load center, both electricity and heating network should be considered to arrive at an optimal location of the HP.
Both size and location of distributed generation (DG) influences the power loss and operational parameters of a multicarrier distribution networks. As analytical methods may not be appropriate to optimally place multiple DGs, Kansal et al. [9] proposed a hybrid of PSO and analytical algorithms. The PSO is used to find the optimal location of DGs which is followed by the analytical method to determine their optimal sizes. Yammani et al. [8], on the other hand, used a variety of genetic algorithm to optimally place and size DGs in electricity network. Pazouki et al. [10] also used genetic algorithm to find out the optimal place and size of a CHP plant in a MES consisting of gas, heat and electricity energy carriers. Heat demands are assumed to be supplied locally either from a CHP or a gas boiler. In all of the papers, minimization of losses in the electricity network is considered as a main objective function even for a DGs that produces both heat and electricity.
Unlike the DGs in the electricity network, the optimal placement of DGs in the heating network is not studied very well.
Marguerite et al. [11] used linear programming to study the optimal dispatch of multi-source DHN while Vesterlund et al. [12] used an evolutionary mixed integer linear programing for thermo-economical optimization of a multi-source DHN. Neither of the optimal placement, nor the coupling with the electricity network are considered in their study. A.
Shabanpour-Haghighi and A. R. Seifi [13] could be the first to propose a teaching-learning based algorithm to solve the optimal power flow of a MES consisting of heat, electricity and gas networks. CHP plants, boilers and compressors were used as a coupling technologies at selected nodes. However, the optimal placement of the technologies is not considered.
The main challenge in studying the optimal placement and sizes of coupling technologies in a MES is the absence of suitable modelling tools that are flexible enough to configure the load flow equations automatically when the location of coupling technologies in a MES is varied. Allegrini et al. [14] reviewed existing software tools that can be used to study urban energy systems. Some of the tools presented, such as EnergyPRO and EnergyPLUS, are used at the prefeasibility stage to do techno-economic analysis, while others, such as TRNSYS are used to do detailed modeling at building/plant level. However, there is no tool suggested which can be used to model the detail network parameters of MES.
On the other hand, different approaches of analyzing MES have been reported in literature. Levihn [15] presented the correlation between electricity price signal and the scheduling of coupling technologies (HPs, CHPs and electric boilers), by analyzing historical data collected from a real DHN. Arnaudo et al. [16] used a merit based techno-economic analysis to compare centralized and distributed HPs in an integrated urban energy system. Ma et al. [7] also used a similar approach to study a MES consisting of electricity and heating energy carriers together with renewables. Geidl and Andersson [17] introduced an energy hub concept to model the conversion, transformation and storage relationship between different energy carriers. Aghamohamadi et al. [18], Carpaneto et al. [19], and Gabrielli et al. [20] followed similar approach to find out the combination of different energy carriers in a single energy hub that gives the least operational cost. Rakipour and Barati [21] also conducted probabilistic optimization in a single energy hub with the consideration of renewable uncertainties. In all of them, the focus was mainly to assess the feasibility of different coupling technologies. The operational parameters of the constituent energy networks and, hence, the associated distribution losses were not considered.
Awad et al [22] and Liu et al [23] took the network parameters into consideration while modeling the load flow problems of heating and electricity networks simultaneously. Liu and Mancarella [24] extended the model presented in Liu et al [23] into an energy system consisting of gas, electricity and district heating network. Shabanpour-Haghighi and Seifi [25] also reported a load flow study for similar MES using empirical formulations rather than the energy hub concept. The node-loop equations used in their hydraulic models, however, needs assumption of pseudo-loop paths. The identification of pseudo loops is not straight forward and makes it difficult to develop a general algorithm for the hydraulic model [26].
They also assumed constant value of the overall heat transfer coefficients for all mass flows in their thermal models which is far from reality. Moreover, only loosely coupled networks are considered in their case studies.
Widl et al. [6] studied technical and economic aspects of coupled electricity and heating networks by running two commercially available tools, DIgSILENT and Dymola, in a co-simulation environment called FUMOLA. A single electric boiler, which runs whenever there is surplus of electricity from the solar PV, is used to loosely couple the two networks. Ayele et al. [27], on the other hand, presented a general and flexible load flow model of a highly coupled MES based on an extended energy hub approach. Its modularity and flexibility is illustrated using electricity and heating networks that are coupled by multiple CHPs.
In summary, only few researchers have made an effort in developing a modular and flexible load flow model of highly coupled MES. Furthermore, optimal placement of coupling technologies, such as HPs is not addressed in literature. The scientific contribution of this paper in filling those gaps are the following: • A highly coupled electricity and heating networks consisting of solar PV, heat-only boiler (HOB), HP and waste to energy CHP plants are considered. The HP-HOB and HP-CHP interactions between electricity, heat and fuel are modelled as modular units of the extended energy hub approach.
• A nested particle swarm optimization (PSO) is proposed to find out the optimal placement and sizes of HPs and a HOB.
• Detailed operational parameters, distribution losses and network constraints of both the heating and electricity networks are considered in the optimization.
The remaining part of this paper is organized as follows. Section 2 presents the methodology (the extended energy hub approach, the HP-HOB and HP-CHP interactions and the nested PSO) in brief. Section 3 describes case studies selected for illustration. Section 4 presents the results and discussion and Section 5 concludes the paper.

The extended energy hub modelling approach
In this paper, the steady state load flow problem of coupled electricity and heating networks are solved using the extended energy hub approach. This approach, the details of which is discussed in [27], is a deterministic steady state model in which the pipe and transmission line physical parameters are assumed to be constant. Generation and load profiles are fed into the model as inputs. At each hub, the local demands, generations and net injections into the multi-carrier network are modeled in modular form using coupling matrices. The power injections for each carrier are, in turn, modeled as a function of each network's physical and operational parameters. Each hub can act as a consumer or a producer depending on its local demand, generation and coupling technology.
A general form of a coupling matrix that describes the relationship between active electricity, reactive electricity, heat and fuel power at each hub in a coupled electricity and heating network is given by equation (1) [27]. Equation (1) is straightforward for modelling solar PV, wind plant, CHP and HOB as they are source of electricity or/and heat networks.
However, HPs are loads for electricity networks and, hence, equation (1) needs modification as it will be discussed in the following subsections.
where represents a coupling coefficient relating generation type with load type ; Pepg, Peqg, Phg, and Pfg are local generations of active electric, reactive electric, heat and fuel powers, respectively while Pep, Peq, and Ph represent the active electric, reactive electric and heat power injections into the network, respectively.

Fuel, heat and electricity interaction in HP-HOB energy hub
The HOB uses fuel to produce heat. Its performance is expressed by its thermal efficiency, . The HP, on the other hand, uses electricity to transfer heat from a colder region to a hotter region. It can also be used to increase the temperature that is initially available so that it can be used for higher temperature application. The source of heat can be sewage, seawater, geothermal, waste heat or ambient air. Its performance is expressed using a coefficient of performance (COP) which is the ratio of heat power transported to the amount of active electric power consumed. Although the COP depends on the type of refrigerant, sink temperature, source temperature and compressor speed control, it can be assumed constant in certain operation conditions. For a 10℃ source temperature, a HP operating in a DHN can have an average COP of 5 [28]. Levihn [15], on the other hand, concluded that a COP of 3.5 for a seawater based large scale HP operating in a temperature range of 65 -115℃ is practically achievable. The effect of COP variation due to partial loading and operating temperatures of a DHN at the distribution level is negligible [29]. Although series connection of HP with HOB or CHP can be used to increase the output temperature level at transmission level (i.e. at sources far from the consumers) [30], parallel connections are preferable to have higher and stable COP at distribution level (i.e. near to the consumers) of the DHN [29]. Apart from the COP, the power factor, ! , is used to calculate the reactive power consumption of the HP. All the HPs considered in this study are assumed to have a COP of 4.0 and a lagging power factor of 0.9.

Fuel, heat and electricity interaction in HP-CHP energy hub
A waste to energy CHP plant is considered in this study. It produces both active electric power and heat from municipal solid waste. The typical municipal solid waste has a lower heating value of 11MJ/kg [31]. Profitability of waste to energy plants highly depends on the gate fee, which varies from country to country between 10$/ton to 100$/ton [32]. A gate fee of 50$/ton of waste is considered in this study. The performance of a waste to energy CHP is expressed by its electrical ( / and thermal ( 0 ) efficiencies. In this study, a waste to energy CHP plant (with an internal combustion engine) that has an electrical efficiency of 31% and thermal efficiency of 57% is assumed [31]. Depending on whether it is running at lagging or leading power factor ( 1 ! ), the CHP plant either consumes or produces reactive power. A lagging power factor of 0.9 is assumed.
The interaction between fuel, electricity and heat in an energy hub that has an HP and a CHP is shown in  (2) and (3), power factor is considered to be positive if it is lagging and negative otherwise.

Equations governing the electricity network
AC power flow model of an electricity network is used to define the per unit active and reactive power injections at any bus k as a function of all nodal voltages and line parameters, as shown in equations (4a) and (4b) respectively [33].
where @ 38 is the voltage angle difference between bus k and bus j; < F8 + GA F8 = H F8 is an element of the network admittance matrix (the details on its formation can be referred from [27]); |6 3 | is voltage magnitude at bus k and N is the total number of buses.

Equations governing the heating network
A thermo-hydraulic model of the heating network is governed by nodal energy balance, mass balance, pressure drop and temperature drop equations. Equation (5) defines the heat power injection into the network from hub k, 3 , as a function of mass flow from the hub, IJ 3 , and its temperature on the supply and return sides, K L3 and K M3 respectively.
where is specific heat capacity of water.
The mass flow from the hub is subjected to continuity of flow at the corresponding node, which is given by equation (6).
In addition, pipe equations are used to relate the pipe flows with the upstream and downstream hydraulic heads as shown in equation (7).
∑ all mass flows into the node 3 = ∑ all mass flows out of the node 3 (6) where IJ 8F is the pipe mass flow from node j (at hydraulic head \ 8 ) to node i (at hydraulic head \ F ) and ] F8 is the corresponding pressure resistance coefficient which (which is further expressed as a function of flow rate and surface roughness [27]).
Thermal energy balance of mixing water at a given node is given by equation (8).
where (IJ 8#ef0 , K 8#ef0 terms are the outgoing mass flow rates and water temperature at a given node j while (IJ 8#Fg , K 8#Fg terms are the incoming mass flow rates and water temperatures at the same node j respectively. The overall heat transfer coefficient between the water in the pipe and the soil is given by equation (9a)  where K r_ gt and K r_L0ƒM0 are outlet and inlet water temperatures, respectively. K e is soil temperature at the surface of the pipe; U is the overall heat transfer coefficient; " % , " | , " … and " † are the inner radius of carrier pipe, outer radius of the carrier pipe, outer radius of an insulation layer and outer radius of an outer jacket, respectively; ℎ r ,ˆ|,ˆ… and ˆ † represent the convective coefficient of water, thermal conductivity of carrier pipe, thermal conductivity of the insulating material and thermal conductivity of the outer jacket. The convective coefficient of water is function of its flow rate and conductivity [34].
The heat transfer coefficient and the temperature difference between the water and the soil, and the mass flow rate are determinant factors for the heat loss. The supply and return temperatures of a DHN in Stockholm, for example, varies from 60 -115℃ and 35 -55℃, respectively [15]. As a result, the relative percentage of heat losses varies throughout a year. The annual average heat loss in Sweden amounts 9.1% of the total heat produced while it reaches up to 21.0% in Lithuania and 27.8% in Korea [35].

Integrated load flow analysis
Once all the equations are formulated in per unit system the overall load flow problem is solved as a single problem using the Newton-Raphson iteration method (the detailed procedures is found in Ayele et al. [27]). The load flow solution enables us to calculate the distribution of losses in both networks, the voltage profile in the electricity network and temperature and mass flow profiles in the DHN.
The electricity consumption by the circulation pumps to overcome the frictional and local pressure losses in the network as well as to guarantee the minimum pressure difference required between the supply and the return sides of a remote consumer substation is given by equation (10). The first term in equation (10) corresponds to the pressure loss in the network while the second term corresponds to the pressure difference required at consumer substations.
where /_ f• is active electric power consumed by circulation pumps on each branch in Watt; " is gravitational acceleration; ∆\ F8 is the frictional head loss on a pipe connecting nodes i and j (assuming straight pipe); :IJ F8 : is magnitude of the mass flow in the pipe and is efficiency of the circulation pump (it is assumed to be 80% in this study). /e‹ is a fraction to take the local pressure losses due to valves and junctions into account while the factor 2 is due to equal pressure drop on both supply and return pipes. ∆\ ‹eg is the hydraulic head difference at the consumer substation in meter while |I ‹eg | is the magnitude of mass flow rate on the primary side of the substation. In this study, the values of /e‹ and ∆\ ‹eg are assumed to be 0.3 and 5.1m (≈ 50ˆ -), respectively [36]. It is also assumed that the total electricity consumption of the circulation pumps is very small when compared to the total electrical demand in the system. Hence, their impact on the operational parameters of the electricity network is neglected.

Optimal placement and sizing of HPs and a HOB
In this study, a novel methodology on optimal placement of HPs and a HOB with the full consideration of both heating and electricity networks is proposed. The objective function is defined by equation (11). The operational cost is the sum of all costs of local generations and imports in both networks including fuels consumed by the coupling devices. It also includes the consumption of circulation pumps. Two scenarios are considered by varying the α value between 0 and 1.
• In Scenario I, the location of HPs is optimized from the point of only heating network without considering the loss in the electricity distribution network (α = 1).
• In Scenario II, both networks are considered in the optimization (α = 0).
IBC™ * š›-oeB>C-n >?oe − • >?oe > šnš=oe›B=Boež ŸB?oe›B ¡oeB>C n>?? ¢ For each candidate optimal solution, the load flow equations (2) A nested PSO is implemented in this study. The outer PSO is dedicated for optimal placement while the inner PSO determines the corresponding economical dispatch of the HPs and other technologies. PSO is first developed by Kennedy and Eberhart [37] in 1995. It tries to find the global best value in analogy to the way a flock of birds scatter and regroup.
The whole group is referred to as a population (swarm) and its individual members are called particles. If there are M variables of optimization, the position of each particle in the swarm is defined by a point as shown in (13a) which is updated at each iteration using equations (13b) and (13c) [37].
Where ¬ F is the current position of particle i; ¯F #g r is the new velocity of particle i, ° is the inertia/damping factor; ¯F e is the current velocity of particle i; = % is a self-accelerating factor; = | is a global accelerating factor; ¬ F# L0 is the best position of the particle in the past; ¬ # L0 is the global best position achieved by the swarm in the past; › % and › | are random numbers between 0 and 1 and ¬ F#g r is the new position of particle i.
The flow chart of the nested PSO is shown in Fig 3, which follows the following steps. iii. The position of each particle in the outer PSO is updated according to equation (13). vi. The position of each particle in the inner PSO are updated according to equation (13).
vii. Integrated load flow is then computed for each particle in the inner PSO loop which is followed by fitness evaluation using equation (11). viii.
The fitness values are compared with the particle's best fitness and the global best fitness values in the inner PSO. Accordingly, the particle's best position and the global best position are updated ix.
Steps (vi)-(viii) are computed for each of the M particles in the inner PSO. x.
If the preset maximum iteration of the inner PSO (MaxIter2) is not reached, step (ix) is repeated.
xi. The global best fitness value returned from the inner PSO loop is used as the particle's fitness value in the outer PSO loop.
xii. This value is used to update the particle's best and the global best positions in the outer PSO. xiii.
Steps (iii)-(xii) are computed for each of the N particles in the outer PSO. xiv.
If the preset maximum iteration of the outer PSO (MaxIter1) is not reached, step (xiii) is repeated.
xv. The best location of the HPs and HOB is contained in the global best position of the outer PSO while their economical dispatch is contained in the corresponding global best position of the inner PSO.
In both of the PSO loops, the global best value is updated after evaluating each particle's fitness and the latest global best value is used to update the position of the next particle.

Set maximum iteration for inner PSO as MaxIter2
and start with Iter2=1, and particle j=1 Fig. 3. Flow chart of the nested PSO applied for optimal placement and sizing of HPs and HOB.

Case studies
Two case studies consisting of heating and electricity networks that are coupled with a waste to energy CHP and HPs In addition solar PV and wind power plants are considered as distributed generations. HOB is also considered as a competitive technology. The optimal placement of the HPs depends the availability of electricity and thermal load. In order to make the optimization problem more challenging, the source of electricity are assumed to be installed at hubs with no thermal demands. The first case study deals with a small part of a distribution network in a densely populated area while the second case study deals with a sparsely populated area connected with a relatively larger network size. The

General assumptions
The following assumptions are made for both case studies.
• As the focus of this study is the distribution losses of existing networks, only operational costs are considered.
• The installation of solar PV panels, wind turbines and waste to energy CHP plants are site specific and it is difficult to relocate them. Their locations are assumed to be fixed at certain hubs.
• Hub1 is considered as a slack hub at which any import/export from/to the neighborhood network takes place.
• The price of electricity and heat from the neighborhood is assumed to be 0.2€/kWh and 0.1 €/kWh, respectively.
• A power purchase agreement of 0.12€/kWh for active electric power is considered. As such practical agreement is not yet fully developed for heating networks, it is neglected.
• A reactive power penalty of 0.02€/kvarh is considered for the imported reactive electric power at the slack hub in excess of 48.5% (corresponding to a 0.9 lagging power factor) of the imported active electric power.
• The price of gas used in the HOB is taken to be 0.113/kWh. The HOB thermal efficiency is 90%.
• Operating cost of the waste to energy CHP plant is assumed to be -0.016€/kWh which corresponds to a gate fee of 50 €/ton of waste with an average heating value 11MJ/kg.
• The thermal and electrical efficiencies of the waste to energy plant are assumed to be 57% and 31%, respectively at a lagging power factor of 0.9.
• The operating cost of solar PV is assumed to be negligible.
• All HPs are assumed to run at a COP of 4.0 and a lagging power factor of 0.9.
• The electricity network is assumed to be balanced three phase system. All transmission lines are assumed to be three phase ACSR Waxwing type with resistance, reactance, susceptance and ampacity of 0.262Ω/km, 0.386Ω/km, 4.31µS/km and 480A, respectively [38]. Double conductors are used between hubs 1 and 2.
• All pipes in the DHN are of DN-100 type with standard insulation layer. The parameters (carrier pipe thickness and outer diameter of 3.6mm and114.3mm, and outer jacket thickness and outer diameter of 3.2mm and of 200mm, respectively) are taken from isoplus ® [39]. Maximum allowed flow rate in the pipe is 9.3kg/s. Thermal conductivity of the carrier pipe, outer jacket and insulation material are 40W/mK, 0.4W/mK and 0.027W/mK, respectively. Surface roughness of the carrier pipe is assumed to be 0.05mm.
• All hubs acting as a heat source are assumed to have a supply temperature of 70℃ and all hubs acting as a heat consumer are assumed to have a return temperature of 50℃ on the network side. A uniform soil temperature of 4.36℃ is assumed around the supply and return pipes of the DHN.
• Both residential and commercial demands are considered. Peak demand of one/two residential buildings and a small hotel located in Idaho, USA are taken as a reference from an Open Energy Information database [40].
An electrical power factor of 0.95 and 0.9 lagging is assumed for residential and commercial demands, respectively.

Specific inputs for Case 1
A small distribution network topology for Case 1 is shown in Fig 4. The load distributions are summarized in Table 1.
The CHP with a maximum of 200kW waste intake capacity is installed at Hub3. A solar PV with 50kW output is connected at Hubs 5. The electricity distribution network is at 0.4kV. Two heat pumps, each of which are rated at 25kWe, and a HOB that can take up to 30kW of gas are considered in the optimization. In order to find their optimal placement in the network and their economical dispatch, a swarm of 5 particles with 15 iterations is considered in the outer PSO, while a swarm with 10 particles and 70 iterations is used in the inner loop.

Specific inputs for Case 2
A relatively larger and sparsely populated distribution network topology for Case 2 is shown in Fig 5. The load distributions are summarized in Table 2. The CHP is assumed to be installed at Hub2 and its maximum waste intake

Optimization results for Case 1
The optimal locations of the two HPs and the HOB together with their operating size and the resulting distribution of losses are summarized in Table 3. Imported powers for each energy carrier are also given in the same table. Although the highest heat demand is at Hub4, the optimal location of the HPs in Scenario I are found to be at Hubs 2 and 6. This is due to the existence of a cheap heat source from the CHP plant at Hub3 which has no other connection other than Hub4.
In scenario II, on the other hand, the HPs are placed at Hub2 instead of Hub3 or Hub5 where there is cheap electricity generation. The reason is the existence of heat demand at Hub2 and due to the fact that larger proportion of the active electricity comes from the slack hub (see Table 3). In both scenarios, the CHP and solar PV are selected at their maximum output followed by HPs. No HOB is selected. Importing heat is also found to be expensive relative to the heat from HPs.
The results for both scenarios shows that the optimal placement of HPs is not aligned with active electricity generation hubs (hubs 1, 3 and 5). They are, rather, placed at some of the hubs where there is heat demand. In scenario I, 54.6% of the total heat generated (290.5kW) is transported in the heating network while 100% of the total active electricity generated (462.38kW) is transported in the electricity network. In the case of scenario II, however, 73.55% of the total heat generated (291.3kW) and 100% of total active electricity generation (459.52kW) is transported in the corresponding networks. As a result, scenario I has lower heat loss (1.6% of the heat demand) than Scenario II (1.9% of the heat demand).
The active electric power loss, on the other hand, is higher (3.9% of the demand including HPs' consumption) for Scenario I than Scenario II (3.2% of the demand including HPs' consumption).
The pumping power required for circulation in scenarios I and II are 0.12kW and 0.16kW respectively, which are negligible relative to the total electricity demand. It is also found that only 0.05% and 0.07% of this pumping energy is used to overcome the friction and local pressure losses in the network respectively while the remaining is used to keep the pressure drop at the consumer substations. The cost of pumping energy for scenarios I and II is 0.024€/h and 0.032€/h which are less significant compared with the corresponding cost of heat loss (0.46€/h and 0.54€/h respectively). This is in agreement with the findings of [12] but in disagreement with the findings of [36]. The possible reasons could be the difference in the pipe parameters (diameter, length etc.) and operational strategies. The operating cost of the whole network is 68.3€/h in scenario I. This value is reduced to 67.68€/h in scenario II which shows nearly 1% saving. This is due to the fact that the price of electricity is twice of the heat if they are compared at the importing level. The difference grows to even four times if the heat from the HPs is considered.  Figure 6 depicts the temperature and mass flow profiles at different hubs. A positive mass flow implies that the hub is acting as a source of heat. Hub6, for example, is acting as a source in Scenario I and as a consumer in Scenario II. When the mass flow from the hub is closer to zero, the water gets more time to lose its heat and, as a results, its temperature gets closer to the soil temperature (for example the return temperature at Hub2 in Scenario I). If the mass flow is zero, as in the case for Hub5 in both scenarios, the steady state temperature of water on both the return and supply pipes becomes equal to the soil temperature.
(a) (b) Fig. 6. a) The mass flow injected into the network from each hub and b) the corresponding supply and return temperatures for the two scenarios of Case 1.  The voltage magnitude profiles, as shown in Fig 8(a), are within the limits of equation (12a). Scenario II has a relatively better voltage profile. The root mean square current through each transmission line are illustrated in Figs 8(b). It looks that the current in Line1-2 is above the acceptable limit of 480A. But, since double conductors are assumed for each phase, the current in each conductor is still in the acceptable range. It also shows that the current in Line2-4 in Scenario I is 484.5A which slightly higher than the acceptable limit. It means that optimizing HPs without the consideration of electricity network constraints (as is the case for Scenario I) could lead to overloading of the electrical distribution network. Table 4 presents the optimal location and operational sizes of the HPs, HOB together with the economical dispatches of CHP and wind plants. It also shows the distribution losses and amount of import for each energy carrier. The full capacity of the CHP and the wind power plants are selected in both scenarios. The remaining electricity is imported from the neighborhood grid while the remaining heat is generated using HPs. Similar to Case 1, HOB and the neighborhood DHN are avoided. Similar to Case 1, hubs with electricity generation (hubs 1, 2 and 26) are not selected for optimal placement of the HPs. This is again due to the absence of thermal loads at those hubs. However, as it can be referred from Table 4 and Fig 5, the HPs in scenario II are closer to the hubs with electricity generation when compared to the HPs in scenario I. In scenario I, 53.03% of the total heat generated (3.22MW) is transported in the heating network while 100% of the total active electricity generated (4.06MW) is transported in the electricity network. In the case of scenario II, however, 93.55% of the total heat generated (3.23MW) and 100% of total active electricity generation (3.98MW) is transported in the corresponding networks.

Optimization results for Case 2
The heat loss in the distribution network for Scenarios I and II are 9.35% and 9.69% of the heat demand, while the active electricity losses are 5.39% and 3.17% of the electricity consumption, respectively. Compared with Case 1, the losses are higher which is partly due to longer branches and larger percentage of power flows in the network. It is important to note that nodes 9 and 10 in Fig 5 are equally optimal as node 8 due to the symmetry in the network dimension and demands at the nodes (see Table 2). This applies to other symmetrical nodes too. The more asymmetry in the network, the faster the PSO gets to the global optimal solution.  Figure 9 shows the nodal mass injection from each hub and the corresponding temperature profiles on the supply and return pipes. The nodal mass flows for both scenarios lookalike for all nodes except where the HPs are installed (Fig   9(a)). Hubs with no heat consumption and production, for example 13, 19 and 26, can be treated as either sources or consumer in the numerical computation depending on the direction of mass flow obtained from the load flow solution.
The mass flow always is zero at these hubs. However, the temperature can vary depending on the assumption. In any case, the heat power injected into the network is zero and the mass and temperature balances at the mixing nodes are always guaranteed. Hub13, for instance, is considered as a source in Scenario I in which its supply temperature is fixed at 70℃. Its return temperature is equal to the outgoing return temperature of the mixing at Node13. On the other hand, Hub13 is treated as a consumer in Scenario II in which case the return temperature is assumed to be equal to the soil temperature as there is no flow. The supply temperature will then be equal to the outgoing supply temperature of the mixing at Node13. Hub19 is treated as a source in both scenarios.
(a) (b) Fig. 9. a) The mass flow injected into the network from each hub and b) the corresponding supply and return temperatures for the two scenarios of Case 2.  In this paper, a nested PSO algorithm, integrated with the load flow model based on an extended energy hub, is proposed to find out the optimal location and sizes of distributed generation and coupling technologies. Two case studies of coupled heating and electricity networks consisting of a CHP, HOB, HPs, Solar PV and wind power plant are considered. The HPs are found to be more economical than the neighborhood DHN and a HOB for the given price data. The results show that optimizing the location and sizes of HPs only from the heating network point of view could lower the heat loss, but sometimes it could overload and increase the loss in the electricity distribution network, especially when HPs are installed at very low voltage distribution networks. It could also lead to unacceptable voltage profiles. On the other hand, an integrated optimization which considers the constraints and distribution losses in both networks results in a lower operational cost for the overall system with the possibility of better voltage profile in the electricity network. It is also observed that the electricity consumption of circulation pumps is not significant compared to the total electricity demand.
However, the cost might become significant compared to the cost of heat lost specially when there are cheap heat sources in the network. It is observed that the frictional pressure losses are equally important as the pressure losses at the consumer substations for larger networks. Hence, it might be worthy to consider the total cost of pumping energy in thermoeconomic optimizations of large networks.
The methodology presented in this paper is a step forward in optimization of integrated energy systems. It can be applied to different MES consisting of different distributed generations and coupling technologies. It can be further complemented by considering time series data for solar PV and wind power generation together with thermal storage units. Once the optimal locations of coupling technologies are identified, additional optimization variables such as temperature and transformer tap setting can be incorporated. The methodology can also be applied for future energy systems at planning phase in which case the investment cost of the network and coupling technologies will be taken into account. Other optimization algorithms can also be applied and compared with.

Acknowledgments
The research presented is performed within the framework of the