Operational optimization of district heating systems with temperature limited sources

Future district heating systems (DHS) will be supplied by renewable sources, most of which are limited in temperature and ﬂow rate. Therefore, operational optimization of DHS is required to maximize the use of renewable sources and minimize (fossil) peak loads. In this paper, we present a robust and fast model- predictive control approach to use the thermal mass of buildings as a daily storage without violating temperature constraints. The novelty of this paper includes two elements. First, the focus on an operational control strategy that explicitly accounts for temperature-limited renewable sources, like a geothermal source. Secondly, the optimization problem is formulated as a (nearly) convex optimization problem, which is required for adoption of model-predictive control in practice. The examples show that the peak heating demand can be reduced by 50%, if the thermal inertia of the buildings is used and the heating setpoints are adapted. Furthermore, the operational optimization ﬁnds the proper balance between beneﬁts of pre-heating using renewable sources with limited capacity and costs of additional heat losses due to pre-heating. 2020 The Authors. Elsevier B.V. creativecommons.org/licenses/by/4.0/).


Introduction
Scarcity of resources and climate change due to increasing level of CO 2 , has led to increased interest in renewable sources of energy. While wind turbines and solar panels have been in the limelight, relatively little attention has been paid to the energy used for heating and cooling. Heating and cooling account for 38% of the primary energy use in the Netherlands and around 50% in the European Union.
Using renewable energy sources for heating or cooling is subject to a spatial, temporal and temperature mismatch between production and demand. The spatial mismatch can be overcome by using a district heating/cooling network. A temporal mismatch can be overcome by making use of heat storage. This can be for example a seasonal storage like an aquifer thermal energy storage. Another option for short-term storage is preheating of buildings (loadshifting), or the installation of a heat buffer. It is a challenge to make optimal use of the different heat sources and storage possibilities. The temperature mismatch can be solved by (fossil) peak supply, heat pumps and by using short-term storage to limit the peak demand. Since most renewable sources are limited in peak supply and temperature, operational optimization techniques are required to maximize the use of renewable sources and minimize (fossil) peak supply.
The development of operational optimization techniques started in the 1990s. These papers focused on achieving a minimum acceptable supply temperature using the heat demand as a boundary condition to the optimization. Benonysson [1] has addressed many of the challenges associated with operational optimization of a district heating network, including time delays due to supply temperature variations, variable revenues from CHP production units and cost functions. Madsen et al [2] investigated stochastic prediction and control methods to minimize heat losses and CHP heat production costs, subject to constraints on the minimum required supply temperature for each consumer and a constraint on temporal temperature gradient. Dahl et al. [3] have investigated the impact of weather-forecast uncertainty and heat demand uncertainty on the operational supply temperature and associated savings. They conclude that the potential benefit of using dynamic uncertainties is larger for district heating systems with smaller pumping capacities.
More recently, model-predictive control (MPC) solutions were developed for a wide range of control problems in district heating systems. Faharani et al. [4] developed a MPC-based day-ahead pro-duction schedule for a smart thermal grid with CHP units and boilers connecting greenhouses with thermal buffer tanks. The MPC was a Mixed-Integer Linear Programming (MILP) formulation, taking into account capacity constraints, start-up costs, fuel costs, electricity sales with variable prices and heat demand uncertainties. Knudsen et al. [5] developed a MPC controller for demand response in domestic hot water preparation. This MPC was formulated as a non-convex continuous quadratic optimization problem with receding time horizon. Claessens en Vanhoudt et al. [6,7] developed a MPC strategy with 3 different short-term thermal storage configurations: central storage, distributed storage and thermal mass storage. The optimized CHP profit for the thermal mass concept was marginally smaller than for the distributed storage concept, but it is considered a promising demand side management strategy, due to limited implementation costs. The pragmatic MPC strategy is based on a three step approach, derived from electrical grid optimization with modifications.
Kensby et al. [8] have shown experimentally that a significant fraction of the daily heat demand on a winter day can be stored in the thermal mass of relatively heavy buildings without loss of thermal comfort. They provide general guidelines for the amount of heat that can be stored via direct load control without indoor temperature measurements. Lesko et al. [9] have developed an iterative approach using multiple MILP solutions to include nonlinearities in the operational optimization of a DHS using the thermal inertia of buildings as one of the short-term storage options in the DHS, although the impact of this thermal storage on the thermal building behaviour is not included in the optimization. Dominkovic [10] investigated the potential of thermal building mass for storage in DHS. He optimised the operational costs using different scenarios for pre-heating and heating cut-off periods, showing the added value of intraday heat storage in the building thermal mass.
Vandermeulen et al. [11] recently discussed the challenges involved in developing novel control strategies to unlock the flexibility in thermal energy storage present in district heating and cooling networks and the thermal inertia of buildings. They conclude that prediction uncertainties, the complex network dynamics and network size need further research to ensure an efficient and effective control of thermal networks.
Many studies address some kind of operational optimization [9,[12][13][14], but the integrated optimization of temperaturelimited sources and demand side management has not yet been addressed, despite its practical relevance for integration of renewable low-temperature sources into 4th Generation District Heating (4GDH) systems. Many authors derive an optimization model that is in the class of Mixed-Integer Non-linear programmes (MINLP) [9,11,13]. Model-predictive control with real-time optimization needs to meet a number of requirements in order to get adopted by operators, as experienced by water management operators [15]. An acceptable and feasible solution must be guaranteed within a certain bounded time, small perturbations in problem formulation must result in small differences between solutions and each solution must be a good solution in terms of objective function value. These requirements lead to the requirement of a convex problem definition, for which the detailed building and network physics must be simplified. We will show in this paper that a reformulation of the problem definition leads to a more convenient optimization problem in the class of continuous convex problems, while still addressing optimal integration of temperature-limited sources, temperature-limited capacity of heating elements, thermal energy storage in the building inertia and acceptable indoor temperatures.
In this paper we will describe a simplified heat network model and show how this model is used in a MPC-based operational optimization of a DHS with several heat sources and multiple buildings. In the following section, the model components will be described. Afterwards, it will be shown how the model is used for optimization. This model is applied to different test cases with increasing complexity to show the validity of the simplified model. The paper ends with a discussion of the results obtained and conclusions.

Network components and their requirements
In order to allow a controller to compute an optimal heating schedule, a model is required to describe the relation between heating inputs and building temperatures. For deterministic, realtime optimization to be possible, this model has to be convex. Convex optimization problems have several properties that are strongly recommended for real-time optimizations. First, the solution is unique, such that the solver time is well predictable. Secondly, subsequent solutions with slightly different boundary conditions with a receding time horizon will results in similar solutions, which is required to gain the trust of network operators in such advanced operational decision support. The convexity requirement holds for the objective function as well as the inequality constraints depending on the decision variables. Equality constraints should be linear to obtain a convex problem formulation. In practice, this boils down to developing the simplest possible model that is still sufficiently accurate.
In the following, simplified heat network component models will be derived, including sensitivity assessments to substantiate underlying assumptions. The following components are required to model a heat network: 1. Building 2. Pipe 3. Heat source 4. Heat exchanger It is not necessary nor desired to include valves and pumps in the optimization for a number of reasons. First, we are primarily interested in the transport of heat. Secondly, pump and valve settings can be derived in the post-processing of the optimization. The third reason not to include pumps and valves in the optimization, is the fact that these components introduce complex nonlinearities that complicate the optimization unnecessarily. Obviously, a priori estimates of operational limits in terms of heat and mass flows are taken into account.
The convexity requirement applies to the relations between the decision variables and Building model as well. The external processes, like the ambient temperature, solar influx and building temperature bounds, can be arbitrary temporal functions.
The components are combined to construct an operational optimization model of a district heating system. We select the heat flow and source supply temperatures as decision variables. Furthermore, the primary return temperature is assumed constant. Using this assumption and these decision variables, the mas flow rates from all sources and to all buildings can be derived in a unique way. Therefore, the mass flow rates can be determined as post-processing and do not need to be included in the optimization problem. The supply temperatures are used in the constraints that limit the maximum supply from a heat source and limit the maximum heat supply to the buildings. By using the heat flow, instead of using the mass flows and temperature separately, a linear equation is obtained. If the temperature and mass flow would be used as optimization variables, then the merging of two flows would result in a non-linear equality constraint, which violates the convexity requirement.
In the following sections the equations for the different components are discussed. The symbols in bold fonts denote decision variables in the optimization model. The other symbols denote dependent or external variables in the optimization model.

Building
The main purpose of the building model is to calculate the air temperature inside the building and the short-term storage in the thermal mass of the building. Therefore, an energy balance for the air inside the building is made. This is given by:

Heat transmission
Heat can be lost through the walls, rooftop and foundation of the building and by ventilation (mechanical as well as free flow). The transmission refers to heat loss through the roof, walls and floors, and can be modelled with: In which:A = Area through which the heat is transported [m 2 ], U = Thermal transmittance [W/Km 2 ], T out = Outside air temperature [K] This model equation ignores the thermal storage in the building envelope in order to reduce the number of variables per building in the optimization. The thermal storage of the building envelope is lumped with the building floors as detailed hereafter.

Solar influx
The solar influx is assumed to heat the internal walls and floors of the building. The air is indirectly heated via the floor and walls. The solar influx through the windows onto the building floors and internal walls is calculated based on the method as described in the ASHRAE handbook [16]. The total solar heat input is given by: For every outside façade this solar flux is calculated to get to the total solar influx into the building. The sun will also heat the outside walls and roof. This will not contribute to the heating of the building, since the time the heat requires for penetration through a wall is long. This can be calculated with help of the Fourier number, which is given by: In which:a = Thermal diffusion coefficient of building wall[m 2 /s], t = Time required to reach a certain depth [s] When the Fourier number is below 0.1, the temperature on the opposite side of the material which is heated, has not increased. Therefore, it can be calculated what time it would take for the tem-perature on the inside of the building to increase when the sun is shining on the building wall. For concrete, the thermal diffusion coefficient ranges from 3.7e to 8 to 8.4e-8 m 2 /s. For a wall of 20 cm thick this gives a time of 13 h. In this time most of the solar heat in the wall has been transferred back to the outside air. It is therefore assumed that solar radiation is only entering the buildings through windows.

Ventilation and infiltration
The final source of heat loss in this model is the heat loss due to ventilation and infiltration of air. The ventilation heat loss can be calculated from: In which:g = Efficiency of heat recovery [À], Q flow = Air flow rate of ventilation [m 3 /s] The air flow rate is determined by how often the entire air volume inside the building is refreshed. This is an input for the model. The value will be based upon the occupancies level of the building. Normally the factor will be about 2 to 3 times per hour. The heat loss due to infiltration through cracks, windows and doors can also be calculated based upon a flow rate of air from the outside to the inside. However, the total sum of this flow will be much less than the flow due to mechanical ventilation, since it will not replace the entire building volume of air within an hour. Therefore, this term will be much less than the heat loss due to the mechanical ventilation and can thus be neglected.
Heat storage in floors and internal walls The inside of the building consists of air, floors and internal walls, representing the thermal inertia. The floors and walls will exchange heat with the air and they are also heated by the solar radiation on the building. The total volume of walls and floors, assumed to be constructed from the same materials, is modelled as one thermal element with a certain thickness and surface area. The surface area A is set to the gross floor and wall area inside the building for convenience. The total thickness of this layer is based on internal wall and floor volume divided by gross floor area. This thermal inertia is then discretized in smaller sublayers, see Fig. 1. As can be seen in this figure, this is a symmetric situation and only half of the floor sublayers need to be modelled. The heat balance for the top layer is: The heat transfer coefficient can be determined with standard Nusselt relations over the relevant temperature range. In order to avoid introduction of non-linearities an average heat transfer coefficient must be applied (e.g. 20 W/m 2 K). The area is the gross floor area of the building. For the other layers the heat balance is given by: And for the middle layer the balance is: The total heat transfer to the air inside the building is given by: The factor of 2 comes from the fact that the air is on both sides of the floor layer.
To determine the minimum required number of layers for the optimization several simulations with different numbers of layers have been performed. For these simulations the air and floor temperature start at 20°C, the ambient temperature outside of the building is 0°C. Heat is lost through the building envelope and no other heat sources or losses are taken into account. Fig. 2 shows the result for the temperatures of the air inside the building and the temperature of the floor sublayers. The air temperature drops faster in the first hour or so, since a temperature difference must be created first, before heat transfer starts from the walls to the air. After approximately 15 h, the temperatures in all layers drop at the same rate. Fig. 3 shows the air temperature at different time intervals as a function of the number of layers used. For every time interval the influence of additional sub layers decreases. Above ten layers there is almost hardly any influence. Below 10 sub layers (5 calculated values due to symmetry.) the heat transfer from the thermal inertia is overestimated, leading to higher air temperatures at fewer sublayers. The time scale which is of interest for this problem is in the order of 12 h. From this perspective it has been chosen to model 10 sublayers in a building using 5 calculation values due to symmetry. If fewer layers would be used, then the thermal mass would respond too fast to changes in the air temperature.

Pipe
The pipe is used to transport the heat from one location to another location. In the current implementation for operational optimization the heat loss and the thermal delays (travel time) are neglected. Both simplifications are required to ensure a convex model, since the heat loss would depend upon the supply temperature and delays would be affected by the optimized solution. The validity of the second simplification is disputable, but defended by the following arguments: 1. The thermal delay becomes relevant only when supply temperatures may vary from hour to hour, which is not recommended from a thermal stress point of view. 2. If supply temperature setpoints vary gradually on a daily timescale, then deviations in heat supply due to this assumption can be absorbed by flow variations. During part-load conditions, a slightly larger flow rate is not critical. The only critical condition occurs when the supply temperature must increase quickly to a full load condition in the DHS.
The annual heat loss from a network is a substantial fraction of the total heat production and depends on the time-averaged supply and return temperatures in the network. The optimization promotes lower supply temperatures. Therefore, the network heat loss will reduce in an optimized solution compared to a reference solution without temperature limitations. The network heat loss reduction is an additional benefit which is not quantified in the optimization, because these losses are neglected.

Temperature-limited heat sources
The heat sources supply heat to the system. The governing equation is: The maximum amount of heat, a source can supply, is limited by the capacity of the source itself. It can also be limited by the temperature. For example, a geothermal source can only supply its maximum capacity at a fixed temperature difference between the in-and outflow. Furthermore the outflow temperature is fixed (T s,max ). This adds a constraint to a temperature-limited source in the form of: in which Q source = Heat supplied by the source [W], Q tots = Total heat delivered by all sources [W], T prim = Maximum of primary supply temperature of the network and maximum source supply temperature to ensure that the ratio is always less than or equal to 1.
[K], T ret = Primary return temperature from buildings, assumed constant and equal to the design value of the DH network [K], DT s = the design temperature difference at the heat exchanger of the source, based on maximum source supply temperature T s,max and the design value of the primary return temperature T ret :- The key assumption for the optimization is the constant return temperature in the DHS. This assumption is realistic, if the hydraulic controls of the secondary systems are properly designed and implemented [17]. This is a non-linear non-convex constraint, since both the primary temperature and the heat supplied by the sources are decision variables. To linearize the equation, it is rewritten to: In this equation only the first term is non-linear. In order to include this, the homotopy approach is used. The constraint is linearized around a nominal point (Q snom ,T pnom ). This gives: The non-linear equation is then multiplied with a factor X and the linearised equation is multiplied with the factor 1 À X ð Þ . The optimization is then started with this X = 0, basically solving the linear constraint. When a solution is found, this factor is increased and another optimization is performed. The starting point for the next optimization is the solution of the previous optimization. This is continued until X = 1 and the solution is found for the complete non-linear problem. The main reason to maintain this non-convex constraint is illustrated in Fig. 4 hereafter. When the DH system needs a supply temperature, exceeding the maximum temperature of the renewable source, then the renewable source will deliver the largest possible fraction of the total heat supply, which corresponds with a point on the nonlinear (blue) line in Fig. 4. Fig. 4 clearly shows that the linearised equation underestimates the maximum supply from the renewable source. We could use a priori knowledge to select a better linear approximation; e.g. a linearization at 80°C, but that would still result in significant underestimation of the renewable source capacity at 70°C or 100°C. Since the optimal solution will always hit the non-linear boundary, when the supply temperature has to exceed the maximum temperature of the renewable source, we  decided to iterate towards this optimum solution using the homotopy approach.
Next to this constraint there is also a constraint on the temporal supply temperature gradient. This is given by: In which dT max = Maximum temperature change per unit of time for the source; typically equivalent to 5 K/h.[K/s]

Heat exchanger with temperature-limited supply
The heat exchanger (HEX) represents the heat emission system of a building by supplying heat from the district heating system to a building. The governing equation is: The amount of heat supplied is limited by the HEX capacity or the radiator capacity. The design temperature inside the buildings can be assumed constant around 293 K. The maximum supply from the heat emission system (Q max ) is determined based on the capacity at a design temperature T d by: This equation simplifies the log-averaged temperature dependency of heat emission systems. It is based on the temperature inside the radiator which is on the secondary side of the heat exchanger. Therefore, an additional step is required to couple temperature at the secondary side to the temperature on the primary side. The heat exchanger is designed to work with a certain temperature difference between the primary and secondary side. If we assume this difference to be constant over the entire temperature range the primary temperature is given by: In which T prim = Temperature at the primary (inlet) side of the heat exchanger [W], DT hex = Design temperature difference between the primary and secondary side of the building HEX. [K] Finally, combining the above equations leads to a temperaturelimited constraint on the supplied heat to the building: Equation (18) is used to limit the supply to the buildings at a reduced primary supply temperature, taking into account the design capacity of the radiators. In fact, this constraint expresses that the mass flow rates on both sides of the heat exchangers cannot exceed the design value. Therefore, the mass flow rates towards the buildings will not exceed the design value.

Optimization model
The above components can be used to create an operational optimization model of a DHS, defined by a set of linear and linearized algebraic-differential equations in the decision variables. The decision variables are the heat supplies from the sources, the heating and cooling demands by the buildings and the supply temperatures per source. This is combined with a linear objective function. It is possible to use different objective functions. The objective function for these cases minimizes the total costs, defined as follows: The cost per unit of power can also vary in time (e.g. cheap by night expensive during daytime). The pumping energy is in the order of 1% of the heat supply [10] and therefore ignored in the objective function. Next to this a set of constraints is added to the model. Limiting the capacity of the sources and the heat supplied to the buildings has already been discussed. Next to this the temperature inside the buildings should stay within certain bounds. During office hours these bounds are rather strict while during night time the bounds are extended to promote peak shaving by using the thermal inertia of the buildings, see Fig. 5. We have used the operational optimization toolbox RTC-Tools [18,19] to model and solve the optimization problems. RTC-Tools can use different solvers, depending on the problem formulation. In this case the continuous non-linear solver IPOPT [20] has been used.

Test case
To study the results of the optimization model as implemented within RTC-Tools a simple test case is used. This test case consists of three office buildings and two sources, which have been derived from a real-life example. Fig. 6 shows an overview of the system. Table 1 and Table 2 show an overview of the main properties of the office buildings. Table 3 shows the properties of the sources.
A maximum supply temperature of 70°C is a typical geothermal source temperature in the Netherlands, extracted from 2 km deep reservoirs. If the renewable source would be a data centre or surface water, then then a central heat pump preferably lifts the temperature to 70°C at most. The constant return temperature in these test cases is 40°C.

Results
Section 4 presents optimization results of the developed modelpredictive control solver. First the time horizon boundary effects are discussed to assess start-up and end-of-horizon effects in the results. Section 4.2 details optimization results for scenarios with increasing complexity: 1) a reference case with constant outside temperature, 2) realistic environmental conditions, 3) impact of ventilation, 4) source capacity constraints, 5) impact of difference source prices. 6) impact of time-varying heating costs and 7) impact of a temperature-limited source.

