Optimal control of district heating networks using a reduced order model

TL;DR

This paper presents a reduced order model for optimal control of district heating networks, significantly decreasing computation time while managing complex dynamics to optimize energy feed-in power.

Contribution

The paper introduces a novel reduced order modeling approach based on system theory for efficient optimal control of large-scale heating networks.

Findings

Reduced model decreases optimization computation time.

Effective control achieved despite complex network dynamics.

Validated on a large-scale heating network with flux direction changes.

Abstract

We study the optimal control of district heating networks using a reduced order model based on a system theoretic description close to the underlying Euler equations. In the presented scenarios, the central task is to limit the maximal feed-in power occurring as a product of control and state variables. The underlying dynamics of heating networks acting as optimization constraints pose the central computational complexity, prohibiting the determination of an optimal control online. The advection of the injected energy density on the network results in an index-1, quadratic in state differential algebraic equation, challenging to reduce. The suggested reduced model decreases the computation time of the optimization significantly. The effectiveness of the presented approach is demonstrated for an existing, large-scale heating network including changes of flux directions.

Tables

Table 1. Table 1: Outline data for the flow part of the heating network presented in fig.  1 .

Edges	Nodes	Loops	Pipelines	Consumers	Pipeline length
1108	770	6	775	333	8676 m

Table 2. Table 2: Description of the considered optimization scenarios. Test cases TC1-TC3 differ by the considered daily mean temperature changing from − 3 ∘ ​ C superscript 3 C -3\,^{\circ}\text{C} to 7.5 ∘ ​ C superscript 7.5 C 7.5\,^{\circ}\text{C} . T min cons , T max net superscript subscript 𝑇 min cons subscript superscript 𝑇 net max T_{\text{min}}^{\text{cons}},T^{\text{net}}_{\text{max}} denote the temperature equivalents of the optimization constraints ( 17 , 16 ).

$T_{\text{min}}^{\text{cons}},T_{\text{max}}^{\text{net}}{/}^{\circ }\text{C}$	$p_{\text{min}}^{\text{h}}$, $p_{\text{max}}^{\text{h}}$, $\overline{\Delta p}/$bar	$\overline{P}$	$T_d{/}^{\circ }\text{C}$	$t_0,t_e/h$	dt/s
75, 110	3.5, 9.1, 2.5	$0.5\left(\overline{G}+G_0\right)$	$-3,3,7.5$	0,72	300

Table 3. Table 3: Optimal controls and runtime comparison of TC1 ( − 3 ∘ ​ C superscript 3 C -3\,\,^{\circ}\text{C} ) for varying number of state variables (DOF) including full and reduced order models. The reference control u ~ ~ 𝑢 \tilde{u} results from the optimal control strategy described in section 5 using the reference discretization FOMR. The feed-in P 𝑃 P is measured by FOMR  based on the control suggested by each coarse model. The last row measures the relative error of outputs y 𝑦 y comparing reference model FOMR, and each coarse model.

	FOM 0	FOM 1	ROM 1	FOMR
runtime opt./s	2765.0	2640.0	608.5	16100.0
runtime sim./s	2685.0	2585.0	557.0	15900.0
# solves DAE	16	9	11	12
DOF	775	1789	180	9538
$J\left(u\right)$	55.9	53.0	52.9	52.6
${\parallel u-\tilde{u}\parallel }_2/{\parallel \tilde{u}\parallel }_2$	8.10$\times {10}^{-3}$	9.48$\times {10}^{-4}$	7.25$\times {10}^{-4}$	0
${\parallel P-\overline{P}\parallel }_{\text{max}}/\overline{P}$	5.41$\times {10}^{-3}$	4.60$\times {10}^{-4}$	4.49$\times {10}^{-4}$	0
${\max }_h{\parallel y^h-y_R^h\parallel }_2/{\parallel y_R^h\parallel }_2$	3.03$\times {10}^{-2}$	1.23$\times {10}^{-3}$	5.28$\times {10}^{-3}$	0

Table 4. Table 4: Optimization results and runtime comparison of TC 2( + 3 ∘ ​ C superscript 3 C +3\,\,^{\circ}\text{C} ) for different spatial discretizations including full and reduced order models. For a detailed explanation we refer the reader to table  3 .

	FOM 0	FOM 1	ROM 1	FOMR
runtime opt./s	2525.0	3210.0	596.2	18300.0
runtime sim./s	2455.0	3146.7	546.2	18100.0
# solves DAE	15	11	11	13
DOF	775	1789	180	9538
$J\left(u\right)$	16.5	14.8	14.8	14.6
${\parallel u-\tilde{u}\parallel }_2/{\parallel \tilde{u}\parallel }_2$	7.10$\times {10}^{-3}$	6.17$\times {10}^{-4}$	5.26$\times {10}^{-4}$	0
${\parallel P-\overline{P}\parallel }_{\text{max}}/\overline{P}$	8.11$\times {10}^{-3}$	6.86$\times {10}^{-4}$	6.80$\times {10}^{-4}$	0
${\max }_h{\parallel y^h-y_R^h\parallel }_2/{\parallel y_R^h\parallel }_2$	2.60$\times {10}^{-2}$	1.25$\times {10}^{-3}$	4.72$\times {10}^{-3}$	0

Table 5. Table 5: Optimization results and runtime comparison of TC 3( + 7.5 ∘ ​ C superscript 7.5 C +7.5\,\,^{\circ}\text{C} ) for different spatial discretizations including full and reduced order models. For a detailed explanation we refer the reader to table  3 .

	FOM 0	FOM 1	ROM 1	FOMR
