Model Predictive Control of District Heating Grids Using Stabilizing Terminal Ingredients


The transformation of fossil fuel-based district heating grids (DHGs) to CO\(_2\)-neutral DHGs requires the development of novel operating strategies. Model predictive control (MPC) is a promising approach, as knowledge about future heat demand and heat supply can be incorporated into the control, operating constraints can be ensured and the stability of the closed-loop system can be guaranteed. In this paper, we employ MPC for DHGs to control the system mass flows and injected heat flows. Following common practice, we derive terminal ingredients to stabilize given steady state temperatures and storage masses in the DHG. To apply MPC with terminal ingredients, it is crucial that the system under control is stabilizable. By exploiting the particular system structure, we give a sufficient condition for the stabilizability in terms of the grid topology and hence, for the applicability of the MPC scheme to DHGs. Furthermore, we demonstrate the practicability of the application of MPC to an exemplary DHG in a numerical case study.


1 Introduction↩︎

District heating grids (DHGs) play an important role in the decarbonization of the heat sector [1]. A DHG is a network of pipes, which are used to transport water as heat energy carrier in between producers, where heat flows are injected, and consumers, where heat flows are extracted with individual temperature requirements. Major advantages of DHGs are the possibility of a decentralized integration of RES-based heat producers, low temperature waste heat sources and thermal energy storages (TESs), which are all needed for the transformation towards RES-based DHGs [1], [2]. This transformation influences the operation of future DHGs, such that conventional operating strategies are not feasible for an RES-based operation of DHGs [3]. This motivates the development of novel operating strategies for RES-based DHGs.

To give a brief overview of the state of the art in the field of operational management of DHGs, we distinguish between energy management strategies (EMSs), which compute control signals by solving an optimization problem iteratively, and control approaches, which ensure the stability of the closed-loop system by a suitable feedback controller.

An example of economically motivated EMS for DHGs can be found e.g. in [4], where the flexibility stemming from sector coupling is used to increase operational efficiency. In [5], suitable discretization schemes to obtain a nonlinear optimization model for DHGs with increased modeling accuracy are investigated. A reduced order model of a DHG is investigated to reduce computational complexity in [6].

Control strategies for the nonlinear hydraulics of DHGs are studied e.g. in [7], where a passivity-based approach is proposed. A temperature controller based on Lyapunov-Krasovskii stability theory is proposed and experimentally validated in [8]. In [9], the mass stored in a TES as well as the system temperatures are controlled in a DHG with a single producer and multiple consumers. Furthermore, passivity-based approaches are proposed in [10] and [11] to control the thermo-hydraulics of DHGs and to control the frequency and temperature in a coupled power and heating grid, respectively.

To combine the advantages from both optimization-based EMS, e.g. a predictive operating behavior as well as constraint satisfaction, and control strategies guaranteeing closed-loop stability, model predictive control (MPC) is a highly promising approach. In MPC, an optimization problem is solved iteratively to obtain a stabilizing control input for each control action [12]. However, none of the above papers addresses an optimization-based controller with assured stability for DHGs.

In [13], a set constraint, so called terminal constraint, and a cost term, so called terminal cost, was added to the nominal optimization problem of a MPC to obtain an asymptotically stable closed-loop behaviour. The terminal constraint is defined via a neighborhood around the desired steady-state set point, called terminal region, and enforces the last entry of each optimal state trajectory resulting from solving the optimization problem to lie in this region. Additionally, the terminal cost penalizes the distance of the last entry of each optimal state trajectory to the desired steady-state set point. Together, terminal region and terminal cost form the terminal ingredients. The approach from [13] is extended to the discrete-time case in [14]. Additionally, [13], [14] provide a framework to derive a terminal cost and a terminal region. However, to apply MPC with terminal ingredients, it is crucial that the linearized system model is stabilizable with respect to the desired steady state. Furthermore, the algorithm to calculate terminal cost and terminal constraint can result in conservative MPC laws [13]. This motivates both the examination of the local stabilizability of a DHG model as well as a numerical investigation of the practicability of MPC with terminal ingredients conducted in the present work.

Driven by the importance of RES-based DHG operation in the transition to a zero-emission heating sector and recognizing the pivotal role of MPC in this context, this paper presents three primary contributions: At first, building upon [5], [10], we develop, for a DHG with multiple producers, multiple distributed storage units and consumers, an ordinary differential equation-based (ODE-based) state model that describes relevant mass and temperature dynamics and which can be used to implement an MPC. Compared to [5], [10], we introduce heat losses from [5] into the modeling approach in [10] and allow for a less restrictive placement of TES which were not considered in [5]. Additionally, multiple temperature layers are considere in this work. Second, we derive a sufficient condition for the local stabilizability of the developed DHG model by exploiting its structural properties. Third, we calculate suitable terminal ingredients for an exemplary DHG and demonstrate the practicability of MPC with terminal ingredients to control DHGs in a case study.

The remainder of the paper is structured as follows. In Section 2, we describe the DHG model. Then, in Section 3, we formulate the optimization problem of the MPC, show the stabilizability property of the system and derive a stabilizing control law. A case study demonstrating the practicability of the approach is presented in Section 4. The presented work is summarized and conclusions are drawn in Section 5.


The set of (positive and negative) real numbers and non-negative integers is denoted by (\(\mathbb{R}_+\) and \(\mathbb{R}_-\)) \(\mathbb{R}\) and \(\mathbb{N}\), respectively. For a finite set \(\mathcal{S}\), \(|\mathcal{S}|\) denotes its cardinality. For a vector \(x \in \mathbb{R}^n\), \(\mathrm{diag}\big({x}\big) \in \mathbb{R}^{n \times n}\) denotes a diagonal matrix having the elements of \(x\) on the main diagonal. The weighted norm \(x^\top Q x\), where \(Q \in \mathbb{R}^{n \times n}\), is denoted by \(\|x\|_Q^2\). Let \(X \in \mathbb{R}^{n \times m}\) denote a \(n \times m\) matrix and let \(\mathcal{P}\) denote a set with \(|\mathcal{P}|=n\), in which each element of \(\mathcal{P}=\{p_1,\ldots,p_n\}\) is associated with one row of \(x\) or \(X\) and \(\mathcal{P}_{\mathrm{sub}} \subset \mathcal{P}\). Then, \((x)_{p \in \mathcal{P}_{\mathrm{sub}}} \in \mathbb{R}^{|\mathcal{P}_{\mathrm{sub}}|}\) and \((X)_{p \in \mathcal{P}_{\mathrm{sub}}} \in \mathbb{R}^{|\mathcal{P}_{\mathrm{sub}}| \times m}\) denote a vector and a matrix comprising only those rows of \(x\) and \(X\), respectively, for which \(p \in \mathcal{P}_{\mathrm{sub}}\), and \((x)_{p_i} \in \mathbb{R}\) denotes the entry of \(x\) that is associated with \(p_i \in \mathcal{P}\). The \(n \times n\) identity matrix is denoted by \(I_n\). Further, \(|X| \in \mathbb{R}^{n \times m}\) is the matrix resulting from \(X\) by taking absolute values of its entries.

2 Thermo-hydraulic model↩︎

