Peak Shaving in District Heating Utilizing Adaptive Predictive Control

Abstract

District heating systems (DHS) are driven by the heat demands of their consumers, with higher demands giving a higher load on the heat production. While heat demands are human-dependent, they contain diurnal behaviors and weather dependencies. The diurnal behaviors contain periods with high demands causing peak loads on the heat production, which is operationally costly. This is especially true for heat pumps, a solution for DHS to include green energy, as the cost depends directly on the needed temperature. This paper presents a formulation of adaptive model predictive control (MPC) for inducing peak shaving on the production load to handle the peak load problem by using the DHS distribution network as a heat storage. It also presents a simulator model to describe the DHS. The MPC was applied to data from a case study of the DHS in Brønderslev, Denmark, showing a peak reduction of around $8\%$.

1. Introduction

The heating consumption of houses, industries, cities, etc., all have diurnal behaviors [1], daily periods where their demand is high and periods where it is low. This results in peak loads on the heat production, where the peak loads are often more expensive to produce due to the higher temperature needs during the peak load as well as the production unit mix used to cover the peak.

District heating systems (DHS) is particularly affected by these peak loads, as the production is far from the consumption, giving a delay between temperature increase and its distribution, inducing higher inflows to fulfill demands, thus raising the supply power and production costs. The importance of this issue is likely to increase as heat pumps are being integrated into DHS as part of the transition to green energy, with their production cost being directly dependent on the needed supply temperature.

A DHS is a network for the distribution of heat to an entire district [2], see Figure 1, utilizing large-scale heat production instead of independent production per consumer, as in the case of gas. DHS is generally controlled through the supply temperature, while the inflows are governed by differential pressure. The heat is distributed through a medium, usually water, using the pipes of the network. While the pipe loses heat to the surroundings, it also stores the water until it reaches the consumer. Given the diurnal nature of the demand, it is possible to forecast when the demand is expected to increase [3].

To ensure that the supply temperature is high enough to satisfy the demand, but not higher than is needed, model predictive control (MPC) methods have previously been utilized, in combination with forecasts of heat demand [4,5]. MPC is a well-studied control method [6,7], capable of dealing with prediction and constraints. It has successfully been applied to multiple fields, such as cement production [8], urban drainage [9], district heating [10], and more [11].

To address the problem of peak load on the production, this paper contributes with the formulation of an MPC design to minimize the peak loads, using the pipe system to store heat energy. The main idea is to use the preheating of the water entering the pipe system prior to the peak load occurring, utilizing the transport delay as an opportunity for heat storage. The stored water in the pipe system then has a higher temperature, from which the coming load can draw power, thus reducing the need for increasing the supply power at the time.

Additionally, a surrogate model for the DHS is formulated and calibrated for the case of the district heating system in Brønderslev, Denmark. Brønderslev is a city of approx. 12,500 people, where the DHS supplies 5000 out of 5500 houses with heating [12]. This model is used to evaluate the proposed controller.

In the following sections, we start with outlining the surrogate model, then the MPC design, and end with the evaluation of the design.

2. Modeling of District Heating

District heating systems consist of a network of pipes. They can in general be considered to consist of two main pipes, one supplying and one returning, and branches with subnetworks as illustrated in Figure 1. The areas in these branches can then be lumped together as single houses, giving a simple lumped model of houses along a pipe.

As the return temperature of each house can be assumed roughly the same, our model assumes the temperature of the return pipe to be the same along the pipe, simplifying the model to a single pipe.

As pipes in district heating are constantly full, we describe the pipe using volume elements in discrete time, such that in each period, the volume leaving the pipe equals the volume entering, as illustrated in Figure 2. Each element is then described by its temperature, volume, and position in the pipe, while being updated as it traverses through the pipe. This gives a discrete description of the stored heat energy and its distribution along the pipe.

The model includes two updates for the elements: (1) the contained volume; as the element passes by houses, flows to the houses are subtracted from the volume by (1). (2) Heat loss to the environment [13]; as an element moves through a pipe section, heat is lost as (2), and illustrated in Figure 3.