Time horizon boundary effects
Before optimizations are performed and analyzed first the effect at the boundaries of the time horizon used for the optimization need to be investigated. An optimization has been performed to investigate the effect of the boundaries for the example system. In this optimization no additional heat flows (solar, ventilation etc.) are taken into account. Furthermore, the outside temperature is kept constant at 0°C. A period of 31 days has been optimized. For each day the same constraint on the temperature inside the buildings is used (Fig. 5). There are no temperature constraints on the source, which makes the model fully linear.
By using these constraints, it is expected that the time series over the successive days will be similar. Only at the beginning and end of the optimization time horizon it is expected that there is a difference. Fig. 7 shows the time series of the heat supplied by both sources. The mean of the 31 days is plotted (black line) together with this mean plus and minus 2 times the standard deviation. Furthermore, the heat supplied at day 1 (red), 2 (blue) and 31 (dashed red line) is shown. The first day clearly shows a different behavior than the other days due to initialization effects. Hence, day 1 should be neglected in the analysis of further results, since the difference is most of the times more than two times the standard deviation. Fig. 8 shows the mean air temperature of building 1 for the 31 days of the optimization together with the air temperature at day 30 and 31. As can be seen the temperature at day 31 shows a large deviation of the mean, especially at the end of the day.  Based on this it is concluded that this day should also be neglected in the results. Therefore, every optimization should be performed with a start-up day and an end day, which are not considered in the results.
The results from this section must be used for long-term optimizations as well, such as scenarios for a climate year. For such studies the optimization is split in periods of 1 week or 1 month, depending on the problem size. The receding time horizon must use at least two overlapping days as concluded from Fig. 7 and Fig. 8: one start-up day and one day to obtain results for the last day of the previous horizon.

