# Neuron to Graph: Interpreting Language Model Neurons at Scale

Alex Foote1, Neel Nanda2, Esben Kran1, Ioannis Konstas3, Shay B. Cohen4, Fazl
Barez1,4,5††footnotemark:
1Apart Research 2Independent 3 School of Mathematical and Computer Sciences
Heriot-Watt University
4 School of Informatics, University of Edinburgh 5 University of Oxford
Significant Contribution. Correspondence to<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Advances in Large Language Models (LLMs) have led to remarkable capabilities,
yet their inner mechanisms remain largely unknown. To understand these models,
we need to unravel the functions of individual neurons and their contribution
to the network. This paper introduces a novel automated approach designed to
scale interpretability techniques across a vast array of neurons within LLMs,
to make them more interpretable and ultimately safe. Conventional methods
require examination of examples with strong neuron activation and manual
identification of patterns to decipher the concepts a neuron responds to. We
propose Neuron to Graph (N$2$G), an innovative tool that automatically
extracts a neuron’s behaviour from the dataset it was trained on and
translates it into an interpretable graph. N$2$G uses truncation and saliency
methods to emphasise only the most pertinent tokens to a neuron while
enriching dataset examples with diverse samples to better encompass the full
spectrum of neuron behaviour. These graphs can be visualised to aid
researchers’ manual interpretation, and can generate token activations on text
for automatic validation by comparison with the neuron’s ground truth
activations, which we use to show that the model is better at predicting
neuron activation than two baseline methods. We also demonstrate how the
generated graph representations can be flexibly used to facilitate further
automation of interpretability research, by searching for neurons with
particular properties, or programmatically comparing neurons to each other to
identify similar neurons. Our method easily scales to build graph
representations for all neurons in a 6-layer Transformer model using a single
Tesla T4 GPU, allowing for wide usability. We release the code and
instructions for use at https://github.com/alexjfoote/Neuron2Graph.

##  1 Introduction

Interpretability of machine learning models is an active research topic [17,
1] and can have a wide range of applications from bias detection [32] to
autonomous vehicles [2] and Large Language Models (LLMs; [11]). The growing
sub-field of mechanistic interpretability aims to understand the behaviour of
individual neurons within models as well as how they combine into larger
circuits of neurons that perform a particular function [22, 20, 16], with the
ultimate aim of decomposing a model into interpretable components and using
this to ensure model safety.

Feature visualisation [23] is a tool for interpreting neurons in image models,
whereby a synthetic input image is optimised to understand a target neuron.
This has significantly aided work interpreting vision models. For example, to
identify multimodal neurons which respond to abstract concepts [16], and to
catalogue the behaviour of all early neurons in Inceptionv1 [21]. Similar
interpretability tools for understanding neurons in LLMs are lacking.
Currently, researchers often look at dataset examples containing tokens upon
which a neuron strongly activates and investigate common elements and themes
across examples to give some insight into neuron behaviour [11, 15]. However,
this can give the illusion of interpretability when real behaviour is more
complex [5], and measuring the degree to which these insights are correct is
challenging. Additionally, inspecting individual neurons by hand is time-
consuming and unlikely to scale to entire models.

To overcome these challenges, we present Neuron to Graph (N2G), which
automatically converts a target neuron within an LLM to an interpretable graph
that captures a neuron’s behaviour. Our method takes maximally activating
dataset examples for a target neuron, prunes them to remove irrelevant
context, identifies the tokens which are important for neuron activation, and
creates additional examples by replacing the important tokens with likely
substitutes using DistilBERT [28]. These processed examples are then given as
input to the graph builder, which removes unimportant tokens and creates a
condensed graph representation. The graph can be visualised to facilitate
understanding the neuron’s behaviour, as well as used to process text and
produce predicted token activations. This allows us to measure the
correspondence between the target neuron’s activations and the graph’s
structure, which provides a direct measurement of the degree to which a graph
captures the neuron’s behaviour. Once built, the graphs can be searched to
quickly identify neurons with particular properties, which could help
facilitate interpretability research.

N2G provides a promising research direction for understanding the function of
individual neurons in language models by building interpretable
representations for every neuron in a model. Our main contributions are four-
fold:

  * •

An input pruning method which removes any context from an input example that
is unnecessary for neuron activation.

  * •

Token saliency computation to identify the importance of each context token
for neuron activation.

  * •

An augmentation technique for pruned inputs to better explore and understand
neuron behaviour by generating varied inputs through predicted token
substitutions.

  * •

A graph-building process that results in a graph representation for each
neuron. The quality of the representation can be automatically measured by
comparing the output of the representation to the real neuron, the graph can
be visualised for human interpretation, and the representations can be
searched and compared to help automate parts of interpretability research.

##  2 Related Work

Neuron analysis is a branch of natural language processing (NLP) research that
investigates the structure and function of neurons within an LLM. It plays an
important role in understanding the inner workings of a model and has the
potential to enable a mechanistic understanding of large models.

Prior work in neuron analysis has identified the presence of neurons
correlated with specific concepts [26]. For instance, [8] explored neurons
which specialised in linguistic and non-linguistic concepts in large language
models, and [30] evaluated neurons which handle concepts such as causation in
language. The existence of similar concepts embedded within models can also be
found across different architectures. [33] and [29] examined neuron
distributions across models and found that different architectures have
similar localised representations of information, even when [10] used a
combination of neuron analysis and visualisation techniques to compare
transformer and recurrent models, finding that the transformer produces fewer
neurons but exhibits stronger dynamics.

There are various methods of identifying concept neurons [15]. In [3], a
method was proposed for identifying important neurons across models by
analyzing correlations between neurons from different models. In contrast, [7]
developed a method to identify concept neurons in transformer feed-forward
networks by computing the contribution of each neuron to the knowledge
prediction. In contrast, we focus on identifying neurons using highly
activating dataset examples. [18] demonstrated how the co-variance of neuron
activations on a dataset can be used to distinguish neurons that are related
to a particular concept. [31] also used neuron activations to train a probe
which automatically evaluates language models for neurons correlated to
linguistic concepts.

One limitation of using highly activating dataset examples is that the
accurate identification of concepts correlated with a neuron is limited by the
dataset itself. A neuron may represent several concepts, and [6] emphasise the
importance of conducting interpretability research on varied datasets, in
order to avoid the “interpretability illusion”, in which neurons that show
consistent patterns of activation in one dataset activate on different
concepts in another. [25] also showed the limitations of datasets in concept
neuron identification. They demonstrated that generating synthetic language
inputs that maximise the activations of a neuron surpasses naive search on a
corpus.

Researchers also used GPT-4 to simulate neurons and predict neuron activation,
and then again use GPT-4 to explain the neuron behaviour [4]. This method
represents an important step towards scalable neuron interpretability,
providing the ability to automatically generate explanations for all neurons
in a Language Model. However, whilst a natural language explanation can be
very useful, it does not enable further automated analysis in the same way a
more structured representation can. For example, searching explanations
requires complex semantic search compared to simple syntactic matching, and
similarly comparing the explanations for multiple neurons is more challenging.

##  3 Methodology

In this section, we first discuss the model architecture we seek to create
neuron interpretations for (§3.1) and then discuss our algorithm for creating
these interpretations (§3.2).

###  3.1 Model Architecture

To test our algorithm, we analyse neurons of a SoLU model [11], an auto-
regressive model from the Transformer family. This model uses the SoLU
activation function, which pushes neurons to be monosemantic by penalising
many neurons in a layer from firing at once. SoLU models have been shown to
have a higher prevalence of neurons that represent a single feature, and are
therefore more easily interpretable. They therefore provide an ideal test-bed
for our model.

Similarly to other auto-regressive Transformer models, the SoLU model
implementation we use [19] includes multi-layered perceptron (MLP) layers
within the Transformer blocks that make up the model. The model has six
blocks, each with one MLP layer containing $3072$ neurons, and was trained on
the Pile [14]. We denote the set of neurons from all MLP layers by
$\mathcal{N}$, where each element in $\mathcal{N}$ is $n_{\ell j}$, indexing
the $j$th neuron in layer $\ell$ of the SoLU model.

During its feed-forward step on a set of tokens, $x=x_{1}\cdots x_{n}$, a
given neuron $n_{\ell j}$ generates an activation over all tokens $x_{1}\cdots
x_{n}$. We denote by $a(i,x,\ell,j)$ the function that returns the activation
of neuron $j$ in layer $\ell$ on input $x$ for the $i$th token $x_{i}$. We
denote by $a(\cdot,x,\ell,j)$ to be $a(|x|,x,\ell,j)$. For integers $c,d$, we
denote by $x_{c\colon d}$ the substring $x_{c}x_{c+1}\cdots x_{d}$.

Given a fixed neuron $n_{\ell j}\in\mathcal{N}$ that we are seeking to create
an interpretation for, we obtain the maximum activating dataset sample from
the Neuroscope resource [19]. This resource provides the top 20 most
activating examples from the training dataset for each neuron in the model.
Each dataset example is a 1024-token string, and we refer to the final
activating dataset sample of strings for any particular neuron as
$\mathcal{X}=\\{x^{(1)},\ldots,x^{(m)}\\}$.

###  3.2 N$2$G: Neuron to Graph Algorithm

Figure 1: Overall architecture of N2G. Activations of the target neuron on the
dataset examples are retrieved (neuron and activating tokens in red). Prompts
are pruned and the importance of each token for neuron activation is measured
(important tokens in blue). Pruned prompts are augmented by replacing
important tokens with high-probability substitutes using DistilBERT. The
augmented set of prompts are converted to a graph. The output graph is a real
example which activates on the token “except” when preceded by any of the
other tokens.

