Management of a district heating network using model predictive control with and without thermal storage

District heating and cooling networks are a key infrastructure to decarbonise the heating and cooling sector. Besides the design of new networks according to the principles of the 4th and 5th generation, operational aspects may significantly contribute to improve the efficiency of existing networks from both economic and environmental standpoints. This article is the second step of a work that aims to exploit the flexibility of existing networks and improve their economic and environmental performance, using the district heating network of Verona as a case study. In particular, the first part of the research demonstrated through numerical simulations that the thermal inertia of the water contained in the pipes can be used to shift the heat production of the generators over time by acting on the flow rate circulating in the network. This article shifts the focus from the heat distribution side to the heat supply. A model predictive control strategy was formulated as a MILP optimization problem to schedule the heat supply of the cogeneration plants, heat pump and gas boilers as a function of heat load, waste heat production and electricity price forecasts. Computer simulations of considered district heating network were carried out executing the optimization with a rolling-horizon scheme over two typical weeks. Results show that the proposed look-ahead control achieves a reduction in the operational costs of about 12.5% and 5.8%, respectively in a middle season and a winter representative week. Increasing the flexibility of the system with a centralized heat storage tank connected to the CHP and HP units, these percentage rise to respectively 20% and 6.3%. In the warmest periods, when the total installed power of the CHP and HP plants is sufficient to supply the entire heat demand during the peak, and the modulation of these plants has a higher impact, the cost reduction related to the additional thermal energy storage is more relevant.


