\section{Introduction} \label{sec:introduction}

Model Order Reduction (MOR) is an indispensable technology in the context of multi-query, real-time or slim computing context,
as full oder models (FOMs) of real  physical processes may be unaffordable.
Many techniques for MOR exist, cf.\ \cite{benner2020model}.
A large class of MOR methods are reduced basis (RB) models, which are projection-based and the projection spaces are typically constructed from samples of the FOM and allow rigorous a-posteriori error estimation by residual analysis~\cite{haasdonk2017reduced}.
As a drawback, the time-integration requires time-marching which may be a bottleneck in case of a fine time resolution.
In recent years, also machine learning (ML) based techniques have entered the field, which enable rapid prediction without expensive time-marching \cite{wang2019nonintrusive}.
Such data-based approaches, however, mostly lack any error guarantees.
Recently, we presented an adaptive and hierarchical RB-ML-ROM framework, which combines the benefits of those two approaches, i.e.~rapid prediction of ML-models including certification from RB methods \cite{haasdonk2022certified}.
In the current presentation, we demonstrate the versatility of the RB-ML-ROM framework, by including efficient multilayer kernel models for the ML part, namely Structured Deep Kernel Networks (SDKN) and two-layered greedy kernel models (2L-VKOGA).

The paper is organized as follows: In \Cref{subsec:rom} and \Cref{subsec:ml} we recall some more background information on reduced order modeling and machine learning.
\Cref{sec:deep_kernels} elaborates on novel deep kernel models and \Cref{sec:certified} reviews the recently introduced certified and adaptive RB-ML-MOR model.
Subsequently, \Cref{sec:num_appl} provides the numerical experiments on the combination of these approaches and finally \Cref{sec:conclusion} summarizes and gives an outlook.
\vspace{-.3cm}

\section{Reduced order modeling} \label{subsec:rom}
\vspace{-.3cm}