The goal of this section is to develop an ODE-based model of a general DHG that can be used to design an MPC feedback law that stabilizes the water mass stored at TESs and the temperatures within the DHG. Then, the model needs to describe the hydrodynamics of water mass stored at TESs, the thermodynamics of water in the DHG components as well as Kirchhoff’s laws for hydraulic networks. To do so, we write-down the balance of mass and the balance of energy of a general DHG that consists of a network of pipes connecting stratified TESs and heat exchangers (HXs). The HXs are associated with either an injected heat flow from a producer or an extracted heat flow from a consumer. In Figure 1, a schematic representation of an exemplary DHG containing two HXs associated with the injected heat flows \(P_{\mathrm{pr},1} \in \mathbb{R}_+\) and \(P_{\mathrm{pr},2} \in \mathbb{R}_+\), three HXs associated with the extracted heat flows \(P_{\mathrm{d},1} \in \mathbb{R}_+\), \(P_{\mathrm{d},2} \in \mathbb{R}_+\) and \(P_{\mathrm{d},3} \in \mathbb{R}_+\) and two TESs interconnected by a network of pipes is shown.

Figure 1: Schematic of an exemplary DHG.

We make the following standard assumptions to derive the DHG model [5], [15], [16].

Assumption 1. (i) A fixed control volume is associated with each component. (ii) All components are completely filled with water at all times. (iii) The water in the components is incompressible. (iv) Gravitational acceleration is neglected. (v) The mass flow through a pipe is one-dimensional. (vi) The internal energy of water is regarded as the main source of energy, i.e., other forms of energy are negligible. (vii) The internal energy of water depends linearly on the temperature of the water.

For modeling TESs, we make the following assumptions [10].

Assumption 2. (i) Forces acting within a TES or between individual water layers of a TES are neglected. (ii) The control volume of a TES can ideally be divided into a hot layer and a cold layer.

Formally, we represent the topology of a DHG by a connected graph \(\mathcal{G}=(\mathcal{V},\mathcal{E})\), where \(\mathcal{V}=\{v_1, \ldots, v_{|\mathcal{V}|}\}\) and \(\mathcal{E}=\{e_1, \ldots, e_{|\mathcal{E}|}\}\) denote the sets of vertices and edges, respectively. The set \(\mathcal{V}\) consists of three subsets, namely the set of vertices associated with pipe intersections \(\mathcal{V}_{\mathrm{p}} \subset \mathcal{V}\), the set of vertices associated with hot layers of a TES \(\mathcal{V}_{\mathrm{h}} \subset \mathcal{V}\), and the set of vertices associated with cold layers of a TES \(\mathcal{V}_{\mathrm{c}} \subset \mathcal{V}\), such that \(\mathcal{V} = \mathcal{V}_{\mathrm{p}} \cup \mathcal{V}_{\mathrm{h}} \cup \mathcal{V}_{\mathrm{c}}\). An edge \(e \in \mathcal{E}\) is associated with a pipe \(\mathcal{E}_{\mathrm{p}} \subset \mathcal{E}\), or a HX. Edges associated with HXs, where heat flow is injected or heat flow is extracted, are collected in the sets \(\mathcal{E}_{\mathrm{pr}} \subset \mathcal{E}\) or \(\mathcal{E}_{\mathrm{d}} \subset \mathcal{E}\), respectively, such that \(\mathcal{E} = \mathcal{E}_{\mathrm{p}} \cup \mathcal{E}_{\mathrm{pr}} \cup \mathcal{E}_{\mathrm{d}}\). For the following modeling, a directed graph is required, which can be obtained by assigning an arbitrary orientation to \(e \in \mathcal{E}\). For an intuitive physical interpretation, we choose the orientation of \(e \in \mathcal{E}\) corresponding to the direction of flow through a pipe associated with this edge. Then, the edge \(e = (v,w) \in \mathcal{E}\) exists, if water flows from \(v \in \mathcal{V}\) to \(w \in \mathcal{V}\) for \(v \neq w\). In Figure 2, the graph representation of the DHG from Figure 1 is shown. The vertex-edge incidence matrix \(B \in \mathbb{R}^{|\mathcal{V}| \times |\mathcal{E}|}\) of \(\mathcal{G}\) is defined element-wise via \[(B)_{i,j}= \begin{cases} 1, & \text{if~e_j=(w,v_i)\in\mathcal{E}}, \\ -1, & \text{if e_j=(v_i,w)\in\mathcal{E}}, \\ 0, & \text{otherwise.} \end{cases}\]

Figure 2: Graph representation of the exemplary DHG from Figure 1 with definitions of associated sets \(\mathcal{E}_{\mathrm{d}}\), \(\mathcal{E}_{\mathrm{pr}}\), \(\mathcal{V}_{\mathrm{h}}\) and \(\mathcal{V}_{\mathrm{c}}\). Dotted rectangles frame vertices of a corresponding TES.

We proceed by introducing some further concepts from graph theory from [17] that are used in the following sections. A path defines a sequence of edges that link a sequence of distinct vertices. Let \(G=(V,E)\) with set of vertices \(V\) and set of edges \(E=\{ \mathrm{e}_1, \ldots, \mathrm{e}_{|E|} \}\) denote a weakly connected graph, i.e., a graph that has an undirected path between every two vertices. Moreover, let a cycle be an undirected sequence of consecutive vertices and edges in which only the first and last vertex are equal. If \(G\) contains cycles, removing all edges \(b \in \mathcal{F} \subset E\) from \(G\) yields a connected subgraph \(S=(V,E \setminus \mathcal{F})\) that does not contain cycles. Then \(S\) is called the spanning tree of \(G\), whereas \(\mathcal{F}\) is called the set of chords. For the example given in Figure 2, \(\{ e_3, e_6, e_7 \}\) would be a valid choice of chords. A weakly connected graph has at least one spanning tree [17]. In what follows, every set of edges \(L_j \subset E\), \(j=1,\ldots,|\mathcal{F}|\), that forms a unique cycle by adding \(b \in \mathcal{F}\) to \(S\) is called a fundamental cycle. The orientation of a fundamental cycle is defined in accordance with the orientation of \(b \in \mathcal{F}\) that forms it. Now, we can define the fundamental cycle matrix \(F \in \mathbb{R}^{|\mathcal{F}| \times |E|}\) of \(G\) element-wise as \[(F)_{i,j}= \begin{cases} 1, & \text{if } \mathrm{e}_j \in L_i \text{ is oriented as } L_i, \\ -1, & \text{if } \mathrm{e}_j \in L_i \text{ is not oriented as } L_i, \\ 0, & \text{if } \mathrm{e}_j \notin L_i \text{, otherwise.} \end{cases}\]

Balance of mass↩︎

We assign a mass to each component of the DHG and a mass flow to each pipe. Thus, vertex \(v \in \mathcal{V}\) has a mass and each edge \(e \in \mathcal{E}\) has a mass and a mass flow. The sign of a mass flow is considered to be positive in the direction of the associated edge. From Assumption 1 it follows that the masse of edges \(e \in \mathcal{E}\) and vertice \({v \in \mathcal{V}_{\mathrm{p}}}\) are constant and collected in the vectors \({(m_{\mathrm{v}})_{v \in \mathcal{V}_{\mathrm{p}}} \in \mathbb{R}_+^{|\mathcal{V}_{\mathrm{p}}|}}\) and \({m_{\mathrm{e}} \in \mathbb{R}_+^{|\mathcal{E}|}}\), respectively. Furthermore, from Assumption 1 it also follows that the total mass of the \(i^{\text{th}}\) TES denoted by \(m_{\mathrm{tes},i} \in \mathbb{R}_+\) is constant, while we allow for a varying ratio of hot to cold water in each TES. Thus, if the vertices \(v_h \in \mathcal{V}_{\mathrm{h}}\) and \(v_c \in \mathcal{V}_{\mathrm{c}}\) are associated with the \(i^{\text{th}}\) TES, the balance of mass reads \[\label{eq:mh43mc61mtes} m_{\mathrm{tes,}i} = (m_{\mathrm{v}})_{v_h}(t) + (m_{\mathrm{v}})_{v_c}(t),\quad \forall t \geq 0.\tag{1}\]