Introduction
District heating networks (DHN) play a key role in the transition towards sustainable cities, thanks to their ability to efficiently provide space heating and domestic hot water to buildings through the use of renewable sources, waste heat sources and cogeneration plants. Since almost half of the final energy use in Europe is related to heating services, the construction and the expansion of efficient DH systems are some of the key points of European energy policy (Directive (EU) 2018). Besides the design of new DH networks and the refurbishment of existing ones according to the principles of fourth generation (Lund et al. 2014) and fifth generation (Buffa et al. 2019) district heating and cooling networks, improving the operational aspects of traditional district heating networks is important to increase their efficiency from both economic and environmental standpoints.
The optimal operation of heat supply units in smart grids is not a new topic in the scientific literature. In this context, several research papers show that Model Predictive Control may be a viable option to efficiently match energy supply and demand. For instance, Fink et al. (2015) proposed two control methods based on mixed-integer linear programming (named global MILP and time-scale MILP) to flatten the power consumption profile of a group of heat pumps supplied by a biogas-fired CHP that also supplies heat to a district heating network. Both methods showed promising performance in flattening the power demand of the heat pumps compared to a standard PI on-off control. Vivian et al. (2020a) showed that a centralized management of heat pumps can significantly reduce the electrical peak load in a residential district using distributed thermal storage tanks for space heating and domestic hot water production.
As far as district heating networks are concerned, different objectives have been pursued such as the peak shaving or the reduction of costs and primary energy consumption. To this end, different flexibility sources can be exploited, namely the thermal inertia of buildings, small thermal storage tanks in the district heating substations, the thermal inertia of the heat carrier fluid in the network pipelines or large thermal storage tanks in heat supply stations. Guelpa (2021) demonstrated that the energy absorbed by the thermal capacitances of these components can be significant, especially during shutdown or setback of the generators, and that a proper management of thermal capacitances for energy saving can be considered. According to Vandermeulen et al. (2018), the contribution of the network pipes is limited compared to that of the buildings' thermal inertia. An example of this kind has been shown by Aoun et al. (2019), who used reduced order models to predict the space heating needs of buildings connected to a DHN and shift their consumption towards off-peak hours maintaining thermal comfort in the indoor environment. Despite this big potential, Demand Side Management programs are still rare in the DH world because DH operators consider the user substations as the point where their operating space ends.
An alternative optimized management at the user substations was proposed by Capone et al. (2019) that was able to reduce the thermal power peak of the network from 20 to 42% depending on the maximum allowed anticipation of the mass flow request during the early morning.
As far as the thermal inertia of the heat carrier fluid is concerned, in a previous paper the authors obtained significant results in terms of thermal load shifting and peak shaving by modulating the circulating mass flow rate via by-pass valves (Vivian et al. 2020b). It was demonstrated that the control parameters (e.g. start of the network pre-charge) vary according to the period considered because in the coldest periods of the year the flow velocity of the heat carrier fluid is higher and therefore the response times of the network as thermal storage are faster than in the midseason. With regard to the optimization of the heat generators, Bavière and Vallée (2018) focused on the design of a model predictive controller that optimizes the distribution of heat by an appropriate scheduling of supply temperatures and differential pressure at the production level. In the early work by Benonysson et al. (1995), the supply temperature was optimized to minimize the operational costs. Guelpa et al. (2017) proposed a model predictive control to minimize the primary energy consumption at the thermal plants through the reduction of the thermal peak request and allow the full exploitation of the cogeneration plants. The importance of the integrated management of cogeneration plants and district heating networks in terms of economic competitiveness and energy saving has been underlined also by Pini Prato et al. (2012). Schweiger et al. (2017) proposed a new framework also based on Modelica, allowing full dynamic simulation and optimization of low temperature district heating system. Verrilli et al. (2017) proposed an optimal control strategy for district heating power plants with TES and they demonstrated its robustness against heat demand forecast uncertainty. Also Vivian et al. (2017) developed an optimal control to manage a low temperature DHN with high share of waste heat. This is able to guarantee a correct hydraulic and thermal balance of the system, and to minimize the operational costs by using a central storage tank and the network pipelines as buffers. Lésko et al. (2018) proposed linear models to include three different sources of energy flexibility (hot water tanks, thermal inertia of both the DH network and the connected buildings) that can be integrated in operational optimization programs of DH systems. Vesterlund et al. (2017) developed a hybrid evolutionary-MILP optimization algorithm and coupled it to thermal network model simulations in order to achieve an optimal control of heat generators in a complex DH network with multiple heat sources. Including complex thermal networks in the constraints of the optimization problem introduces nonlinear phenomena such as mass and heat propagation that are governed by the advection equation. Recently, nonlinear optimization approaches have been proposed by Krug et al. (2020) and by Rein et al. (2020), whereas other researchers preferred to approximate nonlinear phenomena through piecewise linearization techniques (Bordin 2015;Giraud et al. 2017). However, it is unclear whether these simplified models approximate the true state with sufficient accuracy. An alternative to piecewise linearization has been proposed by Mehrmann et al. (2018) for gas networks, where the error is controlled to adapt the space discretization of the physical model in the optimization constraints. The thermal inertia of both district heating network and buildings has been exploited by Li et al. (2020), using an optimization scheduling model. They concluded that the heat storage capacity of the district heating system and aggregated buildings effectively improves the flexibility of the electricity system, allows CHP units to be operated more flexibly and promotes wind power integration, without affecting the comfort requirement of the demand side. Similar results were obtained by Gu et al. (2017). Finally, Saletti et al. (2021) developed a novel approach to consider the thermal capacitance of the connected buildings in a large-scale heating networks. It was include in a Model Predictive Controller that optimizes the management of the production plants and the heat distribution.
In this paper, a mixed integer linear programming problem has been formulated and used to exploit the flexibility offered by the heat carrier fluid in the pipes and by an additional thermal energy storage tank. The optimization problem has been coupled to a detailed district heating network model for the optimal control of the heat supply units in a receding horizon manner assuming ideal forecasts of heat demand, waste heat availability and electricity price. In the simulations the mass flow rate circulating in the network has been iteratively adapted in order to follow the state predicted by the optimal controller.
The objective of this work is to assess to what extent the aforementioned flexibility sources contribute to the reduction of operational costs for the DH operator.
The novelty of the study, with respect to the current state-of-the-art research, is determined by the comparison between the two considered sources of energy flexibility in a real case study and by the method used to exploit the thermal inertia of the heat carrier fluid in the network pipelines. Unlike several studies, where the supply temperature set-point is adjusted, here the control variable is the circulating mass flow rate. The mass flow rate variations have the purpose of anticipating or delaying the heat supply with respect to the situation where heat generators strictly follow the demand of the district heating customers, exploiting the thermal capacitance of the water in the network and minimizing the operational costs.
The paper is organized as follows. In Sect. 2 we present the optimization problem and the finite difference method model of the district heating network, that we used for operational optimization in large-scale DHN. In Sect. 3 we give an overview of the district heating network under consideration. In Sect. 4 we describe the coupling of the optimal control problem with the heat distribution model in the model predictive control (MPC), the method and the main assumptions considered to perform the simulations. In Sect. 5 we present the results of the optimization problem without coupling to the district heating model and the simulation results of the model predictive control. Finally in Sect. 6 we report on the conclusions.

