# Fairness evaluation during the conceptual design of heat grids with quantum
annealers

Kelvin Loh<EMAIL_ADDRESS>0000-0002-8085-6495 TNO - Quantum
TechnologyStieltjesweg 1Delft2628CKthe Netherlands

###### Abstract.

This paper presents a workflow to evaluate the macro scale thermal fairness
for producers in a district heating network during the conceptual design phase
of such a network. It uses the workflow to evaluate two types of proposed
topologies for a future network to be constructed in Europe. The workflow is
also novel in the sense that it has been demonstrated in the paper that within
the implementation, the load balancing step can be solved readily by a quantum
computing system (in particular, a quantum annealer from DWave). The concept
of fairness is also addressed in the workflow and the results show that there
exist an optimum number of producers for a given district heating topology.

district heating, fairness, quantum annealing, qubo, graph partition, load
balancing

††ccs: Mathematics of computing Optimization with randomized search
heuristics††ccs: Applied computing Mathematics and statistics††ccs:
Mathematics of computing Graph coloring††ccs: Applied computing Environmental
sciences††ccs: Hardware Quantum computation

##  1\. Introduction

The concept of using district heating (DH) networks as a method of providing
heat and hot water to households are gaining popularity within Europe. In the
Netherlands, there are now increased initiatives at the local municipalities
(Niessink and Rösler, 2015) and federal government levels to replace gas
networks with 4th generation district heat grids (4GDH). The idea is to use
multiple sources of hot water (such as from geothermal wells, and waste heat
from industrial plants) to provide for the consumers within the network, thus,
reducing or replacing the need of gas as the main source of heating energy.
There have already been many investigations done on the characteristics of a
4GDH ((Rezaie and Rosen, 2012), (Lund et al., 2014), (Lund et al., 2010)).
Amongst the key characteristic is that 4GDH involves strategic and innovative
planning and the integration of DH into the operation of smart energy systems
(Lund et al., 2018). This means that when planning and designing for a DH,
(e.g. during the conceptual design phase) network designers need to take a
holistic approach to come up with robust designs for all the stakeholders
within the DH network. There have been many works on optimization of DH during
the design (Li and Svendsen, 2013) and operation (Sakawa et al., 2001) phases,
however, none has taken into account the concept of fairness during the design
phase (Sarbu et al., 2019).

This paper investigates on the idea of thermal fairness (Bhattacharya et al.,
2016) (in the reference, micro-level fairness is applied to the consumers
instead of producers of the network during operation planning to provide a
fair demand response), which on a macro scale could correspond to a measure of
the tendency for monopolistic business positions to occur within the network.
We claim that the topology of the network itself combined with the consumer
demands in the network would affect the macro fairness of the network (taking
the concept of fairness from computer networks (Denda et al., 2000)), thus
either promoting or discouraging monopoly. The workflow and implementation
presented here could be used in the conceptual design phase for the network
designer to encourage fair and even competition within a DH network. We apply
the concept of fairness to an upcoming DHN in Europe and use actual energy
demands of the consumers in the proposed network in our approach. This paper
also demonstrates the use of a quantum annealer to solve the generic $k$-graph
partitioning problem when we need to split the demands of the consumers fairly
amongst the total number of producers within the grid. We thus claim that the
solution to a balanced $k$-graph is fair for the producers and policy makers
can use this as a measure to either reward or penalize business practices
which deviate from the optimum solution.

###  1.1. DHN Fairness

Typically in the design for 4GDH networks, a monopolistic position of the DHN
operator is always assumed ((Buffa et al., 2019), ((UK), 2018)). This is a
problem often mentioned in the literature (Osman, 2015). The reason this
occurs is because the owner of the network is also responsible for heat
delivery, thus providing no room for competition. In the Dutch heat market
where customers are used to switching energy providers (through electricity
and natural gas networks), the monopolistic position of DHNs is perceived as a
step backwards in terms of fair competitions. In order to protect the
customers of DHNs from the monopolistic position of DHN operators, the
government established the Dutch Heat Act (Osman, 2015). In order to improve
upon this situation further, we need to already include the concept of
fairness during the conceptual design phase of a 4GDH network ((UK), 2018).
This approach is as recommended by the Competition and Markets Authority in
the UK, recommending that network design and build should be improved to
better align the interests of the stakeholders in the network ((UK), 2018).