For each neuron in the model we automatically build an interpretable
representation that aims to capture the target neuron’s behaviour. Our
algorithm is inspired by the notion that the role of a neuron, and hence its
possible interpretation, can be understood by examining the cases in which its
activation levels are high. We, therefore, aim to expand the list of sequences
in $\mathcal{X}$ to related sequences that also generate high activations for
the target neuron $n_{\ell j}$. This expanded list of sequences is then
compiled into a lattice (a form of a directed acyclic graph) where each path
in this lattice is a _minimally_ representative token sequence which results
in high activation for the neuron. We present pseudocode for the algorithm in
the Appendix.

Pruning

The algorithm begins with a pruning step, which aims to extract the key
substring of each maximally activating dataset example. More formally, for
each $x^{(i)}\in\mathcal{X}$ it identifies the pivot token index:

$e_{i}=\arg\max_{k}a(k,x^{(i)},\ell,j).$ (1)

This pivot token has the highest activation among all tokens in $x^{(i)}$ for
the target neuron $n_{\ell j}$, and therefore, is assumed to be the
representative token for interpreting $n_{\ell j}$. Given this pivot token, we
want to find the shortest prior context required to activate the neuron.
Following the identification of $e_{i}$, we find a minimal subsequence of
$x^{(i)}$ that ends in position $e_{i}$ and starts in position $s_{i}$ such
that:

$\frac{e_{i},x^{(i)}_{s_{i}}\cdots
x^{(i)}_{e_{i}},\ell,j)}{a(e_{i},x^{(i)},\ell,j)}\geq 0.5.$ (2)

Concretely, we iteratively add tokens to the left of the pivot token, each
time re-calculating the activation on the pivot token in the new subsequence
until it is at least half 111We used empirical guidance to choose this
threshold. of the activation on the pivot token in the full sequence.

The purpose of this pruning is to remove extraneous information that is
irrelevant for neuron activation, which facilitates the later stages of the
algorithm. After pruning, we end with a set of minimal subsequences
$\mathcal{Y}=\\{y^{(1)},\ldots,y^{(m)}\\}$ where each one is generated from a
single element in $\mathcal{X}$, $y^{(i)}=x^{(i)}_{s_{i}\colon e_{i}}$.

Saliency Identification

Given the set $\mathcal{Y}$ of strings, we follow another saliency
identification step. For each example $y^{(i)}\in\mathcal{Y}$, we measure the
relative importance of the $k$th token $y^{(i)}_{k}$ to the activation of the
target neuron $n_{\ell j}$ on the $k^{\prime}$th token $x_{k^{\prime}}$. We
calculate the value:

$\alpha_{k,k^{\prime}}=1-\displaystyle\frac{a(k^{\prime},y^{(i)},\ell,j)}{a(k^{\prime},\hat{y}^{(i)},\ell,j)},$
(3)

where in this context, ${\hat{y}^{(i)}}$ is the same as the sequence
$y^{(i)}$, except that the $k$th token is replaced with a special padding
token. The relative importance $\alpha_{k,l}$ indicates what would happen to
the target neuron activation on the $k^{\prime}$th token of $y^{(i)}$ if we
perturbed the $k$th token of that sequence. Intuitively, if the activation of
the neuron on token $k^{\prime}$ is much lower in the perturbed sequence than
the original sequence, the token $k$ is providing necessary context that is
important for neuron activation. This information is useful in later stages of
the algorithm, allowing us to identify the key context tokens for a given
neuron, and discard the unimportant tokens that are irrelevant to the neuron’s
behaviour.

Augmentation

We then automatically augment our set of examples $\mathcal{Y}$ with
additional examples. For each example $y^{(i)}\in\mathcal{Y}$, we identify the
tokens that are important for neuron activation on the pivot token by
thresholding the value $\alpha_{k,e_{i}}$ for all tokens at index $k<e_{i}$.
Each of these important tokens is then in turn replaced with related tokens
that are predicted from a helper model, such as BERT [9]. These tokens,
$\mathcal{R}_{k,i}$ are a set of elements in the vocabulary that the helper
model gives high probability to at position $k$ when fed with $y^{(i)}$ masked
at position $k$. Once we have generated a candidate set of augmentations
$\mathcal{A}$, we pass each new example $a^{(i)}$ through the target model and
measure the activation of the target neuron on the pivot token,
$a(e_{i},a^{(i)},\ell,j)$, and enrich $\mathcal{Y}$ with all $a^{(i)}$ where
that pivot activation is more than half of the pivot activation on the
originating sequence $y^{(i)}$:

$\frac{a(e_{i},a^{(i)},\ell,j)}{a(e_{i},y^{(i)},\ell,j)}\geq 0.5.$ (4)

Intuitively, we are exploring the space around each input example to find
other, similar, examples, and retaining the ones that still strongly activate
the target neuron. This allows us to better understand the full extent of a
neuron’s behaviour. We now update $\mathcal{Y}$ to include all augmented
examples.

For our helper model, we chose DistilBERT [28], as its bidirectionality
provides the ability to use both the prior and post context for better token
predictions. In addition, DistilBERT is a compact and efficient variant of
BERT model [9], trained using knowledge distillation, allowing it to run 60%
faster than BERT with similar statistical performance. This makes it an ideal
choice for efficiently predicting substitute tokens.

Graph Building

After completing the previous stages, we have a set of dataset examples
$\mathcal{Y}$ which strongly activate our target neuron. We then aim to build
a compact representation which concisely represents the tokens on which our
neuron activates, as well as the context necessary for activation on these
tokens. To do this, we create a lattice structure which has tokens as nodes,
where each path is constructed from a sequence in $\mathcal{Y}$ and contains
the important context tokens for activation on the final pivot token, which
ends the path.

For each example $y\in\mathcal{Y}$, we iterate over the tokens and retrieve
the importance of that token relative to the pivot token. For each token
$y_{k}$, where $k<e_{i}$, if its importance is above the threshold we add it
to the path as an important token. Otherwise, we add it to the path as a
special ignore token. We then end the path with pivot token. Ignore tokens are
placeholder’s indicating that some token is necessary here, but that the
specific token does not matter - in contrast to an important token, where that
specific token has been identified as key for neuron activation the pivot
token.

After creating the paths, we store them in a trie. We work backwards through
each path, first adding the pivot token as a top layer activating token, then
adding each context token as an important node or an ignore node as required.
At the end of each path we add a special end node, which denotes a complete
path through the trie. On this end node, we store the normalised activation of
the target neuron on the activating token for that path. We compute the
normalised activation by retrieving the maximum observed activation of the
neuron on any token in the training dataset $\mathcal{D}$,
$a_{max}=\max_{i,k}{a(k,x^{(i)},\ell,j)}$. We then compute the normalised
activation as $a(e_{i},y,\ell,j)/a_{max}$ for the example $y$ which formed the
path.

By storing our representation in this way, we can use our representation to
simulate neuron behaviour by predicting token activations. Specifically, we
can pass an input $x$ consisting of tokens $x_{1}\cdots x_{n}$ to the trie,
and it will output a predicted normalised activation for each token $x_{i}$ in
$x$. Starting at the root of the trie, we check if the token $x_{i}$ is in the
root’s child nodes. If it is, we have identified that $x_{i}$ is an activating
token, so we then continue to check if there is the necessary context for
activation. Iterating backwards starting at $x_{i-1}$, we check if the token
matches any of the current node’s children - specifically, whether the token
is either in the current node’s children, or whether the current node has an
ignore token in its children. If at any point we reach a node that has an end
node as a child, we have traced a valid path through the trie, so record the
normalised activation stored on the node. We continue these steps until we run
out of prior tokens, or the current token fails either of the matching
conditions. We then return the maximum observed normalised activation, or 0 if
we did not find any matching paths.

This representation therefore allows us to directly measure the correspondence
between a target neuron’s activations and the predicted activations from the
trie for any input text, which allows us to quantify the accuracy of the
representation. In addition, we are able to visualise the representation to
facilitate human interpretation. To do this, we remove ignore nodes as they do
not contain relevant information, then at each layer in the trie we collapse
identical tokens to form the directed graph structure. We colour context nodes
as blue according to their importance, and activating nodes as red according
to their normalised activation. Brighter colours indicate higher importance or
activation, respectively. Additionally, we denote valid end nodes with a bold
outline.

We refer to this visualisation, as well as the underlying trie, as a neuron
graph, which we can visually inspect or automatically evaluate. Figure 1 shows
an example of a simple graph for a neuron in layer 1 of the model. The neuron
graphs also facilitate new forms of analysis. For example, they provide a
simple searchable structure. We provide the ability to efficiently search
graphs by both activating tokens and context tokens. For example, we could
search for a graph that activates on the token “everyone” when it occurs with
the context token “except”, to retrieve the graph in Figure 1. This search
functionality allows us to quickly look for neurons with a particular
behaviour in which we are interested. In comparison, attempting to search
maximally activating dataset examples in the same way is much less useful, as
they often contain significant amounts of unimportant information which
renders any search on context tokens useless.

##  4 Results and Discussion

We describe a quantitative analysis, comparing our methods to simple look-up
methods (§4.1) and two case studies where we use N$2$G to discover interesting
neuron behaviour (§4.2). We also provide additional case studies in the
Appendix.

###  4.1 Quantitative Results