Models
This paper is focused on the optimized generation of heat in district heating networks. The first Subsection of the present Section outlines in detail the optimization problem used to schedule the heat generators. The second Subsection describes the thermal network model. The coupling of the two models is explained in the Methods Section.

The optimization problem
The optimization problem relies on a Mixed-Integer Linear Programming (MILP) formulation, meaning that all the relations between physical variables (e.g. energy balances) are formulated as linear equality and inequality constraints. This Subsection resumes the optimization variables, the constraints and the cost function to minimize. The values of the variables are discretized and considered constant through the time step. The problem has 27 "physical" independent variables that are summarized in Table 1. If the time horizon of the optimization is set equal to 24 h, there is one schedule of 24 values for each "physical" variable. Therefore, the problem has 27 × 24 = 648 variables (of which 552 are integers).
The network energy balance is considered as a linear constraint as shown in Eq. (1). The index j in the equations refers to the j-th hour. Thus, the following holds true for j = 1…24: where V is the total volume of water enclosed in the network pipelines; Q gb,min is the minimum thermal energy when the gas boiler is in operation; Q wh,j is the available waste heat recovered from an industrial process at j-th hour; Q dem,j is the total Management of a district heating network using model predictive… heat demand at the user substations at j-th hour; UA is the heat loss coefficient of the pipes; T supp,j and T g,j are the network supply temperature setpoint and the ground temperature at j-th hour. By shifting all the unknown variables to the left side Eq. (1) becomes: Including the thermal capacitance of grid (1) means that the network is modelled as an equivalent thermal storage system, which allows the optimization to store heat in the thermal grid when convenient. According to this strategy, the state of charge of the network is indicated by only one value, which can be interpreted as the average temperature of the heat carrier fluid in the return pipes ( T r ). In Eq. (3) there is the energy balance of the storage tank, which is considered connected only to cogeneration units and heat pumps, at the generation side.
where Q chp,nom is the heat production of each cogeneration unit in nominal condition; Q hp the heat production of each heat pump connected to the corresponding CHP unit; T amb the temperature in the supply station where there is the storage tank.
In Eq. (4), the number of operating cogeneration plants is equal to the sum of the on/off status of each CHP unit: In Eq. (5), the number of the operating heat pumps is equal to the sum of the on/ off state of each heat pump.
The optimization problem is characterised by 648 variables and 96 equations. Therefore, the optimization problem has 552 degrees of freedom, which correspond to the optimal schedules of 23 decision variables. Inequalities (6) and (7) express that the heat output of the gas boiler is modulated between Q gb,min and Q gb,max and only if the gas boiler is on ( x b,j = 1). Inequality (8) defines the minimum amount of the total thermal output delivered to the network at each time step using the coefficient k dem , that ranges from 0 to 1. In the first case, there could be moments with no heat delivered to the networks whereas higher values of k dem make the heat production follow the heat demand. Inequalities (9)-(11) define the moment when a startup of the CHP occurs ( x su,j = 1), which is needed to apply a penalty for an excessive number of start-ups. They are repeated 5 times, once for each CHP unit. Inequalities (12)-(14), instead, define the moment when a shut-down of any CHP units occurs ( x sd,j = 1). Therefore, they are also repeated for five times. Ineq. (15) establishes that a heat pump can be switched on only if the corresponding cogeneration engine is running. Inequalities (16)-(17) establish that at the end of the simulation horizon respectively the average temperature of the storage tank and the return temperature of the network are greater or equal to the initial ones. These constraints avoid that at the end of the simulation the energy levels of the tank and of the network are lower than the initial ones, providing energy for free. Furthermore, inequality (18) constrains the heat supplied by the storage tank to the network as a function of the difference between the average temperature inside the tank and the return network temperature, with a heat transfer coefficient (k he ).
The objective function of the optimization aims to minimize the overall operational costs of the system, as given in Eq. (19).
where c gas is the cost of the gas used by the gas boilers and the internal combustion engines; Q chp,nom the heat production of each natural gas cogeneration engine operating in nominal condition; chp the efficiency of each cogeneration engine; p el,sell is the selling price of electricity; is the ratio between the power and heat production of the cogeneration engine in nominal condition; c man,chp the maintenance cost of each cogeneration engine; Q hp,nom the heat production of each heat pump operating in nominal condition;a hp is the excise duty paid for the electrical energy used in the heat pump; COP is the coefficient of performance of the heat pump; c man,hp the maintenance cost of each heat pump; gb is the efficiency of the gas boiler; c man,gb the maintenance cost of the gas boiler; c wh is the cost of the waste heat. This cost can be discounted by the Italian energy efficiency certificates. c TEE is the monetary value of each energy efficacy certificate; TEE is the number of efficiency certificates corresponding to 1 MWh of recovered waste heat. In conclusion, the optimization problem is subjected to Eqs. (2)-(5) and to inequalities (6)-(18) and consists of the minimization of the objective function f .

