Supervisory control design for balancing supply and demand in a district heating system with thermal energy storage

This paper presents a systematic comparison between three alternatives to design the supervisory control layer of a district heating network composed of a waste heat boiler, an electric boiler, a dump, a hot water storage tank, and a set of consumers. The three alternatives are split range control, controllers with different setpoints, and model predictive control. We evaluate the closed-loop performance in the face of time-varying supply and demand, and constant electricity prices. All alternatives were found to give similar performance. Controllers with different setpoints is the easiest to implement, while model predictive control is the most diﬃcult.


Introduction
The storage of thermal energy is an effective solution to the problem of integrating intermittent heat sources such as solar thermal and industrial waste heat into district heating systems ( Lund et al., 2014;Miró et al., 2016;Guelpa and Verda, 2019 ). The most common form of storage currently used in district heating is the accumulation of hot water in tanks due to its low installation cost and high reliability ( Hennessy et al., 2019 ).
The optimal use of energy storage has been studied extensively in the operations research literature ( van de Ven et al., 2013;Harsha and Dahleh, 2014;Zhou et al., 2016;. Despite a few problems with special structure, most energy storage problems relevant in applications do not have a closed-form solution and a numerical solution by dynamic programming is impractical. Therefore, we resort instead to suboptimal, yet effective, policies. Economic model predictive control is an example of such policy ( Ma et al., 2009;Kumar et al., 2020 ). An alternative approach is to design a hierarchy of optimization and control layers that work independently on different time scales ( Skogestad and Postlethwaite, 2005 ). This is known as hierarchical control and is the approach we adopt in this paper.
Within the hierarchical control paradigm, model predictive control (MPC) has become the technique of choice for designing the supervisory control layer in most of the thermal energy storage systems reported in the literature due to its ability to handle constrained multivariable systems by design ( Cole et al., 2012 ). On the other hand, classical advanced PID-based control structures are widely used in practice De Oliveira et al., 2016 ). However, the systematic design and benchmarking of these control structures for thermal energy storage systems have received little attention in the literature.
We consider a district heating system with thermal energy storage. The system is assumed to be already designed and we only consider its operation. We also assume constant electricity prices. This paper extends the work of Zotica et al. (2020) , where a decentralized control structure based on split range control (SRC) and selectors was compared with model predictive control. In this paper, we propose an alternative decentralized control solution based on PI controllers with different setpoints. The advantage of SRC and controllers with different setpoints is that they handle active manipulated variable (MV) constraint switching without explicitly solving an optimization problem.
The main contribution of this paper is a systematic comparison of two decentralized control solutions (SRC and PI controllers with different setpoints) and MPC for designing the supervisory control layer for a district heating system with thermal energy storage. Fig. 1 shows the control hierarchy in a process plant. The control layer receives its setpoints (CV1 s ) from the upper optimization layer, and is divided into an upper supervisory control and a lower regulatory control. The latter handles control on the fastest time scale by controlling variables (CV2) that contribute to stabi- lizing the process, e.g. levels and pressures. It may also include tight control of economic controlled variables (CV1), usually active constraints. The regulatory control layer typically consists of PID controllers. These loops are usually not subject to reconfiguration, and therefore, extra consideration should be given to what happens when a constraint is reached, in particular, when a manipulated variable (MV) becomes saturated. Hence, the upper supervisory control layer should supervise the regulatory control layer. The main focus of this paper is the design of the supervisory control layer. It acts on a slower time scale, and its roles are ( Reyes-Lúa and Skogestad, 2019b ):

Control layers in a process plant
1. Keep the economic controlled variables (CV1) at the setpoints by using as degrees of freedom the setpoints (CV2 s ) to the regulatory layer or unused MVs. 2. Prevent saturation of the MVs in the regulatory layer. 3. Identify and take care of changes in active constraints.
The supervisory layer can be decentralized or centralized. The former implies several independent controllers that do not communicate with each other. Centralized control implies only one multivariable controller (e.g. MPC) that receives all the measurements from the layer below and simultaneously coordinates all the controllers in the regulatory layer ( Skogestad and Postlethwaite, 2005 ).