Optimization results
In the remainder of the results section we will analyse and reflect on the optimization results. The main goal is to show that the optimization model produces results which are understandable and logical. This can serve as a first proof that the optimization model can be used to control a real DHS.
Unconstrained solution (reference case) Before discussing the effect of several parameters on the solution the so-called basic solution is discussed. This is an optimization for one week, without any constraint on the sources and heat supplied to the buildings. The model includes a periodic temperature constraint on the air temperature inside the buildings (dashed lines in Fig. 9) with narrow temperature limits during office hours and wider limits outside office hours. Furthermore, all external and internal heat fluxes to the building have been set to zero, except for the heat loss to the surroundings. The outside temperature is set to 0°C for the entire week to investigate the impact of different parameters. Fig. 9 shows the results for the air temperature inside the buildings, while Fig. 10 shows the total heat delivered by the sources. As expected, the temperature is kept as much as possible at the lower bound. If the temperature would be increased above the lower bound, this would yield additional heat loss and increase the costs. During the night time almost no heat is supplied to the buildings and they cool down. As soon as a building reaches the lower bound some heat is supplied to keep it at this bound. Also visible in Fig. 10 is the high peak heat demand in the morning to compensate for the heat losses during the night. Conclusion: the optimization model yields acceptable and understandable results for this basic case showing typical dynamics in the inside building temperature and heat supply.