District heating network model (NeMo)
The district heating network model NeMo has been already described in the aforementioned paper (Vivian et al. 2018). The topology of the network is described by using graph theory. The network is represented by a set of nodes and oriented branches and an adjacency matrix determines their mutual connections. Once the geometry is established, the pressure and temperature profiles are calculated. When forced convection occurs, the flow velocity of the heat carrier fluid does not depend on the temperature distribution. Therefore, the hydraulic and thermal sub-models can be uncoupled. This allows the calculation of the mass flow rates and the pressures across the network in a first step; then, given the mass flow rates, the energy balance is performed to find out the temperature distribution. The model assumes a slug flow (i.e., one dimensional model) and neglects both the heat conduction in the axial direction and the thermal capacitance of the surrounding ground. The heat transfer in the radial direction considers the convection between the heat carrier fluid and the inner pipe surface, the thermal insulation of the pipe and the thermal resistance of the surrounding ground. Due to the uncompressible nature of the heat carrier fluid, the hydraulic problem can be described using only two equations: the continuity and the momentum equations. NeMo solves Management of a district heating network using model predictive… these equations using the SIMPLE method proposed by Patankar and Spalding (1972). The heat propagation in the network is then described by the energy balance performed on the volume of heat carrier fluid around the nodes of the network. The control volume of the i-th node corresponds to half of the heat carrier fluid volume of all the branches connected to it. Applying the energy balance to the node shown in Fig. 1 leads to Eq. (20): where G are the mass flow rates, V is the volume of heat carrier fluid enclosed in the control volume, Ω is the perimeter of the pipe section, U is the radial heat transmission coefficient from fluid to the ground and T g is the undisturbed ground temperature. The temperature of the branches is then associated to the temperature of the corresponding upwind nodes, according to the well-known upwind scheme. Therefore, Eq. (20) becomes: Equation (21) can be represented in matrix form as: where M and K are the so-called mass matrix and stiffness matrix, respectively. The temperature at the inlet node is fixed (Dirichlet condition) and the missing mass is attributed to the adjacent nodes. The first-order ordinary differential equation (ODE) (22) can be rewritten to give a linear system that can be solved by Gauss elimination methods as shown in Eq. (23) and (24): where T 0 represents the temperature vector with the temperature values of the preceding time-step (initial network temperature at the beginning of the simulation).
In the current version of the model NeMo, the user can choose the resolution method for the transient heat propagation problem between on the ODE solver (Eq. 22) (Shampine and Reichelt 1997) or the linear system solver (Eq. 24). The latter allows the user to set the time-step of the internal solver Δt.