We define \(m_{\mathrm{h}} := (m_{\mathrm{v}})_{v \in \mathcal{V}_{\mathrm{h}}}\) and \(m_{\mathrm{c}} := (m_{\mathrm{v}})_{v \in \mathcal{V}_{\mathrm{c}}}\) to obtain the following system of differential-algebraic equations \[\label{eq:DAE} \begin{bmatrix} 0 \\ \dot{m}_{\mathrm{h}}(t) \\ \dot{m}_{\mathrm{c}}(t) \end{bmatrix} = \begin{bmatrix} (B)_{v \in \mathcal{V}_{\mathrm{p}}} \\ (B)_{v \in \mathcal{V}_{\mathrm{h}}} \\ (B)_{v \in \mathcal{V}_{\mathrm{c}}} \end{bmatrix} q_{\mathrm{e}}(t) = \begin{bmatrix} (B)_{v \in \mathcal{V}_{\mathrm{p}}} \\ (B)_{v \in \mathcal{V}_{\mathrm{h}}} \\ (-B)_{v \in \mathcal{V}_{\mathrm{h}}} \end{bmatrix} q_{\mathrm{e}}(t),\tag{2}\] which collect the mass balance at each vertex of the DHG, where \(q_{\mathrm{e}} \in \mathbb{R}_+^{|\mathcal{E}|}\) denotes a vector collecting the mass flows through all edges. Note that \(\dot{m}_{\mathrm{h}} = -\dot{m}_{\mathrm{c}}\) follows from 2 . Thus, the dynamics of the mass stored in the TESs can completely be determined by the dynamics of \(\dot{m}_{\mathrm{h}}\). Also note, since the mass flows over the edges \(q_{\mathrm{e}}\) linearly depend on each other due to Kirchhoff’s laws, see first line of 2 , they can be expressed by a subset of independent variables [18], which we refer to as the fundamental mass flows describing the entire hydraulic state of the DHG.

For the considered DHG, we proceed to identify the fundamental mass flows following the procedure described in [10]. In hydraulic networks whose topology can be represented by a graph, where a constant mass is associated with all vertices, the fundamental mass flows are those over the chords of its graph [18]. In [10], it is shown for a hydraulic network with variable masses associated with certain vertices of its graph representation, as considered in this work, that the fundamental mass flows are those over the chords of a reduced graph of \(\mathcal{G}\), i.e., a graph that contains the same amount of edges but less vertices than \(\mathcal{G}\). In what follows, let \(G\) denote this reduced graph. To obtain \(G\), we merge the vertices from \(\mathcal{V}\) associated with the same TESs. This yields the reduced graph of \(\mathcal{G}\) representing the same hydraulic network, but assigning a constant mass to all vertices of \(G\). The mass flows over \(\mathcal{F}\) of \(G\) are also the fundamental mass flows of the hydraulic network represented by \(\mathcal{G}\) [10].

To describe the relations between the fundamental mass flows and \(q_{\mathrm{e}}\), let \(q_{\mathrm{c}} \in \mathbb{R}_+^{|\mathcal{F}|}\) denote a vector that collects the mass flows over the chords \(\mathcal{F}\). Then, it is possible to write \(q_{\mathrm{e}} = F^\top q_{\mathrm{c}}\) [10] and obtain the ODE-based balance of mass \[\label{eq:reduced32mass32balance} \dot{m}_{\mathrm{h}}(t) = (B)_{v \in \mathcal{V}_{\mathrm{h}}} F^\top q_{\mathrm{c}}(t).\tag{3}\]

Balance of energy↩︎

From Assumption 1, it follows that the thermodynamic state of each component of the DHG is characterized by its temperature. We assign an average temperature to the respective control volume of each pipe, HX, pipe intersection and TES layer. Thus, any vertex \(v \in \mathcal{V}\) or any edge \(e \in \mathcal{E}\) also has an average temperature. Let \({T_{\mathrm{v}}\in\mathbb{R}^{|\mathcal{V}|}}\) and \({T_{\mathrm{e}}\in\mathbb{R}^{|\mathcal{E}|}}\) denote vectors collecting the vertices’ temperatures and the edges’ temperatures, respectively. Following [5], [10], we assume the incoming temperature of an edge to be equal to the temperature of its source vertex as well as the outgoing temperature of an edge to be equal to its average temperature.

Additionally, we define the mass matrices \(M_{\mathrm{v}} = \mathrm{diag}\big({m_{\mathrm{v}}}\big)\) and \(M_{\mathrm{e}} = \mathrm{diag}\big({m_{\mathrm{e}}}\big)\) collecting the mass of vertices and edges on their diagonals, respectively, and introduce the abbreviation \(B_+ = \frac{1}{2}(|B|+B)\) and \(B_- = \frac{1}{2}(|B|-B)\). Then, it is possible to write the balance of energy at vertices and edges in vector form as follows [10]: \[\label{eq:energy32balance32reduced32form} \begin{align} M_{\mathrm{v}}(t) \dot{T}_{\mathrm{v}}(t) = &-\mathrm{diag}\big({B_+ F^\top q_{\mathrm{c}}(t) + \kappa_{\mathrm{v}}}\big) T_{\mathrm{v}}(t) \\ &+ B_+ \mathrm{diag}\big({F^\top q_{\mathrm{c}}(t)}\big) T_{\mathrm{e}}(t) + \kappa_{\mathrm{v}} T_{\mathrm{a}}, \\ M_{\mathrm{e}} \dot{T}_{\mathrm{e}}(t) =\, &\mathrm{diag}\big({F^\top q_{\mathrm{c}}(t)}\big) B_-^\top T_{\mathrm{v}}(t) \\ &- \mathrm{diag}\big({F^\top q_{\mathrm{c}}(t) + \kappa_{\mathrm{e}}}\big) T_{\mathrm{e}}(t) + \kappa_{\mathrm{e}} T_{\mathrm{a}} \\ &+ P_{\mathrm{pr}}(t) - P_{\mathrm{d}}(t), \end{align}\tag{4}\] where \(\kappa_{\mathrm{v}} \in \mathbb{R}_+^{|\mathcal{V}|}\) and \(\kappa_{\mathrm{e}} \in \mathbb{R}_+^{|\mathcal{E}|}\) denote vectors collecting heat loss coefficients of vertices and edges, respectively, and \({T_{\mathrm{a}} \in \mathbb{R}}\), \({P_{\mathrm{pr}} \in \mathbb{R}_+^{|\mathcal{E}|}}\) and \({P_{\mathrm{d}} \in \mathbb{R}_+^{|\mathcal{E}|}}\) denote ambient temperature, injected and extracted heat flow at the edges, respectively.

State space model↩︎

