An Adjoint Optimization Approach for the Topological Design of

Large-Scale District Heating Networks based on Nonlinear Models

M. Blommaert∗,1,2,3 , Y. Wack1,2,3 , M. Baelmans1,2
PUBLISHED IN: APPLIED ENERGY, Volume 280, 15 December 2020, 116025
DOI:10.1016/j.apenergy.2020.116025
1
Department of Mechanical Engineering, KU Leuven, Celestijnenlaan 300 box 2421, 3001
arXiv:2008.08328v2 [cs.CE] 16 Oct 2020

Leuven, Belgium
2
EnergyVille, Thor Park, Poort Genk 8310, 3600 Genk, Belgium
3
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 parameters
for realistic district heating networks. Present design tools that employ non-linear flow and heat transport models for
topological 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. In
order to avoid a strong computational cost scaling with the network size, we aggregate consumer constraints with a
constraint aggregation strategy. Moreover, to align this continuous optimization strategy with the discrete nature of
topology optimization and pipe size choices, we present a numerical continuation strategy that gradually forces the
design variables towards discrete design choices. As such, optimal network topology and pipe sizes are determined
simultaneously. 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 pumping
power, while keeping the heat supplied to the consumers within a thermal comfort range of 5%. Starting from a
uniform distribution of 15 cm wide piping throughout the network, the novel algorithm finds a network lay-out that
reduces piping investment by 23% 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 variation
in the consumer flow rates of 72%.