Description of the case study
The case study under consideration is the district heating network of Verona Centro Città (Italy), shown in Fig. 2. This DHN operates at constant supply temperature and variable flow rate. Furthermore, the network extends for a length of about 25 km with a volume of 652 m 3 , providing heat for space heating and-in some casesdomestic hot water to 247 users. The user buildings, which are supplied by the network, have a total volume of about 3.2 mm 3 and they need about 70 GWh/year of heat, with a peak load of approximately 38 MW. The blue dots in Fig. 2 represent the substations of the users, while the three green dots correspond to the supply stations. The main one is CCC, that supplies heat produced by five CHP units, heat pumps (supplied by the cooling circuit of the CHP) and three auxiliary gas boilers; CRV recovers heat from a foundry and CSD has other three auxiliary and reserve gas boilers. CRV can always supply heat to the network, when available from the industrial process.
The design values of the heat supply units are summarized in Table 2. A rulebased control has been implemented and fine-tuned over the years to prevent discomfort due to low mass flow circulation and to keep the supply temperature around  . The priority in the production is given to CRV and then to the CHP units. When their combined production is not sufficient to reach the supply setpoint, more heat is supplied by the gas boilers. The simulations also consider a centralized heat storage tank (HS) connected only to the cogeneration units and to the heat pumps, at the generation side of CCC station. The considered thermal energy storage has a volume of 500 m 3 and its actual installation is being currently considered by the DH operator to improve the flexibility of the heat supply.

Methods
The optimization problem was executed iteratively within a district heating network simulation of the thermal grid described in the previous Section. The coupling of the optimal control problem with the heat distribution model in the model predictive control (MPC), allows to (1) assess the benefits of the proposed control strategy on a time interval bigger than the time horizon of the model predictive controller and (2) consider the deviations of the district heating system from the predictions of the controller, thus improving the reliability of the results.

Simulation framework
As shown in Fig. 3, the set of optimization variables that the determine the heat generated at each time-step is then used to determine the total mass flow rate circulating in the network. The mass flow rate modulation, as shown in a previous study (Vivian et al. 2018), allows the DH system operator to charge and discharge the network to pursue an economical objective. The objective is calculated at each time-step by the optimal controller to minimize the operational costs depending on 24-h ahead forecasts of heat demand, waste heat availability and electricity prices. The return network temperature, which is an output of the detailed network simulation, is used Fig. 3 Schematic illustration of the coupling between the optimization problem and the detailed district heating network simulation as a feedback signal from the system and it is updated with hourly time-step. So, the decisions taken by the optimal controller are always linked to a measurable physical variable that represents the "current state" of the network. The mass flow rate setpoint is thus the link between the output of the optimization and the input of the district heating network model.
In reality, the DH operator does not control directly the mass flow rate but only the head of the circulation pumps. Nonetheless, the total circulating mass flow rate could be varied by the operator by opening and closing a limited number of by-pass valves installed in strategic points of the network. The opening of the valves relies on the mass flow rate setpoint, that is currently set through Eq. (25): Equation (25) outlines that the mass flow rate setpoint depends both on the current heat demand Q dem and the current heat supply Q gen . When α = 0 the situation is similar to the current flow control, where the mass flow rate only depends on the users heat demand. In the opposite situation (α = 1), the mass flow rate only depends on the decisions taken by the optimal controller. This setting would lead to moments when the heat carrier fluid is not circulating and all the heat is provided by the network temperature variation. This situation, predicted by the linear model inside the controller, is not realistic and leads to a very unstable operation of the heat generators. Moreover, a still heat carrier fluid would probably lead to discomfort for the users. The value α = 0.3 was found to be a good compromise between the optimal control strategy of the heat generators and the heating needs of the DH customers.