runtime opt./s	3165.0	5930.0	801.8	16400.0
runtime sim./s	3070.0	5795.0	735.2	16200.0
# solves DAE	19	21	15	13
DOF	775	1789	180	9538
$J\left(u\right)$	6.7	5.6	5.6	5.5
${\parallel u-\tilde{u}\parallel }_2/{\parallel \tilde{u}\parallel }_2$	6.16$\times {10}^{-3}$	6.14$\times {10}^{-4}$	9.63$\times {10}^{-4}$	0
${\parallel P-\overline{P}\parallel }_{\text{max}}/\overline{P}$	1.46$\times {10}^{-2}$	1.34$\times {10}^{-3}$	1.37$\times {10}^{-3}$	0
${\max }_h{\parallel y^h-y_R^h\parallel }_2/{\parallel y_R^h\parallel }_2$	2.35$\times {10}^{-2}$	1.57$\times {10}^{-3}$	5.01$\times {10}^{-3}$	0

Table 6. Table 6: Maximal population density determined in the simulation of system operator and Jacobian for different discretizations. DOF denotes the number of differential state variables.

	FOM 0	FOM 1	ROM 1	FOMR
# nonzero entries $A\left(v\right)$	1555	3583	6827	19081
# nonzero entries $J_{\dot{e}}(v,e)$	67837	139922	11211	662200
DOF	775	1789	180	9538

Equations

$$0$$

$${\dot{e}}^{\alpha }$$

$$\sum _{\alpha \in {\delta }^{\text{out}}(N)}{q}^{\alpha }(t)e^{\alpha }(t,0)$$

$$e^{\alpha }(t,0)$$

$$\sum _{\alpha \in {\delta }^{\text{out}}(N)}{q}^{\alpha }(t)$$

$$\sum _{\alpha \in O}\frac{{\lambda }^{\alpha }\rho l^{\alpha }}{2d^{\alpha }}|v^{\alpha }(t)|v^{\alpha }(t)$$

$$\displaystyle\big{[}q^{i}(t)\cdot(e^{i}(t,0)-e_{R})\big{]}_{i\in\mathcal{H}}$$

$$e^s(t,l^s)$$

$$p^s(t,l^s)$$

$$\dot{e}$$

$$y$$

$$0$$

$$A(v)$$

$$B(v)$$

$${\dot{e}}_r$$

$$y_r$$

$$0$$

$$J(u_T^{\kappa })$$

$$u_T^{\kappa }(t)$$

$$e_{i,n}(t)$$

$$p^i(t)$$

$$(u_T^{\kappa }(t)-e_R)\sum _{i\in \mathcal{H}}{q}^i(t)$$

$$D(\dot{e}(t),e(t),e_0,q(t),u_T^{\kappa },G(t))$$

$$P_{\text{hyd}}(t)=\Delta p^s(t)\sum _{i\in \mathcal{H}}\frac{G_i(t)}{e_i(t)-e_R}\le \frac{{\max }_{t\in {\mathcal{T}}_c}(\Delta p^s(t))}{{\min }_{t\in {\mathcal{T}}_c}e(t)}\underset{t\in {\mathcal{T}}_c}{\max }\sum _{i\in \mathcal{H}}{G}_i(t),$$

$$\Delta p^h\equiv p^h-p^s=\rho g(z^s-z^h)-\rho \sum _{i\in K_h}{\lambda }^i\frac{l^i}{2d^i}{v}^i|v^i|,$$

$$\underset{h\in \mathcal{H}}{\max }({\tilde{p}}^h(t))-\underset{h\in \mathcal{H}}{\min }({\tilde{p}}^h(t))\le \overline{\Delta p}\le p_{\text{max}}^{\text{h}}-p_{\text{min}}^{\text{h}}t\in {\mathcal{T}}_d,$$

$$u_p(t)=p_{\text{min}}^{\text{h}}-\underset{h\in \mathcal{H}}{\min }(\Delta p^h(t))t\in {\mathcal{T}}_d.$$

$$\underset{h\in \mathcal{H}}{\max }({\tilde{p}}^h-{\tilde{u}}_p)-\underset{h\in \mathcal{H}}{\min }({\tilde{p}}^h-{\tilde{u}}_p)=\underset{h\in \mathcal{H}}{\max }(\Delta p^h)-\underset{h\in \mathcal{H}}{\min }(\Delta p^h)\le p_{\text{max}}^{\text{h}}-p_{\text{min}}^{\text{h}},$$

$$p^h=p^s+\Delta p^h.$$

$$\underset{h\in \mathcal{H}}{\min }(p^h)=p^s+\underset{h\in \mathcal{H}}{\min }(\Delta p^h)=p_{\text{min}}^{\text{h}}-\underset{h\in \mathcal{H}}{\min }(\Delta p^h)+\underset{h\in \mathcal{H}}{\min }(\Delta p^h)=p_{\text{min}}^{\text{h}}.$$

$$\underset{h\in \mathcal{H}}{\max }(p^h)$$

$$=p_{\text{min}}^{\text{h}}-\underset{h\in \mathcal{H}}{\min }(\Delta p^h)+\underset{h\in \mathcal{H}}{\max }(\Delta p^h)$$

$$\le p_{\text{min}}^{\text{h}}+p_{\text{max}}^{\text{h}}-p_{\text{min}}^{\text{h}}$$

$$=p_{\text{max}}^{\text{h}},$$

$$P$$

$$=(u_T^{\kappa }(t)-e_R)\sum _{i\in \mathcal{H}}\frac{G_i(t)}{e_{i,n_i}(t)-e_R}.$$

$$D(\dot{e}(t),e(t),e_0,q(t),u_T^{{\kappa }_i},G(t))=0,t\in {\mathcal{T}}_d.$$

$$\tilde{f}(x_{\mu },x_{\mu +1},u)=0.$$

$${\partial }_{\kappa }{x}_{\mu +1}$$

Full text

Keywords. energy networks; optimal control; model order reduction; linear time varying system; Galerkin projection; district heating networks; process optimization.

1 Introduction