1 Introduction
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 additionally
integrating renewable heating sources—such as biomass, geothermal, and solar heat—district heating systems provide
a viable solution to reduce the global dependence of heating and cooling supplies on fossil fuels [2]. 4th generation
district heating networks that aim at fully exploiting this district heating potential while minimizing grid losses should
hereby be targeted [3]. The low network operation temperatures of 4th generation networks, however, put a lot of
pressure 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 to
the use of Mixed-Integer Programming (MIP). In combination with the nonlinearity of the underlying district heating
network physics this poses a challenging optimization problem. Given the difficult nature of such problems, many
authors simplify the underlying model to obtain a Mixed-Integer Linear Program (MILP) for the optimal network.
Söderman [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 topology
and pipe sizes of a district heating network. The application of MILP was later expanded to simultaneously optimize an

1
increasing amount of parameters of these networks. Haikarainen et al. [6] simultaneously optimized the topology and
operations, considering different supplier technologies and heat storage facilities. Mazairac et al. considered topology
optimization of multi-carrier networks including gas and electricity [7]. Morvay et al. [8] combined optimal design of
the network with designing the supplied energy mix. In a study to identify an optimal set of consumers to connect to a
district 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 linearized
nature fails to accurately describe the flow pattern in networks with multiple producers or with loops, the heat losses
throughout 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 heating
network models. Given the difficulty of solving a Mixed-Integer Nonlinear Program (MINLP), one popular approach
is to employ heuristic optimization approaches. Tan et al. [10] used a whale optimization metaheuristic to optimize
the daily operation of a distributed energy system, while Ayele et al. [11] employed a particle swarm algorithm to
optimize temperature levels and the size of heat pumps within a district heating network. The most commonly used
heuristic approaches are evolutionary algorithms. Wang et al. [12] optimized the pressure head of circulation pumps
within a multi-source district heating network using a genetic algorithm. Their study includes a detailed nonlinear
hydraulic model of the pipe network, heat sources, and substations. Hirsch et al. [13] performed an optimal design
study of a single heat connection using a genetic algorithm. The authors optimized the design of a long distance heat
transportation system, accounting for a detailed thermal and hydraulic model of the connecting pipeline. In another
approach proposed by Vesterlund et al. [14], a set of producer parameters in a multi-source network is optimized
using a hybrid Genetic Algorithm - MILP strategy. This two layer approach uses a full thermo-hydraulic network
model to optimize the producer parameters for minimal fuel costs in a nested MILP layer. None of the aforementioned
heuristic approaches include the network topology in the optimization. Li and Svendsen [15], conversely, used genetic
algorithms for the optimization of network topology and pipe diameter dimensioning. In their work, they optimized a
network of 10 consumers, while accounting for nonlinear thermo-hydraulic network modelling aspects. Discrete pipe
diameters were also obtained by simple rounding.
Another strategy to cope with the nonlinearities in the network is to solve the MINLP using commercial MINLP
solvers. This approach is particularly popular for solving scheduling and operational planning problems of district
heating 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 unit
commitment problem. Only few authors apply MINLP solvers to network topology and pipe design problems. Marty et
al. applied the DICOPT-GAMS MINLP solver to the simultaneous optimal design of the Rankine Cycle of a geothermal
plant 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 academic
district 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 follow
from the treatment of momentum and energy conservation in the topology optimization problem, the amount of design
variables that can be dealt with is strongly limited by their steep computational cost scaling. This is also evident from
the 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] and
von Rhein et al. [23]. The authors investigate the steep scaling of the topological search space of a 5th generation
District Heating and Cooling system for an exhaustive search. To circumvent this, they propose a minimal spanning
tree heuristic to only explore a subset of the search space. The number of minimal spanning trees to explore for
each possible subset of connected buildings still becomes intractable for districts of increasing size [22]. The design
examples shown in these papers are limited to networks with four potential consumers. The argument of intractable
cost scaling is further supported by the work of Weinand et al. [24]. The authors developed a heuristic to design
district heating system topologies, circumventing the steep cost scaling. Their comparative optimization approach with
CPLEX was not able to exceed 7 settlements.
In contrast, optimal design applications for inherently nonlinear problems in amongst others computational fluid
dynamics (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-based
optimization 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.

2
Also in the design of flow networks, the use of these methods has been explored for topology optimization [29], and
combined 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, not
considering 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 sizing
problem. 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 optimal
design 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 complete
transport model to deal with the nonlinear effects that characterize 4th generation district heating networks. By being
based on the adjoint approach, it second (2) offers good scalability towards large numbers of design variables. The
approach 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 heating
network topology and integrating additional design constraints using aggregation techniques. Since this is a proof of
concept study, we solve an optimization problem that aims at minimizing the cost associated with the installed network
piping, and installed pump capacity and its operation, while meeting the thermal demand of all consumers within a
reasonable tolerance. While this is not a full detailed thermo-economic cost function yet, it contains the CAPEX and
OPEX 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 valid
alternative 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 to
handle discrete design variables. A related challenge here is how to account for costs intrinsically related to topological
choices, such as the cost of the piping works. Secondly, although the adjoint method does not scale with the number of
design 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, since
it 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. In
order to avoid the need of relaxing the discrete optimization problem to a diameter sizing problem, a new method is
introduced inspired by the smoothed Heaviside projection method from structural topology optimization [34], along
with a penalization approach to account for pipe investment costs. This allows for the pipes to be selected from a
discrete 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, and
energy, that governs the nonlinearities in the energy transport, and that includes a flow friction correlation valid for
all 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. Because
of the novelty of the adjoint approach for the thermo-hydraulic optimization of district heating networks, the efficient
implementation of the adjoint sensitivity calculation is also elaborated in C. In addition, the model linearizations that
are needed for the adjoint approach can be exploited to build an efficient Newton solver for this highly nonlinear
model. This approach is detailed in B. The notations that are necessary to describe both the state and adjoint solvers
are introduced in A.

2 A transport-based thermal network model

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 optimization
approach 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 model
equations for each network component. For clarity, we start this section by introducing the notations needed to represent
the thermal network as a mathematical graph.

3
Figure 1: An illustration of a directed graph representation of a district heating network. Producer, consumer, and
internal nodes are denoted by hexagons, squares, and circles, respectively. Dark full arrows represent feed arcs, while
light full arrows represent return arcs. Dotted arrows represent producer (red) and consumer arcs (blue).

2.1 Graph representation

Referring to Figure 1, the thermal network—including pipes, junctions, producers, and consumers—can be represented
by a directed graph G(N, A), with N the set of all nodes and A the set of all arcs in the graph. Physically, the
graph consists of a feed and return network, with arcs of opposite directions. This graph, which includes all potential
connections, is often referred to as the superstructure. We can further subdivide the set of nodes N into three subsets
Np ∪ Nc ∪ Ni = N , denoting the producer, consumer, and internal nodes, respectively. Similarly, the set of arcs A
is partitioned into the subsets Ap ∪ Ac ∪ Ai = 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 transitions
at the producer and consumer ends, but at the same physical location. In what follows, we can now denote a certain
network 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.

2.2 Network Component models

We will consider 4 different categories of thermal network components that need modelling, as indicated with dashed
boundary lines in Figure 1: pipes, pipe junctions, consumers, and producers. We outline their models here one by one.

2.2.1 Pipe model

For pipe flow, the momentum equation relates the gradient of the static pressure p to the viscous friction. We consider
water as the carrier fluid, with a constant density ρ, dynamic viscosity µ, and specific heat capacity cp , thus neglecting
their temperature-dependence. In the example elaborated in section 4, the values will be taken as corresponding to a
water temperature of 60◦ C. We then use the empirical Darcy-Weisbach equation to model the viscous pressure drop
for incompressible flow as a function of the volumetric flow rate qij through a pipe (i, j) with length Lij ,

(pi − pj )gij = qij , ∀ij ∈ Ai . (1)

The flow conductance gij in this expression is based on the correlation of Cheng [35] covering the entire range of flow
regimes from laminar to turbulent, which we modified to behave regularly towards the limit of zero-sized diameters. It
is given by
2f f (1−f3 ) 2f2 (1−f1 )(1−f3 ) 2(1−f1 )(1−f2 )(1−f3 )
g = |q|(f1 −1)(1−f3 ) (d − 2) 3 ΘL1 ΘTS ΘTR ΘC ,

4
 ρ 1
with ΘL = , f1 = 9 ,
16dπµ

Re
1+ ReLT
 
Re 1
ΘTS = 1.8 log10 , f2 = 2 ,
6.8

Re
1+ 320( 2 (2)
d
)
5 2
d π 1
ΘC = , f3 = ,
8ρL d−2 6


1+ 2
 
3.7d
ΘTR = 2 log10 .

Here,  denotes the absolute pipe roughness, d the inner diameter of the pipe and Re the Reynolds number, defined
as Re = 4ρ|q|/(πµd) . The transition from laminar to the turbulent regime is controlled by ReLT = 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 difference
at the node i, at which the flow enters the pipe ij, and θij at the pipe exit. The pipe exit temperature difference θij , due
to heat loss to the environment is given by
 
−Lij
θij = θi exp , ∀ij ∈ Ai , (3)
ρcp |qij |Rt,ij
with Rt,ij the thermal resistance per unit pipe length between the energy carrier and the environment. For a pipe with
an outer isolation casing diameter do,ij that is assumed to be bigger than the inner diameter dij by a fixed ratio, i.e.
do,ij = rdij , the combined thermal resistance of pipe and soil per unit length is [36]

ln(4h/(rdij )) ln r
Rt,ij = + , (4)
2πλg 2πλi
with λi and λg the thermal conductivity of the insulation and the surrounding ground, respectively, and h the depth at
which the pipe is located. The above relations for the pipe temperature drop are validated with experimental data in
Ref. [37].

2.2.2 Pipe junction model

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 qa through the arcs a
connected to the junction, i.e. X X
qa − qa = 0, ∀n ∈ N. (5)
a=(i,n)∈A a=(n,j)∈A

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 depending
on 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
X
(max(qa , 0) θa + min(qa , 0) θn ) −
a=(i,n)∈A
X (6)
(max(qa , 0) θn + min(qa , 0) θa ) = 0, ∀n ∈ N.
a=(n,j)∈A

2.2.3 Consumer model

We now introduce a basic model to estimate the heat transferred to the consumer. We suppose the consumer substation
consists of a bypass and a heating system directly connected to the network (see Figure 2). This leads to two subsets of

5
 Figure 2: Illustration of the consumer model

consumer arcs, namely the set of bypass arcs Acb ∈ Ac and the set of heating system arcs Ach ∈ Ac . Both arcs have a
control valve to regulate the flow. The pressure drop over the consumer heating arc is assumed to be of the form

|qij |qij
pi − pj = ζ 2 , ∀ij ∈ Ach , (7)
αij

with ζ a constant determined from nominal operating conditions [38], and α ∈ [0, 1] the valve control. It is important
to note that the valve control is not trying to model the physical control of the valve, but rather attempts to obtain good
conditioning for the optimization problem, while respecting the lower limit of the pressure drop (αij = 1). Similarly,
the bypass flow is controlled through the relation
qmax,b
qij = βij (pi − pj ) , ∀ij ∈ Acb , (8)
∆pdes,b

with β ∈ [0, 1] the bypass valve control, qmax,b a maximal bypass flow and ∆pdes,b a typical pressure drop over the
bypass.
Conservation of energy in the heating system is given by

ρcp qij (θi − θij ) = Qij , ∀ij ∈ Ach , (9)

with Qij the heat transferred to the house through the heating system. The latter is modelled in the form of the
characteristic equation for radiators [39, 40]:
 nij
θ i − θ ij
Qij = ξij   , (10)
i −θhouse )
ln (θ(θi,j −θhouse )

with θhouse the known temperature difference between the indoor and environment temperature at the house. Values of
the coefficients ξij and nij are tabulated for individual radiators, according to the EN 442-2 norm [41]. The bypass
arcs on the other hand are assumed to be free of heat losses, i.e.

θij = θi ∀ij ∈ Acb (11)

2.2.4 Producer model

In the producer arcs, a fixed input flow qb is imposed as boundary condition for this system of equations. In addition,
a given temperature is imposed for the heat source. This leads to

qa = qb,a , θa = θb,a ∀a ∈ Ap . (12)

To uniquely define the pressures throughout the network, a reference pressure is imposed in one of the producer return
nodes. Since only pressure differences influence the flow solution, the exact choice of reference does not impact the
solution.

6
2.3 Solution procedure
Together with this reference pressure, equations (1-12) jointly lead to a well-defined system that, for given design and
operational parameters, can be solved for the flow rates qij , pressures pi , temperatures θi and θij , and consumer heat
fluxes Qij throughout the network. These variables are called the state variables of the model equations and can be
grouped into a single state vector x. The system is solved by an efficient Newton procedure to obtain second order
convergence. The analytical linearization of the model equations that is required for this also forms the basis for the
adjoint solver (see section 3). The structure of the model equations can be exploited to decouple the solution of flow
rates and pressures at the one hand and temperatures at the other hand. We detail the procedure to efficiently solve this
strongly nonlinear system of equations in B. In order to describe the solver, we need to reformulate the equations in
graph notation as an algebraic system of equations. This is done in A.

3 Combined size and topology optimization of thermal networks

For a given superstructure, we can now use the model presented in section 2 to find an optimal configuration for
a thermal network by optimizing the pipe diameters d, along with the network operation parameters, including the
producer inflows qb , and the consumer controls α and β. We will jointly denote all these design variables with the
design variable vector ϕ.
To formulate the design problem as a mathematical optimization problem, the first crucial step is the choice of
the objective/cost function J(ϕ, x) that is to be minimized. Ideally, one would optimize the network design for the
maximal return on investment or maximal energy efficiency. However, the aim of this paper is rather to demonstrate the
validity of the presented approach. As such, we deliberately choose a cost function that captures only two of the main
contributions, being the cost associated with the installed pump capacity and its operation, and that of the investment
cost related to piping. We thus obtain the following cost function formulation that will be employed for design in a
worst-case scenario: X X
J(ϕ, x) = λP (pi − pj ) qij + C(dij )Lij , (13)
ij∈Ai ij∈Ai

with λP the price of increased pump capacity and operation, and C(d) a function representing the investment cost of
piping placement per unit length as a function of the pipe diameter. Because of the nonlinear dependence of both terms
on 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
5 % range of the desired heat supply Qd,a , i.e.

0.95 Qd,a ≤ Qa ≤ 1.05 Qd,a , ∀a ∈ Ach (14)

Because these constraints depend on the state variables x, they will be referred to as state constraints. They will be
succinctly denoted further using a single generic state constraint vector h(ϕ, x) ≤ 0. Finally, minimal and maximal
values for all design variables are typically imposed, leading to box constraints in the form ϕmin ≤ ϕ ≤ ϕmax . These
constraints are relatively easy to deal with and will be further on only implicitly denoted by an admissible set Φad of
design values, giving ϕ ∈ Φad .

3.1 Optimization problem formulation

We can then write our design problem in the generic optimization problem formulation

min J (ϕ, x) cost function

ϕ∈Φad ,x

s.t. c(ϕ, x) = 0, model equations (15)

h(ϕ, x) ≤ 0, state constraints

with c(ϕ, x) = 0 succinctly denoting the model equations (1-12), as detailed in A.

7
 Since the set of model equations c(ϕ, x) = 0 we defined in section 2 can be solved for a state solution x(ϕ) for
each evaluated network design ϕ, we can use them to eliminate the state variables x in the optimization problem. The
so-called reduced optimization problem then reads

min Jˆ (ϕ)
ϕ∈Φad
(16)
s.t. ĥ(ϕ) ≤ 0, state constraints,

with Jˆ and ĥ the reduced cost function and reduced state constraint vector, respectively. The adjoint optimization
strategy discussed in the remainder of this section gives us an efficient approach to solve this optimization problem.

3.2 Adjoint-based continuous optimization procedure

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 pl = ϕl+1 − ϕl ,

| 1
min ∇Jˆ (ϕl ) pl + p|l Bk pl
pl ∈Φp ={pl | ϕl +pl ∈Φad } 2 (17)
s.t. b l ) + ∇h
h(ϕ b (ϕl )| pl ≤ 0,

where Bk is an approximation of the reduced hessian ∇2 J, ˆ which will be estimated from gradient information in
subsequent optimization iterations with a damped BFGS algorithm [42]. It is clear that the main challenge is thus the
efficient calculation of the cost function and state constraint sensitivities ∇Jˆ and ∇ĥ. For this purpose, the adjoint
approach is used.
Because of the novelty of the adjoint approach in its application to district heating network optimization, we will
elaborate next on how to evaluate these sensitivities using an adjoint approach. For compactness of notation, we will
leave out the operator arguments below if they are clear from the context. In addition, we define the vector function
I = [J, h| ]| that groups the quantities of interest of which we need the sensitivities to solve the quadratic program
|
dÎ
(17) for an optimization step pl . Its sensitivities are then given by the sensitivity matrix ∇Î := dϕ .
Using the chain rule of differentiation, we can decompose this matrix into two contributions, i.e.

dÎ ∂I ∂I dx
= + . (18)
dϕ ∂ϕ ∂x dϕ
It should be noted that while the first term is a quite simple term that only contains the direct dependence of the cost
dx
function on our design variables (zero in our case), the second term features the nonlinear dependence dϕ of the model
dx
equation solution. An expression for dϕ can be found by linearising the model equations c(ϕ, x) = 0 around the
design variables ϕl and model solution x̄l = x(ϕl ) of the current iterate, giving

∂c dx ∂c
=− . (19)
∂x dϕ ∂ϕ
It can be observed that the number of linear algebraic sets of equations that need to be solved to evaluate the complete
dx
sensitivity matrix dϕ scales directly with the number of design variables in ϕ. The adjoint method avoids this in an
elegant way. It can be easily derived by substituting (19) in (18). After rearranging terms we get
"  −|  | #|
dÎ ∂I ∂c ∂I ∂c
= + − .
dϕ ∂ϕ ∂x ∂x ∂ϕ

The term between brackets can be obtained by solving the adjoint equation
 |  |
∂c ∗ ∂I
x =− , (20)
∂x ∂x

8
for the so-called adjoint variables x̄∗l = x∗ (ϕl , x̄l ). The derivative matrix in optimization iteration l then evaluates as
 |  |
∂I ∂c
∇Î (ϕl ) = (ϕl , x̄l ) + (ϕl , x̄l ) x̄∗l . (21)
∂ϕ ∂ϕ
∂c
This latter expression is cheap to evaluate, since ∂ϕ 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 algebraic
system 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 a
scalar-valued cost function, the gradient calculation cost is therefore reduced to a single system solution. In the current
optimization 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 a
variant of the Kreisselmeier-Steinhauser function [33], turning these constraints into a single scalar-valued constraint
ĥks (ϕ, γ) ≤ 0, with
nh
1 X
ĥks (ϕ, γ) = ln eγ[ĥ]i . (22)
γ i=1

In this equation, [ĥ]i represents a specific component of the vector ĥ with length nh , and γ ∈]0, ∞[ is a parameter
that 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 combination
with the adjoint method we effectively reduce the computational costs of the complete gradient matrix evaluation
to approximately two algebraic system solutions, independent of the number of design variables. Thus, a scalable
optimization approach is achieved.

3.3 A numerical continuation strategy for discrete diameter optimization

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 to
perform a true topology optimization we need to include the choice whether or not to put a pipe in a particular arc
of the superstructure. We therefore modify the diameter design variable to either not place a pipe at all (d → 0),
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 = {0, 0.05, 0.125, 0.20} m.1 To achieve this discrete nature without giving
up the benefits of the adjoint approach, we introduce a novel numerical continuation strategy that gradually forces
the continuous optimization method towards these discrete diameter choices. The proposed strategy consists of two
elements: the smoothed projection of the diameters onto the discrete diameter set and the gradual penalization of
intermediate diameter values through the piping cost relation C(d). We will now elaborate these components.

3.3.1 Smoothed multi-projection method for discrete diameter optimization

For the projection of the diameters on the discrete diameter set D, we modify the so-called projection method that is
used 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]