To conclude the thermo-hydraulic modeling, let us summarize equations 3 and 4 into an ODE-based state space model of the form \[\label{eq:overall32state32model} \begin{align} \dot{x}(t) &= f(x,u,d) = M(t)^{-1} \big( A(t) x(t) + E_{\mathrm{u}} u(t) + E_{\mathrm{d}} d(t) \big), \end{align}\tag{5}\] using the state vector \(x=\begin{bmatrix} m_{\mathrm{h}}^\top & T_{\mathrm{v}}^\top & T_{\mathrm{e}}^\top\end{bmatrix}^\top \in \mathbb{R}^n\), the control input vector \(u=\begin{bmatrix} q_{\mathrm{c}}^\top & P_{\mathrm{pr}}^\top\end{bmatrix}^\top \in \mathbb{R}^m\), and the disturbance vector \(d=\begin{bmatrix} P_{\mathrm{d}}^\top & T_{\mathrm{a}} \end{bmatrix}^\top \in \mathbb{R}^p\), where \(n=|\mathcal{V}_{\mathrm{h}}|+|\mathcal{V}|+|\mathcal{E}|\), \(m=|\mathcal{F}|+|\mathcal{E}_{\mathrm{pr}}|\) and \(p=|\mathcal{E}_{\mathrm{d}}|+1\). Then, we define the matrices \[\begin{align} \label{eq:def95M95A} M(t)^{-1}&= \begin{bmatrix} I_{|\mathcal{V}_{\mathrm{h}}|} & 0 & 0 \\ 0 & M_{\mathrm{v}}(t) & 0 \\ 0 & 0 & M_{\mathrm{e}} \end{bmatrix}^{-1},\; \end{align}\tag{6}\] \[\begin{align} A(t) &= \begin{bmatrix} 0 & 0 \\ 0 & \tilde{A}(q_{\mathrm{c}}(t))-\tilde{A}_0 \end{bmatrix},~ \tilde{A}_0={\rm diag}\left(\begin{bmatrix} \kappa_{\mathrm{v}} \\ \kappa_{\mathrm{e}} \end{bmatrix}\right) \end{align}\] \[\label{eq:def95tilde95A} \tilde{A}(q_{\mathrm{c}})=\begin{bmatrix} -\mathrm{diag}\big({B_+ F^\top q_{\mathrm{c}} }\big) & B_+ \mathrm{diag}\big({F^\top q_{\mathrm{c}}}\big) \\ \mathrm{diag}\big({F^\top q_{\mathrm{c}}}\big) B_-^\top & -\mathrm{diag}\big({F^\top q_{\mathrm{c}} }\big) \end{bmatrix},\tag{7}\] \[\begin{align} E_{\mathrm{u}} = \begin{bmatrix} (B)_{v \in \mathcal{V}_{\mathrm{h}}} F^\top & 0 \\ 0 & 0 \\ 0 & I_{|\mathcal{E}|} \end{bmatrix},&~ E_{\mathrm{d}} = \begin{bmatrix} 0 & 0 \\ 0 & \kappa_{\mathrm{v}} \\ -I_{|\mathcal{E}|} & \kappa_{\mathrm{e}} \end{bmatrix}, \end{align}\] to write the dynamics 3 and 4 in the form of 5 .

For the design of the MPC, we consider an explicit Euler discretization of the system dynamics 5 with step size \({\Delta t\in\mathbb{R}_+}\) given by \[\label{eq:euler} \begin{align} x(k+1) &= f_\delta(x(k), u(k), d(k)) \\ &:=x(k) + \Delta t \, f(x(k), u(k), d(k)), \end{align}\tag{8}\] where the discrete-time input \(u(k)\) is given by the continuous time input at time \(k\,\Delta t\) for \(k \geq 0\).

3 Model predictive controller↩︎

We follow [13], [14] and use a quadratic stage cost function, a terminal cost function and a terminal region defined as \[\label{eq:terminal95ingredients} \begin{align} \ell(x,u) &= \|x-\overline{x}\|_Q^2 + \|u-\overline{u}\|_R^2,~ C_P(x) = \|x - \overline{x}\|_P^2, \\ \mathbb{X}_{C_P,\alpha} &= \{x \in \mathbb{R}^n ~:~ C_P(x) \leq \alpha\}, \text{ for some }\alpha > 0, \end{align}\tag{9}\] respectively, where \(Q \in \mathbb{R}^{n \times n}\), \(R \in \mathbb{R}^{m \times m}\), \(P \in \mathbb{R}^{n \times n}\) denote positive definite weight matrices and \({(\overline{x}, \overline{u}, \overline{d})}\) denote a steady state tuple of 8 . Moreover, we collect lower and upper bounds for the states as well as for the control inputs in the vectors \(x_{\mathrm{lb}} \in \mathbb{R}^n\) and \(x_{\mathrm{ub}} \in \mathbb{R}^n\) as well as \(u_{\mathrm{lb}} \in \mathbb{R}^m\) and \(u_{\mathrm{ub}} \in \mathbb{R}^m\), respectively. The vector \(q_{\mathrm{ub}} \in \mathbb{R}_+^{|\mathcal{E}|}\) collects the upper limits for the mass flow over each pipe. Then, we can define a set describing the state constrains and a set describing the control input constraints as follows \[\begin{align} \mathbb{X} &= \{ x \in \mathbb{R}^n~:~ x_{\mathrm{lb}} \leq x \leq x_{\mathrm{ub}} \}, \\ \mathbb{U} &= \{ u \in \mathbb{R}^m~:~ F^\top q_{\mathrm{c}} \leq q_{\mathrm{ub}}, u_{\mathrm{lb}} \leq u \leq u_{\mathrm{ub}} \}, \end{align}\] respectively. Now, we can formulate the corresponding MPC problem with terminal ingredients for DHGs.

Problem 1 (MPC with terminal ingredients for DHG). Consider the discrete-time dynamics 8 . At time \(k_0\), given the initial state \(x_{k_0} \in \mathbb{X}\), and a desired steady state tuple \((\overline{x},\overline{u},\overline{d})\), find a state and control input trajectory \((x_{k_0}^*(\cdot),u_{k_0}^*(\cdot))\) that solves \[\begin{align} \underset{x(\cdot),u(\cdot)}{\min} ~ & \sum_{k=0}^{N-1} \ell(x(k),u(k)) + C_P(x(N)) \\ \text{s.t. } & x(k+1) = f_\delta(x(k),u(k),\overline{d}), \quad k=0,\ldots,N-1,\\ & (x(k), u(k)) \in \mathbb{X} \times \mathbb{U}, \quad k=0,\ldots,N-1,\\ & x(N) \in \mathbb{X}_{C_P,\alpha}, \\ & x(0) = x_{k_0}. \end{align}\]

Note that knowledge of future disturbance values, i.e., perfect forecast, needs to be assumed in order to solve Problem 1. By iteratively solving Problem 1, we obtain an optimal input sequence \(u_{k_0}^*=(u_{k_0}^*(0),u_{k_0}^*(1),\ldots,u_{k_0}^*(N-1))\in \mathbb{U}^{N}\). The MPC closed-loop control law \(\mu:\mathbb{X} \rightarrow \mathbb{U}\) is then obtained by using only the first element of the input sequence \(u_{k}^*(0)\) for the given initial state \(x_{k}\), i.e., \(\mu(x_{k}):=u_{k}^*(0)\). Applying \(\mu(x_{k_0})\) will result in a system state \(x(k_0+1)=x_{k_0+1}\), the time-horizon is moved one step ahead and Problem 1 is solved again with initial value \(x_{k_0+1}\) instead of \(x_{k_0}\).

