\section{Introduction}

\gls{dhn} are considered one of the core technologies to enable carbon-neutral space heating \cite{OECD/IEA2019}. It has the ability to connect a multitude of different renewable heat sources and provide heat to districts and entire cities. In \gls{dhn}s, the typically high upfront investment cost of groundworks and piping is a core decision variable for the feasibility of a development project. Therefore in the planning phase it is crucial to design the pipe routing (network topology), pipe sizing and heat production capacities in an optimal way. This topology optimization problem solves the question of where to place heating network pipes and at what diameter, while accounting for the non-linear nature of the heat transport physics. Together with the binary choice  of pipe placement it is therefore inherently a \gls{minlp}. 

Solving non-linear combinatorial problems directly can be challenging. As a result, the topology optimization problem of heating networks is often linearized, resulting in a linear optimization problem that can be solved efficiently using Mixed Integer Linear Programming (MILP) solvers.  S\"oderman \cite{Soderman2007}, for example, optimized the structure and configuration of a district cooling network using a MILP approach. Dorfner and Hamacher \cite{dorfner2014large} later used a linear approach to optimize the topology and pipe sizes of a DHN. Haikarainen et al. \cite{Haikarainen2014} optimized the topology and operations, while accounting for different production technologies and heat storage. Mazairac et al. \cite{Mazairac2015a} optimized the topology of a multi-carrier network that incorporates gas and electricity  supply. Morvay et al. \cite{Morvaj2016} optimally designed the network while also optimizing the energy mix supplied. Still using a MILP approach, Bordin et al. \cite{bordin2016optimization} optimized the network topology while studying the set of consumers to optimally connect to a DHN. Most recently, Resimont et al. \cite{Resimont2021} used a MILP approach to optimize a city-scale heating network. The transformation of modern \gls{dhn}s towards multi-source, low-temperature networks \cite{Lund2021} invalidates the assumptions of most linear heating network models. Accurate modeling of heat losses and the different feasible temperature levels of supply and demand sites in \gls{dhn}s requires a non-linear representation of the network physics. Solving the resulting \gls{minlp} problem is challenging. In the literature, one approach to solving this \gls{minlp} problem is to use heuristic search algorithms. For example, Li and Svendsen used a genetic algorithm \cite{Li2013}, and Allen et al. deployed both a minimal spanning tree heuristic and a particle swarm algorithm \cite{Allen2022} to the topology optimization of a small DHN.

To fully leverage the potential of mathematical optimization notably two approaches have been successfully applied to the non-linear optimal topology problem of heating networks. Mertz et al. \cite{mertz2016minlp}, solves the full \gls{minlp} resolving the discrete nature of pipe routing choices, using a combinatorial optimization approach. On the other hand, Blommaert et al. \cite{Blommaert2020a} and Wack et al. \cite{wack2022}, solve a relaxed \gls{nlp} ensuring a discrete network topology through penalization. This method is inspired by partial differential equation (PDE) constrained topology optimization, where penalization methods are combined with adjoint optimization to solve problems governed by PDEs. These methods promise to scale well with the problem size, allowing to optimize large-scale problems.

The current bottleneck to a widespread use of optimization methods for the design and topology optimization of \gls{dhn}s are an accurate representation of the network while maintaining scalability of the approach. The goal of this paper is therefore to compare the performance of the combinatorial approach by Mertz et al. \cite{mertz2016minlp} and the relaxed adjoint approach by Wack et al. \cite{wack2022} to non-linear topology optimization of \gls{dhn}s. First, the computational cost scaling of both optimization approaches with the network size is benchmarked to assess whether they are suited to be used for large-scale \gls{dhn} development projects. Both approaches solve non-linear and non-convex optimization problems, so the found optimal network designs can be different local optima. Therefore, in a second step, the found optimal designs are compared. Here, the influence of resolving the discrete nature of the pipe routing problem, as well as the influence of initialization strategies on the found local optima is discussed. This comparison is performed for a simple single heat producer case and, to closer mimic modern heating network problems, for a multi-producer case with different heat production temperatures.

\section{The topology optimization problem of \gls{dhn}s}\label{sec:OptProblem}
To compare both approaches, a topology optimization problem is defined that aligns the problem definitions of Mertz et al.\cite{mertz2016minlp} and Wack et al \cite{wack2022}. This ensures comparability of the results with both approaches. For completeness, the optimization problem definition is repeated here in short, for detailed discussion the reader is referred to the afore-mentioned publications.

In order to represent the design and topology optimization problem of a \gls{dhn}, a set of design variables $\defDesignVar$ is defined, containing the pipe diameters $\ve{\topVarSkalar}$ as well as the normalized producer inflows $\ve{\prodInput} $. To represent the physical state of a given network, a vector of physical variables $\defStateVar$ is defined. It contains the flow rates $\ve{\flow} $, nodal pressures $\ve{\pressure} $ and nodal and pipe exit temperatures $\ve{\dTinf} $. The temperatures $\ve{\dTinf} = \ve{T} - \temperature_\infty$ are defined as the difference between the absolute water temperature $\ve{T}$ and the outside air temperature $\temperature_\infty$. Now the topology optimization problem for a \gls{dhn} can be posed as a generic optimization problem of the form:

\begin{equation} \label{eq:opt}
	\begin{aligned} \min_{\ve{\designVarskal},\ve{x}}& \qquad
		\mathcal{J} \left(\ve{\designVarskal},\ve{x}\right) \\
		s.t.& \qquad \ve{h}(\ve{\designVarskal},\ve{x}) \leq 0.  \\
	\end{aligned}
\end{equation}

Here $\mathcal{J}$ represents the cost function and is defined as the total cost of the project over an investment horizon of $A = 30 ~\mathrm{years}$. $\ve{h}(\ve{\designVarskal},\ve{x})$ defines the set of model and technological constraints. A full definition of the cost function and model constraints for this benchmark can be found in \ref{app:model}. This optimization problem is based on Mertz et al. \cite{mertz2016minlp} and Wack et al. \cite{wack2022} and a detailed description of the underlying costs and models can also be found there. In the problem definition, the set definition in table \ref{tab:sets} is used to reference different parts of the network.

\begin{table}
	\centering
	\caption{Notation used to describe different subsets of the \gls{dhn}.}
	\label{tab:sets}
	\begin{tabular}{ll}
		\hline
		Set	& Symbol  \\
		\hline\hline
		All edges &  $E$ \\  
		Producer input & $\setEdges_\pro$ \\  
		Pipes& $\setEdges_\pipe$ \\
		Feed edges &  $\setEdges_{\mathrm{f}}$  \\
		Heating system & $\setEdges_\rad$  \\
		\hline
	\end{tabular}
\end{table}

\section{Methodology}
To compare the two approaches, the topology optimization problem (described in equation \ref{eq:optProblem}) is solved using two different methods: solving the full \gls{minlp} problem in GAMS using the implementation by Mertz et al. \cite{mertz2016minlp} (referred to as fMINLP{} for simplicity), and using the penalized adjoint optimization approach by Wack et al. \cite{wack2022} (referred to as pNLP{}). Both implementations were simplified to facilitate comparison. The pNLP{} implementation no longer accounts for multiple discrete pipe diameters, while the fMINLP{} implementation no longer optimizes the production side or cascading supply between high and low temperature consumers. The following sections briefly discuss the two methodologies and their different treatment of the binary pipe existence variable.

\subsection*{Combinatorial \gls{minlp} approach: fMINLP{}}
For the comparison, the topology optimization problem will be solved using a combinatorial \gls{minlp} approach. This approach was previously published by Mertz et al. \cite{mertz2016minlp}, and a comprehensive discussion of the approach can be found there. It is implemented in GAMS and uses the \gls{milp} solver CPLEX and the \gls{nlp} solver CONOPT. To pose the optimization problem as a combinatorial problem, the topological choice of pipe placement is represented by a binary variable $\ve{\phi}$. Physical variables defined on these pipes are then coupled to this existence variable using the bigM method ($\mathcal{M}$) \cite{mertz2016minlp}, as e.g. in the flow velocity definition: 

