Modelling and Optimising the Marginal Expansion of an Existing District Heating Network

Although district heating networks have a key role to play in tackling greenhouse gas emissions associated with urban energy systems, little work has been carried out on district heating networks expansion in the literature. The present article develops a methodology to find the best district heating network expansion strategy under a set of given constraints. Using a mixed-integer linear programming approach, the model developed optimises the future energy centre operation by selecting the best mix of technologies to achieve a given purpose (e.g. cost savings maximisation or greenhouse gas emissions minimisation). Spatial expansion features are also considered in the methodology. Applied to a case study, the model demonstrates that depending on the optimisation performed, some building connection strategies have to be prioritised. Outputs also prove that district heating schemes’ financial viability may be affected by the connection scenario chosen, highlighting the necessity of planning strategies for district heating networks. The proposed approach is highly flexible as it can be adapted to other district heating network schemes and modified to integrate more aspects and constraints.

QI y,s,t Income from heat sales (£)

Introduction
Developing or expanding a district heating network (DHN) requires a good implementation strategy to take into account the economical and environmental impacts associated with district heating. Indeed, 12% of the EU27 heat supply is already provided by DHNs [1] and the potential of DHNs for the European future has been widely highlighted [2,3]. As a consequence, systematic design and operating procedures are needed for an efficient use of DHNs. Building a new network involves a deep understanding of the network physics, good planning, design and operation strategies, but expanding a network is even more challenging. Indeed, expanding a network requires taking into account the existing network specifications as well as targeting a good expansion strategy with an appropriate time horizon. In these circumstances, having a good investment and expansion schedule generates optimal revenues while ensuring that energy resources are used efficiently to supply heat in an urban environment.
A large number of academic papers already deal with DHNs modelling techniques and develop methodologies to design networks from scratch. These articles focus either on technology selection and sizing [4,5], on the design of essential network parameters (optimal temperature, pipes selection, etc.) [6][7][8] or on the network operating parameters [9,10]. However, literature dealing with DHN expansion, and especially with the associated modelling aspects, remains scarce and is mainly owned by private companies. Firstly, most of the methodologies developed seek to identify the areas that would be suitable for extension, using specific metrics as well as geographic information systems [11][12][13].
Furthermore, strategic energy masterplans and company reports present interesting methodologies for network expansion. For instance, the London areas of Redbridge and Bunhill-Shoreditch identify the subareas that would be economically viable for an expansion by performing a cost assessment and present a set of possible extensions from a cluster network [14,15]. However, they do not detail the investments or extension phasing strategies. Other feasibility studies may include a more detailed phasing of the extension.
For example, in the London Borough of Islington, the potential future heat sources and the thermal stores associated are sized. An optimum network layout is also determined [16]. Taking into account various key drivers and scenarios, several extension planning strategies are obtained for the area, which highlights the importance of considering an incremental approach. However, the modelling techniques used in the study are not presented either.
In the academic literature, optimal sizing of heat sources and installation phasing in the context of DHNs extension has been considered by Mojica-Velazquez [17]. In such work, an optimisation framework is presented to find an optimum investment schedule, considering several heat technologies and implementing a nonlinear model. Lambert et al. also considered the investments phasing. Applying first a classical method to optimally size the heat network, their multi-stage stochastic programming approach was then used to simulate the optimal marginal expansion of an existing network [18]. Their sequential decisionmaking approach could be further enriched to take into account production facilities, maintenance costs or hydraulics.
Finally, it can be interesting to look at capacity expansion and capacity planning problems in other industrial and scientific domains, especially for power production. For example Saif and Almansoori studied the capacity expansion of a desalination and power supply chain, using the General Algebraic Modelling System (GAMS) software and considering various capacity expansion scenarios. Their optimisation model uses a multi-period mixed integer linear programming (MILP) approach [19]. Henggeler-Antunes et al. also developed a multi-objective MILP approach to provide decision support in the evaluation of power generation capacity expansion policies. Considering expansion costs as well as environmental impacts, their model also takes into account demand side management aspects [20].
Although a very small number of articles detail the approach used, most of the documents dealing with network expansion are owned by private companies [14-16, 21, 22]. They give few details on the modelling and optimisation approaches applied. However, capacity expansion problems, which have been studied in other domains, apply methodologies that may be reused in DHN expansion problems. In these circumstances, the model developed in this paper attempts to address the following questions: 1. Focusing on the energy centre design, is there an optimal expansion strategy to maximise the network's cost savings or decrease its environmental impact?
2. In parallel to capacity expansion analysis, what operating strategy of the heating technologies selected can help achieving the same goals?
The paper is structured as follows: the mathematical formulation applied to design and operate the future energy centre with DHNs, as well as the inputs used in the model are detailed first. The methodology developed is then applied to a case study to determine the best expansion strategy and investment schedule. Finally, the model pertinence and its limitations are discussed.

