Reduction of a district heating model using network decomposition

In this contribution we study model order reduction of nonlinear transport networks for district heating systems. In these, heating energy injected at a centralized power plant is transported to consumers. They are modeled by a hyperbolic differential algebraic system with large state space dimension. The network structure introduces sparse system dynamics, which transform to a dense reduced system leading to unacceptable computational costs [1]. To exploit the benefits of sparsity, sub‐parts of the network are reduced separately in a structure preserving way using Galerkin projection [2]. After introducing the dynamical model, we discuss a decomposition strategy which aims at minimizing the number of observables for each subnetwork. We demonstrate its numerical benefits at an existing large scale heating network.


Parametric model for district heating and reduction by Galerkin projection
The dynamics of heating networks are described by one dimensional incompressible Euler equations. The latter incorporate the advection of the energy density ϕ by a purely time dependent volume flow q(t) in each pipeline, as well as the conservation of momentum. In this, the pressure difference along a pipeline is dominated by frictional forces, allowing to neglect the contribution of acceleration in a quasi-stationary limit. After spatial discretization of the advection equation using the upwind scheme, a state space decomposition is performed allowing to formulate the differential algebraic equation The algebraic constraints (3) result from the coupling of pipelines at nodes of the network and gather three contributions: the conservation of volume, Kirchhoff's second law claiming that pressure differences along network loops sum up to zero, and the volumeflow balance at consumers. It forces the product of volume flow q at consumer stations and observed energy density y to equal the power consumption u H . Hence, (3) collects quadratic constraints in the state variables ϕ, q. The system description is decomposed to a main network ϕ 0 to which the thermal control u T is applied, and c ∈ N subnetworks forming ϕ which receive observablesỹ as artificial inputs. Observables y are measured by both main-and subnetworks, whileỹ is the system state of the main network at attachment points to subnetworks.Ã,B,C represent continued block diagonal matrices. Their blocks A s , B s , C s , s ∈ {1, .., c} describe the dynamics of the subnetworks. Both A s (q), and B s (q) depend on the volume flow q, and can be represented in an affine decomposition We approximate the input-output characteristics of each system s ∈ {0, .., c} by moment matching of the volume flow dependent transfer function H(q, ω) = C(ω1 − A(q)) −1 B(q) in frequency domain. Considering the volume flow q as a timevarying parameter to the advection of the energy density, local Galerkin projections V s e (q e ) are performed for representative volume flow vectors q e . These are picked using a greedy strategy and ultimately combined to a global projection V s using a singular value decomposition for each s ∈ {0, .., c}. Application of V s to (1-3) defines the reduced order model (ROM).
Eliminatingỹ in (1), using (2), one obtains an undecomposed system description Σ equivalent to (1)(2)(3) where the network links are rewritten to the system matrix A(q). It can be shown that Σ and its reduced model Σ r obtained by a Galerkin projection V are Lyapunov stable for an upwind discretization [3]. Forming V by diagonal blocks matching the network decomposition allows to conclude for Lyapunov stability of the reduced, decomposed system Π r , which is seen as follows. 8.39 × 10 −3 10928 Table 1: Maximum relative error δ, runtime and estimated complexity for the reference network RND. The error δ = max i∈{1,..,n(y)} compares the PDE solution vector y to the approximation y r . Results are obtained using MATLAB ® (R2016b) on an Intel ® Xeon ® CPU E5-2670 @ 2.60GHz. DOF denotes the number of finite volume cells for both full and reduced models. The time integration for all models is performed by the MATLAB ® implicit Runge-Kutta scheme ode15s.
By eliminatingỹ, Π r is again equivalent to the reduction of the undecomposed system Σ r , which is known to be Lyapunov stable. Hence, in the formation of a block diagonal V , any Galerkin projection V s , s ∈ {0, .., c} can be applied to each of the subsystems, ensuring Lyapunov stability.

Decomposition algorithm and numerical validation
Towards a beneficial network decomposition, we aim at equalizing the number of observables n(y s ) to be approximated in each subnetwork. To this end, the following heuristic is proposed. In a directed graph G, a recursive tree search is performed. The number of allowed outputs in potential subnets is restricted to {2, .., N }, where N ∈ N is an external parameter. Based on the current root node r, the dimension n(y s ) of the output vector y s in the underlying subnet S s is counted. If n(y s ) ∈ {2, .., N }, and S s contains no outputs which are already clustered, S s is defined to be a new subnet. Since the optimal choice of N is not defined a-priori, the cost function J(ỹ) = n(ỹ) 2 + c s=0 n(y s ) 2 is evaluated to indicate beneficial cluster decompositions. J(ỹ) penalizes both the number of connection points n(ỹ) of main-and subnetworks, and the number of observables in each subnetwork n(y s ) quadratically. For the reference network street (RNS), the cost-minimal choice N = 8 resulting from the variaton of N is visualized, cf. figure 1.
The influence of the proposed decomposition on the simulation runtime is demonstrated for a computationally more demanding reference network district (RND) including n(y) = 333 outputs. Since it includes changes of flux directions which are complex to reduce, we encapsulate them in a single subnet which is not reduced. We then compare three simulation models. The full order model discretized by the upwind scheme with 1 subnet (FOM 1), the FOM with 29 subnets proposed by our algorithm (FOM 29) and the ROM obtained from reducing all 29 subnetworks except for the one including flux changes (ROM 29), cf. table 1. Comparing FOM 1 to FOM 29 shows that a large speed up already results for the unreduced system by network decomposition, by a more efficient operator assembly of A(q). An additional saving is accounted to the reduction by Galerkin projection.
We presented a heuristic for a network decomposition targeting to reduce the number of observables to be approximated in each subnet. We demonstrated that the resulting ROMs allow for a faster computation of the dynamics of heating networks. An additional speed up by the reduction of flux changing networks as well as different decomposition strategies [4] will be investigated in the future.