Optimal control of district heating networks

Standard models for district heating networks usually lead to a high number of coupled nonlinear and algebraic equations, which makes an optimal control of such networks a difficult task. We present an iterative algorithm, in which we separately resolve the high dimensional equality constraints, leaving us with a series of low dimensional quadratic programs with linear inequality constraints, which can be solved by established methods.


Mathematical Network Model And Problem Description
In district heating, thermal energy is generated at one (or multiple) centralized producers, and then transported to consumers and back to the production site via water pipelines. This leads to the fact, that only the pressure levels and the supply temperature can be controlled independently, whereas the mass flows result from the coupling with consumer stations.
We briefly outline our mathematical model for the district heating network. The hydraulic behaviour of water pipelines is modelled by the one dimensional, incompressible Euler equations, where pressure loss due to friction is described by model of the Colebrook and White. The transport of internal energy is described by an unsteady advection equation without diffusion, and a sink term describing energy loss due to wall cooling. Consumer stations and producers are simply described by algebraic constraints, describing power demand and temperature input, respectively. As coupling condition at junctions we choose conservation of mass-and energy flows, where the average junction temperature follows from the perfect mixing assumption. In addition, we formulate box constraints for the controls, as well as lower bounds for consumer supply temperature and pressure gradients.
Classically, such networks are operated in either a pressure or a temperature driven manner, where only one of these quantities is varied actively, whereas the other one is held constant. In the pressure driven case, a quasistationary model yields a sufficient approximation of the networks behaviour, and can be directly embedded into an optimal control problem [1]. If temperature variations are admitted, varying transport velocities and time delays become an essential part of the system's dynamics, making it necessary to consider the fully instationary model.
The problem we are going to examine throughout the remainder of this article slightly differs from usual optimal control tasks in district heating: Given a initial state x(t 0 ) with control u(t 0 ) = u ref , we want to solve the quadratic problem u = arg min where the additional constraint enforces an upper bound for the supply power, which depends on both, the actively controlled supply temperature, and the mass flow rate. Therefore, this approach is neither temperature nor pressure driven.

Numerical Optimization by Reduction to Quadratic-Linear Programs
In this section, we briefly describe an iterative algorithm which efficiently solves the discretized counterpart of the optimal control problem (1). Let ∈ R n·m denote the discretized state after the k th iteration, where each x (k) i ∈ R m represents the spatial discretization of the system's state at time t i = t 0 +i∆ t. We consider the subset of equations describing the transition between two time steps where (3) corresponds to the network's interior dynamics, (4) is the control input, and (5) are the state dependent inequality constraints, including the supply power threshold. The system (3)-(4) is solved using Newtons method, which yields x (k) i . Thereafter, we calculate the sensitivities up to the i th time step by the recursion formula which is a consequence of the implicit function theorem. Here, J i , which can be computed effectively using sparse automatic differentiation. With this information at hand, we provide a first order approximation of the inequality constraints w.r.t. the control inputs u (k) . Together with the discretized objective functional, this results in a quadratic program with linear constraints with respect to the search direction d (k) = u (k+1) − u (k) , which we solve using a dual method by Goldfarb et al. [2].

Numerical Results And Conclusion
We apply our method to a real world network with approximately 1500 pipelines, where the upper power limit has been set to 4.8 MW. The resulting control is shown in figure 1, where it is compared to a constant reference curve. By slightly increasing supply power ahead of consumption peaks, the increased amount of energy within the system, which makes it possible to stay below the threshold. At this point we should notice that the supply power constraint adds a lot of complexity to the optimization problem. In general it is not obvious, whether any given initial state admits a feasible (dynamic) control at all. To circumvent this issue, we need use homotopy methods to obtain a valid initial state, for which we can solve our optimal control problem.