The operation of district heating with heat pumps and thermal energy storage in a zero-emission scenario

With the decarbonisation of electricity generation, large scale heat pumps are becoming an increasingly viable prospect for district heating installations. Heat pumps couple heat demands to an intermittent electricity supply with varying electricity prices with the use of thermal energy storage providing flexibility to avoid peak electricity charges and minimise operating costs. However, the operating strategy for cost minimising in district heating system models is dependent on the size of heat pump and thermal energy storage capacity chosen and its operational conditions. Model predictive control techniques can be used to explore district heating configurations with varying forecast horizons. This study applies optimisation to a district heating operation model simulation to find low cost combinations of heat pump and thermal energy storage sizes. Physics-based representations of a district heating network and thermal energy storage are developed with ground source heat pumps and applied to a district heat load profile with hourly marginal electricity costs derived from a modelled zero-carbon electricity system as a basis for operation. Using a dynamic programming algorithm with different forecast horizons to minimise operational costs, the total costs of combinations of heat pump and thermal energy storage sizes are calculated. The operation at smaller thermal store sizes shows cycling multiple times per day, while at larger sizes these sub-daily cycles are maintained but longer multi-day cycles become more predominant. It was found that thermal energy storage equivalent of around 1% of annual demand is sufficient to minimise operating costs and enables flexibility beyond 4 days. This has important consequences for the electricity system and can facilitate the integration of variable renewable electricity.


Introduction
The Climate Change Act (2008) set a target for the UK to reduce greenhouse gas emissions by 80% from 1990 baseline levels, this was later amended to set a net zero or 100% reduction [1]. The provision of space, water and industrial heat accounts for around 40% of all energy consumption and 40%-50% of greenhouse gas emissions in the UK. To meet the net zero target, emissions from heat and the building sector are required to be near zero [2]. The decarbonisation of electricity generation in the UK is facilitating a shift away from fossil based heat generation to electrified heating with heat pumps [3] in district heating (DH) and consumer systems. Heat pumps (HP) are likely to play a key role in the decarbonisation of heat and their integration into DH systems is a promising application [4,5]. Indeed, Part L of the building regulation now encourages developers to consider the use of heat pumps and connection to DH in order to not exceed target emission rates [6]. The use of thermal energy storage (TES) in DH in combination with heat pumps, offers a potential for the management electricity systems with variable renewables [7]. The market penetration of both heat pumps and DH in the UK is currently very low and there is thus a need for research and analysis on the design and modelling of DH in a decarbonised energy system [8,9]. This paper presents the implementation of a model predictive control (MPC) algorithm to simulate the operation of DH with large scale heat pumps in a decarbonised energy system.

Background
The optimisation of DH operation has been developed here to inform real-time control in actual systems. Optimal control often employs model predictive control (MPC) with dynamic process models, optimising over a finite time horizon. MPCs can be applied to real time operation with the use of feedback loops to adjust processes and predict how a system is likely to respond. They employ an algorithmic optimisation such as: linear programming, mixed integer linear programming (MILP), mixed integer nonlinear programming (MINLP), genetic algorithms and dynamic programming (DP) and are typically computed with commercially available solvers such as GAMS and CPLEX.
The choice of algorithm depends on the system functions being modelled. The processes and components in a DH system can be nonlinear. While nonlinear algorithms allow complex interactions to be modelled, linear algorithms are often faster and scale better but can be limited in applications and care must be used with the formulation as global optimum solutions are not always attainable or indeed known. Nonlinear processes and constraints need to be linearised or approximated when used with linear algorithms and this can affect the accuracy of results.
Much of the literature on operation and control has focused on the optimisation of CHP in DH, often in conjunction with storage. Comprehensive reviews on the use of MPC for DH modelling have been covered by Sameti and Haghighat [10] and Vandermeulen et al. [11]. The literature features many examples of single objective MILP optimisation for a variety of configurations. It has been applied to the scheduling of CHP with TES and boilers by Verrilli et al. [12] and for CHP utilising building thermal inertia by Leśko et al. [13] who use an iterative approach to approximate nonlinear interactions. Moustakidis et al. [14] applied a formulation that used multiple optimisation controls over different time-scales and was used with month-ahead renewable generation forecasts to minimise global costs [15].
Comparatively few studies have used nonlinear (NL) methods to control operation. The thermal inertia of a DH system and buildings was used in conjunction with CHP to improve the utilisation of renewable electricity minimise operational costs using a MINLP [16]. Powell et al. [17] present a MINLP to find optimal charge rates for TES with the objective of minimising costs when participating in the wholesale electricity market over a 24-hr time horizon. Other formulations have studied various storage configuration with electricity spot markets [18,19].
Despite the design problem being similar, there has been little analysis of HPs in a DH context, with control studies focusing on building systems [20]. The use of MPC with HPs in tandem with batteries has been analysed, with a focus on dynamic pricing and grid ancillary services [21,22]. The goal with this analysis is to control the use of HP and TES, controls that are essentially binary decision variables and suited to the use of a DP algorithm. DP is a method that allows a high degree of expression with the system formulation and has been shown to perform better than other NL and LP algorithms when applied to an energy storage system with variable electricity costs [23].