In this contribution an optimal control problem for district heating networks is solved using reduced order models. Heating networks are of particular interest for low-carbon energy supply due to their flexibility in using different sources of energy [1, 2]. The energy density $u_{T}$ injected at a power plant is guided to consumers of different sizes using a network of pipelines referred to as flow network. At the consumers, the local volume flow is regulated using heat exchangers to match the time dependent power consumption $G$ given the currently available energy density $e$. Fig. 1 illustrates an existing large scale network considered in this contribution. Its outline data is supplied in tab. 1. A central aim of operating these networks lies in efficiently planning the input energy density $u_{T}$. It defines the power feed-in $P=(u_{T}-e_{R})\hat{q}$ in combination with the aggregated volume flow $\hat{q}$, and the energy density of the cooled fluid $e_{R}$ entering the plant in the return network. Due to the high transport times from source to consumption in large scale networks, the power feed-in $P(t)$ at each point of time $t$ needs not to match the current power consumption $G$ at the heat exchangers. While both quantities are coupled by conservation of energy within the heating network, there exists an essential optimization potential in distributing $u_{T}$ over time. Injecting a high energy density in times of large volume flows $\hat{q}$ requires firing additional vessels which might be unfavorable for economic and ecological reasons. Vice versa, planning might also enable to use external overcapacities of energy resulting from renewable energies. In the specific application, a waste to heat incineration plant is able to deliver power up to a maximum level at no cost. When exceeding this level for only a short time interval, significant costs result from using gas boilers to cover peaks in the injected power.

In the optimization task, the dynamical energy transport on these complex networks is an essential constraint to ensure a physically relevant control and requires a significant amount of computation time. Since controls of power plants are updated every $15\,$min in application, the corresponding optimization needs to be sufficiently fast. To this end, the usage of reduced order models is central [3, 4, 5]. For network systems, both tasks of defining an optimal control [6], as well as formulating reduced order models are complex. Due to their mathematical complexity and the large benefits, the optimization of energy networks such as electric [7], water [8], gas and heating networks is an active field of research. Concerning the model order reduction of energy systems, different works already exist for gas networks [9, 10], electric networks [11] and other applications [12]. Similarly, many publications focus on the optimization of gas networks [13, 14, 8, 15, 16, 17]. While gas networks are modeled by compressible Euler equations, the transport fluid for heating networks is water in the fluid phase inducing incompressible dynamics. On the one hand this simplifies the equations for the conservation of mass and momentum significantly. On the other hand, for heating networks the conservation of inner energy describing the dynamical transport of thermal energy is the dominant effect and mathematically the most challenging one. It leads to a time-dependent advection of the injected thermal energy yielding a large, dynamically changing delay in the energy transport between source and consumption. Towards the formulation of a reduced order model, advection on network systems is a demanding problem.

Selected works discussing the optimal control of heating networks are mentioned subsequently. In [18], pumping costs resulting from the variation of pressure and massflow are optimized using a reduced order model without reflecting the dynamically changing thermal transport. In [19], a predictive controller is formulated. The advection of thermal energy is reflected by the method of characteristics on each pipeline assuming a constant time delay between source and households. In [20] the transport of thermal energy from source to each consumer is described by virtual single pipelines including a constant time delay from source to consumption points. In this contribution, similar to the approach used for gas networks, a model very close to the underlying Euler equations is formulated, allowing to precisely model the thermal transport dynamics and addressing the central difficulties of heating networks. These are dynamically changing delay times from source to consumers, as well as changes of flux direction. Using a system theoretic approach to model the dynamics allows to use effective tools from model order reduction.

The article is structured as follows. After introducing the model for the energy transport within heating networks in section 2, a reduced order model required to efficiently simulate the relevant outputs of the heating network, is described in section 3. Hereafter, the optimal control problem is stated in section 4, and the numerical determination of an optimal control is discussed in section 5. Finally, the application of the reduced order model for determining an optimal control is numerically investigated in section 6 and compared to different fidelities of full order models.

2 Model for the transport dynamics

In this section, we suggest a system description of the dynamics of heating networks beneficial for the formulation of a reduced order model. The assumptions used in the derivation are the following. Since water is the transport medium, the incompressible limit is used. Furthermore, the pipelines transporting the fluid are assumed to be perfectly isolated avoiding thermal losses. In agreement with the incompressibility assumption, the density is assumed to be constant. Since for heating networks the pressure difference introduced by friction at the pipeline border dominates the pressure difference induced by acceleration, inertial effects are neglected. Finally, the return temperatures of all heat exchangers are assumed to be equal and constant, in line with the available data. Since these return temperatures define the input for the return network, this allows to neglect the return network when modeling heating networks, because it contains no observable relevant for the dynamics. Under these assumptions, the Euler equations of the transported fluid within a pipeline $\alpha\in\mathcal{P}$ read

[TABLE]

where $\mathcal{P}$ denotes the set of pipes including $n_{P}$ elements. While the conservation of mass (1) degenerates to the incompressibility assumption, the conservation of momentum was integrated along the pipeline length $l^{\alpha}$ to obtain (2), introducing pressure and height differences $\Delta p$, $\Delta z$ over pipeline borders. The loss of pressure by friction with the pipeline border is modeled using the Darcy-Weisbach equation [21]. It couples to a dimensionless friction factor $\lambda^{\alpha}$ depending on the properties of pipeline $\alpha$. The change of the friction factor with Reynolds number is neglected, allowing to set $\lambda^{\alpha}$ constant in time within this contribution. The remaining quantities are the pipeline diameter $d^{\alpha}$, the density $\rho$, and the gravitational constant $g$. Eq. (3) describes the advection of the energy density $e$ subject to the velocity $v$. While the energy density at the start of the pipe is defined by the inflow, the velocity is defined by the consumers stations entering as additional network components. Transferring the pipeline model to the network context, the resulting algebraic constraints at the nodes are given by

[TABLE]