Given that the neuron graphs built by the algorithm can be directly used to
process text and predict token activations, we can evaluate the degree to
which they accurately capture the target neuron’s behaviour by measuring the
correspondence between the activations of the neuron and the predicted
activations of the graph on some evaluation text. We conduct all experiments
on the SoLU model described in §3.1.

For each neuron, we take the top $20$ dataset examples from Neuroscope and
randomly split them in half to form a train and test set. We give the training
examples to N$2$G to create a neuron graph, then take the test examples
compute the normalised token activations as described in §3.2. We apply a
threshold to the normalised token activations, defining an activation above
the threshold as a firing of the neuron, and an activation below the threshold
as the neuron not firing. In these experiments we set the threshold to $0.5$.
We then process the test prompts with the neuron graph to produce predicted
token firings. We can then measure the precision, recall, and $F1$ score of
the graph’s predictions compared to the ground truth firings. Building a
neuron graph for every neuron in the model took approximately 48 hours of
processing on a single Tesla T4 GPU.

To ground the results of N$2$G we compare to two simple baselines. The first
is a per-neuron token lookup table. For each neuron, we take the same training
set of 10 dataset examples which we give to the N$2$G algorithm, pass them
through the model and record the activation of the neuron on each token. We
then create a lookup table that maps every token to the maximum observed
activation for any occurrence of that token. We can then predict token
activations on text by outputting the stored activation for each token in the
input (or 0 if we did not see the token when creating the lookup).
Intuitively, this presents a very recall-focused baseline, as it ignores any
context and instead just identifies whether a token has ever been seen to
activate a neuron. The second baseline generalises the token lookup to an
n-gram lookup, following the same process as above but storing the prior $n$
context tokens for each token in the input. We choose $n=5$ for our baseline,
which provides a precision focused baseline that requires specific context for
neuron activation on any given token.

Table 1 shows the average precision, recall, and $F1$ score of the three
methods for each layer of the model, averaged across the neurons in that
layer222$F1$ scores computed on a per-neuron basis then averaged, rather than
computing on the average precision and recall.. For each neuron, we only
compute these statistics on the tokens that caused the neuron to fire. This is
because any given neuron typically fires on very few tokens, so the prediction
problem is very imbalanced - predicting the neuron will never fire would give
very good performance but provide no useful information.

In layers 0 and 1 of the model, the neuron graphs built with N$2$G on average
capture the behaviour of the neurons well, with high recall and decent
precision on the firing tokens. In comparison to the baselines, N$2$G achieves
close to the recall of the high recall baseline, with substantially better
precision. However, the precision baseline achieves significantly higher
precision than N$2$G, suggesting that the method is not fully capturing the
context needed to precisely determine when the neuron will fire.

However, as we progress to deeper layers of the model, the recall and
precision of N$2$G generally decreases. This corresponds to neurons in the
later layers on average exhibiting more complex behaviour that is less
completely captured in the training examples. Specifically, neurons in early
layers tend to respond to a small number of specific tokens in specific,
narrow contexts, whereas later layers often respond to more abstract concepts
represented by a wider array of tokens in many different contexts, which was
similarly observed by [11]. We see similar trends for the recall baseline, and
the decrease in recall in particular suggests that the training examples on
which we fit our models do not fully capture the behaviour of neurons in later
layers. This suggests there could be improvements from collecting a larger set
of dataset examples for each neuron and demonstrates the value of augmentation
to better explore the space of activating inputs.

Interestingly, the recall of the n-gram lookup increases as we go deeper into
the model. This is likely because later layers tend to require longer context
and more specific context, whereas early layers may only require very short
context, less than the 5 tokens used in the baseline. This demonstrates the
importance of the pruning and saliency mechanisms, to adaptively identify the
necessary context for activation on a per neuron basis.

| N$2$G |  |  | Token Lookup |  |  | n-Gram Lookup |  |
---|---|---|---|---|---|---|---|---|---
Layer | Precision | Recall | F1 | Precision | Recall | F1 | Precision | Recall | F1
0 | 0.67 | 0.84 | 0.68 | 0.45 | 0.89 | 0.53 | 0.97 | 0.13 | 0.17
1 | 0.64 | 0.82 | 0.65 | 0.36 | 0.90 | 0.44 | 0.96 | 0.14 | 0.18
2 | 0.56 | 0.75 | 0.55 | 0.29 | 0.85 | 0.37 | 0.93 | 0.14 | 0.17
3 | 0.45 | 0.71 | 0.45 | 0.23 | 0.82 | 0.30 | 0.91 | 0.16 | 0.18
4 | 0.39 | 0.66 | 0.38 | 0.22 | 0.78 | 0.28 | 0.89 | 0.16 | 0.18
5 | 0.38 | 0.68 | 0.37 | 0.20 | 0.81 | 0.27 | 0.86 | 0.26 | 0.28

Table 1:  N2G compared to two baseline. Note the statistics are averaged
across all neurons in a layer, and only computed for the tokens on which the
neuron fired. Best $F1$ scores in bold.

(a) | (b) | (c)
---|---|---
|  |

Figure 2:  (a) Activation trajectory averaged across 100 neurons from each
layer of the model. (b) and (c): Activation and Importance trajectories for a
neuron from layer 4. Points show occurrence of the important “idae” token.
Below, graph for the corresponding neuron from layer 4.

The trend of later layers requiring more context for activation is
demonstrated by Figure 2 a). This shows how the normalised activation of a
neuron on the pivot token in an input changes as we progressively remove
tokens from the prior context. We refer to this as an activation trajectory,
and compute the average trajectory across 100 neurons in each layer of the
model. We observe that as we move from early to later layers the average
activation trajectory becomes shallower, showing that neurons in later layers
on average require a much longer context to activate than neurons in early
layers.

We show a typical example of a neuron from layer 4 of the model in Figure 2
(b). The average trajectories show a smooth decrease in activation as we
remove context tokens, but the more typical shape of any individual neuron is
for there to be distinct jumps in activation as we reach important tokens. The
trajectories start off identically for the first five tokens corresponding to
“category: reptiles”, which are important for activation and therefore
included in the neuron graph at the bottom of the figure. There is then a
longer-range dependency on the token “idae” (which for example, could occur at
the end of the word “Gekkonidae”, the scientific name for Geckos), which can
be of variable length. There are distinct spikes where this token occurs, both
in the activation and importance trajectories, marked with points in the
graphs. This demonstrates how the neuron graph captures the key information
for activation, whilst abstracting away information about the length
dependency between tokens. Additionally, we can see that the importance
measure from the saliency mechanism closely relates to the activation
trajectories, with a spike in importance corresponding to the increase in
activation.

###  4.2 Large Scale Neuron Case Studies

Constructing the neuron graphs for every neuron in a model enables new
workflows for mechanistic interpretability research. They provide a simple
searchable structure which we can use to identify neurons with specific
properties, and can be compared to each other to find similarities between
neurons. Here, we present some case studies demonstrating how these
representations can enable these forms of investigation.

In-context Learning

Figure 3: Neuron graph for an in-context learning neuron that activates on
repeated token sequences.

In-context learning refers to the ability of Language Models to use prior
information in their prompt to better predict later tokens, and is a defining
property of these models. A key mechanism behind in-context learning is the
induction head [24], an attention head which learns to look for repeated token
sequences and increase the probability of predicting later tokens in the
sequence if the beginning of the sequence is seen again. For example, if we
see “abc…ab”, an induction head would contribute to increasing the probability
of predicting the token “c” will reoccur.

We automatically searched the graph representations to identify graphs that
contained any particular token in both the context tokens and the activating
tokens. This enabled us to discover many neurons that activate on repeated
tokens, and in particular we found what appears to be an in-context learning
neuron. Shown in Figure 3, this neuron responds to a wide variety of
apparently unrelated tokens, as long as they are part of a repeating sequence
such as “ab…ab”. This demonstrates how the neuron graph representation can
enable researchers to much more efficiently explore and understand particular
neuron behaviours via search, facilitating mechanistic interpretability
research.

Discovering Neuron Groups

Figure 4: A neuron graph that occurs for a neuron in Layer 1 and a neuron in
Layer 4 of the model. The neurons have identical behaviour, and were
recognised as a similar pair through an automated graph comparison process.

In addition, this representation offers the possibility of automatically
identifying higher level structures within Language Models, such as simple
circuits [13]. We provide an initial step towards this by automatically
identifying neurons with similar behaviours by comparing their
representations. Specifically, for every pair of neurons we measure the
proportion overlap between the pairs’ context tokens, as well as the
proportion overlap between their activating tokens. We can then automatically
identify similar neurons by retrieving pairs with more than 90% overlap for
both context tokens and activating tokens.

We identify 60 pairs of similar neurons within the model by automatically
comparing pairs of neuron graphs. Figure 4 shows an example of a neuron graph
that represents a neuron in layer 1 of the model. Our analysis identified that
an identical graph also represents a neuron in layer 4 of the model.

One possible explanation for the existence of such pairs of neurons is the
spreading out of features due to superposition. Previous research has provided
evidence for superposition, where models use individual neurons to represent
multiple unrelated features. This phenomenon presents a barrier to
interpretability as it makes understanding neuron function more challenging
[11]. Superposition implies that features are spread across multiple neurons
[22], to enable a model to represent more features than it has neurons. As
such, it could provide an explanation for the occurrence of similar neurons.
Alternatively, the model may be incentivized by dropout during the training
process to represent highly important features in multiple neurons, so that
dropping out any particular neuron does not hurt the loss. Further research
into superposition and the existence of neurons with very similar behaviours
would be beneficial to better understand the causes of these phenomena.