\begin{align}\label{eq:optProblem}
	&\gijEdge{m} -\mathcal{M}\left( 1-\gijEdge{\phi} \right)\leq\gijEdge{v} \density \pi \frac{1}{4}\gijEdge{d}^2\\ &\gijEdge{v} \density \pi \frac{1}{4}\gijEdge{d}^2\leq \gijEdge{m} +\mathcal{M}\left( 1-\gijEdge{\phi} \right) , \quad \forall \gi\gj \in \Epipe.
\end{align}
Here $\mathcal{M} \gg \gijEdge{m}$ is chosen in such a way that if a pipe exists $\gijEdge{\phi} = 1$, the bounds on the flow velocity constraint are tight. If $\gijEdge{\phi} = 0$, the bounds of this constraints are relaxed. To force the installation of both feed and return pipes the following constraint is used:
\begin{equation}
	\gijEdge{\phi}= \gEdge{\phi}{j}{i},\quad\forall \gi\gj \in \Epipe.
\end{equation}

The fMINLP{} implementation features additional constraints to facilitate convergence. They are defined in a way to not be active in the final optimal design in order to maintain the same final problem formulation as the pNLP{}. First, the installed capacity $H$ is bound by the sum of consumer demands: 

\begin{equation}
	\sum_{ij\in\setEdges_{\pro,\mathrm{f}}}\gijEdge{H} =\sum_{ij\in\setEdges_{\pro,\mathrm{f}}}\gijEdge{\left(q\dTinf  \right)} \density c_{\mathrm{p}} \leq 1.5 \sum_{kl\in \setEdges_\rad}\gEdge{Q}{k}{l}.
\end{equation}

To facilitate convergence on the momentum equations, the velocity and pressure drop over a pipe are bound by:

\begin{align}
	\gijEdge{v} &\leq 3.5, \\
	\gijEdge{v} &\leq 18.438 \gijEdge{d} + 0.2186, \\
	\gijEdge{v} &\geq v_{\mathrm{min}}\gijEdge{\phi},\\
	\giNode{p} - \gjNode{p}&\leq\SI{200}{\pascal} \quad \forall \gi\gj \in \Epipe.
\end{align}

This constraints are not active in the final optimized design, as they are dominated by the pipe momentum equation (\ref{eq:pipeMomentum}). A fixed pressure drop is set for the heat exchanger

\begin{equation}
	\giNode{p} - \gjNode{p}\geq\SI{2}{\kilo\pascal}, \quad \forall \gi\gj \in \Erad,
\end{equation}

and the exit temperature of the radiator is required to be smaller than the inflow temperature:
\begin{equation}
	\giNode{\dTinf} \geq \gjNode{\dTinf} , \quad \forall \gi\gj \in \Erad.	
\end{equation}

Finally, this optimization strategy uses a sequence of initializations to ensure a stable convergence of the final \gls{minlp} solve, based on the strategy of Marty et al. \cite{marty2018simultaneous}. For completeness, this initialization sequence is briefly described in \ref{sec:Apendix_init}.

\subsection*{Adjoint-based penalization approach: pNLP{}}

The other considered approach is the adjoint-based penalization approach by Wack et al. \cite{wack2022}. Here the combinatorial problem is relaxed, allowing for a continuous pipe placement choice. A near-discrete topology is then enforced by penalization, e.g. by replacing $p_0$ in the investment cost function (equation \ref{eq:Jpipe}) with
\begin{equation}
	\bar{p}_0(\gijEdge{d})= p_0 \left(\frac{1}{(1+\exp(-k\left(\gijEdge{d}-\diameter_{\mathrm{min}}\right))})-1\right)\,,
\end{equation}
or by explicitly penalizing the $d$ in the investment cost function and within the momentum and energy equations as described in Wack et al. \cite{wack2022}. 

Because of this relaxation, the \gls{minlp} reduces to solving a series of \gls{nlp}s. Here, the optimization is initialized from a uniform distribution of pipe diameters.

\section{Benchmarking computational cost - a single producer case}\label{chap:simple}

Modern \gls{dhn}s grow ever bigger and more complex, including multiple heat production sites at different temperatures, which optimal design tools need to be capable to handle. These tools therefore need to scale well with the heating network size in order to be applicable to relevant cases. In this first benchmark, the computational cost scaling of a combinatorial \gls{minlp} approach (fMINLP{}) is compared to solving a relaxed penalized \gls{nlp} (pNLP{}). This is done on a benchmark case containing a single producer. 

\subsection*{Setup}

To be able to compare the computational cost scaling of both optimization approaches with increasing network size an easily scalable \gls{dhn} optimization case is set up. This benchmark case is visualized in figure \ref{fig:simpleSetup}. 

\begin{figure}[h]
	\includegraphics[width=1\columnwidth]{fig01_SingleProducerSetupTextwidth}	
	\caption{Setup of the first benchmark case. Around a central producer, heat consumers (houses) are arranged in a circular pattern. They are connected by possible pipe routes (black lines). Pipe junctions are visualized as black circles. To increase the size of this case, additional segments $s$ are successively added to the outside of the network.}
	\label{fig:simpleSetup}
\end{figure}

Here, a heating network is to be designed around a single heat producer in the center of the network. This producer provides heat at $\SI{70}{\degreeCelsius}$. Around this central producer, houses and possible pipe connections are arranged in a circular manner. All houses have a heat demand of $Q = \SI{15}{\kilo\watt}$, and their heating system characterized by $\heaterCoef=\SI{200}{\watt\per\kelvin^\heaterExp}$ and $\heaterExp=1.2$. To investigate the cost scaling, the size of this circular network is successively increased by adding additional segments $s$. With each segment, additional heat consumers and 5 additional potential pipe connections are added, with the number of potential pipes following a linear scaling: $n\left( s \right) = 5s+13$. In figure \ref{fig:simpleSetup} the addition of segments $s = 10$ and $s=11$ are visualized in red and blue.

Now the benchmark of both approaches is conducted by creating a series of heating networks of increasing size. Segments $s$ are added in steps of 10, starting from $s=0$ up to $s=190$ therefore creating a sequence of heating networks with $(s_k)_{k=0}^{19}, ~s_{k}=10k$ segments. This sequence of optimization problems is solved by both the fMINLP{} implementation as well as the pNLP{}. This optimizations were performed on the same computer (Using a single Intel Xeon 3.20 GHz processor core). 

\subsection*{Comparison of the computational cost}
The optimization was repeated 3 times for each network and the mean runtime until convergence (wall time) was recorded. The wall time for each network size is visualized in figure \ref{fig:simpleScaling} on a a) semi-log graph and a b) log-log graph.

\begin{figure*}[h]
	\centering
	\includegraphics[width=0.736\linewidth]{fig02_SimpleExampleScaling}	
	\caption{Wall time scaling of the two approaches with increasing network size $n$. The scaling of the approach in the fMINLP{} is described well by an exponential fit a), while the scaling of the pNLP{} can be described well by a power fit b). The pNLP{} did not converge for case $n=713$ while the fMINLP{} did not converge for every case $n\geq663$. All non-converged optimizations are excluded from the figure and subsequent analysis}
	\label{fig:simpleScaling}
\end{figure*}

The figures show that the wall time for the fMINLP{} implementation scales faster with the network size then for the pNLP{}. In order to get an understanding of how long it would take both approaches to optimize networks of considerable size, it is useful to understand what function governs their cost scaling. Considering the algorithmic differences of solving a \gls{minlp} in the case of the fMINLP{} vs solving a relaxed \gls{nlp} in the case of the pNLP{}, both an exponential and a polynomial time scaling was tested on both approaches. First, an exponential fit was performed for the scaling of both approaches (see figure \ref{fig:simpleScaling} a)). The wall time scaling of the fMINLP{} can here be described by the function $\walltime_{\mathrm{\gams{}}}\left(  n\right) = 51.43e^{0.012n} \si{\second}$ with an coefficient of determination $R^2=0.87$, while the wall time scaling of the pNLP{} can be described by $\walltime_{\mathrm{\pathopt{}}}\left(  n\right) = 47.31e^{0.006n}\si{\second}$ with a coefficient of determination of $R^2 = -0.01$. From this fit, it can be concluded that the wall time scaling of the fMINLP{} approach can be described well by an exponential function, while the wall time scaling of the pNLP{} approach cannot. Second, a power function was fitted to the wall time scaling of both approaches. Here, the scaling of the pNLP{} can be described with $\walltime_{\mathrm{\pathopt{}}}\left(  n\right) = 0.04n^{1.72}$ with a coefficient of determination of $R^2=0.85$ and the scaling of the fMINLP{} can be described by $\walltime_{\mathrm{\gams{}}}\left(  n\right) = 0.16n^{1.75}$ with $R^2=0.04$. It can be concluded that the wall time scaling of the pNLP{} follows a polynomial function reasonably well, while the fMINLP{} cannot be described with a polynomial function. 