###  1.2. Quantum Technology

As we are using a relatively new technology to solve the classical
combinatorial optimization problem, we briefly mention about the technology.
There are two types of quantum computing implementations (Preskill, 2018). The
first is what we know as a universal gated quantum computer and the other is a
quantum annealer.

We can solve our problem in gated quantum computers using an algorithm known
as Quantum Approximate Optimization Algorithm (QAOA). This is a hybrid
optimization algorithm introduced by (Farhi et al., 2014) to take advantage of
Noisy Intermediate Scale Quantum (NISQ) (Preskill, 2018) devices to
demonstrate that early quantum computers can already be useful for certain
classical combinatorial problems. However, the number of qubits required to
solve our problem is still too large compared to current publicly available
gated quantum computers. As such, in this paper, we focus on using the
(DW_2000Q) quantum annealer system provided by DWave (Systems, 2017).

Quantum annealing (QA) is the quantum counterpart to simulated annealing. It
takes advantage of quantum mechanical fluctuations instead of thermal, to get
to an optimal solution with faster convergence (Kadowaki and Nishimori, 1998).
The D-Wave system minimizes the following Quadratic Unconstrained Binary
Optimization (QUBO) objective function (equation (1)).

(1) $\displaystyle O(\mathbf{Q},\mathbf{x})$
$\displaystyle=\sum_{i}Q_{ii}x_{i}+\sum_{i<j}Q_{ij}x_{i}x_{j}$

By tuning the input of the objective function $\mathbf{Q}$, we change the
system Hamiltonian for the qubits. The output of the anneal is a low-energy
ground state, which consists of an Ising spin for each qubit in the system.
This is the mechanism upon which a quantum annealer is constructed which gives
the ability to quickly solve certain classes of NP-hard complex problems such
as optimization, machine learning and sampling problems (Ushijima-Mwesigwa et
al., 2017).

As we are trying to solve a graph partition problem in the energy domain, we
are only aware of the work (Ajagekar and You, 2019) in energy systems
optimization using quantum computing systems.

##  2\. Theoretical Background

In this section we describe the theoretical foundations of our implementation.
When we want to evenly distribute a set of nodes in a given network, the
problem is known as a graph partitioning problem. We start by formulating this
problem to a QUBO problem for $k$-clusters (Ushijima-Mwesigwa et al., 2017)
such that we can solve this problem with the DWave system.

Let $G(V,E)$ be the graph which connects the consumers where the nodes are the
consumer nodes and edges represent the proposed connection between the
consumer nodes. We would like to partition this graph between $k$ producers.
The network has $n$ nodes and $m$ edges. Note that the producers are assumed
to have infinite capacities and are always capable of meeting the demands of
the network. Another assumption is that the producers are not fixed in any
location.

We start the theoretical foundation by first deriving the QUBO formulation of
our problem.

###  2.1. Derivation of the QUBO

In order to partition a graph concurrently into an arbitrary number of parts,
$k$, we need to formulate the problem as a QUBO, which takes on binary
variables. The decision variables are given by

(2) $\displaystyle x_{i,j}$ $\displaystyle=\begin{cases}1,&\text{if node
}i\text{ is in producer }j\\\ 0,&\text{otherwise}\end{cases}$

The constraint

(3) $\displaystyle\sum_{j=1}^{k}x_{i,j}$ $\displaystyle=1$

for each node $i$ ensures that each node is in exactly one producer. The
second constraint

(4) $\displaystyle\sum_{i=1}^{n}x_{i,j}$ $\displaystyle=\frac{n}{k}$

for $j=1,...,k$ are the balancing constraints for the producer sizes.

$\displaystyle\mathbf{x}_{j}=\left[\begin{matrix}x_{1,j}\\\ x_{2,j}\\\
\vdots\\\ x_{n,j}\end{matrix}\right]$

The number of cut edges across the $k$ producers is given by