Superposition also has implications for attempting to control undesirable
behaviours. One method for preventing a model from exhibiting some negative
behaviour, such as toxic outputs, would be to identify specific neurons that
are responsible for this behaviour and ablate them. Indeed, [27] found such a
sentiment neuron in an LSTM model, and used it to control the sentiment of
generated text. Superposition implies that features such as sentiment or
toxicity may be spread across multiple neurons [11], so ablating any a single
neuron may not significantly affect the feature representation. The ability to
identify similar neurons therefore could facilitate identifying clusters of
neurons that contribute to a negative behaviour, and allow us to efficiently
ablate all necessary neurons. This shows the potential for the N$2$G
representations to facilitate diverse tasks in ML safety.

##  5 Conclusions and Limitations

We presented N$2$G, an approach for converting neurons in LLMs into
interpretable graphs that can be visualised, evaluated and searched. The
degree to which a neuron graph captures the behaviour of a target neuron can
be directly measured by comparing the output of the graph to the activations
of the neuron, making this method a step towards scalable interpretability
methods for LLMs. We find that neuron graphs capture neuron behaviour well for
early layers of the model but only partially capture the behaviour for later
layers due to increasingly complex neuron behaviour, and this problem will
likely become more prominent in larger models. Additionally, we evaluate our
method on a SoLU [11] model to minimise polysemanticity, though it can also be
applied to models with more common activation functions. The greater
prevalence of polysemantic neurons in typical Transformer models could reduce
the ability of N$2$G to fully capture neuron behaviour. Future work could
address these limitations by using more training examples to better cover the
full extent of a neuron’s behaviour, better exploring the input space via more
sophisticated and extensive augmentation, and generalising from exact token
matches to matching abstract concepts, for example, by using token embeddings.

The neuron graphs built by N$2$G allow for new methods of analysis. They can
be searched to identify neurons with interesting properties, which we
demonstrate by finding an in-context learning neuron that responds to repeated
sequences of the form “ab…ab”. They can also be programmatically compared with
each other, which we use to identify pairs of neurons with very similar
behaviours. These case studies show how the graph representation can enable
new forms of interpretability analysis and make some tasks much faster by
automating parts of the discovery process.

There could be significant scope to develop new processes that act on the
neuron representations to assist with interpretability research. For example,
future work could look at using the graph representations to automatically
identify circuits of neurons within LLMs that together perform a particular
task. For example, we could automatically look for instances where multiple
neurons that each respond to a single concept are combined to form a neuron in
a later layer that responds to all of those concepts, by extending the
similarity comparison to identify neurons which have a subset of another
neuron’s behaviour. In general, once we have generated the neuron graphs for
all neurons in a model, these become a valuable resource upon which
researchers can develop new tools for analysis.

##  6 Acknowledgments

We’d like to thank Michelle Lo for providing assistance to the literature
review and Charlotte Siegmann for comments on our previous draft, Mor Geva for
encouraging the idea and Logan Smith for helpful discussions. We are also
grateful to all members of the Cohort group at the School of Informatics for
helpful comments and suggestions.

## References

  * [1] Dario Amodei, Chris Olah, Jacob Steinhardt, Paul Christiano, John Schulman, and Dan Mané.  Concrete problems in ai safety.  arXiv preprint arXiv:1606.06565, 2016.
  * [2] Fazl Barez, Hosein Hasanbeig, and Alessandro Abate.  System III: Learning with domain knowledge for safety constraints.  In NeurIPS ML Safety Workshop, 2022.
  * [3] Anthony Bau, Yonatan Belinkov, Hassan Sajjad, Nadir Durrani, Fahim Dalvi, and James R. Glass.  Identifying and controlling important neurons in neural machine translation.  CoRR, abs/1811.01157, 2018.
  * [4] Steven Bills, Nick Cammarata, Dan Mossing, Henk Tillman, Leo Gao, Gabriel Goh, Ilya Sutskever, Jan Leike, Jeff Wu, and William Saunders.  Language models can explain neurons in language models.  2023\.
  * [5] Tolga Bolukbasi, Adam Pearce, Ann Yuan, Andy Coenen, Emily Reif, Fernanda Viégas, and Martin Wattenberg.  An interpretability illusion for bert.  arXiv preprint arXiv:2104.07143, 2021.
  * [6] Tolga Bolukbasi, Adam Pearce, Ann Yuan, Andy Coenen, Emily Reif, Fernanda B. Viégas, and Martin Wattenberg.  An interpretability illusion for BERT.  CoRR, abs/2104.07143, 2021.
  * [7] Damai Dai, Li Dong, Yaru Hao, Zhifang Sui, and Furu Wei.  Knowledge neurons in pretrained transformers.  CoRR, abs/2104.08696, 2021.
  * [8] Fahim Dalvi, Nadir Durrani, Hassan Sajjad, Yonatan Belinkov, Anthony Bau, and James Glass.  What is one grain of sand in the desert? analyzing individual neurons in deep NLP models.  Proceedings of the AAAI Conference on Artificial Intelligence, 33(01):6309–6317, July 2019.
  * [9] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding, 2018.  cite arxiv:1810.04805Comment: 13 pages.
  * [10] Nadir Durrani, Hassan Sajjad, Fahim Dalvi, and Yonatan Belinkov.  Analyzing individual neurons in pre-trained language models.  In Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP), pages 4865–4880, Online, November 2020\. Association for Computational Linguistics.
  * [11] Nelson Elhage, Tristan Hume, Catherine Olsson, Neel Nanda, Tom Henighan, Scott Johnston, Sheer ElShowk, Nicholas Joseph, Nova DasSarma, Ben Mann, Danny Hernandez, Amanda Askell, Kamal Ndousse, And Jones, Dawn Drain, Anna Chen, Yuntao Bai, Deep Ganguli, Liane Lovitt, Zac Hatfield-Dodds, Jackson Kernion, Tom Conerly, Shauna Kravec, Stanislav Fort, Saurav Kadavath, Josh Jacobson, Eli Tran-Johnson, Jared Kaplan, Jack Clark, Tom Brown, Sam McCandlish, Dario Amodei, and Christopher Olah.  Softmax linear units.  Transformer Circuits Thread, 2022.  https://transformer-circuits.pub/2022/solu/index.html.
  * [12] Nelson Elhage, Tristan Hume, Catherine Olsson, Nicholas Schiefer, Tom Henighan, Shauna Kravec, Zac Hatfield-Dodds, Robert Lasenby, Dawn Drain, Carol Chen, Roger Grosse, Sam McCandlish, Jared Kaplan, Dario Amodei, Martin Wattenberg, and Christopher Olah.  Toy models of superposition, 2022.
  * [13] Nelson Elhage, Neel Nanda, Catherine Olsson, Tom Henighan, Nicholas Joseph, Ben Mann, Amanda Askell, Yuntao Bai, Anna Chen, Tom Conerly, Nova DasSarma, Dawn Drain, Deep Ganguli, Zac Hatfield-Dodds, Danny Hernandez, Andy Jones, Jackson Kernion, Liane Lovitt, Kamal Ndousse, Dario Amodei, Tom Brown, Jack Clark, Jared Kaplan, Sam McCandlish, and Chris Olah.  A mathematical framework for transformer circuits.  Transformer Circuits Thread, 2021.  https://transformer-circuits.pub/2021/framework/index.html.
  * [14] Leo Gao, Stella Biderman, Sid Black, Laurence Golding, Travis Hoppe, Charles Foster, Jason Phang, Horace He, Anish Thite, Noa Nabeshima, Shawn Presser, and Connor Leahy.  The Pile: An 800gb dataset of diverse text for language modeling.  arXiv preprint arXiv:2101.00027, 2020.
  * [15] Mor Geva, Roei Schuster, Jonathan Berant, and Omer Levy.  Transformer feed-forward layers are key-value memories.  arXiv preprint arXiv:2012.14913, 2020.
  * [16] Gabriel Goh, Nick Cammarata, Chelsea Voss, Shan Carter, Michael Petrov, Ludwig Schubert, Alec Radford, and Chris Olah.  Multimodal neurons in artificial neural networks.  Distill, 6(3):e30, 2021.
  * [17] Dan Hendrycks, Nicholas Carlini, John Schulman, and Jacob Steinhardt.  Unsolved problems in ml safety.  arXiv preprint arXiv:2109.13916, 2021.
  * [18] Jesse Mu and Jacob Andreas.  Compositional explanations of neurons.  CoRR, abs/2006.14032, 2020.
  * [19] Neel Nanda.  Neuroscope, 2022.
  * [20] Chris Olah.  Mechanistic interpretability, variables, and the importance of interpretable bases.  Transformer Circuits Thread(June 27). http://www. transformer-circuits. pub/2022/mech-interp-essay/index. html, 2022.
  * [21] Chris Olah, Nick Cammarata, Ludwig Schubert, Gabriel Goh, Michael Petrov, and Shan Carter.  An overview of early vision in inceptionv1.  Distill, 5(4):e00024–002, 2020.
  * [22] Chris Olah, Nick Cammarata, Ludwig Schubert, Gabriel Goh, Michael Petrov, and Shan Carter.  Zoom in: An introduction to circuits.  Distill, 2020.  https://distill.pub/2020/circuits/zoom-in.
  * [23] Chris Olah, Alexander Mordvintsev, and Ludwig Schubert.  Feature visualization.  Distill, 2017.  https://distill.pub/2017/feature-visualization.
  * [24] Catherine Olsson, Nelson Elhage, Neel Nanda, Nicholas Joseph, Nova DasSarma, Tom Henighan, Ben Mann, Amanda Askell, Yuntao Bai, Anna Chen, Tom Conerly, Dawn Drain, Deep Ganguli, Zac Hatfield-Dodds, Danny Hernandez, Scott Johnston, Andy Jones, Jackson Kernion, Liane Lovitt, Kamal Ndousse, Dario Amodei, Tom Brown, Jack Clark, Jared Kaplan, Sam McCandlish, and Chris Olah.  In-context learning and induction heads.  Transformer Circuits Thread, 2022.  https://transformer-circuits.pub/2022/in-context-learning-and-induction-heads/index.html.
  * [25] Nina Poerner, Benjamin Roth, and Hinrich Schütze.  Interpretable textual neuron representations for NLP.  In Proceedings of the 2018 EMNLP Workshop BlackboxNLP: Analyzing and Interpreting Neural Networks for NLP, pages 325–327, Brussels, Belgium, November 2018. Association for Computational Linguistics.
  * [26] Alec Radford, Rafal Jozefowicz, and Ilya Sutskever.  Learning to generate reviews and discovering sentiment.  arXiv preprint arXiv:1704.01444, 2017.
  * [27] Alec Radford, Rafal Józefowicz, and Ilya Sutskever.  Learning to generate reviews and discovering sentiment.  ArXiv, abs/1704.01444, 2017.
  * [28] Victor Sanh, Lysandre Debut, Julien Chaumond, and Thomas Wolf.  Distilbert, a distilled version of bert: smaller, faster, cheaper and lighter, 2020.
  * [29] Ludwig Schubert, Chelsea Voss, Nick Cammarata, Gabriel Goh, and Chris Olah.  High-low frequency detectors.  Distill, 6(1):e00024–005, 2021.
  * [30] Esther Seyffarth, Younes Samih, Laura Kallmeyer, and Hassan Sajjad.  Implicit representations of event properties within contextual language models: Searching for “causativity neurons”.  In Proceedings of the 14th International Conference on Computational Semantics (IWCS), pages 110–120, Groningen, The Netherlands (online), June 2021. Association for Computational Linguistics.
  * [31] Lucas Torroba Hennigen, Adina Williams, and Ryan Cotterell.  Intrinsic probing through dimension selection.  In Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP), pages 197–216, Online, November 2020. Association for Computational Linguistics.
  * [32] Jesse Vig, Sebastian Gehrmann, Yonatan Belinkov, Sharon Qian, Daniel Nevo, Yaron Singer, and Stuart Shieber.  Investigating Gender Bias in Language Models Using Causal Mediation Analysis.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 12388–12401. Curran Associates, Inc., 2020.
  * [33] John Wu, Yonatan Belinkov, Hassan Sajjad, Nadir Durrani, Fahim Dalvi, and James Glass.  Similarity analysis of contextual word representation models.  In Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, pages 4638–4655, Online, July 2020. Association for Computational Linguistics.

