**Integrated Optimal Design and Control of Fourth Generation District Heating Networks with Thermal Energy Storage**

**Bram van der Heijde 1,2,3 , Annelies Vandermeulen 1,2,3 , Robbe Salenbien 1,3 and Lieve Helsen 1,2,\***


Received: 27 May 2019; Accepted: 10 July 2019; Published: 18 July 2019

**Abstract:** In the quest to increase the share of renewable and residual energy sources in our energy system, and to reduce its greenhouse gas emissions, district heating networks and seasonal thermal energy storage have the potential to play a key role. Different studies prove the techno-economic potential of these technologies but, due to the added complexity, it is challenging to design and control such systems. This paper describes an integrated optimal design and control algorithm, which is applied to the design of a district heating network with solar thermal collectors, seasonal thermal energy storage and excess heat injection. The focus is mostly on the choice of the size and location of these technologies and less on the network layout optimisation. The algorithm uses a two-layer program, namely with a design optimisation layer implemented as a genetic algorithm and an optimal control evaluation layer implemented using the Python optimal control problem toolbox called modesto. This optimisation strategy is applied to the fictional district energy system case of the city of Genk in Belgium. We show that this algorithm can find optimal designs with respect to multiple objective functions and that even in the cheaper, less renewable solutions, seasonal thermal energy storage systems are installed in large quantities.

**Keywords:** optimal design; optimal control; district heating; district energy systems; genetic algorithm; seasonal thermal energy storage; renewable energy

#### **1. Introduction**

Our energy system is one of the main contributors to the ever-increasing greenhouse gas (GHG) emissions, which calls for the identification of solutions within this sector. In particular, heating and cooling in residential and service buildings contribute to no less than 40% of the total final energy requirements in Europe [1]. Currently, 75% of the heating and cooling demand in buildings in the European Union (EU) (including industrial processes) is met by purely fossil resources [2], while of the remaining fraction 11% is provided by bio-mass, although these are usually polluting wood-stoves. Another 7% uses nuclear energy as a source through electricity and only the remaining 7% of heating and cooling come from 'truly' renewable sources, such as hydro, wind, solar and geothermal power. This makes the heating and cooling sector an ideal candidate to tackle the problem of both energy demand and GHG emissions in an efficient way. An often suggested solution is that of district heating and cooling systems to provide the heat and cold demand using centralised production with a large potential for reusing residual heat sources. Moreover, the most modern district heating and cooling systems (of the so-called *fourth generation*) allow for the inclusion of many low-temperature sources, thermal energy storage (TES) and prosumers that inject heat or cold surpluses back into the network [3]. Dahash et al. [4] provided a comprehensive overview of large-scale thermal energy storage systems, concluding that although not necessarily the most cost-effective, tank and pit storage systems are often the most practical to install. They also found a clear gap in the research towards system integration of these seasonal storage systems in terms of modelling both accurately and in reasonable time.

In particular, Paardekooper et al. [5] calculated that switching to a large share of district heating in the European energy system—including only established technologies—enables a reduction by 86% of the CO2 emissions compared to the levels of 1990 but also that district heating can cost-effectively provide at least half of the heating demand in 2050 in the 14 countries that are studied in the Heat Roadmap Europe projects. This reduction is the result of a switch to renewable and residual energy sources (R2ES), enabled by the heat transport provided by district heating. Moreover, according to Lund et al. [6], (seasonal) energy storage will be needed in a highly interconnected energy system, namely to bridge the fluctuations in the availability of renewable energy sources. They calculated that thermal energy and fuel storage are by far cheaper technologies than electrical (battery) storage.

Although it is important to know the potential of district heating and cooling systems, particularly in combination with TES systems, these previous studies only describe fourth generation district energy networks in general terms. Hence, to realise these innovative energy networks, three challenges need to be overcome. Firstly, the design of systems with large shares of fluctuating renewable energy sources will be much more complicated than that of present thermal networks. Secondly, identifying and implementing the right control strategy for a given network will be harder as well, due to fluctuating energy sources and large shares of energy storing components in the network, including energy flexibility. The third challenge follows from the first and second, namely the fact that the choice of control strategy will influence the optimal design and vice versa. This calls for an integrated strategy in which control and design are concurrently optimised.

#### *1.1. Previous Studies on District Energy System Design*

Söderman and Pettersson [7,8] made a topology optimisation algorithm for district energy systems (DES). The algorithm was based on a mixed integer linear program (MILP) for a district including thermal and an electric grid. Thermal energy storage was included in the optimisation, too. They limited the problem to eight representative time instances, namely typical daily and nightly operation conditions in the four seasons. Weber [9] integrated the optimisation of both design and control of poly-generation systems in DES with different energy carriers, but without considering TES. Again, the temporal detail remained limited. Weber used a bi-level solution strategy, where a master optimisation (evolutionary algorithm) chose the type, size and location of technologies to be installed in the network. The slave optimisation (mixed integer non-linear program) decides the layout of the network and the operational strategy.

Fazlollahi et al. [10] presented a multi-objective, non-linear optimisation strategy for DES including district heating and poly-generation, but without considering large-scale TES systems. They used a problem subdivision similar to that of Weber, where a master evolutionary algorithm varies the design parameters, and the proposed designs are evaluated by an MILP, which optimises the energy flows during 8 typical periods. Fazlollahi implemented an additional layer for the thermo-economic optimisation, and a post-processing step to assess the emissions of the proposed designs. The optimal results were summarised in Pareto-fronts according to system efficiency, total annual system cost and CO2 emissions. In this study, the district heating supply and return temperatures were varied as a function of the ambient temperature.

Other studies combine the entire optimisation in a single mathematical problem, often an MILP. Dorfner and Hamacher [11] used this strategy to find the optimal lay-out and pipe size of district heating networks in Germany. This study omitted the operational aspect, instead only considering peak loads. Morvaj et al. [12] presented a single optimisation problem integrating design, operation and network layout for an urban energy system with 12 buildings. They considered one representative day for each month, averaging the electric and heat load profiles for a whole year. Falke et al. [13]

presented a similar multi-objective optimisation problem as Fazlollahi and Weber, but they considered a rule-based control flow for the operational layer, as opposed to an optimal control strategy.

In a wider energy system context, Patteeuw and Helsen [14] presented an integrated control and design optimisation algorithm for the design of the space heating and domestic hot water production system for residential buildings in the Belgian energy system, assuming a number of scenarios for the composition of the future electricity system. They used a single-layer MILP optimisation algorithm with representative weeks to reduce the temporal complexity of the optimisation problem. However, they found that this approach is very slow. They suggest that a scenario-based optimisation is more efficient than a full optimisation problem in which the scenario parameters are included as decision variables. This suggestion clearly points in the direction of a two-layered optimisation approach.

Lund and Mohammadi [15] presented a methodology to optimise the choice of insulation standard for pipes in thermal networks. Their method is split in two calculation tools: one to calculate different scenarios of heat loss behaviour in the thermal network, and the other where the energy flows in the larger system are optimised using EnergyPlan. An evolutionary design algorithm was coupled to EnergyPlan as the evaluation problem by Prina et al. [16]. Their focus was on the operation of regional energy systems to find both techno-economically feasible, as well as sustainable energy system designs. In a later step [17], they accounted for the long-term investment planning problem, considering the evolution of the price for different technologies and the remaining value of previously installed systems as they are replaced by more modern technologies.

Bornatico et al. [18] used a particle swarm optimisation (PSO) algorithm to optimise the thermal system of a Swiss single residential building (hence no DES or thermal network was considered). They aimed to optimise the size of a solar heating system, including a solar collector, storage tank and auxiliary power unit. In this study, the system was simulated in *Polysun*, coupled to MATLAB for the PSO implementation. Whether Polysun implements a heuristic or optimal control was not specified. The results of the PSO were compared to a genetic algorithm and the results were found to be similar. Ghaem Sigarchian et al. [19] optimised a hybrid microgrid including solar photovoltaic panels and concentrated solar power collectors, an organic Rankine cycle to convert heat to electricity, electric and thermal energy storage and a gas-fired backup generator. Both design variables (in the PSO) and a variable operation (in HOMER) were considered. The objective function was the energy tariff to be paid by the consumers in the network, which had to be minimised. The fitness evaluation function seems to be an optimal control problem implemented in HOMER, although this is not clearly explained.

In conclusion of the previous work, a clear pattern is that the optimisation algorithm is subdivided in two layers, where one layer is aimed at evaluating the operational aspect of a particular design—the lower layer or *slave algorithm*—and the other focusses on the exploration of the design parameter space—the upper layer or *master algorithm*. While there are subtle variations where for instance the slave algorithm also optimises part of the design variables, this general structure holds for most of the above discussed references. Still, a smaller number of studies use a completely integrated control and design algorithm, with a single layer that optimises both operational variables and design parameters. Clearly, this approach represents only a minority in the discussed studies and is only suited for design problems with a limited size and (temporal) complexity.

#### *1.2. Novelty and Contribution*

The aim of this paper is to develop an integrated design and control optimisation algorithm for future district heating systems with large shares of R2ES and seasonal thermal energy storage. This algorithm is illustrated in a fictitious district heating system for an existing city in Belgium and the design results from the optimisation are studied in detail. Note that the focus is on the methodological contribution, rather than on the absolute numbers resulting from the case study.

Compared to pre-existing studies, this paper uses a two-layer approach, focussing on the integration of a higher-resolution full-year optimal control problem (OCP) as the lower-level optimisation layer, with particular attention paid to the high operation variability of future energy systems with distributed energy resources. In order to do so, a Python toolbox called modesto (see Vandermeulen et al. [20]) is used to set up these OCPs. We use an optimal selection of representative days compatible with seasonal thermal energy storage systems to reduce the calculation time. To our knowledge, this is also the first study in which a concurrent design of TES volumes, pipe diameters and heat generation systems is considered, together with a more detailed model of the district heating system.

#### **2. Methodology**

The optimisation framework is conceived as a two-layer integrated optimal design and control algorithm. In Section 2.1, the heat demand for space heating and domestic hot water is calculated deterministically and used as a fixed boundary condition for the algorithm. The slave optimisation is a linear optimisation which determines the optimal energy flows in the network for a given design, including the TES charging behaviour, implemented in modesto. This layer of the algorithm is explained in Section 2.2. The master optimisation is a genetic algorithm which looks for the optimal combination of design parameters, based on a number of objective functions. This layer is described in Section 2.3. Apart from the implementation, this section also summarises the available design choices for the chosen case study, as well as the considered scenarios.

#### *2.1. Case Study*

The optimisation algorithm is illustrated by means of a fictitious DES for the city of Genk in Belgium, called *GenkNet*. Spread over 9 neighbourhoods, 7775 building models were constructed based on geometric data for single family residential buildings. Although the network configuration and the choice of the connected neighbourhoods are hypothetical, the data with which the building models were constructed are real.

The building models are equivalent resistance-capacitance models based on the TEASER FourElement structure (see Remmen et al. [21]). Assumptions on the building materials and wall thicknesses were based on the building age, which was assumed fixed for all buildings belonging to one neighbourhood. The workflow to derive the building model parameters was developed by De Jaeger et al. [22]. The heat demand resulting from space heating and domestic hot water (DHW) production was simulated using a typical meteorological year for Belgium and stochastic occupant profiles as boundary conditions. The occupant profiles contain the space heating temperature set point and the DHW draw-off for every individual building and were derived using the StROBe toolbox (see Baetens and Saelens [23]). An ideal building heating system (neglecting the effects of the heating system temperatures on the heat injection in the building) was assumed. All buildings were simulated during a full year with a 900 s time step using a minimum energy objective (assuming a fixed cost for heat), after which the heat demand of all buildings belonging to one neighbourhood was summed and modelled as a single demand node in the network. The heat distribution network on the neighbourhood level is omitted, which means an underestimation of the total heat losses in the network. Instead, the neighbourhood is represented as a single node, connected to the backbone network through a single service pipe.

The resulting heat demand for the 7775 buildings amounts to 430.5 GWh per year, with an average energy use intensity of 210.8 kWh/m2 per year. Cooling and electricity demand were not considered in this study.

The neighbourhoods were located alongside a central thermal network backbone, as indicated in Figure 1.

**Figure 1.** Network layout for the fictitious case study in Genk, Belgium. The heat demand of 9 neighbourhoods alongside a central backbone connection was aggregated per neighbourhood. Two additional nodes without heat demand, but with the option to install heat sources and/or thermal energy storage (TES) systems are situated at both ends of the backbone.

#### *2.2. Operational Optimisation*

We have chosen for a full-year OCP with a 2 h time-step for the evaluation of the DES performance with respect to a number of objective functions (see Section 2.3). The OCP optimises heat and mass flows in the network, the operation of the heat production systems and the charging behaviour of the TES systems in order to satisfy the heat demand of the neighbourhoods. This optimisation has a single objective function to minimise the operational cost. The choice for a minimal operation cost objective is based on the habit of operating real systems to maximise their profit. Except for a limited number of experimental systems, systems are seldom operated to minimise energy use or CO2 emissions, unless this is linked to additional economic incentives. The choice for an optimal control strategy as opposed to a simulation based evaluation or a rule-based control strategy is justified by the quantification of the maximum potential of every design. This potential is always reached when we assume an optimal control strategy exists and can be implemented, whereas heuristic control strategies might penalise designs that are harder to control. As such, we arrive at a fair comparison of designs.

This OCP is implemented using the Python toolbox modesto (see Vandermeulen et al. [20]). This toolbox is built on top of Pyomo [24] and implements a library of linear optimisation models for common DES components and communicates with an optimisation solver to find the optimal operation strategy. More details about the used component models can be found in Appendix A. All models were either derived from literature or verified by the authors. The optimisation variables considered in the operational layer are the magnitude of heat and mass flows in all components, the thermal output of heating systems, possible curtailment of heat from the solar thermal collectors and the state of charge of the TES systems.

The solution time of the OCP is reduced by using representative days. Van der Heijde et al. [25] have developed a method to optimally select representative days and to restore the chronology such that the original data is approximated as closely as possible. This method is applied here. Based on a number of input time series, specified by the user, an optimal set of representative days is chosen, after which the algorithm determines for each day of the year which representative day it will be represented by. In this work, the chosen input time series were the aggregated heat demand for all neighbourhoods, the solar radiation on a unit surface area, the ambient temperature and the hourly electricity price. This method furthermore makes sure that seasonal effects in the TES systems are modelled accurately. We have limited the representative day selection to 12 days as this was shown to be sufficient to represent the full-year OCP with acceptable accuracy, see the conclusion made by van der Heijde et al. [25].

#### *2.3. Design Optimisation*