tanh (χη) + tanh (χ(ϕ − η))

ϕ̄ = P (ϕ, η, χ) = (23)
tanh (χη) + tanh (χ(1 − η))

to project the design variable ϕ on a binary decision variable ϕ̄ ∈ {0, 1} that determines whether or not material
should be placed at a certain location in a structure. The variables χ ∈]0, ∞[ and η ∈ [0, 1] are two factors controlling,
respectively, the steepness of the Heaviside approximation and the threshold value above which the variable ϕ is
projected onto to upper limit.
1 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 still

well-defined and their solution numerically well-conditioned.

9
Figure 3: (a) Illustration of the smoothed projection of a design variable ϕ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 here
to be a set of pipe diameters D = {0, 0.05, 0.125, 0.20} m.

We extend this principle here to allow projection on the multiple discrete variables in D. To this end, we construct
an extended projection operation d¯ = Pext (d, χ) that is a piecewise function assembled from rescaled versions of
the projection operation P, as shown in Figure 3. The extended projection then uses the basic projection function
P (d, 0.5, χ) to project inbetween two discrete values, as well as project downwards (i.e. P (d, 0, χ)) for values bigger
than the largest value in D. For a diameter d and discrete design variable set D = {0, d1 , d2 , ..., dn }, the extended
projection operator then becomes
  
d
 d P , 0.5, χ d ≤ d1
 1 d1



 d−dj 
d¯ = Pext (d, χ) = dj + dj+1 − dj P d −d , 0.5, χ dj ≤ d ≤ dj+1 , (24)
   j+1 j
 d−dn
dn + wP w , 0, χ dn ≤ d



with w a parameter set such that d < dn + w for all values of d and where we assume that d is always positive.

3.3.2 Piping cost penalization

In the beginning of this section we introduced C(d) as a function that represents the piping cost as a function
of the pipe diameter d. It should be noted, however, that the only relevant values of this function are those in the
discrete set C = {0, c1 , c2 , . . . , cn } corresponding to the discrete diameters d ∈ D, since a cost can only be provided
for these available pipe sizes. For the values in between, a continuous interpolation is to be provided to account for
the gradient-based optimization approach. We choose this interpolation as to deliberately make ‘intermediate’ values
(d 6∈ D) unfavorable for the optimization. Following this reasoning, we smoothly approximate a piecewise constant
function that features a strong increase in price when the diameter exceeds a value in the set D. Such a function can be
described by the piecewise function
c (d−d )
  
d

 υc 1 P d1 , 0, ω + (1 − υ) nd −d 1 d ≤ d1

  n 1

 d−dj cn (d−d1 )

C(d, υ, ω) = υ cj + cj+1 − cj P d −d , 0, ω + (1 − υ) d −d dj ≤ d ≤ dj+1 . (25)
 j+1 j n 1
( )
    c d−d
 υ c + c P d−dn , 0, ω + (1 − υ) n
 1
dn ≤ d
n w w d −d n 1

Here, cw is an arbitrary high cost value to penalize oversized diameters2 . The variable υ ∈ [0, 1] controls the
interpolation between the piecewise linear- and the piecewise constant function. ω ∈]0, ∞[ controls the steepness of
2 Oversized diameters are prevented with an appropriate choice of the box constraints, represented by the admissible set Φad .

10
 600

450
350

0
0 0.05 0.125 0.2

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 towards
a Heaviside step function (green) by υ, representing the real cost of piping. The grey lines illustrate intermediate steps
in the continuation.

the Heaviside approximation for the piecewise constant function. The process of interpolation and steepness increase
can be seen in Figure 4. Given the ill-posedness that the steep derivatives of this function induce on the optimization
problem, a numerical continuation strategy is used to relax this problem in the initial stages of optimization. This will
be explained next.
3.3.3 Numerical continuation strategy
After including the aggregated constraint, the diameter projection, and the continuous piping cost function C(ϕi , υ, ω),
the optimization problem formulation reads3

min Jˆ (ϕ̄, υ, ω)
ϕ∈Φad

s.t. ĥks (ϕ, γ) ≤ 0, aggregated constraint, (26)

ϕ̄ = Pext (ϕ, χ) ≤ 0, extended projection,

in which the continuation variables υ, ω, γ, and χ control the trade-off between relaxing the discrete problem that is to
be solved to increase the well-posedness of the optimization problem, and the exactness of its approximation. We will
therefore apply numerical continuation to gradually force the design towards the discrete solution we are interested
in. That is, we add an outer loop to the optimization in which the above optimization problem is solved for different
values of the continuation variables, and initialize the optimization procedure with the solution of the last continuation
iteration 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. Using
the chain rule of differentiation, we find that the derivative of the quantities of interest is then given by

dÎ dÎ ∂Pext

= , (27)
dϕ dϕ̄ ∂ϕ
3 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. ϕ̄ = Pext (ϕ, χ), in such a way that only the discrete diameter variable components of ϕ are affected.

11
Figure 5: a) Initial topology and configuration of the network. The green triangles represent heat producers of Tb,2 =
70 °C (left) and Tb,1 = 65 °C (right). Consumers are represented by black circles of varying sizes corresponding to
their respective heat demand (5 kW, 15 kW, 50 kW). b) Total length of all placed discrete pipe diameters in the feed
network.