$$\begin{array}{cc}\hfill q_{house}& =\frac{Q_{house}}{C_p\rho (T_H-T_r)}\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill Q_{loss}& =\frac{2\pi L_p}{ln\left(4\frac{d_p}{D_p}\right)}{\sigma }_p(T_{in}-T_g),T_L=\frac{Q_{loss}}{C_p\rho q_{in}}\hfill \end{array}$$

(2)

where $T_r$ is the return temperature, $T_H$ is the temperature going into the house, $T_g$ is the ground temperature, $T_{in}$ is the ingoing temperature of the section, and $q_{in}$ is the inflow to the section. The remaining parameters are the section length $L_p$, the section depth $d_p$, the section diameter $D_p$, the section’s thermal conductivity ${\sigma }_p$, the heat loss $Q_{loss}$, the house flow $q_{house}$, and the heat demand $Q_{house}$.

In district heating, the pipe dimensions decrease along the network; by splitting the pipe model into sections between houses, the dimensions can be represented as seen in Figure 4.

2.1. Model Dimensioning—Brønderslev

For the particular model discussed in this paper, we utilized data and forecasts from the district heating system in Brønderslev, Denmark. For the dimensioning of said model, we considered the static case, where the flow velocity v is constant throughout the pipe, giving a constant ratio between static flow $q_i$ and cross-area $A_{p,i}$:

$$\begin{array}{c}\hfill v=\frac{q_i}{A_{p,i}}=\frac{q_{i+1}}{A_{p,i+1}}\end{array}$$

(3)

The model was calibrated assuming a total temperature loss $T_{loss}$ of 5 °C from the entrance to the last house, a typical value for a Danish DHS during winter operation. For the calibration, the static supply temperature, return temperature, and heat load were determined as the average from the available data.

Additionally, we assumed that the total transport time $t_{trans}^{Total}$ through the pipe system was 4 h, that there were 10 houses ($N_{house}$) along the pipe (branches), and that they were equally spaced with equal heat demands. Each pipe section was assumed to have equal temperature loss. With an assumed constant velocity, the transportation time, loss, and pipe temperature of each pipe section are given by:

$$\begin{array}{cc}\hfill t_{trans,i}& =\frac{t_{trans}^{Total}}{N_{house}},T_{L,i}=\frac{T_{loss}}{N_{house}},T_{pipe,i}=T_{sup}-\sum _{j=1}^i{T}_{L,j}\hfill \end{array}$$

(4)

Assuming each house has equal heat demand, the flows to each house, inflow, and flow through each pipe section are given by:

$$\begin{array}{cc}\hfill q_{house,i}& =\frac{\frac{1}{N_{house}}{Q}_{load}}{C_p\rho (T_{pipe,i}-T_r)},q^{in}=\sum _{i=1}^{N_{house}}{q}_{house,i},q_{pipe,i}=q^{in}-\sum _{j=1}^{i-1}{q}_{house,j}\hfill \end{array}$$

(5)

With the flows, temperatures, and transport times determined, the dimensions of the pipe model can be computed according to the mentioned assumptions.

$$\begin{array}{ccc}\hfill A_{p,i}& =A_{p,i-1}\frac{q_{pipe,i}}{q_{pipe,i-1}},\hfill & \hfill D_{p,i}=\sqrt{\frac{4}{\pi }{A}_{p,i}},d_{p,i}=d_{p,i-1}+D_{p,i-1}-D_{p,i}\end{array}$$

(6)

$$\begin{array}{ccc}\hfill L_{p,i}& =\frac{q_{pipe}{t}_{trans,i}}{A_{p,i}},\hfill & \hfill S_i=\frac{2\pi L_{p,i}}{ln\left(\frac{4d_{p,i}}{D_{p,i}}\right)},{\sigma }_{p,i}=\frac{C_p\rho q_{pipe,i}}{S_i}\frac{T_{L,i}}{T_{pipe,i}-T_g}\end{array}$$