In the presented network description, $\mathcal{P}$ denotes the set of pipelines, $\mathcal{H}$ the set of heat exchangers or consumers. Each element in $\mathcal{L}$ represents a set of paths forming a network loop. The presented algebraic equations describe the conservation of energy (4) and volume (6) over node $N$, where $\delta^{\text{in}}(N)$ ($\delta^{\text{out}}(N)$) denotes edges entering (exiting) node $N$. Eq. (5), claims that different incoming energy densities instantly mix to a defined outgoing energy density identical for all edges exiting node $N$. Eq. (7) results from claiming continuity of the pressure, what defines pressure levels at each node $N$. Summing the resulting pressure differences along loop $O$ results in (7). As the second of the two flow defining equations, (8) defines how the power consumption $G$ is covered by the difference in energy density at each house and its volume flow. Each of the corresponding heat exchangers ensures by regulation of the volume flow $q$ that each consumer receives the required power consumption $G$ for an arbitrary energy difference within the allowed range. Finally, (9) and (10) set the controls for energy density $u_{T}$ and the pressure level $u_{p}$ at the power plant entering the flow network.

For a numerical treatment of the partial differential algebraic equation (PDAE) (1-9), we perform a spatial discretization of the energy density $e$ using the upwind scheme [22]. To this end, $e_{\alpha,\beta}$ denotes cell $\beta$ on pipeline $\alpha\in\mathcal{P}$. This allows to formulate the differential algebraic equation (DAE) used to simulate the dynamics of the full order model,

[TABLE]

Here $e$ denotes the vector of energy densities of all finite volume cells on all pipelines. The operator $A\in\mathbb{R}^{n\times n}$ couples both finite volume cells locally on a pipeline and globally over every node including conservation of energy. The operator $B\in\mathbb{R}^{n\times 1}$ defines where the thermal input is applied to the network. The output matrix $C\in\mathbb{R}^{o\times n}$ measures observables at consumer stations. The system description (11) is Lyapunov stable [23] with known energy matrix $Q$. Hence, (11) can be reduced in a stability preserving way as described in section 3. Eq. (11) forms a DAE of index 1. For the rest of the contribution, we refer to (11) as full order model (FOM). Different full order models resulting from varying number of finite volume cells are evaluated towards their applicability in determining an optimal control in section 6. This includes the comparison to a reduced order model (ROM) formulated in the following section.

3 Formulation of a reduced order model

Subsequently, the formulation of a reduced order model reducing the computational cost of solving the dynamics of heating networks (11) in time is summarized. A detailed derivation of the reduction technique is given in [23].

As shown in [23], the velocity dependent system operators $A(v),B(v)$ defined in (11) allow for an affine decomposition with respect to the velocity,

[TABLE]

where $B^{v}_{i},A^{v}_{i}$ are time-independent matrices multiplied with the velocities defining the time-dependent weighting function $\gamma:\mathds{R}^{n_{F}^{v}}\rightarrow\mathds{R}^{n_{F}^{v}}$. Here, $n_{F}^{v}$ denotes the number of velocity configurations required to describe the dynamics. The function $\gamma_{i}$ equals the velocity $v^{i}$ if $A^{v}_{i}$ is valid regardless of the flux directions in the network. If in contrast $A^{v}_{i}$ is only valid for a specific flow direction, $\gamma_{i}$ is a nonlinear function. If no changes of the flux direction occur in the network, $\gamma_{i}=v^{i},\;\forall i\in\mathcal{P}$. The determination of all possible flux directions is performed in the offline-phase of the reduction process. This description allows to form a reduced model offline based on projection, where $V,W\in\mathbb{R}^{n\times r}$ denote Petrov-Galerkin projection matrices [24]. It can be shown that Lyapunov stability translates to the ROM, when implying the energy matrix $Q$ in the projection step [25, 26]. For the model of heating networks considered here, the construction of the corresponding energy matrix $Q$ is shown [23]. Specifically, volume conservation (6) defined by the current velocity field ensures Lyapunov stability of the FOM [23]. Hence, it is advantageous to rewrite the affine decomposition (12) in terms of independent volume flows $\tilde{q}$, what automatically ensures the volume conservation. This is possible, since (6) can be solved a-priori based on the network diameters, forming the set of independent volume flows $v=N\tilde{q}$. The resulting reduced order model replaces the full order DAE (21) in finding an optimal control and is given by

[TABLE]

Since system (14) changes dynamically, a global Galerkin projection has to be determined, reflecting all relevant linear submodels $\mathcal{Q}^{e}$ [27]. The latter are determined using a greedy-strategy [28]. Based on the current projection $V$, the linearization $v^{*}\in\mathcal{Q}^{\delta}\supset\mathcal{Q}^{e}$ is determined, exhibiting the largest error in the transfer function between full and reduced order model. Here, $\mathcal{Q}^{\delta}$ is the set of linearizations, for which the deviation is tested. It is determined based on training simulations spanning the range of scenarios for which the ROM should be applied hereafter. For the linearization $v^{*}$, the local Galerkin projection is performed and added to the global projection using a singular value decomposition. This process is repeated, until all linearizations of both full and reduced order model exhibit a relative error of the transfer function below $\Delta$. To obtain the local Galerkin projection, we use a moment-matching technique in frequency space [23]. Similar to IRKA [29], it forms a reduced order model interpolating the transfer function of the full order model for a fixed volume flow field at certain points in frequency space. The designed ROM should be valid for different environmental temperatures and only depend on the network topology. Thus, the time in the offline-phase to form the ROM is not taken into account.

4 Control problem

For the set of discrete points of time $\mathcal{T}_{d}=\{t_{0},...,t_{e}\}$, and the corresponding continuous interval $\mathcal{T}_{c}=[t_{0},t_{e}]$ we search a parameterized control $u_{T}^{\kappa}:\mathcal{T}_{c}\rightarrow\mathds{R}$ of the input energy density, minimizing the following objective

[TABLE]

subject to

[TABLE]