Î
where ddϕ̄ can be efficiently calculated by the adjoint gradient calculation (21) as before and multiplied with the
projection gradient.
4 Demonstration of the optimal design approach on a fictitious design prob-
lem
In this section, we will demonstrate the functioning and scalability of the network optimization algorithm on a fictive
test 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 (Tb,1 = 65 °C, Tb,2 = 70 °C).
160 consumers of 3 different consumer types are distributed throughout the neighbourhood. A summary of the different
consumer characteristics can be found in Table 1. As can be seen in Figure 5, the pipe superstructure for the network
optimization is placed onto a part of the neighborhood’s street grid. In this test case, the total amount of design variables
is 632. To the best of the authors’ knowledge, it is the largest reported application of topology optimization for district
heating networks that is based on a fully-fledged flow and heat transport model.

Table 1: Consumer characteristics

Pa s
Type Qd,a Tin,d ∆Td nd ξ [ KkW
nd ] ζ [1012 m3
]
Dwelling 15 kW 55 °C 20 °C 1.2 0.34 1.5
Renovated dwelling 5 kW 55 °C 20 °C 1.42 0.06 0.14
Commercial demand 50 kW 55 °C 20 °C 1.2 1.23 13.54

Assuming a worst-case situation with an outdoor temperature T∞ = −8 °C, we now try to design the network to
have the least possible investment cost, as well as pump capacity and operation cost. The weight for pump-related
costs was chosen to be λP = 102 Me/kW in the current cost function definition, in order to limit the pressure drops
in the piping to reasonable engineering limits. The design is initialized by distributing a uniform pipe network with a
diameter of d = 0.15 m over the entire superstructure (see Figure 5). Initially the control valves (heating system and
bypass valves, as well as producer inflows) are fully opened (α = 1, β = 1, qb = qb,max ). To start the optimization