Asymptotic stabilization of set points↩︎

In the following, we show that the states of the system 8 will asymptotically converge to a given desired steady state tuple \((\overline{x},\overline{u},\overline{d})\), when the control law \(\mu(\cdot)\) is applied. Recall that for a discrete time system \(x(k+1)=f_\delta(x(k))\), \(x(0)=x_0\) an equilibrium \(\overline{x}\in\mathbb{X}\) is called asymptotically stable if for each \(\varepsilon>0\) there exists \(\eta(\varepsilon)>0\) such that \(\|x_{0}-\overline{x}\|<\eta(\varepsilon)\) implies \(\|x(k)-\overline{x}\|<\varepsilon\) for all \(k\geq 0\) and that \(x(k)\rightarrow\overline{x}\) holds as \(k\rightarrow\infty\). Furthermore, we find it useful to recall the meaning of stabilizability.

Definition 2 (Stabilizability [19]). The linear discrete-time system \(x(k+1)=A_\delta x(k) + B_\delta u(k)\) is called stabilizable if there exists a matrix \(G\), such that \(A_\delta + B_\delta G\) is stable, i.e., all its eigenvalues lie in the unit disc.

It was shown in [14] that the closed-loop system’s equilibrium is asymptotically stable, see also [13], if the following conditions hold.

Condition 3. (i) \(f_\delta\) is continuously differentiable. (ii) \(f_\delta(\overline{x},\overline{u},\overline{d})=\overline{x}\). (iii) The discrete time system has a unique solution for every input sequence and any initial value. (iv) \(\mathbb{U}\) is compact and \(\overline{u}\) lies in the interior of \(\mathbb{U}\). (v) The states are perfectly measured. (vi) The MPC algorithm is initially feasible, i.e., Problem 1 has a solution for the initial state \(x_{k_0}\), (vii) the linearized system of 8 , i.e., \[x(k+1)= \overline{x} + \left.\frac{\partial f_\delta}{\partial x}\right\rvert_{(\overline{x},\overline{u},\overline{d})} (x(k)-\overline{x}) + \left.\frac{\partial f_\delta}{\partial u}\right\rvert_{(\overline{x},\overline{u},\overline{d})}(u(k)-\overline{u}),\] is stabilizable in the sense of Definition 2.

The foregoing Conditions 3.(i) to 3.(v) are trivially fulfilled. Condition 3.(vi) is usually needed in the literature, see e.g. [12], and cannot be proven in general.

We will show that Condition 3.(vii) is satisfied for any DHG that can be described by the system 8 . For this, we need the following technical assumption, which can be shown to hold due to the hydraulic decoupling introduced by TESs in typical DHG configurations and in particular holds for the exemplary DHG from Figure 1 and 2.

Assumption 3. For 8 , \(\ker F(B)_{v\in\mathcal{V}_{\mathrm{h}}}^\top=\{0\}\) holds.

Furthermore, recall [10].

Lemma 1. For the steady state mass flows \(\overline{q}_{\mathrm{c}}\), the matrix \(\tilde{A}\) given by 7 fulfills \(\tilde{A}(\overline{q}_{\mathrm{c}})+\tilde{A}(\overline{q}_{\mathrm{c}})^\top\leq 0\).

Proposition 4. Consider the discretized system 8 with step size \(\Delta t>0\) and Assumption 3. There exists a sufficiently small \(\Delta t\), such that Condition 3 holds.

Proof. We prove that Condition 3.(vii) is satisfied, by constructing a matrix \(G\), such that Definition 2 holds for the system 8 . We define the shifted variables \(\tilde{x}(k):=x(k)-\overline{x}\) and \(\tilde{u}(k):=u(k)-\overline{u}\). Then, the linearized discrete time system from Condition 3.(vii) can be rewritten using 8 as \[\begin{align} x(k+1) \nonumber &=\overline{x} + \frac{\partial f_\delta}{\partial x}\tilde{x}(k)+ \frac{\partial f_\delta}{\partial u}\tilde{u}(k)\nonumber \\ &=x(k)+\Delta t \left(\frac{\partial f}{\partial x}\tilde{x}(k)+ \frac{\partial f}{\partial u}\tilde{u}(k)\right),\label{eq:discrete95with95cnts} \end{align}\tag{10}\] where the derivatives are evaluated at \((\overline{x},\overline{u},\overline{d})\). To compute the derivatives, we define \({A_0 = -\mathrm{diag}\big({\begin{bmatrix} 0 & \kappa_v^\top & \kappa_e^\top \end{bmatrix}^\top}\big)}\), \({e_i\in\mathbb{R}^{|\mathcal{F}|}}\) as the \({i^{\mathrm{th}}}\) canonical unit vector, i.e., \({(e_i)_j=1}\) if \({i=j}\), and zero entries otherwise, and \({A_i=\begin{bmatrix} 0\\ \hat{A}_i \end{bmatrix}}\), where \({\hat{A}_i=\tilde{A}(e_i)}\), is given by 7 . Then, we can write \[f(x,u,\overline{d})=M^{-1}(A_0+\sum_{i=1}^{|\mathcal{F}|} u_i A_i)x + E_u u + E_d \overline{d}.\] Using the steady state condition \(f(\overline{x},\overline{u},\overline{d})=0\), abbreviating \(\hat{A}(\overline{x}):=\begin{bmatrix}\hat{A}_1\overline{x} \ldots \hat{A}_{|\mathcal{F}|}\overline{x}\end{bmatrix}\) and by considering the positive definite matrix \(\overline{M}\), which is obtained from \(M(t)\) given in 6 by using the steady state values for the vertex masses, this leads to \[\begin{align} \frac{\partial f}{\partial x}(\overline{x},\overline{u},\overline{d})&=\overline{M}^{-1}(A_0+\sum_{i=1}^{|\mathcal{F}|} \overline{u}_iA_i),\tag{11}\\ \frac{\partial f}{\partial u}(\overline{x},\overline{u},\overline{d})&=\overline{M}^{-1}\left(\begin{bmatrix}A_1\overline{x} \ldots A_{|\mathcal{F}|}\overline{x}&0\end{bmatrix}+E_u\right)\nonumber \\ &=\overline{M}^{-1}\begin{bmatrix} (B)_{v\in\mathcal{V}_h}F^\top &0\\ \hat{A}(\overline{x}) & \begin{bmatrix} 0\\I \end{bmatrix} \end{bmatrix}. \tag{12} \end{align}\]