The objective function (15) penalizes the temporal variation and the distance to the regularization parameter $\eta_{2}$ of the parameterized control $u_{T}^{\kappa}$. The regularization parameters $\eta_{1},\eta_{2}$ are used to equalize both contributions in the objective which are motivated as follows. Minimizing the temporal variation leads to realistic controls, which can be realized properly by the power plant. By choosing $\eta_{2}$ sufficiently small, the mean value of the control decreases which systematically reduces thermal losses by cooling effects. Decreasing the injected energy density will lead to higher pressures in the network, which are restricted as explained below. Constraints (16-19) are technical restrictions in line with standard operation instructions formulated as optimization constraints. The upper energy limit (16) forces the fluid to be in the liquid phase. Note that energy densities in the network can never exceed the one at the inlet. In (17) a minimal energy density is required for proper operating conditions of the heat exchangers. In addition, the pressure levels of consumption nodes in the flow network are restricted to upper and lower bounds (18,19). Eq. (20) sets an upper bound to the maximal injected power, which avoids the use of additional energy sources. Below this limit, the power plant can supply demands by energy stemming from a waste to heat incineration plant at no costs. Due to the transport time of the injected energy from source to consumers, the control has a delayed effect on the consumer, allowing to influence the temporal distribution of the injected power. Finally, (21) reflects the energy transport (3) along the network restricted by the algebraic coupling conditions (4-10) and the initial state $e_{0}(u_{T}^{\kappa})$, cf. section 2. For full and reduced order models, (21) is replaced by (11) and (14) respectively.

The consumption $G$ is assumed to be known a-priori, which is a typical assumption in the simulation of heating networks. Since cooling effects are neglected and the energy densities in the return network are modeled equal and constant, an open loop control problem results. The return network exhibits a constant energy density $e_{R}$ entering the feed-in power (20) as a parameter.

The required pumping power $P_{\text{hyd}}$ necessary to retain a fluid flow in the network is stemmed by the pumps in the depot. It is bounded above as follows,

[TABLE]

where $\Delta p^{s}$ describes the pressure difference achieved at the source edge in the depot. For typical networks, the maximum pressure difference at the depot is smaller than $10\,$bar. Approximating the energy density by $e\approx\rho c_{p}T$, with material constants described in section 6, a maximum aggregated power consumption of $1\,$MW leads to a corresponding pumping power of $16\,$kW. Thus the pumping power is suppressed by almost two orders of magnitude compared to the thermal power.

Treatment of pressure constraints

Subsequently, we explain how to fulfill the pressure constraints (18, 19) in a simplified manner. It relies on limiting the difference of the maximum and minimum pressure levels measured at all consumption points. This allows to adjust the pressure control $u_{p}:\mathcal{T}_{c}\rightarrow\mathds{R}$ after finding the optimal control of the energy density.

By (2), the pressure difference from source to consumption point $h$ is defined by

[TABLE]

where $K_{h}$ denotes an arbitrary path from the source to consumer $h\in\mathcal{H}$. Although the pipeline velocities change dynamically, diameters, lengths and the height profile on the path from source to each household determine the resulting pressure difference to a large extent. This stabilizes the constraint limiting the maximum pressure difference and motivates the following proposition.

Proposition 1**.**

The pressure constraints (18,19) are satisfied by defining an alternative constraint on the difference of the maximum and minimum pressure realized at all consumption points,

[TABLE]

where $\tilde{p}^{h}$ denotes the pressure level at consumption point $h$ resulting from a simulation with source pressure $\tilde{u}_{p}$. Here, $\overline{\Delta p}$ is the true limit for the pressure difference entering the alternative optimization constraint (24). The control of the pressure level leading to admissible pressures at consumption points is obtained a-posteriori by a time dependent shift,

[TABLE]

Proof.

By limiting the difference between maximal and minimal absolute pressure levels in (24), their value relative to the pressure control is limited as well,

[TABLE]

where $\Delta p^{h}$ denotes the pressure difference from house to source, which is independent of the pressure control at the source by (23). Thus, the pressure level at each consumption point for an arbitrary source pressure $p^{s}$ reads,

[TABLE]

Inserting the suggested control (25) allows to determine the minimum pressure at each consumption point by

[TABLE]

Similarly, the maximum pressure level is limited by

[TABLE]

where the inequality in (29) is obtained by using (26). ∎

Feed-in power and control of energy density

The feed-in power is the central constraint to limit additional costs, since it avoids the usage of additional energy resources. It is defined by

[TABLE]

The control $u_{T}^{\kappa}$ affects the feed-in $P$ (30) in two ways. First, by setting the current input energy density $u_{T}^{\kappa}$ and second, by defining the volume flow which results from the current energy densities at heat exchangers. These in turn equal the control $u_{T}^{\kappa}(\tau)$ at a past time $\tau$. Depending on the current state $e_{h,n_{h}}$ at consumer stations, the input control can both amplify and weaken the feed-in power with regard to the current consumption $G$. In the stationary case $e=u_{0}$, where $u_{0}$ denotes the constant input, the feed-in power is the temporally shifted consumption profile. Hence, it also matches the high characteristic power peaks in the morning and the evening hours. In contrast, by anticipating the expected consumption and the transport time of the injected power, peaks in the injected power can be reduced. Since the determination of an optimal control is initialized with a constant temperature, the red, solid lines in parts (b) of fig. 2, 4, 6 visualize the consumption profile equaling the displayed feed-in.

5 Determination of an optimal control

Subsequently, we discuss the computation of an optimal control for the problem (15-21). The main idea is to eliminate the transport dynamics from the optimization constraints by solving them explicitly and passing the remaining constraints to the MATLAB nonlinear optimization tool fmincon.

Supplying an initial control which satisfies all constraints is an open problem. Hence, the initial parameter set $\kappa_{0}$ generally at least violates the feed-in constraint (20) and the determination of a feasible solution is performed initially. For the current parameter vector $\kappa_{i}$, the transport dynamics described by the DAE (21) are solved along $\mathcal{T}_{d}$. In the solution process, both the trajectory of state variables and their parameter gradients are calculated. This allows to evaluate the true optimization constraints (16-17, 20, 24), and their gradients with respect to the current parameter $\kappa_{i}$. In this step, the limits for the pressure levels (18, 19) are replaced by (24) as described in proposition 1. This allows to focus on the determination of the thermal control $u_{T}^{\kappa_{i}}$ in solving the optimal control problem. The required pumping power resulting from the pressure control $u_{p}$ can be neglected as described in section 4. Solving the dynamics explicitly avoids the large computational cost of passing them as optimization constraints. Values and gradients of the optimization constraints for the current parameter are passed to fmincon using the active set method together with the value and parameter of the objective function (15). A summary of the algorithm used to determine the optimal control is provided in alg. 1.