12
within the heat satisfaction boundaries set by (14), we first optimize the network operation parameters for minimal
consumer heat dissatisfaction, i.e.
2
J(ϕ, x) = (Qa − Qd,a ) , (28)
while fixing the pipe diameters. To choose reasonable consumer heating characteristics, we suppose the consumer
heating systems are designed to provide this heat load at nominal/design conditions, described next. First of all, we
assume the design is based on a desired temperature drop ∆Td and pressure drop ∆pd of 50 kPa over the heating
system. Next, we take the heating system inflow temperature Tin,d and a room temperature of Troom,d = 20 °C. We
can then obtain the heating system characteristics from
 −nij
Tin,d − Tout,d
ξ ≈ Qd  (T −θ )  with Tout,d = Tin,d − ∆Td ,
in,d house
ln (Tout,d −θhouse ) (29)
∆pd ∆pd ∗ (ρcp ∆Td )2
ζ= = ,
qd2 Q2d

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 Tb,1 = 65 °C and Tb,2 = 70 °C.
Consequently, we can directly relate the controlled volumetric inflow qb at the producers to the heat inflow via

Qin = ρcp qb (Tb − Tout,d ), (30)

with Tout,d = 40 °C the design outflow temperature of the network. We can thus directly control the maximal input
heat delivered by the producers with a simple constraint on the inflow qb . Given a total heat demand of 1.77 MW for
the 160 consumers, we can only satisfy the demand if the combined maximal heat provided by both producers is big
enough 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 to
optimize 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·103 , 105 ],
χ ∈ [0, 100], υ ∈ [0, 1] and ω ∈ [0, 100]. The maximal heat input Qin,max,1 = Qin,max,2 = 2 MW of each individual
producer 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 C