(5)
$\displaystyle\frac{1}{2}\left(\sum_{j=1}^{k}\mathbf{x}_{j}^{\text{T}}L\mathbf{x}_{j}\right)$

where $L$ is the Laplacian matrix of the graph.

The cost function is then given by

(6)
$\displaystyle\beta\left(\sum_{j=1}^{k}\mathbf{x}_{j}^{\text{T}}L\mathbf{x}_{j}\right)+\sum_{j=1}^{k}\alpha_{j}\left(\sum_{i=1}^{n}x_{i,j}-\frac{n}{k}\right)^{2}+\sum_{i=1}^{n}\gamma_{i}\left(\sum_{j=1}^{k}x_{i,j}-1\right)^{2}$

where $\beta$, $\alpha_{i}$, and $\gamma_{i}$ are positive penalty constants.

Following from the derivation and the variables as presented by (Ushijima-
Mwesigwa et al., 2017), the QUBO formulation is

(7) $\displaystyle
min_{\mathbf{x}}\left(\mathbf{X}^{\text{T}}(\beta\mathcal{L}+\alpha\mathcal{I}+\mathbf{B}_{\Gamma})\mathbf{X}-(2\Gamma^{\text{T}}+2\frac{n}{k}\alpha\mathbbm{1}_{\mathcal{N}}^{\text{T}})\mathbf{X}\right)$

where
$\mathbf{X}^{\text{T}}=\left[\mathbf{x}_{1}^{\text{T}},\dots,\mathbf{x}_{k}^{\text{T}}\right]$.

As we have noticed, Equation (6) assumes that the demands are equally weighted
and the edge weights are also equal. We can modify the problem in such a way
that we fulfil different demand constraints and different distances between
the consumer nodes. In order to take into account the different demand
constraints, we need to modify the constraint on the equal parts Eq (4) to

(8) $\displaystyle\sum_{i=1}^{n}w_{i}x_{i,j}$
$\displaystyle=\frac{\sum_{i=1}^{n}w_{i}}{k}$

where $w_{i}$ is the weighted demand of the $i$-th consumer. As for the
condition that each cluster should minimize the total distances between the
member nodes, we can modify Equation (5) to

(9)
$\displaystyle\frac{1}{2}\left(\sum_{j=1}^{k}\mathbf{x}_{j}^{\text{T}}\Delta_{\mathcal{L}}\mathbf{x}_{j}\right)$

where $\Delta_{\mathcal{L}}=$ off-diagonal elements of the matrix
$(-IDI^{\text{T}})$, $I\in\mathbb{R}^{n\times m}$ is the incidence matrix of
the network, $D\in\mathbb{R}^{m\times m}$ is a diagonal matrix with entries
$d_{i,i}={i\text{-th edge distance}}$.

The modified cost function is then given by Equation (10).

(10)
$\displaystyle\beta\left(\sum_{j=1}^{k}\mathbf{x}_{j}^{\text{T}}\Delta_{\mathcal{L}}\mathbf{x}_{j}\right)+\sum_{j=1}^{k}\alpha_{j}\left(\sum_{i=1}^{n}w_{i}x_{i,j}-\frac{\sum_{i=1}^{n}w_{i}}{k}\right)^{2}+\sum_{i=1}^{n}\gamma_{i}\left(\sum_{j=1}^{k}x_{i,j}-1\right)^{2}$

###  2.2. Fairness measures and Heat Grid overall KPI

In order to evaluate whether a network design is capable of being fair and
also to judge whether the optimized solution of equation (10) is indeed fair,
we need to first define what we mean by fairness of a network, and how it
applies to the heat grid with a customized KPI. We consider first the equal
distribution of demands in our KPI. For a system which allocates the demands
to $k$ producers such that the $i$-th producer needs to provide $y_{i}$ of the
total network demand, then the intuitive measures for fairness given a demand
distribution would be through the variance, coefficient of variation or the
min-max ratio (Jain et al., 1984). In this paper, we would use the Jain
fairness index, equation (11) as proposed by (Jain et al., 1984).

(11) $\displaystyle f_{\text{jain}}(y)$
$\displaystyle=\frac{\left[\sum_{i=1}^{k}y_{i}\right]^{2}}{k\sum_{i=1}^{k}y_{i}^{2}}$