Extraction of parameter gradients

To estimate the effect of a change in the parameterized control on the relevant outputs of the heating networks, the sensitivities of both the objective function and the constraints with respect to the parameters have to be determined $\forall t\in\mathcal{T}_{d}$. To this end, gradients of both the control and state variables regarding the control parameters have to be extracted from the forward solution of the DAE (21). For the input signals typically applied to heating networks, an implicit time integration of the DAE proved to be beneficial. A general implicit time integration scheme, in which $x_{\mu}$ and $x_{\mu+1}$ denote the state variables at former and future time levels can be written as

[TABLE]

The derivative of the future state variable $\partial_{\kappa}x_{\mu+1}$ is obtained by the derivative of the old state variable $\partial_{\kappa}x_{\mu}$ using the implicit function theorem,

[TABLE]

Hence, based on the sensitivity of the initial state $\partial_{\kappa}e_{\alpha,\beta}(t=0)$, the gradient information can be propagated along the solution of the DAE. Using (33) allows to determine the gradient after solving for the new time step. This is in contrast to many automatic differentiation approaches in which the gradient information has to be tracked during the determination of the future time layer causing additional computational cost.

For the time integration of the DAE in this contribution, the implicit midpoint rule as a second order symplectic integrator is used,

[TABLE]

where dt denotes the time step, and $f$ the dynamical part of the DAE (21).

6 Numerical validation

6.1 Time integration of the forward simulation

The solution of the forward problem (21) within the determination of an optimal control is performed by the system descriptions (11, 14). To solve the DAEs (11, 14) within the time horizon required for the optimal control, the implicit midpoint-rule (34) is used. Full order models are unreduced ($W=V=\mathds{1}$), while for the reduced order model a Galerkin projection is applied. Sparse matrix operations are considered in the full order case. To solve systems (11,14) efficiently, a domain decomposition is performed [23]. Different parts of the network are treated as independent systems, with their linkages moderated by artificial inputs. This allows to form the system operators and solve the nonlinear equations introduced by the implicit time integration scheme efficiently. Full and reduced order models are simulated using the decomposition. Due to the affine system representation, the determination of the Jacobian can be determined analytically for both full- and reduced order models. The simulations presented in the following sections are performed using MATLAB(R) R2016b on an Intel(R) XEON(R) CPU E5-2670 processor @ 2.60GHz. The nonlinear system of equations resulting from both the algebraic equations and the time integration scheme are solved using the MATLAB function fsolve.

6.2 Definition of test scenarios

To demonstrate the effectiveness of the reduced order model, the large-scale heating network presented in fig. 1 is studied. The loops visible in the left part of the network pose a central difficulty since the thermal transport can take different paths at the same time to reach a certain destination in the network. Moreover, changes of flux direction occur, which change the set of possible paths the transported quantity takes dynamically.

The robustness of the reduced model towards its application in the optimization is evaluated for different environmental temperatures defining the consumption behavior for the large scale network. These scenarios cover the relevant mean daily temperatures $\{-3,3,7.5\}\,\,^{\circ}\text{C}$. The interval spanned by $[-3,7.5]\,\,^{\circ}\text{C}$ exhibits a large optimization potential in terms of distributing the feed-in power. For colder or warmer environmental temperatures, either all or none of the energy capacities within the power plant will be used. The lower and upper temperature constraints are in line with standard operation conditions of heating networks. Table 2 presents a detailed description of the optimization scenarios under investigation.

The observation interval, in which the constraints and the objective function are evaluated, is set to $\text{dt}=300\,$s, which is smaller than the typical plant operation interval of $900\,$s. This allows to approximate the underlying dynamics more precisely, while matching the relevant decision interval. The power constraint $\bar{P}\in[G_{0},\bar{G}]$ is chosen within the mean ($G_{0}$) and maximum ($\bar{G}$) daily consumption. While the mean consumption $G_{0}$ naturally poses a lower limit for the maximum injected power, the maximum consumption $\bar{G}$ is an upper limit, since it can always be achieved by a stationary control. Due to the initialization with a stationary solution, during the first period the power restriction is relaxed to the maximum consumption. During this time, the output energy density and the corresponding constraints are shaped by the initial solution and not the control. The friction factor $\lambda$ entering the Darcy-Weisbach law (2) is modeled by the Colebrook-White equation depending on the Reynolds number of the fluid, and the roughness and diameter of each pipeline. A computationally efficient approximation of the frictional model is achieved by setting a time-independent friction factor $\lambda_{i},\;i\in\mathcal{P}$ for every pipeline by a-priori defining a suitable Reynolds number. The difference of maximal and minimal pressure levels at all consumption points in the flow network is given by $\overline{\Delta p}=$2.5 bar and replaces the optimization constraints (19,18). For the gravitational constant we use a value of $g=9.81\,$m s*-2*. Specific heat capacity and density are set to $c_{p}=4.16\,$kJ K*-1* kg*-1*, and $\rho=1000\,$kg m*-3*. To transform energy densities to temperatures presented in the section, we use the approximation

[TABLE]

This allows to transform the constraints for the minimum energy density at consumer stations (17), as well as the maximum energy density in the network (16) to their corresponding temperature values $T_{\text{min}}^{\text{cons}}$, and $T^{\text{net}}_{\text{max}}$ displayed in tab. 2.

The power extraction $G$ is modeled using demand profiles typically employed in the simulation of heating networks [30]. Specifically, each heat exchanger exhibits the power demand

[TABLE]