Methodology
The objective of this paper is to propose a methodology to find the optimal investment schedule for a marginal DHN extension, the buildings considered for the expansion being located at the existing network periphery. To achieve this goal, the model presented hereinafter selects the best mix of technologies to be installed in the energy centre and the moment these technologies should be phased in. It also optimises their operation. The general model structure is presented on Figure 1. Prior to the energy centre design, the future pipe network is sized according to the increase of the demand considered over the project lifetime.
Thus the model combines an optimal network flow and capacity planning problem along with an energy system optimisation approach.
As the spatial network extension mainly depends on the willingness of stakeholders to connect new buildings (soft-engineering constraints), predefined building connection scenarios and suitable pipe layouts were built beforehand. Therefore, the interest of the model proposed is that the future network extension is optimised based on the existing network's capacities and layout as well as on realistic connection scenarios and constraints identified by stakeholders. The model's key inputs include the aggregated diversified heat demand (total heat demand taking into account diversity aspects) as well as the connection strategies identified. The technical and economical characteristics of various heat sources that can be installed in the energy centre (including thermal storage) are also considered. Finally, the existing network's characteristics (existing heat producing units, buildings heat demand and spatial layout), which are case-specific, are included. For simplification purpose, no distinction is made between domestic hot water requirements and space heating requirements within the buildings.
The model developed in this paper is formulated using a mixed-integer linear programming approach. In this context, a specific emphasis is placed on linearising several technologies parameters, as it will be seen in

Spatial network extension
As has been mentioned previously, the spatial network expansion is optimised prior to the energy centre design. The pumping costs required to run the network are also estimated at this stage. The existing pipe layout is already taken into account and possible locations for new pipes are also pre-identified. The pipe layout optimisation is similar to previous work reported in the literature [7,18]. The problem considered is a cost minimisation problem: The network layout is constrained, since pipes can be built in one direction only, they must respect the predefined pipe layout and only one diameter can be chosen for each pipe segment: Hydraulic constraints are kept simple by reusing the work done in [7]. The flow velocity, the heat flow and the pressure drop allowed in the pipes are limited by diameter-related constraints: where klin d is the linear coefficient relating the velocity and the pressure drop, following the approach selected by Pirouti et al. [7]. This approach combines Equations (8) and (9) to express the pressure drop as a linear function of the heat flow velocity, which leads to a slight overestimation but keeps the optimisation structure linear.
where the friction factor f i,j,d is approximated by the Haaland equation: In this equation, k' represents the pipes roughness and Re the Reynolds number associated to the heat flow.
Finally, the heat flow at each node depends on the heat consumption at this node and on the heat losses: where the network losses are estimated according to Li et al. description [8]: k is the average heat exchange coefficient (k = 0.3-0.5 W/m 2 . • C) and β is the heat loss coefficient of the district heating network attachment (β = 0.15). The main assumptions used to model pipes are detailed in Appendix A.
Using the pipes layout optimisation results along with Equation (12), the overall pumping power required to run the network can eventually be calculated from: where η P U M P is the network pumps efficiency.
Although the previous approach overestimates the pumping power required to run the network, it provides the proper sizing for future pipes. This step could be replaced by a more detailed hydraulic analysis, similar to what can be found in feasibility studies. The pumping power estimation is assumed to be reasonable as a basic input. The next paragraphs detail the equations used in the energy centre design model.