Main assumptions
The DH network was simulated over two reference weeks that adequately represent the system behaviour for the entire year: • a winter week (WW), from February 27th to March 5th. The thermal load ranges between 7 and 30 MW t with an average value of 14.15 MW t . It is supplied by the cogeneration plants, which run almost always under full-load, the heat pumps, the heat recovered from the foundry and the gas boilers, which supply the remaining part of the heat load; • a week during middle season (MSW), from the 24th to 30th of April, in which the thermal load ranges from 2.5 to 10 MW t and its average value is equal to 4.8 MW t . The heat demand is provided by a single CHP (2.2 MW t ) always at full load, the heat pump connected to the operating CHP, the waste heat and the gas boilers.
The return temperature in the optimization problem can vary between 55 and 70 °C. That means that the average temperature difference between supply and return line may change between 10 and 25 °C. These values are realistic if compared to the temperature difference fluctuations recorded by the DHN operator at the supply station CCC. The mean temperature of the virtual HS ( T s ) can range between 95 and 75 °C. The aforementioned parameter k dem was set equal to 0.6. The CHP engines, which have a thermal efficiency of 38% in nominal condition, can be switched either on or off, and the same holds true for the corresponding heat pumps, with COP = 4. Therefore, they behave as "improved" CHP units with higher thermal efficiency and lower electrical efficiency compared to stand-alone internal combustion engines. Start-up and shut-down costs were set empirically so that the engines do not start up more than three times a day, as required by the DHS operator. The gas boilers are indicated here as a unique heat generator called GB, with a thermal efficiency equal to 0.9%. So, it can modulate its thermal output from 0.1 MW to 25.5 MW. The simulation time-step was set to 1 h.

Evaluation of the simulation results
The optimization problem has been first run without being coupled to the district heating simulation using 1-week horizon for both reference weeks. These results represent a benchmark scenario (OPT), i.e. what can be ideally achieved if the (simulated) system perfectly follows the heat supply schedules predicted by the optimization. Then, the optimization has been coupled to the detailed DH simulation model as described in Sect. 4.1. In the simulated scenarios the DH is operated with a model predictive controller (MPC) based on the same optimization problem. Both OPT and MPC results consider either the network as the only source of flexibility (net scenario) or the combination of both the network and the additional heat storage (net + hs scenario).
In order to evaluate the strategies, the operational costs were compared to the current ones under the same boundary conditions. The current costs depend on the current regulation explained in the description of the case study. The selling price of electricity was evaluated using the Italian national market price (PUN) profile of year 2017 (www. merca toele ttrico. org). In the WW the average PUN is 48.6 €/MWh and ranges between 30.0 and 66.9 €/MWh with significant daily cycles. A similar trend occurs in the MSW with the minimum price falling to 15 €/MWh. An excise duty a hp = 12.5 €/MWh is imposed to the electrical energy consumed by the heat pumps. The cost for the natural gas used for both the gas boilers and the internal combustion engines was assumed equal to approximately c gas = 32 €/MWh after the conversion. The waste heat recovered by the foundry is valued c wh = 20 €/MWh. This cost is discounted by the Italian energy efficacy certificates, which value is 250 € each. One TEE is obtained recovering about 7.7 MWh of waste energy. Furthermore, the following maintenance costs have been considered for the different generation plants: c man,chp = 7 €/MWh, c man,hp = 0.5 €/MWh, c man,gb = 1 €/MWh.

Results and discussion
In the first Subsection, the strategies adopted by the optimization model are shown and commented, for the case with the storage tank at the generation side. The second Subsection compares the results of the advanced control with those of the optimization model and of the basic control.