The customer specific scaling factor $c_{i}$ represents an estimate for the total daily energy demand. The profile $s_{m}(t,T_{d})$ models the time dependent demand for a given consumption class $m$ and the daily mean environmental temperature $T_{d}$. Every member of the class thus shares the same normalized profile while exhibiting an own specific consumption. The available demand profiles adjust the daily consumption by hourly scaling factors. The latter also depend on $T_{d}$, adjusting the relative weights of each hour in the daily consumption. Based on these hourly values, a spline interpolating the consumption is generated. In the considered network, the largest part of consumers belongs to the same consumption class $s_{0}$. To this end, all consumers are modeled by this class.

Since the typical consumption follows a periodic profile if the daily mean environmental temperature does not change, the control $u_{T}^{\kappa}$ is parameterized by a Fourier series,

[TABLE]

which approximates any control $u\in L_{2}$. Here, the frequency $\omega=2\pi/\theta$ is fixed to the period length $\theta$ of the consumption signal corresponding to 24 h. The remaining Fourier coefficients act as parameters to be optimized.

6.3 Optimal control for TC1

The discussion of the numerical optimization results starts with TC1 simulating a mean daily temperature of $-3\,^{\circ}\text{C}$ with a mean and maximum consumption of $1.64\,$MW, and $2.29\,$MW.

Before comparing runtime and optimal controls obtained by different spatial discretizations, we analyze the optimal control suggested by ROM1. This model is gained by Galerkin projection of the upwind discretization FOM1, cf. tab. 3. Based on a constant initial control of $90\,^{\circ}\text{C}$, avoiding high feed-in peaks forces the control to increase its temporal variation, cf. fig. 2(a). The resulting feed-in power remains below the power constraint for all times cf. fig. 2(b). Since the mean temperature changes only slightly, the resulting total volume flow injected at the power plant also remains on the same level compared to the initial control, cf. fig. 2(a).

The suggested temporal variation of the thermal control which is necessary to limit the feed-in power induces pre-heating effects, visualized in fig. 3. Thus, reflecting the transport delay from source to the sink, in certain time intervals marked in red, the injected feed-in power exceeds the current consumption. More specifically, the maximum temperature level is injected at the power plant before the maximum aggregated consumption occurs at the consumption points. A common feature of the standardized consumption profiles are peak consumption around 6 a.m. and 6 p.m. Indeed, two pre-heating phases can be observed reflecting these two phases of high consumption, cf. fig. 3.

To avoid finite horizon effects, in which the suggested optimal control exploits the energy incorporated in the initial state, the setup is simulated for three periods of 24 hours. Focusing on the injected power coupling to volume flows at consumer stations as state variables, a state close to periodicity is reached quickly.

Optimal controls resulting from different spatial discretizations

The control resulting from ROM1 which is presented in fig. 2 attains the feed-in constraint at several points of time. To check for feasibility of the solution, the suggested optimal controls resulting from both full- and reduced order models are compared and validated using a reference discretization FOMR. The latter is given by an upwind discretization in space with a high number of finite volume cells. This will answer the question, whether a coarse, unreduced model is appropriate for an optimization task as well. To this end, table 3 compares the runtimes and approximation qualities of both full and reduced order models. While FOM0 denotes the minimal upwind discretization, in which each pipeline receives one finite volume cell, FOM1 is a finer discretization, with small approximation errors compared to FOMR. Finally, ROM1 is the reduced order model obtained by reducing the discretized model FOM1. As quality indicators of the optimal control, the objective function, the relative $l_{2}$ error of the control and the maximal relative error of the feed-in constraint are considered. In addition, the maximum relative $l_{2}$ error of all outputs $y$ is measured. Here, $y^{h}_{R}\in\mathcal{H}$ refers to output $h$ determined by the reference discretization. Specifically, $\tilde{u}$ is the optimal control determined by the fine model FOMR, and $\lVert{P-\bar{P}}\rVert_{\text{max}}/\bar{P}$ results from the control suggested by a coarse model simulated using FOMR.

Focusing on TC1 and the suggested objective functions $J(u)$, FOM1, FOMR, and ROM1 converge to a comparable value with FOMR taking the minimum of 52.6. FOM0 deviates clearly to 55.9. Regarding the relative deviation to the reference control, ROM1 shows the best approximation with $7.25\times 10^{-4}$, followed by the full order models FOM1 and FOM0. The relative violation of the feed-in constraint is smallest for ROM1, followed by FOM1 and FOM0, while all models exhibit relative errors below one percent. Regarding the approximation quality of the outputs, FOM0 shows the expected strong diffusion, leading to a maximum relative error of $3.03\times 10^{-2}$. Although this error does not affect the feasibility of the feed-in constraint, it leads to violations of the temperature constraints, measurable in practical applications. In contrast, ROM1 still approximates the outputs with an error of $5.28\times 10^{-3}$.

Focusing on the runtimes for determination of the optimal control, ROM1 allows for a speed-up of 4.5 of the entire optimization compared to the coarsest and thus fastest possible unreduced model FOM0. In addition, the speed up compared to FOM1 amounts to 4.3, while achieving comparable results. The runtime of FOM0 results from a higher number of iterations necessary to fulfill the optimization tolerances. For the determination of the optimal control, resimulating the dynamics for a new control candidate takes the largest computational cost.

Different environmental temperatures

Hereafter, test cases TC2, TC3 are discussed, simulating higher mean daily temperatures. In these scenarios, both mean and maximum power consumption decrease compared to TC1. To achieve comparable power constraints for different test cases, a relative constraint $\bar{P}=G_{0}+0.5(\bar{G}-G_{0})$ for the feed-in is chosen.