Technology sizing and working constraints
Technology sizing and working constraints are a set of rules defining which technologies can be installed in the energy centre and how these technologies should be operated. Firstly, the heat output of each technology is limited by its size: This equation applies to all heating technologies that may be installed in the energy centre: CHP plants, heat pumps, natural gas boilers and biomass boilers.
Once a plant is selected, it is installed in the energy centre and starts working: work p,y,s,t Ex p,y q p,y,s,t M · work p,y,s,t As for the pipe layout design, the variables Ex p,y are set equal to 1 at year 1 for the technologies that are initially installed in the energy centre. Heat production is also constrained by the minimum part load allowed for each technology. As may be observed in the following equation, the active working level of a given technology is constrained but can take continuous values (see also Section 3): Finally, constraints on the minimum number of working hours per year or on the minimum number of consecutive working hours allowed are added for model consistency: Thus, all the previous equations shape the heat demand production for each time period and each technology. If necessary, it is also possible to add several constraints to limit the number of plants that can be installed in the energy centre or to impose a minimum annual number of working hours for a given technology.

Heat and electricity balances
Energy balance equations are divided into heat balance and electricity balance. The heat balance equation states that the aggregated heat demand has to be satisfied at each period, the surplus of heat produced being stored in thermal storage facilities: Electricity balance equations can be split into two major electricity balance equations and several constraints on electricity production. The first electricity balance equation states that the amount of electricity produced by combined heat and power plants (CHP) and the electricity bought from the grid are either used to run heat pumps or sold to the grid: Ec HP p,y,s,t + Es G y,s,t The other electricity balance equation states that the electricity used to run heat pumps can be either provided by the grid or by working CHP plants: The amount of electricity produced by CHP units is linked to their heat output by the following linear equation (see also Section 3 and Appendix A): The amount of electricity needed to run heat pumps depends on the quantity of heat they produce and on their coefficient of performance (which is assumed constant): Ec HP p,y,s,t = q p,y,s,t COP ∀p ∈ HP (p), y, s, t The electricity bought from the grid and used to run heat pumps needs to be smaller than the total amount of electricity bought from the grid: Eb G y,s,t Eb G HP y,s,t Finally the electricity produced by CHP plants constraints the amount of electricity that can be sold to the grid, as well as the amount of electricity that can be used internally to run heat pumps: Ep CHP cHP y,s,t Ep CHP y,s,t

Storage constraints
Thermal storage is modelled using a daily approach. Indeed, daily thermal storage is the most common form of thermal storage employed in DH schemes [1,23] and is easier to model using day types. Consequently, it is assumed that thermal stores resume their initial state at the end of the day: the heat storage level is set to be equal to the storage level at the beginning of the day. This approach is similar to the constraints used by Hakarainen et al. and by Oluleye et al. in their respective models [6,24]. The main equation used to model storage is the hourly heat balance, linking the thermal store level, heat inputs, heat outputs and hourly thermal losses: For the first hour of the day, no heat can be discharged from the thermal stores. The thermal storage level is also constrained by the size of the thermal store:

Costs calculations, resource use and greenhouse gas emissions
The overall cost function Cost Obj, detailed in 2.6, is obtained by combining different cost variables.
First, the annualised investment costs are calculated using fixed and variable investment costs as well as the annuity factor. Considering annualised investment costs allows to spread the investments over the Cost an y = p (C p · Ex p,y + Ex p,y · e p · V ar p ) · An p Annual operation and maintenance costs are calculated with the following equation: Resource use and prices have to be taken into account when calculating the amount of natural gas, electricity and biomass needed by the technologies considered. Since heat pumps and biomass boilers are subsidised through the UK non-domestic RHI scheme (see also Section 3), resource equations also involve subsidies for these technologies: Finally, the heat produced in the energy centre is sold to customers, which is modelled by the following equation: Concerning greenhouse gas (GHG) emissions, as the considered technologies have different contributions, it is possible to split the emissions between the on-site combustion of natural gas in boilers, the on-site combustion of biomass in boilers, the indirect emissions associated with the external electricity use and the emissions due to working CHP plants. Regarding CHP units and their associated emissions, it is necessary to separate heat and electricity outputs to properly allocate the total amount of fuel they use. For this purpose, it was decided to follow the UK Department for Environment, Food & Rural Affairs guidance [25] and to select the standardised boiler displacement method to apportion the total fuel to the CHP scheme and the separate heat and electricity outputs. Other standardised methods like the power station displacement method and the DUKES method could also be used for the modelling (see also Section 3).
The non-dimensional B DISP factor introduced in the previous equation is used to calculate the emissions associated with the electricity produced by CHP plants according to the boiler displacement method (see also Section 3).

Objective function
The energy centre can be designed following either a cost savings maximisation approach or a GHG emissions minimisation principle. If the cost savings maximisation is chosen, the overall objective cost function, which minimises the overall investments, is defined by: Cost Obj = min y s t Boiler costs y,s,t + E res y,s,t − QI y,s,t · N days s + Cost an y + Cost OM y On the contrary, if the GHG emissions minimisation approach is selected, the following function has to be minimised:

Model inputs
Given the model sensitivity to the input data, gathering reasonable data or generating them using justifiable trends or projections is crucial. The assumptions used as model inputs are gathered in the next sections.

Combined heat and power
CHP technical and economical characteristics are based on data from manufacturers [26]. Only natural gas CHP with fixed design sizes are selected for the model. Using CHP technical characteristics, it was noticed that their power consumption as well as their electrical output evolve linearly with their part load, as can be seen in Appendix B. Thus it was decided to consider continuous part loads instead of discrete part loads in the model. Concerning the operation of CHP plants, the following assumptions have been made: 1. The model is allowed to select the same CHP capacity unit more than once for the energy centre in order to avoid oversizing problems.
2. As stated previously, CHP plant gas consumption evolves linearly with part load, which is modelled using coefficients a p and b p . Their electrical output as a function of part load is also linearly modelled using coefficients w p and z p .
3. The minimum part load allowed for CHP plants is set to 70% of the design load [4]. 70 % is used as a minimum as the model aims to have a high generation efficiency. This assumption can be easily modified and the minimum decreased given the previous linearisation.

The investment and maintenance costs used in the model are specific to each CHP (see also Appendix
A for the pre-selected capacities and related costs). The investments costs also include the installation costs and are given in 2016 prices.
5. An interest rate of 7 % and a 15-year lifetime are assumed for CHP units [4,26].
6. CHP plants can not be turned on for less than four consecutive hours [4].

Heat pumps
Given the river proximity in the case study (see also Section 4), centralised water source heat pumps are included in the portfolio of technologies selected. As for CHP technologies, fixed capacities are selected for the model (1 and 2 MW th ). The following assumptions are used for the case study analysed in this paper: 1. Given the network temperature, ammonia heat pumps can be installed in the energy centre. For these pumps and for a given sink temperature, the coefficient of performance evolves linearly with the water source temperature [27].
2. The water source temperature is set at 10 • C, which gives a coefficient of performance of 2.86 [27].
3. Centralised water heat pump minimum part load is set at 20 % and a minimum of two consecutive working hours is selected [4].
4. The specific investment costs considered for heat pumps are 1 936 000 £/MW th , with a maintenance cost factor of 6 % [4]. The investment costs also include the installation costs.