##  Appendix A Neuron Graph Examples

In this appendix we explore some interesting or characteristic behaviours of
the neuron graphs. Polysemanticity, the phenomenon where a neuron exhibits
multiple unrelated behaviours, is one of the current major challenges of
neuron interpretability [12]. When present, polysemanticity often shows up
clearly in the neuron graphs as distinct, disconnected subgraphs. For example,
in Figure 5, there are three separate subgraphs corresponding to three clearly
distinct behaviours. The top subgraph responds to a phrase in Dutch -
variations on _de betrokken_ , where not all tokens in _betrokken_ were
important enough to include in the graph. The middle subgraph responds to a
phrase in English - variations on _a fun, over the top_. The bottom subgraph
responds to a phrase in Swedish - _kollegers berättigade_ , with unimportant
tokens not included. This natural separation of behaviours into separate
subgraphs could potentially make it easier to interpret polysemantic neurons,
but more experimentation would be needed to develop and test this further.

Figure 5: A neuron graph exhibiting polysemanticity, with three disconnected
subgraphs each responding to a phrase in a different language. Figure 6:
Neurons related to programming syntax, specifically import statements. Top -
Neuron graph illustrating import syntax for the Go programming language.
Bottom Left: Neuron graph showing fundamental elements of Python import
syntax. Bottom Right: Neuron graph for a neuron that responds to the imports
of widely-used Python packages.

Neuron graphs appear to work particularly well for “syntactic neurons” that
respond to structural patterns, concisely and simply capturing the structure.
One rich source of syntactic text is programming, suggesting that N$2$G could
particularly useful for analysing models trained to write code. We use the
search capability on the neuron graphs of the SoLU model to identify many
neurons related to import statements in various languages. Figure 6 shows
three examples from different layers of the model. The top and bottom left are
from layer 2, and represent basic syntax in Go and Python respectively. The
bottom right graph is for a neuron in layer 4, and shows a neuron that
responds to imports of common Python packages. The circuits line of research
[22] would suggest that later-layer neurons like the one in layer 4 may be
"composed" of neurons in early layers - for example, a simple way to do this
would be to union over several neurons that each respond to an import of one
of the packages (i.e., the later composed neuron activates if any of the
previous neurons activates). Moving from understanding individual neurons to
understanding circuits of neurons is a crucial step in interpretability
research. The ability to automatically identify similar neurons could be
expanded to identify neurons that have a subset of another neuron’s behaviour,
which could provide a method for discovering simple circuits in language
models. This demonstrates how the flexible representations built by N$2$G
could help facilitate new methods for interpretability research.

##  Appendix B N$2$G Pseudocode

Algorithm 1

1:procedure N$2$G Algorithm Input: Target neuron $n_{\ell j}$, list of
sequences $\mathcal{X}$ Prune, Saliency Identification and Augment Steps:

2: for $x^{(i)}\in\mathcal{X}$ do

3: Compute $e_{i}\leftarrow$ pivot token index

4: Find $y^{(i)}\leftarrow$ minimal sub sequence with activation ratio $\geq
0.5$

5: Form $\mathcal{Y}\leftarrow\\{y^{(1)},\ldots,y^{(m)}\\}$

6: for $y^{(i)}\in\mathcal{Y}$ do

7: Compute relative importance value $\alpha_{i,k}$ for each token in
$y^{(i)}_{k}$

8: Identify tokens with high $\alpha_{i,k}$

9: Obtain replacement tokens $\mathcal{R}_{k,i}$ from helper model Construct
Lattice and Optimise Steps:

10: Combine $y^{(i)}$ and $\mathcal{R}_{k,i}$ to form a lattice of augmented
minimal subsequences

11: Find optimal token combination in the lattice that maximises target neuron
activation $n_{\ell j}$ Output:

12: The optimal lattice of augmented minimal subsequences representing high
activation contexts for the target neuron $n_{\ell j}$

##  Appendix C Broader Impact

It is worth mentioning that the field of mechanistic interpretability is ideal
for understanding the black box nature of neural networks. This approach can
help improve model comprehension and enable researchers to build more
transparent AI systems. However, it is essential to recognize that as models
become more capable, they may also be used for purposes that are not aligned
with societal needs and safety. This presents potential ethical concerns and
underscores the need for responsible development and implementation of such
technologies. Our work helps illustrate in an accessible manner the inner-
workings of neural networks, a step towards aligning their use with
responsible use. Researchers and practitioners must remain vigilant in
addressing these challenges and ensuring that AI advances contribute
positively to society. Additionally, interdisciplinary collaborations and
public discussions can aid in raising awareness and developing robust
strategies to mitigate potential risks.


# Non-linear Topology Optimization of District Heating Networks: A benchmark
of Mixed-Integer and Adjoint Approaches

Yannick Wack 111Department of Mechanical Engineering, KU Leuven,
Celestijnenlaan 300 box 2421, 3001 Leuven, Belgium222EnergyVille, Thor Park,
Poort Genk 8310, 3600 Genk, Belgium333Flemish Institute for Technological
Research (VITO), Boeretang 200, 2400 Mol, Belgium<EMAIL_ADDRESS>Sylvain Serra 444Universite de Pau et des Pays de l’Adour, E2S UPPA, LaTEP,
Pau, France  Martine Baelmans 555Department of Mechanical Engineering, KU
Leuven, Celestijnenlaan 300 box 2421, 3001 Leuven, Belgium666EnergyVille, Thor
Park, Poort Genk 8310, 3600 Genk, Belgium Jean-Michel Reneaume 777Universite
de Pau et des Pays de l’Adour, E2S UPPA, LaTEP, Pau, France  Maarten Blommaert
888Department of Mechanical Engineering, KU Leuven, Celestijnenlaan 300 box
2421, 3001 Leuven, Belgium999EnergyVille, Thor Park, Poort Genk 8310, 3600
Genk, Belgium

###### Abstract

The widespread use of optimization methods in the design phase of District
Heating Networks is currently limited by the availability of scalable
optimization approaches that accurately represent the network. In this paper,
we compare and benchmark two different approaches to non-linear topology
optimization of District Heating Networks in terms of computational cost and
optimality gap. The first approach solves a mixed-integer non-linear
optimization problem that resolves the binary constraints of pipe routing
choices using a combinatorial optimization approach. The second approach
solves a relaxed optimization problem using an adjoint optimization approach,
and enforces a discrete network topology through penalization. Our benchmark
shows that the relaxed penalized problem has a polynomial computational cost
scaling, while the combinatorial solution scales exponentially, making it
intractable for practical-sized networks. We also evaluate the optimality gap
between the two approaches on two different District Heating Network
optimization cases. We find that the mixed-integer approach outperforms the
adjoint approach on a single-producer case, but the relaxed penalized problem
is superior on a multi-producer case. Based on this study, we discuss the
importance of initialization strategies for solving the optimal topology and
design problem of District Heating Networks as a non-linear optimization
problem.

