# Benchmarking Quantum(-inspired) Annealing Hardware on Practical Use Cases

##  Appendix A Datasets

All datasets in this paper are available from
https://github.com/ianmalcolm/annealer-benchmark. We are adding README.md,
ground truth and python scripts to improve the usability of the repo.

##  Appendix B Backends

A brief introduction of the backends are as follows.

D-Wave is a world-leading company that designs and builds quantum annealers.
The Quantum Processing Unit (QPU) in the D-Wave’s quantum annealing system is
a lattice of interconnected qubits. Each qubit is made of a superconducting
loop. D-Wave released the D-Wave 2000Q system in 2017, the QPU in which
employs a Chimera graph architecture Vert et al. (2019), equipped with 2048
qubits and 6016 couplers. Due to the variation in the manufacturing, the
actual working resource in a quantum annealer varies.

Figure 1: Chimera graph architecture Inc. . Circles are qubits. A black
straight lines represents an internal coupler. The red curved lines are the
external couplers.

Fig.1 shows a part of the _Chimera_ graph architecture. Each black circle
represents a qubit. Every eight qubits form a complete bipartite graph. The
connections within a bipartite graph are internal couplers. The connections
across bipartite graphs are external couplers. A connection between two qubits
corresponds to a quadratic term of two decision variables in a QUBO problem.

D-Wave released the Advantage system in 2019. One of the major differences
between these two systems is the graph architecture of their QPUs. The
Advantage system QPU employs a _Pegasus_ Inc.  graph architecture, which has
5640 qubits and 40,484 couplers. Generally speaking, we expect a quantum
annealer to solve larger and denser problems if it has more qubits and
couplers.

The Fujitsu Digital Annealer (DA) Tsukamoto et al. (2017) is a hardware
implementation of an enhanced variant of the Simulated Annealing (SA)
algorithm. The enhancements include but are limited to parallel-trial, dynamic
offset and Parallel Tempering (PT) with Isoenergetic Cluster Moves (ICM) Zhu
et al. (2020). The workflow of DA is shown in Fig.2 Aramon et al. (2019).

Initialise bit arrayUpdateSwapBit flip & evaluateAccept
RejectTerminate?UpdateSwap?ReturnYesNoNoYes

Figure 2: Fujitsu Digital Annealer Workflow Diagram

DA employs Metropolis-Hastings update with Parallel Tempering. In each cycle,
it adjusts the temperature, proposes bit flips, and accepts/rejects the
proposals according to the temperature and periodically swaps configurations
between systems. Apart from the sophisticated techniques in proposing and
selecting updates, the competence of DA also comes from the efficient hardware
implementation that exploits the parallelism in the algorithm. We include DA
in the evaluation.

We include the D-Wave QPU of Chimera and Pegasus architecture, and Fujitsu
Digital Annealer in the evaluation. We also include Simulated Annealing (SA)
by D-Wave, which is an open-source CPU implementation based on the Metropolis-
Hastings algorithm. We include QBSolv from D-Wave, which is a heuristic hybrid
solver that incorporates classical computer and quantum annealer. We configure
QBSolv to use Tabu search and use it as a software baseline in the
experiments. We also include Gurobi 9.1, which represents the SOTA commercial
optimisation solver. When we generate a synthetic dataset for the experiments,
we use Gurobi to find the approximated optimal solutions for reference.

Table 1: Solver settings

Solver |  Chimera |  Pegasus |  DA |  QBSolv |  SA
---|---|---|---|---|---
Implementation |  D-Wave DW_2000Q_61 |  D-Wave Advantage System 1.12 |  Fujitsu Quantum-inspired Computing Digital Annealer On-Premises Service |  D-Wave Ocean Software3 |  D-Wave Ocean Software4
Number of samples |  100 |  100 |  128 replicas |  20 |  100
Annealing process |  Default annealing schedule5 with annealing time ranging from 2 to 2000 $\mu s$, annealing offsets disabled, 1000us programming thermalisation and 0s readout thermalisation.6 |  Same as Chimera |  Parallel tempering, with 128 replicas. The number of iterations ranges from $10^{5}$ to $10^{9}$. Offset increase rate is disabled 7 |  Default solver Tabu search |  Problem specific and geometry beta scheduling. Total number of sweeps ranges from $10^{2}$ to $10^{4}$ 8
Precision of QUBO coefficient |  Rescale fp32 and operates in anologue mode |  Same as Chimera |  Cast fp64 to 64/76-bit integer for quadratic/linear term for problem size below 4096 qubits, or to 16/76-bit integer for problem size over 4096 qubits |  fp32 |  fp32
Postprocess |  Disabled |  Disabled |  Disabled |  N.A. |  N.A.
Random seed |  N.A. |  N.A. |  N.A. |  1234 |  1234
Initial state |  None |  None |  None |  None |  None

  * 1

AWS Braket Device ARN: arn:aws:braket:::device/qpu/d-wave/DW_2000Q_6

  * 2

AWS Braket Device ARN: arn:aws:braket:::device/qpu/d-wave/Advantage_system1

  * 3

https://github.com/dwavesystems/qbsolv

  * 4

https://github.com/dwavesystems/dwave-neal

  * 5

QPU-specific anneal schedules documents:
https://docs.dwavesys.com/docs/latest/doc_physical_properties.html#doc-qpu-
characteristics

  * 6

Complete D-Wave QPU parameters:
https://docs.dwavesys.com/docs/latest/c_solver_parameters.html

  * 7

Complete DA annealing parameters:
https://portal.aispf.global.fujitsu.com/apidoc/da/jp/api-ref/da-qubo-en.html

  * 8

Temperature scheduling for SA https://github.com/dwavesystems/dwave-
neal/blob/a11e477c3a6b3585d75ee8b58be75a4127d0c17c/neal/sampler.py#L281

The settings of the solvers are presented in Table 1. We set the _number of
samples_ (a.k.a number of shots, solutions, or readout) per task #samples=100
for Chimera, Pegasus, SA. The annealing process is crucial to the quality of
solutions. We are following the default annealing schedule provided by D-Wave,
which empirically works on most of problems. We carry out the experiments with
different _annealing time_ , which determines the duration of annealing
schedule, to understand if a problem is sensitive to it. For Chimera and
Pegasus, we change the annealing time from 2 to 2000 $\mu s$. This range
almost cover the full range allowed by the devices, which is 1 to 2000 $\mu
s$.

In the implementation of SA, number of sweeps, or #sweeps, is a hyper-
parameter that is similar to the annealing time in D-Wave annealer. we change
#sweeps from $10^{2}$ to $10^{4}$. $10^{3}$ is the default setting, which lead
to 1 sweep for each annealing temperature beta. #sweeps higher than $10^{3}$
causes multiple sweeps per beta, lower causes skip of some beta.

For DA, we use Parallel Tempering mode with #replicas=128. The #replicas=128
means DA runs 128 copies of a problem instance in parallel, each of the
instances is randomly initialized, at different temperatures. At the end of
annealing, we get 128 samples from DA. This is different from the rest of the
backends, where a problem instance is sampled repeatedly, in a sequential way.
The #replicas corresponds to the parallelism of DA. The upper limit of
#replicas is 128. We fix #replicas=128 to make sure DA has its advantage in
the comparison. The number of iterations, or #iterations, is a hyper-parameter
similar to the annealing time in D-Wave annealer. We are varying #iterations
from $10^{5}$ to $10^{9}$. $10^{5}$ is the lower limit of #iterations.
Experimental results suggest that #iterations higher than $10^{9}$ does not
make big difference in terms of the quality of solutions. Some problems are
simple enough to be solved with fewer #iterations. Having too large
#iterations puts DA in disadvantage in terms of speed.

We use randomly generated initial state for all backends. For SA and QBSolv,
we use random seed=1234 to improve the reproducibility.

In the main text, we only include the results with the best hyper-parameters
in the comparison. By saying best, we refer to the following rules: If a
hyper-parameter outperforms others in terms of objective energy, it is the
best. If all hyper-parameters are similar energy-wise, then the one with
lowest time cost is the best. An optimal hyper-parameter on one problem
setting is not necessarily optimal on other problems settings.

###  B.1 Gurobi

The problem formulation for Gurobi goes like this. For max-cut problem, we
formulate its QUBO form and solve it on Gurobi. For MVC and QAP, we formulate
the objective as QUBO form and implement the problem constraints using Gurobi
constraints. QUBO form is not in favour of Gurobi. There are other smart
classical encoding methods for Gurobi. For QAP for example, instead of using
one-hot encoding, one can use categorical variables to represent locations and
factories and solve the problem more efficiently. But we choose QUBO form on
Gubori because it provides a better comparability between annealing-based
binary solvers and classical heuristic solvers.

The solving time of Gurobi presented in the comparison plot of the main text
is the time for traversing the whole solution space, if not terminated by a
pre-defined timeout threshold. It is possible that a global optima is
discovered way before the solution space is traversed. Gurobi program
maintains a “best-so-far” solution during the search. We plot time-to-solution
for Gurobi, to understand how fast Gurobi approaches global optimal solution.

