Implementing flexibility into energy planning models: Soft-linking of a high-level energy planning model and a short-term operational model

The operation of electric and heat grids alike is complicated due to the dynamic demand, with the increasing penetration of renewable energy sources adding to the problem. In order to improve the integration of variable renewable energy sources, the flexibility of the system needs to be improved. This paper proposed a novel characterization of the short-term energy flexibility, which was further utilized for the district heating capacity extension. The soft-linking of the models includes feedback, but the added computational complexity is kept at a minimum. Compared to the other literature in the field, due to the accurate characterization of the dynamics of the energy flexibility, flexibility is utilized much more frequently. The method was demonstrated for the case of the district heating of Zagreb. Results showed that both capital and operational savings can be achieved by adopting the proposed method. In the best performing scenario, which included the capacity extension planning, the savings of the district heating system were 5.4%. The extensive power exchange in the best performing scenario meant that the flexibility was used to help balancing the power grid as well.


Introduction
Energy systems in many countries across the globe have started transitioning towards sustainable energy supply. Worldwide, the share of renewable energy sources in total energy generation was 25% in 2017 [1]. Hydropower is dominating the renewable energy capacity installed, having a share of 58% of the total installed renewable energy capacity [2]. However, the largest increase in the last decade occurred in the installed capacity of photovoltaics (PV) and wind turbines, both being variable renewable energy generators. Between 2010 and 2017, the installed capacity of wind turbines increased from 0.1 TW to 0.51 TW, while installed PV capacity increased from 0.015 TW to 0.39 TW (an increase of more than 20 times) [2]. Future increase in renewable energy generation capacity is expected to be heavily dominated by variable renewable energy sources [1]. Some countries have already achieved large shares of variable renewable energy generation in total power demand. Denmark is one example, managing to meet more than 40% of total electricity demand by wind generation already in 2015 [3]. Moreover, power demand will significantly increase in the future. According to the International Energy Agency, electricity demand will increase by 40% by 2040 under the business-as-usual scenario [1].
Variable renewable energy sources can be successfully integrated into future energy systems if sufficient flexibility exists in the energy system. There are four main flexibility sources in an energy system import and export of energy over the system boundaries, energy storage, power-to-heat and power-to-gas-to-power technologies, as well as different demand-response mechanisms [4]. The latter flexibility source is the focus of this research paper.
Looking into final energy consumption by sectors in the European Union (EU), 38% of the consumption occurs in the building sector [5]. Contrary to the power sector, the inertia of the heating and cooling systems allows for a temporary mismatch between demand and supply, without causing severe consequences. Taking into account the latter possibility, and successfully implementing controls for the activation of heat generators based on different signals, it is possible to implement demand-response mechanisms in the systems that supply thermal energy to buildings. District heating can make up a large proportion of the building energy consumption, as it is the case in Scandinavia, Eastern and South-eastern Europe and other regions.
There are two main options in district heating systems to increase the flexibility without imposing additional capital investments. Those are to use the thermal mass of buildings or heat capacity of water flowing through a heat distribution grid [6]. This paper focused on developing methods for flexibly utilizing the heat capacity of water flowing though the heat distribution grid.
Although some papers, such as [7], focused on representing flexibility on the building side, there is a lack of papers assessing the potential for utilizing the flexibility of the district heating grid and building's thermal mass for integrated energy planning. One paper has focused on the electrification of the heat supply for an Irish case study [8]. The authors have concluded that by utilizing the inertia of building systems, heating costs could be reduced to the cost of the cheapest alternative technology (gas boilers); however, district heating systems were not considered [8]. Another group of authors have developed an add-on for linear system-level optimization model (Balmorel) that could take building flexibility into account [9]. However, they have used a static thermal heat transfer model of the buildings, resulting in an under-representation of the peak loads that occur after demand-side management events take place. Implementing more detailed building dynamics into holistic energy planning models, either linear optimization ones or simulation ones, can be done via soft-linking or hardlinking of models. Moreover, soft-linking can include feedback between the two models. Soft-linking of a detailed dynamical building simulation model and holistic system-level linear optimization model has been carried out in [10]. It was found that the flexible heat load accounted for 5.5%-7.7% of the total district heating demand [10]. However, soft-linking of the models has not implemented a feedback loop. Due to the latter, a building simulation part of the model calculates the flexibility potential for simulated conditions defined prior to the revelation of the real weather conditions, potentially resulting in inaccuracies of building flexibility representation [10]. Furthermore, the computational complexity of the latter approach does not allow for the implementation of controls in a real-time [10].
Furthermore, a research paper on the utilization of heat capacity of water inside distribution grids has pointed out that the corresponding storage capacity is much smaller compared to the thermal mass of buildings as energy storage [6]. The latter paper has also pointed out that the storage capacity per 1°C is approximately 1000 times larger for the case of the thermal mass of buildings than the storage capacity of water located in the district heating distribution grid. Usually, the temperature of the water inside the distribution grid can vary much more than the thermal mass of buildings; however, the resulting difference in the storage potential is still 100 times. Another research paper has carried out a simulation using the heat capacity of water flowing in the distribution district heating grid [11]. The control strategy consisted of varying the supply temperature on a typical winter day. A typical loading phase started around 2 AM and lasted for 3.5 h, in order to reduce the morning peak load [11]. The result showed that up to 15% of daily peak demand can be moved, increasing the distribution losses by about 0.3% [11]. A review paper on flexibility in district heating has mentioned the possibility of storing the energy by temporarily raising or lowering the temperature of water in the thermal network [12]. However, the authors have raised the concern, although without any results or references provided, that more frequent cycling could result in changing thermal stresses, possibly developing pipe cracks that could lead to pipe failures [12].
Although theoretically the thermal mass of buildings has much larger storage potential than the water in district heating systems, there are two challenges related to controlling the thermal building mass of buildings. The first problem arises from the estimation of the flexibility of many different building archetypes that can be found in a city. showed that detecting typical archetypes and finding proper heat transfer and ventilation models on city-scales is complicated both for households [10] and for commercial buildings [13]. The second problem concerns the controlling of all the different buildings in a district heating system. Usually, the central district heating operators do not control temperature set-points of buildings and changing this could be challenging. On the other hand, differential pressure and supply temperatures in district heating distribution grids are usually controlled by the central operators. Thus, the flexible operation of the distribution grid itself is much easier from the implementation point of view. Likewise, economies of scale make the water in the district heating systems much more attractive than the building mass, since just a few controllers are needed to utilize the district heating water as opposed to a controller for every single household needed to utilize the building mass. Thus, in order to improve the current state-of-the-art in research on district heating flexibility representation in holistic level energy models, and to develop a control approach that could be used in a real-time, soft-linking of energy flexibility and simulation-based energy planning is proposed in this paper. Moreover, the soft-linked models included the feedback from the operational model (flexibility function) to the District heating system model (Fig. 2). This paper presents a continuation of the work carried out in [14], in which energy flexibility was characterized dynamically. The latter research was improved and expanded in this research paper. Furthermore, compared to the approach used in [6], which uses a very small part of a district heating grid for storing heat in heat capacity of water, and [11], which uses a fixed control strategy for starting and stopping the increased load in district heating distribution grid at predefined hours, we utilize the district heating grid for storage much more frequently, but with lesser amplitudes (temperature changes of ± 3.5°C).
This paper improves the current state-of-the-art with the following three novelties: (1) Dynamic characterization of energy flexibility based on differential equations is proposed and used to represent the energy flexibility of district heating systems. (2) An energy system model that focuses on district heating systems is soft-linked with the energy flexibility of the district heating systems, with the objective of minimizing operational and/or total system costs of the district heating system. (3) The energy flexibility in district heating systems is utilized by using the heat capacity of the water in the district heating distribution system itself as storage. In this way, there is no need for complicated, data-intensive and time-consuming representations of different building archetypes for a whole city.
In the following section, the district heating model and the model predictive control are presented, as well as the linking between them. Section 2 consists of three different subsections. The first two subsections are related to the explanation of each of the two soft-linked models. The third subsection presents the developed method for the soft-linking of the models. Section 2 is followed by Section 3, in which the district heating system of the city of Zagreb is described. Sections 4 and 5 economic savings achieved by the model coupling, discuss the reasons for the encountered behaviour and point to the possible future research pathways. Moreover, Section 4 shows the impact of implementing flexibility on both operational, as well as capacity extension planning. The main qualitative and quantitative findings can be found in Section 6.

