=.75mm

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.

Notation↩︎

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.

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

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

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}\]

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.

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.

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

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.

Furthermore, recall [10].

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 implemented³ 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.

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.

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.

Acknowledgment↩︎

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).

References↩︎