Natural gas boilers
Both domestic and non-domestic natural gas boilers are considered in the study. Domestic boilers are solely used in the baseline scenario (see also Section 4). As far as the operation of non domestic boilers is concerned, the following assumptions have been made: 1. It is assumed that non domestic boilers investment costs (including installation and flue evacuation) amount to 121 000 £/MW th with a 18 % maintenance cost factor [4,27].
2. Natural gas boilers can directly generate heat when needed, and no minimum part load is considered.
4. Predefined capacities are considered (from 0.5 to 1.5 MW th ) Concerning domestic natural gas boilers, the following hypotheses are made: 1. An average boiler cost of 41 £/kW th is assumed. This cost was obtained by averaging the boilers costs obtained from manufacturers. An installation cost of 1500-1800 $ per boiler is also assumed [30].

Biomass boilers
As biomass boilers are commonly used in Scandinavian DH schemes [1,31], and since they are subsidised in the UK through the non-domestic RHI scheme [29], biomass boilers are included in the technology portfolio with predefined capacities (1 MW th , 2 MW th ). The following assumptions are used for biomass boilers modelling: 1. The specific investment costs considered for biomass boilers are 133 300 £/MW th with a 16 % maintenance cost factor [27].
2. Since biomass boiler efficiency is considerably reduced as the part load decreases, a minimum allowed part load of 70 % is assumed [1].

A 15
-year lifetime is assumed and estimations of subsidies are based on the UK non-domestic renewable heat incentive tariffs [4,29].

Thermal stores
As explained in Section 2, daily thermal storage is modelled using day types. Simple design considerations including cylindrical thermal diffusion modelling, as well as manufacturers' data are used to estimate thermal stores main properties. The following assumptions are used in the model: 1. Thermal store physical properties are estimated assuming a height-diameter ratio of 3/1 , an insulation thickness of 0.3 m and a fiberglass insulation [1].

2.
A storage price factor of 960 £/m 3 and a dimensionless maintenance factor (expressed as a percentage of the capital costs) of 1 % are assumed [27].
3. An average heat loss percentage of 0.1 %/hour is considered 1 .
3.2. Heat, electricity and resources prices

Heat prices
Heat prices are assumed to be given inputs. As heat prices vary with consumer categories, it is necessary to evaluate the share of the heat demand for each customer type P rop c (public users, new residential customers, etc.). Existing heat prices from previously published feasibility studies are used for this purpose [15]. For example, new residential customers are assumed to pay 8.10 p/kWh heat .