The energy model of the district heating system
A linear optimization district heating model was developed specifically for this research question in order to make it straightforward to soft-link it with the flexibility function. The hourly time resolution was used and the model was created in Python using Gurobi. Python is freeware while a Gurobi license can be obtained for free for academic purposes. The district heating model was based on the hourly heat energy demand, and energy equations were used instead of temperature and/or flow values in order to keep the problem linear. Keeping the problem linear is useful for capacity extension problems when many different investment possibilities exist. The grid layout was not explicitly modelled; it was rather taken into account via heat losses, which were obtained from the district heating operator.
The technologies predefined in the model were gas combined heat and power (CHP) plants, gas boilers, electric boilers, electric heat pumps and a heat accumulator (TES). The model allows for capacity expansion planning, or for operational optimization if the capacities of the plants are bounded in the model. The interaction between heat and power markets is modelled for CHPs, heat pumps and electric boilers. Revenues from electricity sales on day-ahead power market were modelled as income in the district heating system, while the total carbon emission, fuel and operating and maintenance (O&M) costs were all modelled as expenditures of the district heating system. Electricity needed for driving electric boilers and electric heat pumps was priced according to the hourly prices achieved on the day-ahead power market, on top of the fixed distribution and transmission fees.
The objective function consists of the total district heating socioeconomic costs which included fuel costs (including electricity costs), fixed and variable O&M costs, annualized investment and carbon costs, while electricity sales from CHPs were considered as revenue of the system (1). By treating electricity revenues as an income in the district heating system, all the fuel and emission costs were also assigned to the district heating system, avoiding potential uncertainties regarding the assignment of costs to the power and district heating sectors separately.
The model can handle geographically distributed district heating systems. In this case, equations are represented for two separated district heating systems. Variables in the model were heat generation power in each of the district heating systems for each hour, TES charging, discharging and a level of heat stored in TES. For the initial run of the district heating system model, a perfect forecast was assumed and all the hourly heat demand values were known in advance.
Different equality and inequality constraints were used in the model. All the generation from plants within a certain district heating system needed to match the heat demand in each hour (2) and (3).
where = H 1 h. The capacity of each plant needed to be larger than or equal to the generation values for each hour throughout the year (4) and (5).
Storage level in each hour needed to be equal to the storage level of the previous hour plus storage charge minus storage discharge, as described in (6).
All the heat generation variables and storage charge and discharge were related to the corresponding efficiencies, which due to simplicity were not presented in the equations above.
The district heating linear optimization model was soft-linked with the flexibility function, which description can be found in the following subchapter.