Furthermore, since \(\ker F(B)_{v\in\mathcal{V}_h}^\top=\{0\}\) holds by Assumption 3, the right inverse \(W\) of \((B)_{v\in\mathcal{V}_h}F^\top\) exists, i.e., \((B)_{v\in\mathcal{V}_h}F^\top W=I\). Then, for some \(\varepsilon>0\), we define the following state feedback \[\begin{align} \label{eq:feedback} \tilde{u}(k) &=\begin{bmatrix} q_c-\overline{q}_c\\ P_{pr}-\overline{P}_{pr}\end{bmatrix}=G \tilde{x}(k)\\ &=\begin{bmatrix} -\varepsilon W &\varepsilon W(\hat{A}(\overline{x})W)^\top \\ 0&0 \end{bmatrix}\begin{bmatrix} m_{\mathrm{h}}(k)-\overline{m}_{\mathrm{h}} \\ \begin{bmatrix} T_v(k)-\overline{T_v}\\ T_e(k)-\overline{T_e} \end{bmatrix}\end{bmatrix}. \nonumber \end{align}\tag{13}\] Using 10 together with 11 , 12 and 13 , leads to \[\tilde{x}(k+1)=\tilde{x}(k)+\Delta t\overline{M}^{-1}\begin{bmatrix} K_1 &-K_2^\top \\ K_2 & K_3 \end{bmatrix} \tilde{x}(k)=: A_d\tilde{x}(k),\] with \(K_1:=-\varepsilon I\), \(K_2:= -\varepsilon\hat{A}(\overline{x}) W\) as well as \[\begin{align} K_3:=-\sum_{i=1}^{|\mathcal{F}|}\overline{u}_i\hat{A}_i -\tilde{A}_0+\varepsilon (\hat{A}(\overline{x}))W(\hat{A}(\overline{x})W)^\top. \end{align}\] Then, \(K_1\) is negative definite and by choosing \(\varepsilon\) sufficiently small, it follows from [20] that \(K_3\) fulfills \(z^\top( K_3+K_3^\top)z\leq 0\). Consider now \[\begin{align} \label{eq:discrete95lyap} &~~~~A_d^\top \overline{M}A_d-\overline{M} \\ &=2\Delta t\begin{bmatrix} K_1&0\\0& K_3 \end{bmatrix} + \Delta t^2\begin{bmatrix} K_1&K_2^\top\\-K_2& K_3 \end{bmatrix}\overline{M}^{-1}\begin{bmatrix} K_1&-K_2^\top\\K_2& K_3 \end{bmatrix}. \nonumber \end{align}\tag{14}\] If \(\Delta t>0\) is sufficiently small, then it follows again from [20], that 14 is negative definite and therefore the closed-loop system \(\tilde{x}(k+1)=A_d\tilde{x}(k)\) has the Lyapunov function \({V(\tilde{x})=\tilde{x}^\top\overline{M}\tilde{x}}\) with the discrete derivative along the solutions \({\Delta V(\tilde{x}(k))=\tilde{x}(k)^\top(A_d^\top \overline{M}A_d-\overline{M})\tilde{x}(k)}\). Hence, the feedback 13 stabilizes the linearization of 8 in \((\overline{x},\overline{u},\overline{d})\). Therefore, Condition 3 is fulfilled. ◻

Thus, if the DHG fulfills \(\ker F(B)_{v\in\mathcal{V}_h}^\top=\{0\}\) and if the step size \(\Delta t\) is sufficiently small, then the closed-loop system’s equilibrium resulting from applying the control law \(\mu_{N}(\cdot)\) is asymptotically stable.

4 Case Study↩︎

In this section, we demonstrate the practicability of the MPC approach outlined in Section  3 via a numerical case study. For this, the exemplary DHG shown in Figure 1 is modeled as described in Section 2 and numerically implemented3 using julia programming language [21]. In what follows, we first describe the parametrization of the DHG model and define nominal operating conditions. Based on that, we define two steady states with suitable terminal ingredients for each steady state. Finally, we simulate the stabilization of the two setpoints via MPC with terminal ingredients.

The length of each edge is set to \(L = 500 \text{m}\), the total masses of the TESs are set to \(m_{\mathrm{tes,}1} = 50 \cdot 10^3\text{kg}\) and \(m_{\mathrm{tes,}2} = 30 \cdot 10^3\text{kg}\), respectively.The heat loss coefficient is set to \((\kappa_{\mathrm{v}})_v = (\kappa_{\mathrm{e}})_e = 0.2 \frac{\text{kJ}}{\text{K}\cdot\text{s}}\) for \(v \in \mathcal{V}\), \(e \in \mathcal{E}\). Furthermore, we set the density of water, the ambient temperature, the specific heat capacity of water and the pipe friction coefficient to \(\rho = 988.05 \frac{\text{kg}}{\text{m}^3}\), \(T_{\mathrm{a}} = 10.0^\circ\text{C}\), \(c_p = 4.18 \frac{\text{kJ}}{\text{kg}\cdot\text{K}}\), \(\lambda = 0.02\), respectively.

We define two steady state set points \((\overline{x}^{I},\overline{u}^{I})\) and \((\overline{x}^{II},\overline{u}^{II})\) for the the DHG dynamics 8 . In what follows, we use the superscript (\(I\) or \(II\)) to associate the variables to each of the set points. One set point is computed for \({P_{\mathrm{d,}1}^{I} = 1.5 \text{MW}}\), \({P_{\mathrm{d,}2}^{I} = 2.5 \text{MW}}\) and \({P_{\mathrm{d,}3}^{I} = 1 \text{MW}}\) with nominal temperatures at vertex \(v_1\) and \(v_6\) of \({T_{\mathrm{v,}1}^{\mathrm{nom,}I} = 45 ^\circ\text{C}}\) and \({T_{\mathrm{v,}6}^{\mathrm{nom,}I} = 75 ^\circ\text{C}}\), respectively. The pipe diameters are designed using the Darcy-Weisbach equation such that the pressure drop equals \(300~\text{Pa}/\text{m}\) when operated at \((\overline{x}^{I},\overline{u}^{I})\). The second set point is computed for \({P_{\mathrm{d,}1}^{II} = 1.5 \text{MW}}\), \({P_{\mathrm{d,}2}^{II} = 2 \text{MW}}\) and \({P_{\mathrm{d,}3}^{II} = 1 \text{MW}}\) with nominal temperatures at vertex \(v_1\) and \(v_6\) of \({T_{\mathrm{v,}1}^{\mathrm{nom,}II} = 46 ^\circ\text{C}}\) and \({T_{\mathrm{v,}6}^{\mathrm{nom,}II} = 77 ^\circ\text{C}}\), respectively. With the given parameters, the heat losses are 9% and 8% in the stationary states \((\overline{x}^{I},\overline{u}^{I})\) and \((\overline{x}^{II},\overline{u}^{II})\), respectively.

Note that the Conditions 3 (i)-(vi) are satisfied for the described setup and the steady states \((\overline{x}^{I},\overline{u}^{I})\) and \((\overline{x}^{II},\overline{u}^{II})\). Furthermore, Assumption 3 holds, and it can be verified that Proposition 4 is applicable for \(\Delta t = 60\text{s}\). Therefore, the Condition 3 (vii) holds, which means that the exemplary DHG is stabilizable in the sense of Definition 2. Hence, the linear quadratic regulator (LQR) can be used to derive suitable terminal ingredients [13]. The prediction horizon is set to \({N = 60}\). Then, terminal ingredients can be derived via an auxiliary nonlinear optimization problem [13]. Without going into the details of this optimization problem, we point out that due to the nonlinearity of the DHG model, it is not possible to make a statement about the globality of the optimizer without further investigation, which is beyond the scope of this paper. Instead, we solve the problem iteratively with 100 different normally distributed initial values and assume that in this way we can approximate the global solution of the optimizer with sufficient accuracy.

The procedure results in one ellipsoid per set point. Each ellipsoid describing a terminal region (and a terminal cost) has dimension 18 and is defined by the solution \(P^{(\cdot)}\) of the Lyapunov equation for the linearized closed-loop system. We denote each solution by \(P^{I}\) and \(P^{II}\), with corresponding \(\alpha^{I}\) and \(\alpha^{II}\), respectively; see 9 .