Classical decentralized control schemes to switch between active MV constraints
A decentralized supervisory layer that handles MV-MV switching can be designed using split range control, controllers with different setpoints, or valve position control ( Reyes-Lúa and Skogestad, 2019b ). The first two options are considered in this paper, and are described next.

MV-MV switching: split range control
Split range control is a multi-input single-output (MISO) control structure that extends the steady-state operating range for the CV by using a new MV when the initial MV becomes saturated. Split range control with three MVs ( u 1 , u 2 , and u 3 ) and one CV ( y ). v is the internal signal from the controller (C) to the split range block.    2 shows the block diagram for SRC with three MVs ( u 1 , u 2 , and u 3 ) and one CV ( y ). A feedback controller ( C), sends an internal signal ( v ) to the split range block (SR). This returns the values for the physical MVs ( u ). Fig. 3 shows a split range block for three MVs (e.g. u 1 , u 2 , and u 3 ) with different gains to the CV, but the same sign, and therefore three positive slopes ( α).
The split range block can be implemented in different ways, e.g. logic (if-else statements), lookup tables, or functions. In this work we use the last. SRC is commonly presented in process control textbooks ( Bequette, 2002;Seborg et al., 2003 ) and in industrial applications (e.g. Forsman and Adlouni (2018) ). However, despite its widespread use, only recently a systematic tuning procedure has been proposed ( Reyes-Lúa et al., 2019 ). This procedure adjusts the slopes ( α) in the split range block to get the desired controller gain for each MV considering the different dynamic effects of each MV on the CV.

MV-MV switching: controllers with different setpoints
Alternatively, controllers with different setpoints can be used for MV-MV switching ( Reyes-Lúa and Skogestad, 2019a ). The advantages are that the controllers can easily be tuned independently and no logic is needed. The disadvantage is that the setpoint is not constant and this will cause some delay during switching. Fig. 4 shows the block diagram for three independent controllers with different setpoints y s 1 , y s 2 , and y s 3 , that manipulate u 1 , u 2 , and u 3 , respectively, to control y . The order of activating the MVs is given by their selected setpoint difference (see Section 5.3.1 ).

Fig. 5.
Flowsheet of distribution network studied in this work with one waste heat boiler, one electric boiler, one air cooling (dump), and one hot water storage tank supplying hot water through a pipeline network to consumers.

Centralized control
Model predictive control is a unified systematic procedure for controlling constrained multivariable systems commonly used in industrial applications ( Qin and Badgwell, 2003;Mayne, 2014 ). At each sampling time, it uses the current plant measurement as the initial state to solve a finite-horizon open-loop optimal control problem to determine the optimal control sequence. Then, the first control is applied to the plant and the process is repeated at the next time step ( Mayne et al.,20 0 0 ). It handles constraints and interactive processes by design. However, it requires a detailed process model, which may not be available at the plant start-up.
The closed-loop performance of MPC depends firstly on the accuracy of the dynamic model and secondly on the choice of tuning parameters, e.g. weights in the objective function, prediction and control horizon, input rate constraints or constraints back-off ( Lu et al., 2020 ).
Finding the MPC tuning parameters is done by trial and error, heuristics, or optimization of indicators of the closedloop performance, e.g. overshoot, integral square error, robustness, etc. The work by Garriga and Soroush (2010) presents an overview of theoretical and practical guidelines for tuning the controller parameters in an MPC. The work by Lozano Santamaría and Gómez (2016) presents a gradient-based tuning algorithm with application to chemical processes. The work by Lu et al. (2020) presents derivative-free tuning algorithms based on Bayesian optimization techniques. In this work, we use trial and error.