Description of flexibility representation
The description of energy flexibility used here relies on the assumption that the temperature of all water in the district heating network increases whenever the heat generated by the heat plants exceeds the heating requirements of consumers, and vice versa when it is surpassed by the demand. This is modelled by the simple ordinary differential equation (ODE): Here X t is the water temperature at time t scaled between 0 and 1, so that = X 0/1 t means that the water temperature is at the lowest/ highest acceptable level at time t. Y t is the heat generation, while B t is the initial heat demand (before accounting for demand response). Thus, C is the total energy required to change the temperature from the lowest acceptable level to the highest acceptable level. The expected energy demand, Y t , is modelled by the equations where logit is the logistic function (1 exp( )) , that maps real values to the interval (0, 1) and 1 is the indicator function for which means that the consumption is reduced as much as possible, and = 1 t means that it is increased as much as possible. This change is governed by two terms, the first describing the desire to be close to normal operation, defined as having = X µ t . When this is not the case, µ X ( ) t will change the consumption in order to make X t approach µ. The rate at which this happens is decided by the size of greater than 0. The second term describes how price influences the consumption. Most of the time, this is approximately given by ku t , where u t is the price of energy (district heating in this case) normalized between 1 and 1 and k > 0 is the parameter describing how important price is for the flexibility. However, when the state of charge is close to 0 or 1, it cannot be decreased/increased any further, and the appetite for being flexible diminishes gradually, as the state of charge approaches the boundaries. This is described by g, which is a function for which = g (0) 0 and = g (1) 1. An example of this is , in this case, the system is twice as sensitive to low prices < u ( 0) when the state of charge is e.g. 0.8 compared to 0.4. Usually, it is more realistic to assume that the system is equally sensitive to the price as long as the state of charge is not close to the boundaries. In this paper = g x x ( ) 1 exp( ), which for large values of quickly approaches 1 as x increases.
Finally, the expected generation Y t , normalized between 0 and 1, is given by (9) as the baseline demand, B t , modified according to t . Since the consumption is non-negative, and cannot surpass the maximum consumption, the potential increase is scaled by the difference between the baseline consumption and maximum consumption, B 1 t . Similarly, potential decreases are scaled by B t . It is also reasonable to expect that there is a limited flexible capacity. This is given by , where for example = 0.6 means that, if flexibility is available consumption can be decreased by 60% or 60% of unused capacity can be switched on.
Each of the parameters controls a particular part of the flexibility, which is illustrated for different parameter values. For the examples shown here, the parameter values shown in the description of Fig. 1 are used unless otherwise stated. Fig. 1(a) shows the change in demand for block increases in price with varying amplitudes. Unlike linear flexibility functions, the amplitude of the response is not linearly related to the amplitude of the price. In fact, the amplitude of the response for a step increase of 1 is only slightly larger than an increase of 0.5, but the larger increases in price sustain the large response for longer, until the flexibility is exhausted, after which the response quickly drops off. For the smaller increases in price, the response is smaller and quickly decreases, but since this means that the flexibility is not exhausted, the response is sustained to a small extent for a long time. Nevertheless, a larger change in price means that a larger part of the energy flexibility is used. This can be observed when the price drops down to the initial value, where the increase in demand reflects the state of charge. The case with the largest price also has the largest amplitude in the rebound, since in this case the flexibility has almost been depleted, while there is more left in the other cases. Similarly, Fig. 1(b) shows the change in demand for a fixed block increase of price and varying values of , the maximum change in demand. It can be seen how the initial change in demand is proportional to . However, as the total amount of available flexibility is independent of , this flexibility is simply spent faster for larger values of . The rebound is the same, where larger values of allows for faster recharge of the state of charge. In Fig. 1(c) the energy capacity, C , is varied. In this case, the initial change in demand is equal for all cases, but as the flexibility is depleted faster for smaller values of C , the demand returns to the baseline faster in these cases. Similarly, it takes more time to recharge the flexibility for larger values of C . Finally, Fig. 1(d) shows the effect of changing the propensity for utilizing flexibility, by varying . This parameter controls how inclined the system is to have a state of charge close to = µ 0.5 in this case. Consequently, the initial change in demand is equal for all cases, but as the state of charge decreases, those cases where is larger returns to the baseline demand faster. Similarly, larger values of results in a more pronounced rebound effect, since the state of charge, X t , is brought back to µ faster.