In order to give a rough idea of the size of the terminal region, we identify the smallest box containing each ellipsoid and compute its one-dimensional projection on the coordinates corresponding to the states of the TESs. Here, we select all states associated with the TESs as states of particular interest because solely the TESs allow for a variable mass ratio of hot to cold water and have the largest inertia within the DHG due to their large total masses. In Table 1, we display the distance with respect to the associated steady state value representing the projections of the approximate terminal region of TES 1 and TES 2 onto the masses \(\Delta m_{\mathrm{v,}1} \in \mathbb{R}\), \(\Delta m_{\mathrm{v,}6} \in \mathbb{R}\) and temperatures \(\Delta T_{\mathrm{v,}1} \in \mathbb{R}\), \(\Delta T_{\mathrm{v,}6} \in \mathbb{R}\), respectively.

Table 1: Values approximating the terminal region for states associated with TESs.
\(\Delta m_{\mathrm{v,}1}\) \(\Delta T_{\mathrm{v,}1}\) \(\Delta m_{\mathrm{v,}6}\) \(\Delta T_{\mathrm{v,}6}\)
I \(0.75 \text{t}\) \(2.05 ^\circ\text{C}\) \(0.45 \text{t}\) \(0.4 ^\circ\text{C}\)
II \(0.75 \text{t}\) \(3.3 ^\circ\text{C}\) \(0.45 \text{t}\) \(2.07 ^\circ\text{C}\)

In order to assess the closed-loop system’s performance over \(N_{\text{sim}} \in \mathbb{N}\) discrete-time steps in the face of a sudden set point change at the discrete-time step \(k_{\mathrm{step}} \in \mathbb{N}\), we use the calculated set points and define two different controllers, namely MPC 1 and MPC 2. We use the stage cost, terminal cost and terminal region \[\begin{align} \ell_1(x,u)&=\begin{cases} \|x-\overline{x}^I\|_Q^2 + \|u-\overline{u}^I\|_R^2,~\text{if}~k \in \mathcal{I}_1, \\ \|x-\overline{x}^{II}\|_Q^2 + \|u-\overline{u}^{II}\|_R^2,~\text{if}~k \in \mathcal{I}_2, \end{cases} \\ C_1(x)&=\begin{cases} C_{P^I}(x),~\text{if}~k \in \mathcal{I}_1, \\ C_{P^{II}}(x),~\text{if}~k \in \mathcal{I}_2, \end{cases} \\ \mathbb{X}_{C,1}&=\begin{cases} \mathbb{X}_{C_{P^I},\alpha^I},~\text{if}~k \in \mathcal{I}_1, \\ \mathbb{X}_{C_{P^{II}},\alpha^{II}},~\text{if}~k \in \mathcal{I}_2, \end{cases} \end{align}\] respectively, where \(\mathcal{I}_1 = \{0,\ldots,k_{\mathrm{step}}-1\} \subset \mathbb{N}\) and \(\mathcal{I}_2 = \{k_{\mathrm{step}},N_{\mathrm{sim}}\}\subset \mathbb{N}\), to define MPC 1. Thus, for \(k \in \mathcal{I}_1\), MPC 1 aims to stabilize \((\overline{x}^{I},\overline{u}^{I})\). Then, for \(k \in \mathcal{I}_2\), \((\overline{x}^{II},\overline{u}^{II})\) is stabilized. To obtain MPC 2, we use a piece-wise constant reference trajectory \[\label{eq:piece-wise32constant32reference} (\overline{x}(k),\overline{u}(k)) = \begin{cases} (\overline{x}^{I},\overline{u}^{I}), \text{ for } k \in \mathcal{I}_1, \\ (\overline{x}^{II},\overline{u}^{II}), \text{ for } k \in \mathcal{I}_2, \end{cases}\tag{15}\] to define the stage cost function \(\ell_2(k,x,u) = \|x-\overline{x}(k)\|_Q^2 + \|u-\overline{u}(k)\|_R^2\), together with piece-wise constant terminal cost and terminal region \[\label{eq:32C232X95C2} \begin{align} C_2(k,x)&=\begin{cases} C_{P^I}(x),~\text{if}~k \in \mathcal{I}_3, \\ C_{P^{II}}(x),~\text{if}~k \in \mathcal{I}_4, \end{cases} \\ \mathbb{X}_{C,2}(k)&=\begin{cases} \mathbb{X}_{C_{P^I},\alpha^I},~\text{if}~k \in \mathcal{I}_3, \\ \mathbb{X}_{C_{P^{II}},\alpha^{II}},~\text{if}~k \in \mathcal{I}_4, \end{cases} \end{align}\tag{16}\] respectively, where \(\mathcal{I}_3 = \{0,\ldots,k_{\mathrm{step}}-N-1\} \subset \mathbb{N}\) and \(\mathcal{I}_4 = \{k_{\mathrm{step}}-N,\ldots,N_{\text{sim}}\} \subset \mathbb{N}\). It follows, that \(C_2(k,x)\) and \(\mathbb{X}_{C,2}(k)\) suddenly change at time step \(k_{\mathrm{step}}-N\), whereas \(\ell_2(k,x,u)\) changes for \(k \in \mathcal{I}_1 \cup \mathcal{I}_4\). Thus, MPC 1 and MPC 2 only differ within the transition phase, i.e., for \(k \in \mathcal{I}_1 \cup \mathcal{I}_4\). Through the application of MPC 2, we investigate the ability of MPC in the context of DHG to utilize knowledge of future heat flow extraction and temperature set points, i.e., to incorporate a predictive behavior into the closed-loop DHG, during the transition phase.

We simulate the closed-loop DHG for \(3\text{h}\), i.e., for \(N_{\text{sim}} = 180\), and choose \(k_{\mathrm{step}} = 90\). The average time required to solve the optimization problem of MPC 1 and MPC 2 is \(0.12\text{s}\) and \(0.14\text{s}\), respectively. The maximum time required to solve the optimization problem of MPC 1 and MPC 2 is \(0.96\text{s}\) and \(0.42\text{s}\), respectively. The results of the described simulation are shown in Figure 3. It can be observed that the DHG under MPC 1 is steered towards \((\overline{x}^{I},\overline{u}^{I})\) for \(k \in \mathcal{I}_1\) and that all desired set points which are indicated by the black dotted lines are reached before the set point changes. Afterwards, for \(k \in \mathcal{I}_2\), MPC 1 steers the DHG towards \((\overline{x}^{II},\overline{u}^{II})\) so that all set points are reached latest after \(2.75\text{h}\). This illustrates the theoretical results on asymptotic stabilization of constant set points from Section 3.

For MPC 2, it can be observed that the DHG is steered towards \((\overline{x}^{I},\overline{u}^{I})\) for \(k \in \mathcal{I}_3\). As soon as the second steady state lies within the prediction horizon of the MPC, i.e., for \(k \in \mathcal{I}_4\), the terminal region and terminal cost of the optimization problem of MPC 2 change according to 16 . Additionally, with each iteration for \(k \in \mathcal{I}_1 \cup \mathcal{I}_4\) the stage cost function is updated according to 15 . Therefore, we can see that the closed-loop behavior of the DHG already adapts for \(k \in \mathcal{I}_1 \cup \mathcal{I}_4\), which is manifested in the form of a increase of the temperatures of TES 1 and TES 2 as well as an increase of mass stored in TES 1 and TES 2 towards \((\overline{x}^{II},\overline{u}^{II})\). Using MPC 2, the steady state \((\overline{x}^{II},\overline{u}^{II})\) is reached after approximately \(t = 2.5 \text{h}\). Compared to MPC 1, MPC 2 reduces required control input changes in between successive control actions in the required heat flows by up to \(3.3\) times. This indicates that the incorporation of a stage cost of the form of \(l_2(k,x,u)\) is beneficial for smooth control input trajectories.