District heating control problem
Fig. 5 shows the thermal energy storage system analyzed in this work, where the working fluid is hot water. For example, this can be a district heating network supplying hot water to residential households. This system has a direct physical connection between supply and demand, such that hot water can be directly sent to the consumers bypassing the storage tank (flowrate q SP ). There are also other distribution networks, for example, in industrial clusters, that exchange energy only through the storage tank ( Knudsen et  The operational objective of the network in Fig. 5 is to manipulate the inputs u to minimize the electric boiler usage ( u 2 = q EP ) ( Eq. 1a ), while balancing the supply and demand ( Eq. 1b ). Furthermore, the storage capacity constraints ( Eq. 1c ) and the model equations (see Section 4 for details) must be satisfied.
The four inputs (degrees of freedom) in Fig. 5 In the system, heat is mainly supplied by the waste heat boiler ( q S ), which extracts heat by burning waste. This may be viewed as a disturbance ( d ) to the system along with the consumer demand ( q ), i.e. d = [ q S q ] . Note that the inputs u correspond to physical valves in Fig. 5 . Later, we will make use of some transformed inputs (MVs) for the purpose of balancing supply and demand. Based on process insight, the four degrees of freedom can be used to balance supply and demand as follows: 1. Excess supply: charge hot water to storage ( q ST ). 2. Excess supply (when storage is full): dump hot water to air ( q D ). 3. Excess demand: discharge hot water from storage ( q T P ). 4. Excess demand (when storage is empty): use electric boiler ( q EP ).
The switching between these four operating regions must be taken care of by the supervisory control layer. Given that electricity prices are assumed constant, this corresponds with the optimal storage policy. That is, this is the optimal solution to Eq. 1 . If electricity prices were time-varying, it would be optimal to store energy at low prices to be used later when prices are high, but this is beyond the scope of this work.
The main objective of the supervisory control system is to meet the energy demand of the consumers by switching between the four operating regions. To simplify the design of the supervisory control system, we will consider three MVs in this layer, rather than the four physical valves (degrees of freedom) shown in Fig. 5 : MV1: hot water from waste heat boiler ( q SP ). MV2: hot water from storage tank ( q T P ). MV3: hot water from electric boiler ( q EP ).
The motivation of selecting the three MVs is that they are the three suppliers of hot water to the consumers and that these MVs were selected for control design in an actual district heating system. Here, MV1 is the flow in the direct physical connection from the variable supply to the consumers. Note that MV1 does not correspond to a physical valve, but it is indirectly given by the material balance in Eq. 2 .
Here, the supply q S is a disturbance, whereas the charge q ST and dump q D are physical valves. For cases where we want q SP to be smaller than q S (that is, we have excess supply), we first charge the tank ( q ST ) and then, when it is full, we start dumping ( q D ). We will assume that this logic is taken care of by a separate block "charging policy logic", which will be part of the regulatory control system (see Section 5.1 ). As mentioned, the reason for doing this is to simplify the design of the supervisory control system.
In summary, the three MVs are related as follows to the physical inputs ( u ): Note that the maximum values for MV1 = q SP and MV2 = q T P are time-varying. That is, q SP is limited by the hot water from the waste-heat boiler ( q SP ≤ q S ), and q T P = 0 when the storage tank is empty.
The temperature of the hot water produced in the waste heat and electric boiler is kept constant by water injection, not shown in Fig. 5 . This is a common practice in district heating networks. There are also other flows not shown in Fig. 5 , for example, the water supply to the electric boiler and the water return from the air cooling. Actually, in Fig. 5 , it may be better to view the varying water flows q [ m 3 h −1 ] as being energy flows Q [ Jh −1 ]. Because of the assumption of constant temperature, q and Q are directly proportional: Q = kq .
In Zotica et al. (2020) , the information about the heat demand was assumed to be available. In practice, this is not realistic, and we instead control the network pressure ( p), which is proportional to the mass m (see Fig. 5 ) in the pipeline. This is a dynamic variable that couples the supply and demand. Therefore, it is an indirect and reliable measurement of the supply-demand balance in a water distribution network.
In summary, Table 1 shows the three MVs, one CV and the two main disturbances (DVs) for the supervisory control considered in this work.

Process model
The change of mass ( m ) in the pipeline system is given by the mass balance in the network, Eq. 4 . dm dt where, ρ [ kg m −3 ] is the water density, assumed constant.
To model the changes in the network (pipeline) pressure ( p), we consider that it is proportional to the change of network mass ( m ), Eq 5 .
where, m [kg] is the network water mass, m 0 [kg] is the water mass at the reference flow, is the constant compressibility coefficient and p 0 [bar] is the pressure at the reference flow. Hence m = ρV, the compresibility factor takes into account the increase in liquid density and more importantly the increase in pipeline volume by increasing the pressure.
Substituting Eq. 5 into Eq. 4 yields the mass balance expressed in terms of the network pressure, Eq. 6 . dp dt = To model the storage tank inventory, we neglect changes in density ( ρ) and assume no heat losses. The dynamic mass balance for the storage tank is given by Eq. 7 where, V h [ m 3 ] is the volume of the hot water, which must be within the limits V min = 0 and V max .
With constant inlet temperature ( T S ), constant heat capacity ( c P ) and perfect mixing, the tank temperature ( T h ) is constant and equal to the inlet temperature ( T S = T h ). Table 2 shows the model parameters.

Dynamic behaviour
The dynamic behaviour of the model is analyzed from step responses in the disturbances and inputs. MV3 = q EP has the same effect on p as MV1 and MV2 and is not shown. Fig. 6 shows the response for a step increase in MV1 = q SP = 250 m 3 h −1 given by an increase in the available hot water supply (DV1 = q S = 250 m 3 h −1 ). Fig. 6 c shows the network pressure response which is an integrating process as given in Eq. 5 with slope k = 3 . 33 bar m −3 . The hot water volume ( Fig. 6 b) is constant. Note that the time scale is in seconds because the pressure dynamics are fast. Fig. 7 shows the response to an increase in the discharge MV2 = q T P = 250 m 3 h −1 . The hot water volume ( Fig. 9 b) is an integrating process with initial slope k = −1 . Fig. 8 shows the response for a step increase in the available hot water supply (DV1 = q S = 250 m 3 h −1 ) with MV1 = q SP constant. The "charging policy logic" first sends the excess hot water to the storage tank, and once the tank is full, it is dumped ( Fig. 8 a). The time scale is in hours because the storage tank dynamics are slow. Fig. 8 b shows the hot water volume. The network pressure ( Fig. 8 c) is constant because MV1 = q SP is constant.
Finally, Fig. 9 shows the response to an increase in the demand DV2 = q = 250 m 3 h −1 with the MVs constant.

Control system design
We want to implement a control system that optimally switches between the four operating regionsr described in Section 3 .

Regulatory control: charging policy logic
As mentioned before, to simplify the design of the control system we include the charging policy logic in the regulatory control system. The logic is:    1. The storage tank is charged with excess hot water ( q ST ) when the hot water storage is below maximum capacity, Eq. 8 .
2. On the other hand, when the storage tank is full, excess heat is dumped, Eq. 9 .

Alternative 1 for supervisory control: split range control
We first consider split range control to keep the network pressure at the setpoint by using one MV at a time, starting with the cheapest, and switching to the more expensive as demand increases or the availability of the cheap MV decreases. In our case, we first want to use the available hot water from the waste heat boiler (MV1 = q SP ), followed by the hot water stored in the tank (MV2 = q T P ), and lastly the electric boiler (MV3 = q EP ). Fig. 10 shows the SRC implementation.

Tuning parameters for SRC
We follow the procedure of Reyes-Lúa et al. (2019) . We define the normal range for the internal signal v to be from 0 % to 100 % , and we scale the MVs from 0% to 100% , see Fig. 11 . The tuning parameters for SRC are the PI-tunings for the common controller C and the slopes α i in the split range block in Fig. 11 . The slopes α i are used to allow for different controller gains for each MV, however, from the network mass balance ( Eq. 4 ), all MVs have the same effect on the CV. Therefore, the three slopes in the split range block are equal, and we get α SP = α T P = α EP = 3 . Fig. 11 shows the split range block. The controller parameters are obtained by applying the SIMC tuning rules ( Skogestad, 2003 ) for an integrating process (i.e. the network pressure balance Eq. 6 ). We select the closedloop time constant τ C = 10 s, resulting in the proportional gain for each MV is K C,i = 54 and the integral time τ I = 8 s (see Section 7.4 ).
We find the common controller gain K C = K C,i = 54 /α i = 18 . To handle the time-varying availability of MV1 and MV2 we update the PI-controller bias ( u 0 ), such that when one MV is no longer available, the new MV starts from the value of the former. Alternatively, this could have been handled by using a "discharging policy logic" block and only two MVs for the SRC (see Section 7.2 ). Note that we have assumed that we directly manipulate the flows q i , that is, we have assumed that all valves have flow controllers. Without flow controllers, we would have had to use different slopes α i in the split range block ( Fig. 11 ). Fig. 12 shows the control structure with three different PIcontrollers with different setpoints for controlling the network pressure ( p) that uses as degrees of freedom MV1 ( q SP ), MV2 ( q T P ) and MV3 ( q EP ). Similar to SRC, we order the use of MVs based on economics, and we use the cheapest MV first. Therefore, we order the three setpoints SP 1 > SP 2 > SP 3 ( p s SP > p s TP > p s EP ) such that only one MV is actively used at any given time.

Setpoints selection
We select the setpoints order based on physical insight. The process gain from the MVs to the CV is positive ( Fig. 6 ). Therefore the controller gain is positive, and a negative controller error ( p s − p) gives a negative controller output (see Eq. 13 ). Setting the controller bias u 0 = 0 , and considering that the minimum physical limit for the MV is 0, the MV only starts to open when the controller error becomes positive. Specifically, when p ≥ p s T P , MV1 is active, and MV2 and MV3 are fully closed. Once MV1 reaches its maximum limit, supply is smaller than demand, and the network pressure drops. Once the pressure reaches a lower threshold, p s EP < p < p s T P , MV2 becomes the active MV, while MV3 is fully closed. Finally, when the supply is smaller than demand, MV2 becomes saturated at its maximum and the pressure drops. Once it reaches an even lower threshold p < p s EP , MV3 becomes the active MV. This control structure handles the intermittent availability of MV1 and MV2 by design as long as antiwindup with tracking of the plant input is implemented.

Tuning parameters for controllers with different setpoints
As

Alternative 3 for supervisory control: model predictive control
We design the MPC to handle also the charging of the storage tank. Therefore it has the role of the supervisory and regulatory layer previously described. We include q D as the fourth MV, while q ST is calculated from the mass balance ( Eq. 2 ).
The system we are analyzing is somewhat atypical because it has more MVs than CVs. Therefore, tuning the MPC is not straightforward, and we must give careful consideration in setting up the objective function to prioritize the use of MVs. We achieve this by selecting the weights ( ω) in the objective function ( Reyes-Lúa et al., 2018 ). We formulate the optimal control problem with the objective function given in Eq. 10a . We want to maximize discharging the tank ( q ST ), minimize dumping ( q D ), minimize using the electric boiler and keep the network pressure ( p) at its setpoint ( p s ). As mentioned before, the MPC controls the network pressure as an indirect measure of the hot water demand. However, the MPC uses the full model ( Section 4 ), and it requires information about the demand hot water. We solve the optimization problem subject to model Eqs. 10b, 10d and 10c , and operation constraints (Eqs. 10e, 10f, 10g and 10h ). Here, k is the current iteration, N is the number of control intervals, ω i are the weights in the optimization problem, and Eqs.
10b and 10c are discretized versions of the mass balances Eqs. 6 and 7 , respectively. We formulate the MPC problem using CasADi ( Andersson, 2013 ), and we use IPOPT to solve the NLP ( Wächter and Biegler, 2006 ). The tuning parameters were found by trial and error. We use 10 control intervals, 10 min prediction horizon and 1 min sampling time. In practice, the sampling time will need to be smaller because of the fast pressure dynamics. The weights in the optimization function were also found by trial and error and are ω T P = 10 −6 , ω D = 10 −5 , ω EP = 10 −3 and ω p = 10 4 .

Simulation case study
We compare the performance of the three control system alternatives to switch between the four operating regions (see Section 3 ) using the model described in Section 4 .
At the initial state of the system, the tank is half full ( V h = 2500 m 3 ), and the hot water supply from the waste heat boiler is equal to the demand, i.e. q S = q = 500 m 3 h −1 .

Simulations step changes
We perform the following series of step changes (see also Fig. 13 ): Step 1 at time t = 1 h : The hot water supply from the waste heat boiler (DV1) increases from q S = 500 m 3 h −1 to q S = 10 0 0 m 3 h −1 .
Step 2 at time t = 12 h : The hot water demand (DV2) increases from q = 500 m 3 h −1 to q = 10 0 0 m 3 h −1 . Step 3 at time t = 15 h : The hot water supply from the waste heat boiler (DV1) decreases from q S = 10 0 0 m 3 h −1 to q S = 500 m 3 h −1 .

Results analysis
The closed-loop simulation results in Fig. 14 demonstrate that all three alternative control structures can successfully implement optimal operation for this system. Because the simulation time scale is in hours, it is not easy to see some of the differences between the three alternatives for MV3 = q EP and for the network pressure ( p). SRC and MPC perform similarly. In the MPC formulation, the step changes in disturbances happen at the sampling time and there is feedforward action from disturbances q s and q at each sampling time. Therefore, there is less variation in the network pressure p in Fig. 14 i compared to Fig. 14 g. This is shown more clearly at time t = 12 h in Fig. 15 . The pressure response for controllers with different setpoints ( Fig. 14 h) is as expected different from SRC ( Fig. 14 g) and MPC ( Fig. 14 i). However, the response for the tank storage ( Fig. 14 e) and MVs ( Fig. 14 b) is not significantly different because 1) the pressure only has a dynamic effect on the flows, that is, the flow values are independent of the pressure setpoint and 2) pressure dynamics are fast compared to the storage dynamics.
In Fig. 14 a, SRC changes q T P instantaneously because we update the controller bias. However, for controllers with different setpoints (3C), there is a small delay until the new MV takes over, ent setpoints (3C) and MPC. Among the three, the alternative of controllers with different setpoints shows the highest input usage. The total variation for q EP is equal for SRC and MPC, meaning that the electricity cost is equal. Compared to SRC, MPC uses slightly less q SP compensated by using slightly more q T P and marginally dumping more heat q D . This means that in the MPC implementation, the tank is simultaneously charged and discharged. However, the short-term peaks in q EP for controllers with different setpoint (3C) which do not matter much for the integrated cost, see Eq. 1 .