The comparison above highlights the potential exponential solution time of solving full \gls{minlp}s. In the context of \gls{dhn} optimization, this steep cost scaling does not only make design studies slow and expensive, its exponential nature renders optimizations of considerably sized networks intractable and serves as a bottleneck for the optimization of large-scale \gls{dhn}s, containing thousands of potential pipe connections. For example to optimize a network containing $n = 2000$ potential pipe connections, an optimization time of $\walltime_{\mathrm{\gams{}}}\left(  n\right) \approx 43000 ~\mathrm{years}$ would be needed following the observed exponential trend of the \gls{minlp} solution time.
In the context of topology optimization of \gls{dhn}s it is computationally favorable to solve a relaxed \gls{nlp}, as is evident from the polynomial wall time scaling of the pNLP{} implementation. The polynomial scaling of the computational solution time with this approach makes large-scale \gls{dhn} optimization feasible. This is highlighted for the above mentioned example of a network of $n = 2000$, where the necessary optimization time reduces to $\walltime_{\mathrm{\pathopt{}}}\left(  n\right) \approx 5.3 ~\mathrm{hours}$ when using the relaxed \gls{nlp} implementation of the pNLP{}, assuming the observed polynomial trend.

The potentially exponential computational cost scaling of solving full \gls{minlp}s should be taken into account when choosing optimization strategies for \gls{dhn} topologies. It will make the optimization of large-scale \gls{dhn}s intractable and hinder the application of such optimization tools to real-world \gls{dhn} problems. It is therefore advisable to simplify the \gls{minlp} formulation by either linearizing the non-linear model constraints, posing the problem as a \gls{milp}, or by relaxing the integer constraints, reformulating the problem as a \gls{nlp}. This second relaxation is done within the pNLP{} implementation, and while keeping the non-linear model in the overall optimization procedure, its significant speed up in the computational cost is shown in this paper.  

\subsection*{Comparison of found optimal network topologies}

Both approaches solve a non-linear optimization problem. This is necessary because the non-linear heat transport physics can have a big influence on the network design. Non-linearities in the optimization problem, however, can lead to multiple local optima. The proneness of both optimization approaches to these local optima is therefore studied in a next step. Here it can be investigated, if the relaxation of the integer constraints of pipe placement by the pNLP{} has a relevant influence on the cost of the found optimal network design. For this, the optimality gap of both approaches is compared. The optimal total annualized cost for each benchmark step found with both approaches is visualized in figure \ref{fig:simpleCost}.\footnote{To prevent potential remaining modelling differences from skewing the comparison, both found optimal designs where evaluated in a simulation within the pNLP{}.}

\begin{figure}[h]
	\includegraphics[width=1\columnwidth]{fig03_SimpleExampleCost}	
	\caption{Comparison of the optimal total annualized cost for each benchmark step found by the fMINLP{} (orange) and the pNLP{} implementation (blue). To check if the fMINLP{} optimum indeed consistently found a cheaper local optimum, the pNLP{} optimization is run using the fMINLP{} optima as an initialization (green). The pNLP{} did not converge for case $n=713$ while the fMINLP{} did not converge for every case $n\geq663$. All non-converged optimizations are excluded from the figure and subsequent analysis}
	\label{fig:simpleCost}
\end{figure}

It can be seen that in this case, the fMINLP{} implementation (Figure \ref{fig:simpleCost} red) consistently finds a cheaper \gls{dhn} design than the pNLP{} implementation (Figure \ref{fig:simpleCost} blue). To check if this optimal design found by the fMINLP{} is indeed a better local optimum, the optimization is repeated in the pNLP{} and initialized with the optimal design found in the fMINLP{} (Figure \ref{fig:simpleCost} green). As can be seen in figure \ref{fig:simpleCost}, the pNLP{} remains in the optimum of the fMINLP{}, indicating that it is indeed a better local optimum that the pNLP{} was unable to find with the used initialization strategy. 

\begin{figure*}[h]
	\centering
	\includegraphics[width=0.736\linewidth]{fig04_TopologyComparison_figma}	
	\caption{Comparison of the optimal topology found with the fMINLP{} and the pNLP{} for the case of size $n = 113$ and the case of size $n = 463$. The optimal topologies found by the fMINLP{} use a shorter pipe length then the topologies found by the pNLP{}. Red lines represent the pipe network connecting the producer, heat consumers (black dots) and pipe junctions (grey dots). The line thickness indicates the pipe diameter while unused pipe connections are drawn in grey.}
	\label{fig:simpleTopo}
\end{figure*}

To better understand the difference in the optimal network designs causing this cost difference, the optimal network topologies found for network sizes $n = 63$ and $n = 463$ by both approaches are visualized in figure \ref{fig:simpleTopo}. The figure shows, that the optimal topologies of both approaches differ, with the pNLP{} using overall a longer pipe length. This difference in the optimal design can be explained with the different initialization methods of the optimization approaches. To facilitate convergence, the fMINLP{} implementation is initialized with an \gls{milp} solve, that effectively minimizes the overall used pipe length while satisfying the heat demands of the consumers. This initialization promotes topologies with minimal pipe lengths, that turns out to be a good optimal network topology for single producer networks. This optimal topology continues to be an optimum when later solving the full \gls{minlp} as well, as can be seen in figure \ref{fig:simpleTopo} b) and d). The pNLP{} on the other hand, initializing the full \gls{minlp} from a uniform distribution of pipes, finds a different local optimum using an overall larger pipe length (see figure \ref{fig:simpleTopo}). This knowledge of the influence of initializations on the optimal topology, can be used to reduce the proneness to local optima for both optimization approaches. In the following section \ref{sec:2producer} it will be studied how the initialization of both approaches perform on a multi producer case.

\section{The influence of initialization - a two producer case}\label{sec:2producer}
Modern $4^{\text{th}}$ generation \gls{dhn}s often feature multiple producers with different injection temperatures. To compare how both optimization approaches perform when designing such networks, a benchmark case was designed featuring two heat producers. 

\subsection*{Setup}
For this benchmark, a circular scalable test case with two producers with different injection temperatures is used. The layout of this case is shown in figure \ref{fig:hardSetup}. The circular network superstructure for this case is equivalent to the case in section \ref{chap:simple}, representing a potential district to be connected to a \gls{dhn}. Here, two producers are placed on the left and right of this district. The two producers represent different heat production technologies, with the left producer supplying heat at $\SI{70}{\degreeCelsius}$ at high costs ($\npvCost_{\mathrm{hC}} = \SI[per-mode=symbol]{800}{\sieuro\per\kilo\watt}$, $\npvCost_{\mathrm{hO}} = \SI[per-mode=symbol]{8}{\ct\per\kilo\watt\per\hour}$) while the producer on the right supplies heat at $\SI{55}{\degreeCelsius}$ at a lower cost ($\npvCost_{\mathrm{hC}} = \SI[per-mode=symbol]{0}{\sieuro\per\kilo\watt}$, $\npvCost_{\mathrm{hO}} = \SI[per-mode=symbol]{4}{\ct\per\kilo\watt\per\hour}$). The district contains houses with different heat system characteristics. While the houses on the top right, using a modern heating systems, can work with water at $\gtrapprox\SI{50}{\degreeCelsius}$, the rest of the houses require heat at higher temperatures ($\gtrapprox\SI{60}{\degreeCelsius}$) than the low temperature producer can provide. The size of this network can be increased by adding additional rings of houses to this network, representing different sizes of DHN development projects.

\begin{figure}[h]
	\includegraphics[width=1\columnwidth]{fig05_TwoProducerSetupTextwidth}	
	\caption{Setup for the two producer case. This setup features two producers on the left and right of a circular network. The network features old houses that require high temperature heat (red) and newer houses that can satisfy their demand with lower temperatures (green). The heat from the left producer (red) can supply all houses, while the right producer (green) can only supply the houses in the top right quadrant}
	\label{fig:hardSetup}
\end{figure}