Additionally, the successful incorporation of control input constraints can be observed for the injected heat flow \(P_{\mathrm{pr,}1}^{\mathrm{MPC}1}\), \(P_{\mathrm{pr,}2}^{\mathrm{MPC}1}\) and \(P_{\mathrm{pr,}2}^{\mathrm{MPC}2}\) for \(t \geq 1.5 \text{h}\) since the heat flow injection reaches its limit, but does not overshoot. Remaining constraints also hold but are never active.

To summarize, this case study demonstrates the practicability of the MPC approach using terminal ingredients for a realistic small scale DHG for two different steady states. In this process, two MPCs with terminal ingredients are compared that coordinate two heat producers so that three consumers can draw heat from the DHG at desired temperature levels. Both MPCs stabilize given piece-wise constant set points and satisfy constraints. Additionally, the investigated slightly extended version of the MPC from Section 3, i.e. MPC 2, showed a predictive behavior and beneficial performance using piece-wise constant stage cost, terminal cost and terminal region. This allowed us to emphasize the ability to use information about future heat demand and temperature requirements in the context of an MPC.

Figure 3: Masses and temperatures at TESs with corresponding mass flow rates as well as heat flow injections for MPC 1 and MPC 2. Set points are indicated by the dotted lines.

5 Conclusion↩︎

Given the significance of developing operating strategies for RES-based DHGs, we have identified MPC as a promising strategy that effectively amalgamates the benefits of existing optimization-based EMSs with the strengths of established control approaches for DHGs. To this end, we derived an ODE-based model describing the thermo-hydraulics of DHGs containing TESs with variable storage mass ratio as well as Kirchhoff’s laws for hydraulic networks. We proved the stabilizability of the DHG model, which is a necessary system property to apply MPC with terminal ingredients. Via a case study, we demonstrated asymptotic stabilization of different steady state set points of an exemplary DHG.

In future work, we plan to investigate the impact of incorporating a piece-wise constant stage cost function into the MPC scheme and examine the scalability of our approach.


This research received funding from the German Federal Government, the Federal Ministry of Education and Research, and the State of Brandenburg within the framework of the joint project EIZ: Energy Innovation Center (project numbers 85056897 and 03SF0693A).


H. Lund, F. Hvelplund, P. Østergaard, B. Möller, B. V. Mathiesen, D. Connolly, and A. N. Andersen, “Chapter 6 - analysis: Smart energy systems and infrastructures,” in Renewable Energy Systems (Second Edition), second edition ed., H. Lund, Ed.1em plus 0.5em minus 0.4emAcademic Press, 2014.
D. Trier, C. Rothballer, G. Stiff, and B. V. Mathiesen, “Guidelines for the EnergySystemTransition: TheEnergyUnionPerspective - HeatRoadmapEurope,” 2018.
H. Lund, S. Werner, R. Wiltshire, S. Svendsen, J. E. Thorsen, F. Hvelplund, and B. V. Mathiesen, “4th generation district heating: Integrating smart thermal grids into future sustainable energy systems,” Energy, vol. 68, 2014.
M. Rose, C. A. Hans, and J. Schiffer, “A predictive operation controller for an electro-thermal microgrid utilizing variable flow temperatures,” in 22nd IFAC World Congress, Yokohama, Japan, 2023.
R. Krug, V. Mehrmann, and M. Schmidt, “Nonlinear optimization of district heating networks,” Optim. Eng., vol. 22, no. 2, 2021.
M. Rein, J. Mohring, T. Damm, and A. Klar, “Optimal control of district heating networks using a reduced order model,” Optimal Control Applications and Methods, vol. 41, no. 4, 2020.
J. E. Machado, M. Cucuzzella, N. Pronk, and J. M. A. Scherpen, “Adaptive control for flow and volume regulation in multi-producer district heating systems,” IEEE Control Syst. Lett., vol. 6, 2022.
J. Bendtsen, J. Val, C. Kallesøe, and M. Krstic, “Control of DistrictHeatingSystem with Flow-dependent Delays,” IFAC-PapersOnLine, vol. 50, no. 1, 2017.
T. Scholten, C. De Persis, and P. Tesi, “Modeling and control of heat networks with storage: The single-producer multiple-consumer case,” in 2015 European Control Conference, 2015.
J. E. Machado, M. Cucuzzella, and J. M. A. Scherpen, “Modeling and passivity properties of multi-producer district heating systems,” Automatica, vol. 142, 2022.
A. Krishna and J. Schiffer, “A port-Hamiltonian approach to modeling and control of an electro-thermal microgrid,” IFAC-PapersOnLine, vol. 54, no. 19, 2021.
L. Grüne and J. Pannek, Nonlinear Model Predictive Control : Theory and Algorithms. 2nd Edition.1em plus 0.5em minus 0.4emSpringer, 2017.
H. Chen and F. Allgöwer, “,” , vol. 34, no. 10, 1998.
C. Rajhans, S. C. Patwardhan, and H. Pillai, “Discrete time formulation of quasi infinite horizon nonlinear model predictive control scheme with guaranteed stability,” IFAC-PapersOnLine, vol. 50, no. 1, 2017.
R. Borsche, M. Eimer, and N. Siedow, “A local time stepping method for thermal energy transport in district heating networks,” Applied Mathematics and Computation, vol. 353, 2019.
S.-A. Hauschild, N. Marheineke, V. Mehrmann, J. Mohring, A. Moses Badlyan, M. Rein, and M. Schmidt, “Port-Hamiltonian modeling of district heating networks,” Progress in Differential-Algebraic Equations II, 2020.
B. Bollobás, Modern Graph Theory, ser. Graduate Texts in Mathematics.1em plus 0.5em minus 0.4emSpringer, 1998, vol. 184.
C. De Persis and C. S. Kallesoe, “Pressure Regulation in NonlinearHydraulicNetworks by Positive and QuantizedControls,” IEEE Transactions on Control Systems Technology, vol. 19, no. 6, 2011.
C. Heij, A. C. M. Ran, and F. van Schagen, Introduction to Mathematical Systems Theory: Discrete Time Linear Systems, Control and Identification.1em plus 0.5em minus 0.4emCham: Springer International Publishing, 2021.
J. Schiffer, F. Dörfler, and E. Fridman, “Robustness of distributed averaging control in power systems: Time delays & dynamic communication topology,” Automatica, vol. 80, 2017.
J. Bezanson, S. Karpinski, V. B. Shah, and A. Edelman, “Julia: AFastDynamicLanguage for TechnicalComputing,” Sep. 2012.

  1. \(^{1}\)Fraunhofer Research Institution for Energy Infrastructures and Geothermal Systems IEG, 03046 Cottbus, Germany, {max.rose, hannes.gernandt, johannes.schiffer}↩︎

  2. \(^{2}\)Brandenburg University of Technology Cottbus-Senftenberg BTU, 03046 Cottbus, Germany, {machadom, schiffer}↩︎

  3. The source code can be accessed in the following repository:↩︎