Ease of implementation
In terms of ease of implementation, the use of three controllers is the simplest. It allows for using three independently tuned controllers and it avoids the logic needed in SRC. The logic can be avoided because the switching is done based on the output (CV = p) and not the limit on the MV-value. However, it has two disadvantages: 1) somewhat worse dynamic performance (see Table 3 ) and 2) varying setpoint. Because each controller can be tuned independently, we do not need to compromise on the integral time as in SRC. This was not relevant for the process studied because the MVs have the same effect on the CV (see the model in Section 4 ), but it can become important for other systems with different dynamics, for example with different valve sizes.
The MPC controller is by far the most difficult to implement. In addition to requiring a dynamic model, it was also difficult to adjust the tuning parameters (e.g. weights in the objective function, prediction horizon and sampling time) to obtain the desired performance. We selected a short prediction horizon because with larger values, hot water was dumped before the tank was at maximum storage capacity. Formulating the objective function ( Eq. 10a ) was also done by trial and error. In the end, the weights on q T P and q D were added to give information to the MPC about which flowrates to prioritize. Without these, the MPC would use two MVs simultaneously. For example, with no excess heat, it would discharge the storage tank instead of using hot water from the waste heat boiler, and dump the remaining hot water. This was because the MPC does not have information about future demand in its prediction, and it is not aware that it should charge the storage tank. One could also add a penalty for not charging the storage tank in the objective function. However, as with any multiobjective problem, there will be a compromise and the tank will not be charged to maximum capacity.