\subsection*{Comparison of the optimized network designs}
Now this case is optimized using both the fMINLP{} and the pNLP{} implementation. To compare the two approaches, three different networks of increasing sizes where optimized (Case 1 with $n=138$, case 2 with $n = 298$ and case 3 with $n = 618$). The annualized costs of the resulting optima found by both approaches are visualized in figure \ref{fig:hardCost}.\footnote{Again the optima of both approaches where evaluated with a pNLP{} simulation to avoid discrepancies.}

\begin{figure*}[h]
	\centering
	\includegraphics[width=0.736\linewidth]{fig06_HarderCaseCost_h}	
	\caption{Comparison of the optimal total annualized cost for case 1-3 found by the fMINLP{} and the pNLP{} implementation. While the fMINLP{} finds optimal design of lower pipe cost (green), the pNLP{} finds optimal design that save on heating costs (red), ultimately leading to a lower total annualized cost for all three cases.}
	\label{fig:hardCost}
\end{figure*}

The cost comparison shows that in this two producer case, the pNLP{} found an optimal network design with a lower total annualized cost then the fMINLP{} in all three cases. While the optimal network design of the fMINLP{} is again cheaper in pipe investment cost, large savings in heat costs can be made with the proposed optimal design by the pNLP{}. While in case 2 the pipe investment cost of the optimal network design found by the fMINLP{} is $\SI{2.1}{\mega \sieuro}$ lower, the optimal design found by the pNLP{} achieves heat cost savings of $\SI{14.4}{\mega \sieuro}$, ultimately leading to a lower total annualized cost. 

For a better understanding of the difference in the optimal designs, the optimal topologies of both approaches are compared. In figure \ref{fig:hardTopo}, the optimal network topologies of case 1 and 2 optimized by both the pNLP{} and the fMINLP{} are visualized. It can be seen that the optimal topologies found by the pNLP{} connect most new houses to the low temperature producer on the top right, effectively creating two separate networks at different temperatures. This use of cheaper, low temperature heat from the producer on the right leads to the cheaper optimal network designs in comparison to the fMINLP{} that where observed in figure \ref{fig:hardCost}. The initialization steps of the fMINLP{}\footnote{Further described in \ref{sec:Apendix_init}} on the other hand, minimizes pipe length while satisfying the heat demands, therefore favoring an optimal topology for both cases that connects all houses to the high temperature producer on the left. The topology resulting from this initialization remains a local optimum also for the full \gls{minlp} optimization. While saving on pipe investment cost, this optimal network topology found by the fMINLP{}, is more expensive than the optimal topology by the pNLP{} because it relies on buying more high temperature expensive heat from the producer on the left.   

\begin{figure*}[h]
	\centering
	\includegraphics[width=0.736\linewidth]{fig07_HardCase_TopoComparison_figma}	
	\caption{Colored lines represent the pipe network connecting the producer, heat consumers (black dots) and pipe junctions (grey dots). The line thickness indicates the pipe diameter while unused pipe connections are drawn in grey. The line color corresponds to the water temperature in the pipe. While the pNLP{} finds an optimal topology connecting modern houses to the low temperature source, the fMINLP{} implementation favors a single network connecting all houses to the high temperature source.}
	\label{fig:hardTopo}
\end{figure*}