###### keywords:

District Heating Networks, topology optimization , mixed-integer non-linear
optimization , benchmark

††journal: Energy

##  1 Introduction

District Heating Networks (DHN) are considered one of the core technologies to
enable carbon-neutral space heating [1]. It has the ability to connect a
multitude of different renewable heat sources and provide heat to districts
and entire cities. In 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 Mixed
Integer Non-Linear Program (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öderman
[2], for example, optimized the structure and configuration of a district
cooling network using a MILP approach. Dorfner and Hamacher [3] later used a
linear approach to optimize the topology and pipe sizes of a DHN. Haikarainen
et al. [4] optimized the topology and operations, while accounting for
different production technologies and heat storage. Mazairac et al. [5]
optimized the topology of a multi-carrier network that incorporates gas and
electricity supply. Morvay et al. [6] optimally designed the network while
also optimizing the energy mix supplied. Still using a MILP approach, Bordin
et al. [7] optimized the network topology while studying the set of consumers
to optimally connect to a DHN. Most recently, Resimont et al. [8] used a MILP
approach to optimize a city-scale heating network. The transformation of
modern DHN s towards multi-source, low-temperature networks [9] 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 DHN s requires a non-linear representation of the network physics.
Solving the resulting MINLP problem is challenging. In the literature, one
approach to solving this MINLP problem is to use heuristic search algorithms.
For example, Li and Svendsen used a genetic algorithm [10], and Allen et al.
deployed both a minimal spanning tree heuristic and a particle swarm algorithm
[11] 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. [12], solves the full MINLP
resolving the discrete nature of pipe routing choices, using a combinatorial
optimization approach. On the other hand, Blommaert et al. [13] and Wack et
al. [14], solve a relaxed Non Linear Problem (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 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. [12] and the relaxed adjoint approach by Wack et al. [14] to non-
linear topology optimization of 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 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.

##  2 The topology optimization problem of DHN s

To compare both approaches, a topology optimization problem is defined that
aligns the problem definitions of Mertz et al.[12] and Wack et al [14]. 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 DHN, a
set of design variables $\bm{\varphi}=\left[\bm{d},\bm{\gamma}\right]$ is
defined, containing the pipe diameters $\bm{d}$ as well as the normalized
producer inflows $\bm{\gamma}$. To represent the physical state of a given
network, a vector of physical variables
$\bm{\bm{x}}=\left[\bm{q},\bm{p},\bm{\theta}\right]$ is defined. It contains
the flow rates $\bm{q}$, nodal pressures $\bm{p}$ and nodal and pipe exit
temperatures $\bm{\theta}$. The temperatures $\bm{\theta}=\bm{T}-T_{\infty}$
are defined as the difference between the absolute water temperature $\bm{T}$
and the outside air temperature $T_{\infty}$. Now the topology optimization
problem for a DHN can be posed as a generic optimization problem of the form:

$\displaystyle\min_{\bm{\varphi},\bm{x}}$
$\displaystyle\qquad\mathcal{J}\left(\bm{\varphi},\bm{x}\right)$ (1)
$\displaystyle s.t.$ $\displaystyle\qquad\bm{h}(\bm{\varphi},\bm{x})\leq 0.$

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}$.
$\bm{h}(\bm{\varphi},\bm{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 A. This optimization problem is based on Mertz
et al. [12] and Wack et al. [14] and a detailed description of the underlying
costs and models can also be found there. In the problem definition, the set
definition in table 1 is used to reference different parts of the network.

Table 1: Notation used to describe different subsets of the DHN. Set | Symbol
---|---
All edges | $E$
Producer input | $E_{\mathrm{pr}}$
Pipes | $E_{\mathrm{pipe}}$
Feed edges | $E_{\mathrm{f}}$
Heating system | $E_{\mathrm{hs}}$

##  3 Methodology

To compare the two approaches, the topology optimization problem (described in
equation 2) is solved using two different methods: solving the full MINLP
problem in GAMS using the implementation by Mertz et al. [12] (referred to as
fMINLP for simplicity), and using the penalized adjoint optimization approach
by Wack et al. [14] (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.

### Combinatorial MINLP approach: fMINLP

For the comparison, the topology optimization problem will be solved using a
combinatorial MINLP approach. This approach was previously published by Mertz
et al. [12], and a comprehensive discussion of the approach can be found
there. It is implemented in GAMS and uses the Mixed Integer Linear Program
(MILP) solver CPLEX and the NLP solver CONOPT. To pose the optimization
problem as a combinatorial problem, the topological choice of pipe placement
is represented by a binary variable $\bm{\phi}$. Physical variables defined on
these pipes are then coupled to this existence variable using the bigM method
($\mathcal{M}$) [12], as e.g. in the flow velocity definition:

$\displaystyle m_{ij}-\mathcal{M}\left(1-\phi_{ij}\right)\leq
v_{ij}\rho\pi\frac{1}{4}d_{ij}^{2}$ (2) $\displaystyle
v_{ij}\rho\pi\frac{1}{4}d_{ij}^{2}\leq
m_{ij}+\mathcal{M}\left(1-\phi_{ij}\right),\quad\forall ij\in
E_{\mathrm{pipe}}.$ (3)

Here $\mathcal{M}\gg m_{ij}$ is chosen in such a way that if a pipe exists
$\phi_{ij}=1$, the bounds on the flow velocity constraint are tight. If
$\phi_{ij}=0$, the bounds of this constraints are relaxed. To force the
installation of both feed and return pipes the following constraint is used:

$\phi_{ij}=\phi_{ji},\quad\forall ij\in E_{\mathrm{pipe}}.$ (4)

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:

$\sum_{ij\in E_{\mathrm{pr},\mathrm{f}}}H_{ij}=\sum_{ij\in
E_{\mathrm{pr},\mathrm{f}}}\left(q\theta\right)_{ij}\rho c_{\mathrm{p}}\leq
1.5\sum_{kl\in E_{\mathrm{hs}}}Q_{kl}.$ (5)

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

$\displaystyle v_{ij}$ $\displaystyle\leq 3.5,$ (6) $\displaystyle v_{ij}$
$\displaystyle\leq 18.438d_{ij}+0.2186,$ (7) $\displaystyle v_{ij}$
$\displaystyle\geq v_{\mathrm{min}}\phi_{ij},$ (8) $\displaystyle p_{i}-p_{j}$
$\displaystyle\leq$200\text{\,}\mathrm{Pa}$\quad\forall ij\in
E_{\mathrm{pipe}}.$ (9)

This constraints are not active in the final optimized design, as they are
dominated by the pipe momentum equation (20). A fixed pressure drop is set for
the heat exchanger

$p_{i}-p_{j}\geq$2\text{\,}\mathrm{kPa}$,\quad\forall ij\in E_{\mathrm{hs}},$
(10)

and the exit temperature of the radiator is required to be smaller than the
inflow temperature:

$\theta_{i}\geq\theta_{j},\quad\forall ij\in E_{\mathrm{hs}}.$ (11)

Finally, this optimization strategy uses a sequence of initializations to
ensure a stable convergence of the final MINLP solve, based on the strategy of
Marty et al. [15]. For completeness, this initialization sequence is briefly
described in B.

### Adjoint-based penalization approach: pNLP

The other considered approach is the adjoint-based penalization approach by
Wack et al. [14]. 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 16) with

$\bar{p}_{0}(d_{ij})=p_{0}\left(\frac{1}{(1+\exp(-k\left(d_{ij}-d_{\mathrm{min}}\right))})-1\right)\,,$
(12)

or by explicitly penalizing the $d$ in the investment cost function and within
the momentum and energy equations as described in Wack et al. [14].

Because of this relaxation, the MINLP reduces to solving a series of NLP s.
Here, the optimization is initialized from a uniform distribution of pipe
diameters.

##  4 Benchmarking computational cost - a single producer case

Modern 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 MINLP
approach (fMINLP) is compared to solving a relaxed penalized NLP (pNLP). This
is done on a benchmark case containing a single producer.

### Setup

To be able to compare the computational cost scaling of both optimization
approaches with increasing network size an easily scalable DHN optimization
case is set up. This benchmark case is visualized in figure 1.

Figure 1: 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.

Here, a heating network is to be designed around a single heat producer in the
center of the network. This producer provides heat at
$70\text{\,}\mathrm{\SIUnitSymbolCelsius}$. Around this central producer,
houses and possible pipe connections are arranged in a circular manner. All
houses have a heat demand of $Q=$15\text{\,}\mathrm{kW}$$, and their heating
system characterized by $\xi=$200\text{\,}\mathrm{W}\mathrm{/}\mathrm{K}^{n}$$
and $n=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 1 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).

### 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 2 on a a) semi-log graph and a b) log-log graph.

Figure 2: 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\geq 663$. All non-converged optimizations are
excluded from the figure and subsequent analysis

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 MINLP in
the case of the fMINLP vs solving a relaxed 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 2 a)). The wall time scaling of the fMINLP can here be described
by the function
$w_{\mathrm{fMINLP{}}}\left(n\right)=51.43e^{0.012n}$\mathrm{s}$$ with an
coefficient of determination $R^{2}=0.87$, while the wall time scaling of the
pNLP can be described by
$w_{\mathrm{pNLP{}}}\left(n\right)=47.31e^{0.006n}$\mathrm{s}$$ 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
$w_{\mathrm{pNLP{}}}\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 $w_{\mathrm{fMINLP{}}}\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 MINLP s. In the context of 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 DHN s, containing
thousands of potential pipe connections. For example to optimize a network
containing $n=2000$ potential pipe connections, an optimization time of
$w_{\mathrm{fMINLP{}}}\left(n\right)\approx 43000~{}\mathrm{years}$ would be
needed following the observed exponential trend of the MINLP solution time. In
the context of topology optimization of DHN s it is computationally favorable
to solve a relaxed 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 DHN optimization feasible. This is
highlighted for the above mentioned example of a network of $n=2000$, where
the necessary optimization time reduces to
$w_{\mathrm{pNLP{}}}\left(n\right)\approx 5.3~{}\mathrm{hours}$ when using the
relaxed NLP implementation of the pNLP, assuming the observed polynomial
trend.