The Gurobi program is running on a server, which is equipped with an Intel
Core i9-10900X CPU, 128GB DDR4 memory and 128GB HDD swap memory. In our
experimental settings, Gurobi always instantiate 20 threads and occupies
almost all CPU time slices, DDR4 memory and swap memory, until it finishes the
search. We pay efforts to minimise the interference from irrelevant processes
on the same server, because some problem instances in this paper are really
pushing the server toward limit and cause Gurobi to terminate with
“OUT_OF_MEMORY” errors.

##  Appendix C max-cut

###  C.1 Pegasus-like max-cut problems

Problem generation We start the problem generation by first checking the
working status of D-Wave QPU. The actual architecture of a QPU could be
different from the description in the official documentation. This is because
there is variation in manufacturing that puts some resources of the QPU in
non-working conditions. In our case, we are accessing D-Wave from AWS Braket
Gonzalez (2021). There are 5436 qubits and 37440 couplers available on D-Wave
Advantage System 1.1. By accessing DWaveSampler.target_structure, one can find
the graph, $G_{\textrm{pegasus}}=(V_{\textrm{pegasus}},E_{\textrm{pegasus}})$
representing the architecture of the target QPU.

Based on $G_{\textrm{pegasus}}$, we generate random graphs with a specified
number of nodes. Given $G_{\textrm{pegasus}}$ and a specific number of nodes
$n$, we randomly choose $\left|V_{\textrm{pegasus}}\right|-n$ nodes from
$G_{\textrm{pegasus}}$, remove these nodes and corresponding edges. The
randomness in the graph generation is controlled by a random seed for good
reproducibility. The resulting graph is a sub-graph of $G_{\textrm{pegasus}}$
with $n$ nodes. We can map this graph onto QPU without a resort to minor
embedding effort.

In the pegasus-like max-cut experiment, we have ten random graphs of such.
$\left|V\right|$ of these graphs ranges from 543 to 5430.

At the late stage of work, AWS Braket D-Wave Advantage System 1.1 is retired,
and is replaced with 4.1 and 6.1. Although they follow exactly the same
architecture, the generated problems based on 1.1 are not applicable to 4.1
and 6.1, because some qubits or couplers that were available on 1.1 are not
necessarily in working condition on 4.1 or 6.1.

(a) Pegasus

(b) DA

(c) SA

Figure 3: max-cut, Pegasus-like problems. The error bar of energy and timing.
X axis is the number of nodes in a graph. The primary Y axis is energy
normalised to that of Gurobi. The secondary Y axis corresponds to dotted
curves and represents timing in seconds, in log scale. The circle of an error
bar represents mean energy, while the top and bottom represent $\pm$ standard
deviation. We shift the error bars a little bit to improve readability and
avoid overlap. Since max-cut is a maximisation problem, higher energy is
better.

Figure 3 shows the performance of Pegasus, DA and SA on max-cut Pegasus-like
problems. For D-Wave Pegasus, figure 3(a) suggest that annealing time of 2000
$\mu s$ achieves higher energy and outperforms shorter annealing time.
Adiabatic theory Born and Fock (1928); Tong et al. (2005) suggests that longer
annealing time will ensure a slower adiabatic process and a better result. In
practice McGeoch (2014); Hauke et al. (2020), D-Wave annealers operate in non-
adiabatic mode. An open quantum system de-coheres and suffers from thermal
noises, which sometimes could also improve the performance of a quantum
system. In our case, we cannot distinguish if the improvement comes from the
quantum mechanics or the thermal noises. We include annealing time of 2000
$\mu s$ for D-Wave Pegasus in the main text for comparison between solvers.

Figure 3(b) shows the performance of DA on the Pegasus-like max-cut problem.
The runs with iterations larger than $10^{6}$ are relatively consistent in
terms of the variance of energy. Further increase of the iterations does not
make prominent difference. We present the results from #iterations=$10^{6}$ in
the main text and compare it with other solvers. DA with $10^{5}$ iteration on
large problem size leads to $\pm 10\%$ variance in energy. We are not going
investigate the reason, since in-depth analysis of a solver and its underlying
algorithm is out of the scope of this work.

Figure 3(c) shows the performance of SA on the Pegasus-like max-cut problem.
The performance are relatively consistent on problems of different size. But
higher number of sweeps provides significantly better results. We include the
results of #sweeps=$10^{4}$ in the main text for comparison.

Figure 4: Gurobi on max-cut Pegasus-like problem. Time to solution. X axis is
the time in seconds in log scale. Y axis is the energy normalised to optimal
solution (found by Gurobi)

Figure 4 shows the time-to-solution of Gurobi on max-cut Pegasus-like
problems. Although some large problem instances, e.g. the one with 5430 nodes,
take quite a long time before the termination of the optimisation program
(because of reaching pre-configured timeout), Gurobi can find very promising
solutions within roughly 100 seconds.

Please note that a point in the figure marks the time Gurobi spends to
discover the promising solution at that moment. Gurobi is not sure if the
discovered the solution is globally optimal. It has to spend some more time to
traverse the whole solution space with some smart pruning techniques to make
sure it is not possible to find any better solutions. The time for traversing
the solution space is reported in Table 2.

Table 2: Gurobi time-to-traverse on Pegasus-like max-cut problems #nodes | Time (s) | #nodes | Time (s)
---|---|---|---
543 | 0.041288 | 3258 | 36000.05
1086 | 0.25132 | 3801 | 36000.03
1629 | 14.66663 | 4344 | 36000.41
2172 | 430.4116 | 4887 | 36000.05
2715 | 1131.357 | 5430 | 36000.05

To collect the intermediate results of Gurobi, we have to pass a callback to
Gurobi program, which introduce extra time cost. For a problem instance that
last for 36000 seconds, the time spent on the callback is over 1000 seconds.

###  C.2 Chimera-like max-cut problems

The chimera-like graphs are randomly generated, following the D-Wave Chimera
architecture. The generation is similar to that of the Pegasus-like max-cut
problem. Please refer to section C.1 for the details of the problem
generation.

(a) Chimera

(b) DA

(c) SA

Figure 5: max-cut, Chimera-like problems. The error bar of energy and timing.
X axis is the number of nodes in a graph. The primary Y axis is energy
normalised to that of Gurobi. The secondary Y axis corresponds to dotted
curves and represents timing in seconds, in log scale. The circle of an error
bar represents mean energy, while the top and bottom represent $\pm$ standard
deviation. We shift the error bars a little bit to improve readability and
avoid overlap. Since max-cut is a maximisation problem, higher energy is
better.

Figure 5 shows the performance of D-Wave Chimera and the Digital Annealer on
the max-cut Chimera-like problems. Through figure 5(a) we understand that a
longer annealing time, i.e. 2000 $\mu s$, leads to the best performance in
energy. But 2000 $\mu s$ does not give Chimera prominent advantage over that
of 200 $\mu s$. This result coincides with our observation with the Pegasus
architecture. The performance of quantum and classical annealers in terms of
energy on max-cut Chimera-like problems are generally improved given higher
hyper-parameters. But for DA, the room of improvement is limited.

For Chimera, we only include the results of $200\mu s$ for comparison between
solvers. For DA and SA We only include the results of #iterations=$10^{6}$ and
#sweeps=$10^{4}$ for comparison between solvers.

(a) Timing

(b) Normalised energy

Figure 6: Max-cut, Chimera-like problems. Comparison between solvers. a) X
axis is the number of nodes. Y axis is the timing of the solvers, in log
scale. b) X axis is the number nodes. Y axis is the energy normalised to that
of Gurobi. We include mean energy (dotted lines) and best energy (solid lines)
to represent the statistics of energy.

Figure 6 shows the performance of different solvers on max-cut Chimera-like
problems. From figure 6(a) we understand that Chimera is the fastest solver
for the most of the time. Gurobi is faster than classical annealing-based
solvers and a little bit slower than Chimera. In terms of energy, according to
figure 6(b) Chimera is the most promising solver among other annealing-based
solvers. The gap between Chimera and DA is widening, as the problem size
increases.

Figure 7: Gurobi on max-cut Chimera-like problems. Time to solution. X axis is
the time in seconds in log scale. Y axis is the energy normalised to optimal
solution (found by Gurobi)

Figure 7 shows the performance of Gurobi on max-cut Chimera-like problems. In
this figure, most of the promising solutions can be found within 0.1 seconds.
The time-to-traverse information is listed in Table 3

Table 3: Gurobi time-to-traverse on Chimera-like max-cut problems #nodes | Time (s) | #nodes | Time (s)
---|---|---|---
204 | 0.088003159 | 1224 | 0.047677994
408 | 0.050300837 | 1428 | 0.087327957
612 | 0.054126978 | 1632 | 0.083067894
816 | 0.057450056 | 1836 | 0.103567839
1020 | 0.067600012 | 2040 | 0.116798878

###  C.3 Connectivity-varied max-cut problems