Soft-linking of the models
Schematic representation of the district heating planning and operational models can be seen in Fig. 2. Initially, the district heating model was run in order to detect the hourly shadow prices in the district heating system, i.e. the price of generating one more unit of district heat. The shadow prices were then used as an input to the flexibility function, which output was new hourly district heat demand (Flexible Demand in Fig. 2). The initial district heating planning model was then updated with new heat demand time series. The second run of the district heating model was the one that generated the final results.
One should note here that q n t , from (2) and q s t , from (3) are equal to B t in the original run but equal to Y t (7) in the second run of the District heating system model. Besides the shadow price of the district heating system, those variables were used to interlink the models.

District heating system used in the case study
A district heating system located in the city of Zagreb, Croatia was chosen for a case study. Zagreb has 800,000 inhabitants. The district heating system consists of two separated grids, one located in the northwest, and the other in the south-east. Supply plants are locally called EL-TO Zagreb (north-west) and TE-TO Zagreb (south-east). Both of those locations consist of several different heat and/or heat and power generators. Approximately 28% of households are connected to one of the district heating systems [15]. All the currently operating district heating plants are driven by natural gas. The complete overview of the plants and their technical parameters can be seen in Table 1.
Economic parameters of the heat generation plants used in this paper can be seen in Table 2.
In Table 3 the annual heat generation and demand for the district heating networks are shown and it can be noted that generation and demand have been steady in the last three years. Zagreb District heating operator provided the info presented in Table 3, as well as the hourly pattern of heat generation [19]. The hourly pattern of heat generation was not allowed to be publicly shared and that is the reason why it is not presented here. Yearly pattern and total heat generation for the year 2017 has been used in the case study carried out in this paper. Furthermore, the Zagreb District heating operator provided the information that the total volume of water in the distribution grid equals to 48,800 m 3 .
Although Croatia started its power day-ahead market in February 2016, called Cropex, and its data has been used for the hourly dayahead electricity prices [20]. The latter also means that the methods presented in this paper are transferable to any country that has developed day-ahead power markets. The average electricity price on Cropex day-ahead power market for the year 2017 was 52.04 EUR/MWh (with a standard deviation of 24.05 EUR/MWh).
As gas prices for non-household consumers were very volatile between different years in the period 2013-2018, an average between those years for non-household consumers for the case of Croatia was used in the case study, which equalled to 34 EUR/MWh of gas [21]. The carbon price was assumed to be 50 €/t. Opposite to the business-economic approach, the socio-economic approach used in this paper does not include business taxes in the model, as those are only internal transfers within the society [22]. However, negative externalities such as climate change effects should be internalized in the socio-economic approach and thus, the CO2 emission price was adopted in the model. In the time of writing of this paper, transmission and distribution tariffs for non-household consumers were complex, and several different tariffs existed in Croatia. However, the difference between different tariffs was not large, so the average price of 20 €/MWh of consumed electricity was used for the transmission tariff and a price of 20 €/MWh was also used for the distribution tariff (totalling at 40 €/MWh of electricity consumed). The latter numbers were added on top of the day-ahead electricity price for electric boilers and electric heat pumps, in scenarios where those two technologies were installed.