Effect of environmental parameters
In the previous optimization the effect of the solar influx and changes in the ambient temperature in time were not considered. Both can be taken into account by using climate design data. Dutch climate data for the first week of January and the first week of May have been used. Fig. 11 shows the ambient temperature for these weeks. Fig. 12 shows the solar influx to all the buildings, taking the cloud coverage into account. Fig. 13 and Fig. 14 show the results for the first week of January, showing marginal differences between these results and the one of the basic case. The main difference is that the amount of heat required now changes for the different days, due to the variation in ambient temperature. During the weekend the solar influx into the buildings is visible in the temperature results. This heats the top sub-layer of the floor, which in turn heats the air inside the building. This is visible around hour 130 in Fig. 13 at that time the heat supplied by the sources is zero, while there is an increase of building temperature, see Fig. 14. Fig. 15 and Fig. 16 show the results for the first week of May. The required heat reduces a lot, since the ambient temperature is much higher and there is a high solar influx into the buildings.
The temperature inside the buildings is almost always against the upper bound, while in January it was at the lower bound. To keep the temperature within the bounds, cooling must be supplied to the buildings, which can be seen in Fig. 17. Despite the large cooling demands, there is still some heat demand at the end of the night to maintain the inside temperature within the bounds.
Conclusion, the optimization can cope with changes in environmental conditions and still produces understandable results and correct dynamics in the thermal energy supply and temperatures.
The impact of a schedule of internal heat gains for electricity consumption and occupancy can be included in the optimization in a similar way, since the internal temperature is linearly related to these internal heat gains following equation (1).
Effect of ventilation In the previous optimization ventilation was not included. This effect is taken into account with equation (5). An efficiency of 65% is used, which is an average efficiency for a heat recovery unit. The air flow is based on air change rate, which denotes how many times the total air volume in the building is refreshed within one hour. During office hours this factor is about 2-3, for other hours it is between 0-1.5. Fig. 18 and Fig. 19 show the results. As can  be seen the total amount of heat required is more than in the base case with the environmental parameters turned off (Fig. 10), the maximum required power rises from 3.6 to 4.2 MW.
Effect of limited source capacity In the previous section the focus was on the building. In this section we will focus on sources with limited capacity. As a first step the maximum heat from the sources is limited to 2 MW. Fig. 20 shows the temperature inside the buildings and Fig. 21 shows the amount of heat supplied by the sources. As can be seen it reduces the morning peak to 2.0 MW. During the week days building 1 is preheated first, this is followed by building 2; building 3 is kept at the lowest possible temperature. This order is found by the solver, because this order minimizes the overall heat losses due to insulation and size of the buildings. This can also be seen in Fig. 9 of the reference scenario, where the decrease in temperature during night time is the smallest for building 1, then building 2 and finally building 3.
Next to this, the solver proposes to heat building 1 to the maximum possible value during the weekend. This pre-heating will store energy in the concrete floors of the buildings. These will start to provide heat to the building when it cools down again. Conclusion: adding bounds to the source capacity successfully forces the optimization solver to find heat supply profiles with load shifting and peak shaving.
Effect of heating costs To investigate the effect of the cost function two sources are included: one source with a low price, but with a maximum capacity of 1.0 MW (representing a renewable source) and a second source with a price which is 10 times higher. This source has enough capacity to deliver the peak demand (typical for gas boiler). Fig. 22 and Fig. 23 show the result. As can be seen source 1 (the cheapest source) is the preferred source and is almost constantly supplying its maximum power, while source 2 only supplies power at mornings to heat the buildings to the lower bounds of the allowable temperature range. Furthermore, building 1 is preheated much earlier than in the previous cases due to the lower costs of the low capacity source. Hence, the optimization finds the proper balance between the load shifting for pre-heating and the overall  heating costs. The maximum heating power in this scenario is just over 2 MW and half the required heating power for the scenario with a heat source with ample capacity.
Effect of time varying heating costs Fig. 24 shows the variation of the costs as function of the time. Source 1 is cheap during the night time and more expensive during daytime, source 2 has a constant price, lower than the daytime price of source 1. Both sources have unlimited capacity. Fig. 25 and Fig. 26 show the results. Source 1 preheats building 1 and 2 during the night. Furthermore, during daytime source 1 is switched off and source 2 is supplying, as expected. Fig. 26 illustrates that variable prices lead to excessive supply variations when the price is low. These dynamic variations are not recommended for many renewable sources. It is concluded that the optimization model can cope with varying costs function.
Effect of temperature-limitations of sources and heat emission systems In this final section the effect of temperature-limited sources (equation (12)) and the rate of change of the primary supply temperature (equation (14)) are included. Also the temperature limitation of heat emission systems at the building level is included as given by equation (18). Source 1 is bound by capacity and temperature as given in Table 3. These capacity and temperature limitations are characteristic for renewable heat sources like a geothermal source, surface water or low-temperature datacenter heat. The unit cost of renewable source 1 is assumed 10 times lower than the unit cost of source 2. This cost could represent the CO 2 emission, the actual operational costs or a combination of these parameters. Fig. 27 shows the results for the temperature within the 3 buildings. Fig. 28 shows the heat supplied by the 2 sources and the maximum amount of heat that source 1 could supply due to the temperature limitation. Fig. 29 shows the heat supplied to the buildings. Furthermore, this figure shows the maximum possible heat supply at the optimized supply temperature, due to the constraint of the heat emission system at building level. Fig. 29 shows that the heat supplied to the building 3 is almost always at the maximum possible. Finally, Fig. 30 shows the primary temperature of the district heating system. This minimum supply temperature is just sufficient, such that all buildings together use source 1 as much as possible; see Fig. 28. It is concluded that all degrees of freedom for optimization are fully used.  These results are compared with a heat supply based on a conventional heating curve. This curve is given by the following relation (temperatures in°C): With the above equation the primary temperature is no longer a decision variable and equation (11) is no longer non-linear. Using this equation yields the primary temperature as a function of the time as given in Fig. 30. Fig. 31 shows how much heat can be supplied from the renewable source 1 for the 2 operational control strategies: optimized control versus heat curve control. It is emphasized that renewable source 1 is limited in supply temperature and power (2 MW at 70°C, Table 3). The heat supplied for the optimized case is much more at a stable level, while for the case with the heating curve, the day-night pattern in the heat demand is more visible and the supply from the renewable source 1 is much smaller. When the primary temperature is not fixed, this can be reduced, which enables source 1 to supply a large baseload heat even when source 2 is not available. Since the cost for source 1 are much lower it results in an even spread of the heat supplied by source 1 as can be seen in Fig. 31.