The index, $f_{\text{jain}}(y)\in\mathbb{R}_{[0,1]}$, has the value of 1 when
the total network demands are being equally distributed amongst the producers
(thus, the solution is 100% fair). As the disparity increases, the fairness
index decreases, and the solution which favors a few producers would be
indicated by an index which is close to 0.

In addition to the demand being equally distributed, we need to take into
account the distances of the pipes connecting the demand nodes to each
producer. Ideally, each producer groups should have minimal pipe lengths (i.e.
minimal edge lengths) connecting the demand nodes so as to minimize
transportation costs due to a reduction in heat loss. Thus, we also need to
evaluate the distances term which contributes to the KPI. We let $d_{i,j}$ as
the physical distance of the $i$-th path, assigned to the $j$-th producer. In
other words, the $i$-th path is the shortest path which connects two nodes
that are assigned to producer $j$. We then come up with equation (12) to
evaluate the quality of the solution based on the distance criteria alone.

(12) $\displaystyle D$
$\displaystyle=1-\frac{\sum_{j=1}^{k}\sum_{i=1}^{n_{k}}d_{i,j}}{\sum_{i=1}^{N_{c}}d_{i}}$

where $n_{k}$ is the total number of node combinations assigned within the
producer set $k$, and $N_{c}$ is the total number of path combinations for a
complete graph consisting of the demand nodes. The distance index is scaled by
the maximum possible distance covered by all possible combinations of paths so
that the resulting index, $D\in\mathbb{R}_{[0,1]}$. In effect, smaller
clusters within each producer groups would cause the value of $D$ to tend to
1.

Since both $y$ and $d$ from equations (11) and (12) depend on the solution $X$
of equation (10), the KPI to evaluate the solution for a given network is then
simply given by equation (13)

(13) $\displaystyle KPI\left(\mathbf{X}\right)$ $\displaystyle=\alpha_{j}\cdot
f_{\text{jain}}\left(\mathbf{X}\right)+(1-\alpha_{j})\cdot
D\left(\mathbf{X}\right)$

where $\alpha_{j}\in\mathbb{R}_{[0,1]}$ represent the relative importance
factor between the KPI terms. Throughout the paper, the presented cases use an
equal weighting (i.e. $\alpha_{j}=0.5$) unless stated otherwise.

##  3\. Implementation

In this section we describe the implementation in more detail. First we start
by describing how we determine the weighting factor of each node as they need
to represent the actual heat demands inside the proposed network. After that
we present the two proposed networks for which we need to evaluate their
fairness indices with different number of producers in the network. We next
discuss the solvers for which we can obtain the optimized solutions of
equation (10).

###  3.1. Demand nodes

In order for us to represent the actual heat demands in the network, we start
by transforming the energy demands profile of the year 2015 (with $Nt$
discrete number of timesteps) for each node to a static form. We take the max
of the annual demand of each node. This is because the network has to be able
to handle the worst case scenario of providing for the maximum energy demands
from the users of the network. We take the max of the annual demand for each
node as it represents the highest possible energy demand which a particular
node would expect to be delivered by the network producers. We then scale it
by the sum over all the consumers since we would like to have a weight factor
which sums to 1. Let $w_{i,t}$ be the time-dependent demand profile of the
$i$-th consumer node, at timestep $t$ of the profile. The fraction used as
input when solving for equation (10) is then given by equation (14).

(14) $\displaystyle w_{i}$ $\displaystyle=\frac{\max_{t\in
Nt}\left(w_{i,t}\right)}{\sum_{i=1}^{n}\max_{t\in Nt}\left(w_{i,t}\right)}$

The effect of applying equation (14) to the demand profiles is shown in figure
1. The normalized profile in figure 1(a) is translated to figure 1(b) for each
consumer nodes.

(a) Normalized energy demands per node for 2015

(b) Weighted energy demands per node for 2015

Figure 1. Demands profile per node

###  3.2. Topology types

