An adjoint optimization approach for the topological design of large-scale district heating networks based on nonlinear models
AAn Adjoint Optimization Approach for the Topological Design ofLarge-Scale District Heating Networks based on Nonlinear Models
M. Blommaert ∗ , , , , Y. Wack , , , M. Baelmans , PUBLISHED IN: APPLIED ENERGY, Volume 280, 15 December 2020, 116025DOI:10.1016/j.apenergy.2020.116025 Department of Mechanical Engineering, KU Leuven, Celestijnenlaan 300 box 2421, 3001Leuven, Belgium EnergyVille, Thor Park, Poort Genk 8310, 3600 Genk, Belgium Flemish Institute for Technological Research (VITO), Boeretang 200, 2400 Mol, Belgium ∗ corresponding author, e-mail: [email protected] Abstract
This article deals with the problem of finding the best topology, pipe diameter choices, and operation parametersfor realistic district heating networks. Present design tools that employ non-linear flow and heat transport models fortopological design are limited to small heating networks with up to 20 potential consumers. We introduce an alterna-tive adjoint-based numerical optimization strategy to enable large-scale nonlinear thermal network optimization. Inorder to avoid a strong computational cost scaling with the network size, we aggregate consumer constraints with aconstraint aggregation strategy. Moreover, to align this continuous optimization strategy with the discrete nature oftopology optimization and pipe size choices, we present a numerical continuation strategy that gradually forces thedesign variables towards discrete design choices. As such, optimal network topology and pipe sizes are determinedsimultaneously. Finally, we demonstrate the scalability of the algorithm by designing a fictitious district heating net-work with 160 consumers. As a proof-of-concept, the network is optimized for minimal investment cost and pumpingpower, while keeping the heat supplied to the consumers within a thermal comfort range of . Starting from auniform distribution of
15 cm wide piping throughout the network, the novel algorithm finds a network lay-out thatreduces piping investment by and pump-related costs by a factor of 14 in less than an hour on a standard laptop.Moreover, the importance of embedding the non-linear transport model is clear from a temperature-induced variationin the consumer flow rates of . The decarbonization of heating supplies is a crucial step in the fight against climate change. With district heating systems, large shares of presently untapped waste heat could be harnessed for use in space heating [1]. By additionallyintegrating renewable heating sources—such as biomass, geothermal, and solar heat—district heating systems providea viable solution to reduce the global dependence of heating and cooling supplies on fossil fuels [2]. that aim at fully exploiting this district heating potential while minimizing grid losses shouldhereby be targeted [3]. The low network operation temperatures of 4th generation networks, however, put a lot ofpressure on the (re-)design of district heating networks.Recently, a number of researchers have explored the use of numerical optimization to obtain optimal network topologies , dimensioning, and operation parameters. The discrete nature of many design variables naturally leads tothe use of
Mixed-Integer Programming (MIP). In combination with the nonlinearity of the underlying district heatingnetwork physics this poses a challenging optimization problem. Given the difficult nature of such problems, manyauthors simplify the underlying model to obtain a Mixed-Integer Linear Program (MILP) for the optimal network.S¨oderman [4], for example, used a MILP approach for the optimal structure and configuration of a district cooling net-work including cold media storage. Later, Dorfner and Hamacher [5] used a similar approach to optimize the topologyand pipe sizes of a district heating network. The application of MILP was later expanded to simultaneously optimize an1 a r X i v : . [ c s . C E ] O c t ncreasing amount of parameters of these networks. Haikarainen et al. [6] simultaneously optimized the topology andoperations, considering different supplier technologies and heat storage facilities. Mazairac et al. considered topologyoptimization of multi-carrier networks including gas and electricity [7]. Morvay et al. [8] combined optimal design ofthe network with designing the supplied energy mix. In a study to identify an optimal set of consumers to connect to adistrict heating network, Bordin et al. [9] optimized the topology using a MILP approach. These tools can, however,not accurately capture some important intrinsically nonlinear features of 4th generation networks. Their linearizednature fails to accurately describe the flow pattern in networks with multiple producers or with loops, the heat lossesthroughout the networks, or the temperature-dependent consumer heat transfer effectiveness.More recent papers therefore shift towards at least partially capturing the nonlinearities in the district heatingnetwork models. Given the difficulty of solving a Mixed-Integer Nonlinear Program (MINLP), one popular approachis to employ heuristic optimization approaches. Tan et al. [10] used a whale optimization metaheuristic to optimizethe daily operation of a distributed energy system, while Ayele et al. [11] employed a particle swarm algorithm tooptimize temperature levels and the size of heat pumps within a district heating network. The most commonly usedheuristic approaches are evolutionary algorithms. Wang et al. [12] optimized the pressure head of circulation pumpswithin a multi-source district heating network using a genetic algorithm. Their study includes a detailed nonlinearhydraulic model of the pipe network, heat sources, and substations. Hirsch et al. [13] performed an optimal designstudy of a single heat connection using a genetic algorithm. The authors optimized the design of a long distance heattransportation system, accounting for a detailed thermal and hydraulic model of the connecting pipeline. In anotherapproach proposed by Vesterlund et al. [14], a set of producer parameters in a multi-source network is optimizedusing a hybrid Genetic Algorithm - MILP strategy. This two layer approach uses a full thermo-hydraulic networkmodel to optimize the producer parameters for minimal fuel costs in a nested MILP layer. None of the aforementionedheuristic approaches include the network topology in the optimization. Li and Svendsen [15], conversely, used geneticalgorithms for the optimization of network topology and pipe diameter dimensioning. In their work, they optimized anetwork of 10 consumers, while accounting for nonlinear thermo-hydraulic network modelling aspects. Discrete pipediameters were also obtained by simple rounding.Another strategy to cope with the nonlinearities in the network is to solve the MINLP using commercial MINLPsolvers. This approach is particularly popular for solving scheduling and operational planning problems of districtheating networks. Examples can be found in Deng et al. [16] who employ it for a nonlinear scheduling optimization,and Zheng et al. [17] who optimize the operational planning of different heat sources in an integrated energy system.Merkert et al. [18] optimized the optimal operation of district heating grids while exploring globalization of the unitcommitment problem. Only few authors apply MINLP solvers to network topology and pipe design problems. Marty etal. applied the DICOPT-GAMS MINLP solver to the simultaneous optimal design of the Rankine Cycle of a geothermalplant and a small heating network topology with up to 8 consumers [19], however, using a simplified network model.The same MINLP solver was also applied by Mertz et al. to the topological design and pipe sizing for an academicdistrict heating network study case [20] and later to the design of a small existing network with 19 consumers [21].While only Li and Svendsen [15], Marty et al. [19], and Mertz et al. [20] deal with the nonlinearities that followfrom the treatment of momentum and energy conservation in the topology optimization problem, the amount of designvariables that can be dealt with is strongly limited by their steep computational cost scaling. This is also evident fromthe relatively small design problems that are considered in these papers.The problem of computational cost scaling is further explored in a topology analysis by Allen et al. [22] andvon Rhein et al. [23]. The authors investigate the steep scaling of the topological search space of a 5th generationDistrict Heating and Cooling system for an exhaustive search. To circumvent this, they propose a minimal spanningtree heuristic to only explore a subset of the search space. The number of minimal spanning trees to explore foreach possible subset of connected buildings still becomes intractable for districts of increasing size [22]. The designexamples shown in these papers are limited to networks with four potential consumers. The argument of intractablecost scaling is further supported by the work of Weinand et al. [24]. The authors developed a heuristic to designdistrict heating system topologies, circumventing the steep cost scaling. Their comparative optimization approach withCPLEX was not able to exceed 7 settlements.In contrast, optimal design applications for inherently nonlinear problems in amongst others computational fluiddynamics (CFD)—such as shape optimization of airfoils [25], topology optimization for flow [26] and heat transfer[27] problems, and even the heat exhaust optimization in nuclear fusion reactors [28]—exploit efficient adjoint -basedoptimization techniques. These methods are continuous optimization methods based on adjoint gradient calculations,which exhibit a computational cost scaling independent of the number of design variables in the optimization problem.2lso in the design of flow networks, the use of these methods has been explored for topology optimization [29], andcombined shape and topological optimization [30]. And most recently, in the context of district heating networks,Pizzolato et al. [31] applied an adjoint-based procedure to the robust hydraulic optimization of the topology, notconsidering thermal aspects. Around the same time, the authors of this paper published a brief first study of an adjoint-based topology optimization approach for district heating networks, emphasizing the importance of combined thermal-hydraulic modelling [32]. Both papers relaxed the topology optimization problem to that of a continuous pipe sizingproblem. Moreover, neither of these papers consider thermo-economic optimal design of district heating networks.The contribution of this paper is to introduce a scalable adjoint-based approach for thermo-economical optimaldesign of district heating network topologies that simultaneously optimizes discrete pipe size and operational parame-ters. In contrast to simplified linearized approaches used in MILP optimization, it is first of all (1) based on a completetransport model to deal with the nonlinear effects that characterize 4th generation district heating networks. By beingbased on the adjoint approach, it second (2) offers good scalability towards large numbers of design variables. Theapproach therefore is suitable for large scale district heating network design. Third (3), it improves upon the methodol-ogy introduced by the authors of this paper [32] by handling the discrete nature of design choices in the district heatingnetwork topology and integrating additional design constraints using aggregation techniques. Since this is a proof ofconcept study, we solve an optimization problem that aims at minimizing the cost associated with the installed networkpiping, and installed pump capacity and its operation, while meeting the thermal demand of all consumers within areasonable tolerance. While this is not a full detailed thermo-economic cost function yet, it contains the CAPEX andOPEX contributions of a return-on-investment analysis that challenge the optimization approach most.Several difficulties have to be overcome within this paper before the adjoint method can be considered a validalternative to the above-mentioned MIP methods. First of all, the adjoint optimization method is continuous by nature.Since the design of topologies is discrete (to put or not to put a pipe somewhere), the method needs to be modified tohandle discrete design variables. A related challenge here is how to account for costs intrinsically related to topologicalchoices, such as the cost of the piping works. Secondly, although the adjoint method does not scale with the number ofdesign variables, it does scale directly with the number of so-called ‘state’ constraints. These are constraints on flow,pressure, or temperature variables that depend on the solution of a network simulation themselves. And lastly, sinceit is a gradient-based optimization method, it is likely that the optimization result is a local and not a global optimum.The question therefore arises whether a useful solution can be found at all.To deal with the state constraints on the thermal energy delivery, we propose the use of constraint aggregation [33],in order to maintain a computational cost scaling that is independent of the number of consumers in the problem. Inorder to avoid the need of relaxing the discrete optimization problem to a diameter sizing problem, a new method isintroduced inspired by the smoothed Heaviside projection method from structural topology optimization [34], alongwith a penalization approach to account for pipe investment costs. This allows for the pipes to be selected from adiscrete set of available sizes, next to the selection of the optimal topology.The paper is outlined as follows. In section 2, we present a model based on conservation of mass, momentum, andenergy, that governs the nonlinearities in the energy transport, and that includes a flow friction correlation valid forall relevant pipe flow regimes. Next, in section 3, we introduce the optimization problem and its solution approach,including the constraint aggregation and the new projection method. Finally, in section 4, we show the correct func-tioning of the procedure on a thermal network design problem of realistic size, demonstrating its scalability. Becauseof the novelty of the adjoint approach for the thermo-hydraulic optimization of district heating networks, the efficientimplementation of the adjoint sensitivity calculation is also elaborated in C. In addition, the model linearizations thatare needed for the adjoint approach can be exploited to build an efficient Newton solver for this highly nonlinearmodel. This approach is detailed in B. The notations that are necessary to describe both the state and adjoint solversare introduced in A. In this paper, we will start from a physics-based thermal network model, building on conservation of mass, momentum,and energy as summarized in a previous work by the authors of this paper [32]. In the next section the optimizationapproach will then be built upon this state-of-the-art model. The network will be designed for a worst-case situation,so that a steady-state model can be adopted. In this section, we will give an overview of the complete set of modelequations for each network component. For clarity, we start this section by introducing the notations needed to representthe thermal network as a mathematical graph . 3igure 1: An illustration of a directed graph representation of a district heating network. Producer, consumer, andinternal nodes are denoted by hexagons, squares, and circles, respectively. Dark full arrows represent feed arcs, whilelight full arrows represent return arcs. Dotted arrows represent producer (red) and consumer arcs (blue).
Referring to Figure 1, the thermal network—including pipes, junctions, producers, and consumers—can be representedby a directed graph G ( N, A ) , with N the set of all nodes and A the set of all arcs in the graph. Physically, thegraph consists of a feed and return network, with arcs of opposite directions. This graph, which includes all potentialconnections, is often referred to as the superstructure . We can further subdivide the set of nodes N into three subsets N p ∪ N c ∪ N i = N , denoting the producer, consumer, and internal nodes, respectively. Similarly, the set of arcs A is partitioned into the subsets A p ∪ A c ∪ A i = A , denoting the producer, consumer, and internal arcs, respectively.The internal arcs hereby represent a geometrical pipe connection, while the other arcs rather represent state transitionsat the producer and consumer ends, but at the same physical location. In what follows, we can now denote a certainnetwork node as n ∈ N and we can denote a directed arc going from node i to node j as ( i, j ) ∈ A , or succinctly as ij ∈ A or even more compactly as a ∈ A . We will consider 4 different categories of thermal network components that need modelling, as indicated with dashedboundary lines in Figure 1: pipes, pipe junctions, consumers, and producers. We outline their models here one by one.
For pipe flow, the momentum equation relates the gradient of the static pressure p to the viscous friction. We considerwater as the carrier fluid, with a constant density ρ , dynamic viscosity µ , and specific heat capacity c p , thus neglectingtheir temperature-dependence. In the example elaborated in section 4, the values will be taken as corresponding to awater temperature of ◦ C. We then use the empirical Darcy-Weisbach equation to model the viscous pressure dropfor incompressible flow as a function of the volumetric flow rate q ij through a pipe ( i, j ) with length L ij , ( p i − p j ) g ij = q ij , ∀ ij ∈ A i . (1)The flow conductance g ij in this expression is based on the correlation of Cheng [35] covering the entire range of flowregimes from laminar to turbulent, which we modified to behave regularly towards the limit of zero-sized diameters. Itis given by g = | q | ( f − − f ) ( d − (cid:15) ) f Θ f (1 − f )L Θ f (1 − f )(1 − f )TS Θ − f )(1 − f )(1 − f )TR Θ C , ith Θ L = ρ dπµ , f = 11 + (cid:16) ReRe LT (cid:17) , Θ TS = 1 . (cid:18) Re . (cid:19) , f = 11 + (cid:18) Re ( d (cid:15) ) (cid:19) , Θ C = d π ρL , f = 11 + (cid:0) d − (cid:15) (cid:15) (cid:1) , Θ TR = 2 log (cid:18) . d(cid:15) (cid:19) . (2)Here, (cid:15) denotes the absolute pipe roughness, d the inner diameter of the pipe and Re the Reynolds number, definedas Re = 4 ρ | q | /( πµd ) . The transition from laminar to the turbulent regime is controlled by Re LT = 2720 .Next, we consider the thermal aspects of an insulated pipe installed underground. First, we introduce θ = T − T ∞ as the temperature difference with the outside air temperature T ∞ . Let us consider θ i to be this temperature differenceat the node i , at which the flow enters the pipe ij , and θ ij at the pipe exit. The pipe exit temperature difference θ ij , dueto heat loss to the environment is given by θ ij = θ i exp (cid:18) − L ij ρc p | q ij | R t,ij (cid:19) , ∀ ij ∈ A i , (3)with R t,ij the thermal resistance per unit pipe length between the energy carrier and the environment. For a pipe withan outer isolation casing diameter d o,ij that is assumed to be bigger than the inner diameter d ij by a fixed ratio, i.e. d o,ij = rd ij , the combined thermal resistance of pipe and soil per unit length is [36] R t,ij = ln(4 h/ ( rd ij ))2 πλ g + ln r πλ i , (4)with λ i and λ g the thermal conductivity of the insulation and the surrounding ground, respectively, and h the depth atwhich the pipe is located. The above relations for the pipe temperature drop are validated with experimental data inRef. [37]. All nodes in the network represent pipe junctions. The flow in these junctions is governed by conservation of mass.For the incompressible flow under consideration, this reduces to conservation of the flow rate q a through the arcs a connected to the junction, i.e. (cid:88) a =( i,n ) ∈ A q a − (cid:88) a =( n,j ) ∈ A q a = 0 , ∀ n ∈ N. (5)Similarly, conservation of the convected energy will give us relations for the temperatures in the nodes. In the junction,perfect mixing of the incoming flows is supposed to obtain outgoing flows at the node temperature. Note that dependingon the sign of the flow rate q in the directed arcs connected to the node, the flow will either enter or leave the junction.Energy conservation can thus be formulated as (cid:88) a =( i,n ) ∈ A (max( q a , θ a + min( q a , θ n ) − (cid:88) a =( n,j ) ∈ A (max( q a , θ n + min( q a , θ a ) = 0 , ∀ n ∈ N. (6) We now introduce a basic model to estimate the heat transferred to the consumer. We suppose the consumer substationconsists of a bypass and a heating system directly connected to the network (see Figure 2). This leads to two subsets of5igure 2: Illustration of the consumer modelconsumer arcs, namely the set of bypass arcs A cb ∈ A c and the set of heating system arcs A ch ∈ A c . Both arcs have acontrol valve to regulate the flow. The pressure drop over the consumer heating arc is assumed to be of the form p i − p j = ζ | q ij | q ij α ij , ∀ ij ∈ A ch , (7)with ζ a constant determined from nominal operating conditions [38], and α ∈ [0 , the valve control. It is importantto note that the valve control is not trying to model the physical control of the valve, but rather attempts to obtain goodconditioning for the optimization problem, while respecting the lower limit of the pressure drop ( α ij = 1 ). Similarly,the bypass flow is controlled through the relation q ij = β ij q max,b ∆ p des,b ( p i − p j ) , ∀ ij ∈ A cb , (8)with β ∈ [0 , the bypass valve control, q max,b a maximal bypass flow and ∆ p des,b a typical pressure drop over thebypass.Conservation of energy in the heating system is given by ρc p q ij ( θ i − θ ij ) = Q ij , ∀ ij ∈ A ch , (9)with Q ij the heat transferred to the house through the heating system. The latter is modelled in the form of thecharacteristic equation for radiators [39, 40]: Q ij = ξ ij θ i − θ ij ln ( θ i − θ house )( θ i,j − θ house ) n ij , (10)with θ house the known temperature difference between the indoor and environment temperature at the house. Values ofthe coefficients ξ ij and n ij are tabulated for individual radiators, according to the EN 442-2 norm [41]. The bypassarcs on the other hand are assumed to be free of heat losses, i.e. θ ij = θ i ∀ ij ∈ A cb (11) In the producer arcs, a fixed input flow q b is imposed as boundary condition for this system of equations. In addition,a given temperature is imposed for the heat source. This leads to q a = q b,a , θ a = θ b,a ∀ a ∈ A p . (12)To uniquely define the pressures throughout the network, a reference pressure is imposed in one of the producer returnnodes. Since only pressure differences influence the flow solution, the exact choice of reference does not impact thesolution. 6 .3 Solution procedure Together with this reference pressure, equations (1-12) jointly lead to a well-defined system that, for given design andoperational parameters, can be solved for the flow rates q ij , pressures p i , temperatures θ i and θ ij , and consumer heatfluxes Q ij throughout the network. These variables are called the state variables of the model equations and can begrouped into a single state vector x . The system is solved by an efficient Newton procedure to obtain second orderconvergence. The analytical linearization of the model equations that is required for this also forms the basis for theadjoint solver (see section 3). The structure of the model equations can be exploited to decouple the solution of flowrates and pressures at the one hand and temperatures at the other hand. We detail the procedure to efficiently solve thisstrongly nonlinear system of equations in B. In order to describe the solver, we need to reformulate the equations ingraph notation as an algebraic system of equations. This is done in A. For a given superstructure, we can now use the model presented in section 2 to find an optimal configuration fora thermal network by optimizing the pipe diameters d , along with the network operation parameters, including theproducer inflows q b , and the consumer controls α and β . We will jointly denote all these design variables with thedesign variable vector ϕ .To formulate the design problem as a mathematical optimization problem, the first crucial step is the choice ofthe objective/cost function J ( ϕ , x ) that is to be minimized. Ideally, one would optimize the network design for themaximal return on investment or maximal energy efficiency. However, the aim of this paper is rather to demonstrate thevalidity of the presented approach. As such, we deliberately choose a cost function that captures only two of the maincontributions, being the cost associated with the installed pump capacity and its operation, and that of the investmentcost related to piping. We thus obtain the following cost function formulation that will be employed for design in aworst-case scenario: J ( ϕ , x ) = λ P (cid:88) ij ∈ A i ( p i − p j ) q ij + (cid:88) ij ∈ A i C ( d ij ) L ij , (13)with λ P the price of increased pump capacity and operation, and C ( d ) a function representing the investment cost ofpiping placement per unit length as a function of the pipe diameter. Because of the nonlinear dependence of both termson the pipe diameters, this choice of cost function illustrates well the advantages of the adjoint optimization strategy.Note that especially the first term is difficult to accurately represent using a linearized approximation.Simultaneously, we impose that in this worst-case situation the heat supplied to the consumers remains within a range of the desired heat supply Q d ,a , i.e. . Q d ,a ≤ Q a ≤ . Q d ,a , ∀ a ∈ A ch (14)Because these constraints depend on the state variables x , they will be referred to as state constraints . They will besuccinctly denoted further using a single generic state constraint vector h ( ϕ , x ) ≤ . Finally, minimal and maximalvalues for all design variables are typically imposed, leading to box constraints in the form ϕ min ≤ ϕ ≤ ϕ max . Theseconstraints are relatively easy to deal with and will be further on only implicitly denoted by an admissible set Φ ad ofdesign values, giving ϕ ∈ Φ ad . We can then write our design problem in the generic optimization problem formulation min ϕ ∈ Φ ad , x J ( ϕ , x ) cost function s.t. c ( ϕ , x ) = 0 , model equations h ( ϕ , x ) ≤ , state constraints (15)with c ( ϕ , x ) = 0 succinctly denoting the model equations (1-12), as detailed in A.7ince the set of model equations c ( ϕ , x ) = 0 we defined in section 2 can be solved for a state solution x ( ϕ ) foreach evaluated network design ϕ , we can use them to eliminate the state variables x in the optimization problem. Theso-called reduced optimization problem then reads min ϕ ∈ Φ ad ˆ J ( ϕ ) s.t. ˆ h ( ϕ ) ≤ , state constraints , (16)with ˆ J and ˆ h the reduced cost function and reduced state constraint vector, respectively. The adjoint optimizationstrategy discussed in the remainder of this section gives us an efficient approach to solve this optimization problem. One possible strategy to solve the nonlinear optimization problem (16) is by using a
Sequential Quadratic Programming (SQP) strategy. One hereby linearises the constrained optimization problem in optimization iteration l to a linearly-constrained quadratic program for the step size p l = ϕ l +1 − ϕ l , min p l ∈ Φ p = { p l | ϕ l + p l ∈ Φ ad } ∇ ˆ J ( ϕ l ) (cid:124) p l + 12 p (cid:124) l B k p l s.t. (cid:98) h ( ϕ l ) + ∇ (cid:98) h ( ϕ l ) (cid:124) p l ≤ , (17)where B k is an approximation of the reduced hessian ∇ ˆ J , which will be estimated from gradient information insubsequent optimization iterations with a damped BFGS algorithm [42]. It is clear that the main challenge is thus theefficient calculation of the cost function and state constraint sensitivities ∇ ˆ J and ∇ ˆ h . For this purpose, the adjointapproach is used.Because of the novelty of the adjoint approach in its application to district heating network optimization, we willelaborate next on how to evaluate these sensitivities using an adjoint approach. For compactness of notation, we willleave out the operator arguments below if they are clear from the context. In addition, we define the vector function I = [ J , h (cid:124) ] (cid:124) that groups the quantities of interest of which we need the sensitivities to solve the quadratic program(17) for an optimization step p l . Its sensitivities are then given by the sensitivity matrix ∇ ˆ I := d ˆ I d ϕ (cid:124) .Using the chain rule of differentiation, we can decompose this matrix into two contributions, i.e. d ˆ I d ϕ = ∂ I ∂ ϕ + ∂ I ∂ x d x d ϕ . (18)It should be noted that while the first term is a quite simple term that only contains the direct dependence of the costfunction on our design variables (zero in our case), the second term features the nonlinear dependence d x d ϕ of the modelequation solution. An expression for d x d ϕ can be found by linearising the model equations c ( ϕ , x ) = 0 around thedesign variables ϕ l and model solution ¯ x l = x ( ϕ l ) of the current iterate, giving ∂ c ∂ x d x d ϕ = − ∂ c ∂ ϕ . (19)It can be observed that the number of linear algebraic sets of equations that need to be solved to evaluate the completesensitivity matrix d x d ϕ scales directly with the number of design variables in ϕ . The adjoint method avoids this in anelegant way. It can be easily derived by substituting (19) in (18). After rearranging terms we get d ˆ I d ϕ = ∂ I ∂ ϕ + (cid:34) − (cid:18) ∂ c ∂ x (cid:19) − (cid:124) (cid:18) ∂ I ∂ x (cid:19) (cid:124) (cid:35) (cid:124) ∂ c ∂ ϕ . The term between brackets can be obtained by solving the adjoint equation (cid:18) ∂ c ∂ x (cid:19) (cid:124) x ∗ = − (cid:18) ∂ I ∂ x (cid:19) (cid:124) , (20)8or the so-called adjoint variables ¯ x ∗ l = x ∗ ( ϕ l , ¯ x l ) . The derivative matrix in optimization iteration l then evaluates as ∇ ˆ I ( ϕ l ) = (cid:18) ∂ I ∂ ϕ ( ϕ l , ¯ x l ) (cid:19) (cid:124) + (cid:18) ∂ c ∂ ϕ ( ϕ l , ¯ x l ) (cid:19) (cid:124) ¯ x ∗ l . (21)This latter expression is cheap to evaluate, since ∂ c ∂ ϕ can be derived analytically and evaluated. The main compu-tational cost is therefore in the solution of the adjoint equation (20). It can be observed that the number of algebraicsystem solutions needed to solve the adjoint equation is now independent of the number of design variable entries in ϕ , but rather scales with the number of quantities of interest in I . For an unconstrained optimization problem with ascalar-valued cost function, the gradient calculation cost is therefore reduced to a single system solution. In the currentoptimization problem formulation, an additional system solution is needed for each state constraint.We therefore choose to limit the number of state constraints drastically by aggregating all these constraints using avariant of the Kreisselmeier-Steinhauser function [33], turning these constraints into a single scalar-valued constraint ˆ h ks ( ϕ , γ ) ≤ , with ˆ h ks ( ϕ , γ ) = 1 γ ln n h (cid:88) i =1 e γ [ˆ h ] i . (22)In this equation, [ˆ h ] i represents a specific component of the vector ˆ h with length n h , and γ ∈ ]0 , ∞ [ is a parameterthat controls the trade-off between smoothness and numerical conditioning of the approximation on the one hand,and exactness of the approximation on the other. In practice, numerical continuation is used to increase it’s value γ in different stages of the optimization as explained in section 3.3. Using this aggregated constraint in combinationwith the adjoint method we effectively reduce the computational costs of the complete gradient matrix evaluationto approximately two algebraic system solutions, independent of the number of design variables. Thus, a scalableoptimization approach is achieved. The algorithm presented thus far accounts only for continuous optimization variables, such as the pipe diameters.This could, in theory, lead to a whole range of diameters, with some of them perhaps of negligible size. In order toperform a true topology optimization we need to include the choice whether or not to put a pipe in a particular arcof the superstructure. We therefore modify the diameter design variable to either not place a pipe at all ( d → ),or to choose from a number of standard pipe sizes. We then end up with a discrete design variable d ∈ D , with D a discrete set of diameters, e.g. D = { , . , . , . } m. To achieve this discrete nature without givingup the benefits of the adjoint approach, we introduce a novel numerical continuation strategy that gradually forcesthe continuous optimization method towards these discrete diameter choices. The proposed strategy consists of twoelements: the smoothed projection of the diameters onto the discrete diameter set and the gradual penalization ofintermediate diameter values through the piping cost relation C ( d ) . We will now elaborate these components. For the projection of the diameters on the discrete diameter set D , we modify the so-called projection method that isused in the field of structural topology optimization to obtain binary variables with a continuous optimization method[34]. This technique uses a smooth approximation of a Heaviside function with described by [43] ¯ ϕ = P ( ϕ, η, χ ) = tanh ( χη ) + tanh ( χ ( ϕ − η ))tanh ( χη ) + tanh ( χ (1 − η )) (23)to project the design variable ϕ on a binary decision variable ¯ ϕ ∈ { , } that determines whether or not materialshould be placed at a certain location in a structure. The variables χ ∈ ]0 , ∞ [ and η ∈ [0 , are two factors controlling,respectively, the steepness of the Heaviside approximation and the threshold value above which the variable ϕ isprojected onto to upper limit. Note that we in practice substitute the discrete diameter “0” in this set with a small but non-zero diameter at which the model equations are stillwell-defined and their solution numerically well-conditioned. ϕ i onto a binary design variable ¯ ϕ i . (b)Illustration of the extended smoothed Heaviside projection of ϕ i onto a variable ¯ ϕ i onto some discrete set, chosen hereto be a set of pipe diameters D = { , . , . , . } m.We extend this principle here to allow projection on the multiple discrete variables in D . To this end, we constructan extended projection operation ¯ d = P ext ( d, χ ) that is a piecewise function assembled from rescaled versions ofthe projection operation P , as shown in Figure 3. The extended projection then uses the basic projection function P ( d, . , χ ) to project inbetween two discrete values, as well as project downwards (i.e. P ( d, , χ ) ) for values biggerthan the largest value in D . For a diameter d and discrete design variable set D = { , d , d , ..., d n } , the extendedprojection operator then becomes ¯ d = P ext ( d, χ ) = d P (cid:16) dd , . , χ (cid:17) d ≤ d d j + (cid:0) d j +1 − d j (cid:1) P (cid:16) d − d j d j +1 − d j , . , χ (cid:17) d j ≤ d ≤ d j +1 d n + w P (cid:16) d − d n w , , χ (cid:17) d n ≤ d , (24)with w a parameter set such that d < d n + w for all values of d and where we assume that d is always positive. In the beginning of this section we introduced C ( d ) as a function that represents the piping cost as a functionof the pipe diameter d . It should be noted, however, that the only relevant values of this function are those in thediscrete set C = { , c , c , . . . , c n } corresponding to the discrete diameters d ∈ D , since a cost can only be providedfor these available pipe sizes. For the values in between, a continuous interpolation is to be provided to account forthe gradient-based optimization approach. We choose this interpolation as to deliberately make ‘intermediate’ values( d (cid:54)∈ D ) unfavorable for the optimization. Following this reasoning, we smoothly approximate a piecewise constantfunction that features a strong increase in price when the diameter exceeds a value in the set D . Such a function can bedescribed by the piecewise function C ( d, υ, ω ) = υc P (cid:16) dd , , ω (cid:17) + (1 − υ ) c n ( d − d ) d n − d d ≤ d υ (cid:16) c j + (cid:0) c j +1 − c j (cid:1) P (cid:16) d − d j d j +1 − d j , , ω (cid:17)(cid:17) + (1 − υ ) c n ( d − d ) d n − d d j ≤ d ≤ d j +1 υ (cid:16) c n + c w P (cid:16) d − d n w , , ω (cid:17)(cid:17) + (1 − υ ) c n ( d − d ) d n − d d n ≤ d . (25)Here, c w is an arbitrary high cost value to penalize oversized diameters . The variable υ ∈ [0 , controls theinterpolation between the piecewise linear- and the piecewise constant function. ω ∈ ]0 , ∞ [ controls the steepness of Oversized diameters are prevented with an appropriate choice of the box constraints, represented by the admissible set Φ ad . Figure 4: Evolution of the penalization on the pipe investment cost C . The optimization starts on a linear function(blue) that gets pushed towards a piecewise linear function (orange) by increasing ω . Next the cost is pushed towardsa Heaviside step function (green) by υ , representing the real cost of piping. The grey lines illustrate intermediate stepsin the continuation.the Heaviside approximation for the piecewise constant function. The process of interpolation and steepness increasecan be seen in Figure 4. Given the ill-posedness that the steep derivatives of this function induce on the optimizationproblem, a numerical continuation strategy is used to relax this problem in the initial stages of optimization. This willbe explained next. After including the aggregated constraint, the diameter projection, and the continuous piping cost function C ( ϕ i , υ, ω ) ,the optimization problem formulation reads min ϕ ∈ Φ ad ˆ J ( ¯ ϕ , υ, ω ) s.t. ˆ h ks ( ϕ , γ ) ≤ , aggregated constraint , ¯ ϕ = P ext ( ϕ , χ ) ≤ , extended projection , (26)in which the continuation variables υ , ω , γ , and χ control the trade-off between relaxing the discrete problem that is tobe solved to increase the well-posedness of the optimization problem, and the exactness of its approximation. We willtherefore apply numerical continuation to gradually force the design towards the discrete solution we are interestedin. That is, we add an outer loop to the optimization in which the above optimization problem is solved for differentvalues of the continuation variables, and initialize the optimization procedure with the solution of the last continuationiteration to overcome the problem of ill-conditioning.Note that the cost function and state constraints are evaluated in (26) based on the projected design variable. Usingthe chain rule of differentiation, we find that the derivative of the quantities of interest is then given by d ˆ I d ϕ = d ˆ I d ¯ ϕ ∂ P ext ∂ ϕ , (27) It is understood that the extended projection operator ¯ d = P ext ( d, χ ) can be generalized to operate on the complete vector of design variables,i.e. ¯ ϕ = P ext ( ϕ , χ ) , in such a way that only the discrete diameter variable components of ϕ are affected. T b , =70 ° C (left) and T b , = 65 ° C (right). Consumers are represented by black circles of varying sizes corresponding totheir respective heat demand ( ,
15 kW ,
50 kW ). b) Total length of all placed discrete pipe diameters in the feednetwork.where d ˆ I d ¯ ϕ can be efficiently calculated by the adjoint gradient calculation (21) as before and multiplied with theprojection gradient. In this section, we will demonstrate the functioning and scalability of the network optimization algorithm on a fictivetest case. Yet, we choose a somewhat realistic test case that aims at designing a medium-sized district in Waterschei,Genk. We take a case set-up with two producers at different inflow temperature levels ( T b , = 65 ° C , T b , = 70 ° C ).160 consumers of 3 different consumer types are distributed throughout the neighbourhood. A summary of the differentconsumer characteristics can be found in Table 1. As can be seen in Figure 5, the pipe superstructure for the networkoptimization is placed onto a part of the neighborhood’s street grid. In this test case, the total amount of design variablesis 632. To the best of the authors’ knowledge, it is the largest reported application of topology optimization for districtheating networks that is based on a fully-fledged flow and heat transport model.Table 1: Consumer characteristics Type Q d , a T in , d ∆T d n d ξ [ kWK n d ] ζ [10
12 Pa sm ] Dwelling
15 kW 55 ° C 20 ° C 1 . .
34 1 . Renovated dwelling ° C 20 ° C 1 .
42 0 .
06 0 . Commercial demand
50 kW 55 ° C 20 ° C 1 . .
23 13 . Assuming a worst-case situation with an outdoor temperature T ∞ = − ° C , we now try to design the network tohave the least possible investment cost, as well as pump capacity and operation cost. The weight for pump-relatedcosts was chosen to be λ P = 10 M e / kW in the current cost function definition, in order to limit the pressure dropsin the piping to reasonable engineering limits. The design is initialized by distributing a uniform pipe network with adiameter of d = 0 .
15 m over the entire superstructure (see Figure 5). Initially the control valves (heating system andbypass valves, as well as producer inflows) are fully opened ( α = 1 , β = 1 , q b = q b , max ). To start the optimization12ithin the heat satisfaction boundaries set by (14), we first optimize the network operation parameters for minimalconsumer heat dissatisfaction, i.e. J ( ϕ , x ) = ( Q a − Q d , a ) , (28)while fixing the pipe diameters. To choose reasonable consumer heating characteristics, we suppose the consumerheating systems are designed to provide this heat load at nominal/design conditions, described next. First of all, weassume the design is based on a desired temperature drop ∆ T d and pressure drop ∆ p d of
50 kPa over the heatingsystem. Next, we take the heating system inflow temperature T in,d and a room temperature of T room , d = 20 ° C . Wecan then obtain the heating system characteristics from ξ ≈ Q d T in,d − T out,d ln ( T in,d − θ house )( T out,d − θ house ) − n ij with T out,d = T in,d − ∆ T d ,ζ = ∆ p d q d = ∆ p d ∗ ( ρc p ∆ T d ) Q , (29)using the coefficient n found in Table 1 as a typical value for heating systems (see Ref. [39]).We suppose the producers deliver heat to the network at a constant temperature of T b , = 65 ° C and T b , = 70 ° C .Consequently, we can directly relate the controlled volumetric inflow q b at the producers to the heat inflow via Q in = ρc p q b ( T b − T out,d ) , (30)with T out,d = 40 ° C the design outflow temperature of the network. We can thus directly control the maximal inputheat delivered by the producers with a simple constraint on the inflow q b . Given a total heat demand of .
77 MW forthe 160 consumers, we can only satisfy the demand if the combined maximal heat provided by both producers is bigenough to cover both this demand and the heat losses on the feed line.After the aforementioned warm-starting approach, we apply the optimization procedure described in section 3 tooptimize the layout and pipe dimensions of this network, along with the producer inflow and consumer valve settings.A continuation strategy with 20 iterations was used for updating the continuation variables in the aggregated constraint,the smoothed multi-projection and the investment cost penalization, increasing them within the ranges γ ∈ [5 · , ] , χ ∈ [0 , , υ ∈ [0 , and ω ∈ [0 , . The maximal heat input Q in,max,1 = Q in,max,2 = 2 MW of each individualproducer was chosen high enough to supply the entire network. We consider a set of available pipe diameter choices D and their corresponding cost for road works and piping C as tabulated in Table 2.Table 2: Discrete diameters D available and their respective investment cost CD [m] C [ e / m]
46 min . The optimization was run on a single CPU of a standard laptop with an Intel Core i7 processor (1.9 GHz). Adecrease in the piping investment cost by
23 % from . e to . e with respect to the uniform pipe distributiondesign of Figure 5 was achieved (i.e. the design obtained after exclusively optimizing the network operation variables).It is understood that this uniform design features .
15 m pipes throughout the superstructure and that the optimizedoperation parameters satisfy the consumer heat demand up to a very high tolerance. The pump operation cost wasreduced by a factor of 14 with respect to this ‘uniform’ pipe design. The resulting network topology is given inFigure 6a. It can be directly observed that larger pipe diameters are provided for high demand consumers, due totheir greater flow needs. Moreover, the resulting pipe diameters all have discrete values thanks to the smoothed multi-projection strategy presented in this paper. The total length of pipe installed for each of the different pipe size optionscan be seen in Figure 6b. It should be noted that the topology optimization algorithm explicitly decided not to installany piping in about of the superstructure.The heat demands of all consumers are met with a producer inflow Q in , = 0 .
92 MW and Q in , = 1 .
09 MW .Minimizing the pump operational cost (13) the algorithm is able to push the supplied consumer heat Q a close to thelower limit of the thermal comfort range of the consumer (14). It is intuitive that the cost of the network designwill indeed be lowest if these constraints are only marginally satisfied. Indeed, through the increase of the constraint13igure 6: a) Optimal network topology and configuration. The green triangles represent heat producers of T b , = 70 ° C (left) and T b , = 65 ° C (right). Consumers are represented by black circles of varying sizes corresponding to theirrespective heat demand ( ,
15 kW ,
50 kW ). Placed pipes are drawn as red lines. Removed pipes are drawn in grey.The resulting topology shows a clear separation of the two neighborhoods. b) Total length of all placed discrete pipediameters in the feed network.smoothness parameter γ in the numerical continuation strategy, the algorithm achieves a maximum deviation of theconsumer heat loads from the lower limit of .
008 % while conservatively satisfying (14).It can be seen in Figure 7a that the network is split into a low temperature- and a high temperature network. Thetemperature differences throughout the network illustrate the importance of basing the network design on a fully-fledged non-linear model for flow and energy modelling. In order to satisfy the heat demand of consumers in thecolder parts of the network, the flow rate through the consumers heating system has to be increased (as can be seenin Figure 7b). The increase of flow rate in this test case reaches up to
72 % when comparing the flow rate for thehighest inflow temperatures of T in = 69 . ° C to the lowest inflow temperatures of T in = 56 . ° C , for the consumertype ‘dwelling’. It should be noted that this increase in flow rate influences the cost balance between investment andpump operation cost, ultimately driving the pipe design towards bigger diameters. We presented a methodology for the automated design of district heating network topologies on the basis of a nonlinearmodel of flow and heat transport. Because the computational cost of the adjoint sensitivity calculations are intrinsicallyindependent of the number of considered design variables, the methodology is expected to scale well to large heatingnetworks.This has been showcased with the design of a fictitious district heating network with a scale well beyond anyother applications of nonlinear topological design tools. On a fictitious design problem of a medium-sized district inWaterschei, the method shows to correctly retrieve a discrete network design that is able to reduce the piping investmentby 23% and the cost related to installed pump capacity and operation by a factor of 14 compared to the initial uniformpipe distribution. While further extensions of the cost function towards the direct optimization of the net present valuecould still improve the relevancy of the resulting designs, they can easily be incorporated in the current methodologywithout impacting the computational cost.The computational cost in this proof-of-concept design remained well below an hour on a standard laptop. The14igure 7: a) Average pipe temperature throughout the network. A clear separation in a low temperature network anda high temperature network can be observed. b) Temperature dependency of the flow rate required by the consumersheating system to satisfy the heat demand. In this figure only ‘dwellings’ and ‘renovated dwellings’ are considered.aggregation of additional design constraints is an important step to achieve this. In the design problem this approachwas tested for a constraint that enforces the heat delivered at each single consumer to deviate less than from its heatdemand for worst-case conditions. The constraint aggregation method successfully reduced them to a single constraint,while the delivered consumer heat is found to deviate less than .
008 % from the solution that is expected if imposedin a non-aggregated way.The added value of considering a full nonlinear model in the design clearly follows from the observation that theconsumer flows vary up to
72 % , only because of temperature variations throughout the network. These temperature-induced flow variations in turn influence the optimal pipe design and operational parameters, and are thus important tocapture already in the design phase. Linearized models, in contrast, fail to describe the flow pattern in networks withmultiple producers or with loops, as well as the heat losses throughout the network and the temperature-dependentconsumer heat transfer effectiveness. Nonlinear model-based design is therefore expected to play an increasingly-important role in the design of next-generation networks.Now that the validity of the method has been established, an important next step is the detailed benchmark of thenovel method and its scaling to that of other mixed-integer non-linear programming methods. It is important to notethat gradient-based optimization methods like the one presented in this paper could be more sensitive to find localoptima than other mixed-integer non-linear programming methods. It should therefore be thoroughly investigated fora number of practical design cases how good the resulting designs perform with respect to those resulting from otherdesign approaches. Finally, it should be mentioned that several uncertainties play a crucial role in the design of thermalnetworks, such as producer unavailabilites, possible future network extensions, or potential pipe failures. To be of trulypractical value, optimal design methods for district heating networks should eventually account for these uncertainties.
In this paper, we propose a scalable adjoint-based optimization strategy to determine the optimal network topology,discrete pipe diameters, and operational parameters on the basis of a fully-fledged non-linear flow and energy trans-port model. Building on recently introduced adjoint approaches for thermal network design with continuous designvariables, we introduce a numerical continuation strategy that combines a smoothed projection of the diameters with a15enalization of intermediate diameter values to gradually converge the optimization algorithm towards discrete diam-eter values. As such, an adjoint-based optimal design of network topology and discrete pipe diameters is enabled forthe first time.The scalability of the method is demonstrated with the design of the –to the best of the authors’ knowledge– largesttopological network design application reported thus far on the basis of a non-linear transport model. This is enabledhere through the intrinsic design-independent scaling of the adjoint method if combined with the suggested constraintaggregation approach. Moreover, by building the optimization on the basis of a nonlinear physics model, heat lossesand the temperature-dependent effectiveness of consumer heating systems can be dealt with already in the design phase.The strong variations in local temperatures and flow rates that are found in the test case, confirm the potential impactof these nonlinearities on the optimal design of next-generation thermal networks.Future research should aim at setting up a thorough benchmark with standard mixed-integer nonlinear programmingtools and developing a methodology to account for uncertainties in the operation conditions.
Acknowledgements
Maarten Blommaert is a postdoctoral research fellow of the Research Foundation - Flanders (FWO) and the FlemishInstitute for Technological Research (VITO), and Yannick Wack is a doctoral research fellow of VITO.
A Thermal network equations in vector form
To describe the discrete implementation of the solver for both model equations and adjoint equations, we first need toreformulate the equations (1-12) in vector form. Consider a graph with n nodes and m directed arcs. Let us start bydefining the following vectors for flow rates, pressures, temperatures and consumer heat fluxes: p = p i , i ∈ N, q = q ij , ij ∈ A, θ = (cid:32) θ n θ a (cid:33) , with θ n = θ i , i ∈ N and θ a = θ ij , ij ∈ A, Q = Q ij , ij ∈ A ch . (31)Using the same subscripts as used before for the subsets, we introduce vectors related to only one of the subsets A i , A c , A ch , A cb , A p of the arcs A in a similar way, e.g. q i = q ij , ij ∈ A i denotes the flow in the internal arcs. Next,we introduce the node-pipe incidence matrix A = ( a , . . . , a m ) ∈ R n × m that contains information on which nodes are connected by a directed arc. That is, each column a i contains “1” at thestart node i of arc ij , and a “-1” at the end node j . The continuity condition in the pipe junctions (5) is then easilytranslated to vector form as A q = 0 . The pressure drop over all arcs on the other hand is given by ∆ p = A (cid:124) p . Alsohere we can easily extend the definition to incidence matrices that include only a number of rows related to a specificsubset of arcs, e.g. A i = ( A ) i,j , i ∈ N and j ∈ A i .Using this notation, we can formulate the state equations related to the hydraulic flow succinctly in vector formas H ( ϕ , y ) = 0 , with y = ( pq ) the hydraulic state vector and ϕ the vector of design variables (defined further insection 3). That is, the vector H ( ϕ , y ) = A q A i (cid:124) p − R i ( ϕ , q ) q i A c (cid:124) p − R c ( ϕ , q ) q c = 0 , (32)then contains the residuals of node continuity, pipe momentum, and consumer arc momentum balances, respectively.The diagonal matrices R i ( ϕ , q ) and R c ( ϕ , q ) hereby represent the nonlinear flow resistance for internal and consumerarcs, respectively, as can easily be extracted from relations (1), (7), and (8). Note that we leave out the trivial producer16onditions and reference pressure condition for simplicity of notation, although they are in principle an integral part ofthe vector equation H ( ϕ , y ) = 0 .To obtain a similar vector equation for energy conservation, we introduce the arc inflow and outflow matrices A in = max( A diag ( sgn ( q )) , A out = min( A diag ( sgn ( q )) , , (33)with “sgn” the elementwise sign function and “diag” an operator that maps some argument vector d with elements d i to a diagonal matrix D with diagonal elements ( D ) i,i = d i . Energy conservation in the pipe junctions, pipes, consumerheating arcs, and bypass arcs can then be succinctly written as a vector equation E ( ϕ , y , θ ) = diag ( A in | q | ) θ n + A out diag ( | q | ) θ a A i , in (cid:124) F ( ϕ , y ) θ n − θ a ρc p diag ( q ch ) ( A ch , in (cid:124) θ n − θ a,ch ) − ξ (cid:16) A ch, in (cid:124) θ n + θ a,ch − θ house (cid:17) ◦ n A cb , in (cid:124) θ n − θ a,cb = 0 , (34)with F ( ϕ , y ) a diagonal matrix that has the fractional reduction of θ through the pipes on its diagonal, as can bededuced from (3).The hydraulic transport equation H ( ϕ , y ) = 0 and energy transport equation E ( ϕ , y , θ ) = 0 together form the state equations of the thermal network that can be solved for the hydraulic state y and thermal state θ , as will be furtherelaborated in the next subsection. In an even more compact form, we can thus represent the thermal network model asthe vector equation c ( ϕ , x ) = (cid:32) H ( ϕ , y ) E ( ϕ , y , θ ) (cid:33) = 0 , (35)with x = (cid:0) yθ (cid:1) = (cid:16) pqθ (cid:17) its state variables . B Implementation of the state solver
To efficiently deal with the nonlinear character of the presented model equations, we employ a Newton scheme toiteratively solve the model equations (35) for an update of the state variables. That is, in each iteration “ k ”, we obtainan update ∆ x k = (cid:16) ∆ y k ∆ θ k (cid:17) from ∂ H ∂ y ( ϕ , y k ) 0 ∂ E ∂ y ( ϕ , y k , θ k ) ∂ E ∂ θ ( ϕ , y k , θ k ) (cid:20) ∆ y k ∆ θ k (cid:21) = − (cid:20) H ( ϕ , y k ) E ( ϕ , y k , θ k ) (cid:21) (36)Note that while the solution of the hydraulic flow equations H ( ϕ , y ) = 0 does not depend on the temperatures θ ,the energy advection in the thermal equations E ( ϕ , y , θ ) does give rise to a flow dependency. We can therefore solve c ( ϕ , x ) = 0 for its solution ¯ x = x ( ϕ ) in two steps:1. Solve the hydraulic equations first iteratively for the hydraulic state ¯ y using updates ∆ y k from ∂ H ∂ y ( ϕ , y k )∆ y k = − H ( ϕ , y k ) . (37)2. Then plug this solution into the energy equations and solve for the thermal state ¯ θ using ∂ E ∂ y ( ϕ , ¯ y , θ k )∆ θ k = − E ( ϕ , ¯ y , θ k ) . (38) The operator ( a ) ◦ b denotes the Hadamard exponential , which simply represents the pointwise exponent of a vector a to the power b . ∂ H ∂ y and ∂ E ∂ θ . These deriva-tives can be obtained by analytical linearization of (32) and (34) with respect to the state variables. Though this is asomewhat cumbersome work, the approach induces fast second-order convergence. Moreover, the approach is partic-ularly interesting in combination with the adjoint approach, as the model linearization is needed again for the adjointsolver. C Implementation of the adjoint solver
The crucial step in the adjoint gradient calculation is the solution of the adjoint equation (20). Let us therefore elaborateon the efficient implementation of the adjoint solver, given the thermo-hydraulic model at hand. Rewriting the adjointequation in terms of H and E reveals a similar block-triangular structure as in the state equations: (cid:18) ∂ H ∂ y (cid:19) (cid:124) (cid:18) ∂ E ∂ y (cid:19) (cid:124) (cid:18) ∂ E ∂ θ (cid:19) (cid:124) (cid:20) y ∗ θ ∗ (cid:21) = − (cid:18) ∂ I ∂ y (cid:19) (cid:124) (cid:18) ∂ I ∂ θ (cid:19) (cid:124) , (39)with x ∗ = (cid:16) y ∗ θ ∗ (cid:17) .Block-Gauss elimination therefore leads to the following two-step approach to solve for the adjoint variables y ∗ and θ ∗ :1. Solve the adjoint thermal equations first iteratively for the adjoint thermal state θ ∗ from (cid:18) ∂ E ∂ θ (cid:19) (cid:124) θ ∗ = − (cid:18) ∂ I ∂ θ (cid:19) (cid:124) . (40)2. Then plug this solution into the adjoint hydraulic equations and solve for the adjoint hydraulic state y ∗ using (cid:18) ∂ H ∂ y (cid:19) (cid:124) y ∗ = − (cid:18) ∂ I ∂ y (cid:19) (cid:124) − (cid:18) ∂ E ∂ y (cid:19) (cid:124) θ ∗ . (41)The adjoint equations thus show a similar sequential structure as the state equations, though in reversed order. In-formation on the cost function dependence can be as such thought of as propagating backwards through the modelequations. Since these equations follow from linearization around ( ϕ , ¯ x ) , the matrices on the left-hand side of theadjoint equations are obtained by simple transposition of the matrices from the last iteration of the model solver ((37)and (38)). Because of the linear nature of the adjoint equations, a complete gradient evaluation can thus be performedin only a fraction of the time needed for a model evaluation, independent of the number of design variables. References [1] U. Persson, B. M¨oller, and S. Werner, “Heat roadmap europe: Identifying strategic heat synergy regions,”
EnergyPolicy , vol. 74, pp. 663–681, Nov. 2014.[2] S. Werner, “International review of district heating and cooling,”
Energy , vol. 137, pp. 617–631, Oct. 2017.[3] H. Lund, S. Werner, R. Wiltshire, S. Svendsen, J. E. Thorsen, F. Hvelplund, and B. V. Mathiesen, “4th GenerationDistrict Heating (4GDH),”
Energy , vol. 68, pp. 1–11, Apr. 2014.[4] J. S¨oderman, “Optimisation of structure and operation of district cooling networks in urban regions,”
Appliedthermal engineering , vol. 27, no. 16, pp. 2665–2676, 2007.[5] J. Dorfner and T. Hamacher, “Large-Scale District Heating Network Optimization,”
IEEE Transactions on SmartGrid , vol. 5, pp. 1884–1891, July 2014. 186] C. Haikarainen, F. Pettersson, and H. Sax´en, “A model for structural and operational optimization of distributedenergy systems,”
Applied Thermal Engineering , vol. 70, no. 1, pp. 211–218, 2014.[7] W. Mazairac, R. Salenbien, and B. de Vries, “Towards an optimal topology for hybrid energy networks,” in
Proceedings of the 22nd EG-ICE International Workshop, 13-16 July 2015, Eindhoven, The Netherlands , pp. 1–10, 2015.[8] B. Morvaj, R. Evins, and J. Carmeliet, “Optimising urban energy systems: Simultaneous system sizing, operationand district heating network layout,”
Energy , vol. 116, pp. 619–636, 2016.[9] C. Bordin, A. Gordini, and D. Vigo, “An optimization approach for district heating strategic network design,”
European Journal of Operational Research , vol. 252, no. 1, pp. 296–307, 2016.[10] Y. Tan, X. Wang, and Y. Zheng, “Modeling and daily operation optimization of a distributed energy systemconsidering economic and energy aspects,”
International Journal of Energy Research , vol. 42, no. 11, pp. 3477–3495, 2018.[11] G. T. Ayele, M. Mabrouk, P. Haurant, B. Laumert, B. Lacarri`ere, and M. Santarelli, “Exergy analysis and thermo-economic optimization of a district heating network with solar-photovoltaic and heat pumps,” in , pp. 1947–1959, Institute of Thermal Technology, 2019.[12] H. Wang, H. Wang, H. Zhou, and T. Zhu, “Modeling and optimization for hydraulic performance design in multi-source district heating with fluctuating renewables,”
Energy Conversion and Management , vol. 156, pp. 113–129,2018.[13] P. Hirsch, M. Grochowski, and K. Duzinkiewicz, “Decision support system for design of long distance heattransportation system,”
Energy and Buildings , vol. 173, pp. 378–388, 2018.[14] M. Vesterlund, A. Toffolo, and J. Dahl, “Optimization of multi-source complex district heating network, a casestudy,”
Energy , vol. 126, pp. 53–63, 2017.[15] H. Li and S. Svendsen, “District heating network design and configuration optimization with genetic algorithm,”
Journal of Sustainable Development of Energy, Water and Environment Systems , vol. 1, no. 4, pp. 291–303, 2013.[16] N. Deng, R. Cai, Y. Gao, Z. Zhou, G. He, D. Liu, and A. Zhang, “A minlp model of optimal scheduling for adistrict heating and cooling system: A case study of an energy station in tianjin,”
Energy , vol. 141, pp. 1750–1763,2017.[17] X. Zheng, G. Wu, Y. Qiu, X. Zhan, N. Shah, N. Li, and Y. Zhao, “A minlp multi-objective optimization modelfor operational planning of a case study cchp system in urban china,”
Applied Energy , vol. 210, pp. 1126–1140,2018.[18] L. Merkert, K. Listmann, and S. Hohmann, “Optimization of thermo-hydraulic systems using multiparametricdelay modeling,”
Energy , vol. 189, p. 116125, 2019.[19] F. Marty, S. Serra, S. Sochard, and J.-M. Reneaume, “Simultaneous optimization of the district heating networktopology and the organic rankine cycle sizing of a geothermal plant,”
Energy , vol. 159, pp. 1060–1074, 2018.[20] T. Mertz, S. Serra, A. Henon, and J.-M. Reneaume, “A minlp optimization of the configuration and the design ofa district heating network: Academic study cases,”
Energy , vol. 117, pp. 450–464, 2016.[21] T. Mertz, S. Serra, A. Henon, and J. M. Reneaume, “A minlp optimization of the configuration and the design ofa district heating network: study case on an existing site,”
Energy Procedia , vol. 116, pp. 236–248, June 2017.[22] A. Allen, G. Henze, K. Baker, and G. Pavlak, “Evaluation of low-exergy heating and cooling systems and topologyoptimization for deep energy savings at the urban district level,”
Energy Conversion and Management , vol. 222,p. 113106, 2020. 1923] J. von Rhein, G. P. Henze, N. Long, and Y. Fu, “Development of a topology analysis tool for fifth-generationdistrict heating and cooling networks,”
Energy conversion and management , vol. 196, pp. 705–716, 2019.[24] J. M. Weinand, M. Kleinebrahm, R. McKenna, K. Mainzer, and W. Fichtner, “Developing a combinatorial optimi-sation approach to design district heating networks based on deep geothermal energy,”
Applied energy , vol. 251,p. 113367, 2019.[25] A. Jameson, L. Martinelli, and N. Pierce, “Optimum Aerodynamic Design Using the Navier-Stokes Equations,”
Theoretical and Computational Fluid Dynamics , vol. 10, no. 1-4, pp. 213–237, 1998.[26] T. Borrvall and J. Petersson, “Topology optimization of fluids in stokes flow,”
International journal for numericalmethods in fluids , vol. 41, no. 1, pp. 77–107, 2003.[27] T. E. Bruns, “Topology optimization of convection-dominated, steady-state heat transfer problems,”
InternationalJournal of Heat and Mass Transfer , vol. 50, pp. 2859–2873, July 2007.[28] M. Baelmans, M. Blommaert, W. Dekeyser, and T. Van Oevelen, “Achievements and challenges in automatedparameter, shape and topology optimization for divertor design,”
Nuclear Fusion , vol. 57, no. 3, p. 036022, 2017.[29] A. Klarbring, J. Petersson, B. Torstenfelt, and M. Karlsson, “Topology optimization of flow networks,”
ComputerMethods in Applied Mechanics and Engineering , vol. 192, no. 35, pp. 3909–3932, 2003.[30] A. Evgrafov, “Simultaneous optimization of topology and geometry of flow networks,”
Structural and Multidis-ciplinary Optimization , vol. 32, no. 2, pp. 99–109, 2006.[31] A. Pizzolato, A. Sciacovelli, and V. Verda, “Topology optimization of robust district heating networks,”
Journalof Energy Resources Technology , vol. 140, no. 2, p. 020905, 2018.[32] M. Blommaert, R. Salenbien, and M. Baelmans, “An adjoint approach to thermal network topology optimization,”in
Proceedings of the 16th International Heat Transfer Conference, 10-15 August 2018, Beijing, China , 2018.[33] G. J. Kennedy and J. E. Hicken, “Improved constraint-aggregation methods,”
Computer Methods in AppliedMechanics and Engineering , vol. 289, pp. 332–354, June 2015.[34] J. K. Guest, J. H. Pr´evost, and T. Belytschko, “Achieving minimum length scale in topology optimization usingnodal design variables and projection functions,”
Int. J. Numer. Meth. Engng. , vol. 61, pp. 238–254, May 2004.[35] N.-S. Cheng, “Formulas for friction factor in transitional regimes,”
Journal of Hydraulic Engineering , vol. 134,no. 9, pp. 1357–1362, 2008.[36] D. D’Eustachio, “Criteria for thermal insulation for use on underground piping,” ASTM International, 1957.[37] B. van der Heijde, M. Fuchs, C. Ribas Tugores, G. Schweiger, K. Sartor, D. Basciotti, D. M¨uller, C. Nytsch-Geusen, M. Wetter, and L. Helsen, “Dynamic equation-based thermo-hydraulic pipe model for district heatingand cooling systems,”
Energy Conversion and Management , vol. 151, pp. 158–169, Nov. 2017.[38] M. Pirouti, A. Bagdanavicius, J. Ekanayake, J. Wu, and N. Jenkins, “Energy consumption and economic analysesof a district heating network,”
Energy , vol. 57, pp. 149–159, Aug. 2013.[39] , ch. HVAC systems and equipment. Atlanta, GA: American Society of Heating, Re-frigeration and Air Conditioning Engineers, 2012.[40] J. Grote, Karl-Heinrich; Feldhusen,
DUBBEL: Taschenbuch f¨ur den Maschinenbau . Springer-Verlag, 2014.[41] “EN 442-2:2014; radiators and convectors – Part2: Test methods and rating,” standard, European Committee forStandardazation, Brussels, Belgium, Oct. 2014.[42] J. Nocedal and S. Wright,
Numerical Optimization . Springer, 2006.[43] F. Wang, B. S. Lazarov, and O. Sigmund, “On projection methods, convergence and robust formulations intopology optimization,”