Discussion
The optimization model gives understandable and promising results for all scenarios with increasing complexity. However, there are still some points of discussion in the applicability of the model. The first one is the correctness of the model, especially for the building. The buildings in the test cases have been modelled as single zone models. These models can be extended to multiple zones without violating the conclusions of this paper. Furthermore, the structure of the building model is similar to the RC-modelling approach adopted and validated by Bacher and Madsen [21]. Bacher's model includes three heat capacity terms for the internal air, the envelope and the radiators. Our model includes six heat capacity terms for the internal air and five thermal storage layers, representing the floor and walls. It was shown that a 5-layer ther-  mal storage model is physically realistic. Therefore, we expect this building model to be accurate enough for practical applications.
Nevertheless, it is recommended to validate the indoor temperature of the building model with more detailed building models and measurements. Once this simplified building model is sufficiently validated, it is straightforward to extend the optimization to a large number of buildings by aggregating similar buildings with similar orientation in the same neighborhood to a single aggregated building.
Another point of discussion is the air temperature constraint inside the building as given in Fig. 5. In some cases, the air temperature is heated towards the upper temperature bound outside working hours. This is done to store heat in the floor of the buildings. However, such a high temperature might not be acceptable for the plants and other green elements inside the building. Such a wide temperature bandwidth is more suitable for office buildings than for residential buildings. This constraint can be adapted accordingly. Limiting the maximum air temperature will reduce the bandwidth for optimization and thus will result in reduced flexibility and higher costs. These air temperature constraints can be set by the user and thus are not a limiting factor for the applicability of the optimization model. The model can even be used to show the effect of these constraints on the operational costs.
The advantage of using this kind of optimization for the control of a district heating system has been shown in the last example.
Here it was shown that a capacity and temperature-limited source (e.g. geothermal) could be utilized much better, when the optimization was used instead of a fixed heating curve. However, to further investigate the advantage a detailed case study needs to be performed of an existing district heating system.
A hidden benefit of this optimization approach is related to the network heat losses. The network heat loss reduction, due to lower supply temperatures, is an additional benefit which is not quantified in the optimization, because the network heat losses are neglected in the optimization.
A final point of attention is the time used for the optimization runs. If this optimization needs to be used to control the system, the run time for the optimization should be much smaller than             the prediction horizon. At this moment the full optimization problem including the non-linear constraint takes 104 s on an Intel core I7, 2.4 GHz. It should be noted that this was only running the optimization for one week and for a simple system. This time is well within the time boundary for practical applications with hourly updates.