Flexibility parameters
Since the district heating networks are not currently operating with price-responsive control it is not possible to estimate the parameters from data. However, qualified guesses can be made, based on the physical properties of the district heating network. First of all, the total amount of water is approximately 48,800 m 3 , amounting to the heat capacity of 56.8 MWh/K. In this paper, it is assumed that a change of ± 3.5 K in the temperature of the district heating water can be handled by the local heat exchangers receiving the district heating water. The parameter C is scaled by the maximum demand, 1057 MWh. This means that the total range of allowed temperatures is 7 K, which is obtained by having = C 7·56.8/1057 0.4. It is expected that all of the heat generators can be flexible, in which case = 1. The parameter was set to 4, a medium sized value, as it is preferable to have the temperature of the district heating water close to the business as usual temperature, but at the same time the exact temperature is not important. Similarly, there are no fatal consequences of having the temperature leave the temperature interval, and so = 50, which is a large value that allows the temperature to get very close to the boundaries before the price is ignored to keep the temperature within the boundaries. The price responsiveness is controlled by k, where large values makes the system more responsive to price changes. It was chosen to have = k 5, so that the largest possible effect of price is 5, when the price is −1 or 1 which should be compared to = 4, which means that the largest possible effect of the state of charge is = = µ 4·0.5 2, when the state of charge is 0 or 1. The largest Having the largest effect of price 2.5 times larger than the effect induced by the state of charge, means that prices are expected to be 2.5 times more important than efficiency losses, which seems reasonable. Of course, without data to support the parameter estimates, they are nothing but qualified guesses, and thus a sensitivity analysis was carried out to analyze the effect of the parameter estimates on the objective value of the energy planning model. For , values between 10 and 100 were tested with effects below * Although it is a back pressure cogeneration plant, the DH operator reported that it is malfunctioning and thus, only operating in the heat generation mode. ** Capacities installed in the Ele_DH scenario. No capacity installed in the Basic scenario. In the Cap_Ext scenario, the installed capacities of electric boilers and heat pumps were a result of the optimization. If needed, consult Table 4 for the description of the scenarios.

Table 2
Economic parameters of plants used in the case study [18].   Table 4 Description of three different scenarios used in the paper.