The potentially exponential computational cost scaling of solving full MINLP s
should be taken into account when choosing optimization strategies for DHN
topologies. It will make the optimization of large-scale DHN s intractable and
hinder the application of such optimization tools to real-world DHN problems.
It is therefore advisable to simplify the MINLP formulation by either
linearizing the non-linear model constraints, posing the problem as a MILP, or
by relaxing the integer constraints, reformulating the problem as a 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.

### 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 3.111To prevent potential remaining
modelling differences from skewing the comparison, both found optimal designs
where evaluated in a simulation within the pNLP.

Figure 3: 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\geq 663$. All non-converged optimizations are excluded from
the figure and subsequent analysis

It can be seen that in this case, the fMINLP implementation (Figure 3 red)
consistently finds a cheaper DHN design than the pNLP implementation (Figure 3
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 3 green). As can be seen in
figure 3, 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.

Figure 4: 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.

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 4. 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 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 MINLP
as well, as can be seen in figure 4 b) and d). The pNLP on the other hand,
initializing the full MINLP from a uniform distribution of pipes, finds a
different local optimum using an overall larger pipe length (see figure 4).
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 5 it will be studied how the
initialization of both approaches perform on a multi producer case.

##  5 The influence of initialization - a two producer case

Modern $4^{\text{th}}$ generation 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.

### 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 5. The circular network superstructure for this case is equivalent to
the case in section 4, representing a potential district to be connected to a
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 $70\text{\,}\mathrm{\SIUnitSymbolCelsius}$ at
high costs
($C_{\mathrm{hC}}=$800\text{\,}\mathrm{\text{€}}\text{/}\mathrm{kW}$$,
$C_{\mathrm{hO}}=$8\text{\,}\mathrm{ct}\text{/}\text{(}\mathrm{kW}\text{\,}\mathrm{h}\text{)}$$)
while the producer on the right supplies heat at
$55\text{\,}\mathrm{\SIUnitSymbolCelsius}$ at a lower cost
($C_{\mathrm{hC}}=$0\text{\,}\mathrm{\text{€}}\text{/}\mathrm{kW}$$,
$C_{\mathrm{hO}}=$4\text{\,}\mathrm{ct}\text{/}\text{(}\mathrm{kW}\text{\,}\mathrm{h}\text{)}$$).
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$50\text{\,}\mathrm{\SIUnitSymbolCelsius}$$, the rest of
the houses require heat at higher temperatures
($\gtrapprox$60\text{\,}\mathrm{\SIUnitSymbolCelsius}$$) 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.

Figure 5: 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

### 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 6.222Again the optima of both approaches
where evaluated with a pNLP simulation to avoid discrepancies.

Figure 6: 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.

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
$2.1\text{\,}\mathrm{M\text{€}}$ lower, the optimal design found by the pNLP
achieves heat cost savings of $14.4\text{\,}\mathrm{M\text{€}}$, 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 7, 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 6. The initialization steps of the fMINLP333Further
described in B 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
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.

Figure 7: 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.

The results of the benchmarking study indicate that the initialization
strategy for the topology optimization problem of 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.
[15].

##  6 Conclusion

In this paper, we compared the performance of two different approaches to the
non-linear topology optimization of DHN s. While the approach by Mertz et al.
[12] solves a full MINLP, resolving the discrete nature of pipe routing
choices using a combinatorial optimization approach (fMINLP), the approach by
Wack et al [14] solves a relaxed 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 DHN’s, solving the full
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 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 DHN optimization
feasible. We highlighted the scaling difference of both approaches by showing
that for a DHN featuring $n=2000$, an optimization time of
$w_{\mathrm{fMINLP{}}}\left(n\right)\approx 43000~{}\mathrm{years}$ would be
needed following the observed exponential trend of the MINLP solution time in
the fMINLP implementation. The necessary optimization time reduces to
$w_{\mathrm{pNLP{}}}\left(n\right)\approx 5.3~{}\mathrm{hours}$ when using the
relaxed 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 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 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 MINLP does not necessarily lead to better optimal
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 DHN s as a 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
MILP approaches to the topology optimization of DHN’s.

## 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: 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: https://doi.org/10.5281/zenodo.7434451

## 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).

## CRediT authorship contribution statement

Yannick Wack: Conceptualization, Formal analysis, Software, Visualization,
Writing – original draft. Sylvain Serra: Conceptualization, Formal analysis,
Software, Writing - Review & Editing. Martine Baelmans: Conceptualization,
Writing - Review & Editing. Jean-Michel Reneaume: Conceptualization, Writing -
Review & Editing. Maarten Blommaert: Conceptualization, Software, Writing -
Review & Editing, Supervision.

## Compliance with ethical standards

#### Conflict of interest

The authors declare that they have no conflict of interest

#### Funding

The authors did not receive support from any organization for the submitted
work.

#### Ethical approval

This article does not contain any studies with human participants or animals
performed by any of the authors.

##  Appendix A Detailed definition of the topology optimization problem for
DHN s

In this section, the underlying cost function $\mathcal{J}$ and the set of
model and technological constraints $\bm{h}(\bm{\varphi},\bm{x})$ for this
benchmark are defined in detail.

###  A.1 Cost function

The objective function is defined as the total cost of the project over an
investment horizon of $A=30~{}\mathrm{years}$:

$\displaystyle\mathcal{J}\left(\bm{\varphi},\bm{x}\right)$
$\displaystyle=f_{\mathrm{CAP}}\left(J_{\mathrm{pipe,CAP}}\left(\bm{\varphi}\right)+J_{\mathrm{h,CAP}}\left(\bm{x}\right)\right)$
$\displaystyle+f_{\mathrm{OP}}\left(J_{\mathrm{h,OP}}\left(\bm{x}\right)+J_{\mathrm{p,OP}}\left(\bm{x}\right)\right)\,,$
(13)

with

$\displaystyle f_{\mathrm{CAP}}=$
$\displaystyle\left(1+e_{\mathrm{a}}\right)^{A},$ (14) $\displaystyle
f_{\mathrm{OP}}=$
$\displaystyle\frac{1-\left(1+e_{\mathrm{a}}\right)^{A}\left(1+e_{\mathrm{i}}\right)^{A}}{1-\left(1+e_{\mathrm{a}}\right)\left(1+e_{\mathrm{i}}\right)}\,,$
(15)

assuming a discount rate $e_{\mathrm{a}}=0.04$, and an energy inflation rate
$e_{\mathrm{i}}=0.04$ [12]. The invest cost of piping $J_{\mathrm{pipe,CAP}}$
is approximated with a linear interpolation of the catalogue cost of
commercially available pipe diameters and the trench cost:

$J_{\mathrm{pipe,CAP}}\left(\bm{\varphi}\right)=\sum_{ij\in
E_{\mathrm{pipe}}}\left(p_{1}d_{ij}+p_{0}\right)L_{ij}\,.$ (16)

With the interpolation coefficients
$p_{1}=$1976.3\text{\,}\mathrm{\text{€}}\text{\,}{\mathrm{m}}^{-2}$$,
$p_{0}=$301.4\text{\,}\mathrm{\text{€}}\text{\,}{\mathrm{m}}^{-1}$$. The
investment cost for building heat production plants is calculated using

$J_{\mathrm{h,CAP}}\left(\bm{x}\right)=\frac{\rho
c_{\mathrm{p}}}{F\eta_{\mathrm{pr}}}\sum_{ij\in
E_{\mathrm{pr},\mathrm{f}}}\left(C_{\mathrm{hC}}\,q\,\Delta\theta\right)_{ij}\,,$
(17)

with the capacity price of heat production
$C_{\mathrm{hC}}=$800\text{\,}\mathrm{\text{€}}\mathrm{/}{\mathrm{k}\mathrm{W}}$$,
the capacity factor $F=0.33$ and assuming an efficiency of
$\eta_{\mathrm{pr}}=0.9$. The operational heat cost is calculated using

$J_{\mathrm{h,OP}}\left(\bm{x}\right)=\frac{\rho
c_{\mathrm{p}}}{\eta_{\mathrm{pr}}}8760\frac{$\mathrm{h}$}{$\mathrm{yr}$}\sum_{ij\in
E_{\mathrm{pr},\mathrm{f}}}\left(C_{\mathrm{hO}}\,q\,\Delta\theta\right)_{ij}\,,$
(18)