Conclusions and recommendations
This paper describes and analyses an operational modelpredictive control optimization model for a district heating system with capacity-and temperature-limited sources. The operational optimization also takes the influence of the supply temperature on the heat emission systems into account.
Whereas many authors have proposed operational optimizations as mixed-integer non-linear problems, we managed to define a suitable model-predictive control optimization that meets all requirements for practical applicability. The reformulation of the problem definition leads to a convenient optimization problem in the class of continuous (nearly) convex problems, while still addressing optimal integration of temperature-limited sources, temperature-limited capacity of heating elements, thermal energy storage in the building inertia and acceptable indoor temperatures. The operational optimization finds the proper balance between benefits of pre-heating using renewable sources with limited capacity and costs of additional heat losses due to pre-heating.
The test case shows that the peak supply reduces significantly from 4 MW (Fig. 14) to 1.8 MW (Fig. 28). Furthermore, the supply temperature is minimised such that the cheap temperaturelimited source can deliver most of the heat demand in continuous operation. An additional benefit of the reduced supply temperature is the reduction in network heat loss. Further benefits for the district heating network include reduced variations in network flows.
The analysis of the building model performance, using concrete floor parameters, shows that at least 5 floor layers must be included in the building model to properly account for the thermal inertia of buildings with concrete floors.
In addition to operational applications, the developed modelpredictive control solver can be used to quantify the heat flexibility of districts over a wide range of annual scenarios, including variable pricing and limitations in temperature and capacity of renew-  able heat sources. The annual scenarios must use a receding time horizon with two overlapping days as detailed in section 4.1.
It is concluded from the test case that the optimization model gives optimal results for a wide range of realistic scenarios. For future research it is recommended to study an existing district heating system. Furthermore, a comparison of the results of the single zone building model with results from a more complex model and measurements should be made.

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.