General control structure for balancing supply and demand
A general solution for balancing variable supply and variable demand is shown in Fig. 16 . This is an inventory ( m ) control problem with two MVs (MV s and MV d ) and two DVs ( d 1 and d 2 ), and The idea is that the adjustable supply MV s should be used when the variable demand ( d 2 ) is larger than the base load supply ( d 1 ), and the adjustable demand (MV d ) should be used when d 1 > d 2 .
One should normally not use MV s and MV d simultaneously. The values of MV s and MV d are set by a feedback controller that controls the inventory m that indirectly measures the imbalance between supply and demand (i.e. pressure in our case), and the MV s -MV d switching is taken care of by, for example, SRC . Fig. 17 shows the split range block At the split value ( v * ), d 1 = d 2 and the variable demand balances the variable supply. Alternatively, we may use two controllers for MV s and MV d with two different setpoints for the inventory m .
Next, we need to decide on how to implement MV s and MV d using the physical inputs ( u 1 , u 2 , u 3 and u 4 , in our case). From Fig. 16 we have MV s = q from storage + q electric boiler = u 2 + u 3 (12a) MV d = q to storage + q dump = u 1 + u 4 (12b) Since our objective is to minimize the use of electric boiler (and minimize dump to store when possible), it becomes clear what we should do. If MV s is active, we first use hot water from storage ( u 3 ) until the storage tank is empty, and then use the more expensive electric boiler ( u 2 ). If MV d is active, we first send excess hot water to storage ( u 4 ) and when the storage tank is full, send it to dump ( u 1 ). This corresponds to the four operation regions in Section 3 . In practice, this may be implemented using a "charging policy" logic for MV s and a "discharging policy" for MV d

Setpoint difference for three controllers
We select the setpoints for the three controllers in Section 5.3.2 by trial and error. This choice is a trade-off. A smaller setpoint difference may result in having more than one MV active at a given time, while a larger setpoint difference results in a larger delay until the next MV activates.

PI Controller tuning
The PI-controller is given in Eq. 13 .
where u 0 is the controller bias, and K C is the proportional gain and τ I is the integral time constant, derived from applying the SIMC tuning rules ( Eq. 14 ) ( Skogestad, 2003 ).
where k is the initial slope of the step response, θ is the time delay, τ is the time constant, and τ C is the desired closed-loop time constant. The PI-controllers in this paper are for integrating process, when τ → ∞ .

Conclusion
The use of inventory (pressure) control is an effective way of balancing supply and demand for the district heating system. For the case of constant electricity prices, optimal operation for the system studied is easy to identify based on physical insight (see Section 3 ). In this work, we compare three alternative control implementations (split range control, controllers with different setpoints, and model predictive control) to handle MV-MV switches. The closed-loop simulation results in Fig. 14 show that all control structures successfully switch between the four operating regions to balance supply and demand. However, MPC requires careful tuning to obtain the desired performance, making it more difficult to implement than the decentralized solutions.

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.