D [m] 0 0.032 0.065 0.1 0.15 0.2 0.3 0.4
C [e/m] 0 2202 2218 2258 2448 2461 2665 2922

For this test case with 632 design variables, the optimization procedure, including warm-start, converged after
46 min. The optimization was run on a single CPU of a standard laptop with an Intel Core i7 processor (1.9 GHz). A
decrease in the piping investment cost by 23 % from 88.8 Me to 68.7 Me with respect to the uniform pipe distribution
design 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 0.15 m pipes throughout the superstructure and that the optimized
operation parameters satisfy the consumer heat demand up to a very high tolerance. The pump operation cost was
reduced by a factor of 14 with respect to this ‘uniform’ pipe design. The resulting network topology is given in
Figure 6a. It can be directly observed that larger pipe diameters are provided for high demand consumers, due to
their 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 options
can be seen in Figure 6b. It should be noted that the topology optimization algorithm explicitly decided not to install
any piping in about 3 km of the superstructure.
The heat demands of all consumers are met with a producer inflow Qin,1 = 0.92 MW and Qin,2 = 1.09 MW.
Minimizing the pump operational cost (13) the algorithm is able to push the supplied consumer heat Qa close to the
lower limit of the thermal comfort range of the consumer (14). It is intuitive that the cost of the network design
will indeed be lowest if these constraints are only marginally satisfied. Indeed, through the increase of the constraint