The results of the benchmarking study indicate that the initialization strategy for the topology optimization problem of \gls{dhn}s has a significant impact on the optimal designs found for both the single- and two-producer cases. This dependence on the initialization should always be taken into account when solving non-convex optimal design problems for DHNs that use non-linear models or cost functions. This sensitivity to the initialization could be mitigated by using globalization strategies, such as running multiple optimizations from different initializations, as demonstrated by Marty et al. \cite{marty2018simultaneous}.


 \section{Conclusion}
 In this paper, we compared the performance of two different approaches to the non-linear topology optimization of \gls{dhn}s. While the approach by Mertz et al. \cite{mertz2016minlp} solves a full \gls{minlp}, resolving the discrete nature of pipe routing choices using a combinatorial optimization approach (fMINLP{}), the approach by Wack et al \cite{wack2022} solves a relaxed \gls{nlp} ensuring a discrete network topology through penalization (pNLP{}). By performing a benchmark of both approaches, we showed that in the context of topology optimization of \gls{dhn}'s, solving the full \gls{minlp} leads to an exponential scaling of the computational cost with the network size. This steep cost scaling does not only make design studies slow and expensive, its exponential nature renders optimizations of large networks containing thousands of potential pipe connections intractable. On the other hand, solving the topology optimization problem as a relaxed \gls{nlp}, as is done by the pNLP{}, maintains a polynomial scaling of the computational cost. This polynomial scaling of the optimization time makes large-scale \gls{dhn} optimization feasible. We highlighted the scaling difference of both approaches by showing that for a \gls{dhn} featuring $n=2000$, an optimization time of $\walltime_{\mathrm{\gams{}}}\left(  n\right) \approx 43000 ~\mathrm{years}$ would be needed following the observed exponential trend of the \gls{minlp} solution time in the fMINLP{} implementation. The necessary optimization time reduces to $\walltime_{\mathrm{\pathopt{}}}\left(  n\right) \approx 5.3 ~\mathrm{hours}$ when using the relaxed \gls{nlp} implementation of the pNLP{}, assuming the observed polynomial trend.
 
 Posing the optimization problem as a non-linear optimization problem has the benefit of an accurate representation of the \gls{dhn} physics. Non-linearities in the optimization problem can lead to multiple local optima. We showed in this paper how the initialization of different optimization approaches can therefore lead to different optimal designs. For a \gls{dhn} optimization case containing a single producer, the fMINLP{} implementation consistently found cheaper optimal network designs than the pNLP{} approach. This difference in optimal design was explained by the different initialization strategies, notably the fMINLP{} implementation focusing on minimizing pipe length in the initialization. The resulting optimal topologies featuring minimal pipe lengths, prove beneficial for single producer networks. For a two producer case on the other hand, this initialization strategy showed to be a limitation and the pNLP{} approach consistently found cheaper optimal network topologies than the fMINLP{} implementation. The comparison of the optimality gap between both approaches shows that resolving the binary nature of the pipe routing problem by solving a full \gls{minlp} does not necessarily lead to better optimal \gls{dhn} topologies, while significantly increasing the computational cost. 
 
 The benchmark also showed that the proneness to local optima and the consequent importance of initialization strategies should be taken into account when deciding to solve the optimal design problem of \gls{dhn}s as a \gls{nlp}. This proneness to local optima can be mitigated with globalization strategies like e.g. running multiple optimizations from different initializations. Although used in literature, in this paper heuristic algorithms and linearized models were omitted. Future work should complement this comparison by benchmarking the scalability of said heuristics (e.g. genetic algorithms) and \gls{milp} approaches to the topology optimization of \gls{dhn}'s. 
 

 \section*{Data Availability}
 A data-set including the structure, input parameters and optimization results of the heating networks used in both benchmark cases of this paper is available at the following link: \url{https://doi.org/10.48804/DO1BRQ}. The optimization results can be replicated using the methodology and formulations described in this paper.
 
 Additionally a small case generator was written that can be used to recreate all benchmark cases of this paper. The repository can be found at the following link: \url{https://doi.org/10.5281/zenodo.7434451}
 
 
 \section*{Acknowledgements}
A special thank you to Valentin Gressot for his support in aligning the software base of both approaches for the benchmark.

Yannick Wack is funded by the Flemish institute for technological research
(VITO). 

\section*{CRediT authorship contribution statement}
\textbf{Yannick Wack:} Conceptualization, Formal analysis, Software, Visualization,
Writing – original draft. \textbf{Sylvain Serra:} Conceptualization, Formal analysis, Software,
Writing - Review \& Editing. \textbf{Martine Baelmans:} Conceptualization, Writing - Review \& Editing.  \textbf{Jean-Michel Reneaume:} Conceptualization, Writing - Review \& Editing.  \textbf{Maarten Blommaert:} Conceptualization, Software, Writing - Review \& Editing, Supervision. 

\section*{Compliance with ethical standards}
\subsubsection*{Conflict of interest}
The authors declare that they have no conflict of interest
\subsubsection*{Funding}
The authors did not receive support from any organization for the submitted
work.
\subsubsection*{Ethical approval}
This article does not contain any studies with human participants or animals
performed by any of the authors.


\section{Introduction}

\gls{dhn} are considered one of the core technologies to enable carbon-neutral space heating \cite{OECD/IEA2019}. It has the ability to connect a multitude of different renewable heat sources and provide heat to districts and entire cities. In \gls{dhn}s, the typically high upfront investment cost of groundworks and piping is a core decision variable for the feasibility of a development project. Therefore in the planning phase it is crucial to design the pipe routing (network topology), pipe sizing and heat production capacities in an optimal way. This topology optimization problem solves the question of where to place heating network pipes and at what diameter, while accounting for the non-linear nature of the heat transport physics. Together with the binary choice  of pipe placement it is therefore inherently a \gls{minlp}. 

Solving non-linear combinatorial problems directly can be challenging. As a result, the topology optimization problem of heating networks is often linearized, resulting in a linear optimization problem that can be solved efficiently using Mixed Integer Linear Programming (MILP) solvers.  S\"oderman \cite{Soderman2007}, for example, optimized the structure and configuration of a district cooling network using a MILP approach. Dorfner and Hamacher \cite{dorfner2014large} later used a linear approach to optimize the topology and pipe sizes of a DHN. Haikarainen et al. \cite{Haikarainen2014} optimized the topology and operations, while accounting for different production technologies and heat storage. Mazairac et al. \cite{Mazairac2015a} optimized the topology of a multi-carrier network that incorporates gas and electricity  supply. Morvay et al. \cite{Morvaj2016} optimally designed the network while also optimizing the energy mix supplied. Still using a MILP approach, Bordin et al. \cite{bordin2016optimization} optimized the network topology while studying the set of consumers to optimally connect to a DHN. Most recently, Resimont et al. \cite{Resimont2021} used a MILP approach to optimize a city-scale heating network. The transformation of modern \gls{dhn}s towards multi-source, low-temperature networks \cite{Lund2021} invalidates the assumptions of most linear heating network models. Accurate modeling of heat losses and the different feasible temperature levels of supply and demand sites in \gls{dhn}s requires a non-linear representation of the network physics. Solving the resulting \gls{minlp} problem is challenging. In the literature, one approach to solving this \gls{minlp} problem is to use heuristic search algorithms. For example, Li and Svendsen used a genetic algorithm \cite{Li2013}, and Allen et al. deployed both a minimal spanning tree heuristic and a particle swarm algorithm \cite{Allen2022} to the topology optimization of a small DHN.

To fully leverage the potential of mathematical optimization notably two approaches have been successfully applied to the non-linear optimal topology problem of heating networks. Mertz et al. \cite{mertz2016minlp}, solves the full \gls{minlp} resolving the discrete nature of pipe routing choices, using a combinatorial optimization approach. On the other hand, Blommaert et al. \cite{Blommaert2020a} and Wack et al. \cite{wack2022}, solve a relaxed \gls{nlp} ensuring a discrete network topology through penalization. This method is inspired by partial differential equation (PDE) constrained topology optimization, where penalization methods are combined with adjoint optimization to solve problems governed by PDEs. These methods promise to scale well with the problem size, allowing to optimize large-scale problems.

The current bottleneck to a widespread use of optimization methods for the design and topology optimization of \gls{dhn}s are an accurate representation of the network while maintaining scalability of the approach. The goal of this paper is therefore to compare the performance of the combinatorial approach by Mertz et al. \cite{mertz2016minlp} and the relaxed adjoint approach by Wack et al. \cite{wack2022} to non-linear topology optimization of \gls{dhn}s. First, the computational cost scaling of both optimization approaches with the network size is benchmarked to assess whether they are suited to be used for large-scale \gls{dhn} development projects. Both approaches solve non-linear and non-convex optimization problems, so the found optimal network designs can be different local optima. Therefore, in a second step, the found optimal designs are compared. Here, the influence of resolving the discrete nature of the pipe routing problem, as well as the influence of initialization strategies on the found local optima is discussed. This comparison is performed for a simple single heat producer case and, to closer mimic modern heating network problems, for a multi-producer case with different heat production temperatures.

\section{The topology optimization problem of \gls{dhn}s}\label{sec:OptProblem}
To compare both approaches, a topology optimization problem is defined that aligns the problem definitions of Mertz et al.\cite{mertz2016minlp} and Wack et al \cite{wack2022}. This ensures comparability of the results with both approaches. For completeness, the optimization problem definition is repeated here in short, for detailed discussion the reader is referred to the afore-mentioned publications.

In order to represent the design and topology optimization problem of a \gls{dhn}, a set of design variables $\defDesignVar$ is defined, containing the pipe diameters $\ve{\topVarSkalar}$ as well as the normalized producer inflows $\ve{\prodInput} $. To represent the physical state of a given network, a vector of physical variables $\defStateVar$ is defined. It contains the flow rates $\ve{\flow} $, nodal pressures $\ve{\pressure} $ and nodal and pipe exit temperatures $\ve{\dTinf} $. The temperatures $\ve{\dTinf} = \ve{T} - \temperature_\infty$ are defined as the difference between the absolute water temperature $\ve{T}$ and the outside air temperature $\temperature_\infty$. Now the topology optimization problem for a \gls{dhn} can be posed as a generic optimization problem of the form:

\begin{equation} \label{eq:opt}
	\begin{aligned} \min_{\ve{\designVarskal},\ve{x}}& \qquad
		\mathcal{J} \left(\ve{\designVarskal},\ve{x}\right) \\
		s.t.& \qquad \ve{h}(\ve{\designVarskal},\ve{x}) \leq 0.  \\
	\end{aligned}
\end{equation}

Here $\mathcal{J}$ represents the cost function and is defined as the total cost of the project over an investment horizon of $A = 30 ~\mathrm{years}$. $\ve{h}(\ve{\designVarskal},\ve{x})$ defines the set of model and technological constraints. A full definition of the cost function and model constraints for this benchmark can be found in \ref{app:model}. This optimization problem is based on Mertz et al. \cite{mertz2016minlp} and Wack et al. \cite{wack2022} and a detailed description of the underlying costs and models can also be found there. In the problem definition, the set definition in table \ref{tab:sets} is used to reference different parts of the network.

\begin{table}
	\centering
	\caption{Notation used to describe different subsets of the \gls{dhn}.}
	\label{tab:sets}
	\begin{tabular}{ll}
		\hline
		Set	& Symbol  \\
		\hline\hline
		All edges &  $E$ \\  
		Producer input & $\setEdges_\pro$ \\  
		Pipes& $\setEdges_\pipe$ \\
		Feed edges &  $\setEdges_{\mathrm{f}}$  \\
		Heating system & $\setEdges_\rad$  \\
		\hline
	\end{tabular}
\end{table}

\section{Methodology}
To compare the two approaches, the topology optimization problem (described in equation \ref{eq:optProblem}) is solved using two different methods: solving the full \gls{minlp} problem in GAMS using the implementation by Mertz et al. \cite{mertz2016minlp} (referred to as fMINLP{} for simplicity), and using the penalized adjoint optimization approach by Wack et al. \cite{wack2022} (referred to as pNLP{}). Both implementations were simplified to facilitate comparison. The pNLP{} implementation no longer accounts for multiple discrete pipe diameters, while the fMINLP{} implementation no longer optimizes the production side or cascading supply between high and low temperature consumers. The following sections briefly discuss the two methodologies and their different treatment of the binary pipe existence variable.

\subsection*{Combinatorial \gls{minlp} approach: fMINLP{}}
For the comparison, the topology optimization problem will be solved using a combinatorial \gls{minlp} approach. This approach was previously published by Mertz et al. \cite{mertz2016minlp}, and a comprehensive discussion of the approach can be found there. It is implemented in GAMS and uses the \gls{milp} solver CPLEX and the \gls{nlp} solver CONOPT. To pose the optimization problem as a combinatorial problem, the topological choice of pipe placement is represented by a binary variable $\ve{\phi}$. Physical variables defined on these pipes are then coupled to this existence variable using the bigM method ($\mathcal{M}$) \cite{mertz2016minlp}, as e.g. in the flow velocity definition: 

\begin{align}\label{eq:optProblem}
	&\gijEdge{m} -\mathcal{M}\left( 1-\gijEdge{\phi} \right)\leq\gijEdge{v} \density \pi \frac{1}{4}\gijEdge{d}^2\\ &\gijEdge{v} \density \pi \frac{1}{4}\gijEdge{d}^2\leq \gijEdge{m} +\mathcal{M}\left( 1-\gijEdge{\phi} \right) , \quad \forall \gi\gj \in \Epipe.
\end{align}
Here $\mathcal{M} \gg \gijEdge{m}$ is chosen in such a way that if a pipe exists $\gijEdge{\phi} = 1$, the bounds on the flow velocity constraint are tight. If $\gijEdge{\phi} = 0$, the bounds of this constraints are relaxed. To force the installation of both feed and return pipes the following constraint is used:
\begin{equation}
	\gijEdge{\phi}= \gEdge{\phi}{j}{i},\quad\forall \gi\gj \in \Epipe.
\end{equation}

The fMINLP{} implementation features additional constraints to facilitate convergence. They are defined in a way to not be active in the final optimal design in order to maintain the same final problem formulation as the pNLP{}. First, the installed capacity $H$ is bound by the sum of consumer demands: 

\begin{equation}
	\sum_{ij\in\setEdges_{\pro,\mathrm{f}}}\gijEdge{H} =\sum_{ij\in\setEdges_{\pro,\mathrm{f}}}\gijEdge{\left(q\dTinf  \right)} \density c_{\mathrm{p}} \leq 1.5 \sum_{kl\in \setEdges_\rad}\gEdge{Q}{k}{l}.
\end{equation}

To facilitate convergence on the momentum equations, the velocity and pressure drop over a pipe are bound by:

\begin{align}
	\gijEdge{v} &\leq 3.5, \\
	\gijEdge{v} &\leq 18.438 \gijEdge{d} + 0.2186, \\
	\gijEdge{v} &\geq v_{\mathrm{min}}\gijEdge{\phi},\\
	\giNode{p} - \gjNode{p}&\leq\SI{200}{\pascal} \quad \forall \gi\gj \in \Epipe.
\end{align}

This constraints are not active in the final optimized design, as they are dominated by the pipe momentum equation (\ref{eq:pipeMomentum}). A fixed pressure drop is set for the heat exchanger

\begin{equation}
	\giNode{p} - \gjNode{p}\geq\SI{2}{\kilo\pascal}, \quad \forall \gi\gj \in \Erad,
\end{equation}

and the exit temperature of the radiator is required to be smaller than the inflow temperature:
\begin{equation}
	\giNode{\dTinf} \geq \gjNode{\dTinf} , \quad \forall \gi\gj \in \Erad.	
\end{equation}

Finally, this optimization strategy uses a sequence of initializations to ensure a stable convergence of the final \gls{minlp} solve, based on the strategy of Marty et al. \cite{marty2018simultaneous}. For completeness, this initialization sequence is briefly described in \ref{sec:Apendix_init}.

\subsection*{Adjoint-based penalization approach: pNLP{}}

The other considered approach is the adjoint-based penalization approach by Wack et al. \cite{wack2022}. Here the combinatorial problem is relaxed, allowing for a continuous pipe placement choice. A near-discrete topology is then enforced by penalization, e.g. by replacing $p_0$ in the investment cost function (equation \ref{eq:Jpipe}) with
\begin{equation}
	\bar{p}_0(\gijEdge{d})= p_0 \left(\frac{1}{(1+\exp(-k\left(\gijEdge{d}-\diameter_{\mathrm{min}}\right))})-1\right)\,,
\end{equation}
or by explicitly penalizing the $d$ in the investment cost function and within the momentum and energy equations as described in Wack et al. \cite{wack2022}. 

Because of this relaxation, the \gls{minlp} reduces to solving a series of \gls{nlp}s. Here, the optimization is initialized from a uniform distribution of pipe diameters.

\section{Benchmarking computational cost - a single producer case}\label{chap:simple}

Modern \gls{dhn}s grow ever bigger and more complex, including multiple heat production sites at different temperatures, which optimal design tools need to be capable to handle. These tools therefore need to scale well with the heating network size in order to be applicable to relevant cases. In this first benchmark, the computational cost scaling of a combinatorial \gls{minlp} approach (fMINLP{}) is compared to solving a relaxed penalized \gls{nlp} (pNLP{}). This is done on a benchmark case containing a single producer. 

\subsection*{Setup}

To be able to compare the computational cost scaling of both optimization approaches with increasing network size an easily scalable \gls{dhn} optimization case is set up. This benchmark case is visualized in figure \ref{fig:simpleSetup}. 

\begin{figure}[h]
	\includegraphics[width=1\columnwidth]{fig01_SingleProducerSetupTextwidth}	
	\caption{Setup of the first benchmark case. Around a central producer, heat consumers (houses) are arranged in a circular pattern. They are connected by possible pipe routes (black lines). Pipe junctions are visualized as black circles. To increase the size of this case, additional segments $s$ are successively added to the outside of the network.}
	\label{fig:simpleSetup}
\end{figure}

Here, a heating network is to be designed around a single heat producer in the center of the network. This producer provides heat at $\SI{70}{\degreeCelsius}$. Around this central producer, houses and possible pipe connections are arranged in a circular manner. All houses have a heat demand of $Q = \SI{15}{\kilo\watt}$, and their heating system characterized by $\heaterCoef=\SI{200}{\watt\per\kelvin^\heaterExp}$ and $\heaterExp=1.2$. To investigate the cost scaling, the size of this circular network is successively increased by adding additional segments $s$. With each segment, additional heat consumers and 5 additional potential pipe connections are added, with the number of potential pipes following a linear scaling: $n\left( s \right) = 5s+13$. In figure \ref{fig:simpleSetup} the addition of segments $s = 10$ and $s=11$ are visualized in red and blue.

Now the benchmark of both approaches is conducted by creating a series of heating networks of increasing size. Segments $s$ are added in steps of 10, starting from $s=0$ up to $s=190$ therefore creating a sequence of heating networks with $(s_k)_{k=0}^{19}, ~s_{k}=10k$ segments. This sequence of optimization problems is solved by both the fMINLP{} implementation as well as the pNLP{}. This optimizations were performed on the same computer (Using a single Intel Xeon 3.20 GHz processor core). 

\subsection*{Comparison of the computational cost}
The optimization was repeated 3 times for each network and the mean runtime until convergence (wall time) was recorded. The wall time for each network size is visualized in figure \ref{fig:simpleScaling} on a a) semi-log graph and a b) log-log graph.