(7)

where $q_{pipe,0}=q^{in}$ and $T_{pipe,0}=T_{sup}$.

2.2. Model Computation

As discussed, the model consisted of a series of volume elements describing the state along the pipe. Each element tracked its current/previous position, its temperature, volume, and its initial temperature $T_0$ when entering the current section.

At each time step, the simulator determined the needed inflow, by computing the heat consumption of the houses based on their demands $Q_H$ and the pipe and supply temperatures. A new element was then added to the pipe, filling the pipe volume, while empty ones were removed. The heat loss of each element was then computed based on the element’s translation through the pipe. The details of the consumption and heat loss are discussed in the following parts, where $t_s$ is the sampling time of the discrete model.

2.2.1. Consumption

The consumption of a given house was computed by subtracting the needed volumes from the nearest pipe elements at or upstream of the house. Starting with the house at the end of the pipe, the consumption was computed independently for each house. From each element, only the volume part that was upstream of a house was used. If the positions of the element and the house are given in volume along the pipe, then the available needed volume of the element to fulfill the heat demand of the house is given by

$$\begin{array}{c}\hfill V_H:=min(\frac{Q_H}{C_p\rho (T_{el}-T_r)},V_{el},V_H^{pos}-V_{back})\end{array}$$

(8)

where $V_H^{Pos}$ is the house position, $V_{back}$ is the back position of the element (upstream end), $Q_H$ is the current demand of the house, and $V_{el}$ and $T_{el}$ are the volume and temperature of the element, respectively. Correspondingly, the element’s volume and the demand of the house are then reduced:

$$\begin{array}{cc}\hfill V_{el}& :=V_{el}-V_H,Q_H:=Q_H-\frac{V_H}{t_s}{C}_p\rho (T_{el}-T_r)\hfill \end{array}$$

(9)

If the demand is unfulfilled (not zero), the procedure is repeated for the next element upstream of the house. When demand is fulfilled after, e.g., n elements, the flow and temperature supplied to the house for satisfying the demand are given by

$$\begin{array}{c}\hfill q_H:=\frac{1}{t_s}\sum _{i=1}^n{V}_{H,i},T_H:=\frac{{\sum }_{i=1}^N{T}_{el,i}{V}_{H,i}}{{\sum }_{i=1}^n{V}_{H,i}}\end{array}$$

(10)

The necessary inflow is then obtained after repeating the above computation for each house:

$$\begin{array}{c}\hfill q_{in}:=\sum _{i=1}^{N_{house}}{q}_H\end{array}$$

(11)

2.2.2. Computations of Losses

The heat loss of the elements can then be determined from the element’s change in position. As each pipe section has different dimensions, the heat losses are computed in parts, starting from the previous position $V_{pos}^{past}$ to the nearest passed house $V_{H,i}^{pos}$, then to the next passed house until the current position $V_{pos}$ is reached.

$$\begin{array}{c}\hfill V_{pos}^{past}\le V_{H,i}^{pos}\le V_{pos}\end{array}$$

(12)

If $V_{part}$ is the end of the considered part of section i (the smallest of $V_{H,i}^{pos}$ and $V_{pos}$), then the procedure is given by Algorithm 1.

Algorithm 1 Heat Loss of Element.

Require: i = section # of past position

$V_{part}:=min(V_{H,i}^{pos},V_{pos})$

while $V_{pos}^{past}<V_{pos}$ do

(1) Determine the moved amount:

$\Delta v:=V_{part}-V_{pos}^{past}$

(2) Compute the time spent and the moved length:

$\Delta t:=\frac{\Delta v}{V_{el}}{t}_s\Delta L:=\frac{\Delta v}{A_i}$

(3) Determine the heat loss, where $T_0$ is the temperature when entering that section.

$Q_{Loss}:=\frac{2\pi \Delta L}{log\left(4\frac{d_i}{D_i}\right)}{\sigma }_i(T_0-T_g)$