13
Figure 6: a) Optimal network topology and configuration. The green triangles represent heat producers of Tb,2 = 70 °C
(left) and Tb,1 = 65 °C (right). Consumers are represented by black circles of varying sizes corresponding to their
respective heat demand (5 kW, 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 pipe
diameters in the feed network.

smoothness parameter γ in the numerical continuation strategy, the algorithm achieves a maximum deviation of the
consumer heat loads from the lower limit of 0.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. The
temperature 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 the
colder parts of the network, the flow rate through the consumers heating system has to be increased (as can be seen
in Figure 7b). The increase of flow rate in this test case reaches up to 72 % when comparing the flow rate for the
highest inflow temperatures of Tin = 69.5 °C to the lowest inflow temperatures of Tin = 56.9 °C, for the consumer
type ‘dwelling’. It should be noted that this increase in flow rate influences the cost balance between investment and
pump operation cost, ultimately driving the pipe design towards bigger diameters.

5 Discussion
We presented a methodology for the automated design of district heating network topologies on the basis of a nonlinear
model of flow and heat transport. Because the computational cost of the adjoint sensitivity calculations are intrinsically
independent of the number of considered design variables, the methodology is expected to scale well to large heating
networks.
This has been showcased with the design of a fictitious district heating network with a scale well beyond any
other applications of nonlinear topological design tools. On a fictitious design problem of a medium-sized district in
Waterschei, the method shows to correctly retrieve a discrete network design that is able to reduce the piping investment
by 23% and the cost related to installed pump capacity and operation by a factor of 14 compared to the initial uniform
pipe distribution. While further extensions of the cost function towards the direct optimization of the net present value
could still improve the relevancy of the resulting designs, they can easily be incorporated in the current methodology
without impacting the computational cost.
The computational cost in this proof-of-concept design remained well below an hour on a standard laptop. The

14
Figure 7: a) Average pipe temperature throughout the network. A clear separation in a low temperature network and
a high temperature network can be observed. b) Temperature dependency of the flow rate required by the consumers
heating 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 approach
was tested for a constraint that enforces the heat delivered at each single consumer to deviate less than 5 % from its heat
demand 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 0.008 % from the solution that is expected if imposed
in a non-aggregated way.
The added value of considering a full nonlinear model in the design clearly follows from the observation that the
consumer 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 to
capture already in the design phase. Linearized models, in contrast, fail to describe the flow pattern in networks with
multiple producers or with loops, as well as the heat losses throughout the network and the temperature-dependent
consumer 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 the
novel method and its scaling to that of other mixed-integer non-linear programming methods. It is important to note
that gradient-based optimization methods like the one presented in this paper could be more sensitive to find local
optima than other mixed-integer non-linear programming methods. It should therefore be thoroughly investigated for
a number of practical design cases how good the resulting designs perform with respect to those resulting from other
design approaches. Finally, it should be mentioned that several uncertainties play a crucial role in the design of thermal
networks, such as producer unavailabilites, possible future network extensions, or potential pipe failures. To be of truly
practical value, optimal design methods for district heating networks should eventually account for these uncertainties.

6 Conclusions and outlook

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 design
variables, we introduce a numerical continuation strategy that combines a smoothed projection of the diameters with a

15
penalization 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 for
the first time.
The scalability of the method is demonstrated with the design of the –to the best of the authors’ knowledge– largest
topological network design application reported thus far on the basis of a non-linear transport model. This is enabled
here through the intrinsic design-independent scaling of the adjoint method if combined with the suggested constraint
aggregation approach. Moreover, by building the optimization on the basis of a nonlinear physics model, heat losses
and 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 impact
of 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 programming
tools 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 Flemish
Institute 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 to
reformulate the equations (1-12) in vector form. Consider a graph with n nodes and m directed arcs. Let us start by
defining the following vectors for flow rates, pressures, temperatures and consumer heat fluxes:

p = pi , i ∈ N,
q = qij , ij ∈ A,
!
θn (31)
θ= , with θn = θi , i ∈ N and θa = θij , ij ∈ A,
θa
Q = Qij , ij ∈ Ach .

Using the same subscripts as used before for the subsets, we introduce vectors related to only one of the subsets
Ai , Ac , Ach , Acb , Ap of the arcs A in a similar way, e.g. qi = qij , ij ∈ Ai denotes the flow in the internal arcs. Next,
we introduce the node-pipe incidence matrix

A = (a1 , . . . , am ) ∈ Rn×m

that contains information on which nodes are connected by a directed arc. That is, each column ai contains “1” at the
start node i of arc ij, and a “-1” at the end node j. The continuity condition in the pipe junctions (5) is then easily
translated to vector form as Aq = 0. The pressure drop over all arcs on the other hand is given by ∆p = A| p. Also
here we can easily extend the definition to incidence matrices that include only a number of rows related to a specific
subset of arcs, e.g. Ai = (A)i,j , i ∈ N and j ∈ Ai .
Using this notation, we can formulate the state equations related to the hydraulic flow succinctly in vector form
as H(ϕ, y) = 0, with y = ( pq ) the hydraulic state vector and ϕ the vector of design variables (defined further in
section 3). That is, the vector  
Aq
H(ϕ, y) = Ai | p − Ri (ϕ, q)qi  = 0, (32)
 
|
Ac p − Rc (ϕ, q)qc
then contains the residuals of node continuity, pipe momentum, and consumer arc momentum balances, respectively.
The diagonal matrices Ri (ϕ, q) and Rc (ϕ, q) hereby represent the nonlinear flow resistance for internal and consumer
arcs, respectively, as can easily be extracted from relations (1), (7), and (8). Note that we leave out the trivial producer

16
conditions and reference pressure condition for simplicity of notation, although they are in principle an integral part of
the vector equation H(ϕ, y) = 0.
To obtain a similar vector equation for energy conservation, we introduce the arc inflow and outflow matrices

Ain = max(A diag(sgn(q)) , 0)

(33)
Aout = min(A diag(sgn(q)) , 0) ,

with “sgn” the elementwise sign function and “diag” an operator that maps some argument vector d with elements di
to a diagonal matrix D with diagonal elements (D)i,i = di . Energy conservation in the pipe junctions, pipes, consumer
heating arcs, and bypass arcs can then be succinctly written as a vector equation4
 
diag(Ain |q|) θn + Aout diag(|q|) θa
Ai,in | F(ϕ, y)θn − θa
 

E(ϕ, y, θ) =  ◦n  = 0, (34)
 
|
 
ρcp diag(qch ) (Ach,in | θn − θa,ch ) − ξ Ach,in θ2n +θa,ch − θhouse 
 