\begin{figure*}[h]
	\centering
	\includegraphics[width=0.736\linewidth]{fig02_SimpleExampleScaling}	
	\caption{Wall time scaling of the two approaches with increasing network size $n$. The scaling of the approach in the fMINLP{} is described well by an exponential fit a), while the scaling of the pNLP{} can be described well by a power fit b). The pNLP{} did not converge for case $n=713$ while the fMINLP{} did not converge for every case $n\geq663$. All non-converged optimizations are excluded from the figure and subsequent analysis}
	\label{fig:simpleScaling}
\end{figure*}

The figures show that the wall time for the fMINLP{} implementation scales faster with the network size then for the pNLP{}. In order to get an understanding of how long it would take both approaches to optimize networks of considerable size, it is useful to understand what function governs their cost scaling. Considering the algorithmic differences of solving a \gls{minlp} in the case of the fMINLP{} vs solving a relaxed \gls{nlp} in the case of the pNLP{}, both an exponential and a polynomial time scaling was tested on both approaches. First, an exponential fit was performed for the scaling of both approaches (see figure \ref{fig:simpleScaling} a)). The wall time scaling of the fMINLP{} can here be described by the function $\walltime_{\mathrm{\gams{}}}\left(  n\right) = 51.43e^{0.012n} \si{\second}$ with an coefficient of determination $R^2=0.87$, while the wall time scaling of the pNLP{} can be described by $\walltime_{\mathrm{\pathopt{}}}\left(  n\right) = 47.31e^{0.006n}\si{\second}$ with a coefficient of determination of $R^2 = -0.01$. From this fit, it can be concluded that the wall time scaling of the fMINLP{} approach can be described well by an exponential function, while the wall time scaling of the pNLP{} approach cannot. Second, a power function was fitted to the wall time scaling of both approaches. Here, the scaling of the pNLP{} can be described with $\walltime_{\mathrm{\pathopt{}}}\left(  n\right) = 0.04n^{1.72}$ with a coefficient of determination of $R^2=0.85$ and the scaling of the fMINLP{} can be described by $\walltime_{\mathrm{\gams{}}}\left(  n\right) = 0.16n^{1.75}$ with $R^2=0.04$. It can be concluded that the wall time scaling of the pNLP{} follows a polynomial function reasonably well, while the fMINLP{} cannot be described with a polynomial function. 