In contrast to TC1, the optimal control suggested for TC2 by ROM1 decreases the thermal control $u_{T}^{\kappa}$ compared to the initial control at $90\,^{\circ}\text{C}$. The corresponding average volume flow increases. This happens at the expense of increasing pumping costs which can be neglected as described in section 4. Since the total power consumption for TC2 is smaller, the resulting increase in the volume flow does not exceed the level of TC1, allowing for identical pressure differences within the required constraints. In addition to a decreased mean value of the thermal control, also the temporal variation can be decreases. This leads to a smaller objective function of 14.6. Concerning error measures, the same observations discussed for TC1 apply. The speed-up of ROM1 compared to FOM0(FOM1) results in a factor of 4.2(5.4).

For the last scenario TC3 simulating a mean daily temperature of $7.5\,^{\circ}\text{C}$, the thermal control decreases in average value and temporal variation even further. The reduced consumption allows to increase the injected volume without violating pressure constraints. The injected flow temperature does now approaches the lower limit defined at $75\,^{\circ}\text{C}$. As observed for the other test cases, the approximation of the feed-in power is remarkably precise even for large deviations in the approximations of the outputs. For the coarsest discretization FOM0, the relative error of the feed-in constraint results in $1.46\times 10^{-2}$. To illustrate this observation for TC3 , fig. 5(a) shows the output with the largest relative error compared to the reference discretization FOMR. The robustness of the feed-in constraint towards errors in the state-space approximation is unexpected by (30), in which the feed-in depends on the volume flows at households defined by the thermal outputs $y$. Two explanations can be supplied for this effect. First, the feed-in depends on the sum of volume flows over all consumer points, allowing approximation errors to cancel. Furthermore, as the upwind discretization is a conservative finite volume scheme, the total stored energy is preserved on every discretization level.

In contrast to the large deviations observed for FOM0, ROM1 displayed in fig. 5(b) exhibits smaller errors in the output approximation. Specifically, it mainly deviates around the discontinuity at $t=54$ h resulting from changes in the flux direction occurring in the dynamical simulation visible for the reference discretization. Again, the error in the feed-in constraint resulting from an imprecise approximation of flux changes is small by two reasons. First, only few consumers are affected by changing flux directions. Second, a change of flux directions is associated with the volume flow tending to zero, forcing the implied power to be zero as well. Hence, the absolute error in the power approximation remains small.

6.4 Speed-up of the reduced order model

As discussed above, resimulating the dynamics is the central computational cost in determining an optimal control. One cause for the speed-up of the ROM is the few number of simulations of the forward problem to satisfy the tolerances for constraints and the objective function. The second cause is the speed-up resulting for a single solution of the DAE (21), which is discussed subsequently. For the implicit midpoint rule with $n_{t}$ time steps, the computational cost $c_{sim}$ splits into the following parts,

[TABLE]

where $n_{imp}$ is the number of iterations to solve for the upcoming timestep, $c_{h}$ is the cost of solving the algebraic equations (6-8) defining the flow field, and $c_{f}$ is the cost of evaluating the differential part of the DAE. Finally, $c_{J}$ denotes the cost to determine the Jacobian of the DAE with respect to energy densities. Since only the thermal transport is reduced, $c_{h}$ is identical for both full and reduced order models. The cost for the evaluation of the ODE scales with the number of entries implied in the system operators $A(v),A_{r}(v)$ defined in (11,14). Using the upwind discretization, $A(v)\in\mathbb{R}^{n\times n}$ is a sparse matrix. For typical networks, the number of nonzero entries can be limited by $3n$. In contrast, the system operator $A_{r}(v)\in\mathbb{R}^{r\times r}$ resulting from a Galerkin projection is dense with $r^{2}$ nonzero entries. As a consequence, the number of reduced states needs to be significantly smaller to reduce the computational cost. For the coarse discretization sufficient for the determination of an optimal control, this degree of reduction is barely possible. The key saving in applying the reduced order model stems from the computation of the Jacobian matrix. Based on the system operator description and the definition of the DAE (11)

[TABLE]

the Jacobian reads

[TABLE]

Since by the algebraic equation (8) each velocity depends on the energy densities at households, $J_{\dot{e}}(v,e)$ carries significantly more non-zero entries in the full order case than $A(v)$, cf. tab. 6. In contrast, since the reduced order operator $A_{r}(v)$ is already densely populated by Galerkin projection, the number of non-zero entries does not increase for the reduced Jacobian.

7 Conclusions

In this contribution, we discussed the optimal control of district heating networks utilizing a reduced order model (ROM). The suggested optimal controls resulting from minimizing the temporal variation of the control successfully limit the maximum feed-in power to the average of mean and maximum total consumption. In addition, practically relevant constraints on temperature and pressure are included reproducing realistic operation conditions. For the presented scenarios, this allows to avoid the usage of additional, unfavorable sources of energy. The quadratic in-state DAE of index 1 is split into a thermal and a flow defining part, allowing to describe the thermal part as a linear, parameter-varying system with the velocity field acting as the parameter. Using a greedy strategy, relevant velocity configurations are implied in a global Galerkin projection forming a Lyapunov stable ROM. While the ROM approximates both relevant state variables and gradient information sufficiently fine for the determination of an optimal control, it speeds up the entire optimization phase by at minimum a factor 4, compared to even coarse levels of upwind discretizations used as full order models. For distinct test scenarios we observe even higher speed-ups of 7.3. Thus, the ROM gaps the bridge towards the determination of an optimal control within an online planning. The effectiveness of the ROM is demonstrated for an existing large scale network in which different pipelines change their flux direction dynamically. Runtime and approximation quality are studied for multiple real world scenarios including varying daily mean temperatures. This allows to apply the presented model to other networks and operation conditions relevant in practice.

For further research, we study the benefits of the reduced order model in a feedback control, in which additional advantages might result from its reduced state space dimension. Furthermore, a comparison to direct optimization approaches will be interesting, in which the transport dynamics of the network directly appear as optimization constraints, avoiding to resimulate the network dynamics in each iteration.

Acknowledgments

We acknowledge the financial support by the German Federal Ministry BMBF within the project EiFer (Förderkennzeichen: 05M18AMB), and by the German Federal Ministry BMWi within the project DYNEEF (Förderkennzeichen: 03ET1346B).