Methodology
The District Heating Model (DiHeM) by the authors consists of a simplified representation of DH (as shown in Fig. 1) with an exogenously set DH heat load of 50 GWh per year with a constant 12% distribution loss factor and electricity price signals based on a zero emission scenario from Siddiqui et al. [24]. All flow conditions and outputs are calculated at an hourly resolution using mass and energy balances. The TES is modelled as a pressure connected, stratified water tank. We assume the tank hot water temperature T h , is the same as the flow temperature and stratified into two layers, with the cold layer temperature, T c , dependent on the DH return temperature, T r . Thermal losses are modelled assuming steel construction and 300 mm of insulation [25]. The mass of hot water, M h,i , in the tank at any time determines the state of charge and the flow rate,Ṁ TES,i , measured in magnitude and direction, enable the state of charge to be calculated in the next timestep, M h,i+1 . Fourth generation DH systems anticipate flow temperatures of 70 • C and below with average return temperatures circa 20 • C [26].
This model assumes a fixed flow-supply temperature, T f , of 70 • C, which is considered on the upper limit of what can be considered 'low temperature DH'. We select a minimum return temperature of 15 • C and peak return temperature of 50 • C, calculated hourly using a linear relationship with the DH load. The return flow is then either recirculated after being reheated by the heat pump or stored in the TES cold layer.
While the time it takes for water to flow through the network is not instantaneous, based on typical flow velocities, it is assumed that the flow returns to the HP within the hour [27]. The thermal storage of the distribution network has not been modelled. However, much of the thermal capacity of the system is contained in the building fabric which has been accounted for in the DH load profile and using an appropriate diversity factor.
The possible operating state of the DHN's components over an hour is defined as one of three discrete modes: 1. The HP is off, there is no electrical input and the heat load is met entirely by discharge from the TES.
2. The HP is on and covers the entire heat load.
3. The HP is at full capacity and covers the entire heat load using the residual spare capacity to charge the TES. In the event where heat load is partially met as either the HP is of insufficient capacity or there is insufficient discharge capacity in the TES, a penalty equivalent to delivering the remaining heat load using an electric boiler (or HP with COP=1) is applied. In mode 1, mass flow rate of discharge from the TES is equal to the mass flow rate of the DH. The returning flow re-enters and mixes with the TES on the cold side. In mode 3 we assume that the HP heat output Q HP , is at the maximum HP capacity. The mass flow rate through the TES,Ṁ TES , is taken in the negative direction with inflow charging the TES and T m is the temperature of the mixture of the return flow and TES cold discharge at temperature T c .
We model the COP as a fraction η HP = 65% of the Carnot cycle efficiency, and assume a ramp rate of 20% of installed capacity per minute and a minimum output of 25% of capacity [28]. A heat exchanger between the HP condenser and the DH flow is required to raise its temperature from the HP heat exchanger inlet to T f . The heat pump condenser temperature T cond is a constant 75 • C and must then be higher than T f to allow for losses across the heat exchanger which we assume is a counterflow heat exchanger. The log-mean temperature across the heat exchanger is then used to calculate the COP: We assume a ground source heat pump, though the model is easily adapted to other sources such as air, water or sewage. Ground temperatures exhibit less locational variance and small seasonal variance compared to air and other sources. Beardsmore and Cull [29] have given the calculation of temperatures at depth from periodic surface heating. Temperatures were calculated at an assumed depth of 3 m and the thermal properties of soil were estimated as per the data in Busby [30]. This is then used to calculate HP COPs. The COPs show seasonal variance as well as hourly variance according to T r which is determined by DH load. The mean modelled COP across all weather is 4.3 and corresponds well when compared to surveyed operational and modelled HPs in this temperature range [4,31,32].