The comparison above highlights the potential exponential solution time of solving full \gls{minlp}s. In the context of \gls{dhn} optimization, this steep cost scaling does not only make design studies slow and expensive, its exponential nature renders optimizations of considerably sized networks intractable and serves as a bottleneck for the optimization of large-scale \gls{dhn}s, containing thousands of potential pipe connections. For example to optimize a network containing $n = 2000$ potential pipe connections, an optimization time of $\walltime_{\mathrm{\gams{}}}\left(  n\right) \approx 43000 ~\mathrm{years}$ would be needed following the observed exponential trend of the \gls{minlp} solution time.
In the context of topology optimization of \gls{dhn}s it is computationally favorable to solve a relaxed \gls{nlp}, as is evident from the polynomial wall time scaling of the pNLP{} implementation. The polynomial scaling of the computational solution time with this approach makes large-scale \gls{dhn} optimization feasible. This is highlighted for the above mentioned example of a network of $n = 2000$, where the necessary optimization time reduces to $\walltime_{\mathrm{\pathopt{}}}\left(  n\right) \approx 5.3 ~\mathrm{hours}$ when using the relaxed \gls{nlp} implementation of the pNLP{}, assuming the observed polynomial trend.

The potentially exponential computational cost scaling of solving full \gls{minlp}s should be taken into account when choosing optimization strategies for \gls{dhn} topologies. It will make the optimization of large-scale \gls{dhn}s intractable and hinder the application of such optimization tools to real-world \gls{dhn} problems. It is therefore advisable to simplify the \gls{minlp} formulation by either linearizing the non-linear model constraints, posing the problem as a \gls{milp}, or by relaxing the integer constraints, reformulating the problem as a \gls{nlp}. This second relaxation is done within the pNLP{} implementation, and while keeping the non-linear model in the overall optimization procedure, its significant speed up in the computational cost is shown in this paper.  

\subsection*{Comparison of found optimal network topologies}

Both approaches solve a non-linear optimization problem. This is necessary because the non-linear heat transport physics can have a big influence on the network design. Non-linearities in the optimization problem, however, can lead to multiple local optima. The proneness of both optimization approaches to these local optima is therefore studied in a next step. Here it can be investigated, if the relaxation of the integer constraints of pipe placement by the pNLP{} has a relevant influence on the cost of the found optimal network design. For this, the optimality gap of both approaches is compared. The optimal total annualized cost for each benchmark step found with both approaches is visualized in figure \ref{fig:simpleCost}.\footnote{To prevent potential remaining modelling differences from skewing the comparison, both found optimal designs where evaluated in a simulation within the pNLP{}.}

\begin{figure}[h]
	\includegraphics[width=1\columnwidth]{fig03_SimpleExampleCost}	
	\caption{Comparison of the optimal total annualized cost for each benchmark step found by the fMINLP{} (orange) and the pNLP{} implementation (blue). To check if the fMINLP{} optimum indeed consistently found a cheaper local optimum, the pNLP{} optimization is run using the fMINLP{} optima as an initialization (green). The pNLP{} did not converge for case $n=713$ while the fMINLP{} did not converge for every case $n\geq663$. All non-converged optimizations are excluded from the figure and subsequent analysis}
	\label{fig:simpleCost}
\end{figure}

It can be seen that in this case, the fMINLP{} implementation (Figure \ref{fig:simpleCost} red) consistently finds a cheaper \gls{dhn} design than the pNLP{} implementation (Figure \ref{fig:simpleCost} blue). To check if this optimal design found by the fMINLP{} is indeed a better local optimum, the optimization is repeated in the pNLP{} and initialized with the optimal design found in the fMINLP{} (Figure \ref{fig:simpleCost} green). As can be seen in figure \ref{fig:simpleCost}, the pNLP{} remains in the optimum of the fMINLP{}, indicating that it is indeed a better local optimum that the pNLP{} was unable to find with the used initialization strategy. 

\begin{figure*}[h]
	\centering
	\includegraphics[width=0.736\linewidth]{fig04_TopologyComparison_figma}	
	\caption{Comparison of the optimal topology found with the fMINLP{} and the pNLP{} for the case of size $n = 113$ and the case of size $n = 463$. The optimal topologies found by the fMINLP{} use a shorter pipe length then the topologies found by the pNLP{}. Red lines represent the pipe network connecting the producer, heat consumers (black dots) and pipe junctions (grey dots). The line thickness indicates the pipe diameter while unused pipe connections are drawn in grey.}
	\label{fig:simpleTopo}
\end{figure*}

To better understand the difference in the optimal network designs causing this cost difference, the optimal network topologies found for network sizes $n = 63$ and $n = 463$ by both approaches are visualized in figure \ref{fig:simpleTopo}. The figure shows, that the optimal topologies of both approaches differ, with the pNLP{} using overall a longer pipe length. This difference in the optimal design can be explained with the different initialization methods of the optimization approaches. To facilitate convergence, the fMINLP{} implementation is initialized with an \gls{milp} solve, that effectively minimizes the overall used pipe length while satisfying the heat demands of the consumers. This initialization promotes topologies with minimal pipe lengths, that turns out to be a good optimal network topology for single producer networks. This optimal topology continues to be an optimum when later solving the full \gls{minlp} as well, as can be seen in figure \ref{fig:simpleTopo} b) and d). The pNLP{} on the other hand, initializing the full \gls{minlp} from a uniform distribution of pipes, finds a different local optimum using an overall larger pipe length (see figure \ref{fig:simpleTopo}). This knowledge of the influence of initializations on the optimal topology, can be used to reduce the proneness to local optima for both optimization approaches. In the following section \ref{sec:2producer} it will be studied how the initialization of both approaches perform on a multi producer case.

\section{The influence of initialization - a two producer case}\label{sec:2producer}
Modern $4^{\text{th}}$ generation \gls{dhn}s often feature multiple producers with different injection temperatures. To compare how both optimization approaches perform when designing such networks, a benchmark case was designed featuring two heat producers. 

\subsection*{Setup}
For this benchmark, a circular scalable test case with two producers with different injection temperatures is used. The layout of this case is shown in figure \ref{fig:hardSetup}. The circular network superstructure for this case is equivalent to the case in section \ref{chap:simple}, representing a potential district to be connected to a \gls{dhn}. Here, two producers are placed on the left and right of this district. The two producers represent different heat production technologies, with the left producer supplying heat at $\SI{70}{\degreeCelsius}$ at high costs ($\npvCost_{\mathrm{hC}} = \SI[per-mode=symbol]{800}{\sieuro\per\kilo\watt}$, $\npvCost_{\mathrm{hO}} = \SI[per-mode=symbol]{8}{\ct\per\kilo\watt\per\hour}$) while the producer on the right supplies heat at $\SI{55}{\degreeCelsius}$ at a lower cost ($\npvCost_{\mathrm{hC}} = \SI[per-mode=symbol]{0}{\sieuro\per\kilo\watt}$, $\npvCost_{\mathrm{hO}} = \SI[per-mode=symbol]{4}{\ct\per\kilo\watt\per\hour}$). The district contains houses with different heat system characteristics. While the houses on the top right, using a modern heating systems, can work with water at $\gtrapprox\SI{50}{\degreeCelsius}$, the rest of the houses require heat at higher temperatures ($\gtrapprox\SI{60}{\degreeCelsius}$) than the low temperature producer can provide. The size of this network can be increased by adding additional rings of houses to this network, representing different sizes of DHN development projects.

\begin{figure}[h]
	\includegraphics[width=1\columnwidth]{fig05_TwoProducerSetupTextwidth}	
	\caption{Setup for the two producer case. This setup features two producers on the left and right of a circular network. The network features old houses that require high temperature heat (red) and newer houses that can satisfy their demand with lower temperatures (green). The heat from the left producer (red) can supply all houses, while the right producer (green) can only supply the houses in the top right quadrant}
	\label{fig:hardSetup}
\end{figure}