Electricity prices
Hourly electricity import and export prices are generated for each time period using the existing and forecasted wholesale electricity prices as well as all the factors included in the electricity cost breakdown (existing and projected).
1 Losses have been estimated for various thermal store loads (from 1 % to 100 %) and averaged to obtain a representative loss percentage per hour. Since the heat losses obtained with this approach are close from one load to another, averaging them is a reasonable assumption. Although all the factors are included in the import electricity cost breakdown, only the commodity, transmission loss multiplier (TLM), the line loss factor (LLF) and the DUoS commodity are used to calculate the export electricity prices (see also

Resource prices
Resource prices include both natural gas and biomass costs. Natural gas prices are assumed to vary on a yearly basis. Their evaluation is based on DECC 2015 projections, the central price scenario being selected for the model [35]. For example, the projections assume a price of 16.39 /MWh for the first year of the project. A static price of 31 £/MWh heat produced is assumed for biomass [27].

Fuel allocation to CHP plants electricity output
As explained in 2.5, it was chosen to apportion the total fuel to the CHP scheme and to the separate heat and electricity outputs following the boiler displacement method, which is a standardised method used in the UK industry for energy reporting [25]. This method assumes that the heat generated by a CHP unit displaces the heat raised by a boiler with an efficiency of 81% on a gross calorific value basis. To determine the GHG emissions, the method considers that the boiler uses the same fuel mix as the actual fuel mix used by the CHP plant (which is natural gas in that case). Mathematically, the CHP electricity carbon factor can be calculated with the following equation: reasonable assumption, since the coefficients obtained for all CHP units are relatively close one to another [33].
Although other fuel allocation methods could be used to apportion the fuel emissions to heat and electricity outputs, one needs to emphasise that the boiler displacement method leads to the least biased results in the model. Indeed, the other methods tend to maximise the electricity sales, which displaces the emissions produced by the DH scheme on-site and artificially reduces the emissions produced. As it is not a sustainable way to proceed, it was decided to select the boiler displacement method instead.

Resource contribution to greenhouse gas emissions
Although natural gas and biomass are assumed to have a constant carbon factor (respectively 0.185 tCO 2(eq) /MWh [35] and 0.016 tCO 2(eq) /MWh [27]), the electricity grid carbon factor varies dynamically and continuously throughout the year, as result of how the UK fuel mix is monitored. Such variations are included in the model. To generate hourly dynamic factors, the following data are combined: • Forecasted electricity carbon footprints for the next years to come, the number of years selected depending on the time horizon chosen for the problem. DECC reference projections [43]  Multiplying both data, a dynamic carbon footprint profile is obtained and used as input for the model, one electricity grid carbon factor being obtained for each time interval considered over the project lifetime.

Case study
In the next sections, the optimisation program developed is applied to a real district heating network.
The existing network chosen supplies heat to 22 buildings using 2.5 kilometres of pipes. In the existing energy centre, two CHP plants, two natural gas boilers and two thermal stores ensure the continuous supply of heat to the connected dwellings. The extension considered for this network includes 31 buildings preidentified, their connection requiring the installation of 3.5 kilometres of new pipes (see also Figure 2). A

12-years horizon time is chosen for the case study.
For the area considered, the energy centre is designed following either a cost savings maximisation approach or a GHG emissions minimisation principle, as stated in 2.6. A hard urban environment is assumed to evaluate the costs of laying out the new pipes [21]. This assumption arises from talking to local DHNs experts. Since the energy centre is close to a river, centralised water source heat pumps are included in the model. Various connection scenarios are considered and a baseline scenario, in which no marginal extension occurs, is also included in the study.

Heat demand analysis and scenarios
Given the time horizon chosen and the problem size, it is necessary to reduce the number of time periods considered, as it has already been done in many papers dealing with district heating networks [4][5][6]24].
For this purpose, two typical days per month (weekday and weekend) and hourly heat demand data are considered. Since the project covers a 12-years horizon period, this eventually gives 6912 time intervals.
Consumer archetypes are then used to generate consumer profiles for each node [27]. These profiles are adjusted to the aggregated annual heat demand and to each building property (type, area, etc.) and the diversified aggregated demand is finally obtained by applying standard diversity curves [27]. The average heat loss percentage obtained from the pipe layout optimisation is also applied to this demand. Since the connection patterns mainly depend on the willingness of stakeholders to connect, buildings are then gathered into ten clusters, each cluster being connected to the network at a given time. Gathering loads into clusters is mainly based on geographic considerations, each cluster being connected to the network with a main transmission pipe and several distribution pipes.
Following thermal demand data discretisation and the gathering of loads into clusters, connection scenarios are established. The demand analysis is exclusively used to build these scenarios, the connection rapidity and its amplitude varying from one scenario to another. GHG emissions reduction criteria were not considered to build the scenarios, although the emissions reduction potential associated with each scenario is accounted in the model.
Finally three real connection scenarios were built for the model [33]: • A baseline scenario in which the existing network is not expanded and thus customers considered for the expansion continue using domestic natural gas boilers for heating purposes.
• Scenario 1 considers an incremental connection of buildings, from small clusters with a low heat demand to big clusters with a high heat demand. The heat demand evolution is slow at the beginning and then increases rapidly.
• Scenario 2 also considers an incremental connection of buildings, clusters with a high heat demand being connected first. The heat demand increases rapidly at the beginning to almost reach a plateau after a few years.