Problem generation We use the “dense_gnm_random_graph” method from networkx,
i.e., the open-source python library, to generate random graphs. We set
$n=145$, as this is the largest complete graph that can be mapped onto a
Pegasus QPU. We sweep the averaged degree to create graphs with different
connectivity. For example, to generate a graph with an averaged degree of 100,
we set $m=\mathrm{floor}(145\times 100\div 2)$, which represents the total
number of edges in the graph.

(a) Pegasus

(b) DA

(c) SA

Figure 8: Max-cut, connectivity-varied problems. The error bar of energy and
timing. X axis is the averaged degrees, in log scale. The primary Y axis is
energy normalised to that of Gurobi. The secondary Y axis corresponds to
dotted curves and represents timing in seconds, in log scale. The circle of an
error bar represents mean energy, while the top and bottom represent $\pm$
standard deviation. We shift the error bars a little bit to improve
readability and avoid overlap. Since max-cut is a maximisation problem, higher
energy is better.

Figure 8 shows the performance of the solvers on the max-cut Chimera-like
problems. Through figure 8(a) we understand that a longer annealing time, i.e.
2000 $\mu s$, leads to the best performance in energy, when the averaged
degree is below 12. But with higher averaged degrees, the advantage disappear.
DA is not sensitive to its hyper-parameter. SA is more sensitive to its hyper-
parameter, compared with DA.

For Pegasus, we only include the results of $2000\mu s$ in the main text for
comparison between solvers. For DA and SA We only include the results of
#iterations=$10^{6}$ and #sweeps=$10^{4}$ in the main text for comparison
between solvers.

Figure 9: Gurobi on max-cut connectivity-varied problem. Time to solution. X
axis is the time in seconds in log scale. Y axis is the energy normalised to
optimal solution (found by Gurobi)

Figure 4 shows the time-to-solution of Gurobi on max-cut connectivity-varied
problems. Gurobi can find very promising solutions within roughly 10 seconds.
The time-to-traverse information is listed in Table 4.

Table 4: Gurobi time-to-traverse on connectivity-varied max-cut problems #degrees | Time (s) | #degrees | Time (s)
---|---|---|---
1 | 0.004468 | 17 | 19396.73
2 | 0.034891 | 18 | 36000.02
3 | 0.051368 | 19 | 36000.02
4 | 0.252771 | 20 | 36000.01
5 | 0.319412 | 30 | 36000.03
6 | 3.782623 | 40 | 36000.01
7 | 4.946026 | 50 | 36000.11
8 | 7.411871 | 60 | 36000.13
9 | 42.18519 | 70 | 36000.1
10 | 224.94 | 80 | 36000.33
11 | 348.8753 | 90 | 36000.26
12 | 674.4734 | 100 | 36000.18
13 | 1015.575 | 110 | 36000.79
14 | 1766.733 | 120 | 36000.79
15 | 3245.935 | 130 | 36000.14
16 | 9300.061 | 140 | 36000.14

##  Appendix D mvc

###  D.1 MVC Pegasus-like problems

We reuse the graphs generated for Section C.1 and formulate mvc problem based
on these graphs.

(a) Error bar of energy and timing

(b) Probability of feasibility

Figure 10: SA on Pegasus-like MVC problems. a) The error bar of energy and
timing. X axis is the number of nodes in a graph. The primary Y axis is energy
normalised to that of Gurobi. The secondary Y axis corresponds to dotted
curves and represents timing in seconds, in log scale. The circle of an error
bar represents mean energy, while the top and bottom represent $\pm$ standard
deviation. We shift the error bars a little bit to improve readability and
avoid overlap. Since MVC is a minimisation problem, lower energy is better. b)
Probability of feasibility. Y axis represents the percentage of samples that
meet the constraints of the problem.

Figure 10 shows the performance of SA on Pegasus-like MVC problems. For SA,
large #sweeps translates to better performance and longer time cost. This is
consistent across a range of problem size and hyper-parameter settings.

(a) Error bar of energy and timing

(b) Probability of feasibility

(c) Constraint violation

Figure 11: D-Wave Pegasus on Pegasus-like MVC problems. a) The error bar of
energy and timing. X axis is the number of nodes in a graph. The primary Y
axis is energy normalised to that of Gurobi. The secondary Y axis corresponds
to dotted curves and represents timing in seconds, in log scale. The circle of
an error bar represents mean energy, while the top and bottom represent $\pm$
standard deviation. We shift the error bars a little bit to improve
readability and avoid overlap. Since MVC is a minimisation problem, lower
energy is better. b) Probability of feasibility. Y axis represents the
percentage of samples that meet the constraints of the problem. c) Constraint
violations over a range of problem sizes and hyper-parameters. Each mark is a
solution. X axis represents the percentage of violated constraints in a
problem. Y axis represents pseudo energy, normalised to that of Gurobi.

Figure 11 shows the performance of D-Wave Pegasus on Pegasus-like MVC
problems. From figure 11(a) we understand that longer annealing time usually
produces the best solutions. From figure 11(b) we can also see that longer
annealing time produces the best probability of feasibility. Although Pegasus-
like MVC problems have exactly the same graph topology as those of Pegasus-
like max-cut problems, the optimal setting of the annealing time is on the
different end of the allowed range. We reach to the conclusion that the
optimal value of annealing time is closely related to the problem settings.
However, hyper-parameter optimisation is not the focus of this paper. We only
pick the hyper-parameters with the most promising results to be included for
comparison between solvers.

Figure 11(c) aims to compare between feasible and infeasible solutions, over a
range of problem sizes and hyper-parameter settings. If you have a look at the
first row, where the annealing time is 2 $\mu s$, as the problem size is
getting larger, the cluster of solutions is moving away from 0% violation. The
pseudo energy of the solutions are lowering, meaning that the feasibility is
traded for lower pseudo energy. If you have a look at the last row, where the
annealing time is 2000 $\mu s$, the situation is the other way around. The
pseudo energy is higher on larger problems, The lower pseudo energy is traded
for higher feasibility.

We know that quantum mechanics is not the only driving power for problem
optimisation in a quantum annealer Hauke et al. (2020). Other things, like
thermal noise, also play important roles in the optimisation process. It would
be interesting if we can identify which part is dominating the performance in
Pegasus-like max-cut and Pegasus-like MVC problems, and find a guideline on
the proper setting of the annealing time.

(a) Error bar of energy and timing

(b) Probability of feasibility

(c) Constraint violation

Figure 12: DA on Pegasus-like MVC problems. The setting of the plots is the
same as that of figure 11

Figure 12 shows the performance of DA on Pegasus-like MVC problems. Figure
12(a) and 12(b) suggest that when #iterations is small, e.g. $10^{5}$, there
is large variance in energy and low probability of feasibility. Figure 12(c)
suggest that small #iterations, like $10^{5}$, cannot effectively explore the
solution space and leave some of the solution infeasible. the feasibility is
traded for low pseudo energy.

For Pegasus, we only include the results of $2000\mu s$ for comparison between
solvers. For DA and SA We only include the results of #iterations=$10^{7}$ and
#sweeps=$10^{4}$ for comparison between solvers.

(a) Timing

(b) Normalised energy

Figure 13: Energy and timing of solvers on Pegasus-like MVC problems. We
choose $T_{anneal}=2000\mu s$ for D-Wave Pegasus, #iterations=$10^{5}$ for DA
and #sweeps=$10^{4}$ for SA. The setting of the figure is the same as that of
figure 6.

Figure 13 shows the energy and timing of the solvers on Pegasus-like MVC
problems. D-Wave Pegasus is faster than the other annealing-based solvers, but
is worst in terms of the energy of solutions. DA is faster than SA, and is the
best among all annealing-based solvers in terms of energy.

(a) Probability of feasibility

(b) Constraint Violation

Figure 14: Feasibility/violation of solvers on Pegasus-like MVC problems. The
hyper-parameter setting is the same as that in figure 13. a) X axis is the
number of nodes. Y axis is percentage of samples that meet the constraints of
the problem. b) constraint violation on pegasus_node5430. X axis is the
percentage of constraints being violated. Y axis is the pseudo energy,
normalised to that of Gurobi.

Figure 14 shows the probability of feasibility, as well as the constraint
violation, of solvers on Pegasus-like MVC problems. All solvers except Pegasus
can find feasible solutions effectively on all problem instances. The curves
of probability of feasibility are mostly overlapped with that of Gurobi in
figure 14(a). Figure 14(b) suggests that there is not much violation in the
samples of Pegasus. One can fix the violations and get solutions of quality
that is close to SA.

Figure 15: Gurobi time-to-solution on Pegasus-like MVC problems

Figure 15 shows the Gurobi time-to-solution plot. Gurobi can find promising
solutions within a few seconds. The Gurobi time-to-traverse on Pegasus-like
MVC problems is listed in Table 5