with the unit price of heat
$C_{\mathrm{hO}}=$0.06\text{\,}\mathrm{\text{€}}\mathrm{/}{\mathrm{k}\mathrm{W}\mathrm{h}}$$.
The operational cost of pumps at the heat production sites is computed with

$J_{\mathrm{p,OP}}\left(\bm{x}\right)=\frac{8760}{\eta_{\mathrm{pump}}}\frac{$\mathrm{h}$}{$\mathrm{yr}$}\sum_{ij\in
E_{\mathrm{pr},\mathrm{f}}}C_{\mathrm{pO},ij}\left(p_{j}-p_{a}\right)q_{ij}\,,$
(19)

with the electricity price
$C_{\mathrm{pO}}=$0.11\text{\,}\mathrm{\text{€}}\mathrm{/}{\mathrm{k}\mathrm{W}\mathrm{h}}$$
and a pump efficiency of $\eta_{\mathrm{pump}}=0.7$.

###  A.2 DHN model constraints

For all pipe junctions in the network conservation of mass and energy is
assumed.

#### Pipe model

The momentum equations in the pipes are modelled using the Blasius friction
factor $f$ and assume a singular pressure drop of $30\%$ [12]:

$(p_{i}-p_{j})=\frac{100}{70}f_{ij}\frac{8\rho
L_{ij}}{d_{ij}^{5}\pi^{2}}\lvert q_{ij}\rvert q_{ij}\,,\quad\forall ij\in
E_{\mathrm{pipe}},$ (20) $\mathrm{with}\quad
f_{ij}=0.3164\left(Re\right)^{-\frac{1}{4}}\,,\quad\forall ij\in
E_{\mathrm{pipe}}.$ (21)

The energy equation over pipes are modelled accounting for the thermal
conductivity of the pipe insulation
$\lambda_{\mathrm{i}}=$0.03\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$$
and the surrounding soil
$\lambda_{\mathrm{g}}=$1.4\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$$,

$\theta_{ij}=\theta_{i}\exp{\left(\frac{-L_{ij}}{\rho c_{\mathrm{p}}\lvert
q_{ij}\rvert{U_{ij}}}\right)},\quad\forall ij\in E_{\mathrm{pipe}},$ (22)
$U_{ij}=\frac{\ln(4h/(rd_{ij}))}{2\pi\lambda_{\mathrm{g}}}+\frac{\ln{r}}{2\pi\lambda_{\mathrm{i}}},$
(23)

assuming an insulation ratio $r=1.4$ and a pipe depth
$h=$0.4\text{\,}\mathrm{m}$$ [13].

#### Consumer model

In the consumer arc, conservation of momentum is assumed. Conservation of
energy in the heating system leads to

$\rho c_{\mathrm{p}}q_{ij}(\theta_{i}-\theta_{ij})=Q_{ij},\quad\forall ij\in
E_{\mathrm{hs}},$ (24)

with $Q_{ij}$ the heat transferred to the house through the heating system.
The latter is modelled with the characteristic equation for radiators [16]
using the $LMTD$ approximation by Chen [17]:

$Q_{ij}=\xi_{ij}\left(LMTD\left(\theta_{i}-\theta_{\textrm{house}},\theta_{ij}-\theta_{\textrm{house}}\right)\right)^{n_{ij}}.$
(25) $\displaystyle\mathrm{with}\quad$ $\displaystyle
LMTD\left(\Delta\theta_{\mathrm{A}},\Delta\theta_{\mathrm{B}}\right)$ (26)
$\displaystyle\approx\left(\Delta\theta_{\mathrm{A}}\Delta\theta_{\mathrm{B}}\left(\frac{\Delta\theta_{\mathrm{A}}+\Delta\theta_{\mathrm{B}}}{2}\right)\right)^{\frac{1}{3}}\,.$
(27)

Here, $\theta_{\textrm{house}}$ is the inside temperature of the house, and
$\xi$ and $n$ are heating system specific coefficients.

#### Producer model

At the producer, heat is injected with an input flow $\gamma$ at a fixed
temperature $\Theta$:

$q_{ij}=\gamma_{ij},\quad\theta_{ij}=\Theta_{ij}\quad\forall ij\in
E_{\mathrm{pr},\mathrm{f}}.$ (28)

### Additional state constraints

To ensure that the heat demand $Q_{\mathrm{d},ij}$ of every consumer is met,
the following constraint is defined:

$\displaystyle Q_{ij}-Q_{\mathrm{d},ij}\geq 0\,,\quad\forall ij\in
E_{\mathrm{hs}}.$ (29)

##  Appendix B Initialization strategy of the combinatorial approach
(fMINLP)

The initialization strategy for the MINLP approach is crucial for the overall
solution process. The initial values for the variables and constraints in the
MINLP model can greatly affect the solution’s speed, accuracy, and
feasibility. Figure 8 shows the specific initialization strategy used in the
fMINLP implementation.

Figure 8: Initialization strategy of the combinatorial approach (fMINLP). The
optimal network topology $\phi$ determined by the MIP is used to initialize
the following NLPs (red). The design variables $\bm{\varphi}$ and physical
variables $\bm{x}$ are optimized in a sequence of NLPS (blue).

The initialization strategy for the combinatorial fMINLP implementation
consists of the following steps. In the first step (MIP), a MILP model is
solved to determine an initial network topology. For simplicity, the non-
linearities present in the system are not considered. The optimal network
topology of the MIP step is used in $\mathrm{NLP_{1}}$,$\mathrm{NLP_{2}}$ and
the MINLP as the initial network topology (Visualized by the red arrows in
figure 8). The first NLP model ($\mathrm{NLP_{1}}$) is designed to consider
some of the non-linearities that were ignored in the initial MIP model. In
this step, the characteristic equation of a radiator is used in a simplified
form in order to reduce the computational complexity of the optimization
problem. The second NLP model ($\mathrm{NLP_{2}}$) takes into account the full
set of non-linearities present in the system. Finally, the full MINLP is
solved including all constraints and binary variables.

## References

  * [1] OECD/IEA, “Renewables 2019 Analysis and forecast to 2024,” tech. rep., 2019\.
  * [2] J. Söderman, “Optimisation of structure and operation of district cooling networks in urban regions,” Applied Thermal Engineering, vol. 27, no. 16 SPEC. ISS., pp. 2665–2676, 2007.
  * [3] J. Dorfner and T. Hamacher, “Large-scale district heating network optimization,” IEEE Transactions on Smart Grid, vol. 5, no. 4, pp. 1884–1891, 2014.
  * [4] C. Haikarainen, F. Pettersson, and H. Saxén, “A model for structural and operational optimization of distributed energy systems,” Applied Thermal Engineering, vol. 70, pp. 211–218, sep 2014.
  * [5] W. Mazairac, R. Salenbien, and B. De Vries, “Towards an optimal topology for hybrid energy networks,” EG-ICE 2015 - 22nd Workshop of the European Group of Intelligent Computing in Engineering, 2015.
  * [6] B. Morvaj, R. Evins, and J. Carmeliet, “Optimising urban energy systems: Simultaneous system sizing, operation and district heating network layout,” Energy, vol. 116, pp. 619–636, dec 2016.
  * [7] C. Bordin, A. Gordini, and D. Vigo, “An optimization approach for district heating strategic network design,” European Journal of Operational Research, vol. 252, no. 1, pp. 296–307, 2016.
  * [8] T. Résimont, Q. Louveaux, and P. Dewallef, “Optimization Tool for the Strategic Outline and Sizing of District Heating Networks Using a Geographic Information System,” Energies, vol. 14, p. 5575, sep 2021.
  * [9] H. Lund, P. A. Østergaard, T. B. Nielsen, S. Werner, J. E. Thorsen, O. Gudmundsson, A. Arabkoohsar, and B. V. Mathiesen, “Perspectives on fourth and fifth generation district heating,” Energy, vol. 227, p. 120520, jul 2021.
  * [10] H. Li and S. Svendsen, “District heating network design and configuration optimization with genetic algorithm,” Journal of Sustainable Development of Energy, Water and Environment Systems, vol. 1, no. 4, pp. 291–303, 2013.
  * [11] A. Allen, G. Henze, K. Baker, G. Pavlak, and M. Murphy, “An optimization framework for the network design of advanced district thermal energy systems,” Energy Conversion and Management, vol. 266, 2022.
  * [12] T. Mertz, S. Serra, A. Henon, and J.-M. M. Reneaume, “A MINLP optimization of the configuration and the design of a district heating network: Academic study cases,” Energy, vol. 117, pp. 450–464, 2016.
  * [13] M. Blommaert, Y. Wack, and M. Baelmans, “An adjoint optimization approach for the topological design of large-scale district heating networks based on nonlinear models,” Applied Energy, vol. 280, no. May, p. 116025, 2020\.
  * [14] Y. Wack, M. Baelmans, R. Salenbien, and M. Blommaert, “Economic topology optimization of district heating networks using a pipe penalization approach,” Energy, vol. 264, p. 126161, 2023.
  * [15] F. Marty, S. Serra, S. Sochard, and J.-M. Reneaume, “Simultaneous optimization of the district heating network topology and the Organic Rankine Cycle sizing of a geothermal plant,” Energy, vol. 159, pp. 1060–1074, 2018.
  * [16] Ashrae, “Radiators,” in ASHRAE handbook, ch. Radiators, p. 231, 2014.
  * [17] J. J. Chen, “Comments on improvements on a replacement for the logarithmic mean,” jan 1987.