We investigate two types of abstract topologies as these were the ones
considered by the designers of the network. One is a tree-like topology
(figure 2(a)) which connects the consumers while the other adds a few more
connections to make it a circular or ring-like topology (figure 2(b)). It is
important to note that as this methodology evaluates fairness of a network
topology given the demands, we could also use the workflow presented here
within a topology optimizer which produces the topologies for evaluation and
have a more robust final design of a network.

(a) Tree-like topology

(b) Circular-like topology

Figure 2. Topologies considered for evaluation

###  3.3. Solvers

Since the optimization problem with equation (10) is a QUBO problem, we can
attempt to solve this in a quantum annealer which has been described in more
detail in section 1.2. We implemented the QUBO matrix formed by equation (10).
This matrix is then provided to the qbsolv (M. Booth and Roy, 2017) solver
from D-Wave Systems to connect to their QPU machine with the automatic
embedding to the system through the EmbeddingComposite (Systems, 2017)
function.

We also used the METIS (Karypis and Kumar, 1998) graph partitioning solver to
compare a classical solution with the one obtained from qbsolv.

###  3.4. Workflow

The workflow to evaluate fairness is shown in algorithm 1. The results of the
workflow output can help inform the user on the fairness of the evaluated
topology, $G(V,E)$.

input : The topology, $G(V,E)$, to be evaluated

input : Weighted demands from consumers, $w$

input : Maximum number of producers to be considered, $N_{k}$

for  _$k\leftarrow 1$ to $N_{k}$_ do

$\mathbf{X}_{k}\leftarrow$ QUBOSolver ($k$, $w$, $G(V,E)$) /* QUBOSolver
solves eq (10) using qbsolv or METIS */

kpiArray [$k$] $\leftarrow$ KPI($\mathbf{X}_{k}$) /* KPI evaluations using eq
(13) */

end for

output : kpiArray

Algorithm 1 Workflow used to evaluate the topology KPIs

##  4\. Results

In this section we provide the results of our proposed topologies on the
fairness and overall KPI. We also discuss on the quality of the solutions
produced by both solvers.

Figure 3 show the different fairness indices when we increase the number of
producers as computed by the solvers for different topologies. We notice from
the Jain fairness index (fig 3(a)) that for all the topologies, the index
reduces as we increase the number of producers. As for the distance index (fig
3(b)), it indicates that as the number of producers increase in the network,
the energy transportation costs reduce, thus increasing the index value.

The overall KPI indicates the quality of the solution provided by the solvers
and we can see from figure 3(c) that the optimum quality of the solution
occurs when the number of producers in the network (for both topologies) is
approximately 5. As we change the weighting factor $\alpha_{j}$ of the overall
KPI (equation (13)), we would switch the overall KPI between that of fully
fair (fig 3(a)) to minimization of cluster distance (fig 3(b)).

(a) Jain fairness index

(b) Distance index

(c) Overall fairness KPI

Figure 3. Fairness indices between multiple topologies, number of producers,
and solvers

Another observation from figure 3(c) is that the circular-like topology
produces better quality assignments (figure 4) compared to the tree-like
topology (figure 5) from both solvers. As for the comparison between the
different solvers, we do not notice any significant difference between the
solutions produced by both types of solvers, and this is expected since the
number of consumer nodes in the network is relatively small, and from
(Ushijima-Mwesigwa et al., 2017), we would not notice much added advantage
until we get to networks of size, $n\gtrapprox 200$. However, the same
implementation and workflow can be used for more topologies and network sizes
in the future.

(a) Node assignments for each producer (Map)

(b) Node assignments for each producer

Figure 4. qbsolv solution for $k=5$ and a circular type topology

(a) Node assignments for each producer (Map)

(b) Node assignments for each producer

Figure 5. qbsolv solution for $k=5$ and a tree type topology

##  5\. Conclusions

From the previous section, we have seen that a heat grid network can have an
optimum number of producers. Adding more producers would reduce the fairness
or the overall quality of the assignments. Thus, another conclusion which we
can draw from the results is that during the conceptual design phase, the
fairness concept should be introduced to help provide a robust network for the
future. Trivially, adding more connections between nodes should provide an
increased fairness for all the producers within the network. However, this
does not take into account the original costs of laying the connections in the
first place. We also note that it is up to the user to interpret the KPI and
how it can be used to alleviate the monopolistic problem of DHNs. The KPI just
serves as a quantifiable measure of fairness which can be tuned.