Table 5: Gurobi time-to-traverse on Pegasus-like MVC problems #nodes | Time | #nodes | Time
---|---|---|---
543 | 0.002709 | 3258 | 23.13427
1086 | 0.014983 | 3801 | 42.36129
1629 | 0.095287 | 4344 | 699.0575
2172 | 0.864816 | 4887 | 135.1517
2715 | 2.327012 | 5430 | 49.6165

###  D.2 MVC connectivity-varied problems

We reuse the graphs generated for Section C.3 and formulate mvc problem based
on these graphs.

(a) Error bar of energy and timing

(b) Probability of feasibility

Figure 16: SA on Pegasus-like MVC problems. The setting of the plots is the
same as that of figure 11, except that we skip the plot of constraint
violation, as the solutions returned by SA are all feasible.

Figure 16 shows the performance of SA on connectivity-varied MVC problems. For
SA, large #sweeps translates to better performance and longer time cost. This
is consistent across a range of problem size and hyper-parameter settings.

(a) Error bar of energy and timing

(b) Probability of feasibility

(c) Constraint violation

Figure 17: D-Wave Pegasus on connectivity-varied MVC problems. The plot
setting is the same as that in figure 11

Figure 17 shows the performance of D-Wave Pegasus on connectivity-varied MVC
problems. From figure 17(a) we understand that longer annealing time usually
produces the best solutions, but the annealing time of 20 $\mu s$ produces
better probability of feasibility. From figure 11(b) we known that the
difference between 20$\mu s$ and 2000$\mu s$ is trivial. This observation
suggest that optimality and feasibility can be detached.

In figure 11(c) first row last column sub-plot, the marker spreads widely and
are mostly close to the bottom of the plot. This suggest that shorter
annealing time ensure a lower (pseudo) energy, but is more likely to produce
infeasible solutions. This is opposite to the sub-plot of last row first
column, where markers has higher (pseudo) energy but closer to the 0%
violation.

(a) Error bar of energy and timing

(b) Probability of feasibility

(c) Constraint violation

Figure 18: DA on connectivity-varied MVC problems. The plot setting is the
same as that in figure 11

Figure 18 shows the performance of DA on connectivity-varied MVC problems.
This set of plots are similar to those in figure 12, except that the problem
size in terms of number of nodes is much smaller and trivial to DA. DA can
solve these problems more quickly.

For Pegasus, we only include the results of $2000\mu s$ in the main text for
comparison between solvers. For DA and SA We only include the results of
#iterations=$10^{5}$ and #sweeps=$10^{2}$ in the main text for comparison
between solvers.

Figure 19: Gurobi time-to-solution on connectivity-varied MVC problems

Figure 19 shows the Gurobi time-to-solution plot. Gurobi can find promising
solutions within a second. The Gurobi time-to-traverse on Pegasus-like MVC
problems is listed in Table 6

Table 6: Gurobi time-to-traverse on connectivity-varied MVC problems #degrees | Time | #degrees | Time
---|---|---|---
1 | 0.0009 | 17 | 0.5694
2 | 0.0014 | 18 | 0.4882
3 | 0.0061 | 19 | 0.9114
4 | 0.0128 | 20 | 0.4998
5 | 0.0252 | 30 | 0.9279
6 | 0.0806 | 40 | 0.9257
7 | 0.0722 | 50 | 3.0588
8 | 0.1579 | 60 | 0.3693
9 | 0.5501 | 70 | 0.1204
10 | 0.967 | 80 | 1.1727
11 | 0.4539 | 90 | 1.0164
12 | 0.4778 | 100 | 0.4307
13 | 0.4948 | 110 | 0.9458
14 | 1.6128 | 120 | 0.4478
15 | 0.523 | 130 | 0.1147
16 | 0.4513 | 140 | 0.0948

###  D.3 MVC Benchmark from DIMACS 10th Challenge