Optimization results during two representative days
Figures 4 and 5 show the forecasts (heat demand, waste heat availability and PUN), the heat supplied by the production plants to the network, the schedules calculated by the optimization, the return temperature of the water inside the pipes ( T r ) and the   Fig. 5 Simulation results for a typical winter day: a heat demand and heat supply; b heat supply mix; c average return temperature; d electricity price temperature inside the storage tank ( T s ) for two representative days,-i.e. Monday of both middle season (MSW) and winter weeks (WW). As far as the results of the typical middle season day are concerned, Fig. 4a compares the users heat demand ( Q dem ) with the thermal energy supplied to the district heating network ( Q supp ). As it can be seen during the first 5 h, the return water temperature decreases or increases depending on whether the energy, supplied by WH, GB and the heat storage tank (HS), is respectively lower or higher than the heat demand. From hour 4 to 8, when the electricity price starts to rise, the CHP increases the supplied thermal energy. Part of the energy produced by the cogeneration units and HP is stored in the tank, increasing T s , when CHP + HP is greater than HS, as it can be seen in Fig. 4b, c. In the same way the network is charged, because Q supp is higher than Q dem . From hour 8 to 14 the CHP production decreases in response to the reduction of the electricity price (Fig. 4d). From hour 14 to 15 the energy previously stored in the tank is used, and only WH and HS supply heat. From hour 15 the electricity price starts to increase again, the CHP increases the supplied thermal energy and the pattern seen from hour 4 to 15 is repeated. As regard the level of heat storage in the tank, at the end of the simulation period the average temperature inside the tank T s is equal to the initial one (Fig. 4c), so the constraint defined by inequality (16) is respected.
During the first 3 h of the typical winter day, the network is charged due to the higher thermal input in the network compared to the heat demand. Accordingly, in these hours, the water temperature inside the pipes increases as shown in Fig. 5c. Furthermore, in the first 4 h the electricity price-see Fig. 5d-is high enough to produce thermal power with a CHP unit, as shown in Fig. 5b. Instead, from hour 4 to hour 5 the electricity price is low, so the heat demand is satisfied by the GB. From hour 7 to hour to 23 all the cogeneration engines are running at full capacity and the GB are modulated to follow the energy demand. About the thermal power supplied by the storage tank, it mostly follows the production of CHP + HP (Fig. 5b), except from hour 6 to 10 in which it is lower, determining an increase of the temperature inside the tank ( T s )-see Fig. 5c-, and in the last 2 h in which is higher.

Simulation results of the model predictive control
The most critical part of the simulation framework is the coupling between the optimization model and the district heating network model, to form the model predictive controller (MPC). In general, the MPC output differs from that of the optimization model (OPT). This is because in the optimization model the network is considered as a one-capacitance model, whereas in the simulation the network behaviour is accurately represented using the detailed model NeMo. In order to match the outputs of the MPC to those of the optimization as much as possible, the mass flow rate has to be controlled as specified in Sect. 4.1. The most significant differences in the outputs of the two models are represented by the return temperatures ( T r ), shown in Fig. 6a, and the thermal power supplied by the storage tank to the network ( Q hs ) (Fig. 6b). Figure 6a, b refer to the representative middle season week. The difference between the simulation results compared to the solution of the optimization problem determines an increase of the operational cost as shown in Table 3. This is due to the oversimplification of the network thermal behaviour in the optimization constraints, which leads to errors in the state estimation, especially in the middle season week.
As mentioned above, the cost values determined by the optimization model in the different cases can be considered as a benchmark. The cost results of the MPC are between these benchmarks and the operational costs in the basic control (BAS).  In order to analyse more in detail the improvement determined by the advanced control, the scheduling of the basic control (Fig. 7a) is compared with those of MPC, considering respectively the thermal inertia of the network (the heat carrier fluid inside the network pipelines), in Fig. 7b, and of both the network and the storage tank (Fig. 7c). The period considered is the MSW.
The strategy adopted by the advanced control, Fig. 7b, is to supply heat with a single cogeneration unit and to increase the thermal energy production from CHP only when the value of PUN (the red line in Fig. 7) is higher than a threshold value (about 45 €/MWh). This determines a higher modulation of CHP units compared with the basic control case (Fig. 7a). When the electricity price exceeds the threshold value for a very short period of time (hour 30 and 128), extra cogeneration units don't turn on, to avoid an excessive number of start-ups and shut-downs in a limited period of time. In Fig. 7c, considering the thermal inertia of both the network and the storage tank, the thermal load profile is characterised by higher peak than in Fig. 7b. This is due to the greater thermal storage capacity of the system, which allows to produce more thermal power with CHP and HP, and to accumulate the production in excess in the storage tank, when the electricity price is high. On the other hand, when the electricity price is low, the stored energy can be supplied by the storage tank (blu dotted line in Fig. 7c). The energy production of the heat supply units and the corresponding share within the heat supply mix, which characterize the different type of controls, are shown in Tables 4 and 5, respectively for the MSW and WW.
In the MSW, the strategy adopted by the model predictive control allows a higher share of the heat produced by the CHP units compared to the basic control (from 31.8 to 49%). As a result, in this last case the MPC allows a reduction of heat produced by gas boiler (53.3% in the BAS, 32% in MPC under "net" scenario). The remaining share of the heat demand is covered by the heat pumps (10%) and by the waste heat sources (9%). Considering also the storage tank, the share of the heat produced by CHP is 37%, so we have an increase compared to the basic control and a decrease with respect to the "net" scenario.
On the other hand, in the WW the advanced control determines an increase of heat produced by GB and a reduction of the share supplied by CHP (from 62% in the BAS to 45% in the MPC). This is because, differently from the basic control, in this case the CHP plants are not in operation when the electricity price is lower than a threshold value, and when it is high enough the total installed power of the CHP plants is not sufficient to supply the entire heat demand. Here too the addition of the storage tank determines a decrease of this share to 41.5%. The heat supplied by the HP is characterised by the same trend of the CHP in both "net" and "net + hs" scenarios.
During the simulated periods, the proposed model predictive controller allows a reduction in the operational costs of 12.5% and 5.8% in the MSW and in the WW, respectively. This cost reduction is mainly due to the cost-optimal scheduling of the heat supply units, but also to a reduction in the heat supplied to the customers as shown in Tables 4 and 5. Considering also a greater thermal storage capacity of the system, determined by the tank, the cost savings increase up to 20% in the MSW and to 6.3% in the WW. So, in the representative middle season week, in which the total installed power of the CHP and HP plants is sufficient to supply the entire heat demand during the peak, the modulation of CHP and HP units has a higher impact, and the percentage increase of the cost reduction determined by the HS is higher.