The workflow can also be adapted to be more time dynamic by performing the
fairness optimization step in a quasi-static manner (i.e. by evaluating the
solution at each timestep for the nodes instead of currently taking the
maximum of the demand profile).

On a technological note, the implementation is suitable to be solved in a
quantum annealer through qbsolv, and the advantages could be promising if the
network to be evaluated consist of a large number of edges and vertices. We
could also implement the same QUBO problem in gated quantum computers and use
the QAOA hybrid algorithm to solve the optimization problem.

##  6\. Acknowledgements

The author would like to acknowledge the useful discussions he had with his
fellow colleagues, Paul Egberts, and Can Tumer at the Heat Transfer and Fluid
Dynamics (HTFD) department in TNO, regarding topology optimization and
district heating network fairness. The author would also like to acknowledge
the support from the department manager, Jeremy Veltin, in allowing for this
project to be completed.

## References

  * (1)
  * Ajagekar and You (2019) Akshay Ajagekar and Fengqi You. 2019.  Quantum computing for energy systems optimization: Challenges and opportunities.  _Energy_ 179 (2019), 76 – 89.  https://doi.org/10.1016/j.energy.2019.04.186
  * Bhattacharya et al. (2016) Saptarshi Bhattacharya, Vikas Chandan, Vijay Arya, and Koushik Kar. 2016.  Thermally-fair Demand Response for District Heating and Cooling (DHC) Networks. In _Proceedings of the Seventh International Conference on Future Energy Systems_ _(e-Energy ’16)_. ACM, New York, NY, USA, Article 15, 11 pages.  https://doi.org/10.1145/2934328.2934343
  * Buffa et al. (2019) Simone Buffa, Marco Cozzini, Matteo D’Antoni, Marco Baratieri, and Roberto Fedrizzi. 2019\.  5th generation district heating and cooling systems: A review of existing cases in Europe.  _Renewable and Sustainable Energy Reviews_ 104 (2019), 504 – 522.  https://doi.org/10.1016/j.rser.2018.12.059
  * Denda et al. (2000) Robert Denda, Albert Banchs, and Wolfgang Effelsberg. 2000\.  The Fairness Challenge in Computer Networks. In _Quality of Future Internet Services_ , Jon Crowcroft, James Roberts, and Mikhail I. Smirnov (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, 208–220.
  * Farhi et al. (2014) Edward Farhi, Jeffrey Goldstone, and Sam Gutmann. 2014\.  A Quantum Approximate Optimization Algorithm.  _arXiv e-prints_ 1411, 4028 (Nov 2014).  arXiv:quant-ph/1411.4028
  * Jain et al. (1984) R. Jain, D. Chiu, and W. Hawe. 1984.  _A Quantitative Measure of Fairness and Discrimination for Resource Allocation in Shared Computer Systems_.  DEC Research Report TR-301. Digital Equipment Corporation, Maynard, MA, USA. 38 pages.  ftp://ftp.netlab.ohio-state.edu/pub/jain/papers/fairness.htm
  * Kadowaki and Nishimori (1998) T. Kadowaki and H. Nishimori. 1998.  Quantum annealing in the transverse Ising model.  _Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics_ 58, 5 (1998), 5355–5363.  www.scopus.com Cited By :529.
  * Karypis and Kumar (1998) George Karypis and Vipin Kumar. 1998.  A Fast and High Quality Multilevel Scheme for Partitioning Irregular Graphs.  _SIAM J. Sci. Comput._ 20, 1 (Dec. 1998), 359–392.  https://doi.org/10.1137/S1064827595287997
  * Li and Svendsen (2013) H. Li and S. Svendsen. 2013\.  District heating network design and configuration optimization with genetic algorithm.  _Journal of Sustainable Development of Energy, Water and Environment Systems_ 1, 4 (2013), 291–303.  www.scopus.com Cited By :22.
  * Lund et al. (2010) H. Lund, B. Möller, B.V. Mathiesen, and A. Dyrelund. 2010\.  The role of district heating in future renewable energy systems.  _Energy_ 35, 3 (2010), 1381 – 1390.  https://doi.org/10.1016/j.energy.2009.11.023
  * Lund et al. (2014) Henrik Lund, Sven Werner, Robin Wiltshire, Svend Svendsen, Jan Eric Thorsen, Frede Hvelplund, and Brian Vad Mathiesen. 2014.  4th Generation District Heating (4GDH).  _Energy_ 68, C (2014), 1–11.  https://doi.org/10.1016/j.energy.2014.02.
  * Lund et al. (2018) Henrik Lund, Poul Alberg Østergaard, Miguel Chang, Sven Werner, Svend Svendsen, Peter Sorknæs, Jan Eric Thorsen, Frede Hvelplund, Bent Ole Gram Mortensen, Brian Vad Mathiesen, Carsten Bojesen, Neven Duic, Xiliang Zhang, and Bernd Möller. 2018\.  The status of 4th generation district heating: Research and results.  _Energy_ 164 (2018), 147 – 159.  https://doi.org/10.1016/j.energy.2018.08.206
  * M. Booth and Roy (2017) S.P. Reinhardt M. Booth and A. Roy. 2017.  _Partitioning Optimization Problems for Hybrid Classical/Quantum Execution_.  D-Wave Technical Report Series 14. D-Wave Systems. 1–9 pages.  http://github.com/dwavesystems/qbsolv
  * Niessink and Rösler (2015) Robin Niessink and Hilke Rösler. 2015.  _Development of Heat Distribution Networks in the Netherlands_.  Technical Report. ECN.  https://www.ecn.nl/publicaties/PdfFetch.aspx?nr=ECN-E--15-069
  * Osman (2015) Nejat Osman. 2015\.  _Barriers to district heating development in the Netherlands: a business model perspective_.  Master Thesis. University of Twente.  https://essay.utwente.nl/73054/1/Osman_BA_BMS.pdf
  * Preskill (2018) John Preskill. 2018\.  Quantum Computing in the NISQ era and beyond.  _arXiv e-prints_ 1801, 00862 (Jan 2018).  arXiv:quant-ph/1801.00862
  * Rezaie and Rosen (2012) Behnaz Rezaie and Marc A. Rosen. 2012.  District heating and cooling: Review of technology and potential enhancements.  _Applied Energy_ 93 (2012), 2 – 10.  https://doi.org/10.1016/j.apenergy.2011.04.020 (1) Green Energy; (2)Special Section from papers presented at the 2nd International Enery 2030 Conf.
  * Sakawa et al. (2001) M Sakawa, K Kato, S Ushiro, and M Inaoka. 2001\.  Operation planning of district heating and cooling plants using genetic algorithms for mixed integer programming.  _Applied Soft Computing_ 1, 2 (2001), 139 – 150.  https://doi.org/10.1016/S1568-4946(01)00014-X
  * Sarbu et al. (2019) Ioan Sarbu, Matei Mirza, and Emanuel Crasmareanu. 2019\.  A review of modelling and optimisation techniques for district heating systems.  _International Journal of Energy Research_ (05 2019).  https://doi.org/10.1002/er.4600
  * Systems (2017) D-Wave Systems. 2017\.  _D-Wave Initiates Open Quantum Software Environment_.  Press release. D-Wave Systems.  https://www.dwavesys.com/press-releases/d-wave-initiates-open-quantum-software-environment
  * (UK) (2018) Competition & Markets Authority (UK). 2018\.  _Heat networks market study: Final report_.  Technical Report. UK Government.  https://assets.publishing.service.gov.uk/media/5b55965740f0b6338218d6a4/heat_networks_final_report.pdf
  * Ushijima-Mwesigwa et al. (2017) Hayato Ushijima-Mwesigwa, Christian F. A. Negre, and Susan M. Mniszewski. 2017.  Graph Partitioning using Quantum Annealing on the D-Wave System.  _arXiv e-prints_ 1705, 03082 (May 2017).  arXiv:quant-ph/1705.03082