Operating algorithm
The dynamic optimisation algorithm (DOA) is adapted from an open source operational optimisation library -Prodyn [33] which utilises a DP algorithm to determine the optimal control sequence given an objective function.
Atabay [34] describes how a model function such as the DH system must be discretised when using DP. The DOA system model is applied to DiHeM with discrete decision variables, U i , comprising of the three operating modes for each time step to calculate the state of the TES, X i , at the next timestep, X i+1 = f(X i , U i ). The TES is divided into discrete levels, X, of size 1 MWh for all store sizes, so the number of levels depends on the size of TES. The algorithm computes the cost J i of going from X i to X i+1 when a decision U i is made at a timestep for each possible decision variable and TES, then state J i = g(X i , U i ). A penalty cost is applied if the hourly heat load is not met, equivalent to the electricity cost of a HP with a COP of 1. The algorithm then works forward from a defined initial state X 0 , over N timesteps to finds the sequence, v, which minimises total costs, J, over all timestep in the time horizon.
A specific starting and ending charge state of X must be set. If the ending charge state is left undefined, the algorithm will always choose to empty the storage at the end as the optimum solution. We assume that a short term (up to seven day) DH demand look ahead is possible with reasonably high accuracy and beyond this there would be a good indication of relative conditions [35]; [36]. Similarly, while day ahead electricity system prices are available, with the demand and generation forecast short term electricity cost projections can be made [37] though of course this will depend on the configuration and operation a future electricity system. The hourly operation over which the DOA lookahead is used is then be restricted to 5 days (120 h).
With a small TES, knowledge of future operating conditions beyond a few days is of little use. With a large TES, projections of conditions a week or further in advance are desirable to optimally utilise the capacity. The required rolling DOA (RDOA) lookahead duration is estimated based on the HP and TES capacity. The duration of storage in the TES, D, can be estimated from the mean winter loadL, D TES = TES max /L. The TES level after a given number of hours cannot be more than the HP is able to attain in that time. The average time to recharge the TES from empty during the winter season, D CHR , can be estimated by D CHR = TES max / (HP max −L). If D TES is less than D CHR , the lookahead time horizon is then equal to D TES otherwise it is D CHR .
To define the final charge state X N , we assume that the final charge state must be sufficient to cover high cost demand in a projected period equal to D TES beyond the lookahead time horizon. Where the projected period is beyond the 7-day horizon of hourly forecasts, then the demand is extrapolated from the demand during the lookahead period.
After running the RDOA, the optimal control for the period D TES is returned but only the first 24 h are used, and the computation moves forward by one day to recalculate.

Results
All simulations presented here assume a HP capacity equal to peak heat demand (22 MW th ) to ensure security of supply and avoid the possibility of load not being met and penalty costs. The RDOA simulation is started during midsummer to the following summer with a starting TES = 0. TES levels using RDOA for various TES capacities are shown in Fig. 2 alongside the electricity costs representing marginal generation costs for a highly renewable, zero emission scenario.
The smallest TES size shown, 1 0.2% shows constant cycling throughout the year. The regular daily charge/ discharge frequency is maintained at 1% TES and the full capacity of the TES is frequently utilised in the heating season. Full TES capacity is sparsely utilised at capacities beyond 2% TES and the frequency dynamics decrease with larger TES sizes. The full charge/discharge cycles of the TES at larger sizes correspond directly to periods of high electricity prices.
With a HP capacity equal to peak demand, and 1% TES, charging the TES from empty to full takes between 1 and 2 days. With a 2% TES, this increases to 2 to 4 days. The doubling of TES capacity from 1% to 2% has little effect on the broader operating patterns in comparison to the operating pattern at smaller sizes. With a high forecast accuracy in the sub 2-day period, the RDOA operation of a 1% TES can be practically recreated in 'real world' conditions with a high level of confidence.  The TES operation can be split into sub daily cycles and multi day cycles. Fig. 3shows a frequency analysis of TES cycles with a 1% TES. The frequencies graph for 0.1% is similar to 1% except that the smallest cycle length is 1 cycle per day with none of the super-diurnal cycles shown in the inset. The TES never holds a full charge for longer than a day with a small TES, often cycling multiple times per day, in response to daily load variation. The 0.1% TES also exhibits a stronger twice daily cycle. The 1% TES also demonstrates these sub-diurnal cycles but also has longer multi-day cycles shown in Fig. 2 inset on a logarithmic scale with a 3-day and 7-day cycle showing strong signals.
Running the simulation for multiple annual periods with varying TES sizes reveals that electricity costs rapidly decline with increasing TES size then plateau between 1%-2% depending on the simulated period (average 1.3%). TES capacities beyond this are redundant in the modelled prices of the assumed high renewable system.

Conclusions
This paper has argued the need to study the impact of large-scale heat pumps in district heating systems, and their potential impact on the electricity system and consequently how they may be operated in conjunction with TES to provide flexibility. To analyse this an MPC algorithm using DP optimisation has been applied to a district heating system model which includes physics-based representations of heat pumps and TES. The system was simulated with variable electricity costs derived from a zero-emission scenario.
Minimum electricity costs being achieved at 1% TES suggests that if this capacity is adopted there is significant operational flexibility to be gained using DH with HP and TES. With the ability to shift demand by over 4 days at this TES capacity electricity peak loads can be reduced, and DH can facilitate the integration of variable renewable electricity.
Future work includes a simultaneous optimisation of the heat pump and TES to minimise capital expenditure and operating costs. An extension of this would be to analyse how the demand from district heating in turn impacts upon the electricity system and consequently prices in the case where HP-DH penetration is high enough to have a marginal impact on the electricity system.

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.