(4) Apply the heat loss and update the temperature:

$\Delta q:=\frac{A_i\Delta L}{\Delta t},T_{el}:=T_{el}-\frac{Q_{loss}}{C_p\rho \Delta q},T_0:=T_{el}$

(5) Update the positions (past, part and section # i ):

$V_{pos}^{past}:=V_{part},V_{part}:=min(V_{H,i+1},V_{pos})$

$i:=i+1$

end while

When the heat losses have been applied, then the diffusion between elements can be computed for each element:

$$\begin{array}{c}\hfill T_{el}:=(1-r)T_{el}+r\frac{T_{el,up}+T_{el,down}}{2}\end{array}$$

(13)

where r is the mixing rate (a mixing rate of $r=0.1$ was used throughout this work), with $T_{el,up}$ and $T_{el,down}$ being the temperature of the elements up- and downstream of the specific element.

2.2.3. Model Simulations

To test the model, simulations using historical data were used. As flow was computed by the model, we could compare the computed flow to the historical flow to evaluate the model. The used data of the heat demands were from the utility company’s measurements of the overall demand/consumption for the period of the data. In Figure 5, the flow of two simulations (dimensioning with 0 °C and 5 °C of total temperature loss) is given next to the historical flow. It can be seen that the 0 °C simulation is the closest match to the historical data, while the 5 °C induces slightly higher flows. Both simulations occasionally experience negative overshoots, as the approximation using the lumped model diverges from the historical data, due to the sampling time and the size of the elements.

As mentioned earlier, the 5 °C of total temperature loss was used in the dimensioning of the model in the later sections; this was chosen to ensure the simulation with MPC would have a temperature loss to consider, given the heat load data of the system included not just the actual heat load of houses but also the heat losses.

3. MPC Design

For the MPC design over a prediction horizon $N_H$, a dynamic description is needed to describe the relationship between the control variables and the objectives. As this work focused on peak shaving of supply power, we utilized the definition of the supply power $P_{in,t}$ as a base for the design:

$$\begin{array}{c}\hfill P_{in,t}=C_p\rho q_{in,t}(T_{s,t}-T_{r,t})\end{array}$$

(14)

depending on both the inflow and the supply temperature. The inflow in return relied on the heat demands and the house temperatures throughout the pipe:

$$\begin{array}{c}\hfill C_p\rho q_{in,t}=\sum _{i=1}^{N_{house}}\frac{Q_{H,i,t}}{T_{H,i,t}-T_{r,t}}\end{array}$$

(15)

and the house temperatures were similarly affected by past loads and operations, thus coupling the supply powers at different times.

3.1. System Identification and Models

In the simulations, we used historical data and forecasts from the DHS in Brønderslev. The available data included the supply temperature, return temperature, and total heat load (losses + demand), as well as forecasts of load and supply temperature. The data and forecasts were provided by ENFOR, using methods for heat load forecasts [3] and supply temperature forecasts [14].

As the data and forecasts described the system externally, and the detailed system described above was highly nonlinear, given the house temperatures’ interconnection, we applied an adaptive approach to keep the design model simple, utilizing system identification (SysId) through the onlineforecasting package in R [15,16]. More specifically, we used recursive fitting [17] to obtain a prediction model based on historical forecasts and data, updating the model parameters when new data were available.

In this work, we tested five model structures for the MPC design, see (16)–(20),

$$\begin{array}{cc}\hfill \mathrm{Model}1:& P_{in,t+k|t}={\theta }_{1,k}lp\left(\frac{T_{s,t+k|t}-T_{r,t+k|t}}{{\tilde{T}}_{H,t+k|t}-T_{r,t+k|t}},a_1\right)Q_{+k|t}\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \mathrm{Model}2:& P_{in,t+k|t}={\theta }_{1,k}lp\left(T_{s,t+k|t}-T_{r,t+k|t},a_1\right)Q_{t+k|t}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \mathrm{Model}3:& P_{in,t+k|t}={\theta }_{1,k}lp\left(T_{s,t+k|t},a_1\right)Q_{t+k|t}+{\theta }_{2,k}lp\left(T_{r,t+k|t},a_2\right)Q_{t+k|t}\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \mathrm{Model}4:& P_{in,t+k|t}={\theta }_{1,k}lp\left(\frac{T_{s,t+k|t}-T_{r,t+k|t}}{{\tilde{T}}_{H,t+k|t}-T_{r,t+k|t}},a_1\right)+{\theta }_{2,k}lp\left(Q_{t+k|t},a_2\right)\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \mathrm{Model}5:& P_{in,t+k|t}={\theta }_{1,k}lp\left(T_{s,t+k|t}-T_{r,t+k|t},a_1\right)+{\theta }_{2,k}lp\left(Q_{t+k|t},a_2\right)\hfill \end{array}$$

(20)

where for the k-step ahead, $Q_{t+k|t}$ is the heat load, ${\tilde{T}}_{H,t+k|t}$ is an estimate/forecast of the average house temperature, $a_i$ is the parameter for the first-order unit-gain low-pass function $lp(x,a_i)$, and ${\theta }_{i,k}$ is the model coefficients determined through the recursive fitting.

The estimate ${\tilde{T}}_{H,t+k|t}$ was based on the simple assumption that the supply temperature lost a fixed amount of degrees per house passed, making each house temperature equal to a scaled past supplied temperature:

$$\begin{array}{c}\hfill {\tilde{T}}_{H,t+k|t}=\frac{1}{N_{house}}\sum _{i=1}^{N_{house}}{0.9815}^i{T}_{s,t+k-i|t}\end{array}$$

(21)

Each model and its motivation can be summarized as: (1) Model 1 approximates the physical structure of the system. (2) Model 2 is simplified but preserves the multiplicative structure. (3) Model 3 adds more degrees of freedom to model 2. (4) Model 4 is simplified but approximates the fractional structures. (5) Model 5 is the simpler linear model.

From the above discussion, we can see that all of our models take the same general low-pass form, see (22). For the k-step, the low-pass form can be rewritten to a linear form with respect to the control variables on the form:

$$\begin{array}{cc}\hfill P_{in,t+k|t}& =\sum _{i=1}^n{\theta }_{i,k}lp\left(\frac{T_{i,t+k|t}-V_{i,t+k|t}}{d{V}_{i,t+k|t}},a_i\right)Q_{i,t+k|t}\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & =\sum _{i=1}^n\left(\frac{{\theta }_{i,k}}{{\theta }_{i,0}}\frac{a_i^k}{n}\frac{Q_{i,t+k|t}}{Q_{i,t-1}}{P}_{in,t-1}+{\theta }_{i,k}\sum _{j=0}^k{a}_i^{k-j}(1-a_i)\frac{Q_{i,t+j|t}}{d{V}_{i,t+j|t}}(T_{i,t+j|t}-V_{i,t+j|t})\right)\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill & =F_{1,k}{T}_{1,k}+{\Sigma }_k,T_{1,k}={[T_{1,t|t},T_{1,t+1|t},\cdots ,T_{1,t+k|t}]}^T\hfill \end{array}$$

(24)

where $T_{1,k}$ is the control variables, $T_{i,k}$, $V_{i,t}$, and $d{V}_{i,t}$ are the remaining variables inside the low-pass filter of (16)–(20), and $Q_{i,t}$ is a scaling variable on the low-pass filter. n is the number of terms in the model. Now, let $F_1$ and ${\Sigma }_0$ be the matrix and vector describing the entire power prediction:

$$\begin{array}{cc}\hfill P_{in,N_H-1}& =F_1{T}_{1,N_H-1}+{\Sigma }_0={[P_{in,t|t},P_{in,t+1|t},\cdots ,P_{in,t+N_H-1|t}]}^T\hfill \end{array}$$

(25)

giving the models a linear description over the prediction horizon with regard to the control variables.

3.2. Cost and Objectives

With the models defined, the MPC design can be addressed. The main part of the design was the definition of the optimum through the objectives in the cost function. The objectives were (1) inducing peak shaving through the reference-tracking of the power, (2) keeping the supply temperature low for the production cost, (3) keeping the production steady, through a minimum variation of the supply temperature. The objectives were all given as quadratic costs:

$$\begin{array}{cc}\hfill J_t& =\underset{T_s}{min}\sum _{i=0}^{N_H-1}\left(\left|\right|P_{in,t+k|t}-R_{p,k}{\left|\right|}_{H_y}^2+\left|\right|T_{s,t+k|t}{\left|\right|}_{H_u}^2+\left|\right|\Delta T_{s,t+k|t}{\left|\right|}_{H_{\Delta u}}^2\right)\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill & =\underset{T_s}{min}\left|\right|P_{in,N_H-1}-R_p{\left|\right|}_{H_y}^2+\left|\right|T_{s,N_H-1}{\left|\right|}_{H_u}^2+\left|\right|\Delta T_{s,N_H-1}{\left|\right|}_{H_{\Delta u}}^2\hfill \end{array}$$

(27)

where $R_p$ is the power reference, defined as the mean forecasted power within the horizon:

$$\begin{array}{c}\hfill R_p=mean({\Sigma }_0+F_1{T}_{s,forecast})\end{array}$$

(28)

The idea behind the first objective was to shift the top of the supply peaks to periods where the load from heat demand was below the reference, essentially averaging out the loads. Physically, this is done by storing up heat energy in the pipe system; the flow then distributes the heat along the pipe to be ready, increasing the temperature before higher demands occur. Applying the model structures, the cost function can be written as:

$$\begin{array}{cc}\hfill J_t& =\underset{T_s}{min}{T}_{s,N_H-1}^T{H}_0{T}_{s,N_H-1}+h_0{T}_{s,N_H-1}\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill H_0& =F_1^T{H}_Y{F}_1+H_U+I_d^T{H}_{\Delta U}{I}_d\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill h_0& ={({\Sigma }_0-R_p)}^T{H}_Y{F}_1+{\left(I_0{T}_{s,-1}\right)}^T{H}_Y{I}_d\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill I_0& =\left[\begin{array}{c}-1\\ 0\\ ⋮\\ 0\end{array}\right]I_d=\left[\begin{array}{cccc}1& 0& \cdots & 0\\ -1& 1& \ddots & ⋮\\ ⋮& \ddots & \ddots & 0\\ 0& \cdots & -1& 1\end{array}\right]\hfill \end{array}$$

(32)

The MPC based on each model was tuned visually to induce peak shaving. The tuned weights can be observed in

Table 1. The visual tuned weights of each control design.

Model	${\mathit{H}}_{\mathit{Y}}$	${\mathit{H}}_{\mathit{U}}$	${\mathit{H}}_{\mathsf{\Delta }\mathit{U}}$	Horizon ${\mathit{N}}_{\mathit{H}}$
1	1000	0.01	0	6
2	1000	0.01	0	6
3	1000	0.30	0.02	6
4	1000	0.25	0.02	6
5	1000	0.15	0.1	8

.

3.3. Constraints

The second part of MPC design is the constraints. For these controllers, the operational limits of the supply temperature regarding the temperature and its change were formulated as the constraints of the MPC:

$$\begin{array}{c}\hfill T_{s,forecast}\le T_{s,N_H-1}\le 85{}^{\circ }\mathrm{C}\end{array}$$

(33)

$$\begin{array}{c}\hfill -5{}^{\circ }\mathrm{C}\le \Delta T_{s,N_H-1}\le 5{}^{\circ }\mathrm{C}\end{array}$$

(34)

where the minimum supply temperature was chosen as the forecasted supply temperature, to ensure a guaranteed base performance of the system. The maximum supply temperature was chosen as 85 °C, while the change limits were chosen as a 5 °C limit.

The constraints can be expressed purely in terms of the supply temperatures (no change) as

$$\begin{array}{cc}\hfill G_u{T}_{s,N_H-1}& +G_{u_{-1}}{T}_{s,-1}\le g\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill G_u=\left[\begin{array}{c}{I}_{N_H}\\ -I_{N_H}\\ I_d\\ -I_d\end{array}\right]G_{u_{-1}}& =\left[\begin{array}{c}0\\ 0\\ I_0\\ -I_0\end{array}\right]g=\left[\begin{array}{c}85{}^{\circ }\mathrm{C}{\mathbf{1}}_{N_H}\\ -T_{s,forecast}\\ 5{}^{\circ }\mathrm{C}{\mathbf{1}}_{N_H}\\ 5{}^{\circ }\mathrm{C}{\mathbf{1}}_{N_H}\end{array}\right]\hfill \end{array}$$

(36)

where $T_{s,-1}$ is the past supply temperature.

4. Results

To evaluate the designed MPCs and the approach to peak shaving, we used simulations based on the model discussed in Section 2, using demand and data from a historical scenario. Each model case was given the exact same scenario and data: an initial 70-day data set for initializing the model spanning 1 March 2021 to 9 May 2021, then simulating the adaptive MPC over a period of 10 days, 10 May 2021 to 19 May 2021. For comparison, a base simulation was performed using the historical inputs and the system model. As the data were given in 1 h increments, the simulation and MPCs utilized a 1 h sampling time as well. The used data on the demands had a trend with the highest loads between 5 to 9 AM, with another small high around 18 to 19. The daily variation was around 0.4 MW, see

Table 2. Demand trends: the hourly trends of the heat demand, showing the mean, standard variation, minimum, and maximum values for each hour in MW.

Hour	0	1	2	3	4	5	6	7	8	9	10	11
mean	3.49	3.49	3.53	3.60	3.79	4.32	4.51	4.61	4.47	4.30	4.07	3.89
std. var.	0.30	0.29	0.28	0.28	0.31	0.39	0.33	0.27	0.22	0.26	0.38	0.38
min	3.20	3.12	3.24	3.28	3.32	3.71	3.99	4.14	4.18	3.97	3.64	3.29
max	4.08	3.99	4.01	4.11	4.24	5.03	4.99	5.00	4.69	4.69	4.71	4.54
Hour	12	13	14	15	16	17	18	19	20	21	22	23
mean	3.81	3.68	3.59	3.57	3.66	3.72	3.83	3.82	3.67	3.66	3.58	3.53
std. var.	0.42	0.39	0.43	0.39	0.39	0.37	0.33	0.31	0.33	0.28	0.33	0.31
min	3.08	3.05	2.94	2.99	3.20	3.30	3.48	3.50	3.24	3.40	3.25	3.23
max	4.56	4.36	4.26	4.18	4.29	4.42	4.44	4.43	4.27	4.17	4.13	4.13

.

4.1. Visual Evaluation

In Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10, the performance of the MPCs is shown alongside historical and/or base simulation performance, each figure consisting of four graphs: (1) The top graph illustrates the supply and return temperatures evolution over time, both simulated, forecasted, and historical supply temperature is displayed. (2) The inflow over time is shown in the top-middle graph for the historical inflow, and the simulations of the MPC and the baseline. (3) The bottom-middle graph displays the supply power of the simulations, as well as the utilized historical load of the system. (4) The bottom graph shows how the stored energy in the system changes over time, for the simulations with the MPC and the baseline.

Model 1’s performance is shown in Figure 6. We can observe the supply temperature only diverges from the forecast in order to charge the system with stored energy, temporally increasing supply power. These charging periods happen prior to the peak loads occurring, resulting in the following power peak being reduced/shaved in comparison to the base. Importantly, the new peaks from the charged power are significantly smaller than the base supply peaks. A similar reduction and behavior are seen in the flow.

In Figure 7, the performance of model 2 is shown. It performs similarly to model 1 with an increased supply temperature for storing heat energy in the pipe system, though the amount of shaving is larger for model 2.

Furthermore, model 3 also has a similar performance to model 1, as seen in Figure 8, the general difference being a higher size of the charged power peaks than in model 3, occasionally comparable to the base supply peaks.

In model 4’s case, shown in Figure 9, the charging of the system can be observed, though the peak shaving occurs only occasionally. The charging time with a higher supply temperature also tends to be longer.

The performance of model 5 is shown in Figure 10, showing peak shaving, with substantially longer charging periods, extending into the peak period.

Given the visual performances of the MPCs, model 2 would be the MPC design which performed best. The high-induced supply power peaks excluded model 3 and 4 as the new peaks ran counter to the purpose of peak shaving. Model 5 was likewise unsuitable as its long period of charging indicated higher production costs. Between models 1 and 2, model 2 was preferable with its slightly larger peak shaving.

4.2. KPI-Based Evaluation

Key performance indexes (KPI) can reduce performance to a few numbers. In this evaluation, we focused on the following list of KPIs, with each model’s KPI shown in

Table 3. KPI results for the simulation of the five models.

Model #	1	2	3	4	5	Ranking Best → Worst
Added temperature (°C)	0.853	0.991	0.814	1.088	1.538	3, 1, 2, 4, 5
	1.289%	1.497%	1.230 %	1.644%	2.322%
Added energy (MWh)	6.472	6.309	6.281	7.218	7.961	3, 2, 1, 4, 5
	0.609 %	0.591%	0.594 %	0.679%	0.749%
Shaved peak power 1-sided (MW)	−0.3411	−0.4303	−0.3347	−0.2460	−0.3036	2, 1, 3, 5, 4
Shaved peak power 2-sided (MW)	−0.1299	−0.1936	−0.1277	−0.0595	−0.1395	2, 1, 3, 5, 4

.

• Added temperature: the average amount of degrees added to the supply temperature, in comparison to the forecast.
• Added energy: the additional energy applied to the entrance of the system
• Shaved peak power: the average difference of max peak powers at each event between base and MPC. Each event lies within a 12 h interval (12 AM–12 PM or 12 PM–12 AM).
  • One-sided: only model peaks in the interval occurring prior to the maximum base peak are considered
  • Two-sided: all model peaks in the interval are considered.

Since the performance of models 1–3 were relatively similar, the best model then depended on the importance given to each KPI. The increase in energy consumption was fairly even around 0.6%. Models 1 and 3 performed better keeping the temperature low, while model 2 was better at peak shaving. Overall, model 2 was the best performer when prioritizing peak shaving, with its 0.43 MW average reduction, roughly $8\%$ of the general power peak or $4\%$ for the two-sided peak shaving.

The overall KPI results show the best results were from the models preserving the multiplicative structure, models 1, 2, and 3, outperforming models 4 and 5 in every KPI.

Based on both the KPI and the visual analysis, model 2 was the best choice for a design model, with model 1 being a close alternative, providing a trade-off between peak shaving and low-temperature importance.

5. Conclusions

In this paper, we proposed formulations of adaptive model predictive control (MPC) for handling the peak load problem of supply power in district heating, through pipe energy storage. A model was similarly presented, for the purpose of simulating district heating systems approximated as a single pipe with distributed loads.

Several MPC formulation designs were proposed and applied to a case study of the district heating system in Brønderslev, Denmark. It was found that a simple physics-inspired model, named model 2, had the best performance among the designs, providing around $8\%$ of peak shaving in supply power without inducing unwarranted behaviors, such as new large power peaks.

Future Work

A continuation of the work could include a real-world application to an actual district heating system (DHS) or application to data from other DHS than the Brønderslev case. As model 1 was close to model 2 in performance, an updated design description of the house temperatures or the inclusion of internal pipe temperature measurements in the design could be added.