The design parameter values are varied in the upper layer. While we attempted to keep the slave optimisation linear, both to guarantee a global optimum and to limit the calculation time, non-linear effects do appear in the master optimisation problem. These non-linearities include: investment costs, which can vary with the size of the installed system; discrete decision variables, such as pipe diameters; and the calculation of the state-dependent heat loss from thermal storage tanks. This last phenomenon is caused by the use of a TES model in which the heat loss depends on the actual state of charge. On the other hand, this loss fraction also depends on the size of the TES system, which would render the optimisation problem bi-linear (see the derivation by Vandewalle and D'haeseleer [26]). To avoid this extra non-linearity, the design variable (namely the size of the TES systems) is treated as a constant parameter in the OCP and it is varied in the master optimisation.

We implemented the design optimisation algorithm as a genetic algorithm in Python using the DEAP (Distributed Evolutionary Algorithms in Python) toolbox [27]. The algorithm uses the NSGA-II (Non-dominated Sorting Genetic Algorithm II) selection operator [28]. Crossovers are handled by a simulated binary crossover operator, and for mutation, a polynomial bounded operator is used. Moreover, the genetic algorithm features a small probability of entirely reinitialising some parameters, which is a variation on the mutation operator. The Pareto-optimal solutions of all generations are stored in a *Hall of Fame*. Every new generation is initialised based on all non-dominated individuals, taken over all previous generations. As such, previous optimal solutions cannot be lost in the course of the evolution. A single optimisation run features 100 generations with 60 individuals. Each newly generated individual has a 95% mutation probability and a 70% crossover probability.

Every candidate design is evaluated as an instance of modesto with a minimal cost objective (see Section 2.2). Infeasible optimisation problems result in a high penalty objective value, such that these designs are not selected in the next generation. With modesto, the optimal control trajectory for all energy and mass flows in the network during one year is computed, such that the operational cost is minimal. The workflow of this genetic algorithm is illustrated in Figure 2.

**Figure 2.** Flow chart illustrating the different steps in the genetic algorithm.

#### 2.3.1. Objective Functions

The design was evaluated for optimality based on three objectives, namely the annual primary energy imported from outside the district energy system—used for the heat pumps, geothermal heating plant, possibly auxiliary heating needed for the production of DHW and network pumping power—, the total annualised costs and the CO2 emissions. The total annualised cost *ca* consists of the annual operation cost *cop*, *<sup>a</sup>* as calculated by modesto, increased by the annualised investment *Ia* and fixed annual maintenance cost *cmaint*, *<sup>a</sup>* :

$$\mathcal{L}\_a = I\_a + \mathcal{L}\_{op,a} + \mathcal{L}\_{main,a}.\tag{1}$$

The annualised investment cost *Ia* is calculated using the capital recovery factor:

$$I\_a = I\_{tot} \frac{(1+i)^\tau i}{(1+i)^\tau - 1}. \tag{2}$$

where *Itot* stands for the total investment, and *Ia* for the annualised investment. *τ* denotes the economic lifetime of the technology and *i* is the interest rate, taken as 3%, assuming a public investment on a long term. This is in line with multiple studies, such as that of Möller and Werner [29], Nussbaumer and Thalmann [30] and Steinbach and Staniaszek [31]. The annualisation calculation using the capital recovery factor assumes that every component is replaced by an identical system at the end of its lifetime, and at the same investment. As such, variations in technology prices over time are neglected.

For maintenance (*cmaint*, *<sup>a</sup>*), the annual contribution is estimated to be a fixed fraction of the initial investment. This fraction, as well as the typical economic life time of various technologies, were derived from the EnergyPlan Cost Database [32]. All economic input data for the different technologies considered in this work is summarised in Appendix B.

In order to get a better grasp of the orders of magnitude of the objective functions—that is, annualised total costs, primary energy import and CO2 emissions—we scale them with respect to the total annual heat demand of all neighbourhoods for space heating and DHW. This total heat demand amounts to 430.5 GWh per year. The resulting scaled variables are called the levelised cost of heat

(LCOH, expressed in EUR/kWh), the primary energy import share (PEIS, in %) and the CO2 intensity (in kg CO2/kWh).

#### 2.3.2. Design Choices

The *GenkNet* case has a total of 9 neighbourhoods, 1 industrial node and 1 node with additional heat generation systems but no demand. The design exercise is left very open; the optimiser has to choose how many renewable resources for heat generation are installed at every node, where the TES systems are installed and how large they should be, and how much backup power is needed. The available design choices for the TES systems are listed in Table 1, the solar thermal collector (STC) arrays in Table 2, and the backup heat pumps and geothermal heating plant in Table 3. The maximum volume corresponds to the largest pit and tank TES systems currently found in literature. The maximum STC area corresponds to the available south-oriented roof area of the buildings in the neighbourhoods, however without accounting for previously installed systems, such as PV panels. At node *ThorPark*, a larger area is assumed to be available for the installation of an STC array.

**Table 1.** Available design choices for TES systems at the different nodes in GenkNet. All numbers are expressed in m3.


**Table 2.** Available design choices for the installed STC array area per node. All nodes except *ThorPark* consider the total available South-oriented rooftop area of the considered buildings in that neighbourhood. All numbers are expressed in m2.


**Table 3.** Design choices for the nominal power of the central heat generation systems. The power is expressed in MW. The abbreviations "geo" and "hp" stand for geothermal heating plant and air source heat pump, respectively.


In addition, the pipe diameters are also design decision variables. For the available diameters, the reader is referred to Tables A2 and A3. The smaller diameter pipes are implemented as twin pipes (up to DN 200), the larger pipe diameters as compound pipes. The investment costs for these pipes are discussed in Appendix B.3. The available diameters were derived from IsoPlus [33]. It is also an option to install no pipe at all at a specific network edge; this choice is represented by a 0 m diameter pipe.

All scenarios share the same network layout, as shown in Figure 1. Whereas the choice for the size of the heat pumps and the geothermal heating plant is handled by the optimisation algorithm, the availability of excess heat is fixed at 10 MW, constantly available throughout the year. Whether this resource is utilised or not is a question of the decision of the district heating connection between the node *GenkZuid* and the rest of the network, that is, is it worth the investment to make a connection to the industrial area from the city or not.

#### 2.3.3. Scenarios

While most of the boundary conditions are fixed for the design algorithm, two of them are varied discretely and deterministically to establish their influence on the results, leading to a number of scenarios. The first boundary condition is the nominal temperature level in the network. Four options are available:


Hence, most scenarios use a 20 K nominal Δ*T* with rather low supply temperatures, whereas the last scenario uses medium-high temperatures with a 40 K nominal Δ*T*.

The second scenario parameter is the cost of heat from the industrial excess heat source in the most southern node of the network. The heat prices considered are:


These excess heat costs are substantially higher than the ones discussed by Doraˇci´c et al. [34], but they are chosen to be in line with the expected cost for an industrial company that needs to invest in a connection of its processes to a district heating system.

The different combinations of the two scenario parameters lead to a total of 16 optimisation runs to be performed. However, the focus will be on how the scenarios deviate from the reference scenario—55/35 ◦C with 15 EUR/MWh excess heat cost—leading to a total of 7 scenarios to be studied in detail.

#### **3. Results**

The emphasis of this paper is on the methodological contribution, namely the integrated design and control optimisation algorithm. The results presented in this section should mostly be interpreted as a proof of concept, rather than in absolute numbers.

#### *3.1. Reference Case Results*

As mentioned before, the case with a 55/35 ◦C temperature regime and an excess heat cost of 15 EUR/MWh is chosen as the reference case. This section shows a selection of visualisations of the optimal design results in order to make more sense out of the large amounts of output data. Firstly, we focus on the higher level, using only design parameters and yearly aggregated outcomes, but in a later stage we will also zoom in on the results on smaller time scales.

#### Genetic Algorithm Outcome

The solutions resulting from the genetic design algorithm, ranked by the objectives of LCOH and PEIS are shown in Figure 3. On the one hand, this figure shows the spread of all considered designs that turned out to be feasible in terms of their LCOH and PEIS values. The blue dots show the Pareto-optimal solutions, that is, the solutions that dominate the solution space. The three black markers respectively show the solutions with the lowest PEIS, the one with an approximately average PEIS measured over all non-dominated solutions and the one with the maximal PEIS. These three specific solutions will be treated in more detail in the next subsection. Note that the minimal and maximal PEIS solutions also represent the respective maximal and minimal LCOH solutions.

**Figure 3.** Scatter plot of all genetic design algorithm solutions for the reference design case.

Zooming in on the Pareto-optimal solutions only, we can plot the contributions of operational costs, investment and maintenance costs to the LCOH. This cost breakdown is presented in Figure 4, where we see that the contribution of the annualised investment takes up the largest share of the total annual costs of the system. As expected, as a larger amount of energy is imported from outside the network, the operational costs (representing in this case the cost of electricity use) increase linearly.

**Figure 4.** Different contributions to the total annual LCOH for the reference case, namely that of annual maintenance, operational cost and annualised investment cost. The annualised investment has the largest contribution.

The contribution of the different technologies is shown in Figure 5, showing indeed that the higher investment for the lower PEIS solutions is largely caused by a larger installed amount of TES and STC systems. The bar chart furthermore shows that the investment in auxiliary heating plants (i.e., the heat pumps and the geothermal plant) combined remains more or less constant. This does not mean that the installed auxiliary power remains constant as well, given the different prices per unit of thermal power for the two technologies. Finally, Figure 5 shows that the investment in the transport pipes remains more or less constant, too. On closer inspection (not shown here), we find that the lower PEIS solutions have marginally larger investments in the network, corresponding to wider pipes on average.

**Figure 5.** Bar chart showing the contribution of different technologies to the levelised annualised investment costs. The cost of the auxiliary domestic hot water (DHW) boilers has been omitted since it is the same for all solutions shown.

Figure 6 shows the distribution of the chosen sizes for the STC and TES systems in the network, split by network node and *energy range*. The energy range is a categorisation of the Pareto-optimal solutions based on their PEIS. The range of PEIS values for all Pareto-optimal solutions is split in three equal parts, denoting high, mid and low energy import. The plot shows three violin plots per neighbourhood and parameter, indicating the approximate distribution, as if the design parameters were distributed following some probability function. The plot gets wider where there is a denser distribution of observations for that parameter. The actual parameter values are indicated by the black dots inside the violin plots, which shows that the distribution is in fact very sparse, with a few dense spots here and there.

**Figure 6.** Distribution of solar and storage systems per network node. The Pareto-optimal solutions are split into high, mid and low energy import solutions. The left graph shows the maximum STC area per neighbourhood as the grey diamond.

The solutions with lower PEIS are expected to have a larger share of solar energy, which is only possible if enough STC area is installed, as well as TES volume. For example, *Boxbergheide*, *ThorPark*, *WaterscheiGarden* and *ZwartbergNEast* show this evolution of decreasing STC area with increasing imported energy. For the other neighbourhoods, the spread stays largely the same, although in most cases the average shifts to lower values. Only for *OudWinterslag*, no such shift can be distinguished. Another interesting observation is that more often than not, almost the maximum available area is exploited for the installation of STC systems.

Looking at the TES volumes, such a trend is not so easy to find. In the case of *OudWinterslag*, *TermienWest*, *WaterscheiGarden*, *Winterslag* and *ZwartbergSouth*, we can see that at least the low energy range has the largest storage volumes. However, if we compare this graph with the map (Figure 1), one interesting trend is that larger storage tanks tend to be installed as close as possible to locations where auxiliary heating plants (i.e., heat pumps and geothermal) are available. This is more obvious for the neighbourhoods of *ZwartbergNWest* and *Winterslag-Boxbergheide*. *ZwartbergNEast* seems to act as the storage hub for *ThorPark*, given the smaller distance between these nodes than the distance between

*Waterschei Garden* and *ThorPark*. The reason for the predominantly large volume at *ZwartbergNEast* can be explained by the limited storage potential at the *ThorPark* node.

When the distribution of the nominal power of the production systems is plotted per neighbourhood, we find that the variety in the solutions is rather low. Figure 7 shows that the size of the heat pump in *ZwartbergNWest* is positively correlated with the energy import range and the opposite is true for the geothermal heating plant. The heat pump at *Winterslag* has an almost constant nominal power.

**Figure 7.** Spread of the nominal power of the chosen production systems, again subdivided according to the range of the energy import from outside the network. The heat pump at *Winterslag* is nearly always the same size, at *ZwartbergNWest* the size increases with the range of energy import and the geothermal heating plant decreases with increasing energy range.

A final plot of interest is the spread of the sum of all storage volumes compared to the total STC area in the network. This graph is shown in Figure 8, and a sigmoid distribution appears clearly. For lower STC areas, the storage volume appears to increase linearly with the STC area, until around <sup>150</sup> × <sup>10</sup><sup>3</sup> <sup>m</sup>2, where the storage volume starts increasing very rapidly. As soon as the maximum storage volume is reached, the STC area can still increase, but mostly at a higher cost without much reduction in terms of import share.

**Figure 8.** Distribution of total STC area compared to the total storage volume. The scatter points are coloured according to the LCOH, with darker values representing a lower LCOH than the lighter values. The larger the dot, the larger the import share (and thus the larger the grid dependence of the network). Note that the axes are limited by the minimum and maximum allowed design sizes.

#### *3.2. Detailed Study of Highlighted Reference Solutions*

In Figure 3, three particular solutions are indicated with a black marker. Two solutions respectively at the upper and lower end of the Pareto front were chosen, and one solution that is closest to the average PEIS of all Pareto-optimal solutions. For these designs, we investigate the heat sources in the network and the energy storage levels in more detail.

#### 3.2.1. Maximum LCOH, Minimum PEIS

The first highlighted solution is the one with the lowest primary energy import from outside the network, but the highest LCOH. In the graph of the heat flow rates (the upper subplot in Figure 9), the black line indicates the part of the heat demand of the neighbourhoods which is supplied by the district heating network. We see that the excess heat production is used as a base load throughout most of the year, except for the Summer period, where most of the heat demand is met by a combination of solar thermal power and discharging of the TES systems.

**Figure 9.** Full year time series for minimum PEIS design.

The geothermal plant is runner-up after the excess heat injection. Its operation is fixed by the installed nominal heating power and the fixed operation schedule throughout the year. Additional gaps are filled mostly by the heat pumps and by discharging the storage systems. The STC panels inject heat whenever they can, and surpluses are stored for later use. This is visible from the part below the zero power level, which indicates the charging behaviour of all tanks combined.

When the net power graph is studied in a shorter time range (see Figure 9), we can clearly see that power exceeding the neighbourhood heat demand is stored in the TES systems.

Around 40 GWh of TES is installed, and a clear seasonal charging pattern is obvious from the lower subplot in Figure 9. At the start of March, the storage tanks are completely empty, and during Spring and Summer they are gradually charged with energy until the maximum storage level is reached in November. Then, the storage is quickly depleted until the cycle repeats. A result of the genetic algorithms optimisation is that the full storage volume is used. Solutions where the charged energy profile does not use the full available range are dominated by more efficient designs.

An interesting result is that most storage tanks are used with a more or less similar charging pattern, witnessed by the evenly spread contributions of different systems. Only the smaller *ThorPark* TTES system is overshadowed by the other systems.

Figure 10 zooms in on a winter, spring and summer week to show a more detailed image of the heating power flows in the system. In the winter week, we see that the geothermal plant is working continuously, largely assisted by the heat pumps at their full power.

**Figure 10.** Detailed heating power plot of the minimum PEIS design for a winter, spring and summer week.

The remaining heat demand is met by discharging the TES systems. During the spring week, we see the intermittent behaviour of the geothermal plant (which is fixed by the off- and on-periods as determined by the user), but also the high power coming from the STC systems. Note that the TES tanks are mostly charging and only limited TES discharge is needed to fill some gaps in the heat demand. The heat pump is only on for higher demand levels. Note that the maximum amount of excess heat is being injected during both weeks (winter and spring). Finally, the summer week shows a very low heat demand, and while all "auxiliary" conversion systems are off, the demand is met by discharging the storage and power from the STC systems. Note how the surplus of solar heat is charged to the TES systems.

#### 3.2.2. Intermediate LCOH and PEIS

The intermediate solution (see Figure 11) shows a rather different picture: whereas the previous solution was characterised by heat flow rate peaks over 200 MW, here the peak powers are more limited. The positive peaks do not exceed the maximum heat demand, which seems to suggest that the network has been designed to accommodate the maximum heat demand and not more. A slightly smaller geothermal system seems to be suggested by the power graph, but this is compensated by larger heat pumps. The lower positive peaks show that a smaller STC area has been installed.

**Figure 11.** Full year time series for intermediate PEIS design.

The size of the TES systems is also considerably smaller, peaking at energy levels of 15 GWh. The seasonal behaviour is still apparent, but the charging behaviour is more gradual than that of the minimum PEIS solution. In addition, the faster charging and discharging cycles on top of the seasonal behaviour suggest that in this system, weekly and diurnal storage cycles have a more important role to play, contributing to the lower operational cost of this system. Finally, only four large storage systems are installed here, with the other TES systems playing only a marginal role. Further study will have to prove whether these TES systems will be completely removed if we leave the genetic algorithm to search even more generations with larger populations, or they are actually needed and worth investing in. Whereas the previous solution was characterised by similar charging patterns, the large storage tanks' storage behaviour is clearly shifted in time. To illustrate this, *Boxbergheide* and *TermienEast* only start charging in June, whereas *ZwartbergNEast* already starts filling up in March and remains at a more or less constant level from May until end of December.

Figure 12 again shows a more detailed picture of the heat flows. The behaviour looks largely the same as in Figure 10, the difference being in the scale of the graphs, which shows that the maximum power peaks are considerably lower than those in the minimum PEIS design. This is a result of the smaller STC systems installed here.

**Figure 12.** Detailed heating power plot of the intermediate PEIS design for a winter, spring and summer week.

#### 3.2.3. Minimum LCOH, Maximum PEIS

The final highlighted design is that with a minimal LCOH (see Figure 13), but with the largest demand for primary energy from outside the network. Again, there are clear differences with regard to the previously discussed solutions: in this case, the excess heat injection is almost continuously at the maximum power level, whereas it was pushed out by solar power in the previous solutions. In this case, the installed STC area is very limited, which is witnessed by the absence of large positive power peaks.

**Figure 13.** Full year time series for maximum PEIS design.

However, the TES volume is hardly smaller than that of the intermediate solution. It appears that a combination of the almost constant excess heat supply with the limited solar power injection makes it economically interesting to still have a decent amount of storage available in the network. In comparison to the intermediate solution, the average SoC is lower, which means smaller heat losses from the storage systems. Again, even though there is a clear seasonal charging pattern, there is an emphasis on shorter charging cycles to minimise operational costs.

Figure 14 shows the behaviour of the maximum PEIS system on a smaller time scale. Note the substantially lower solar thermal peaks compared to Figures 10 and 12. Another evident difference is the continuous injection of excess heat into the system, even during summer. Finally, we see that the heat pump is used more often during the spring week and even during the summer week.

**Figure 14.** Detailed heating power plot of the maximum PEIS design for a winter, spring and summer week.

#### 3.2.4. Comparison of Aggregated Results

To conclude this section on the three highlighted design solutions, Figure 15 summarises the heat delivery by the different heat sources in the network and compares these number to the total demand and storage losses. The difference between those two columns (sources vs. demand and losses) only differ by the amount of heat lost in the network. This figure makes clear that these network heat losses only make up a very limited fraction of the total heat balance. Although this result seems alarming at first, it can be explained by the fact that only the transmission network is modelled and the distribution network is omitted. Typically, these transmission pipes are much more efficient because of their wider diameter and the relatively low loss surface compared to the high mass flow rates that flow through these pipes. In addition, further analysis of the charging and discharging behaviour in these three highlighted solutions shows that an annual round-trip efficiency between 83% and 90% is achieved with the current storage model.

**Figure 15.** Summary of heat delivery for the three highlighted solutions.

#### *3.3. Influence of Network Temperatures*

This section compares the results of different temperature scenarios to the reference scenario with a 55/35 ◦C temperature regime. Figure 16 shows the differences between the resulting Pareto-optimal solutions with respect to the LCOH and PEIS objective functions.

**Figure 16.** Comparison of the Pareto-optimal solutions for the different temperature scenarios. The excess heat cost is fixed at 15 EUR/MWh.

There is an obvious shift on the PEIS-axis, correlated with the average network temperature. This shift can most probably be mainly attributed to an increase in the COP (Coefficient of Performance) of the heat pumps and the geothermal plant with lower average network temperatures. The heat losses will also diminish with lower average temperatures, but given the limited influence of these losses, even more so with regard to the primary energy import into the network, it is expected that their role in the PEIS shift is of minor importance.

The 75/35 ◦C case is the only scenario with a 40 ◦C temperature difference, and this translates into a lower average LCOH. This is a result of the smaller pipes required for the same heating power transport compared to a system with a lower temperature difference, and by a more efficient utilisation of the same TES volume or STC area. The 20 ◦C Δ*T* scenarios are characterised by very similar LCOH ranges, but greatly differ in PEIS, where the reduced heat loss from the network and the increased COP of heat pump-based conversion systems are thought to be the main reason for these differences.

Considering the gas price in Belgium is currently between 0.036 and 0.039 EUR/kWh depending on the type of customer (Belgian Commission for Electricity and Gas Regulation CREG [35])—including the commodity price and the network cost, excluding taxes and levies, and disregarding the costs for a gas boiler or the loss of efficiency of a realistic heating system—we see that the current business as usual situation is already close to the LCOHs encountered in the future DES design. Compared to those numbers alone, already many of the proposed DES designs are actually cheaper than individual gas-based heating. At least, the system costs are clearly in the same order of magnitude. To make an honest comparison, the cost of the distribution networks (not considered in this study), but also the investment and maintenance of gas boilers would have to be included in the analysis. Hence, these results need to be interpreted carefully.

Figure 17 shows the Pareto fronts in terms of CO2 intensity and LCOH for the different temperature scenarios. The carbon intensities for the network vary between 0.02 and 0.045 kg/kWh. For comparison, the specific CO2 emissions for combustion of natural gas is 0.2 kg/kWh [36]. Even without accounting for the efficiency of a natural gas boiler, it can be seen that the CO2 emissions of the studied designs are considerably lower. Clearly, Figure 17 is almost identical to Figure 16, except for the different *x*-axis. We refer to Appendix A.2.5 for a description of the time variation of the CO2 intensity of the electricity consumption. Another important cost to include would be the current carbon tax. Moreover, this type of analysis could be used to determine which carbon tax levels would be needed to push the market towards more renewable technologies.

**Figure 17.** CO2 intensity versus LCOH Pareto fronts for different temperature scenarios.

Figure 18 shows the relationships between the installed TES and STC sizes and their influence on LCOH and PEIS. In this figure, we see a more or less linear correlation for all temperature levels except the reference scenario. The 75/35 ◦C scenario has a slightly smaller TES volume distribution, but the spread on the STC area is more or less the same.

**Figure 18.** Correlation plot between total TES volume and STC area for GenkNet for different network temperatures. All plots have a fixed excess heat cost of 15 EUR/MWh.

#### *3.4. Influence of Cost of Excess Heat*

The differences in the Pareto fronts of scenarios with different excess heat cost—but with fixed network temperatures of 55/35 ◦C—are shown in Figure 19. The differences between the fronts are very subtle, and on the first sight they seem to be mostly coinciding. On closer inspection, the front with 20 EUR/MWh for excess heat is usually on top, indicating a slightly higher LCOH for the same PEIS, although in the higher LCOH range this is not always the case. In addition, the Pareto front of this highest excess heat cost stretches out further to the lower right than the other fronts. This means that lower cost solutions can be reached by importing more energy when the excess heat cost is higher, however the differences are very small.

**Figure 19.** Comparison of the Pareto-optimal solutions for the different excess heat cost scenarios. The network temperatures are fixed at 55/35 ◦C.

These results suggest that the excess heat cost is of minor importance to the solution, at least for the cost range studied here. Indeed, the cost of excess heat is small compared to the total levelised costs encountered in the solution space. We expect to see larger differences when the excess heat cost approaches or exceeds the current system LCOHs. Moreover, the importance of the investment to connect the industrial excess heat sources is not to be forgotten. Because the availability of excess heat requires a district heating connection with a substantial investment cost, we expect that the point at which excess heat is no longer included in the optimal design will be already encountered at a lower excess heat price than the current system LCOH levels. Currently this tipping point at which excess heat is no longer chosen as a source of heat has not been encountered yet.

The main investment cost consideration influenced by the price of excess heat is the investment in the connection between the backbone and the node *GenkZuid* at which the excess heat is available. Therefore, we study the chosen diameter of this connection. Figure 20 shows the diameter for different energy ranges and excess heat costs, and it is clear that for the highest excess heat cost, the installed diameters are substantially smaller. Furthermore we see that the low energy import solutions always have a wider connection with the excess heat node than the ones with a higher energy range. Strangely, again we see an increasing trend in the diameters of the excess heat connection for the three lowest excess heat prices, which appears to be contradictory to the intuition that one would invest less to get access to a more expensive commodity, however also operational costs play.

**Figure 20.** Influence of the excess heat cost on the diameter of the connection between the excess heat injection node and the rest of the network.

In summary, the cost of excess heat does not have a very clear influence on the design choices, except for the diameter of the excess heat node connection.

#### *3.5. Calculation Times*

For the seven calculated scenarios, the average calculation time for an entire genetic algorithm run (with 100 generations of 60 individuals) took 17.7 h. The algorithm was run on a Dell Precision 7920 Tower workstation with two Intel Xeon Silver 4116 processors (2.1 GHz–3.0 GHz Turbo, 12 cores) and 64 GB 2666 MHz DDR4 RAM. Hence, the average operational optimisation (a single evaluation) took 127.4 s, incorporating the fact that 12 processes were run simultaneously. Because the solver usually identifies an infeasible design much quicker, the average modesto calculation time for a feasible design is expected to be higher. However, a large spread on the calculation times was observed, where the problem was typically harder to solve as the system design got closer to its feasibility boundaries—that is, when the components were less oversized.

The fastest genetic algorithm run was performed in 12.0 h, the slowest in 24.9 h.

#### **4. Discussion**

The results have mostly shown that the genetic algorithm is able to select optimised designs using modesto as the evaluation core to calculate the objective function values and using representative days to reduce the calculation time. However, finding generalisable conclusions or rules of thumb based on these results is difficult. Even if we could find general rules based on the presented results, more data for different kinds of networks and additional scenarios would be needed in order to really generalise. Moreover, because of the heuristic nature of a genetic algorithm, we cannot prove the optimality of the results. Nonetheless, the convergence of the Pareto-optimal solutions is an indication that at least a local, and in the best case a global optimum is approached.

In addition, further study would be needed to increase the level of detail of the optimisation problem; for example, currently we don't account for the efficiency of substations that deliver the heat from the network to the end-users. Also the distribution networks in the separate neighbourhoods are not modelled at this moment; this means the investment is considerably underestimated, and also the heat demand might increase due to extra heat losses.

The omission of the building thermal flexibility from the model is expected to have a smaller influence, given the presence of seasonality in all storage SoC profiles, even for the solution with the highest share of primary energy import. Still, it would be interesting to see if the possibility of using the already available sources of energy flexibility (see Vandermeulen et al. [37,38] for a comprehensive discussion of this topic) in thermal networks would influence the design choices made.

One of the objectives of this thesis was to optimise the location of TES systems in a network. From the highlighted solutions and the distribution of the storage sizes, a possible conclusion would be

that TES systems are preferably installed as close as possible to the heat sources in the network—that is, the heat pumps and geothermal plant. Of course, this result is biased since the STC systems also inject heat. But in general, the conclusion that heat should be stored as close as possible to where it will be generated is intuitive as it minimises transport losses.

The simultaneous optimisation of the network pipe diameters adds an interesting dimension to the results. Whereas smaller diameters are mentioned as one of the advantages of future, smart thermal networks (see Lund et al. [3]), the design results with maximal renewable and residual energy shares show that the pipe diameters actually become larger to enable large peak flow due to the redistribution of stored energy in the network and due to the high power injection by the solar thermal collectors.

A final remark must be made about the convergence of the genetic algorithm. The comparison of the solutions with different excess heat cost showed some counter-intuitive results, which could be explained by a possible incomplete convergence of the genetic algorithm. Upon closer investigation of the evolution of the hypervolume indicator (not included in the results of this paper), we see that the increase of the indicator flattens, indicating (near) convergence of the algorithm. Still, sometimes a sudden jump in the indicator—usually because of an accompanying jump in the Pareto front—is observed. Hence, it is possible that the most optimal Pareto front has not been reached yet. To be sure of the convergence, the genetic algorithm should be run with an even higher number of generations. On the other hand, the trade-off must be made between the optimality of the outcome, compared to the additional computation time needed to reach it. A similar trade-off is seen in the solution of MILPs, where an optimality gap is allowed to greatly reduce the solution time, at the expense of a slightly weaker certainty about the optimality of the solution.

#### **5. Conclusions and Further Recommendations**

This paper describes and illustrates an integrated design and control optimisation algorithm. For the evaluation of a single design, a Python toolbox named modesto, which implements a full year linear optimal control problem is used, in conjunction with an optimised representative days method to reduce the temporal complexity of the evaluations. By comparing designs using an optimal control problem, they can be compared objectively with regard to their maximum potential, and the comparison becomes independent of the efficiency of chosen control strategies. This optimal control problem is used as the evaluation core of a genetic algorithm, in which the design parameters are varied. As such, a two-layer design optimisation algorithm is constructed, where the lower layer optimises the operational aspect of the proposed networks, whereas the upper layer varies the design parameters in order to find the optimum with respect to a number of predefined objective functions.

This optimisation algorithm is applied to a fictive district heating network for the city of Genk in Belgium, and the influence of the network temperatures and the excess heat cost is investigated using a number of scenarios. It is shown that the design algorithm is able to efficiently find optimal solutions with respect to multiple simultaneous objectives, and that the proposed systems are competitive with individual natural gas-based heating systems in terms of levelised cost of heating, and even outperform this gas-based non-collective reference in terms of CO2 emissions. However, we were unable to derive clear rules of thumb as to where the thermal energy storage must be installed. A clear result, though, is that seasonal energy storage will be crucial for future energy systems, as ample storage volumes are selected even for the cheapest solutions with the highest primary energy import into the network.

Suggestions for further research would be to add more modelling detail, for example by modelling the heat losses in the neighbourhood distribution networks explicitly (currently, only the transport pipes are modelled without considering the additional heat losses and investment in the neighbourhood distribution grids). Modelling the entire distribution network in detail however would probably result in very complex optimisation problems. Therefore, further study of how a homogeneous distribution network with distributed injection of renewable heat (STC panels on the rooftops of every building) could be aggregated, is required. Secondly, the modesto framework allows for easy addition of more energy conversion and thermal energy storage models. As such, the connection between different energy carriers (i.e., natural gas and electricity) could be strengthened with the addition of for example cobined heat and power (CHP) systems. This would also add a larger variation on the CO2 objective, which is currently very strongly linked to the primary energy import into the network. As a last improvement to the optimisation toolbox, also the network temperatures could be varied, but in a deterministic way. The supply temperature could be varied as a function of the ambient temperature (heating curve), and as a pre-processing step this would maintain the linearity of the optimisation problem. However, this would require a modification of the energy storage models, which currently rely on fixed high and low temperature levels.

**Author Contributions:** Conceptualisation, B.v.d.H., A.V., R.S. and L.H.; methodology, B.v.d.H. and A.V.; software, B.v.d.H. and A.V.; verification, B.v.d.H.; writing—original draft preparation, B.v.d.H.; writing—review and editing, B.v.d.H., A.V., R.S. and L.H.; visualisation, B.v.d.H.; supervision, R.S. and L.H.

**Funding:** This research was funded by VITO grant numbers 1510475 and 1710760.

**Acknowledgments:** The authors want to acknowledge the work of Ina De Jaeger for the curation and analysis of the GIS database of Genk and for providing the building model parameters. We are grateful for the feedback on the results by Jan Diriken and Tijs Van Oevelen. Furthermore, we would like to thank Vincent Reinbold for the fruitful discussions that have contributed to this work. The work of Bram van der Heijde and Annelies Vandermeulen is funded through VITO PhD Scholarships.

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

#### **Nomenclature**


#### **Appendix A. Operational Optimisation Models**

This section provides a summary of model equations used. Because of the modular structure and the complex problems considered, giving a comprehensive overview of all optimisation variables and constraints is not possible.

#### *Appendix A.1. Network Models*

The heat losses from the district heating network are calculated (see van der Heijde et al. [39]) based on a nominal supply and return temperature. These nominal temperatures are also used by all other components in the DES, such as TES and conversion systems. To maintain the model's linearity on the one hand, and to avoid infinite temperature differences at low mass flow rates (in the case of nearly mass flow-independent heat losses) on the other, the heat losses from the pipes are made a linear function of the mass flow rate. The nominal heat loss level is reached around the mass flow rate that corresponds to a pressure gradient of 80–100 Pa/m.

Pumping losses in the network are modelled with a set of linear inequality constraints. These linear segments interpolate between equidistant points on the actual pumping power curve (third-degree function of mass flow rate). Because of the inequality constraints, only branched networks can be currently modelled, and the pumping energy must be represented in the operational objective function, such that high deviations from the inequality constraints are penalised. In this case, the (nonnegative) cost for the electricity to drive the pumps is a part of the operational cost.

As such, the following model equations can be written for a pipe model:

$$
\dot{Q}\_{\rm in} = L d\dot{q} + \dot{Q}\_{\rm out} \tag{A1}
$$

$$
\dot{m}c\_p \Delta T \le \dot{Q}\_{\text{in}} \le \dot{m}c\_p (1 + \gamma) \Delta T \tag{A2}
$$

$$
\dot{m}c\_p\Delta T \le \dot{Q}\_{out} \le \dot{m}c\_p(1+\gamma)\Delta T,\tag{A3}
$$

where *γ* denotes the allowed variation on the temperature difference in the network to allow heat losses without violating the energy balance. We limit *γ* to 5%, but the user can specify this parameter as needed. *Q*˙ *in* and *Q*˙ *out* are the respective heat flow rates at the inlet side and outlet side of the pipe system. Δ*T* is the design temperature difference between the supply and return side, neglecting any temperature differences between the in- and outlet side. *L* is the length of the modelled segment. Note that *Q*˙ is defined positive when heat is transported from the inlet to the outlet of the pipe, and vice versa. Heat losses *dq*˙ are positive when heat is lost from the pipe to the surrounding.

The linearised heat loss model is implemented as:

$$d\dot{q}(\dot{m}) = d\dot{q}\_{\text{nom}} \frac{|\dot{m}|}{0.7233 \dot{m}\_{\text{max}}},\tag{A4}$$

where

$$d\dot{q}\_{\text{nom}} = \frac{T\_{\overline{v}} + T\_{\overline{r}} - 2T\_{\overline{a}}}{R\_{\text{s}}}.\tag{A5}$$

The value of *Rs* is listed in Tables A2 and A3. 0.7233 is a factor that influences the slope of the heat loss curve, such that the nominal heat losses occur in the region around a pressure drop of 100 Pa/m.

#### *Appendix A.2. Thermal Energy Conversion Components*

The energy conversion units models are based on the Energy Hub concept (see Geidl et al. [40,41] or Evins et al. [42]). The conversion of one energy form to another is modelled using a predetermined efficiency factor. While the nominal network supply and return temperature to calculate transport heat losses are fixed, a variation on the temperature difference between supply and return *γ* is allowed to be able to account for heat losses in the network. This implies that also at the production side, the relation between the design temperature difference, mass flow and heat flow rates cannot be imposed strictly, to avoid non-linear equations.

Hence, the following inequality constraints on the mass flow rate and the temperature difference form a general representation of an energy conversion system in the thermal network:

$$
\dot{m}c\_p\Delta T \le \dot{Q} \le \dot{m}c\_p(1+\gamma)\Delta T,\tag{A6}
$$

where *m*˙ is the mass flow rate injected into the district heating system's supply and extracted from the return side by the heat generating system and *Q*˙ is the thermal power injected into the network by the conversion system. This mass flow rate is heated by at least Δ*T*, and at most (1 + *γ*) · Δ*T*. These two inequalities make sure that the conversion unit is not able to inject heat into the network at zero mass flow rate, which would imply an infinite temperature difference. In addition, the maximum heat injection is limited to the nominal power of the conversion system, which is a design parameter:

$$
\dot{Q} \le \dot{Q}\_{\text{max}} \tag{A7}
$$

This means that even if a larger temperature difference than the nominal Δ*T* is used, the thermal output of the unit cannot be higher than its nominal power. There is no explicit bound on the mass flow rate, other than that following out of Equation (A6) and out of the mass flow rate limits of the pipes connected to the network node in which the conversion system is embedded. All conversion systems are assumed to be able to modulate their output 100%, unless stated differently.

#### Appendix A.2.1. Heat Pumps

The heat pumps considered in this study are large air source heat pumps. Its COP (Coefficient of Performance) is pre-calculated based on the variation of the ambient temperature, and based on the network temperatures:

$$\text{COP} = \frac{\text{Q}}{\text{W}} = \eta\_{\text{C}} \frac{T\_v}{T\_v - T\_a} \,\text{}\tag{A8}$$

assuming a non-ideal Carnot cycle with a relative efficiency *η<sup>C</sup>* = 0.6 (compared to the ideal Carnot cycle). In this equation, *W*˙ is the mechanical power supplied to the heat pump in the form of electricity, and *Tv* and *Ta* are the respective supply and ambient temperature. The Carnot efficiency, although seemingly optimistic, was chosen in line with data from the Danish Energy Agency [43] for large heat pumps.

#### Appendix A.2.2. Solar Thermal Collectors

The STC model assumes a steady-state heat delivery as a function of the mass flow rate and the solar irradiance only, and was derived from norm EN 12975-2 [44]:

$$\dot{Q}\_{\rm out}(t) = A \cdot \left( \eta\_0 \dot{Q}\_{\rm sal}(t) - a\_1 \left( T\_m - T\_a(t) \right) - a\_2 \left( T\_m - T\_a(t) \right)^2 \right) \tag{A9}$$

In this equation, *Q*˙ *out* and *Q*˙ *sol* represent the heat output and the solar irradiance on the unit surface respectively. The collector surface area is represented by *A*. *η*<sup>0</sup> is the base efficiency and *a*<sup>1</sup> and *a*<sup>2</sup> are temperature dependence parameters. Manufacturers of STC panels measure these values according to the norm mentioned before. This paper assumes Arcon Sunmark HT-SolarBoost 35/10 flat-plate collectors with an *η*<sup>0</sup> base efficiency of 0.839, an *a*<sup>1</sup> value of 2.46 W/(m<sup>2</sup> K) and an *a*<sup>2</sup> factor of 0.0197 W/(m<sup>2</sup> K2) [45]. *Ta* is the ambient temperature, whereas *Tm* represents the mean panel temperature. We assume *Tm* to be the average of the supply and return temperature of the network, which would correspond with a linear temperature increase along the collector pipe. When the heat output would become negative according to Equation (A9), it is set to be exactly 0W.

#### *Energies* **2019**, *12*, 2766

Solar panels are assumed to be oriented South at a 40◦ tilt angle. The solar panels studied in the design optimisation are either installed in an empty field in rows (no shading effects accounted for), or mounted on the south-oriented parts of the roof of the buildings in the selected neighbourhoods. The incident solar radiation on a tilted unit area was simulated in Dymola [46] using typical meteorological year data for Brussels. The effect of the incidence angle of solar radiation on the panel on the transmission and reflection of the solar irradiance is neglected for simplicity.

#### Appendix A.2.3. Geothermal Heating Plant

The Danish Energy Agency [43] considers a geothermal plant for district heating in combination with an electric heat pump that assists in extracting as much heat as possible from the ground. The water is pumped up from the ground extraction well and passes through a heat exchanger to preheat the water from the return side of the district heating system. The water in the geothermal loop is then cooled by the heat pump evaporator before it is pumped back into the injection well. The heat pump condenser in its turn injects heat into the district heating loop until to the desired supply temperature is reached.

A geothermal well is designed to be operated continuously and therefore it can hardly change the output power. We predetermine the period during which the geothermal plant is (in)active, usually in the summer months, in this case the plant is shut down from the end of May until the end of September (day 150 and 270 of the year).

The coefficient of performance of the system is determined based on the temperatures of the geothermal well doublet, combined with the design temperatures of the network and a pinch temperature difference of 5 K in the heat exchanger. This leads to a non-linear system of equations, which is solved as a pre-processing step and the resulting COP is used in modesto. We will not treat the non-linear system of equations in further detail here, as they can be derived from the system lay-out in Reference [43].

#### Appendix A.2.4. Industrial Excess Heat

Industrial excess heat is treated as a simple heat source with a fixed cost per unit of energy. It is assumed that the industry suppliers cover the entire investment and pumping costs and that these expenses are reflected in the excess heat price. Furthermore, we assume a constant availability of the nominal heating power, without the obligation to buy all of it. Down-periods for maintenance of the industrial processes are not accounted for.

#### Appendix A.2.5. Electricity Considerations

The electricity used by the heat pump components is assumed to be extracted from the electricity grid. We are using Belgian electricity grid data, namely the BELPEX day-ahead prices for the year 2014. The primary energy factor (PEF) and CO2 intensity for the Belgian grid have been calculated by Vandermeulen and Vandeplas [47]. We have calculated the hourly average profile as a representation and assume the same profile is valid for every day of the year. The time series of these profiles (CO2 intensity, PEF and electricity price) have been summarised in Figure A1. Note that the electricity price is incorporated in the representative days selection algorithm and hence the hourly average profile is not used in the optimisation framework. The reason for this choice is that the electricity cost has a large influence on the operational objective function of the optimisation, here the real yearly variation has been chosen instead of a daily aggregated version.

**Figure A1.** Time series and hourly average of PEF, CO2 intensity and electricity cost.

The fact that these time series are not realistically correlated with the weather data is disregarded here.

#### *Appendix A.3. Thermal Energy Storage*

This paper considers two types of TES, namely tank and pit thermal energy storage (TTES and PTES). Both are modelled in line with the model used by Vandewalle and D'haeseleer [26], namely a perfectly stratified model (using the network supply and return temperature as the respective high and low TES temperatures), including state-dependent heat losses to the surroundings. The state of charge must be between 0% and 100%, and the (dis)charge heat flow is unconstrained.

#### Appendix A.3.1. Tank TES Systems

TTES systems have the advantage of better insulation over PTES systems, but due to their construction they are limited in size. In addition, the heat losses can be reduced even further by choosing an advantageous tank aspect ratio, minimising the tank surface with respect to the volume. We assume that the tank is constructed with a 0.5 m thick concrete shell, surrounded by 0.3 m of insulating material [48]. The concrete has a thermal conductivity of 1.63 W/(m K) and foam glass gravel is chosen as insulation, with a thermal conductivity of 0.095 W/(m K). Also, we assume that the wall and insulation thickness are the same for all of the tank walls. For TTES, we simply assume a cylindrical shape. Hence, the dependence of the heat losses is perfectly linear in the SoC. To minimise the surface area of the cylinder, we choose the aspect ratio *h*/*D* = 1.

Seasonal TTES systems are generally partly buried (or bermed) underground, which leads us to assume the average of the ambient and ground temperature as the boundary condition for the side walls of the tank. The bottom only sees the undisturbed ground temperature as a boundary condition, whereas the top of the tank is exposed to the ambient temperature.

#### Appendix A.3.2. Pit TES Systems

A comprehensive description about design of PTES systems can be found in Sørensen et al. [49]. Sorknæs [50] suggests a pit side and bottom wall conductivity of 0.5 W/(m·K) in case of sand, with an equivalent soil layer thickness of 2 m before the soil reaches the undisturbed ground temperature. The top insulation has a thickness of 0.24 m with a thermal conductivity of 0.104 W/(m·K). Because of the width variation at different depths, the relation between the SoC and the heat losses is no longer perfectly linear. However, it can be assumed to be linear without too large deviations. This is the result of the trapezoidal cross-section of the pit, where the base is narrower than the top edge. In this paper, a fixed shape of the storage pit is assumed, where the top width is 9 times the height, and the inclination of the sides is exactly 1 in 2 (or 26.57◦), such that the bottom width becomes 5 times the height. Assuming a square top and bottom surface, the maximum considered pit volume of 200,000 m<sup>3</sup> has a height of approximately 16 m, whereas the top and bottom width are 144 m and 80 m respectively.

For the boundary temperature at the bottom of the pit, the undisturbed ground temperature *Tg* is used. The walls of the pit are assumed to be exposed to the average of the ground temperature and the outside air temperature, given the rather limited depth, including the part of the pit above the ground. The top cover is exposed to the ambient temperature only.

#### **Appendix B. Techno-Economic Data**

This section summarises data about costs, efficiency and lifetime of the used technologies.

#### *Appendix B.1. Thermal Energy Conversion Systems*

Table A1 summarises the used unit prices and investment analysis characteristics for the integrated optimal design and control algorithm. Note that cost figures per nominal power are always expressed considering the thermal power output.


**Table A1.** Characteristics of used conversion technologies.

#### *Appendix B.2. Thermal Energy Storage Systems*

The unit cost of large storage systems varies with their size, hence the used cost data is represented graphically in Figure A2. The data points are derived from Schmidt and Miedaner [53]. The economic lifetime of TES systems is assumed to be 20 y, the fixed maintenance cost is 0.70 %/y with respect to the initial investment cost [32]. The economic lifetime and maintenance cost are the same for PTES and TTES.

**Figure A2.** Cost evolution in function of volume for TES tanks and TES pits, using a linear interpolation between known investment costs for real systems. Data derived from Schmidt and Miedaner [53]. Note that the lower graph is representing the same data as the upper plot, but focussing on the lower volume range of the TES tank data.

#### *Appendix B.3. Thermal Network Pipes*

The available sizes of thermal network pipes and their dimensions are summarised in Tables A2 and A3. For the explanation of the different radii *r* and their respective accompanying diameters *d*, see Figure A3. The symmetric thermal resistance *Rs* is calculated using the derivation by Wallentén [54] and using the conclusions of van der Heijde et al. [39].

**Figure A3.** Schematic representation of the types of double pipes considered in this study with clarification of the dimensions used to calculate the thermal resistance of compound and twin pipes. The grey shaded area indicates insulation material, the red and blue shaded areas represent the water in the supply and return pipes, respectively.


**Table A2.** Dimensions of twin pipes and resulting symmetrical thermal resistance. All distances are expressed in m, the thermal resistance per unit length *R* in K m/W. All data has been retrieved from Isoplus catalogues [33].

*do* represents the outer diameter of the pipe wall, inside the insulation. *dcc* is the distance between the walls of the two pipes, such that *D* = *dcc*+*do* <sup>2</sup> . *s* is the pipe wall thickness.

**Table A3.** Dimensions of compound single pipes and resulting symmetrical thermal resistance. All distances are expressed in m, the thermal resistance per unit length *R* in K m/W. All data has been retrieved from Isoplus catalogues [33].


For compound pipes, *dcc* is the distance between the outer edge of the insulation jackets of the separate pipes, or *D* = *dc*+*dcc* <sup>2</sup> .

In the study of Ahlgren et al. [55], a distinction is made between inner and outer city areas. Netterberg et al. [56] seem to have consulted a similar data source. They found a linear investment cost regression:

$$I\_{inner} = 2.18 \cdot d\_{DN} + 308.46 \text{ [EUR/m]} \tag{A10}$$

$$d\_{outer} = 1.8596 \cdot d\_{DN} + 230.5 \text{ [EUR/m]} \text{ }, \tag{A11}$$

where *Iinner* and *Iouter* denote the respective investment cost per unit length for inner and outer city areas and *dDN* is the nominal diameter of the pipe in mm.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **An Approach to Study District Thermal Flexibility Using Generative Modeling from Existing Data**

**Camille Pajot 1,***∗***,†, Nils Artiges 1,2,***∗***,† , Benoit Delinchant 1,***<sup>∗</sup>* **, Simon Rouchier <sup>2</sup> and Frédéric Wurtz 1,***<sup>∗</sup>* **and Yves Maréchal <sup>1</sup>**


Received: 30 July 2019; Accepted: 18 September 2019; Published: 24 September 2019

**Abstract:** Energy planning at the neighborhood level is a major development axis for the energy transition. This scale allows the pooling of production and storage equipment, as well as new possibilities for demand-side management such as flexibility. To manage this growing complexity, one needs two tools. The first concerns modeling, allowing exhaustive simulation analyses of buildings and their energy systems. The second concerns optimization, making it possible to decide on the sizing or control of energy systems. In this article, we analyze, in the case of an existing residential neighborhood, the ability to study by modeling and optimization tools two scenarios of energy flexibility of indoor heating. We propose in particular a method allowing to rely on a varied set of data available to build the various models necessary for optimization tools or dynamic simulation. A study was conducted to identify the neighborhood's flexibility potential in minimizing *CO*<sup>2</sup> emissions, through shared physical storage, or storage in the building envelope. The results of this optimization study were then compared to their application to the virtual neighborhood by simulation.

**Keywords:** district scale; demand-side management; flexibility; MILP; *CO*<sup>2</sup> emissions; heat pump; ETL; data management

#### **1. Introduction**

#### *1.1. Energy Planning and Flexibility at the District Scale: Solutions and Issues*

To fight climate change, many energy transition policies are emerging around the world [1]. With the ambition to achieve a successful transition from fossil fuels to low-carbon production, the share of renewable energies into the energy mix increases. It is well known that buildings represent more than a third of global energy consumption, 40% of *CO*<sup>2</sup> emissions, and much more in urban areas [2]. Besides, the integration of diverse renewable energy sources in cities is a major step to achieve sustainability objectives [3]. This diversity of solutions increases the complexity of urban planning, both for design and retrofit, when one has economical and energy efficiency in mind.

To cope with the energy landscape complexity, several works of research led to software developments towards energy planning. Theses energy planning tools target different time scales (time step and range) and space scales (local to global). Among them, we can cite for example:

• **MODEST:** The MODEST Energy System Optimisation Model aims to compute how energy demand should be covered at the lowest possible cost, using a model of energy networks suitable

at regional and national scales [4]. MODEST uses linear programming (LP) to minimize capital and operation costs. The methodology uses a flexible time division to provide simulation results on both short and long time ranges.


All these models are great for testing and validating energy policies, energy landscape modifications at a wide scale, as well as studying medium- or long-term associated ecological and economic impacts. However, deep integration of intermittent renewable energies in the electrical network induces variability at the production side which could jeopardize the energy systems stability [12]. This phenomenon could be avoided by increasing the flexibility of consumption through demand-side management strategies, i.e., synchronizing the consumption with power production [13,14]. This area is more and more studied and especially applied to buildings whose consumption represents more than 55% of global electricity demand [2]. This raises a need for energy planning tools more suitable at a regional and medium scale (i.e., cities and districts). Many tools exist for this purpose. The reader can refer to the following reviews for an extensive overview: [15–17]. Among these, one can cite:


Generally, one can observe that such energy planning tools can be differentiated by the following criteria (see [25] for an extensive tools review):


Besides, two types of strategies are mainly investigated for demand-side management: electrical appliances that can be shifted [26,27] and thermal loads that can be modulated (with or without storage system) [28]. The principle of heating load modulation without any storage system consists of using the internal mass of buildings as energy storage. Thus, a building can be over-heated when consumption is needed and under-heated when production is lower. To describe this behavior, Panão et al. introduced the concept of Building as Battery (BaB) and illustrated it on residential buildings with photovoltaic panels [29]. Although many studies only focus on the BaB, some others address the challenge of modulating the heating load with Thermal Energy Storage (TES) [30]. The evaluation of the flexibility is often realized thanks to simulation results [28,31].

Therefore, to tackle thermal flexibility at a district scale, the tool to use must be characterized by the following:


From our current knowledge, only Oemof-Solph and the open-source tool OMEGAlpes (Optimization ModEls Generation As Linear Programming for Energy Systems) we are developing in our team seem to comply with such requirements.

OMEGAlpes is dedicated to the generation of linear optimization problems for energy systems [32]. It allows quickly building Mixed-Integer Linear Programs (MILP) to design and manage multi-carrier energy systems. OMEGAlpes models are based on energy flows and energy units allowing to quickly study numerous cases by setting and gathering elementary models. Big optimization problems (hundreds of decision variables) can be quickly solved at the district scale due to linear models.

Oemof is more oriented towards interfaces between complementary tools and is currently less complete than OMEGAlpes on the model side, which led us to pursue our developments on OMEGAlpes for our studies on thermal flexibility in districts.

A final issue in the process of a flexibility study is the good choice of (building) models and their parameter values. This aspect is a key point in the field of Urban District Energy Modeling (UBEM) and is discussed below.

#### *1.2. Objectives and Paper Structure*

In this study, we aim to present a methodology to study the flexibility potential of a district that can be obtained by the heating modulation. In this study case, heating systems are decoupled from the domestic hot water because of different temperature levels. Thus, only heating load modulation is addressed. We show:


The methodology is illustrated for a new residential district heated by a groundwater source. Located in Grenoble (France), the district is composed of 16 buildings outside-insulated with 11 floors on average. All buildings construction were initiated after 2010 and are designed according to the French energy policy for buildings (RT2012)—30% energy performance objective (30% more efficient than the RT2012, corresponding here to 50 kWh·m−2·year−<sup>1</sup> primary energy consumption). A simplified representation of the district is shown in Figure 1 and an overview in Figure 2. More details on geometric and physical parameters are given further in Section 2.3 Table 1.

**Figure 1.** Illustration of the district thermal fed by a groundwater source.

**Figure 2.** Grenoble "Cambridge Eco-district" overview.

The goal of the study case is to quantify the reduction of *CO*<sup>2</sup> emissions that can be obtained through flexibility on the district heating load. First, we present how one can use existing data to produce a suitable district model in OMEGAlpes. Then, we describe models used in OMEGAlpes and how one can use this tool to study two thermal flexibility scenarios. The heating load modulation thanks to thermal energy storage is addressed and then compared to the Building as a Battery (BaB) concept. Finally, results obtained by optimization are confronted with simulation results obtained with simulation models.

#### **2. Eco-District Modeling Based on Available Data**

The first step in our methodology is to build a dynamic thermal model of our district suitable for MILP optimization. In this section, we explain why proper data management is important for district model generation and how a data management tool can be used alongside UBEM tools. Then, we describe the data used and the models generated for OMEGAlpes.

#### *2.1. Data Variability and Heterogeneity*

Contemporary cities, and particularly since the development of the concept of "smart cities", expose more and more data for various users and applications. The data exposition is fostered by various actors, with corresponding privacy levels. For example, one can access city-related data through the following sources:


The main difficulties encountered by the engineer in obtaining city-scale data are as follows:


Since data accessibility is more relevant to organizational problems and policies, the scientific challenge in exploiting these data is more related to their variability and heterogeneity. Indeed, four main axes of heterogeneity are observed: quantity, granularity, structure and semantics (see Figure 3). Therefore, to use these data for modeling purposes, one needs to apply appropriate tools and techniques to handle this heterogeneity.

The problem of data management is commonly encountered in the IT industry. In numerous application fields (online sales, social networks, advertising, etc.), developers have to handle various and sometimes unreliable data, encoded in different formats and databases. Generally, data are stored in different and specific databases (data warehouses) and not exploitable *as is*. One must then develop a middleware layer to extract the data, transform them and load them to client databases and interfaces to comply with clients' needs.

Our problem here is quite similar: we want to extract district-related data from various origins and encoded on various files, pre-treat these data and store them in models able to perform simulation and/or optimization tasks. Consequently, this approach of Extraction–Transformation–Loading (ETL) seems well adapted to the problem of district modeling from existing data. The reader can find general information about ETL techniques in [37].

**Figure 3.** Data heterogeneity in districts.

#### *2.2. District Buildings Modeling*

Detailed dynamic thermal modeling of districts is more encountered in UBEM simulators tools such as CitySim [38], City Energy Analyst [39], TEASER [40] or CityBES [41], than in energy planning tools. The reader can refer to [42] for an extensive review. Considering the complexity to model a whole district, the amount of data potentially required and the nature of available data, some simplifying strategies are often used:


Among UBEM tools, TEASER is particularly adapted to the generation of low order models with few input data. This Python package developed at the University of Aachen can generate a simple "archetype" model of a building with a minimum of five parameters and can involve statistical databases for data enrichment. The generated models are Python objects translated in Modelica using IBPSA annex 60 or Aixlib libraries [40].

We build here "four walls elements (i.e., interior walls, exterior walls, floor plate, and roof) SingleFamilyDwelling" archetypes models using TEASER according to the IWU (Institut Wohnen und Umwelt—Institute for Housing and Environment) topology issued from the EPISCOPE project [36]. "SingleFamilyDwelling" corresponds to the archetype's data enrichment method. At the moment, TEASER mostly supports archetypes issued from studies of the German stock. This is not an issue here

since we only test our methodology, but the implementation of French archetypes should be necessary for more accurate results. The resulting model is a RC reduced-order thermal model corresponding to the "AixLib.ThermalZones.ReducedOrder.RC.FourElements" component of the AixLib Modelica Library [43], as depicted in Figure 4. The corresponding Modelica simulation is further used as a reference for the virtual district.

**Figure 4.** Four elements building model generated by TEASER [43].

A second model is required to perform optimizations. Thus, we developed a simpler linear building model in OMEGAlpes. Furthermore, too many model parameters can be counterproductive for an early-stage study. Therefore, we implemented a simplified RC-model of the Swiss SIA2044 norm in OMEGAlpes meeting our main needs: it can be easily built with few data, all the equations are linear (described Figure 5), and the Swiss building structure is very similar to the French one.

**Figure 5.** RC model according to SIA 2044.

This RC-model is composed of a thermal capacity C*m* and five thermal resistances:


*Energies* **2019**, *12*, 3632

Parameters in TEASER models are translated to the OMEGAlpes model such that global thermal transfer coefficient (U) values and thermal capacitance are preserved.

#### *2.3. Generation of the Building Stock Model from Existing Data*

In our residential district study case, the following data are available:


For OMEGAlpes building models, we need the following data:


All the data required for OMEGAlpes building models can be deduced from generated TEASER models. Then, the UBEM generation is processed according with the workflow summarized in Figure 6.

**Figure 6.** Workflow for district model generation.

In this workflow, one parses all the data files first. For buildings with RSET files, all the required parameters to build TEASER and OMEGAlpes models, except for emissivity, absorbtivity, and transmittance, are present. For other buildings, there are only enough data to generate TEASER archetypes. For some of them, the floor area is not directly available and one has to find corresponding polygons in the land registry file to estimate them. To complete missing data for OMEGAlpes models, one can extract parameters generated by TEASER in archetypes. For each parameter parsing, injection or extraction, one has to deal with different formulations and units. The generated dataset used to build all OMEGAlpes building models is summarized in Table 1.


**Table 1.** District main input parameters statistics for OMEGAlpes.

To apply this workflow, we developed a specific Python package to ease file parsing, data manipulation (with an intermediate SQLite database) and district model generation, with modularity in mind (for further data integration). The architecture of this tool is summarized in Figure 7.

**Figure 7.** Python software architecture for data handling and model generation.

Such an approach is close to the ETL methodology. As ETL processes are well suitable for UML modeling [44], the choice of an Object-Oriented Programming language such as Python is appropriate. Besides, the support of Python by the scientific community eases the development of the "Transform" part.

#### **3. Modeling of the Optimal Planning Problem**

#### *3.1. Optimal Planning of the District Heating Systems with OMEGAlpes*

As already mentioned, two flexibility approaches to manage the district heating systems are addressed through the study case:

• The first one consists in designing an energy system composed of a heat pump and thermal energy storage to minimize the *CO*<sup>2</sup> emissions of a fixed district heating load.

• The second study case also aims to minimize the *CO*<sup>2</sup> emissions of the buildings' heating load, but thanks to flexibility through building envelopes. In this case, specific building models dedicated to the optimization should be used to estimate how the load can be modulated.

In both studies, we aim to estimate the possibility to decrease the *CO*<sup>2</sup> emissions by designing and operating the system. The studies were conducted during two weeks in January, which usually represent the coldest period and are critical for the power system. Thus, the design of the system can be significant for the entire heating period. It is important to notice here that our goal is not to predict energy needs and *CO*<sup>2</sup> emissions for an entire year, but to be closer to the operation. Therefore, focusing on two weeks allows us to anticipate the possibility to pursue our work with a model predictive control approach thereafter.

To define the OMEGAlpes optimization models, a graphical formalism was defined to represent the energy units and power flows (see Figure 8).

**Figure 8.** OMEGAlpes formalism for energy system modeling.

Let us introduce Study Cases 1 and 2; the results are detailed in Section 4.

#### *3.2. Study Case 1: Flexibility through Thermal Energy Storage (TES)*

The first study case deals with energy flexibility provided by a Thermal Energy Storage (TES) to minimize the *CO*<sup>2</sup> emissions of the district heating load. The energy system studied is composed of the district heated by geothermal groundwater through a heat pump and thermal energy storage to provide demand-side management. The goal of this study case is to design the whole supplying system (heat pump and storage). To do so, we used three OMEGAlpes units to model the energy system: the district heating load, the heat pump, and the thermal energy storage, as shown in Figure 9.

**Figure 9.** Modeling of the energy system of the first study case (heat pump, district heating load and thermal energy storage) according to OMEGAlpes formalism.

#### 3.2.1. Estimation of the District Heating Load

In this first study case, only the flexibility provided by the storage energy system is addressed so that the district heating load cannot be modulated and is thus an input of our optimization problem. To estimate the thermal needs of the district, we relied on results from a first optimization obtained

with OMEGAlpes which can be considered as a temperature regulation simulation. All buildings were modeled as described in the previous section and set with standard occupancy schedules obtained by TEASER and a temperature set-point of 20 ◦C. The objective of the optimization is to minimize the sum of the over-heating and the result (see Figure 10) is taken as the dynamic thermal consumption of the district P*dist*(t). In this figure, we can notice that, during the days, the district heating load is very low. This could be explained by the high insulation of the buildings which require low consumption and thus can benefit from occupancy and solar gains to cover their needs.

**Figure 10.** District heating load during a two-week period in January obtained by optimization.

#### 3.2.2. Modeling of Heat Pump

Composed of new residential buildings, the district can be heated by low temperatures around 35 ◦C. With a groundwater temperature around 15 ◦C, this specificity allows the heat pump to reach a high Coefficient Of Performance (COP) of 5. Moreover, the temperature of the groundwater is assumed to be invariant so that we can consider the COP to be a constant equal to 5. Therefore, the heat pump is modeled by the relations between the thermal power provided by the groundwater *Pthermin* (*t*), the electrical power consumed by the heat pump *Pelec*(*t*), the thermal power delivered *Pthermout*(*t*) and the COP, as described in Equation (1).

$$\begin{cases} P\_{therm\_{in}}(t) + P\_{elcc}(t) = P\_{therm\_{out}}(t) \\ P\_{therm\_{out}}(t) = \text{COP} \* P\_{elcc}(t) \end{cases} \tag{1}$$

In this study case, a trade-off was chosen between different levels of accuracy of the whole energy system modeling according to the uncertainties relating to the occupants' behaviors. Indeed, as we aim to estimate orders of magnitude of the *CO*<sup>2</sup> emissions reduction obtained by heating flexibility, the modeling of the heating systems is very simplified. For further studies, a deeper level of modeling could be needed to provide a more accurate estimation.

3.2.3. Modeling of Thermal Energy Storage (TES)

Multiple types of thermal energy storage systems are used in the literature to smooth building thermal needs. However, the most widespread technology used remains water tanks for their simplicity and low costs.

The power stored to the TES *Pstor*(*t*) is defined as the difference between the charging power *Pc*(*t*) and the discharging power *Pd*(*t*) as described by Equation (2).

$$\begin{cases} P\_{\text{stor}}(t) = P\_{\text{c}}(t) - P\_{\text{d}}(t) \\ P\_{\text{c}}(t) \ge 0 \\ P\_{\text{d}}(t) \ge 0 \end{cases} \tag{2}$$

Moreover, the relation between the energy contained in the water tank *e*(*t*) and the charging and discharging powers is defined by Equation (3). The storage capacity *Cstor* is defined as the maximal value of *e*(*t*).

$$\begin{cases} e(t+dt) = e(t) \* (1 - a\_{sd}) + (P\_c \* \eta\_c - \frac{P\_d}{\eta\_d}) \* dt \\ e(t\_0) = e(t\_f) \\ e(t) \le \mathbb{C}\_{stor} \end{cases} \tag{3}$$

where


In this study, a stratified storage system is considered called thermocline storage whose management is more complex than traditional storage (more details can be found in [45]). Indeed, in our case, we assumed that the storage has to be fully charged at least once per five days to optimally operate. The first step to model this constraint is to define a variable to indicate if the storage is fully charged. To do so, a binary variable was introduced: *is\_soc\_max(t)* which equals 0 when the state of charge is lower than 100% and 1 when the storage is fully charged. The definition of this indicator was realized thanks to Equation (4), where *Cstor* is the storage capacity, *estor*(*t*) is the energy contained in the storage at the time t and is taken equal to 10<sup>−</sup>3.

$$\begin{cases} \mathbb{C}\_{\text{star}} \* (1 + \text{is\\_soc\\_max}(t) - \epsilon) \ge \mathfrak{e}\_{\text{star}}(t) \\ \mathbb{C}\_{\text{star}} \* \text{is\\_soc\\_max}(t) \le \mathfrak{e}\_{\text{star}}(t) \end{cases} \tag{4}$$

Then, our constraint can be expressed thanks to a sliding window including five days. Let *tcycl* be the time step corresponding to the end of the first five-day period; the constraint of at least one full charge during five days is defined by Equation (5).

$$\forall t \in \left[t\_{\text{cycl}}; t\_f\right], \sum\_{k=t-t\_{\text{cycl}}}^{t} \text{is\\_soc\\_max}(k) \le 1\tag{5}$$

#### 3.2.4. Modeling of *CO*<sup>2</sup> Emissions of the District Heating Load

In this study case, the *CO*<sup>2</sup> emissions of the district heating load (*EmCO*<sup>2</sup> ) come from the electrical consumption of the heat pump. Fed by the French power system, the heat pump emissions vary dynamically according to the French grid *CO*<sup>2</sup> emissions rate (*emCO*2,*rate*(*t*), see Figure 11).

Thus, the *CO*<sup>2</sup> emissions of the district heating load can be calculated by Equation (6), so that changing the heat pump operation could lead to *CO*<sup>2</sup> reduction, which we tried to achieve thanks to thermal energy storage in this study case.

$$\mathcal{E}m\_{\text{CO}\_2} = \sum\_{t=t\_0}^{t=t\_f} P\_{\text{clcc}}(t) \* \mathcal{E}m\_{\text{CO}\_2, rate}(t) \tag{6}$$

**Figure 11.** *CO*<sup>2</sup> emissions rate of the French power system during a two-week period in January 2018.

#### 3.2.5. Energy System Design Parameters

As explained above, the objective is to minimize the *CO*<sup>2</sup> emissions of the district heating load. To do so, we considered a groundwater source heat pump coupled with the thermal energy storage that we aim to design. Three parameters are optimized:


#### *3.3. Study Case 2: Flexibility through Heating Loads Modulation (BaB)*

In this second study case, thermal flexibility is provided by the building envelopes. Each building is modeled individually so that the district load can be deduced by the addition of each building heating load. Thus, the thermal load of the district can be directly modulated without any external thermal energy storage (see Figure 12).

**Figure 12.** Modeling of the energy system of the second study case (16 buildings and 16 heat pumps) according to OMEGAlpes formalism.

#### 3.3.1. Estimation of the District Heating Load

To guarantee the occupants' thermal comfort, the operative temperature is constrained to be higher or equal to 20 ◦C. Thus, the building can be over-heated by moments to store heat into the buildings while keeping thermal comfort.

The heating load can be calculated for each time step according to the thermal RC model available in OMEGAlpes and presented in the previous section. Besides constraining the operative temperature, the boundaries conditions are the same as before. These internal gains (from occupancy and weather) are applied to the nodes a, c and m (see Figure 5). More details about the model can be found in [46].

3.3.2. Modeling of Heat Pumps and *CO*<sup>2</sup> Emissions of the District Heating Load

The configuration of the district is slightly changed since each building is fed by its heat pump. Each heat pump is designed with an over-sizing (around +66%) according to the reference heating need of being able to use flexibility. The total maximal electrical consumption allowed to feed all the heat pumps was set to 500 kW.

The modeling of the *CO*<sup>2</sup> emissions is similar to the previous study so that each heat pump emits according to its electrical consumption. However, the objective to minimize the *CO*<sup>2</sup> emissions is global.

In this study case, the energy system is designed before running the optimization. Therefore, the minimization of the *CO*<sup>2</sup> emissions is based on finding an optimal operation of all the heat pumps of the districts.

#### **4. Results**

This section is divided into two main subsections:


#### *4.1. Optimization*

The study cases presented in this paper are realized for a time step of 10 min for two weeks in January. For the first one, each optimization problem generated is composed of 38k variables (28k continuous and 10k binaries) for 61k constraints. The resolution was launched on an Intel bicore i5 2.4 GHz CPU with the Gurobi solver so that the optimization problem was solved within less than 10 s on average for 192 optimizations. The corresponding results are detailed Section 4.1.2.

The second study case consists of a single resolution since only one configuration is studied. The associated optimization problem is composed of 1211k variables (1100k continuous and 111k binaries) for 1263k constraints. The resolution was launched on the same Intel bicore i5 2.4 GHz CPU with the Gurobi solver and the optimization problem was solved within 23 min. The dynamic results are detailed Section 4.1.3.

#### 4.1.1. Flexibility Potential

To evaluate the gains obtained by optimization, the first step is to estimate the maximal reduction in *CO*<sup>2</sup> emissions that can be achieved. A simple way to evaluate this maximum is to allow shifting each 10-min power slot of the load curve to minimize the *CO*<sup>2</sup> while consuming the same energy during the two weeks. In this case, the *CO*<sup>2</sup> emissions can be reduced to a maximum of 22% keeping the current heat pump maximal power consumption (300 kW). Of course, in this case, the building comfort (internal temperature) is not guaranteed. This potential of 22% *CO*<sup>2</sup> savings is used to compare the next results.

The dynamic result of this naive estimation is shown Figure 13. We can notice that some *CO*<sup>2</sup> emissions reduction can be obtained by anticipating or removing the consumption for a few hours (short-term flexibility), but the longer-term variation of the *CO*<sup>2</sup> levels lead to a need for longer-term flexibility.

**Figure 13.** Theoretical maximal flexibility to minimize *CO*<sup>2</sup> emissions by shifting 10-min values from the reference profile.

#### 4.1.2. Study Case 1: Flexibility through Thermal Energy Storage (TES)

In this study case, three elements were designed: the storage capacity (*Cstor*), the storage self-discharge coefficient (*αsd*) and the maximal electrical power consumed by the heat pump (*Pmax elec* ). Results are shown in Figure 14 for the three self-discharge coefficients studied (0.125%, 0.25% and 0.5%). The *CO*<sup>2</sup> emissions reduction obtained in each configuration is drawn according to the storage capacity and the maximal electrical power consumed by the heat pump.

Regarding the storage capacity, we can notice that 100 kWh is too small to reduce the *CO*<sup>2</sup> emissions regardless of the two other parameters. For larger capacities (≥1 MWh), the impact on the *CO*<sup>2</sup> emissions reduction begins to be noticeable and is strongly correlated with the self-discharge coefficient and the maximal power consumed by the heat pump.

For a storage capacity under 2 MWh, we can notice that the reduction in *CO*<sup>2</sup> emissions are lower than 3.5% for all designs. However, in the case of a TES of 48 MWh with 0.125% of self-discharge and a maximal electrical power of the heat pump of 2500 MW, we manage to reach 20% reduction in *CO*<sup>2</sup> emissions of the district heating load. Knowing that the average daily heating consumption of the district during the period is 8 MWh, a 4 MWh storage corresponds to 12 h of consumption while a capacity of 48 MWh corresponds to six days.

**Figure 14.** *CO*<sup>2</sup> emissions reduction according to the design of the heat pump and energy storage system.

Two phenomena happen when designing the energy system to use flexibility to reduce *CO*<sup>2</sup> emissions:


Although it seems reachable to have a strong impact on the *CO*<sup>2</sup> emissions of the district heating load with a big TES and a heat pump with high electrical power needs, this design choice leads to other problems. Indeed, choosing the kind of heat pumps means to increase the electrical power peaks and could report *CO*<sup>2</sup> emissions decreases from the heating side to increases at the electrical one. Over-sizing the heat pump should thus be carefully considered taking this effect in mind. Moreover, a 48 MWh water tank is expensive and takes a lot of space so that it is not an ideal solution. Nevertheless, it could be very interesting to deeply consider the level of insulation that can have an important impact.

Finally, using a building's envelope as storage should be investigated.

#### 4.1.3. Study Case 2: Flexibility through Heating Loads Modulation

In this case, the *CO*<sup>2</sup> reduction is low (0.5%). Although the flexibility potential was previously estimated to 22%, it mainly relies on medium- and long-term flexibility. However, with the Building as Battery (BaB) concept, the flexibility addressed in our case can be defined as short term. Indeed, we can notice on the optimization results (Figure 15) that no energy is shifted for more than one day. High consumption peaks allow profiting from low-*CO*<sup>2</sup> rate periods but the energy cannot be stored in the long run. Indeed, with external wall insulation systems (EWIS), buildings envelopes form a relatively small storage capacity, while the *CO*<sup>2</sup> variability is more long-term in this case.

Moreover, over-heating the buildings leads to an increase in district energy consumption (0.9%), so that the environmental gain due to shifting the heating load is reduced. However, the mean operative temperature goes from 20.3 ◦C to 20.4 ◦C, i.e., an increase of 0.1 ◦C (0.4%). Therefore, the increase of the consumption induces a better thermal comfort while reducing *CO*<sup>2</sup> emissions. Many studies achieve a greater reduction by allowing over- and under-heating [47], i.e., by using both energy flexibility and sobriety, while we choose to focus on the impact of flexibility only.

**Figure 15.** Optimal flexibility to minimize *CO*<sup>2</sup> emissions by allowing over-heating.

With a maximal electrical consumption of heat pumps of 500 kW, this scenario can be compared to those with a single 500 kW heat pump. The reduction of *CO*<sup>2</sup> emissions obtained by buildings as storage is similar to results with a TES with a capacity between 100 kWh and 1 MWh, whatever the self-discharge coefficient.

Finally, the reference scenario is the result of an optimization problem by providing minimal energy to maintain thermal comfort. Thus, when over-heating is not compensated by the *CO*<sup>2</sup> diminution, the *CO*<sup>2</sup> emissions minimization corresponds to the energy minimization.

For this reason, we need to compare a simulation reference profile to the flexible scenario to see if the improvement is preserved with a standard controller. Besides, comparing OMEGAlpes results with simulation can lead to investigating optimal control robustness during application.

#### *4.2. Simulation*

Experimenting flexibility scenarios computed by OMEGAlpes directly on the real district is a complex task. Before performing these tests, we started simulation studies to validate our approach and/or identify the main issues towards implementation. Since we used TEASER as an intermediate for building OMEGAlpes models, Modelica models are also ready to simulate for each building. First, we can compare heating needs and operative temperature profiles (defined as a mean between air and radiant temperatures) for both modeling approaches. To do so, we used two specific control scenarios:


computed power shifts on Modelica models, we applied the operative temperatures computed for the flexibility scenario as a new setp oint profile.

In the first case of constant temperature set point, we obtained the results presented in Figure 16 (buildings mean operative temperatures and the sum of all buildings consumption).

Constant setpoint temperature scenario

**Figure 16.** Comparison between Modelica and OMEGAlpes—Constant temperature set point scenario.

The first obvious difference stands in energy needs. Modelica models generated by TEASER have a more important consumption for the same comfort criteria and are therefore probably less insulated. Besides, temperature peaks are shifted between models. These shifts can be due to the differences in control strategies, and/or in computed internal gains despite using the same scenarios and weather files. Consequently, even with the same data sources, it appears hard to obtain identical dynamic behaviors between different modelers. Further effort must be invested to preserve global building characteristics during model translating (in our case, model simplification). Besides, this also suggests the use of a model calibration phase before implementing any model-based control strategy.

If we consider the application of the flexibility scenario in Figure 17, the consumption differences between both models are visually less important, except for higher spikes for Modelica models, but dynamics of operative temperature are sill very present (more inertia to go down for Modelica models).

We also compared performance results between constant temperature set point and computed flexibility temperatures on Modelica models only, to see if it also leads to improvement despite model discrepancies (see Figure 18).

Here, the dynamics induced by the flexibility scenario are very noticeable. As for OMEGAlpes results, we observe power shifts and spikes inducing heat storage in buildings envelopes. Unfortunately, performance is not preserved, since both energy consumption and *CO*<sup>2</sup> emissions are worsened (Table 2).

This implies that the flexibility command computed here is not robust to the model discrepancies we are facing. Therefore, the robustness of flexibility scenarios towards modeling uncertainties is certainly a key research topic before real-life integration if we do not want optimized scenarios to be counterproductive. This is also true during the early stage of design where low order models are used to investigate energy scenarios.

Flexibility scenario - OMEGAlpes vs Modelica

**Figure 17.** Comparison between Modelica and OMEGAlpes—flexibility scenario.

Modelica simulations - Constant setpoint VS Flexibility scenarios

**Figure 18.** Modelica simulations—comparison between control scenarios.

**Table 2.** Performance indices with flexibility scenario to optimize *CO*<sup>2</sup> emissions against constant temperature setpoint.


#### **5. Conclusions**

Based on an ETL (Extract–Transform–Load) method, we have initiated a tool based on the heterogeneous available data of buildings at the district scale, which can generate the necessary data for optimization and simulation models. In particular, we have applied this method to the case of data available on a new residential district composed of 16 residential buildings. This work makes it possible to identify the important parameters for different modeling tools at the neighborhood scale, to extract them from available data, or to estimate them when they are not available.

We then carried out two flexibility studies, based on the OMEGAlpes tool, which requires modeling in MILP formulation. The first study analyzed the design of a heat pump (nominal power) and storage (capacity and self-discharge factor) to desynchronize the production of heat and use to heat buildings. It appears that a large investment is necessary to try to reach the maximum potential (estimated at 22%), which relies in particular on long-term flexibility (more than a week). The second study relied on thermal storage via the building envelope. This zero investment solution is, therefore, a potential alternative to the previous case. However, the results obtained below 1% show that the low storage capacity of these residential buildings does not allow addressing flexibility considering a *CO*<sup>2</sup> variability during several days.

Finally, the tool we developed for data processing at the neighborhood scale allowed us to easily set up a validation process. Thus, we have transmitted the flexibility results in the simulation model using the Modelica AixLib library. The results show a predicted performance degradation compared with optimization results. On very small gains (<1%) obtained by the upward flexibility (temperature > 20 ◦C), it even presents negative performance in energy and *CO*2. This lack of robustness to modeling assumptions reinforces the idea that a tool to generate different levels of modeling based on available data will be indispensable for future studies related to robustness in optimization.

**Author Contributions:** Supervision, B.D., S.R., F.W. and Y.M.; Writing—original draft, C.P. and N.A.; and Writing—review and editing, C.P., N.A. and B.D.; C.P. and N.A. contributed equally to this work.

**Funding:** This work was partially supported by the ANR project ANR-15-IDEX-02.

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

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Combined Optimal Design and Control of Hybrid Thermal-Electrical Distribution Grids Using Co-Simulation**

**Edmund Widl 1,\* , Benedikt Leitner <sup>1</sup> , Daniele Basciotti 1, Sawsan Henein 1, Tarik Ferhatbegovic <sup>1</sup> and René Hofmann 1,2**


Received: 28 February 2020; Accepted: 8 April 2020; Published: 15 April 2020

**Abstract:** Innovations in today's energy grids are mainly driven by the need to reduce carbon emissions and the necessary integration of decentralized renewable energy sources. In this context, a transition towards hybrid distribution systems, which effectively couple thermal and electrical networks, promises to exploit hitherto unused synergies for increasing efficiency and flexibility. However, this transition poses practical challenges, starting already in the design phase where established design optimization approaches struggle to capture the technical details of control and operation of such systems. This work addresses these obstacles by introducing a design approach that enables the analysis and optimization of hybrid thermal-electrical distribution systems with explicit consideration of control. Based on a set of key prerequisites and modeling requirements, co-simulation is identified as the most appropriate method to facilitate the detailed analysis of such systems. Furthermore, a methodology is presented that links the design process with the implementation of different operational strategies. The approach is then successfully applied to two real-world applications, proving its suitability for design optimization under realistic conditions. This provides a significant extension of established tools for the design optimization of multi-energy systems.

**Keywords:** design optimization; control and operation; multi-carrier energy systems; co-simulation

#### **1. Introduction**

The joint design and integrated operation of electrical distribution grids and district heating systems promises to exploit hitherto unused synergies for increasing efficiency and flexibility. The envisaged goal is to achieve an increase of the hosting capability of electrical distribution networks for renewable energy sources (RES), while simultaneously reducing greenhouse gas emissions and primary energy use of district heating systems. For instance, in electrical distribution grids, the integration of photovoltaic (PV) and wind generation into the existing infrastructure has severe consequences on the power quality. At the same time, district heating networks struggle to replace carbon intensive heat plants with economically feasible combined heat and power plants and other renewable heat sources. In case network planning and operation are done appropriately, the diverse storage technologies, the different time constants and the diverse constraints regarding demand and generation can complement each other.

However, there is a lack of tools and methods for designing such integrated energy systems, especially in view of a detailed validation of proposed control and operation schemes. The work

presented here introduces a design approach that addresses the technical challenges of both domains at the same time, including their dynamic interaction at network level as well as local and high-level control (see Figure 1). This enables the exploitation of synergies in the operation of hybrid thermal-electrical distribution system in an optimal way. The presented approach relies on co-simulation, which enables the coupling of domain-specific simulators and multi-purpose tools in a way that allows to combine multi-physics simulations and optimization procedures. This approach is motivated by the criteria established in Section 2, and backed up by reviewing the available tools and methods according to the state-of-the-art for simulation and optimization (Section 3). A description of the modeling approach and the utilized simulation tools is presented in Section 4. Based on this, a novel methodology is presented in Section 5 that links design constraints to suitable operational strategies and optimization methods. Finally, the applicability of the presented approach is demonstrated in two real-world applications in Section 6.

**Figure 1.** Scope of the design approach proposed in this work.

#### *Scope and Main Contributions*

This paper aims at developing an integrated optimal design and control framework for coupled district heating and electrical distribution networks, extending the scope of traditional design tools for multi-energy systems. The simulation and optimization framework is illustrated for designing storage and thermal-electric appliances in two case studies, i.e., an industrial and a rural area. Note that the focus of this work lies on the methodological contribution rather than on the case study results. Compared to existing studies, this work explicitly includes detailed models for the district heating network, the electric distribution grid as well as low- and high-level control implementations already in the design stage. Thus, the method considers the impact of different levels of control and operation on the optimal system design. The coupling of heterogeneous modeling paradigms and tools is established via a co-simulation approach and the design optimization relies on the use of meta-heuristics.

This work addresses experts with backgrounds in district heating, electric distribution, energy storage deployment, control theory, design optimization and (co-)simulation alike. The multi-disciplinary nature of the proposed design approach (see also Figure 1) requires a comprehensive introduction and presentation of the work to make it understandable and useful for readers and experts coming from these different fields.

In summary, this paper contributes to the research field by presenting a simulation-based design optimization approach that: (i) is based on a fully dynamic thermal-hydraulic district heating and electrical distribution network model; and (ii) explicitly includes closed-loop control implementations and, thus, leads to the optimal design being dependent on operational aspects. The method is used in two case studies that exhibit different control complexities, i.e., model predictive and rule-based, and design spaces, i.e., a finite set of allowed solutions and an infinite set of possible designs.

#### **2. Prerequisites for the Design Optimization of Hybrid Distribution Systems**

#### *2.1. Determining Factors of Hybrid Distribution Systems*

Hybrid networks are realized through the physical interconnection and joint operation of electrical and thermal distribution grids. From a technical perspective, this is accomplished with the help of *coupling points*, i.e., devices which directly or indirectly enable the exchange of energy across carrier domains. For the proposed design approach, a (preliminary) *technical system layout* for both the electrical network and the thermal network topology is required from the domain experts, which already includes the type and position of the coupling point(s). The considered *degrees of freedom* in this work are typically related to the sizing of components (e.g., storage capacities or power ratings) or controller set-points (e.g., gains or thresholds). A specific technical system layout together with a specific set of numerical values for theses degrees of freedom is referred to as *system configuration* in this work.

The proposed design optimization approach explicitly considers operational and control aspects at different levels. In general, there are control systems at the process level (e.g., for heat pumps or transformers) that are designed to ensure that objectives are achieved locally (e.g., valve positions and tap changer). However, the complexity of control schemes increases drastically when individual processes are combined to larger systems and new (common) control/optimization targets are defined. In such a case, a higher-level control instance—referred to as *operational strategy* in this work—is required that governs the local processes in compliance with relevant system constraints (compare with Figure 2).

**Figure 2.** Overview of local control (process level) versus operational strategy (system level).

#### *2.2. Key Prerequisites and Modeling Requirements*

This work focuses on the optimization of technical design aspects, assessing the impact of system configurations and operational strategies on the system's performance. This requires not only the quantification of performance indicators on the system level, but also validation down to the component level, in order to avoid infeasible results. Hence, any valid design approach for hybrid thermal-electrical distribution networks has to fulfill the *key prerequisites* (KP) presented in Table 1.

In general, the most viable approaches in this context are simulation-based methods, as they allow the qualitative and quantitative examination of new concepts in a comparably fast and inexpensive way before deployment. However, the accurate and detailed modeling of all involved networks including their coupling points and control is of utmost importance. When deployed in a joint simulation, the (sub-)models representing the individual domains have to comply with the *modeling requirements* (MR) presented in Table 1.

It should be noted that some of the criteria defined above do not necessarily apply for design optimization approaches focusing on socioeconomic goals, especially regarding the required level of detail (spatial and temporal).

**Table 1.** Key prerequisites and modeling requirements for the assessment of hybrid networks.


#### **3. State of the Art Regarding the Simulation and Optimization of Hybrid Distribution Systems**

#### *3.1. Simulation of Coupled Heat and Power Networks*

Previous work on multi-carrier energy systems was mainly focused on determining an optimal mix of energy sources as well as conversion and storage technologies. Comprehensive overviews were given by Mancarella et al. [1] and van Beuzekom et al. [2], who independently concluded that existing approaches are not suited for detailed technical assessments that are required for network infrastructure planning and operation. Consequently, they argued that there is a need for new tools and methods in this regard. One of the main reasons is that established tools for system design, such as EnergyPLAN [3] or HOMER [4], do not offer models that are detailed enough to evaluate the effects of local controls on the process level, thus disagreeing with KP1, MR2 and MR3. Unfortunately, simulation tools that focus on the technical evaluation of energy systems on the process level are by themselves also not suited in this context [5], as they are typically concerned with just a single energy-related domain. Anything outside their direct focus is usually taken into account only implicitly or using simplified models, thus disagreeing with KP1, MR1, MR2 and MR3.

A potential solution is provided by multi-domain modeling languages. For instance, there exist several libraries for the Modelica [6] language that target energy-related domains, such as power systems [7,8] or buildings [9], which in combination allow the modeling of hybrid thermal-electrical distribution systems (and thus satisfying KP1 and MR1). It has also been demonstrated how such models can be utilized for standard optimization approaches (satisfying KP2), see for instance [10]. However, domain experts (thermal, electrical, controls) are often not trained to use such tools, which is complicating compliance with KP3. Furthermore, since Modelica focuses on modeling physical systems, the implementation of controls and operational strategies with a high algorithmic effort, e.g., model predictive control, can become difficult, thus complicating compliance with MR3. Most of the above arguments are also valid for similar languages and tools, e.g., MATLAB/Simulink/Simscape [11]. Hence, in the context of simulating hybrid thermal-electrical distribution networks, the usability of approaches relying exclusively on multi-domain modeling languages is in practice (still) limited, or at least the associated effort is (still) considerably high.

An alternative solution is offered by co-simulation approaches, which overcome these limitations by enabling the coupling of domain-specific simulators and multi-purpose tools (thus satisfying MR2 and MR3). This allows domain experts to use the most appropriate tools according to the state of the art for their respective domain, including advanced optimal control schemes. On the one hand, this guarantees an adequate and precise representation of the individual domains (thus satisfying KP1 and MR1) [12]. On the other hand, it facilitates the participation of and interaction between experts from different domains, whose expertise is often closely linked to specific tools (thus satisfying KP3). This advantage has led to various research activities in many energy-related domains, e.g., buildings [13,14], power systems [15–18] and hybrid distribution networks [12,19].

Given these considerations, co-simulation has been chosen to analyze systems according to the criteria defined in Section 2 for the work presented here. Even though other approaches—and especially approaches based on multi-domain modeling languages—are expected to become mature and flexible enough in the future, the advantages of coupling domain-specific simulators and multi-purpose tools in a co-simulation still prevail.

Within this context, established tools for the design and optimization of multi-energy systems can be considered as important guides for an overall design process, as their results should be used as starting point for a detailed technical evaluation.

#### *3.2. Design Optimization of Coupled Heat and Power Networks*

Available literature on design optimization of hybrid thermal-electric networks is mainly based on the energy hub concept introduced by Geidl et al. [20]. Related work focuses on determining optimal design of energy hubs such as selection and sizing of coupling units and storages [21]. However, applied models rely on many simplifications and are not able to cover technical aspects relevant in power networks (e.g., voltages and reactive power), in district heating (e.g., temperature propagation and pressures), and also in individual components (e.g., temperature stratification in thermal storages). Thus, such an approach contradicts the identified modeling requirements MR2 as well as MR3.

One possible solution, targeted in this paper, is to utilize detailed (co-)simulation setups for design optimization. Even though simulations are very well suited to characterize a given system design using a system performance measure, their application in the context of design optimization is more challenging. From a mathematical programming point-of-view, the simulation-based design leads to objective functions that must, in general, be considered non-linear, multi-modal and discontinuous [22]. To make matters more complicated, the evaluation of the objective function is computationally expensive (minutes to hours or even days per evaluation), depending on the complexity of the simulation model. With the number of possible design variables being high and the range of corresponding input parameter values being huge or even infinite, it becomes infeasible to perform this search by hand or by brute-force. Thus, this problem class requires to either reduce the allowed solution space and/or to use efficient black-box optimization algorithms, where finding a global optimum within finite time is not guaranteed [23].

Nevertheless, simulation-based design optimization is a frequently used technique as it allows the use of detailed models. The possibility to use high-fidelity models for optimization is especially relevant when targeting systems-of-systems, such as hybrid thermal-electric networks in this work, and, thus, satisfies KP2. In energy-related research, applications to building design are most frequent. The dynamic simulation tool IDA ICE together with a multi-objective optimization was used to find cost-optimal energy performance renovation measures for educational buildings by Niemela et al. [24]. Delgarm et al. [25] used the building energy simulation program EnergyPlus and a multi-objective particle swarm optimization to find the building specifications that minimize annual energy consumption. Many more similar studies exist and Nguyen et al. [26] provided an extensive summary of building related simulation-based optimization methods. Recent work uses simulation-based optimization in the context of district heating system design. Wang et al. [27] optimized the hydraulic design of variable-speed pumps in multi-source district heating networks using static hydraulic simulation and a genetic optimizer. Van der Heijde et al. [28] tried to find the cost optimal location and size of thermal storage tanks in district heating networks using

dynamic thermal-hydraulic simulation combined with optimal control and a genetic multi-objective optimization algorithm.

#### *3.3. Summary and Conclusion Regarding the Relevant State of the Art*

Based on the review of the relevant state of the art above, it follows that established optimization models do not provide the level of detail required for assessing hybrid thermal-electrical distribution grids, especially in view of closed-loop control and operation. In practice, currently only the co-simulation of domain-specific technical models can provide the desired degree of accuracy. However, even though these technical models are commonly used by domain experts for detailed assessments, they are typically not designed nor intended to be used for optimization. The literature gives examples of how co-simulation and this class of technical models can be used for optimization, but there exists—to best knowledge of the authors—no general methodology to support this process.

In this context, this work is the first to introduce simulation-based optimization for coupled district heating and electric network simulation combined with closed-loop control. A methodology to assist the co-simulation-based design of coupled heat and power networks is presented by providing conceptual guidance and a proof-of-concept implementation to the above mentioned problem setting. The focus lies on the utilized technical models (see Section 4) and their integration into suitable optimization approaches (see Section 5), not on the formulation of specific optimization algorithms.

Table 2 summarizes this situation in terms of the KPs and MRs identified above (+, full compliance; ◦, compliance with effort; −, insufficient compliance).



#### **4. Simulation Approach for Hybrid Networks**

This section presents the approach used in this work for the detailed simulation of coupled heat and power networks including control. The overall model is highly complex, exhibits non-linear behavior and has no closed algebraic formulation. On the one hand, this rules out the utilization of the most commonly used optimization approaches (e.g., LP or MILP). On the other hand, this additional complexity is unavoidable for analyzing the technical details of control and operation in such networks. In view of the general optimization methodology presented in Section 5, the presented approach can be regarded as a representative example of co-simulation approaches used for technical assessments. As such, it highlights the differences between the class of simulation models targeted by this work and the class of models typically used for optimization.

#### *4.1. Co-Simulation Environment*

A co-simulation approach enables the coupling of the different modeling paradigms, i.e., a transient thermal-hydraulic model, a quasi-static power flow model and time-discrete advanced control models. Thus, the influence of time-discrete *advanced control systems*, e.g., using rule-based control or model predictive control (MPC), on the dynamic *physical system*, i.e., the electric and the district heating (DH) network including supplies, coupling units and consumers, can be studied. This enables the assessment of hybrid thermal-electrical distribution grids with appropriate spatial and temporal resolution for relevant use cases [12].

The modeling activities and environments used in this work are split into the control model and the physical system model, with the latter only including low-level control, e.g., PID controllers. The assessment method is based on modeling tools according to the state of the art for each domain that are presented in more detail in the following sections.

Within the context of this work, the FUMOLA environment [29] has been used as co-simulation environment. FUMOLA is specifically designed to support the features offered by the Functional Mock-up Interface (FMI) specification [30], which defines a standardized application programming interface (API) and model description for both co-simulation and model exchange. FMI was selected as it is a mature, non-proprietary specification, developed by both academia and industry. The FMI standard enables the exchange and extension of tools and methods for the different domains and, thus, makes the approach highly versatile and extensible, especially in selecting the most appropriate method for advanced control system implementation.

#### *4.2. District Heating Network Model*

Thermal networks (including producers, network, thermal storages and consumers) are modeled in Modelica/Dymola [6,31] with the help of the DisHeatLib library [32] that is built upon the IBPSA library [33]. These open-source libraries include models for the most relevant components of district heating networks, considering bi-directional mass flows, heat transport delays, detailed substation and storage models and other thermo-hydraulic aspects that are highly relevant in heat networks. In addition, it provides models of local controllers and interfaces to electric networks. In summary, this modeling approach captures transient thermal and quasi-static hydraulic network behavior. The most relevant models are shortly presented in the following.

All DH pipes are modeled using a plug flow approach. The outlet temperature *T*out and, thus, the heat loss of a fluid parcel passing a pipe is described as:

$$T\_{\rm out} = T\_{\rm g} + \left(T\_{\rm in} - T\_{\rm g}\right)e^{-\frac{r}{R^{\rm C}}} \tag{1}$$

It depends only on the initial temperature *T*in, residence time *τ*, undisturbed ground temperature *T*g calculated using the Kusuda equation [34], thermal resistance of the pipe *R* and heat capacity of the water in the pipe *C*. The pressure drop mass flow correlation along the pipe is given by

$$
\dot{m} = \text{sgn}(\Delta p) k \sqrt{|\Delta p|} \tag{2}
$$

where *k* is the constant flow coefficient calculated for nominal conditions using the Colebrook equation for turbulent flow in rough pipes [35]. Details about implementation and experimental validation can be found in [36]. Heat exchangers in the DH substations are modeled with a variable effectiveness using a number of transfer units approach [37]. Valves are modeled using the above pressure drop and flow rate correlation with the flow coefficient *k* depending on the opening control signal.

Thermal energy storages are modeled using a vertically discretized multi-node approach to account for stratification and buoyancy [38]. Heat pumps are modeled using a Carnot-efficiency-based approach, electric heaters use a constant efficiency and gas and biomass boilers as well as combined heat and powers (CHPs) use heat generation dependent efficiencies. The main district heating supply unit is modeled as an ideal heat and differential pressure source and with no limits on maximum/minimum power or ramp rate. A fixed supply temperature is assumed for all generators.

A full list of model formulations and implementations can be found in the open-source Modelica libraries DisHeatLib and IBPSA.

#### *4.3. Electrical Distribution Grid Model*

Electrical distribution grids (including producers and consumers) are modeled with DIgSILENT PowerFactory [39], an engineering tool targeting primarily professional users. The quasi-static assessment of the electrical networks covers the most important features, such as voltage fluctuations and time-varying loads and generation. DIgSILENT PowerFactory's simulation interface has been extended to enable a series of consecutive power flow calculations in a co-simulation [40].

The power flow equations for node *a* in an *N*-node system are given in complex form:

$$\frac{S\_a^\*}{V\_d^\*} = \sum\_{b=1}^N Y\_{ab} V\_b \tag{3}$$

where *S*∗ *<sup>a</sup>* and *V*<sup>∗</sup> *<sup>a</sup>* denote the complex conjugate apparent power and voltage at node *a*, respectively, *Yab* denotes the bus admittance matrix and *Vb* denotes the complex voltage at node *b*. This results to 2*n* equations for the 4*n* unknowns, i.e., voltage, active/reactive power and voltage angle.

Coupling units are modeled as PQ buses, where active and reactive power is known, using the average power consumption/generation from the dynamic DH network model over the synchronization interval Δ*t*qs as input. Transformer units connect the low-voltage electrical networks to an external grid, i.e., modeled as a slack bus that determines the voltage and phase at the connection point.

#### *4.4. Operation and Control Models*

General-purpose tools such as MATLAB or Python can be integrated into the co-simulation with the help of FMI-compliant interfaces [41,42]. After each synchronization interval Δ*t*ctrl, they are called with the latest simulation outputs (measurement data), based on which they calculate and return new control setpoints that are then fed back to the physical models (feedback loop). This facilitates the implementation of a potentially large range of different types of algorithmic approaches for system-level controllers, including rule-based or model-predictive control schemes (see below). Local controllers on the process level are implemented within their respective subsystem models.

#### 4.4.1. Rule-Based Operational Strategies

Rule-based control generally relies on a list of *i* deterministic rules formulated as logical and/or algebraic expressions, e.g., a hysteresis controller that issues on/off signals. These rules constitute the knowledge base to determine a control action *ψ<sup>i</sup>* for a set of inputs (measurement data) corresponding to a certain system state *φ<sup>i</sup>* .

$$(\phi\_1 \to \psi\_1) \land \dots \land (\phi\_i \to \psi\_i) \tag{4}$$

Constructing this list of rules relies on expert knowledge. The reasoning depends on the specified operational goals and is often highly case-specific.

#### 4.4.2. Model-Based Operational Strategies

Model-based control algorithms, in comparison, rely on knowledge about the dynamic system behavior **x**˙(*t*) and the impact of control actions **u** ∈ **U** to govern the overall system along an optimal trajectory. The model can be used within an optimal control scheme to determine the control action that satisfies all system dynamics and yields an optimal performance metric *J*:

$$\begin{aligned} \min\_{\mathbf{u} \in \mathbf{U}} J &= \int\_{t\_s}^{t\_f} L(\mathbf{x}(t), \mathbf{u}(t)) dt \\ \text{s.t. } \dot{\mathbf{x}}(t) &= f(\mathbf{x}(t), \mathbf{u}(t)) \end{aligned} \tag{5}$$

where *L* is a cost function, *ts* is the start and *tf* is the final time of the control horizon. The continuous time optimal control problem is often transformed into a time discrete version where the dynamics are represented by a (linear) state-space system. At runtime, feedback from the system (measurement data) and predictions of disturbances can be used to provide safe operation and optimal performance with respect to the operational goals. The use of model-based optimal control schemes is often labor intense, especially if model identification is not automated, requires a suitable optimization algorithm and solving the mathematical programming problem might be time consuming.

#### **5. Design Optimization and Control of Hybrid Networks**

This section describes the proposed design optimization framework for hybrid thermal-electrical distribution networks, focusing on the *technical assessment* of such systems for the purpose of network planning and operation. It is based on a simulation-based optimization approach that utilizes the detailed coupled heat and power network simulation including closed-loop control, presented in the previous section. The general resulting design optimization problem is given by

$$\mathbf{z}^\* = \operatorname\*{argmin}\_{\mathbf{z} \in \Omega} c\left(\mathbf{g}(\mathbf{z})\right) \tag{6}$$

where the function *g*(**z**) involves one (co-)simulation run for a specific system configuration **z** and **Ω** describes the solution space, i.e., the set of all possible and allowed system configurations. The goal of the design optimization process is to find the system configuration **z**∗ that minimizes the objective function. Due to the high computational burden involved in executing one call of *g*(·), it is important to a priori reduce the number of possible system configurations. Hence, the proposed design approach assumes that most basic design decisions have already been reached using either expert know-how and/or established design optimization tools, e.g., the choice of conversion and storage technologies has been made employing mixed-integer linear programming techniques.

The *objective function c*(·) relates results from one simulation run *g*(**z**) for a specific system configuration **z** to the optimization criterion, such as costs or technical key-performance indicators. Thus, the objective function maps certain technical and/or economical aspects of the overall system to a numerical (scalar) value. In the case of multi-carrier energy systems, objective functions typically relate aspects of the overall system that are traditionally treated by separate engineering domains, e.g, total energy imports for both heat and power. Furthermore, objective functions may evaluate effects that result from dynamic interactions between the subsystems, especially synergies among production, consumption and storage and their impact on network operation.

#### *5.1. Influence of Operational Strategies on Optimal Design*

The objective function of a given system configuration is highly dependent on the performance of the respective operational strategy and implemented control. Hence, to yield a small value for the objective function, the *operational goals* (e.g., the use of local PV generation for heat pump operation) should be in-line with the design *optimization targets* (e.g., sizing of heat pump to increase PV self-consumption). To this end, the operational goals and design optimization targets need to be translated into an appropriate control implementation.

In this work, two categories of control schemes, i.e., rule-based and model-based, are considered in more detail in Section 4.4. From the point of view of design optimization, these two categories of operational strategies serve very different purposes:


This difference comes from the fact that model-based optimal operational strategies are—by design—able to guide the system evolution in accordance to the defined optimization targets. In contrast, a rule-based approach is always heuristic, such that there is a priori no such guarantee. In this case, both a different set of rules or a different system configuration could yield a performance improvement with respect to the objective function.

#### *5.2. Design Optimization Approaches*

The design process needs to adapt to the specific circumstances of any given project, especially constraints regarding available design options. For instance, there may be economical or legal restrictions or technical constraints (especially due to already existing infrastructure). In practice, this has a strong influence on the number of possible system configurations (see, for instance, the applications in Section 6). Hence, the following two complementary approaches are introduced, taking into account the size of the solution space and the implemented control scheme:


In considerations of the above, the choice between OCS and HPS represents a trade-off among implementation effort, computational complexity and usability of the operational strategies for design optimization. Figure 3 visualizes this trade-off in terms of the size of the solution space, the computational complexity of the utilized operational strategy and the implementation effort for the associated sub-tasks (i.e., controller development and candidate selection). A combination of OCS and HPS, i.e., metaheuristic design optimization on top of an optimal control scheme, for a large number of possible system configurations, although preferable and theoretically possible, is not considered due to the high computational effort involved.

Figure 4 visualizes the different optimization approaches of OCS and HPS. In OCS, the optimizer is an integral part of the operational strategy implementation, guiding the system evolution towards optimal operational performance. In HPS, the optimizer is separated from the co-simulation and design optimization is achieved by repeatedly executing the co-simulation with different parameters. In this context, another benefit of using co-simulation becomes apparent, because the choice between both approaches has virtually no impact on the modeling of the thermal and electrical domain, as long as the operational strategy is implemented as an individual and exchangeable component in the co-simulation.

computaonal complexity of operaonal strategy

**Figure 3.** Comparison of the Optimal Control Scan (OCS) and the Heuristic Parameter Scan (HPS), indicating higher (↑) and lower (↓) implementation effort for specific sub-tasks.

**Figure 4.** Schematic view of the relation between co-simulation and optimizer for: (**a**) the Optimal Control Scan (OCS); and (**b**) the Heuristic Parameter Scan (HPS).

#### *5.3. Optimization Algorithms, Decision Variables and Objective Functions*

The choice between OCS and HPS is first and foremost conceptual, providing a guideline for implementing the actual optimization. For both cases, a large variety of optimization algorithms exists that can be applied, thanks to the flexibility of co-simulation approaches. In the case of HPS, any metaheuristic population-based optimization algorithm can be applied that can take the results from individual co-simulation runs as black-box input, such as PSO [43], Differential Evolution [44] or PSwarm [45]. In the case of OCS, the co-simulation approach enables for instance the integration of existing toolboxes for model-predictive control, LP or MILP for MATLAB or Python. For example implementations—without loss of generality—refer to Section 6.

Using co-simulation of domain-specific models as basis for system assessment allows including more detailed technical information for decision variables compared to traditional optimization approaches. However, the choice between OCS and HPS does have practical implications for the selection of decision variables and objective functions. In the case of HPS, the decision variables used by the optimal control instance have to be based on the measurement data of the current and previous simulation time steps. For HPS, in contrast, the optimizer has not only access to all measurement data but also the full set of results of each completed co-simulation run. This means that overall key performance indicators (e.g., total energy saving or yearly local self-consumption) can be included in the objective function. This distinction is key for understanding the conceptual and qualitative difference regarding the optimality of results in OCS and HPS. Nevertheless, from a quantitative and practical point of view, both approaches yield (near) optimal results.

#### **6. Proof-of-Concept Applications**

The applicability and usefulness of the proposed design approach is demonstrated with the help of two real-world applications. Both applications focus on the design optimization of hybrid thermal-electrical grids from a technical perspective, aiming at exploiting synergies between the networks by mutual support during operation. Hence, no monetary optimizations are applied, i.e., neither investment nor operational costs are considered.

Both application examples aim at local consumption of excess PV generation using heat pumps or electric boilers. However, the proposed design approach is not limited to this kind of applications and could also be applied to applications with a focus on short-term and/or long-term storage, peak shaving, or others.

#### *6.1. OCS Example Application: Suburban Industrial Area*

#### 6.1.1. Technical System Layout

The system is located in a suburban industrial park area, comprising multiple office buildings and industrial facilities at four adjacent sites, with a total yearly electrical demand of 1 GWhel and a total yearly thermal demand of 2.3 GWhth. Figure 5 presents a schematic of the technical system layout.

**Figure 5.** Schematic overview of the hybrid system layout for the suburban industrial area (GB, gas boiler; BB, biomass boiler; GWHP, ground water heat pump; WWHP, waste water heat pump; PV, photovoltaic module; CHP, combined heat and power plant). The components targeted by the design optimization process are highlighted in orange.

The on-site low voltage network consists of 0.6 km of cables and 0.6 km of overhead lines, and has five medium-size PV systems with a total installed capacity of 272 kWpel connected to it. The heat network connects the thermal generators and demand sites, which makes it possible to share heat between them. Three of the sites use previously installed gas boilers (GBs) with a total nominal capacity of 2.44 MWth, whereas the fourth site uses a CHP plant with a nominal capacity of 950 kWth, all four connected to thermal buffer tanks. Moreover, a biomass boiler (BB) with 950 kWth nominal capacity is feeding the heat network.

The main task in this application was to add and size a ground water heat pump (GWHP) and a waste water heat pump (WWHP), representing the *degrees of freedom* in the design process. These heat pumps and the CHP are the system's *coupling points* between the networks.

#### 6.1.2. Operational Strategy

The foremost goal of the operational strategy is to maintain an operational temperature between 80 and 95 ◦C in all thermal buffers, in order to guarantee that the thermal demand can be fulfilled at all times at an admissible temperature level. At the same time, the following *optimization targets* should be considered:


These operational targets can be translated into an operational strategy, which was implemented as model-based optimal control (compare with Section 4.4.2), based on a linear optimization problem formulation. It prioritizes the heat sources according to the following scheme:


At the same time, the controller keeps track of all operational constraints (matching of generation and demand, production thresholds, network constraints, etc.). This guarantees the optimal operation for a given system configuration.

#### 6.1.3. Design Optimization Strategy

For optimizing the system design, the heating capacities of the two heat pumps are the *degrees of freedom*. However, only a limited number of realistic sizing options has been identified beforehand by the owner, based on experience and expertise from domain experts as well as spatial and budget constraints. For instance, since WWHPs are more efficient than GWHPs, operation of the WWHP is prioritized, whereas the GWHP is only turned on if even more electricity from PV is available and heat is needed. Therefore, all selected configurations foresee a bigger WWHP in combination with a smaller GWHP. Furthermore, given the operational strategy explained above, which also aims at a minimization of the electricity consumption from the external grid, the maximal practical size of the heat pumps is limited by the maximal PV production. These considerations led to three potential system configurations (see Table 3), which differ in the sizing of the heat pumps.


With the limited number of system configurations and the possibility to translate the operational strategy into an optimal controller scheme, the OCS is the natural choice of *optimization approach* for this application (see Figure 6).

**Figure 6.** Schematic overview of the OCS workflow applied to the system design of the suburban industrial area.

#### 6.1.4. Results

As expected, the operational strategy is affected by the seasonal variations (temperature, solar irradiance), effectively resulting in seasonal operational modes that prefer different sources at different times of the year. These significant differences necessitate the evaluation of the system performance on a yearly basis, in order to provide a good basis for the choice of the optimal design.

Figure 7 summarizes the performance of the three considered system configurations in terms of yearly total energy production for the heat pumps (*E*HP tot ), the biomass boiler (*E*bio tot ) and the CHP (*E*CHP tot ). In all three cases, the thermal energy production of the gas boilers is around 230 MWhth/a. The figure shows that Configuration C—which has the largest WWHP—fails to exploit the full potential of the PV generation, resulting in the smallest total energy production of all configurations. This is due to the fact that the heat pump must not be operated below 80% of its maximum capacity, which in turn leads to a high threshold for turning it on in Configuration C. Configuration A—with the smallest WWHP—slightly falls behind Configuration B, which provides in this regard the best compromise as it exhibits the largest total generation from the heat pumps. Furthermore, in Configuration C, the WWHP's high operational threshold causes not only an increase of biomass-based generation, but actually leads to an overcompensation at the cost of the CHP-based thermal generation compared to the other configurations.

Analysis of the detailed results from the co-simulation shows that Configuration B has the most favorable impact on the electrical system in view of integrating the on-site PV production. For instance, it shows the most significant improvement of the voltage band usage compared to non-hybrid system layouts by reducing the time and amount the maximum voltage band exceeds the 10% threshold.

In conclusion, even though the improvements from the system level point of view are limited due to the actual available PV generation and other practical constraints (e.g., economical aspects of increasing the CHP capacity), Configuration B shows overall the best thermal and electrical system performance. It optimally exploits the on-site PV generation via the heat pumps (*E*HP tot = 225 MWhth/a) while at the same time providing a modest reduction of fuel-based generation from the CHP and the biomass boiler (Δ*E*fuel tot = −211 MWhth/a). Furthermore, the design goal of minimizing on-site CO2 is successfully met through the preference for the heat pumps, the biomass boiler and the CHP over the gas boiler (see Figure 8).

**Figure 7.** Comparison of the yearly thermal energy production of heat pumps, biomass boiler and CHP for the three considered system configurations (Configuration A, Configuration B and Configuration C).

**Figure 8.** Yearly load duration curve.

#### *6.2. HPS Example Application: Rural Residential Area*

#### 6.2.1. Technical System Layout

The system comprises a rural low-voltage electric grid with a total cable length of about 6.5 km. This network connects residential, commercial and agricultural customers, summing up to a total number of around 110 customers with a total annual power demand of around 775 MWhel and PV systems that are feeding a total of 384 MWhel per year.

The district heating network structure is typical for rural areas with low heat demand density (linear density of about 520 kWhth/a/m and peak load demand of about 2.0 MWth). An outdoor-temperature dependent heating curve is used to set the supply temperature between 90 (winter) and 70 ◦C (summer). Return and supply pipes connect around 60% of the buildings in the area with a base heat generation plant responsible to keep the differential pressures at the consumer substations above a minimum.

Additionally, three electric boilers, i.e., electric heaters combined with thermal storage tanks, are installed as *coupling points* between the networks. Adequate locations for the boilers were identified using a simulation-based sensitivity analysis in an early planning phase. These locations are especially prone to over voltage problems resulting from PV generation, an issue that might be eased by the active conversion of excess power generation. The volume of the thermal storage tanks and the capacities of the electric heaters are chosen as *degrees of freedom* of this system. The system layout is illustrated in Figure 9.

**Figure 9.** Schematic overview of the hybrid system layout for the rural residential area, depicting the electrical distribution network and the district heating network (EH, electric heater; TES, thermal energy storage). The components targeted by the design optimization process are highlighted in orange.

#### 6.2.2. Operational Strategy

Decentralized electricity generation from PV systems of local prosumers pose challenges to the electrical distribution grid. Especially in times when PV production is high and electrical consumption is low, the upper voltage limit in the network can be exceeded and lines or transformers can be overloaded. To mitigate some of these problems, the operational strategy foresees to utilize excess power from PV overproduction. The local power grid is regarded as virtual power plant (VPP) and a supervisory controller tries to use as much of the excess power, i.e., negative residual load, locally via the installed electric heaters. Based on these ideas, a suitable operational strategy was implemented in Python using the following rule-based scheme (compare with Section 4.4.1):


#### 6.2.3. Design Optimization Strategy

The high number of possible system configurations and the use of a rule-based operational strategy makes HPS the most appropriate *optimization approach* for this application. Figure 10 illustrates the overall simulation-based optimization procedure.

The design optimization of the coupling units, i.e., finding optimal sizes for the defined *degrees of freedom*, rests upon the following *optimization targets*:


These targets are combined into a single objective function using appropriate weighting factors. The six-dimensional solution space is reduced by introducing bound constraints for the degrees of freedom, avoiding extremely high and low tank volumes and electric heater capacities.

The implementation relies on the use of a dedicated open-source tool for simulation-based optimization [46], which allows parallelization on multiple machines and the use of optimization algorithms specialized for efficient black-box optimization. In this case, the freely available PSwarm solver [45] is used that combines particle swarm optimization and pattern search for efficient global optimization. The pattern search relies on a coordinate search method that is responsible for local convergence, whereas the population-based particle swarm algorithm performs a global search enabling the exploration of the whole design space. The stopping criterion is based on a maximum number of co-simulation runs, i.e., objective function calls.

**Figure 10.** Schematic overview of the HPS workflow applied to the system design of the rural residential area.

#### 6.2.4. Results

The convergence of the objective function value over the number of co-simulation runs, i.e., assessed system configurations, is shown in Figure 11. The objective function values are shown relative to a reference solution reflecting the current status quo, i.e., without any active coupling between the networks using the electric heaters. The optimization algorithm is able to outperform the reference scenario within only a few simulation runs. After around 450 runs, the algorithm converges to a minimum and stops after it exceeds the maximum number of objective function calls. Although this might only be a local minimum, the found minimum objective function is around 5.7% lower compared to the uncoupled case. The corresponding optimal system configuration exhibits electric heaters with a total capacity of 260 kWel and thermal storage tanks with a total volume of 20 m2.

The impact on the two networks for the uncoupled case and the optimal design case in terms of load duration curves is shown in Figure 12. A total of around 220 MWh electricity is converted in the optimal design case using electric boilers, illustrated by the colored areas. It can be seen that electric heaters in combination with thermal storage tanks are able to use a significant share of the excess power generation from PV by converting it into heat that is fed to the district heating network. Due to the seasonality of PV generation, district heating is mainly affected in low heat demand times, i.e., summer.

**Figure 11.** Convergence of objective function shown relative to the reference scenario without coupling between the networks, i.e., no electric boilers.

**Figure 12.** Load duration curves for the electric network (**a**) and the district heating network (**b**) for the uncoupled and the optimal design system configurations.

#### **7. Conclusions**

This work presents a simulation-based design approach for hybrid thermal-electrical distribution grids. The approach addresses the technical challenges of both domains while at the same time emphasizing their mutual control and operation, in order to exploit hitherto unused synergies in production, storage and consumption. As such, this approach is a significant extension of established tools for the design optimization of multi-energy systems.

Co-simulation is recommended for the technical assessment as a means to bridge the gap between single-domain simulation tools and the multi-domain target of investigation, providing a viable and practical approach to involve experts from different domains. A novel methodology is presented that links complementary design approaches to suitable operational strategies and optimization methods according to the state of the art. This enables the exploitation of synergies in the control and operation of hybrid thermal-electrical distribution systems in an optimal way. Furthermore, an implementation of the proposed approach based on state-of-the-art tools is presented. Finally, the applicability of the presented approach is demonstrated in two real-world applications.

Future work could apply this approach to other network-related multi-energy applications, such as design of short-term and long-term storage, peak shaving and others. A method for representative days selection for coupled heat and power networks could avoid the need for full-year simulations and, thus, reduce computational time. The use of computationally efficient model-based

optimal control models or faster models for the physical system could enable the combination with meta-heuristic design optimization in a tolerable run time.

**Author Contributions:** conceptualization, E.W. and B.L.; methodology, E.W. and B.L.; software, E.W., B.L., D.B., S.H. and T.F.; investigation, E.W., B.L., D.B., S.H. and T.F.; writing—original draft preparation, E.W. and B.L.; writing—review and editing, D.B., S.H. and R.H.; visualization, E.W. and B.L.; supervision, R.H. and E.W.; project administration, E.W. and B.L.; and funding acquisition, E.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** Parts of this work were funded by the Austrian research funding association (FFG) within the scope of the research project *OptHySys—Optimierung Hybrider Energienetze und -Systeme* (project #848778). Parts of this work were carried out within the research project *SmILES—Smart Integration of Energy Storages in Local Multi Energy Systems for maximising the Share of Renewables in Europe's Energy Mix* and received funding from the European Commission's Horizon 2020 Research and Innovation Programme under the Grant Agreement No. 730936.

**Acknowledgments:** The authors want to acknowledge support provided by the cooperation college *Smart Industrial Concept* (SIC!).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Vehicle-To-Grid for Peak Shaving to Unlock the Integration of Distributed Heat Pumps in a Swedish Neighborhood**

#### **Monica Arnaudo \* , Monika Topel and Björn Laumert**

Royal Institute of Technology, 11428 Stockholm, Sweden; monika.topel@energy.kth.se (M.T.); bjorn.laumert@energy.kth.se (B.L.)

**\*** Correspondence: monica.arnaudo@energy.kth.se

Received: 17 February 2020; Accepted: 2 April 2020; Published: 3 April 2020

**Abstract:** The city of Stockholm is close to hitting the capacity limits of its power grid. As an additional challenge, electricity has been identified as a key resource to help the city to meet its environmental targets. This has pushed citizens to prefer power-based technologies, like heat pumps and electric vehicles, thus endangering the stability of the grid. The focus of this paper is on the district of Hammarby Sjöstad. Here, plans are set to switch from district heating to heat pumps. A previous study verified that this choice will cause overloadings on the electricity distribution grid. The present paper tackles this problem by proposing a new energy storage option. By considering the increasing share of electric vehicles, the potential of using the electricity stored in their batteries to support the grid is explored through technical performance simulations. The objective was to enable a bi-directional flow and use the electric vehicles' (EVs)' discharging to shave the peak demand caused by the heat pumps. It was found that this solution can eliminate overloadings up to 50%, with a 100% EV penetration. To overcome the mismatch between the availability of EVs and the overloadings' occurrence, the minimum state of charge for discharging should be lower than 70%.

**Keywords:** vehicle-to-grid; heat pumps; integrated energy systems

#### **1. Introduction**

Sweden has adopted the ambitious goal of becoming a zero carbon emission society by 2045 [1]. Within this context, the energy sector has attracted a special attention because it accounts for more than 70% (including transport) of total greenhouse gas emissions [2]. When looking at the internal electricity mix of this country [3,4], electricity is promoted as a clean resource and thus a promising option to achieve the 2045 target. As a consequence, in Stockholm, citizen-driven initiatives work to promote the installation of distributed residential heat pumps (HPs) and the adoption of electric vehicles (EVs) [5]. A real example is represented by Hammarby Sjöstad, which is a residential neighborhood constituted by multi-apartment buildings. This district is currently connected to the city's district heating (DH) network, but projects are already approved for some customers to switch to their own domestic HP.

On an opposite side, this electrification trend is challenged by a structural problem that impacts the power grid system in Stockholm. In fact, the electricity transmission to the city is close to its maximum power capacity. This means that the distribution system operators (DSOs) are being forced to impose power limits on the distribution grid [6]. Thus, the connection of new loads, like HPs and EVs, could soon be impossible. If capacity expansion investments are to be avoided or delayed, alternative solutions should be quickly found.

A promising option is offered by energy storage [7]. This technology can shift loads in time by charging when surplus energy supply is available and discharging during periods of peak demand. This can be done on both the heat and on the electricity sides, depending on the type of energy storage technology used.

For this paper, the focus was set on the electric power side by considering the electricity storage potential in the batteries of EVs. In this context, an EV is regarded not only as a load to the power grid but also as an extended capacity for supply. This is based on the vehicle-to-grid (V2G) concept. V2G technology enables a bi-directional flow between an EV's battery and the grid [8–11]. Thus, both the charging and discharging of electrical energy are allowed.

The background of this study was linked to a real initiative called "Charge at Home" ("Ladda Hemma" in the original Swedish). Within this context, several housing associations in Hammarby Sjöstad were encouraged and assisted in the process of installing EV charging infrastructures in their buildings. The study presented in this paper explores the potential of enabling a V2G bi-directional flow at these stations. By assuming that distributed domestic HPs are installed in Hammarby Sjöstad, the objective was to use the electricity stored in the EVs' batteries to cover the peak power demand generated by these HPs. The main challenges were the availability of EVs parked at the stations and the level of charge of their batteries.

The importance of this study was stressed by a previous assessment that the installation of distributed HPs in Hammarby Sjöstad will cause overloadings on the local electricity distribution grid [12]. This was found by assuming that DSOs would impose a maximum loading limit equal to 100% for each grid's cable. The possibility of using the thermal mass of the buildings as thermal energy storage (TES) by controlling the thermostats of the apartments was also studied. Despite the implementation of this solution, grid overloadings were still detected. The present work builds on this last outcome by verifying how much V2G technology can further contribute to balance grid overloadings.

Several papers have analyzed the possibility of using V2G for peak shaving and avoiding grid upgrades. Most of them [13–18] proposed optimal control algorithms by looking at costs and emissions. Some of them also took battery degradation [19,20], market composition [21], and incentives [22] into account. The studies of [23,24] combined the capacity provision market with the spinning reserve. The authors of [25] showed the advantage of combining a V2G peak shaving strategy with a domestic battery and a photovoltaic (PV) system.

Despite the broad perspective, none of these papers considered the potential synergies and impacts between the power, transport, and heating sectors.

The authors of [26] implemented V2G with a time-of-use tariff. This study also showed the impact of considering the indoor temperature of buildings as a constraint on the supply dispatch of a micro-grid. However, their perspective was limited to one building.

The energy hub presented in [27] included the heating and cooling demands of a residential community. How the implementation of V2G can impact electricity and cooling prices was shown. However, impact and synergies among these energy carries were not considered from a technical standpoint.

Given these literature gaps, one objective of this paper was to propose a district-level perspective on potential synergies among the heating, electricity, and transport sectors. To approach this sector-coupling problem, the co-simulations of dedicated models [12,28] and stochastic profiles were combined. Thus, the technical performances of technologies belonging to these different sectors were linked.

Here, the focus was on the case of a multi-apartments neighborhood located in Stockholm, Hammarby Sjöstad. By taking into account that distributed HPs can overload the local grid (reference scenario from [12] as a main objective, this paper aimed at estimating the potential of V2G to alleviate this problem (V2G integration scenario). This was done by using the energy stored in the EV's batteries to shave the HPs' peak power demand.

The remainder of the paper is as follows. Firstly, the case study is presented in Section 2. Secondly, the two scenarios and the corresponding simulation models are described in Sections 3 and 4, respectively. Finally, the results are discussed in Section 5.

#### **2. Case Study**

The case study presented in this paper corresponded to a specific area in Hammarby Sjöstad, as shown in Figure 1. All the buildings in this area are connected to a single medium-to-low voltage transformer substation. The conceptual illustration of the case study area in Figure 1 sketches the current type of energy infrastructure in the neighborhood. The multi-apartment buildings are connected to a DH network for space heating and domestic hot water purposes. The installed capacity is about 1.2 MW for a total of 94,555 m2 of heated area. Electricity is provided by a low voltage distribution grid (400 V) that covers a capacity of about 1.6 MW.

**Figure 1.** The case study area: a conceptual illustration and a satellite view BB: building block; DH: district heating.

Within the scope of the scenarios later described in Section 3, the buildings are grouped in 22 building blocks (BB) according to the energy declaration documents of the corresponding housing associations [29]. As reported in Table 1, the BBs are identified through their heated area.


**Table 1.** Square meters of heated area for each BB.

The choice of a limited area was motivated by the focus of the study, which was to show the potential of peak power demand management by using V2G technology. This could be easily replicated to show the impact on the whole neighborhood (or beyond).

#### **3. Scenarios**

#### *3.1. Reference Scenario*

The reference scenario for the present work was taken from a previous study presented in [12]. Given the interest of a few housing associations (Section 1), the focus was to assess the technical

feasibility of replacing a DH network with domestic-distributed HPs. Since the sole installation of the distributed HPs was shown to cause overloadings on the electricity distribution grid, a demand response solution was further assessed. The thermal mass of the buildings was used as TES to alleviate the grid from the overloadings caused by the HPs. A positive contribution was shown with a reduction of the heat demand dependence from a DH of 6%, compared to the case without the thermal mass control. Furthermore, up to 50% fewer overloadings were detected at the HP level. However, it was also highlighted that a full disconnection from the DH network was not possible given the current infrastructure capacity. These outcomes [12] represented the reference scenario for the present study.

#### *3.2. V2G Integration Scenario*

In the present paper, the reference scenario was extended by assessing the potential of V2G technology as a way to further support the electricity distribution grid. A simplified illustration of this concept is shown in Figure 2. For the sake of simplicity, a BB is represented as one building, though it can also correspond to a group of buildings. The V2G charging/discharging station is located upstream of each HP installation. The EVs' charging/discharging patterns were based on stochastic profiles. These were generated according to the assumptions described in Section 4.

**Figure 2.** Conceptual illustration of the case study within the vehicle-to-grid (V2G) scenario. HP: heat pump; EV: electric vehicle.

The EVs' parking stations are located in the BBs, in line with the "Charge at Home" initiative (Section 1). In this scenario, as illustrated in Figure 2, both charging and discharging processes—and thus V2G—were enabled in each BB. The aim was to cover the grid overloadings by performing peak shaving with the electricity stored in the EVs' batteries. This study considered a 100% penetration of the EVs. Furthermore, the EVs' charging was assumed to happen overnight with no overloadings. More specific assumptions are described in Section 4.

#### **4. Models**

Different combinations of models were used to represent the two scenarios described in Section 2. These models are discussed in detail in the following sub-sections.

#### *4.1. Models for the Reference Scenario*

This sub-section briefly highlights that, from a modeling and simulation perspective, a co-simulation method [28,30] was used to assess the interaction between the loading of the electricity distribution grid, the heat supply from HPs and a DH network, and the indoor temperature of residential buildings. Figure 3 sums up the modeling tools that were used for each system involved in the scenarios. The scheme also shows the main heat (red lines), power flows (black lines), and control signals (blue dashed lines).

**Figure 3.** Co-simulation scheme for the reference scenario (EB: electrical backup; TH: thermostat).

The co-simulation environment allowed us to capture two feedback loop signals. On one hand, the loading status of the grid was checked to decide whether the HPs could be operated. An overloading signal was generated when 100% of the capacity was reached. On the other hand, the indoor temperature of the BBs was monitored to decide if the thermostat set point could be changed to load/unload the thermal mass. In this case, the indoor temperature was to be maintained within the range of 20 ± 0.5 ◦C.

The following models were implemented for this reference scenario:


These models were used within a one year simulation with an hourly time step. A detailed explanation, parameters, variables, and assumptions for all the models are provided in [12].

#### *4.2. Models for the V2G Integration Scenario*

Concerning the V2G integration scenarios, two models were combined:


Regarding the first model, a critical day was selected (10th January) from the reference scenario, together with the identified critical cables. With reference to this day, the overloading power for each critical cable was obtained by subtracting the power corresponding to a 100% loading of each critical cable from the actual loading. For the hours that did not present overloadings, the overloading power was set to zero.

In this study, the analysis was extended to all cables, including the parallel ones. Figure 4 shows the one-line diagram of the electricity distribution grid of the case study. The four red-highlighted cables are the ones that were identified as critical in the reference scenario. Two of these cables corresponded to bundled parallel cables. The corresponding labels are listed as the following:


**Figure 4.** One-line diagram of the case study 400 V grid.

The stochastic driving patterns of the EVs were generated by means of an in-house, Python-based model. Table 2 summarizes the main inputs to the model, with the related references for the assumptions. These values were considered valid for the selected location. The number of cars per person referred to statistical data of a typical Stockholm neighborhood. To determine the density of inhabitants, it was reasonably assumed that two people could live in 48 m2. The level of penetration of EVs was related to a future case, in line with the current trend of an increasing EV market share [34]. These three inputs were used by the model to determine the number of EVs per each BB.


**Table 2.** Inputs to the in-house EVs' driving profile model. Soc: state of charge.

From further statistical data, it was possible to assume average driving kilometers and speeds. The state of charge (SoC) of an EV, which refers to the charge level of its battery, was constrained by charging and discharging limits. An EV had to be charged if its SoC was below 50%, while charging was prevented if its SoC was over 90%. This is mainly due to technical performance reasons of the battery. Discharging to the grid was enabled when an EV's SoC was over 80%. Since this represented a key parameter for the present study, a sensitivity analysis was performed for values between 70% and 90% of the minimum SoC for discharging. The number of charging stations for each BB and the number of trips per day referred to the context of a residential neighborhood. This was also in line with the initiative "Charge at Home," presented in Section 1. Thus, the two trips represented a way to work and a way back to home. The potential charging/discharging of an EV could happen after these two trips were completed each day.

Assumed initial hourly profiles for the share of cars starting a trip at a certain hour were taken as a further input to the model. These profiles were randomized using the parameters' ranges presented in Table 2.

Given all these inputs (Table 2), the in-house Python model performed three main steps:


The main outputs of the model were:


Finally, the EVs' available power profiles were combined with the overloading power profiles taken from the reference scenario. The objective was to estimate the potential of V2G to alleviate the grid overloadings caused by the installation of the distributed HPs. Thus, a peak power shaving strategy was implemented by using the electricity stored in the EV's batteries. This was done by overlapping the two hourly profiles along the selected critical day. In particular, a V2G discharging service was activated when the following combination of events occurred:


The validity of the results was bound to the assumption that the EV's driving patterns did not present large differences from one day to the other. This was considered reasonable for a residential district, like Hammarby Sjöstad.

Since it was expected that the assumption for the minimum SoC for discharging played a relevant role, a sensitivity analysis was conducted for minimum SoC values equal to 70%, 75%, 80%, 85%, and 90%.

#### **5. Results**

Within the scope of the present paper, the results are presented and discussed in Figure 4 and in the remainder of this section.

Figure 4 highlights the cables that were found to be still critical within the reference scenario. This means that, despite the utilization of the thermal mass of the BBs as a TES device, not all the grid overloadings, caused by the HPs, could be compensated for [12]. This was mainly due to the constraints in terms of thermal mass capacity of each BB and indoor temperature comfort. The latter was assumed to 20 ± 0.5 ◦C. By introducing V2G technology, according to the models discussed in Section 4.2, the objective was to explore the technical potential of EVs' batteries to provide further capacity to the grid. In this way, it could be shown how V2G can help covering these remaining overloadings generated by the HPs' demand peak.

Concerning the V2G integration scenario, Table 3 shows the number of EVs that were assigned to each BB. This outcome was based on the parameters presented in Table 2, regarding the number of cars per person, the number of people per BB's heated square meters, and the share of EVs in the studied area.


**Table 3.** Number of EVs per each BB.

In terms of electricity discharging to the grid, the number of EVs reported in Table 3 shows the total availability of EVs when no constraints about trips and SoC were considered. However, within the V2G integration scenario, the EVs were firstly required to be parked at their BBs' stations after completing two daily trips (to and from work). Secondly, their SoC had to be above 80%. Based on these constraints, Figure 5 illustrates the EVs' availability when back from work, i.e., after the second daily trip. When looking at the number of EVs parked (top graph in Figure 5), it could be verified that the cars were driven back from work along the day, and, around 6 p.m., most of them were parked at home. However, since an SoC higher than 80% was also required, the actual, much lower, EVs' availability is shown in the bottom graph of Figure 5. For example, at the charging station of BB16, at 6 p.m., only 25% of the parked EVs could discharge electricity to the grid.

**Figure 5.** EVs' availability at each BB station during the selected critical day.

Figure 6 illustrates the impact of the EVs' availability for V2G on the grid's overloadings. In this figure, the reference scenario (without V2G integration) and the V2G integration scenario are compared over the selected critical day. The comparison is done by showing the percentage overloading performance for each critical cable. The loading performance below 100% is not shown in order to set the focus on the overloadings only.

**Figure 6.** *Cont.*

**Figure 6.** Cables overloading performance day before and after the integration of V2G.

As it can be noticed in Figure 6, the overloadings were concentrated during the late afternoon and during the evening. This is reasonable since the studied area is mainly a residential one, so people are at work during the day hours.

The integration of V2G helped to alleviate the overloadings by shaving the peak demand generated by the HPs. All the cables, except Llugnw6 and Llugnw7, fully benefited from this new system with no overloadings left. The other cables remained critical, especially before 6 p.m., when most of the cars were still away from the parking stations. This means that, with a 100% level of EV penetration in the neighborhood, the discharge of electricity from EVs' batteries could cover the amplitudes of the critical power peaks. The time matching between overloadings and the EVs' availability remained a challenge.

From the perspective of the main objective of the case study, Figure 7 shows the relation between the cables and each BB and, thus, each distributed HP unit. Considering that this analysis was related to a critical day (representative according to the assumption in Section 4.2), it can be concluded that the integration of V2G has the potential to fully enable the installation of HPs in BB6, BB19, BB20, and BB21. Concerning the other BBs, further measures should be implemented. For example, if a grid upgrade is to be either avoided or delayed, a different V2G operation logic should be tested. The implementation of smart control devices based on forecasting methods is a suggestion for future work.

**Figure 7.** Cables overloadings over the selected critical day.

Finally, in order to assess the relevance of the assumption for the minimum SoC for discharging, a sensitivity analysis is presented for values equal to 70%, 75%, 80%, 85%, and 90%. It was found that the availability of cars with a minimum SoC of 70% allowed for the covering of all the overloadings during the selected critical day. However, this is expected to increase the following charging requirements, which should be further tested. On the opposite side, no cars were available when the minimum SoC was set to 90%, so none of the overloadings could be balanced. As an example, Figure 8 shows the results of this sensitivity analysis at 5 p.m. during the selected day for all the critical cables. At this specific hour, a minimum SoC of 75% already solved the overall criticality. However, the cable Llgnw7 was still overloaded by about 8% at 1 p.m.

**Figure 8.** Cables' overloadings at 5 p.m. for different minimum SoC values for discharging.

As a more general conclusion, it can be stated that the approach presented in this paper can play a relevant role as a decision support tool for city planners, energy utilities, and engaged citizens. A district-level perspective linked to a sector-coupling approach (heating, electricity, and transport) can unlock the implementation of innovative technologies like active thermal mass and V2G.

#### **6. Conclusions**

In this paper, a district level perspective was applied on an integrated energy infrastructure problem where the electricity, heating, and transport sectors were interconnected.

Hammarby Sjöstad, a residential neighborhood in Stockholm, was selected as a relevant case. A previous study assessed that the plan of installing distributed domestic HPs will overload the local electricity distribution grid. The criticality of this situation can be improved, but not be solved, by using the thermal mass of the buildings as TES.

In the present work, V2G was presented as a new solution to further support the grid. This technology enables a bi-directional flow of electricity between the EVs' batteries and the grid. The aim was to explore the potential of V2G to perform peak power demand shaving by discharging electricity to the grid. The objective was to compensate for the overloadings caused by the HPs.

The technical performance simulation over a representative critical day demonstrated that, with a 100% penetration of EVs, the available power for discharging could cover all the overloadings' amplitudes. However, because of the time mismatch between the cars' availability and the need for balancing, only some cables could be completely relieved. This means that only the buildings connected to these cables could install residential HPs.

This conclusion is strongly dependent on the minimum SoC set for discharging. The study was conducted with a minimum SoC value of 80%. It was further shown that lowering this parameter to 70% could help solving the overall overloading problem. However, this is expected to have an impact on the following charging requirements, which should be further studied.

As future research, a model predictive control logic should be investigated in order to solve the time mismatch challenge. Furthermore, the validity of the results should be tested against potential daily EVs' driving patterns variations. A sensitivity analysis regarding the penetration level of EVs and their SoC boundaries for charging is also suggested.

Finally, the approach used in this study shows that city planners, energy utilities, and engaged citizen can benefit from taking a district-to-city level perspective on integrated energy systems.

**Author Contributions:** Conceptualization, M.A., and M.T.; methodology, M.A.; software, M.A. and M.T.; validation, M.A. and M.T.; formal analysis, M.A.; writing—original draft preparation, M.A.; writing—review and editing, M.T., M.A. and B.L.; visualization, M.A.; supervision, M.T. and B.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