Results
As explained in Section 2 and Figure 1, the model key outputs mainly include the cost savings achieved and the investments required in each scenario, as well as the mix of plants selected by the optimisation program. Moreover, as the model optimises the energy centre operation, the contribution of each plant to the global heat or electricity production and the GHG emissions associated are returned for each time period.
Thus, for each day type and scenario, heat production profiles, electricity production patterns and GHG emissions curves can be built using the model outputs. For example, Figure 3 shows the heat production profile obtained in the first scenario in case of a cost savings maximisation, at year 5, for a weekday in January.   If a cost optimisation methodology is chosen to design and operate the future energy centre, it is first interesting to observe that both expansion scenarios have a positive cost balance: both expansion scenarios are financially viable (see Table 2). Although they both generate revenues from year 1 (positive costs balance), it can be noticed that the second scenario has higher revenues in the early years, whereas the first scenario's revenues are greater during the late years (see Figure 4). Despite this last observation, the first scenario's expenses are much higher than those required in the second scenario. Both expansion scenarios' expenses are also lower than those required in the baseline scenario (no expansion). Finally, the cost savings associated with the second scenario are 30 % higher than in the first scenario. Apart from the percentage, this result is not sensitive to the time horizon chosen in this case study.
Looking closer at the mix of plants selected by the optimisation program (see also Appendix C), it can be observed that if a lower heat production capacity is installed in the second scenario, the thermal storage capacity installed is the same in both expansion scenarios. Thermal stores are installed earlier in the second scenario, since the heat demand increases a lot during the early years. It is also interesting to observe that the total heat production capacity is oversized in both expansion scenarios. Concerning the technologies themselves, although heat pumps are subsidised through the RHI scheme, the model never installs them in the energy centre. On the other hand, in both expansion scenarios, the model prefers investing in biomass boilers than in natural gas boilers or additional CHP units. Most of the times, it also prefers using these biomass boilers instead of the existing natural gas boilers, which highlights the influence of the biomass boilers' RHI subsidies on the mix of plants chosen. Finally, if natural gas boilers are initially present in the energy centre, the heat they may produce is displaced by thermal stores (summer, mid-seasons) and they are scarcely and sporadically used.
Concerning the energy centre operation, installing thermal stores flattens the global heat production profiles.
In  When a GHG emissions minimisation approach is chosen to design and operate the future energy centre, it is interesting to observe that the investment schedule changes compared to the previous model (see Figure   5). If the annual expenses decrease in both expansion scenarios, the revenues generated by the district heating scheme are also reduced significantly. Indeed, the electricity residue (balance between electricity sales and electricity purchases), which was significant in the previous model is almost equal to zero and sometimes even negative. Therefore, only heat sales generate revenues and both expansion scenarios cost savings have decreased considerably (see Table 3). While the second scenario, with a positive cost balance, is still financially viable, the first scenario is not attractive to invest in any more, as its cost balance is negative. However, when considering GHG emissions, the second scenario's emissions are higher than the In conclusion, if a cost optimisation approach is selected to design and operate the future district heating network, the second scenario should be chosen, as its cost savings are 30% higher compared to the first scenario. On the contrary, if the energy centre is designed following a GHG emissions minimisation principle, the first scenario should be chosen, since its emissions are the lowest at year 12. However, in that case, only the second scenario seems financially viable, although this scenario has the highest GHG emissions. One should emphasise that the results obtained in this section are case specific and do not include carbon costs or credits. They may differ depending on the expansion studied and on the network's original properties.
However, as it will be further discussed, the pertinence of several modelling aspects is already highlighted by this example.