In this paper, we consider parametrized parabolic PDEs on a bounded Lipschitz-domain~$\Omega\subset\mathbb{R}^m$ over the time interval~$[0,T]$ with~$T>0$.
Further, let~$\mathcal{P}\subset\mathbb{R}^p$ be the set of parameters.
We denote by~$V_h$ a Hilbert space associated to a finite element discretization of~$\Omega$, such that~$V_h\subset H^1(\Omega)\subset L^2(\Omega)$.
For a parameter~$\mu\in\mathcal{P}$, we consider the semi-discrete problem of finding~$u_h(\mu) \in L^2(0,T; V_h)$ with $\partial_t u_h(\mu) \in L^2(0,T; V_h')$ such that
\begin{align}
	\langle\partial_t u_h(\mu), v_h\rangle + a\big(u_h(\mu), v_h; \mu\big) &= l(v_h; \mu) &&\forall \ v_h \in V_h \ \text{in } [0,T],\label{equ:fom-problem}\\
	u_h(0; \mu) &= u_0(\mu), &&\notag
\end{align}
where~$u_0(\mu)\in V_h$ is the initial condition.
In addition, we are particularly interested in output trajectories, i.e.~we consider an output functional~$f_h\colon\mathcal{P}\to L^2(0,T)$ defined as~$f_h(t;\mu)\coloneqq s(u_h(t;\mu))$ for some~$s\in V_h'$.
For every~$\mu\in\mathcal{P}$, $l(\,\cdot\,;\mu)\in V_h'$ is assumed to be a linear and continuous functional and~$a(\,\cdot\,,\,\cdot\,; \mu)\colon V_h \times V_h \to \mathbb{R}$ denotes a continuous and coercive bilinear form.
\par
Since solving the problem above for many different values of the parameter is costly, we are looking for a reduced order model that is faster to evaluate.
To this end, the reduced basis (RB) method projects the problem~\eqref{equ:fom-problem} onto a well-chosen reduced space~$V_\text{RB}\subset V_h$ with~$\dim V_\text{RB}\eqqcolon N_\text{RB}\ll N_h\coloneqq\dim V_h$ (for a general introduction to RB methods, see~\cite{haasdonk2017reduced}).
One thus searches for an approximate solution~$u_\text{RB}(t;\mu)\in V_\text{RB}$ of~$u_h(t;\mu)\in V_h$ by solving the same semi-discrete problem as in~\eqref{equ:fom-problem} with test functions from the space~$V_\text{RB}$, i.e.
\begin{align}
	\langle\partial_t u_\text{RB}(\mu), v_\text{RB}\rangle\hspace*{-2pt} +\hspace*{-1pt} a\big(u_\text{RB}(\mu), v_\text{RB}; \mu\big) &= l(v_\text{RB}; \mu) &&\forall \ v_\text{RB} \in V_\text{RB}  \  \text{in } [0,T] ,\label{equ:rom-problem}\\
	\langle u_\text{RB}(0;\mu),v_\text{RB}\rangle &= \langle u_0(\mu),v_\text{RB}\rangle \hspace*{-5pt}&&\forall \ v_\text{RB} \in V_\text{RB}.\notag
\end{align}
To obtain an efficient reduced order model, we pose the standard assumption that the bilinear form~$a$ and the right hand side~$l$ are parameter separable.
In that case, many of the required quantities can be precomputed during the so-called \emph{offline}-phase.
In the subsequent \emph{online}-phase, the reduced solution for a new parameter can be determined independently of the dimension~$N_h$ of~$V_h$.
\par
To assess the error between~$u_h$ and~$u_\text{RB}$, one can consider the residual by inserting~$u_\text{RB}$ into Eq.~\eqref{equ:fom-problem}.
The norm of the residual can be computed without actually solving~\eqref{equ:fom-problem}.
Further, the corresponding error estimator~$E_\text{RB}(u;\mu)$ for arbitrary~$u\in V_\text{RB}$, introduced in~\cite{grepl2005aposteriori}, can be evaluated efficiently during the online phase if~$a$ and~$l$ fulfill the parameter separability assumption.
Since the functional~$s$ is assumed to be linear and continuous, the error estimator for the state~$u_\text{RB}$ also results in an estimator for the error in the output functional, i.e.
\[
\lVert f_h(\mu)-s(u_\text{RB}(\mu))\rVert_{L^2(0,T)}\leq\lVert s(\,\cdot\,;\mu)\rVert_{V_h'}E_\text{RB}(u_\text{RB}(\mu);\mu).
\]
We emphasize at this point that the error estimator computes an estimate based on the parameter~$\mu$ and coefficients with respect to the reduced basis.
However, these coefficients do not have to be the result of solving the problem in~\eqref{equ:rom-problem}.
\vspace{-.3cm}

\section{Machine learning for regression} \label{subsec:ml}
\vspace{-.3cm}

In the case of regression, ML procedures typically assume that input data $X_N \coloneqq \{x_1, \dots, x_N \} \subset \mathbb{R}^d$ and corresponding output data $Y_N \coloneqq \{ y_1, \dots, y_N \} \subset \mathbb{R}^b$ is given.
The task is to ''learn`` a function which approximates the mapping $x_i \mapsto y_i$ and generalizes properly to unseen novel inputs $x$.
One mathematically well established approach for this task is to employ kernel methods.
These kind of methods can be understood as applying linear regression to a nonlinear transformation of the data, which is implicitly given by the kernel~$k$, which is in our case a symmetric function $k\colon \mathbb{R}^d \times \mathbb{R}^d \rightarrow \mathbb{R}$.
A well known representer theorem then states that $f^*$ can be expanded as
\begin{align}
	\label{eq:kernel_model}
	s_{X_N} = \sum_{j=1}^N \alpha_j k(\cdot, x_j).
\end{align}
In order to obtain sparse models, greedy methods can be employed to select a suitable subset of centers $X_n \subset X_N$ such that the kernel expansion from Eq.~\eqref{eq:kernel_model} is smaller.
An algorithm for this, which will be leveraged later on, is the vectorial kernel orthogonal greedy algorithm (VKOGA) \cite{santin2021kernel}. \\
Another lass of ML methods is deep learning, where probably neural networks are the most prominent example.
Due to their deep structure, feature learning is possible which can be beneficial especially for medium to high dimensional datasets. \\
In order to make the best out of both approaches, the idea is to combine deep learning methods such as neural networks with the benefits of ''shallow`` kernel methods.
In the following \Cref{sec:deep_kernels}, we review two recently introduced algorithms in this direction:
SDKNs \cite{wenzel2021universality} as well as 2L-VKOGA \cite{wenzel2023data}.
\vspace{-.3cm}

\section{Deep kernel models} \label{sec:deep_kernels}
\vspace{-.3cm}
Both, the SDKN and the 2L-VKOGA approach, are based on a deep kernel representer theorem \cite{bohn2019representer},
which allows to put kernel methods into a multilayer setup in a mathematical principled way: \\
\textbf{SDKNs} \cite{wenzel2021universality} make use of special classes of kernels, which are used alternately in a structured way to derive a powerful setup:
First, matrix-valued linear kernels are used, i.e.\ $k_\text{lin}(x, z) = \langle x, z \rangle_{\mathbb{R}^{d_0}} \cdot I_{d_1}$ with $I_{d_1} \in \mathbb{R}^{d_1 \times d_1}$ being the identity matrix.
Second, single-dimensional kernels are employed that act on single components $x^{(i)}, z^{(i)}$ separately as $k_s(x, z) = \diag(k(x^{(1)}, z^{(1)}), \dots, k(x^{(d)}, z^{(d)}))$.
Hereby, the one dimensional kernel $k\colon \mathbb{R} \times \mathbb{R} \rightarrow \mathbb{R}$ is e.g.\ the Gaussian kernel.
The SDKN model is set up using the linear kernels for every odd layer and single-dimensional kernels for every even layer.
Thus, the odd layers essentially allow for linear combinations of the inputs, while the even layers can be compared to activation functions of neural networks.
However, due to being kernel models that can be optimized instead of fixed activation functions, they allow for more flexibility and thus enable a potentially faster optimization. %
As analyzed in \cite{wenzel2021universality}, these SDKNs enjoy universal approximation properties in various limit cases.\\
The \textbf{2L-VKOGA} \cite{wenzel2023data} combines the use of an optimized two-layered kernel with greedy kernel algorithms, where we leverage the VKOGA algorithm \cite{santin2021kernel}: %
Radial basis function kernels are a popular class of kernels which are given as $k(x, z) = \Phi(\Vert x - z \Vert)$
for some radial basis function $\Phi\colon \mathbb{R}^d \rightarrow \mathbb{R}$,
e.g.\ the Gaussian kernel $k(x, z) = \exp(-\Vert x - z \Vert^2)$.
These RBF kernels usually make use of a shape parameter $\varepsilon > 0$ via $k_\varepsilon(x, z) = \Phi(\varepsilon \Vert x - z \Vert)$,
which can be used to adopt the shape of the kernel to the dataset. %
The two-layered kernel approach generalizes this scale parameter by first allowing for a $d \times d$ matrix $A$ instead of only a scalar valued shape parameter~$\varepsilon$, i.e.\ $k_A(x, z) \coloneqq \Phi(\Vert A(x-z) \Vert)$,
and secondly providing an ML-based gradient descent optimization, which allows to optimize the kernel (i.e.\ the matrix $A$) efficiently.
In particular, this setup can be framed within the deep kernel representer theorem \cite{bohn2019representer} and thus be seen as a two-layered kernel. %
The optimized kernel $k_A$ can subsequently be used with standard kernel algorithms.
Here we use it in conjunction with greedy algorithms, and in particular with the VKOGA algorithm:
In view of Eq.~\eqref{eq:kernel_model}, the VKOGA algorithm selects a suitable subset $X_n \subset X_N$ of the dataset, such that the resulting kernel model consists only of $n \ll N$ summands and is thus sparse. %
For our current application, we extend the 2L-VKOGA approach that was only introduced for scalar valued target data \cite{wenzel2023data} to vectorial data, i.e. $Y_N \subset \mathbb{R}^b$.
Therefore, we consider for the kernel optimization the $\ell^2$-norm over the errors from the single dimensions.
As the VKOGA algorithm is intrinsically already vector-valued, both can be combined.

\vspace{-.3cm}
\section{Certified and adaptive RB-ML-ROM} \label{sec:certified}
\vspace{-.3cm}

We summarize the adaptive model hierarchy as specified in~\cite{haasdonk2022certified}.
Instead of describing the model hierarchy in an abstract way as done in~\cite{haasdonk2022certified}, we make direct use of the RB model from \Cref{subsec:rom} and the deep kernel methods from \Cref{sec:deep_kernels}.
\par
The RB-ML-ROM is built as a hierarchy consisting of three models, a full-order model (FOM), a reduced basis reduced order model (RB-ROM), and a machine learning reduced order model (ML-ROM).
As discussed, the error estimator of the RB-ROM (see \Cref{subsec:rom}) can be used to certify solutions from the RB-ROM as well as from the ML-ROM.
The guiding assumption is that the ML-ROM is the fastest model in the hierarchy.
A query to the RB-ML-ROM now proceeds as follows:
For a new parameter the ML-ROM is queried and its prediction accepted if its error estimation is sufficiently accurate, i.e.~the error is below a desired tolerance~$\varepsilon>0$.
Otherwise, the slightly more expensive RB-ROM is queried and its error is again validated using the error estimator.
Only if the RB-ROM is also inaccurate, one has to fall back to the costly FOM, which is assumed to be accurate enough for every parameter.
\par
Not only the prediction of the hierarchical model is adaptive, but the submodels themselves are adaptively updated as follows: When querying the adaptive RB-ML-ROM for the first time for a certain parameter~$\mu\in\mathcal{P}$, the RB-ROM and the ML-ROM are initialized as empty models, i.e.~they can only return zero as the solution.
Consequently, most likely, the FOM is called for the parameter~$\mu$ and the FOM output trajectory~$f_h(\mu)$ is returned.
As intermediate quantity the FOM solution trajectory~$u_h(\mu)$ is computed, which will be used to set up the RB-ROM.
In our concrete setting, the RB-ROM is constructed using the hierarchical approximate proper orthogonal decomposition (HaPOD)~\cite{himpe2018hierarchical} applied to data obtained from querying the FOM.
Whenever solving the FOM is necessary (due to the RB-ROM and ML-ROM not being accurately enough for the current parameter), a new FOM trajectory becomes available and the RB-ROM can be improved using this new information.
Finally, the ML-ROM is built upon the updated RB-ROM and is based on the same reduced basis and error estimator as the RB-ROM.
More precisely, using an ML algorithm such as the deep kernel models described in \Cref{sec:deep_kernels}, we learn the coefficients of the RB solution with respect to the reduced basis as a function of the parameter, see also~\cite{wang2019nonintrusive}.
The required training data is generated as a byproduct when the RB-ROM is queried.
However, only those RB solutions that fulfill the accuracy requirements (evaluated by the error estimator of the RB-ROM) are added to the training set.
Recall that the error estimator of the RB-ROM takes as input a parameter~$\mu\in\mathcal{P}$ and arbitrary reduced coefficients.
In particular, these coefficients do not have to stem from the RB-ROM, but can also be the result of an ML prediction.
This fact enables the direct application of the error estimator to certify the ML-ROM independent of the actual ML model used.
\par
For the ML-ROM, two different approaches are possible which differ in the way they treat the time component.
The \emph{random-access-in-time} method passes the time instance along with the parameter as an input to the machine learning model.
To obtain the whole solution trajectory, the machine learning model is therefore queried for all different time instances.
In contrast, the \emph{time-vectorized} approach returns the whole time sequence for a given parameter, i.e.~the output contains the reduced coefficients for all time steps at once.
\vspace{-.3cm}

\section{Numerical application} \label{sec:num_appl}

\vspace{-.3cm}
As an application for the deep kernel methods in the framework of the adaptive RB-ML-ROM model, we leverage the Monte Carlo example introduced in \cite[Section 4.2]{haasdonk2022certified}.
Likewise, we can also use the NN used there as a baseline method for comparison.
The quantity of interest is the average (in time) temperature~$\overline{f}\colon \mathcal{P} \rightarrow \mathbb{R}$ within a given room of a building, represented by the domain $\Omega$.
The temperature can be computed by solving the heat equation and is depending on 28 diffusion parameters for walls, doors and heaters.
In order to determine the average and variance of this output quantity over a given probability density~$\rho$ on the 28 dimensional parameter space, we rely on a Monte Carlo estimation
\begin{align*}
	\text{E}[\overline{f}] \eqqcolon \int_\mathcal{P} \rho(\mu)\overline{f}(\mu) ~ \mathrm{d}\mu \approx \frac{1}{N_\text{mc}} \sum_{\mu \in \mathcal{P}_\text{mc}} \overline{f}(\mu),
\end{align*}
and similarly for the variance. \\
For the numerical experiment, the same setup as in \cite{haasdonk2022certified} was used, i.e.\ the full order model consists of 321206 DoFs and 999 timesteps and is thus quite expensive.
Therefore, the use of the adaptive model introduced in \Cref{sec:certified} is desirable, which is applied with a tolerance of $\varepsilon = 5 \cdot 10^{-2}$.
For the machine learning models, the SDKN followed the setup of the NN, i.e.\ four layers with each an inner dimension of 128 while using the Gaussian kernel for its single-dimensional kernel layers.
The 2L-VKOGA used a quadratic Matérn kernel and 10 epochs with a batch size of 32 and an initial learning rate of $5 \cdot 10^{-3}$ for the Adam optimizer for the optimization of the two-layered kernel.
For the greedy algorithm, the so called $f$-greedy selection criterion was used up to an expansion size of 500 centers.
The NN and SDKN were (re)trained as soon as the data of 200 new parameter values were evaluated, the VKOGA and 2L-VKOGA models were (re)trained as soon as $40$ new data samples were available.
As soon as the performance of the ML-ROM models was sufficiently good (use of ML model at least in 60\% of the queried parameter values), no more retrainings were performed.
The results of the runs of the adaptive model are collected in \Cref{tab:overview_results} as well as \Cref{fig:adaptive_model_results}:
\Cref{tab:overview_results} summarizes the key performance numbers for the different machine learning models:
In the column ''Training ML-ROM``, the training set size $n_\mu$ and the training duration $t$ for the ML-ROMs are listed:
It can be seen that the VKOGA models are particularly fast, which is due to their smaller training batch sizes of $40$ instead of $200$, such that they are trained already earlier.
Of particular interest is the column ''ML-ROM ratio``, which shows the ratio of ML-ROM vs.\ ROM evaluations:
For the baseline NN, this ratio increases from $0\%$ (before training) to up to $87\%$ after $3$ training stages.
The shallow VKOGA model does not seem to be accurate at all at any time, its ML-ROM ratio never exceeds 25\% -- and is thus stopped after $2000$ queried parameter evaluations.
For the SDKN, the ML-ROM ratio is already $60.3\%$ after the first training,
thus surpassing the NN performance after two trainings.
Its final ratio of $94.7\%$ is large enough such that no further retraining is required anymore.
The deep 2L-VKOGA model achieves a performance of $86.7\%$ already after the very first training, which included only $40$ evaluated parameter values.
The reason for this superb performance (especially in contrast to the shallow VKOGA model) can be seen in the optimized kernel.
A singular value analysis of the matrix~$A$ of the kernel $k_A$ reveals that (quite independent of the specific run,
i.e.\ random choice of the parameter $\mu$) there is one large singular value of $2.5$, while the second largest one is already only $0.35$ and the remaining ones smaller or equal than $0.35$.
This furthermore highlights the amenability of the 2L-VKOGA for interpretability of the ML model.
The last column ''Prediction models`` shows the number of FOM, ROM and ML-ROM evaluations (which is related to the ratios from the previous column) and the mean prediction times of these models.
The reason for this is that the size of the ML models (in terms of the number of parameters) is approximately the same, namely approximately 80000 parameters for the NN and SDKN and approximately 100000 parameters for the VKOGA models.
In terms of total runtime the NN took 48:57h, the SDKN 48:23h, the 2L-VKOGA 46:58h and the VKOGA was stopped after 54:01h (for only 2000 parameters). \\
The source code to carry out the experiments and generate the figures can be found online.\footnote{\url{https://gitlab.mathematik.uni-stuttgart.de/wenzeltn/paper-2023-deep-kernel-for-rb-ml-rom}}%

\vspace{-.7cm}

\begin{table}
	%
	\centering
	\begin{tabular}{ c || c c | c || c c c }
		\toprule
		\multirow{3}{*}{ML-Model} & \multicolumn{2}{c|}{Trainings of ML-ROM} & \hphantom{i}\multirow{3}{*}{\makecell{ML-ROM\\ratio}}\hphantom{i} & \multicolumn{3}{c}{\hphantom{i}Num.\ of model evaluations\hphantom{i}} \\
		\cline{2-3}\cline{5-7}
		& \hphantom{i}\makecell{Param.\ at which\\training started}\hphantom{i} & \hphantom{i}\makecell{Training\\time (s)}\hphantom{i} & & Model & \makecell{Num.\\evals.} & \makecell{Average\\time (s)} \\
		\midrule \midrule
		\multirow{4}{*}{NN} & & & 0.0\% & FOM & 12 & $234.9$ \\
		& $200$ & $536.3\hphantom{0}$ & 36.8\% & RB-ROM & 876 & $6.7$  \\
		& $516$ & $3919.1$ & 49.1\% & ML-ROM & 4112 & $3.2$ \\
		& $908$ & $6515.1$ & 92.9\% & & & \\ \hline
		\multirow{4}{*}{VKOGA} & & & 0.0\% & FOM & 12 & $235.9$ \\
		& $40$ & $97.8$ & $\leq 25\%$ & RB-ROM & 1789 & $8.6$  \\
		& $\cdots$ & & $\cdots$ & ML-ROM & 199 & $4.2$ \\
		& $1999$ & $3572.3$ & $\leq 25\%$ & & & \\
		\midrule \midrule
		\multirow{3}{*}{SDKN} & & & 0.0\% & FOM & 12 & $233.0$ \\
		& $200$ & $3350.2$ & 60.4\% & RB-ROM & 613 & $6.7$  \\
		& $703$ & $8272.9$ & 94.7\% & ML-ROM & 4375 & $3.3$ \\ \hline
		\hphantom{i}\multirow{3}{*}{2L-VKOGA}\hphantom{i} & & & 0.0\% & FOM & 12 & $235.9$ \\
		& $40$ & $110.3$ & 86.7\% & RB-ROM & 685 & $8.7$  \\
		& & & & ML-ROM & 4303 & $4.1$ \\
		\bottomrule
	\end{tabular}
	\caption{Overview of the results using the adaptive model for different kind of machine learning models:
		Top two boxes: Baseline NN as well as a shallow kernel VKOGA model (\Cref{subsec:ml}).
		Bottom two boxes: The deep kernel approaches SDKN and 2L-VKOGA (\Cref{sec:deep_kernels}).}
	\label{tab:overview_results}
\end{table}

\begin{figure}[hbt!]
	\centering
	\input{Figures/SDKN_200_time_to_output_1000.tex}\hspace*{-0.3cm}
	\input{Figures/VKOGA_20_time_to_output_1000.tex}
	\caption{
		Visualization of the required time (vertical axis) to build and evaluate the adaptive model $M_\text{adapt}$ for each $\mu$ chosen by the random sampling (horizontal axis).
		Left: SDKN as ML-ROM. Right: 2L-VKOGA as ML-ROM.
	}
	\label{fig:adaptive_model_results}
\end{figure}

\vspace{-.3cm}
\section{Conclusion} \label{sec:conclusion}
\vspace{-.3cm}

In this paper we reviewed the certified RB-ML-ROM model and the two deep kernel models SDKN and 2L-VKOGA, whereby we extended the second one for the use with vectorial target data.
These deep kernel models where then used for the machine learning part within the RB-ML-ROM model, where they outperformed the standard NN baseline as well as shallow kernel models. \\

	%
		%
	%
	%
	%
	%
	%
	%
	%
	%
	%

\vspace{-.3cm}
\subsubsection{Acknowledgements}
Funded by BMBF under the contracts 05M20PMA and 05M20VSA.
Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2044 --390685587,
Mathematics Münster: Dynamics--Geometry--Structure, and EXC 2075 --390740016.
We acknowledge the support by the Stuttgart Center for Simulation Science (SimTech).

\bibliographystyle{splncs04}


\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.