The graphs in the experiment are extracted from DIMACS 10th Challenge
(https://www.cc.gatech.edu/dimacs10/index.shtml). Authors with interests can
go to Github and find useful course works on the MVC benchmark in CSE-6410
from Georgia Tech. Students of this course have made a few comprehensive
benchmarks of heuristic methods on these MVC problems. A few examples are
listed below:

  * •

https://github.com/Z-Jiang/CS-6140

  * •

https://github.com/sangyh/minimum-vertex-cover

  * •

https://github.com/ChujieChen/Minimum-Vertex-Cover

  * •

https://github.com/xwave7/minimum-vertex-cover

  * •

https://github.com/arjunchint/Minimum-Vertex-Cover

(a) Error bar of energy and timing

(b) Probability of feasibility

Figure 20: SA on DIMACS10 MVC problems. The setting of the plots is the same
as that of figure 11, except that we skip the plot of constraint violation, as
the solutions returned by SA are all feasible.

Figure 20 shows the performance of SA on DIMACS 10th Challenge MVC problems.
For SA, large #sweeps translates to longer time cost, but make no big
difference in terms of energy.

(a) Error bar of energy and timing

(b) Probability of feasibility

(c) Constraint violation

Figure 21: D-Wave Pegasus on DIMACS10 MVC problems

Figure 21 shows the performance of D-Wave Pegasus on DIMACS 10th Challenge MVC
problems. As the topology of the graphs in the dataset varies largely, we
cannot observe clear trend from figure 21(a) and 21(b). But we can at least
understand that Pegasus is able to find very limited number of feasible
solutions, when the problem size is relatively. Annealing time as short as 2
$\mu s$ can hardly produce feasible solutions for graphs with over 100 nodes.
This observation match with the previous synthetic datasets.

From figure 11(c) we understand there is limited amount of constraint
violation in the solutions, as most of clusters are very close to the 0%
violation line. With a naive fixation step, one can save values from those
broken solutions.

(a) Error bar of energy and timing

(b) Probability of feasibility

(c) Constraint violation

Figure 22: DA on DIMACS10 MVC problems

Figure 22 shows the performance of DA on DIMACS 10th Challenge MVC problems.
This set of plots are similar to those in figure 18, except that the problem
size in terms of number of nodes is much smaller and trivial to DA. DA can
solve these problems more quickly.

For Pegasus, we only include the results of $2000\mu s$ in the main text for
comparison between solvers. For DA and SA We only include the results of
#iterations=$10^{6}$ and #sweeps=$10^{2}$ in the main text for comparison
between solvers.

Figure 23: Gurobi time-to-solution

Figure 23 shows the Gurobi time-to-solution plot. Gurobi can find promising
solutions within a second.

The “dimacs10-netscience” problem describes a network of co-authorships in the
area of network science. It has 1461 nodes and 2742 edges. The more detailed
description of this problem can be found via
http://konect.cc/networks/dimacs10-netscience/.

We use the benchmark from DIMACS 10th Challenge in 2012, because we believe it
closely follows the trend of requirements nowadays. However, most of the
previous works on heuristic MVC methods are based on the benchmark from DIMACS
2nd Challenge in 1992. We evaluate this outdated benchmark as well.

###  D.4 MVC Benchmark from DIMACS 2nd Challenge

The MVC dataset from DIMACS 2nd Challenge
http://archive.dimacs.rutgers.edu/pub/challenge/graph/benchmarks/volume/Clique/
is one of the most widely adopted Wang et al. (2019); Cai et al. (2013); Li et
al. (2020) benchmarks for evaluating MVC related algorithms, despite of the
fact that it is a 30-ish years old benchmark. We do not include quantum
annealer in this experiment because most of the problem instances in the
benchmark is too large or dense for quantum annealers. We set
#iterations=$10^{8}$ for DA and #sweeps=$10^{3}$ for SA.

Table 7: Performance on MVC benchmark from DIMACS 2nd Challenge. The performance is shown in Energy/Timing (s)/$P_{f}$ format. Name | #nodes/#edges | Optimal | da | qbsolv | sa
---|---|---|---|---|---
brock200_2 | 200 / 9876 | 188 | 189/27.033/1 | 189/4.29/1 | 192/2.213/1
brock200_4 | 200 / 13089 | 183 | 192/27.025/1 | 192/1.384/1 | 193/2.398/1
brock400_2 | 400 / 59786 | 371 | 392/27.316/1 | 392/3.149/1 | 395/10.057/1
brock400_4 | 400 / 59765 | 367 | 393/27.317/1 | 393/2.544/1 | 395/10.093/1
brock800_2 | 800 / 208166 | 776 | 790/27.846/1 | 790/8.71/1 | 793/34.695/1
brock800_4 | 800 / 207643 | 774 | 790/27.834/1 | 790/9.48/1 | 792/35.435/1
C125.9 | 125 / 6963 | 91 | 121/26.924/1 | 121/0.665/1 | 122/1.308/1
C250.9 | 250 / 27984 | 206 | 245/27.107/1 | 245/1.659/1 | 245/4.678/1
C500.9 | 500 / 112332 | 443 | 495/27.468/1 | 495/7.144/1 | 496/18.693/1
C1000.9 | 1000 / 450079 | 932 | 994/28.183/1 | 995/14.26/1 | 996/74.469/1
C2000.5 | 2000 / 999836 | 1984 | 1983/32.724/1 | 1985/42.001/1 | 1989/178.126/1
C2000.9 | 2000 / 1799532 | 1920 | 1994/32.946/1 | 1995/55.965/1 | 1995/318.634/1
C4000.5 | 4000 / 4000268 | 3982 | 3982/49.541/1 | 3985/147.304/1 | 3988/730.413/1
DSJC500.5 | 500 / 62624 | 487 | 487/27.457/1 | 487/4.353/1 | 490/10.871/1
DSJC1000.5 | 1000 / 249826 | 985 | 985/28.144/1 | 986/11.126/1 | 990/43.229/1
gen200_p0.9_44 | 200 / 17910 | 156 | 195/27.028/1 | 195/1.01/1 | 196/3.105/1
gen200_p0.9_55 | 200 / 17910 | 145 | 195/27.03/1 | 195/1.33/1 | 196/3.109/1
gen400_p0.9_55 | 400 / 71820 | 345 | 392/27.314/1 | 392/2.408/1 | 392/11.949/1
gen400_p0.9_65 | 400 / 71820 | 335 | 393/27.303/1 | 393/2.653/1 | 394/12.024/1
gen400_p0.9_75 | 400 / 71820 | 325 | 394/27.315/1 | 394/2.599/1 | 395/11.878/1
hamming8-4 | 256 / 20864 | 240 | 240/27.113/1 | 240/1.098/1 | 243/3.69/1
hamming10-4 | 1024 / 434176 | 984 | 1004/28.229/1 | 1004/13.756/1 | 1008/71.61/1
keller4 | 171 / 9435 | 160 | 156/26.995/1 | 156/0.808/1 | 159/1.797/1
keller5 | 776 / 225990 | 749 | 745/27.805/1 | 745/8.611/1 | 754/38.329/1
keller6 | 3361 / 4619898 | 3302 | 3298/47.664/1 | 3298/159.914/1 | 3318/835.81/1
MANN_a27 | 378 / 70551 | 252 | 375/27.277/1 | 375/2.402/1 | 375/11.581/1
MANN_a45 | 1035 / 533115 | 690 | 1032/30.363/1 | 1032/16.445/1 | 1032/88.17/1
MANN_a81 | 3321 / 5506380 | 2221 | 3318/47.809/1 | 3318/172.374/1 | 3318/971.542/1
p_hat300-1 | 300 / 10933 | 292 | 261/27.179/1 | 261/1.868/1 | 271/2.329/1
p_hat300-2 | 300 / 21928 | 275 | 273/27.181/1 | 273/1.816/1 | 281/4.106/1
p_hat300-3 | 300 / 33390 | 264 | 291/27.168/1 | 291/2.674/1 | 293/5.786/1
p_hat700-1 | 700 / 60999 | 689 | 635/27.798/1 | 637/4.187/1 | 660/11.511/1
p_hat700-2 | 700 / 121728 | 656 | 651/27.807/1 | 653/6.29/1 | 668/21.366/1
p_hat700-3 | 700 / 183010 | 638 | 690/27.708/1 | 690/8.831/1 | 692/30.694/1
p_hat1500-1 | 1500 / 284923 | 1488 | 1413/31.966/1 | 1417/13.831/1 | 1449/52.967/1
p_hat1500-2 | 1500 / 568960 | 1435 | 1438/31.598/1 | 1440/21.961/1 | 1462/100.844/1
p_hat1500-3 | 1500 / 847244 | 1406 | 1488/31.725/1 | 1489/29.186/1 | 1492/147.645/1

##  Appendix E qap

###  E.1 TinyQAP problems

Problem generation We start the problem generation by randomly generating $n$
points in a two-dimensional space using “np.random.rand” from NumPy with a
random seed of 1234. The $n$ points are saved for reserved usage. We then
generate the distance adjacency matrix from the list of points. For the
generation of the flow matrix, we generate $F_{n\times n}=(m+m.T)/2$, where
$m$ is an $n\times n$ matrix generated by the same random method. We set the
diagonal of $F_{n\times n}$ to be zero to eliminate loops.

We evaluate Chimera, Pegasus, DA and SA on this dataset. Specifically, we
restrict Chimera on problem of size below 8, since this is the largest QAP it
can handle.

(a) Error bar of energy and timing

(b) Probability of feasibility

(c) Constraint violation

Figure 24: D-Wave Chimera on Tinyqap. The plot setting is the same as that in
figure 11

Figure 24 shows the performance of D-Wave Chimera on Tinyqap problems. Problem
of size larger than 8 is not included because they are beyond the capacity of
Chimera. From figure 24(a) we understand that longer annealing time does not
have clear advantage over shorter ones in terms of energy as well as
feasibility. This is also supported by figure 24(b), which suggest that
Chimera can hardly find any feasible solutions when the problem size is larger
than 5. Figure 24(c) shows that there is high level of constraints violation.
The sub-plots in the upper right corner and lower left corner suggest there
could be as much as 50%-100% constraints being violated in a solution. This
translate to more efforts in fixing the broken constraints, compared with
those solutions in MVC problems.

(a) Error bar of energy and timing

(b) Probability of feasibility

(c) Constraint violation

Figure 25: D-Wave Pegasus on Tinyqap. The plot setting is the same as that in
figure 11

Figure 25 shows the performance of D-Wave Pegasus on Tinyqap problems.
Comparing with figure 24 we known that Pegasus has no clear advantage either
in terms of energy of in feasibility. The solutions by Pegasus have high level
of constraint violations, which also similar to that of Chimera.

(a) Error bar of energy and timing

(b) Probability of feasibility

(c) Constraint violation

Figure 26: DA on Tinyqap. The plot setting is the same as that in figure 11

Figure 25 shows the performance of DA on Tinyqap problems. Figure 26(a)
suggest that higher #iterations can produce better results. This is consistent
across all problem size. Figure 26(b) suggest that DA can find feasible
solutions in most cases. #iterations as low as $10^{5}$ only suffers a few
percent loss in probability of feasibility. Figure 26(c) suggest that the
infeasible cases has much higher level of constraint violations, compared with
DA’s performance on MVC problems.

(a) Error bar of energy and timing

(b) Probability of feasibility

(c) Constraint violation

Figure 27: SA on Tinyqap. The plot setting is the same as that in figure 11

Figure 27 shows the performance of SA on Tinyqap problems. Figure 27(a)
suggest that higher #sweeps can produce slightly better results. This is
consistent across all problem size. Figure 27(b) suggest that SA generally has
difficulty in finding feasible solutions. Higher #sweeps improves probability
of feasibility. On all MVC problems SA can always find feasible solutions, but
figure 27(c) suggest that the infeasible cases have up to 30% constraint
violations.

For D-Wave annealer, we only include the results of $200\mu s$ for Chimera and
$2000\mu s$ for Pegasus in the main text for comparison between solvers. For
DA and SA We only include the results of #iterations=$10^{8}$ and
#sweeps=$3\times 10^{3}$ in the main text for comparison between solvers.

Figure 28: Gurobi time-to-solution on TinyQAP problems

Figure 28 shows the Gurobi time-to-solution plot. Gurobi can find promising
solutions within a few seconds. The Gurobi time-to-traverse on TinyQAP is
listed in Table 8

Table 8: Gurobi time-to-traverse on TinyQAP problems Size | Time | Size | Time
---|---|---|---
3 | 0.025959 | 8 | 1.3042
4 | 0.037391 | 9 | 2.283769
5 | 0.055867 | 10 | 30.85997
6 | 0.095988 | 11 | 217.3864
7 | 0.52411 | 12 | 1443.717

###  E.2 Tai Benchmark from QAPLIB

Tai Benchmark is available from https://coral.ise.lehigh.edu/data-
sets/qaplib/qaplib-problem-instances-and-solutions/#Ta

(a) Error bar of energy and timing

(b) Probability of feasibility

(c) Constraint violation

Figure 29: DA on Tai of QAPLIB. The plot setting is the same as that in figure
11

Figure 29 shows the performance of DA on Tai problems from QAPLIB. Problem of
size larger than 80 is not included because they are beyond the capacity of
DA. From figure 29(a) we understand that higher #iterations produces better
solutions in terms of energy. This is also supported by figure 29(b), which
suggest that lower #iterations has difficulty in finding feasible solutions
Figure 29(c) also support our observation in figure 29(b).

(a) Error bar of energy and timing

(b) Probability of feasibility

(c) Constraint violation

Figure 30: SA on Tai of QAPLIB. The plot setting is the same as that in figure
11

Figure 30 shows the performance of SA on Tai problems from QAPLIB. From figure
30(a) and 30(b) we understand that higher #sweeps produces better solutions in
terms of energy and feasibility. SA generally has difficulty in finding
feasible solutions when problem size is over 60. Figure 30(c) also support our
observation in figure 30(b).

For Pegasus, we only include the results of $2000\mu s$ in the main text for
comparison between solvers. For DA and SA We only include the results of
#iterations=$10^{8}$ and #sweeps=$3\times 10^{3}$ in the main text for
comparison between solvers.

We do not run Gurobi over Tai dataset because we use ground truth from QAPLIB.

##  Appendix F Warehouse

###  F.1 Warehouse Assignment Heuristics

Various storage assignment policies have been developed. For example,
randomized storage policy (RSP) Petersen (1997), assign items to storage
randomly. Dedicated storage policy (DSP) Fumi et al. (2013) always assigns
items to their own locations in the warehouse. This means every item has its
dedicated location. This policy then relies on order pickers’ memory power to
minimise the travel distance. Cube-per-order Index (COI) Malmborg and
Bhaskaran (1990) assigns items according to their ratio between volume and
popularity. For the interest of this project, we consider all items to have
identical volume, so this reduces to arranging items according to popularity.
Order-oriented-swapping (OOS) Mantel et al. (2007) It starts with a random
assignment and runs a certain number of iterations. In each iteration, two
items in the assignment are swapped. The distance is recalculated and the new
assignment is accepted if it reduces the distance, and is accepted with some
probability otherwise. Effectively, this is the technique of Simulated
Annealing (SA). class-based storage policy (ABC) Petersen et al. (2004),
groups SKUs into classes according to popularity and allocates a dedicated
area for each class. The number of popularity classes and the number of items
in each class are arbitrary. The QAP formulation bridges the gap between the
warehouse assignment problem and various sophisticated solvers. Following
this, there are successful attempts in implementing QAP Ruijter (2007) for
warehouse assignment and favourable results are observed.

###  F.2 Computing platforms

_Annealing-based computers_ We include _Chimera_ Quantum Annealer (QA) and
Fujitsu Digital Annealer (DA) in the experiment. The hyper-parameters are the
same as that in Section B, except we use $\\#\text{shots}=10000$ for QA to
increase the chance of finding feasible solutions. When problem size is above
90, we use our decomposition heuristic plus DA to solve the problem. Only
annealer:solve() in Algorithm 6 runs on the annealing hardware. The rest of
the algorithm runs on the CPU of a server. The server has an Intel Xeon Silver
4116 at 2.1GHz with 128GB DDR4 memory and 960GB SSD storage.

_Baselines_ We include three classical methods based on CPU, i.e. _QBSolv_
from D-Wave, Simulated Annealing Raidl et al. (2020) (SA) and Tabu Gasquez
(2015), as baselines. SA (python) and Tabu (C++) are specialised for
permutation optimization, meaning that the neighbouring function for QAP is
implemented in programming language, such that constraints do not have to be
incorporated in the objective function. We also include random solutions in
the comparison, which serves as the bottom line of a QAP solver.

###  F.3 Warehouse Dataset

Block-structural QAP instances of sizes 8, 90, 180, 270, 3600, and 8100 are
synthesised using randomly generated order sets. In addition, for objective
2), input data distribution is perturbed to provide modified versions of
datasets of size 270. This is to match the findings in Tsige (2013) that
interaction-based methods such as QAP work well when 80% of ordered items are
concentrated in 20% of the SKUs. In other words, there is a small set of
commonly ordered products. The following assumptions are made about the order
sets as well as the shape of the respective warehouses. WH-270b below is for
objective 2), whereas others are for objective 1).