Remarks on simulation results
The three main results and the related discussions are remarked in this section.
• In the optimization problem, the oversimplification of the network thermal behaviour results in an underestimation of the operational cost in comparison to the MPC, which considers a detailed district heating network model. This is true especially in the representative middle season week, when the cost result of the optimization problem is 8.5% lower than that evaluated by the MPC, in the scenario with the storage tank. • The cost reduction determined by MPC (exploiting the thermal inertia of the network) compared to the basic control is more relevant in the MSW (− 12.5%) than in the winter one (− 5.8%). This is in accordance with the results of a previous research article (Vivian et al. 2020b), in which it was demonstrated that in the coldest period the flexibility potential of the network is lower. • By adding a centralized storage tank, the cost savings increase much more in the MSW (20%) compared to WW (6.3%). The result is related to the case study under consideration, in which the total installed thermal power of the cogeneration engines and heat pumps, connected to the centralized storage tank, is limited in comparison to the peak load, during WW. In this period, the CHP and HP plants can't produce more energy compared to the basic control, when the electricity price is high. So, a small amount of energy is accumulated in the storage tank and his flexibility potential is not fully exploited.

Conclusions
The current paper presents a MILP optimization problem to schedule the heat supply units of the district heating network of Verona (Italy) based on heat demand, waste heat availability and electricity price forecasts. Coupling the optimization to the detailed district heating network model based on finite difference methods allowed repeating the scheduling with a rolling horizon scheme. The objective of the proposed control strategy was to minimize the operational costs for the DH operator using (1) only the thermal inertia of the water enclosed in the network pipelines as thermal buffer and (2) also using an additional thermal storage tank that could be possibly installed in the main supply station. The simulation in two reference weeks shows that the oversimplification of the network thermal behaviour in the optimization constraints leads to errors in the state estimation, which turns into higher costs compared to those initially predicted by the optimization, especially in the middle season week.
Comparing the results of the model predictive control to those of the current control strategy indicates a reduction in the operational costs of the system in the middle season and winter weeks of 12.5% and 5.8%, respectively. Increasing the flexibility of the system with an additional heat storage tank connected to the cogeneration units and the heat pumps these percentages increase to respectively 20% and 6.3%, respectively. Therefore, the benefit related to the additional thermal energy storage is more relevant in the warmest periods, when the total installed power of the CHP and HP plants is sufficient to supply the entire heat demand during the peak, and the modulation of these plants has a higher impact. Alternative scenarios could pursue an environmental objective instead of an economic one as shown here. Future work will try to improve the coupling between the optimization and the DH simulation, and it will introduce a detailed model of the thermal storage tank.