\subsection*{Comparison of the optimized network designs}
Now this case is optimized using both the fMINLP{} and the pNLP{} implementation. To compare the two approaches, three different networks of increasing sizes where optimized (Case 1 with $n=138$, case 2 with $n = 298$ and case 3 with $n = 618$). The annualized costs of the resulting optima found by both approaches are visualized in figure \ref{fig:hardCost}.\footnote{Again the optima of both approaches where evaluated with a pNLP{} simulation to avoid discrepancies.}

\begin{figure*}[h]
	\centering
	\includegraphics[width=0.736\linewidth]{fig06_HarderCaseCost_h}	
	\caption{Comparison of the optimal total annualized cost for case 1-3 found by the fMINLP{} and the pNLP{} implementation. While the fMINLP{} finds optimal design of lower pipe cost (green), the pNLP{} finds optimal design that save on heating costs (red), ultimately leading to a lower total annualized cost for all three cases.}
	\label{fig:hardCost}
\end{figure*}

The cost comparison shows that in this two producer case, the pNLP{} found an optimal network design with a lower total annualized cost then the fMINLP{} in all three cases. While the optimal network design of the fMINLP{} is again cheaper in pipe investment cost, large savings in heat costs can be made with the proposed optimal design by the pNLP{}. While in case 2 the pipe investment cost of the optimal network design found by the fMINLP{} is $\SI{2.1}{\mega \sieuro}$ lower, the optimal design found by the pNLP{} achieves heat cost savings of $\SI{14.4}{\mega \sieuro}$, ultimately leading to a lower total annualized cost. 

For a better understanding of the difference in the optimal designs, the optimal topologies of both approaches are compared. In figure \ref{fig:hardTopo}, the optimal network topologies of case 1 and 2 optimized by both the pNLP{} and the fMINLP{} are visualized. It can be seen that the optimal topologies found by the pNLP{} connect most new houses to the low temperature producer on the top right, effectively creating two separate networks at different temperatures. This use of cheaper, low temperature heat from the producer on the right leads to the cheaper optimal network designs in comparison to the fMINLP{} that where observed in figure \ref{fig:hardCost}. The initialization steps of the fMINLP{}\footnote{Further described in \ref{sec:Apendix_init}} on the other hand, minimizes pipe length while satisfying the heat demands, therefore favoring an optimal topology for both cases that connects all houses to the high temperature producer on the left. The topology resulting from this initialization remains a local optimum also for the full \gls{minlp} optimization. While saving on pipe investment cost, this optimal network topology found by the fMINLP{}, is more expensive than the optimal topology by the pNLP{} because it relies on buying more high temperature expensive heat from the producer on the left.   

\begin{figure*}[h]
	\centering
	\includegraphics[width=0.736\linewidth]{fig07_HardCase_TopoComparison_figma}	
	\caption{Colored lines represent the pipe network connecting the producer, heat consumers (black dots) and pipe junctions (grey dots). The line thickness indicates the pipe diameter while unused pipe connections are drawn in grey. The line color corresponds to the water temperature in the pipe. While the pNLP{} finds an optimal topology connecting modern houses to the low temperature source, the fMINLP{} implementation favors a single network connecting all houses to the high temperature source.}
	\label{fig:hardTopo}
\end{figure*}

The results of the benchmarking study indicate that the initialization strategy for the topology optimization problem of \gls{dhn}s has a significant impact on the optimal designs found for both the single- and two-producer cases. This dependence on the initialization should always be taken into account when solving non-convex optimal design problems for DHNs that use non-linear models or cost functions. This sensitivity to the initialization could be mitigated by using globalization strategies, such as running multiple optimizations from different initializations, as demonstrated by Marty et al. \cite{marty2018simultaneous}.


 \section{Conclusion}
 In this paper, we compared the performance of two different approaches to the non-linear topology optimization of \gls{dhn}s. While the approach by Mertz et al. \cite{mertz2016minlp} solves a full \gls{minlp}, resolving the discrete nature of pipe routing choices using a combinatorial optimization approach (fMINLP{}), the approach by Wack et al \cite{wack2022} solves a relaxed \gls{nlp} ensuring a discrete network topology through penalization (pNLP{}). By performing a benchmark of both approaches, we showed that in the context of topology optimization of \gls{dhn}'s, solving the full \gls{minlp} leads to an exponential scaling of the computational cost with the network size. This steep cost scaling does not only make design studies slow and expensive, its exponential nature renders optimizations of large networks containing thousands of potential pipe connections intractable. On the other hand, solving the topology optimization problem as a relaxed \gls{nlp}, as is done by the pNLP{}, maintains a polynomial scaling of the computational cost. This polynomial scaling of the optimization time makes large-scale \gls{dhn} optimization feasible. We highlighted the scaling difference of both approaches by showing that for a \gls{dhn} featuring $n=2000$, an optimization time of $\walltime_{\mathrm{\gams{}}}\left(  n\right) \approx 43000 ~\mathrm{years}$ would be needed following the observed exponential trend of the \gls{minlp} solution time in the fMINLP{} implementation. The necessary optimization time reduces to $\walltime_{\mathrm{\pathopt{}}}\left(  n\right) \approx 5.3 ~\mathrm{hours}$ when using the relaxed \gls{nlp} implementation of the pNLP{}, assuming the observed polynomial trend.
 
 Posing the optimization problem as a non-linear optimization problem has the benefit of an accurate representation of the \gls{dhn} physics. Non-linearities in the optimization problem can lead to multiple local optima. We showed in this paper how the initialization of different optimization approaches can therefore lead to different optimal designs. For a \gls{dhn} optimization case containing a single producer, the fMINLP{} implementation consistently found cheaper optimal network designs than the pNLP{} approach. This difference in optimal design was explained by the different initialization strategies, notably the fMINLP{} implementation focusing on minimizing pipe length in the initialization. The resulting optimal topologies featuring minimal pipe lengths, prove beneficial for single producer networks. For a two producer case on the other hand, this initialization strategy showed to be a limitation and the pNLP{} approach consistently found cheaper optimal network topologies than the fMINLP{} implementation. The comparison of the optimality gap between both approaches shows that resolving the binary nature of the pipe routing problem by solving a full \gls{minlp} does not necessarily lead to better optimal \gls{dhn} topologies, while significantly increasing the computational cost. 
 
 The benchmark also showed that the proneness to local optima and the consequent importance of initialization strategies should be taken into account when deciding to solve the optimal design problem of \gls{dhn}s as a \gls{nlp}. This proneness to local optima can be mitigated with globalization strategies like e.g. running multiple optimizations from different initializations. Although used in literature, in this paper heuristic algorithms and linearized models were omitted. Future work should complement this comparison by benchmarking the scalability of said heuristics (e.g. genetic algorithms) and \gls{milp} approaches to the topology optimization of \gls{dhn}'s. 
 

 \section*{Data Availability}
 A data-set including the structure, input parameters and optimization results of the heating networks used in both benchmark cases of this paper is available at the following link: \url{https://doi.org/10.48804/DO1BRQ}. The optimization results can be replicated using the methodology and formulations described in this paper.
 
 Additionally a small case generator was written that can be used to recreate all benchmark cases of this paper. The repository can be found at the following link: \url{https://doi.org/10.5281/zenodo.7434451}
 
 
 \section*{Acknowledgements}
A special thank you to Valentin Gressot for his support in aligning the software base of both approaches for the benchmark.

Yannick Wack is funded by the Flemish institute for technological research
(VITO). 

\section*{CRediT authorship contribution statement}
\textbf{Yannick Wack:} Conceptualization, Formal analysis, Software, Visualization,
Writing – original draft. \textbf{Sylvain Serra:} Conceptualization, Formal analysis, Software,
Writing - Review \& Editing. \textbf{Martine Baelmans:} Conceptualization, Writing - Review \& Editing.  \textbf{Jean-Michel Reneaume:} Conceptualization, Writing - Review \& Editing.  \textbf{Maarten Blommaert:} Conceptualization, Software, Writing - Review \& Editing, Supervision. 

\section*{Compliance with ethical standards}
\subsubsection*{Conflict of interest}
The authors declare that they have no conflict of interest
\subsubsection*{Funding}
The authors did not receive support from any organization for the submitted
work.
\subsubsection*{Ethical approval}
This article does not contain any studies with human participants or animals
performed by any of the authors.