Table 9: Warehouse dataset Name | Problem Size | #subsets | # rows | #columns
---|---|---|---|---
WH-8 | 8 | 1 | 4 | 2
WH-90 | 90 | 1 | 45 | 2
WH-180 | 180 | 2 | 45 | 4
WH-270 | 270 | 3 | 45 | 6
WH-270b | 270 | 3 | 45 | 6
WH-3600 | 3600 | 40 | 45 | 80
WH-8100 | 8100 | 90 | 45 | 180

In this dataset, the smallest size of the problem instance is 8, which is the
maximum problem size directly handled by chimera. 90 is the biggest problem
size supported by DA. Problem size bigger than 90 is divided into groups of 90
for decomposition. For example, the order set of size 3600 can be solved as 40
sub-QAPs, each of size 90. All except size 8 warehouses are assumed to
comprise columns of size 45, and there is an aisle running in between columns.

To ensure a fair comparison between the heuristic and the software library,
critical parameters are set as identical. In particular, both are set to run
for 10,000,000 iterations, which is the determining factor of the speed of any
simulated annealing algorithm. Note that in the heuristic, our hardware has to
perform 10,000,000 iterations for each of the sub-QAPs. Other annealing
parameters such as maximum temperature and temperature interval are left to
the default settings.

The results obtained for each of the datasets are average values over a
certain number of runs. The software library runs 3 times for each dataset.
The number of runs for the decomposition heuristic varies across datasets due
to overhead. For WH-8, WH-90, WH-180, and WH-270 the heuristic is run thrice,
so the readings are averages of 3. For WH-3600 and WH-8100 the heuristic is
run once.

###  F.4 Formal description

We can think of the the distance matrix and frequency matrix as complete
graphs. The nodes of the graphs represent items and locations themselves.
Formally,

$G=(V(G),E(G))\;|V(G)|=n$

is the graph of items.

$H=(V(H),E(H))\;|V(H)|=n$

is the graph of locations. Define the edge weights of $G$ to be interaction
frequencies, and those of $H$ to be distances. Formally,

$f:E\left(G\right)\rightarrow\mathbb{R}\qquad f(i,j)=f_{ij}$
$d:E\left(H\right)\rightarrow\mathbb{R}\qquad d(i,j)=d_{ij}$

An assignment is given by a bi-jection $\varphi:V(G)\rightarrow V(H)$. The
goal is to solve the following minimisation:

$\min_{\varphi}\sum_{(i,j)\in E(G)}f(i,j)d(\varphi(i),\varphi(j))$ (1)

One can check that the objective function of QAP and 1 are equivalent.

###  F.5 Decomposition

Next, we describe the decomposition formally and provide proof of the
theoretical boundary. In order to minimise the travelling distance of an order
picker, intuitively, we want to maximise the interaction frequencies of items
within a block. This objective can be formally described as follows:

$\text{Maximise}\quad\sum_{i,j=1}^{n}\sum_{l=1}^{k}f_{ij}x_{il}x_{jl}$ (2)

Equation 2 denotes the sum of interaction frequencies among items within their
respective subsets. $(x_{ij})$ is an $n\times k$ decision matrix. $x_{ij}$
denotes whether item $i$ goes to subset $j$. The objective comes with the
following constraints:

$\displaystyle\sum_{l=1}^{k}x_{il}=1\quad\forall i,1\leq i\leq n$ (3)
$\displaystyle\sum_{i=1}^{n}x_{il}\leq s\quad\forall l,1\leq l\leq k$ (4)

Equation 3 means an item can belong to exactly one subset. Note that each
subset can have at most $s=\frac{n}{k}$ items. Therefore, equation 4 means
each subset must not exceed its capacity $s$.

Note that equation 3 is actually another QAP with a particular decision matrix
(but with different constraints). Thus, it can be translated into a graph
formulation. In the graph formulation, $(x_{ij})$ is equivalent to a function
$g:V(G)\to K$, which maps an element of $V(G)$ into its subset, such that
$g(i)=a\iff x_{ia}=1$. This function is well-defined due to constraint
equation 3. Equation 2 can be converted to the equivalent graph formulation as
follows:

$\displaystyle\sum_{i,j=1}^{n}\sum_{l=1}^{k}f_{ij}x_{il}x_{jl}$
$\displaystyle=\sum_{l=1}^{k}\sum_{\begin{subarray}{c}i,j=1\\\
x_{il}=x_{jl}=1\end{subarray}}^{n}f_{ij}$
$\displaystyle=\sum_{l=1}^{k}\sum_{\begin{subarray}{c}i,j=1\\\
g(i)=g(j)=l\end{subarray}}^{n}f_{ij}$
$\displaystyle=\sum_{\begin{subarray}{c}(i,j)\in E(G)\\\
g(i)=g(j)\end{subarray}}f(i,j)$ (5)

which, intuitively, is the sum of flows within subsets. This is the setup for
the theorem below, which states that with a solution to partitioning problem,
an optimal solution to the overall QAP objective can be constructed.

######  Theorem 1.

Let $d(x,y)=\begin{cases}\delta&c(x)=c(y)\\\ M&c(x)\neq c(y)\end{cases}$ for
some positive constants $\delta,M$. Then there exists a $\varphi_{0}$ for
which the warehouse objective is achieved, and $\varphi_{0}$ can be
constructed from a solution to partition objective denoted by $g_{0}$.

###### Proof.

Suppose $\varphi:V(G)\to V(H)$ is some assignment. Then the warehouse
objective can be written as:

$\displaystyle\qquad\sum_{(i,j)\in E(G)}f(i,j)d(\varphi(i),\varphi(j))$
$\displaystyle=\sum_{\begin{subarray}{c}(i,j)\in E(G)\\\
c(\varphi(i))=c(\varphi(j))\end{subarray}}f(i,j)d(\varphi(i),\varphi(j))$
$\displaystyle\qquad+\sum_{\begin{subarray}{c}(i,j)\in E(G)\\\
c(\varphi(i))\neq c(\varphi(j))\end{subarray}}f(i,j)d(\varphi(i),\varphi(j))$
$\displaystyle=\delta\left(\sum_{\begin{subarray}{c}(i,j)\in E(G)\\\
c(\varphi(i))=c(\varphi(j))\end{subarray}}f(i,j)\right)+M\left(\sum_{\begin{subarray}{c}(i,j)\in
E(G)\\\ c(\varphi(i))\neq c(\varphi(j))\end{subarray}}f(i,j)\right)$
$\displaystyle=\delta\mathcal{A}+M(\mathcal{F}-\mathcal{A})$
$\displaystyle=M\mathcal{F}+(\delta-M)\mathcal{A}$

where $\mathcal{F}=\sum_{(i,j)\in E(G)}f(i,j)$, which is the sum of
interaction frequencies between all items, and
$\mathcal{A}=\sum_{\begin{subarray}{c}(i,j)\in E(G)\\\
c(\varphi(i))=c(\varphi(j))\end{subarray}}f(i,j)$ is the sum of interaction
frequencies between items within partitions. Since $\delta-M<0$, the original
objective is minimised if and only if $\mathcal{A}$ is maximised.

Note that $\varphi:V(G)\to V(H)$ and $c:V(H)\to K$, therefore we can view
$c\circ\varphi:V(G)\to K$ and $\mathcal{A}=\sum_{\begin{subarray}{c}(i,j)\in
E(G)\\\ c\circ\varphi(i))=c\circ\varphi(j))\end{subarray}}f(i,j)$. Then
$g_{0}:V(G)\to K$ is by definition the solution to

$\max_{g}\sum_{\begin{subarray}{c}(i,j)\in E(G)\\\
g(i)=g(j)\end{subarray}}f(i,j)$

which means

$\sum_{\begin{subarray}{c}(i,j)\in E(G)\\\
g_{0}(i)=g_{0}(j)\end{subarray}}f(i,j)\leq\mathcal{A},\quad\forall\varphi$

Now we need to construct a $\varphi$ such that
$c(\varphi(i))=c(\varphi(j))\iff g_{0}(i)=g_{0}(j)$. This is equivalent to
saying that after splitting items into subsets according to $g_{0}$, assign
the items such that items within the same subset are in the same interval. The
relative position of items within an interval does not matter. Therefore,
there are many possible answers. ∎

The above method assumes a strong condition that $(d_{ij})$ only has two
distinct values, $\delta$ for locations within an interval and $M$ in between
intervals. This is a simplification of the block structure and $(d_{ij})$, in
reality, is usually more complex.

However, the concept of intervals can be generalised such that when $\delta$
and $M$ are non-constant, a function $c$ can still be defined on $V(H)$ such
that the sum of distances within intervals is minimised. The motivation for
such construct is that when a warehouse does not exhibit a clear-cut column
structure, one can still think of an abstract “column” as a group of close
locations. Note that this definition automatically specialises to the previous
definition of $c$ when $\delta$ and $M$ are constant. This intuition can be
formally expressed as

$\text{Minimise}\quad\sum_{i,j=1}^{n}\sum_{l=1}^{k}d_{ij}x_{il}x_{jl}$ (6)

Subject to:

$\displaystyle\sum_{l=1}^{k}x_{il}=1\quad\forall i,1\leq i\leq n$ (7)
$\displaystyle\sum_{i=1}^{n}x_{il}\leq s\quad\forall l,1\leq l\leq k$ (8)

$(x_{ij})$ is an $n\times k$ decision matrix. $x_{ij}$ denotes whether item
$i$ goes to location group $j$. Note that this graph partitioning problem can
be thought of as the dual to equation 2, with maximisation changed to
minimisation. In effect, items that are frequently ordered together will be
assigned to locations that are closer together.

Now both the set of locations and the set of items have been divided into
subsets of equal size. Between subsets of locations, the distances are
maximal, and between subsets of items, the interaction frequency is minimal.
The final step is to produce a bijection from the set of subsets of items to
the set of subsets of items, and for each pair of subsets in the bijection
$(items,locs)$, there is a sub-QAP of size $O(\sqrt{n})$, and therefore a QUBO
of size $O(n)$, where $n$ is the total number of items.

Note that this is subject to $n$ being a perfect square; in practical
situations where $n$ is not a perfect square, compromise has to be made in
either finding the nearest smaller perfect square and do the optimisation on
the smaller set of items and locations only, or use integer divisors of $n$
other than $\sqrt{n}$, and correspondingly deal with sub-QAPs of sizes other
than $\sqrt{n}$. For example, if there are $n=3600$ locations, it is possible
to divide it into 40 groups of 90, where each sub-QAP will have $n=90$.

###  F.6 QAP-to-QUBO conversion

A QAP has to be converted to QUBO before it is solvable by a quantum annealer.
This is also true of each of the sub-QAPs generated by the decomposition
procedure . Here the procedures of QAP-to-QUBO conversion is described.
Readers familiar with square penalty may safely skip this section.

For the linear constraint of the QAP objective, there is one constraint for
each $i$. Thus, for each $i$ the quadratic penalty term
$P(\sum_{k=1}^{n}x_{ik}-1)^{2}$ is added to the original objective function,
where $P$ is some positive constant to be determined a-posteriori. Note that
this penalty term, when expanded, is also a quadratic form, and thus can be
represented as an addition to the QUBO matrix $Q$.

We can express this constraint in a more general form. For the decision vector
$\bar{x}$, locate the position of $x_{ik}$ and instantiate a 0-1 bit-vector
with length $n$, with the corresponding positions having 1 and the other
positions having 0. For all $n$ constraints, collect the respective bit-
vectors as rows of a $n^{2}\times n^{2}$ matrix $A$ (note that not all rows of
A has to be used as a constraint, in which case the row is left as zero). All
constraints are therefore in the form $A\bar{x}=\bar{b}$, where $\bar{b}$ is
just a length $n^{2}$ vector of all 0’s except $n+n$ 1’s (corresponding to the
$n+n$ constraints we have) in the corresponding rows of $A$. To convert the
constraints into a sum of squares we take the inner product

$(Ax-b)^{T}(Ax-b)$

We can embed the penalty $P$ easily into $A$ and $b$ by scaling. Subsequently,
the expanded QUBO coefficient matrix can be obtained by matrix manipulation.
It could be shown that the final coefficient matrix to be augmented is
$(A^{T}A-2D)$, where $D=diag(b^{T}A)$. The complete algorithm is described
below.

Algorithm 1 Conversion from QAP to QUBO

1://Given: FindIndex$(i,j)$ that maps quadratic index of $x_{ij}$ to its
corresponding position in the decision vector $\bar{x}$

2:procedure convert($F$,$D$,$P$)

3: $Q=\textrm{ZeroMatrix}(n^{2}\times n^{2})$

4: for $i$ in $[1..n]$ do

5: for $j$ in $[1..n]$ do

6: for $k$ in $[1..n]$ do

7: for $l$ in $[1..n]$ do

8: $x_{ik}=\textrm{FindIndex}(i,k)$

9: $x_{jl}=\textrm{FindIndex}(j,l)$

10: $Q[x_{ik}][x_{jl}]=F[i][j]\times D[k][l]$

11: end for

12: end for

13: end for

14: end for

15: $Q^{\prime}=$AtoQ(preparematrixA(), preparevectorB(), $P$)

16: return $Q+Q^{\prime}$

17:end procedure

Algorithm 2 Prepare Maxtrix A

1:

2:procedure preparematrixA

3: $A=\textrm{ZeroMatrix}(n^{2}\times n^{2})$

4: $i=1$

5: for $r$ in $[1..n]$ do

6: for $k$ in $1..n$ do

7: $idx=\textrm{FindIndex}(i,k)$

8: $A[r][idx]=1$

9: end for

10: $i=i+1$

11: end for

12: $k=1$

13: for $r$ in $[(n+1)..2n]$ do

14: for $i$ in $1..n$ do

15: $idx=\textrm{FindIndex}(i,k)$

16: $A[r][idx]=1$

17: end for

18: $k=k+1$

19: end for

20: return $A$

21:end procedure

Algorithm 3 Prepare Vector B

1:

2:procedure preparevectorB

3: $b=\textrm{ZeroVector}(n^{2})$

4: for $r$ in $[1..2n]$ do

5: $b[r]=1$

6: end for

7: return $b$

8:end procedure

Algorithm 4 Convert Matrix A to QUBO

1:

2:procedure AtoQ($A$,$b$,$P$)

3: Scale each row of $A$ by $\sqrt{P}$

4: Scale each entry of $b$ by $\sqrt{P}$

5: $D=b^{T}A$

6: return $A^{T}A-2D$

7:end procedure