Discussion
The results obtained in the previous section bring into light the influence modelling hypothesis and design strategies can have on the mix of plants finally selected. Depending on the considerations made and the objective function selected to perform the optimisation, the results vary significantly from one model to another.
Overall, cost savings maximisation models have shown that anticipating the capacity expansion as well as oversizing the heat production capacity can be beneficial. Indeed, it could seem appropriate to invest in smaller capacities that would better match the heat requirements. However, this would lead to a lower revenue and the energy centre would have less operating flexibility (the same model run with a mix of plants manually preselected always leads to lower cost savings in both expansion scenarios). In addition, the influence of the technologies already installed in the energy centre has to be emphasised. While the model reuses them when a cost optimisation is performed, they are almost never used when the model tries to minimise GHG emissions, as the existing technologies run on fossil fuels. Consequently, if the energy centre was designed from scratch, without any existing technology and knowing the connection scenario, the model would certainly invest in different technologies and higher revenues would be obtained. Furthermore, cost optimisation models have shown that despite subsidies, it is more interesting to invest in CHP units and export the electricity they produce to the grid. Nevertheless, electricity export rates would need to be studied carefully before making any decision, as rates are usually influenced by many third parties in the energy market. In both expansion scenarios, the model tends to invest first in CHP plants , biomass boilers being installed later. Moreover, the effect of subsidies depends on the connection strategy, a large increase of the demand in the early years stimulating the installation of biomass boilers (Scenario 2). Including thermal stores in the model also enhances the utilisation of biomass boilers. Heat pumps remain too expensive despite subsidies. The optimisation program never selects a heat pump when a cost optimisation is performed.
Finally, the importance of thermal stores for the financial viability of district heating schemes has been demonstrated. Although heat production profiles are overall flattened by the installation of thermal stores, the model tends to use them as extra-loads to maximise the amount of electricity produced by CHP plants while storing the surplus of heat produced.
Greenhouse gas emissions modelling has also proven to be a major challenge for designing the energy centre while ensuring low emissions on-site. Choosing an appropriate method to apportion CHP fuel to heat and electricity output when GHG emissions are minimised avoids potentially biased results. However, it is important to remember that the methodology chosen places the carbon benefits on the supply side, and only considers a local impact. The time horizon chosen for a given case study and its potential influence on the final results needs also to be carefully considered. In all scenarios, biomass boilers and heat pumps are used to supply base demand. The model naturally invests in biomass boilers, as it only considers on-site emissions. One should not forget that the whole biomass supply chain may have a bad environmental impact, which was not taken into account in the present model. The model also chose heat pumps, the technology becoming particularly interesting from year 7 or 8, when the electricity grid carbon factor has significantly decreased. In any case, a proper life cycle assessment would be required for a more detailed GHG emissions analysis.
Finally, it has been observed that the strategy chosen for connecting the heat clusters strongly influences the mix of plants selected as well as the energy centre operating strategy. In the cost optimisation models, it has been seen that the second scenario, in which big clusters are connected earlier, generates higher and quicker revenues than the first scenario. In addition, lower investment costs are required for this second scenario. When GHG emissions are minimised, it has been observed that the connection strategy chosen differently impacts the financial viability of both expansion scenarios.

Conclusion
The objective of this paper was to find an optimum strategy to expand an existing district heating network, focusing on the energy centre design and on its operation. Applied to a case study, the model presented in this article successfully selected and operated an appropriate mix of technologies in the energy centre, either maximising the network's cost savings or minimising its greenhouse gas emissions. The model outcomes reemphasise the need of using optimisation programs to study district heating network expansions. They also highlight the necessity of carefully phasing the network expansion with effective building connection strategies. The methodology developed can be applied to other networks considering a marginal expansion. Further work should consider to include a sensitivity analysis on electricity, heat and resource prices. Concerning the demand, a sensitivity analysis could also be performed to include demand profile variations. The demand scenarios could also be combined into a single optimisation problem using a stochastic approach. It is observed that the cost of laying out transmission pipes for a hard urban environment evolves linearly with pipe diameter, which has been extrapolated to all pipe diameters [21].

Appendix B. Assumptions used to model CHP plants
The following curve displays CHP plant electrical output as a function of part load. Similar trendlines are obtained for CHP plant power consumption as a function of part load.  As explained in 3.1.1, the model is allowed to install the same CHP capacity unit twice in the energy centre (letters a and b are used to distinguish two CHP units having the same capacity).