Acb,in | θn − θa,cb

with F(ϕ, y) a diagonal matrix that has the fractional reduction of θ through the pipes on its diagonal, as can be
deduced 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 further
elaborated in the next subsection. In an even more compact form, we can thus represent the thermal network model as
the vector equation !
H(ϕ, y)
c(ϕ, x) = = 0, (35)
E(ϕ, y, θ)
p
y

with x = θ = q 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 to
iteratively solve the model
 equations (35) for an update of the state variables. That is, in each iteration “k”, we obtain
∆yk
an update ∆xk = ∆θk from

∂H
 
(ϕ, yk ) 0    
 ∂y  ∆yk H(ϕ, yk )
 ∂E ∂E  ∆θ = − (36)
k E(ϕ, yk , θk )
(ϕ, yk , θk ) (ϕ, yk , θk )
∂y ∂θ

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 ȳ using updates ∆yk from

∂H
(ϕ, yk )∆yk = −H(ϕ, yk ). (37)
∂y

2. Then plug this solution into the energy equations and solve for the thermal state θ̄ using

∂E
(ϕ, ȳ, θk )∆θk = −E(ϕ, ȳ, θk ). (38)
∂y
4 The operator (a)◦b denotes the Hadamard exponential, which simply represents the pointwise exponent of a vector a to the power b.

17
This significantly reduces the size of the linear systems that need to be solved.
∂H ∂E
It should be noted that the application of this Newton scheme requires the derivatives and . These deriva-
∂y ∂θ
tives can be obtained by analytical linearization of (32) and (34) with respect to the state variables. Though this is a
somewhat 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 adjoint
solver.

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 elaborate
on the efficient implementation of the adjoint solver, given the thermo-hydraulic model at hand. Rewriting the adjoint
equation in terms of H and E reveals a similar block-triangular structure as in the state equations:
 |  |   | 
∂H ∂E  ∗ ∂I
 ∂y |  y∗ = −  ∂y |  ,
 ∂y   
 (39)
 ∂E  θ  ∂I 
0
∂θ ∂θ
 ∗
with x∗ = yθ∗ .
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
 |  |
∂E ∗ ∂I
θ =− . (40)
∂θ ∂θ

2. Then plug this solution into the adjoint hydraulic equations and solve for the adjoint hydraulic state y ∗ using
 |  |  |
∂H ∂I ∂E
y∗ = − − θ∗ . (41)
∂y ∂y ∂y

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 model
equations. Since these equations follow from linearization around (ϕ, x̄), the matrices on the left-hand side of the
adjoint 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 performed
in only a fraction of the time needed for a model evaluation, independent of the number of design variables.

References
[1] U. Persson, B. Möller, and S. Werner, “Heat roadmap europe: Identifying strategic heat synergy regions,” Energy
Policy, 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 Generation
District Heating (4GDH),” Energy, vol. 68, pp. 1–11, Apr. 2014.
[4] J. Söderman, “Optimisation of structure and operation of district cooling networks in urban regions,” Applied
thermal engineering, vol. 27, no. 16, pp. 2665–2676, 2007.
[5] J. Dorfner and T. Hamacher, “Large-Scale District Heating Network Optimization,” IEEE Transactions on Smart
Grid, vol. 5, pp. 1884–1891, July 2014.

18
 [6] C. Haikarainen, F. Pettersson, and H. Saxén, “A model for structural and operational optimization of distributed
energy 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, operation
and 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 system
considering 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ère, and M. Santarelli, “Exergy analysis and thermo-
economic optimization of a district heating network with solar-photovoltaic and heat pumps,” in 32nd Interna-
tional Conference on Efficiency, Cost, Optimization, Simulation and Environmental Impact of Energy Systems,
ECOS 2019, 23 June 2019 through 28 June 2019, 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 heat
transportation 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 case
study,” 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 a
district 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 model
for 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 multiparametric
delay modeling,” Energy, vol. 189, p. 116125, 2019.

[19] F. Marty, S. Serra, S. Sochard, and J.-M. Reneaume, “Simultaneous optimization of the district heating network
topology 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 of
a 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 of
a 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 topology
optimization for deep energy savings at the urban district level,” Energy Conversion and Management, vol. 222,
p. 113106, 2020.

19
[23] J. von Rhein, G. P. Henze, N. Long, and Y. Fu, “Development of a topology analysis tool for fifth-generation
district 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 numerical
methods in fluids, vol. 41, no. 1, pp. 77–107, 2003.
[27] T. E. Bruns, “Topology optimization of convection-dominated, steady-state heat transfer problems,” International
Journal 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 automated
parameter, 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,” Computer
Methods 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,” Journal
of 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 Applied
Mechanics and Engineering, vol. 289, pp. 332–354, June 2015.
[34] J. K. Guest, J. H. Prévost, and T. Belytschko, “Achieving minimum length scale in topology optimization using
nodal 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üller, C. Nytsch-
Geusen, M. Wetter, and L. Helsen, “Dynamic equation-based thermo-hydraulic pipe model for district heating
and 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 analyses
of a district heating network,” Energy, vol. 57, pp. 149–159, Aug. 2013.
[39] 2012 ASHRAE Handbook, 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ür den Maschinenbau. Springer-Verlag, 2014.
[41] “EN 442-2:2014; radiators and convectors – Part2: Test methods and rating,” standard, European Committee for
Standardazation, 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 in
topology optimization,” Structural and Multidisciplinary Optimization, vol. 43, pp. 767–784, June 2011.

20