In order to enable quantum annealing of a QAP, the problem must first be
converted into a Quadratic Unconstrained Binary optimisation (QUBO) form,
which effectively means constraints have to be subsumed as part of the
objective function in some way. It is customary to do this by encoding the
constraints as quadratic penalty terms which augment the objective function.
It can be shown that for QUBOs, the optimal solution to the augmented
objective function also minimises the original objective function.

###  F.7 Exterior penalty method

We use exterior penalty method to iteratively solve a QUBO as a series of
QUBO’s with increasing penalty weights. Such is called the exterior penalty
method, whose algorithm is described below. Readers familiar with this method
in optimisation may safely skip this section.

The algorithm below starts with a random permutation matrix as the initial
solution. Note that the procedure is equipped with a $isValid()$ function that
tests if a solution satisfies constraints.

Algorithm 5 Exterior Penalty Method

1:procedure solve($F$,$D$,$\alpha_{0}$,$\beta$)

2: $\alpha=\alpha_{0}$

3: solution $=$ randomPermutationMatrix$(n\times n)$

4: repeat

5: solution $=$ annealer.solve(convert($F,D,\alpha$), solution)

6: $\alpha:=\alpha\times\beta$

7: until isValid(solution)

8:end procedure

###  F.8 Overall procedure

The overall procedure could be abstractly outlined in Algorithm 6.

Algorithm 6 Solving QAP with decomposition heuristics

1:procedure run($F$,$D$,$n$)

2:$\triangleright$ n-#items/locations

3:$\triangleright$ k-#groups

4: $s=\sqrt{n}$ $\triangleright$ assume group size is $\sqrt{n}$

5: $Fs,Ds,FtoD=\textrm{decompose}(F,D,s)$ $\triangleright$ FtoD is the map
between subsets and intervals

6: $\alpha_{0}=10000$

7: $\beta=1.5$

8: $\textrm{subsolutions}=[]$

9: for $i$ in $[1..k]$ do

10: $\textrm{items}=F[i]$

11: $\textrm{locations}=D[FtoD[i]]$

12: $F^{\prime}=\textrm{emptyMatrix}(s\times s)$

13: $D^{\prime}=\textrm{emptyMatrix}(s\times s)$

14: for $p,q$ in $items$ do

15: $F^{\prime}[\textrm{localindex}(p)][\textrm{localindex}(q)]=F[p][q]$

16: end for

17: for $k,l$ in locations do

18: $D^{\prime}[\textrm{localindex}(k)][\textrm{localindex}(l)]=D[k][l]$

19: end for

20:
$\textrm{subsolutions.append}(\textrm{solve}(F^{\prime},D^{\prime},\alpha_{0},\beta))$

21: end for

22: combineSubsolutions(subsolutions)

23:end procedure

Note that $\alpha_{0}$ and $\beta$ are hyper-parameters for penalty weight.
Instead of statically determined, we iteratively adjust the penalty weight for
the warehouse assignment problem to find the best solutions.

###  F.9 Discussion on improving the solutions

An important assumption in the comparison of the quality of solutions is that
QAP is a good way of modelling warehouse assignments. With that assumption in
mind, We observed that OOS gives a shorter picking distance than DA. This is
to be expected since OOS is manually tuned and optimised for WH-270b, whereas
QAP and its heuristic on DA are completely generic, modulo the initial
estimation of penalty weights. A better solution to the QAP on DA will then
yield a better solution for the warehouse assignment problem.

_Model for warehouse assignment_ The current QAP model for warehouse
assignment is the primitive one from Mantel et al. (2007). There are many
holes to plug in its assumptions. In particular, the standard QAPs need to be
solved with standard tools that give reliable solutions, which in turn can be
simulated a-posteriori to ascertain their effectiveness in modelling the
warehouse assignment problem. The problem of modelling warehouse assignment is
also more subtle. There are other measures of interaction frequency other than
the basic one, each of different numerical ranges and effectiveness, as
reviewed in Kofler (2014). This is another aspect that carries the potential
for improvement.

## References

  * Vert et al. [2019] Daniel Vert, Renaud Sirdey, and Stephane Louise.  On the limitations of the chimera graph topology in using analog quantum computers.  In _Proceedings of the 16th ACM international conference on computing frontiers_ , pages 226–229, 2019.
  * [2] D-Wave Systems Inc.  D-wave qpu architecture: Topologies.  https://docs.dwavesys.com/docs/latest/c_gs_4.html.  Accessed: 2021-09-01.
  * Tsukamoto et al. [2017] Sanroku Tsukamoto, Motomu Takatsu, Satoshi Matsubara, and Hirotaka Tamura.  An accelerator architecture for combinatorial optimization problems.  _Fujitsu Sci. Tech. J_ , 53(5):8–13, 2017.
  * Zhu et al. [2020] Zheng Zhu, Chao Fang, and Helmut G. Katzgraber.  borealis —a generalized global update algorithm for boolean optimization problems.  _Optimization Letters_ , 14(8):2495–2514, 2020\.
  * Aramon et al. [2019] Maliheh Aramon, Gili Rosenberg, Elisabetta Valiante, Toshiyuki Miyazawa, Hirotaka Tamura, and Helmut G Katzgraber.  Physics-inspired optimization for quadratic unconstrained problems using a digital annealer.  _Frontiers in Physics_ , 7:48, 2019.
  * Gonzalez [2021] Constantin Gonzalez.  Cloud based qc with amazon braket.  _Digitale Welt_ , 5(2):14–17, 2021.
  * Born and Fock [1928] Max Born and Vladimir Fock.  Beweis des adiabatensatzes.  _Zeitschrift für Physik_ , 51(3):165–180, 1928.
  * Tong et al. [2005] DM Tong, K Singh, Leong Chuan Kwek, and Choo Hiap Oh.  Quantitative conditions do not guarantee the validity of the adiabatic approximation.  _Physical review letters_ , 95(11):110407, 2005\.
  * McGeoch [2014] Catherine C McGeoch.  Adiabatic quantum computation and quantum annealing: Theory and practice.  _Synthesis Lectures on Quantum Computing_ , 5(2):1–93, 2014.
  * Hauke et al. [2020] Philipp Hauke, Helmut G Katzgraber, Wolfgang Lechner, Hidetoshi Nishimori, and William D Oliver.  Perspectives of quantum annealing: Methods and implementations.  _Reports on Progress in Physics_ , 83(5):054401, 2020.
  * Wang et al. [2019] Luzhi Wang, Shuli Hu, Mingyang Li, and Junping Zhou.  An exact algorithm for minimum vertex cover problem.  _Mathematics_ , 7(7):603, 2019.
  * Cai et al. [2013] Shaowei Cai, Kaile Su, Chuan Luo, and Abdul Sattar.  Numvc: An efficient local search algorithm for minimum vertex cover.  _Journal of Artificial Intelligence Research_ , 46:687–716, 2013.
  * Li et al. [2020] Ruizhi Li, Shuli Hu, Shaowei Cai, Jian Gao, Yiyuan Wang, and Minghao Yin.  Numwvc: A novel local search for minimum weighted vertex cover problem.  _Journal of the Operational Research Society_ , 71(9):1498–1509, 2020.
  * Petersen [1997] Charles G. Petersen.  An evaluation of order picking routeing policies.  _International Journal of Operations & Production Management_, 17(11):1098–1111, 1997.
  * Fumi et al. [2013] Andrea Fumi, Laura Scarabotti, and Massimiliano M Schiraldi.  Minimizing warehouse space with a dedicated storage policy.  _International Journal of Engineering Business Management_ , 5:21, 2013.
  * Malmborg and Bhaskaran [1990] Charles J Malmborg and Krishnakumar Bhaskaran.  A revised proof of optimality for the cube-per-order index rule for stored item location.  _Applied Mathematical Modelling_ , 14(2):87–95, 1990.
  * Mantel et al. [2007] Ronald J Mantel, Peter C Schuur, and Sunderesh S Heragu.  Order oriented slotting: a new assignment strategy for warehouses.  _European Journal of Industrial Engineering_ , 1(3):301–316, 2007.
  * Petersen et al. [2004] Charles G. Petersen, Gerald R. Aase, and Daniel R. Heiser.  Improving order-picking performance through the implementation of class-based storage.  _International Journal of Physical Distribution & Logistics Management_, 34(7):534–544, 2004.
  * Ruijter [2007] Herwen Ruijter.  Improved storage in a book warehouse: Design of an efficient tool for slotting the manual picking area at wolters-noordhoff.  Master’s thesis, University of Twente, 2007.
  * Raidl et al. [2020] Günther Raidl et al.  pymhlib.  https://github.com/ac-tuwien/pymhlib, 2020.
  * Gasquez [2015] David Gasquez.  qap.  https://github.com/davidgasquez/qap, 2015.
  * Tsige [2013] Mersha T Tsige.  Improving order-picking efficiency via storage assignments strategies.  Master’s thesis, University of Twente, The Netherlands, February 2013\.
  * Kofler [2014] Monika Kofler.  _Optimising the storage location assignment problem under dynamic conditions/eingereicht von: Monika Kofler_.  PhD thesis, Linz, 2014.