Basic Electrified District Heating (Ele_DH) Capacity extension (Cap_Ext)
All the heating plants data and inputs were adopted as described in the District heating system used in the case study subsection ( 0.05%, while values between 1 and 10 were tested for both k and with effects of at most 0.25% and 0.01%, respectively.

Description of scenarios
Three different scenarios were carried out in the paper, as described in Table 4. The Electrified District Heating (Ele_DH) scenario represented a move to transition the heat sector towards electrification, in order to reduce CO2 emissions and help integrating more variable renewable energy sources in the power sector. For this specific case study, the electric boilers and the heat pumps replaced gas plants older than the year 2000.
The Capacity Extension (Cap_Ext) scenario was created to present the possibility of utilizing the flexibility function when carrying out a capacity extension planning. A potential district heating grid extension was not modelled, as there are currently no plans for doing that in the region. The latter possibility could be added to the model if there would be available input on the heat demand density of the potentially newly connected districts. The Zagreb district heating operator assumes that newly built buildings within the area of the current district heating grid will be connected to it and that this new demand will cover a reduced demand originating from the improved energy efficiency of buildings, keeping the overall district heating demand flat. The Cap_Ext scenario showed the potential for reducing the needed capacity when expanding the district heating system in comparison to the original expansion plan.
In the Cap_Ext scenario, out of the original plants operating in the DH system, only newly built boilers and CHP units were left. More specifically, DH north had installed two CHP generators with a heat capacity of 10.25 MW each and one gas boiler of 116 MW. DH south had two CHP plants with a heat generation capacity of 80 MW each, as well as the heat accumulator. The model allowed new investments in electric boilers, heat pumps, gas boilers, CHP units and the extension of the heat accumulator capacity. It is worth noting that the initial capacity of plants in this scenario was insufficient, so the model had to invest in the additional capacity of heat generation technologies.
The original expansion plan corresponded to the first run of the district heating optimization model, while the second run of the district heating optimization model (after the flexibility function implementation) corresponded to the updated expansion plan (if needed, consult the scheme presented in Fig. 2 for the sake of clarity).
In all the runs of all the scenarios, the same amount of heat was needed in each hour.

Results
The results section starts with the overview of the installed capacities, followed by the presentation of differences between different scenarios in terms of generated energy from specific sources. The results section further continues with the presentation of the overall economic results of different scenarios and finally, it ends up with a detailed representation of operational changes during the chosen two days. The representation of the detailed operational changes explicitly shows the influence of the operational model on the energy planning model.
The resulting capacities of heat generation plants can be found in Table 5. As described earlier, in Ele_DH scenario, two back pressure units that were operating only as gas boilers were replaced with an electric boiler and a heat pump in the DH north. In DH south, four gasdriven heat boilers were replaced by a heat pump and an electric boiler with a capacity of 135 MW each. Comparing the results of the latter scenario in the original run and a run with flexibility included (last two rows in Table 5), one can notice that by using the flexibility function, needed extension capacities of gas boilers and electric boilers dropped in the case of DH north. In the case of DH south, the needed capacity of TES storage reduced by 2.5%.
An example of changes in plant operation in DH north in the Ele_DH scenario, with and without flexibility, can be seen in Fig. 3. One can notice that the gas boilers with flexibility reduced the operation time under the high capacity and shifted that generation to the periods of lower operation capacity. On the other hand, heat pumps increased their operation at the maximum rated capacity, from 1811 h to 2904 h. CHPs also slightly increased their operation under the maximum rated capacity, from 4885 to 4943 h.
Heat generation from different plants in different scenarios can be seen in Table 6. In the Basic scenario, flexibility utilization increased the operation of CHP plants and reduced the operation of gas boilers, reducing the overall socio-economic costs of the DH system.
In the Ele_DH scenario, the utilization of the flexibility increased the use of the CHPs, electric boilers and the heat pumps, while it decreased the use of the gas boilers. Reduced operation of the gas boilers was especially emphasized, as they generated 16.7% less heat in the case with the flexibility.
In the Cap_Ext scenario, the operation of CHP units increased significantly, while the use of gas and electric boilers was reduced correspondingly. Moreover, the operation of heat pumps also reduced in this case.
Summarizing the heat generation results as seen in Table 6, one can note that utilizing the flexibility, the CHP generation increased in all the scenarios, between 0.3% and 3.2%. On the other hand, gas boilers operation reduced in all the scenarios following the implementation of the flexibility. Finally, heat pumps and electric boilers increased their generation in the Ele_DH scenario, while they decreased their generation in the Cap_Ext scenario. The reason for the latter behaviour is that the Cap_Ext scenario had much larger capacities of CHPs than in the Ele_DH scenario, and this additional CHP capacity took over some of the generation of electric boilers and heat pumps.
Economic results of different scenarios are presented in Table 7. One can note that DH north achieved significantly larger savings compare to the DH south, which was the case in all the scenarios. In the basic scenario, the socio-economic savings were relatively low. This scenario was dominated by gas boilers and gas CHPs, without any heat generation plants driven by electricity. The latter shows that in gas dominated district heating system, steady intra-year gas prices reduce the possibility to utilize price arbitrage.
In the Ele_DH scenario, savings were more significant, especially in the case of the DH north. Total socio-economic costs reduced by 2.8% in the DH north. For the case of DH south, the total socio-economic costs of the DH system reduced by 0.6%.
The largest savings occurred in the capacity extension planning scenario, i.e. the Cap_Ext scenario. The objective value of the DH south improved only marginally; however, in the DH north, the savings amounted to the significant 15.5%. The most dominant factor in the improved performance of the DH north was the larger income achieved from electricity sales, i.e. revenues from selling electricity increased by 7.5%.
In the Cap_Ext scenario, the economic results showed that the

Table 6
Heat generation from different plants in three different scenarios, including the relative difference before and after utilizing flexibility.
Total CHP production (GWh) Total gas boilers production (GWh) Total electric boilers production (GWh) Total heat pump production (GWh) needed extension capacity can be reduced if the flexibility function is utilized before the decision on the capacity extension is made. The change in capacities, as seen in the last two columns in Table 5, resulted in absolute savings of EUR 0.55 million. The adoption of the flexibility resulted in a reduced need for 1 MW of gas boilers, marginally reduced capacity of electric boilers, marginally increased capacity of heat pumps, and reduced capacity of heat accumulator for 181 MW. One should note that the latter cost saving is not directly observable in Table 7, as results presented in Table 7 include only annualized investment costs and not the total investment costs. The lowest overall cost of the energy system was achieved in the Cap_Ext scenario, which was expected as there was no excess capacity in the district heating system. The socio-economic cost of the district heating system in the Cap_Ext with flexibility scenario was 30% lower than the Ele_DH with flexibility scenario. It is worth explaining here that the negative values for the Operational costs with electricity sales income that occurred in the DH north in the Cap_Ext were the consequence of electricity income being larger than the operational costs of the district heating system.
In Figs. 4 and 5 the heat production before and after applying the flexibility function is shown for the Southern and Northern grid respectively. The discharge of the heat accumulators is shown such that negative values correspond to filling the storage. The shadow price associated with marginal changes in consumption and the electricity price is shown as well. For DH South the most noticeable change is during the morning of 11-Jan where the electricity prices are very high and the peak demand is shaved so that the electric boiler is not needed, and the heat pumps are used less. Apart from this, the flexibility is mostly being used to change when the heat storage is being charged, so that this happens mostly when the marginal cost is low.
Generally, CHP's were the cheapest source of heat due to high electricity prices in this period, while heat pumps and electric boilers were used for the peak load. An exception was during the night between 11-Jan and 12-Jan, when the electricity price was relatively low, resulting in the heat pumps being cheaper than the CHP's. For DH North the effect of the flexibility is again seen to shave the peaks where the electric boiler was used during the morning of 11-Jan. This was accomplished mainly by increasing the demand during the night hours, where demand was low, as seen during the night between 11-Jan and 12-Jan. It is also clear how the high electricity prices during this period make the gas boilers cheaper than the heat pumps during the day of 11-Jan, 09:00 12-Jan and 20:00 12-Jan.
In Fig. 6, the accumulated savings for each of the district heating grids are shown. As already mentioned, the largest savings are obtained for DH North even though the overall demand is larger for DH South. As the heat demand is larger during the winter, in the heating season, it is no surprise that this is also where most of the total savings occur. Especially for DH South, there are almost no savings at all outside of the heating season, except for one event at the start of August. Except for that event, the savings due to the energy flexibility are steady, with a continuously varying slope that peaks during the coldest months; January, February and December.

Table 7
Economic results from different scenarios (million euros), including the relative difference before and after utilizing flexibility. * In the objective function of the optimization, the goal was to minimize the total costs of the district heating system, which included revenues from selling electricity as an income (the objective function is represented in (1)).

Discussion
The discussion section starts with a brief recap of the main economic results, followed by a discussion on the reasons for the achieved savings. It continues with a discussion on the influence of flexibility utilization on the optimal operation of heat accumulators. Furthermore, the section continues with the discussion on the operational constraints in the system, as well as a discussion on the possibility for an automatic parametrization of the flexibility function. Limitations of this study, comparison with the existing literature and possible future research directions are discussed at the end of the section.
In this paper, energy planning and operational models were soft coupled in order to assess the potential of utilizing water stored in district heating grids as flexible thermal energy storage. The case study carried out to analyze the methods of this paper confirmed the benefits of this approach. All three scenarios showed that significant socioeconomic savings can be achieved by flexible operation of the heat generation in the district heating systems. The total socio-economic savings ranged from 0.4 to 3.1 MEUR/year. In relative terms, economic savings were between 0.4% and 5.4%. The mentioned savings are achieved in the district heating sector. However, feeding more electricity when the electricity prices are high (due to high electricity demand or low availability of electricity supply), and feeding less electricity to the power grid when the electricity prices are low (due to oversupply of cheap electricity or low electricity demand) helped balancing the power system. Nonetheless, quantitatively describing those implicit savings was out of the scope of this paper.
Most of the economic savings occurred from the combination of lower operating costs and/or higher revenues from electricity sales when the flexibility was utilized. The behaviour was different in the DH north and DH south in different scenarios. For the case of Cap_Ext scenario, increased revenues from electricity sales were by far the largest contributor towards reduced overall socio-economic costs in the DH North, while in the DH south, the largest contributor towards the overall socio-economic savings originated from reduced operating costs. The behaviour described here points to the conclusion that both reduction in operating costs, as well as an increase in revenues from better-performing electricity sales, contribute towards the overall economic improvement reached by utilization of the flexibility. This means that the biggest advantage of energy flexibility in district heating grids  is that it supports the power grids.
It is also interesting to take a closer look at the behaviour of the heat accumulator. Following the implementation of flexibility, the utilization of the heat accumulator decreased by 12%, 12% and 19% in the three scenarios, respectively. Furthermore, the reduced installed capacity of the heat accumulator in the Cap_Ext scenario upon the implementation of the flexibility indicates competition between the different flexibility options. The final point that goes along the same line is that the operational savings in the DH South were always lower than in the DH North, which was especially emphasized in the Cap_Ext scenario. In the latter scenario, after the system optimized the initial size of the thermal accumulator, only slight savings were achieved by including flexibility. The competition between different flexibility sources was already tackled and discussed in the literature in a detailed way [23], which needs to be taken into account when discussing the potential of a single source of flexibility.
One should note that by adopting the coefficient C = 0.4 in this paper, we are varying temperatures only up to ± 3.5 K from the original temperature values. In this way, neither a significant additional thermal stress on the piping nor a complicated control action demand on substations would be imposed. Current substations in district heating are already oversized, as they are designed for the critical conditions. Furthermore, the current controllers on the end-user side should already be able to adapt the flow to accommodate the slight temperature changes on the primary side of the system. Due to the latter, the approach used in this paper should be easy to implement.
Contrary to the approach used here, the parameters of the flexibility function should be estimated based on data, when possible. Indeed, if a district heating company would implement the control presented here, the data required for estimating the parameters would become available from their daily operation, and the parameters could adaptively be updated to reflect the flexibility [24]. Furthermore, with enough data available the representation could be improved, using appropriate machine learning techniques.
One of the strengths of the approach used in this paper is the easy implementation of the flexibility function with automatic estimation of the parameters. This implies that the same approach is valid for district heating systems with different characteristics. Also, following the digitalization of the district heating sector that has started gaining a stronger foothold, data-driven solutions will be easy to implement in many different district heating systems. Furthermore, this paper showed that the implemented soft-coupling of models with feedback can be used both for operation and capacity extension planning. In addition, the approach adopted in this study requires only slight temperature variations compared to the base operation, which could be implemented solely by central district heating operators, while the enduser substations should be able to adapt automatically to the changing conditions.
Limitations of the study include the potential reluctance of the centralized district heating operators to change the business-as-usual operation. Furthermore, district heating systems that have a direct exchange of heat, which can still be found in Eastern and Southeastern Europe, could have problems even adopting to the slight temperature variations. Moreover, the approach used in this case study is tailored for the district heating systems that have constant inertia, which means that flow velocities should be kept constant. The latter is usually only the case for winter operation of the district heating systems when flows are kept constant and temperatures are varying. Fortunately, this is exactly where most of the savings were obtained, so the inaccuracy related to non-constant flows should be very limited.
Comparing the results of this study with the literature, one can notice that our implementation of controls is novel compared to the other studies and no direct comparison can be obtained. In [11], the authors simulated a district heating system with the temperature control implemented such that the district heating network was preloaded (between 2 AM and 5:30 AM) in order to reduce the morning peak consumption that occurred at 6 AM in their simulated case [11]. Their temperature limiter was significantly alternated, from maximum supply temperature of 95°C to maximum supply temperature of 115°C. The authors achieved a peak load reduction of up to 15% [11]. In our case, the temperature variation in comparison to the original case was only ± 3.5°C and heat demand pattern was alternated very frequently based on the provided price signal.
In [6], a similar step function system was implemented; however, with preheating of the grid. The implementation of the increase in artificial penalty signal was fixed in the middle of the considered week. The authors concluded that the flexibility of the thermal mass of buildings is two orders of magnitude larger than the heat capacity of water located in the district heating grid [6]. However, having a completely different implementation of the penalty signal, and not calculating consequences on the economic performance of a district heating supplier makes it incompatible for comparison with the case study carried out in this paper.
Up to the knowledge of the authors, this is the first paper that focused on the implementation of the flexibility control, using economic performance as a penalty signal. The results of the paper showed that this set-up can be used both for operational planning, as well as for capacity extension planning. Moreover, penalty signals could easily be tailored in a way that it targets heat loss reduction or CO2 emissions reduction.
Furthermore, the approach taken in this research opens up several possible avenues for future research. One option is to adopt the same approach towards the representation of the flexibility of buildings. For this, the exact same model could be used; however, obtaining enough data for parameterization of the flexibility function of buildings could be challenging. Another problem that could arise with the demand-response utilization in buildings is the control of the substations and the individual thermostatic valves within the apartments, as those are usually not under the central control of district heating operators [12]. On the positive side, the potential thermal mass as storage of buildings is much larger than that of water in district heating pipes [6]. The second option is to apply the same approach to other energy markets, i.e. intra-day and balancing markets. The latter could be done with slight changes compared to the approach taken in this study, i.e. the one based on the day-ahead market operation. The third option could be to show that by tailoring the flexibility function in a different way, the soft-linked models could directly reduce peak demands in district heating systems over the year. The latter approach is different compared to the solution of this study that was minimizing the socio-economic costs of the system, which also included increased revenues from electricity sales.
Some further research possibilities would be to test the solution in different set-ups, such as low-temperature district heating grids, as well as district heating systems of varying sizes. Furthermore, the possibility to test the solution in a detailed dynamic simulation model would be beneficial. It would also be of interest to adopt the coupled model to district heating systems that have more variable heat generation sources in its portfolio.

Conclusions
This paper presented soft coupling of the district heating system optimization model, which took the electricity system into account via day-ahead electricity prices, with the flexibility function that represented the heat demand response following a penalty signal. The penalty signal was obtained by computing the shadow price of consuming heat for each time step.
The method of this study proposed slightly varying supply temperature in the distribution grid compared to the baseline operation, i.e. ± 3.5 K, in order to use the district heating distribution grid as storage. Results of the case study showed that significant socio-economic savings could be achieved. The savings ranged from 0.36 MEUR to 3.1 MEUR. The savings in relative terms ranged from 0.4% to 5.4%. Moreover, the scenario that included capacity extension possibility showed that the heat accumulator extension capacity could be reduced by 6%, still meeting all the district heating demand. Furthermore, the influence of different parameters needed for the flexibility function on the heat demand behaviour was shown in order to make it easier to comprehend how the function would behave in other district heating systems. It was further debated that the parameter estimation in the district heating grids could be estimated in an automated way as operational data becomes available, thus improving the functionality. The latter is possible as the approach used in this paper can be implemented in a centralized way from district heating operators.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.