# Unbalanced penalization: A new approach to encode inequality constraints of
combinatorial problems for quantum optimization algorithms

J. A. Montañez-Barrera  Corresponding author: J. A. Montañez-Barrera;
<EMAIL_ADDRESS>Jülich Supercomputing Centre, Institute for
Advanced Simulation, Forschungszentrum Jülich, 52425 Jülich, Germany  Dennis
Willsch  Jülich Supercomputing Centre, Institute for Advanced Simulation,
Forschungszentrum Jülich, 52425 Jülich, Germany  A. Maldonado-Romo  Centro de
investigación en Computación, Instituto Politécnico Nacional, México
Entropica Labs, 186b Telok Ayer Street, Singapore 068632  Kristel Michielsen
Jülich Supercomputing Centre, Institute for Advanced Simulation,
Forschungszentrum Jülich, 52425 Jülich, Germany  AIDAS, 52425 Jülich, Germany
RWTH Aachen University, 52056 Aachen, Germany

###### Abstract

Solving combinatorial optimization problems of the kind that can be codified
by quadratic unconstrained binary optimization (QUBO) is a promising
application of quantum computation. Some problems of this class suitable for
practical applications such as the traveling salesman problem (TSP), the bin
packing problem (BPP), or the knapsack problem (KP) have inequality
constraints that require a particular cost function encoding. The common
approach is the use of slack variables to represent the inequality constraints
in the cost function. However, the use of slack variables considerably
increases the number of qubits and operations required to solve these problems
using quantum devices. In this work, we present an alternative method that
does not require extra slack variables and consists of using an unbalanced
penalization function to represent the inequality constraints in the QUBO.
This function is characterized by larger penalization when the inequality
constraint is not achieved than when it is. We evaluate our approach on the
TSP, BPP, and KP, successfully encoding the optimal solution of the original
optimization problem near the ground state cost Hamiltonian. Additionally, we
employ D-Wave Advantage and D-Wave hybrid solvers to solve the BPP, surpassing
the performance of the slack variables approach by achieving solutions for up
to 29 items, whereas the slack variables approach only handles up to 11 items.
This new approach can be used to solve combinatorial problems with inequality
constraints with a reduced number of resources compared to the slack variables
approach using quantum annealing or variational quantum algorithms.

Keywords

QUBO; inequality constraints; knapsack; bin packing; travel salesman problem;
quantum optimization; VQA; QAOA; quantum annealing; D-Wave Advantage; D-Wave
Hybrid; combinatorial optimization.

##  I Introduction

The exploration of quantum enhancement in combinatorial optimization problems
has been extended over several years. There are three main reasons for this;
first, many combinatorial optimization problems can be encoded in
Hamiltonians, the ground state being the optimal solution of the problem [1,
2]; second, they are commonly hard to solve and have practical applications
[3]; and third, the quantum algorithms to solve these problems need few
resources and can be tested on current state-of-the-art quantum hardware [4,
5].

The usual approach for encoding combinatorial optimization problems on quantum
processing units (QPUs) is to transform them in their quadratic unconstrained
binary optimization (QUBO) representation and obtain the cost Hamiltonian
after a change of variables. In the QUBO encoding, the constraints of the
combinatorial optimization problems are added to the cost function as
penalization terms along with the objective function.

However, we are currently in a noisy intermediate-scale quantum (NISQ) stage
and QPUs are restricted in terms of qubits number, noise level, and circuit
depth [6]. The main consequence is that numerous problems solvable by quantum
algorithms such as Shor’s algorithm [7, 8] remain unpractical to test on real
quantum hardware. Despite this fact, solving combinatorial optimization
problems is possible on current hardware. For example, on gate-based QPUs,
Variational Quantum Algorithms (VQA) [9] are applicable even with the above
limitations thanks to their short depth and resistance to noise [10]. Besides,
quantum annealers, using the quantum annealing (QA) principle [11, 12, 13],
are another quantum technology used to solve combinatorial optimization
problems [14, 15, 16]. In this respect, QA shows some advantages under certain
conditions compared to classical annealing [17, 18, 19].

Of the aforementioned VQA algorithms, the most studied for combinatorial
optimization problems is the Quantum Approximate Optimization algorithm (QAOA)
[20]. Even if it is still unclear whether QAOA will give any advantage
compared with other algorithms [21], its simplicity makes it of great interest
for analytical and practical purposes [4]. In the simplest QAOA version, the
cost Hamiltonian of a combinatorial optimization problem is encoded in a
parametric unitary gate along with a “mixer”, a second parametric unitary gate
that does not commute with the first unitary gate. In this context, the
parameters are adjusted to minimize the expectation value of the cost
Hamiltonian using a classical optimizer. Multiple approaches for solving
combinatorial optimization problems using QAOA or QA can be found in the
literature, for example, in logistics [22], finance [23, 5, 24, 25], energy
[26], communications [27], automotive industry [28], traffic signaling [29],
among others.

From the different sets of problems that can be solved using quantum hardware,
the ones with inequality constraints require an extra number of variables to
get their QUBO representation. For example, the BPP and the TSP have many
inequality constraints that increase with the number of items and cities,
respectively. The usual approach implemented in software development kits
(SDKs) such as Qiskit [30] or D-Wave Ocean [31] is to use slack variables to
encode the inequality constraints. Still, such an approach increments the
number of qubits and connections needed to solve these problems.

In this paper, we propose an alternative method to encode inequality
constraints in the QUBO formulation of combinatorial optimization problems. In
this new heuristic encoding method, inequality constraints are encoded using
an unbalanced penalization formula. This formulation adds larger penalization
when the inequality constraint is not fulfilled than when it is. We test our
method on the TSP, the BPP, and the KP using OpenQAOA [32], a python-based
library developed by Entropica Labs, and the Jülich universal quantum computer
simulator (JUQCS) [33, 34, 35] for problems with up to 43 qubits. In the
majority of the cases, the optimal solution of the original combinatorial
optimization problem is encoded in the ground state of the cost Hamiltonian of
their respective QUBO. However, in some cases, the optimal solution is in the
vicinity of the ground state. For such cases, our method still provides a good
approximation given the whole set of eigenvalues.

Additionally, we test the performance for finding solutions for the BPP using
the D-Wave Hybrid solver and the quantum annealer D-Wave Advantage 5.3 System
JUPSI, a quantum computing system developed by D-Wave Systems Inc that has
more than 5000 qubits located in Jülich, Germany. Our findings demonstrate
that our method outperforms existing approach in terms of both the quantity
and quality of solutions obtained for the BPP. Specifically, using the
unbalanced penalization approach, we achieve remarkable results on both
systems. With the D-Wave Advantage system, we obtain solutions for up to 7
nodes, while with the D-Wave Hybrid solver, we achieve solutions for up to 29
nodes. In contrast, the traditional slack variables approach only handles
solutions for up to 7 and 11 nodes for the D-Wave Advantage and D-Wave Hybrid,
respectively.

The rest of the paper is organized as follows. Sec. II provides a description
of the implementation of the combinatorial optimization problems using QUBOs
and an overview of the TSP, BPP, and KP. Section III presents a description of
the unbalanced penalization approach and a new metric called coefficient of
performance (CoP) to study the efficiency of QAOA to solve combinatorial
optimization problems. Section IV presents the results and discussion.
Finally, Section V provides some conclusions. The source code for the new
method setup and implementation on the TSP, BPP, and KP can be found at
https://jugit.fz-juelich.de/qip/unbalanced-penalizations-qubo.

##  II Quadratic unconstrained binary optimization

The set of combinatorial problems that can be represented by the QUBO
formulation are characterized by functions of the form

$f(\mathrm{x})=\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}q_{ij}x_{i}x_{j},$ (1)

where $n$ is the number of variables, $q_{ij}\in\mathbb{R}$ are coefficients
associated to the specific problem, and $x_{i}\in\\{0,1\\}$ are the binary
variables of the problem. Note that $x_{i}x_{i}\equiv x_{i}$ and
$q_{ij}=q_{ji}$ in this formulation. Therefore, the general form of a
combinatorial optimization problem solvable by QPUs is given by the cost
function

$f(\mathrm{x})=\sum_{i=1}^{n-1}\sum_{j>i}^{n}q_{ij}x_{i}x_{j}+\sum_{i=1}^{n}q_{ii}x_{i},$
(2)

and equality constraints given by

$\sum_{i=1}^{n}c_{i}x_{i}=C,\ c_{i}\in\mathbb{Z},$ (3)

and inequality constraints given by

$\sum_{i=1}^{n}l_{i}x_{i}\leq B,\ l_{i}\in\mathbb{Z}.$ (4)

To transform these problems into the QUBO formulation the constraints are
added as penalization terms. In this respect, the equality constraints are
included in the cost function using the following penalization term

$\lambda_{0}\left(\sum_{i=1}^{n}c_{i}x_{i}-C\right)^{2},$ (5)

where $\lambda_{0}$ is a penalization coefficient that should be chosen to
guarantee that the equality constraint is fulfilled and $C$ is a constant
value. In the case of inequality constraint, the common approach is to use a
slack variable [36, 37]. The slack variable, $S$, is an auxiliary variable
that makes a penalization term vanish when the inequality constraint is
achieved,

$B-\sum_{i=1}^{n}l_{i}x_{i}-S=0.$ (6)

Therefore, when Eq.(4) is satisfied, Eq.(6) is already zero. This means the
slack variable, $S$, must be in the range $0\leq S\leq
B-\min_{x}\sum_{i=1}^{n}l_{i}x_{i}$. To represent the $slack$ variable in
binary form, the slack is decomposed in
$N=\lfloor\log_{2}(\max_{x}B-\sum_{i=1}^{n}l_{i}x_{i})+1\rfloor$ binary
variables:

$S=\sum_{k=0}^{N-1}2^{k}s_{k},$ (7)

where $s_{k}$ are the slack binary variables. Then, the inequality constraints
are added as penalization terms by

$\lambda_{1}\left(B-\sum_{i=1}^{n}l_{i}x_{i}-\sum_{k=0}^{N-1}2^{k}s_{k}\right)^{2}.$
(8)

Combining Eq.(2) and the two kinds of constraints Eq.(5) and Eq.(8), the
general QUBO representation of a given combinatorial optimization problem is
given by

$\min_{x}\left(\sum_{i=1}^{n-1}\sum_{j>i}^{n}c_{ij}x_{i}x_{j}+\sum_{i=1}^{n}h_{i}x_{i}+\lambda_{0}\left(\sum_{i=1}^{n}q_{i}x_{i}-C\right)^{2}+\lambda_{1}\left(B-\sum_{i=1}^{n}l_{i}x_{i}-\sum_{k=0}^{N-1}2^{k}s_{k}\right)^{2}\right).$
(9)

Following the same principle, more constraints can be added and note that
after some manipulations, Eq.(9) can be rewritten in the form of Eq.(2). The
last step to represent the QUBO problem on QPUs is to change the $x_{i}$
variables to spin variables $z_{i}\in\\{1,-1\\}$ by the transformation
$x_{i}=(1-z_{i})/2$. Hence, Eq.(2) represented in terms of the cost
Hamiltonian model reads

$H_{c}(\mathrm{z})=\sum_{i=1}^{n-1}\sum_{j>i}^{n}q_{ij}(1-z_{i})(1-z_{j})/4+\sum_{i=1}^{n}q_{ii}(1-z_{i})/2.$
(10)

###  II.1 Traveling salesman problem

In the TSP, a traveler has to stop by a set of cities, finding the shortest
route to visit every city once and return to the starting point. Different
formulations of the TSP exist, but in this work, we focus on the Dantzig-
Fulkerson-Johnson (DFJ) formulation [38], which is the linear programming (LP)
version of the TSP. Note that while there are other formulations of this
problem that do not require inequality constraints, e.g., [1], the DFJ
formulation is used here to test the capabilities of the unbalanced
penalization method to encode multiple inequality constraints. In this
formulation, the problem is given by

$\min\sum_{i=1}^{n}\sum_{j\neq i,j=1}^{n}c_{ij}x_{ij},$ (11)

subject to the set of constraints,

$\sum_{i=1,i\neq j}^{n}x_{ij}=1\qquad\forall j=1,...,n,$ (12)

$\sum_{j=1,j\neq i}^{n}x_{ij}=1\qquad\forall i=1,...,n,$ (13)

$\sum_{i\in Q}\sum_{j\neq i,j\in Q}x_{ij}\leq|Q|-1\ \ \forall
Q\subsetneq\\{1,...,n\\},|Q|\geq 2,$ (14)

$x_{ij}\in\\{0,1\\}\qquad\forall i=1,..,n\qquad\forall j=1,..,n,$ (15)

where $n$ is the number of cities (nodes), $c_{ij}$ is the distance between
the cities $i$ and $j$, $x_{ij}$ is a binary variable that represents if the
path from the city $i$ to the city $j$ is used or not, and $Q$ represents a
sub-tour on a set of cities. From the above equations, Eq.(11) is the cost
function for the distance traveled, Eq.(12) and Eq.(13) restrict the traveler
to take only one path to enter and one to leave a city, and Eq.(14) are the
inequality constraints that prohibit sub-tours in the solution. Of all the
inequality constraints analyzed in this work, constraint Eq.(14) is the most
expensive in terms of the number of slack variables. Fig. 1 (a) shows the
number of qubits needed to represent the TSP for a different number of cities.

For this work, we randomly place $n$ cities on a 50x50 grid with $c_{ij}$
being the Euclidian distance between cities. We vary from 2 to 7 cities and
for problems up to 5 cities, we generate 10 random problems to test the
generalization of our method.

Figure 1:  Number of qubits needed to solve (a) the TSP for different number
of cities, (b) the BPP for different number of items, and (c) the KP for
different number of items. The solid line represents the number of variables
of the problem and the dashed line represent the variables of the problem plus
the slack variables needed to represent the inequality constraints of the
problem for a different number of nodes.

###  II.2 Bin packing problem

In the BPP, a collection of items must be stored in the minimum possible
number of bins. In this case, the items have a weight associated and the bins
are restricted to carry a maximum weight. This problem has many real-world
applications such as loading trucks with a weight restriction [39], container
scheduling [40], or design of FPGA chips [41]. In terms of complexity, the BPP
is an NP-hard problem and its formulation is given by

$\min\sum_{j=1}^{m}y_{j},$ (16)

subject to:

$\sum_{i=1}^{n}w_{i}x_{ij}\leq By_{j}\qquad\forall j=1,...,m,$ (17)

$\sum_{j=1}^{m}x_{ij}=1\qquad\forall i=1,...,n,$ (18)

$x_{ij}\in\\{0,1\\}\qquad\forall i=1,..,n\qquad\forall j=1,..,m,$ (19)
$y_{j}\in\\{0,1\\}\qquad\forall j=1,..,m,$ (20)

where $n$ is the number of items (nodes), $m$ is the number of bins, $w_{i}$
is the i-th item weight, $B$ is the maximum weight of each bin, $x_{ij}$ and
$y_{j}$ are binary variables that represent if the item $i$ is in the bin $j$,
and whether bin $j$ is used or not, respectively. From the above equations,
Eq.(16) is the cost function to minimize the number of bins, Eq.(17) is the
inequality constraint for the maximum weight of a bin, Eq.(18) is the equality
constraint to restrict that an item is only in one of the bins, and Eqs.(19)
and (20) means that $y_{i}$ and $x_{ij}$ are binary variables.

Fig. 1 (b) shows the number of variables needed to solve the BPP when the
number of bins is equal to the number of items. The weights are chosen
randomly with values between 4 and 20 and the maximum bin weight equal to 20.
Without losing generality, two further simplifications are made, the first
assigns the first item to the first bin $x_{11}=1$ and therefore $x_{1j}=0\
\forall\ j\in\\{2,...,m\\}$. Second replaces $y_{j}=1\ \forall\ j\in\
\\{1,2,...,N^{bin}_{min}\\}$ with
$N^{bin}_{min}=\lceil(\sum_{i=1}^{n}s(i)\rceil)/B$ because the minimum of bins
required is known.

###  II.3 Knapsack problem

In the KP, a set of items with associated weights and values should be stored
in a knapsack. The problem is to maximize the value of the items (nodes)
transported in the knapsack. The KP is restricted by the maximum weight the
knapsack can carry. The KP is the simplest nontrivial integer programming
model with binary variables, only one constraint, and positive coefficients.
It is formally defined by

$\max\sum_{i=1}^{n}p_{i}x_{i},$ (21)

$\sum_{i=1}^{n}w_{i}x_{i}\leq W,$ (22)

where $n$ is the number of items, $p_{i}$ and $w_{i}$ is the value and weight
of the $ith$ item, respectively, $x_{i}$ is the binary variable that
represents whether the $ith$ item is in the knapsack or not, and W is the
maximum weight that the knapsack can transport. Fig. 1 (c) shows the number of
qubits needed to solve the KP for a different number of items. The set of
problems is created with weights and values selected randomly ranging between
1 and 63 for the values, between 1 and 127 for the weights, and knapsack
maximum weight is equal to 70% of the sum of the weights of all the items.

An estimate of the number of slack variables needed to represent the QUBO of
the TSP, the BPP, and the KP is shown in Fig. 2. In the TSP, the number of
slack variables increases exponentially with the number of cities added. In
the BPP, it requires many slack binary variables increasing proportionally to
the number of bins of the problem. Lastly, in the KP, the number of slack
variables needed is constant and depends on the maximum weight allowed in the
knapsack.

Figure 2:  The number of slack variables needed to solve the TSP (diamonds),
the BPP (circles), and the KP (triangles) based on the number of nodes. The
solid lines are guides to the eye.

##  III Unbalanced penalization

The unbalanced penalization approach we propose is an approximation for
including inequality constraints of combinatorial optimization problems in the
QUBO representation of the cost function. This method does not require
additional slack variables, and the inequality constraint of Eq. (4) given by

$h(\mathrm{x})=B-\sum_{i}l_{i}x_{i}\geq 0,$ (23)

is replaced by a penalization term in the QUBO formulation that makes larger
penalization for negative terms and smaller for positive. The idea of the
unbalanced penalization function comes from the shape of the exponential decay
curve, $f(x)=e^{-h(\mathrm{x})}$. In this function, positive values of
$h(\mathrm{x})$, given by Eq. (23), make $f(\mathrm{x})\approx 0$ while
negative values make it grow exponentially, as shown in Fig. 3. However, the
exponential function cannot be encoded as a QUBO penalization term. Therefore,
we consider an expansion of the exponential function to quadratic order. The
penalization term is given by

$e^{-h(\mathrm{x})}\approx 1-h(\mathrm{x})+\frac{1}{2}h(\mathrm{x})^{2}.$ (24)

In general, we modify Eq. (24) to include free parameters to be adjusted for
the different kinds of problems. Therefore, Eq. (24) is rewritten by

$\displaystyle\zeta(x)$ $\displaystyle=-\lambda_{1}h(x)+\lambda_{2}h(x)^{2}$
(25)
$\displaystyle=-\lambda_{1}\left(B-\sum_{i}l_{i}x_{i}\right)+\lambda_{2}\left(B-\sum_{i}l_{i}x_{i}\right)^{2},$
(26)

where $\lambda_{1,2}$ are a set of multipliers that can be tuned for the
specific problem and the constant term is removed because the position of the
cost function’s minimum is independent of the constant term. Note that this
approach can be extended to $g(\mathrm{x})=B-\sum_{i}l_{i}x_{i}\leq 0$, in
that case $f(x)=e^{g(x)}$. The new cost function based on Eq. (9) and the
unbalanced penalization approach is given by

$\min_{x}\left(\sum_{i=1}^{n}\sum_{j\neq
i}^{n}c_{ij}x_{i}x_{j}+\sum_{i=1}^{n}h_{i}x_{i}+\lambda_{0}\left(\sum_{i=1}^{n}q_{i}x_{i}-C\right)^{2}-\lambda_{1}\left(B-\sum_{i}l_{i}x_{i}\right)+\lambda_{2}\left(B-\sum_{i}l_{i}x_{i}\right)^{2}\right),$
(27)

Figure 3:  Comparison between $e^{-h(\mathrm{x})}$ and the unbalanced function
$1-h(\mathrm{x})+\frac{1}{2}h(\mathrm{x})^{2}$.

To tune the $\lambda_{0,1,2}$ parameters of Eq. (27) for the TSP, the BPP, and
the KP, we use the Nelder-Mead optimization method to explore the set of
$\lambda_{0,1,2}$ values that brings the optimal solution of the original
combinatorial optimization problems as close as possible to the ground state
of the cost Hamiltonian, Eq. (10), of their QUBO representation. For the TSP
we use a random problem with 4 cities, for the BPP we use a random problem
with 4 items, and for the KP we use a problem with 13 items. The set of values
$\lambda_{0,1,2}$ found by the Nelder-Mead optimization method for each
problem are presented in Table 1. After we find the best set of values for the
specific problem, we test if those values generalize to other random cases. It
is a crucial aspect of the unbalanced penalization method that the set of
$\lambda_{0,1,2}$ of a random case generalizes for different and/or larger
problems. The results on these aspects of the unbalanced penalization method
are presented in Sec.IV.

Table 1: Parameters $\lambda_{0,1,2}$ for the TSP, BPP, and KP used to translate the combinatorial optimization problems into the QUBO representation using the unbalanced penalization approach (Eq. (27).) For all the equality constraints of each problem, we use the same $\lambda_{0}$ and for the inequality constraints the same $\lambda_{1,2}$. |  $\lambda_{0}$ |  $\lambda_{1}$ |  $\lambda_{2}$
---|---|---|---
TSP |  38.2584 |  18.2838 |  57.0375
BPP |  20.5198 |  7.2949 |  0.8583
KP |  - |  0.9603 |  0.0371

### Quantum Annealing

Quantum annealing is a computational technique designed to tackle optimization
problems by leveraging quantum mechanical effects in the search for low-energy
states of Ising Hamiltonians [42, 43]. By employing this method, the system
can explore a broader solution space and potentially find good quality
solution more efficiently than classical optimization approaches.

In order to validate the effectiveness of our proposed method on real
hardware, we conducted experiments using the D-Wave Advantage and the D-Wave
Hybrid Quantum-Classical solver to obtain results for the BPP. To ensure a
fair comparison, we selected random BPP problems with the same settings as
described in Section II.2, without any simplifications.

For the experiments using the D-Wave Advantage, we conducted 20 trials on
random problem instances ranging from 3 to 8 nodes (items). Each trial
consisted of 5000 samples. For the hybrid solver, we performed 20 experiments
with a single run, allowing a maximum execution time of 3 seconds. The problem
instances for the hybrid solver ranged from 5 to 31 nodes (items).

In both cases, we utilized the $\lambda_{0,1,2}$ parameters as presented in
Table 1. It is important to note that our objective was not to fine-tune the
hyperparameters for either solver, but rather to demonstrate that under
similar conditions, unbalanced penalization yields significantly better
solutions in terms of both quantity and quality. Results are presented in Sec.
IV.

### Quantum Approximate Optimization Algorithm

To test the difference between the slack variables approach vs. the unbalanced
penalization encoding on different random problems of the TSP, the BPP, and
the KP, we use the QAOA with one layer. It is important to acknowledge that,
while the impact of penalization terms of the constraints remains uncertain in
terms of their influence on the performance of the Quantum Approximate
Optimization Algorithm (QAOA)[44], we employ QAOA to demonstrate the
contrasting probability distributions resulting from the two approaches. In
this case, the cost Hamiltonian, $H_{c}$, obtained from the QUBO formulation,
is translated into a parametric unitary gate given by

$U(H_{c},\gamma)=e^{-i\gamma H_{c}},$ (28)

where $\gamma$ is a parameter to be optimized. A second unitary operator
applied is

$U(B,\beta)=e^{i\beta X},$ (29)

where $\beta$ is the second parameter that must be optimized and
$X=\sum_{i=1}^{n}\sigma_{i}^{x}$ with $\sigma_{i}^{x}$ the Pauli-x quantum
gate applied to qubit $i$. The general QAOA circuit is shown in Fig. 4. Here,
$R_{X}(\theta)=e^{-i\frac{\theta}{2}\sigma_{x}}$, $p$ represents the number of
repetitions of the unitary gates Eqs. 28 and 29 with each repetition having
separate values for $\gamma_{p}$ and $\beta_{p}$, and the initial state is a
superposition state $|+\rangle^{\otimes n}$. The case of QAOA with $p=1$ is
suitable for visualization of the landscape and in Sec. IV, we show a
comparison between both encodings and even though they are similar, the
addition of extra qubits for the slack variables approach has a great impact
on the probability of finding the optimal solution of the original
combinatorial optimization problem.

Figure 4: Schematic representation of QAOA for $p$ layers. The parameters
$\gamma$ and $\beta$ for each layer are the ones to be optimized. QAOA is used
to compare the slack variables method vs. the unbalanced penalization method
for the encoding of the inequality constraints.

### Coefficient of performance

Additionally, we introduce the coefficient of performance (CoP), a quantifier
that expresses how good a quantum algorithm is to give a specific solution
compared to random guessing over a set of solutions. In this perspective, we
use the CoP to determine how likely it is to find the optimal solution for the
TSP, the BPP, and the KP on the local minima of the QAOA algorithm for $p=1$.
The CoP is given by

$CoP=\frac{P(^{*}\mathrm{x})}{P_{R}},$ (30)

where $P(^{*}\mathrm{x})$ is the probability of getting a bit string
${}^{*}\mathrm{x}$ using a quantum algorithm, e.g, QAOA and $P_{R}=1/2^{n}$ is
the probability of a random guessing with n equal to the number of qubits
involved in the problem.

##  IV Results

### Optimal solution distribution

To determine the applicability of our method, we sorted the energy eigenvalues
of the cost Hamiltonian for the different combinatorial optimization problems.
The parameter of interest is the position of the energy eigenvalue that
describes the optimal solution of the original combinatorial optimization
problems within the list of all sorted energy eigenvalues of the cost
Hamiltonian of their QUBO representation. If the position of the optimal
solution is far from the ground energy of the cost Hamiltonian, the unbalanced
penalization encoding is not working as expected. Otherwise, if the energy
eigenvalue corresponding to the optimal solution is close to the ground state
energy given the total number of eigenvalues, the unbalanced encoding is
advantageous.

Figure 5 shows the energy eigenvalues of 10 random cases for the TSP with 5
nodes using the unbalanced penalization encoding approach. The eigenvalues are
sorted in terms of their energy from minimum to maximum with the x-axis
representing the eigenvalue position (number) of eigenvalues, and the y-axis
its energy. In the inset, showing the lowest 25 energy eigenvalues, the big
circles represent where the optimal result is positioned compared to all the
other eigenvalues (small circles) of the cost Hamiltonian. Here, the worst
eigenvalue position is 7 out of $2^{20}$ = 1048576 eigenvalues. This means
that in the worst scenario the optimal solution is among the lowest 0.00057%
of all eigenvalues.

Figure 5:  Eigenvalues distribution for the TSP with 5 cities (20 qubits) for
10 randomly generated problems using the unbalanced penalization method. The
inset shows the lowest 25 energy eigenvalues. The big circles are the optimal
solutions for the random problem and the small circles are the different
eigenvalues of the cost Hamiltonian of its QUBO representation. Note that each
eigenvalue is degenerate with multiplicity two, because the problem is
symmetric, e.g., one clockwise solution and another anti-clockwise.

From Fig. 5 it is seen that the unbalanced penalization approach does not
ensure that the optimal solution is the lowest eigenvalue of the cost
Hamiltonian, but in all the cases analyzed, it is very close to it. For
example, in Fig. 6, this feature is illustrated for one of the random problems
of Fig. 5. Here, there are a couple of eigenvalues that do not fulfill all the
restrictions of the problem but have lower energy than the optimal solution.
In the end, the unbalanced penalization encoding is a trade-off between adding
slack variables and ensuring the ground state of the cost Hamiltonian
corresponds to the optimal solution of the combinatorial optimization problem,
or reducing the number of variables but encoding the optimal solution in the
vicinity of the ground state. In many cases, there is no preference for
encoding the optimal solution in the lowest eigenvalue of the cost
Hamiltonian. For instance, QAOA brings the expectation value of a cost
Hamiltonian to a region of overall minimum energy, so we can expect that the
optimal and sub-optimal solutions occur with reasonable probability in regions
of minimal energy.

Figure 6:  The first 25 sorted eigenvalues of the Hamiltonian $H_{c}$
corresponding to the worst case shown in Fig. 5. The big circle corresponds to
the optimal solution of this problem and the small circles are the other
eigenvalues. The solution shown in the red inset is the ground state, which is
an invalid solution of the TSP. The solution in the green inset is the optimal
solution.

Fig. 7 shows the energy eigenvalues sorted for the cost Hamiltonian of the BPP
for 10 random problems with 5 items. In the inset, showing the 40 lowest
eigenvalues, the optimal solutions of the cost Hamiltonians are depicted as
big circles and the small circles are the other eigenvalues of the cost
Hamiltonian. In this case, there is an increasing number of degeneracies
explained by the symmetries of the problem (exchange items of one bin with the
others). To keep the same number of qubits (21) for all the problems, we
select 10 random cases that require more than 3 bins to store the items. Here,
the worst case is in the position 36 out of $2^{21}$=2097152 eigenvalues which
means it is within $0.00124\%$, of all eigenvalues.

Figure 7:  Eigenvalue distribution for the BPP with 5 items (21 qubits) for 10
randomly generated problems using the unbalanced inequality-constrained
penalization method. The inset shows the 40 lowest energy eigenvalues. The big
circles are the optimal solutions for the random problem and the small circles
are the different eigenvalues for the cost Hamiltonian of its QUBO
representation. Note that each eigenvalue has a different degeneracy, because
the problem has some symmetries.

Finally, Fig. 8 shows the eigenvalue distribution for 10 random cases of the
KP with 21 items (21 qubits). In this case, there are no degeneracies, but the
possible solutions are close to each other. The inset shows the first 50
eigenvalues. In this case, for one of the random cases the optimal solution is
located at position 49 out of $2^{21}$=2097152 eigenvalues, which means it is
within $0.0023\%$ of the first eigenvalues. In this case, the parameters
$\lambda_{0,1,2}$ were optimized for a random case with 13 items, which
suggests great generalization capabilities of the current method.

Figure 8:  Eigenvalue distribution for the KP with 21 items (21 qubits) for 10
randomly generated problems using the unbalanced penalization method. The
inset shows the 50 lowest eigenvalues. The big circles give the position of
the optimal solutions for the random problems and the small circles depict the
position of the different eigenvalues of the cost Hamiltonian. In this case,
there are no degeneracies but the energies are close to each other.

### Quantum Annealing

Figure 9 illustrates the results obtained using the D-Wave Advantage. The
upper x-axis in this figure, as well as in subsequent plots, represents the
logical variables required to represent each problem size, as per the
formulation provided in Section II.2. Across all cases, the unbalanced
penalization approach consistently generates more optimal solutions compared
to the slack variables approach. Although the optimal solution cannot be
guaranteed to be the ground state of the Hamiltonian, more optimal solutions
are found for the unbalanced penalization method.

For the scenario involving 8 nodes, neither of the two methods were able to
find valid solutions. However, it is important to acknowledge that the number
of statistical samples used in these experiments might not have been
sufficient to draw conclusions about the limit beyond which no further
solutions are found on the quantum annealer.

Figure 9:  D-Wave Advantage success rate of finding optimal solutions for
different cases of the BPP ranging from 3 to 8 nodes (items). The 20 small
circles (squares) for each problem size consist of 5000 samples for the
unbalanced (slack) method.

Fig. 10 presents the results of the same experiments depicted in Fig. 9,
focusing on the quantity of valid solutions obtained. Consistent with the
observations in Fig. 9, the employment of the unbalanced penalization method
enhances the likelihood of discovering valid solutions compared to the slack
approach. However, it is important to note that in the case of 7 nodes, where
the slack approach appears to exhibit superior performance, more samples of
valid solutions (than could be produced by the machine) would be needed to
draw a solid conclusion.

Figure 10:  D-Wave Advantage success rate of finding valid solutions for
different cases of the BPP ranging from 3 to 8 nodes (items).

In order to investigate the viability of finding valid solutions for larger
problem sizes, we employed the D-Wave Hybrid solver. We conducted experiments
using problem sizes ranging from 5 to 31 nodes (items), with intervals of 2,
and executed 20 random problem instances for each case. Each problem instance
was solved using a single run on the D-Wave Hybrid solver, with a maximum time
limit of 3 seconds.

Figure 11 depicts the probability of finding optimal solutions across the 20
random problems for each node case. Similar to the observations with the
D-Wave Advantage, the unbalanced penalization method consistently outperforms
the slack approach. It successfully discovers optimal solutions for nearly all
cases, with the exception of 19, 23, 27, and 31 nodes. Conversely, the slack
approach fails to find optimal solutions after 7 nodes.

Figure 11:  D-Wave Hybrid solver success rate of finding optimal solutions for
different cases of the BPP ranging from 5 to 31 nodes (items). Each point
consists of a single execution of 20 random problems with a maximum time for
the hybrid solver of 3 seconds.

Figure 12 illustrates the probability of finding valid solutions using both
the unbalanced penalization method and the slack variables approach on the
D-Wave Advantage. Consistent with the findings in Figure 11, the unbalanced
penalization method demonstrates superior performance compared to the slack
variables approach.

Up to 17 nodes, the unbalanced penalization method consistently discovers
valid solutions, whereas the slack variables approach only manages to find
valid solutions up to 11 nodes. Beyond 17 nodes, the probability of finding
valid solutions gradually declines for the unbalanced penalization method,
with no solutions being found for the 31-node case. This behavior can be
attributed to the fact that the maximum time limit remains constant throughout
the experiments.

Even though increasing the time will impact positively on improving the
performance of both methods, it is important to reiterate that the primary
objective of this section is not to provide high-quality solutions, but rather
to evaluate the performance of both methods under similar conditions.

Figure 12:  D-Wave Hybrid solver success rate of finding valid solutions for
different cases of the BPP ranging from 5 to 31 nodes (items). Each circle
(square) represents the mean value of 20 random problems with a maximum time
for the hybrid solver of 3 seconds for the unbalanced (slack) method.

### QAOA landscape

In this section, the results of the cost Hamiltonian energy landscape using
the QAOA algorithm with one layer are presented for the TSP, the BPP, and the
KP. The slack variables method and the unbalanced penalization method are used
for the cost Hamiltonian encoding. For visualization purposes, the energy
landscape for all cases is limited to the region of
$\gamma\in(-\pi/\max\\{q_{ij}\\}$, 0) where $q_{ij}$ is given in Eq. (10) and
$\beta\in(-\pi/2,0)$.

Fig. 13 shows the energy landscape of the QAOA algorithm for a TSP with 4
cities. This is equivalent to a 17 and 12 qubits problem for the slack
variables and the unbalanced penalization encodings, respectively. The
probabilities for finding the optimal solutions of the TSP are 0.02% and
0.227% for the slack variables and the unbalanced encodings, respectively.
This shows a clear advantage of the unbalanced encoding method that increases
the probability more than 10 times compared to the probability obtained with
the slack variables encoding. This characteristic is mainly attributed to the
difference in the number of qubits needed to represent the problem and thus
the size of the set of possible solutions. For the slack variables encoding
there are 131072 ($2^{17}$) possible solutions while for the unbalanced
penalization encoding there are 4096 ($2^{12}$) possible solutions.

Figure 13:  QAOA energy for $p=1$ for the cost Hamiltonian $\langle
H_{c}\rangle$ of the TSP with 4 cities with (a) slack variables encoding and
(b) unbalanced penalization encoding. The minimum energy is indicated with a
white dot. The probability that this minimum corresponds to the optimal
solution of the TSP problem is P(*x) = 0.02% and P(*x)=0.227% for encoding (a)
and (b), respectively.

Fig. 14 shows the energy landscape of the QAOA algorithm for a BPP with 3
items, which is equivalent to a 6 and 19 qubits problem for the unbalanced and
slack encodings, respectively. The probability of finding the optimal solution
of the BPP in the local minima of QAOA is 0.001% for the slack variables
encoding and 10.66% for the unbalance penalization method. Here, the optimal
solution is 10000 times more likely to be found using the unbalanced
penalization approach than using the slack variables approach.

Figure 14:  Same as the Fig. 13 for the BPP with 3 items. The probability that
this minimum corresponds to the optimal solution of the BPP problem is P(*x) =
0.001% and P(*x)=10.661% for encoding (a) and (b), respectively.

Finally, Fig. 15 shows the energy landscape of the QAOA algorithm for a KP
with 10 items, which is equivalent to a 18 and 10 qubits problem for the slack
variables and the unbalanced penalization encoding, respectively. Here, the
probability to find the optimal solution of the KP is 0.001% and 0.346% for
the slack variables and the unbalanced penalization encoding. This means
finding the solution with the unbalanced penalization method is 346 times more
likely than with the slack variables encoding.

Figure 15: Same as the Fig. 13 for the KP with 10 items. The probability that
this minimum corresponds to the optimal solution of the KP problem is P(*x) =
0.001% and P(*x)=0.346% for encoding (a) and (b), respectively.

### The coefficient of performance

Figure 16:  Coefficient of performance for (a) the TSP, (b) the BPP, and (c)
the KP for different problem sizes. The error bars represent the standard
deviation over 5 different random cases for problems with less than 22 qubits.
For problems with more than 21 qubits the data points represent 1 random case.

Figure 16 shows the CoP for different problem sizes for the TSP, BPP, and KP
using the unbalanced penalization method. All the cases show the probability
of finding the optimal solution for each problem using QAOA with 1 layer in
the minimum point of the energy landscape region $\beta\in(-\pi/2,0)$ and
$\gamma\in(-\pi/\max\\{q_{ij}\\},0)$. At this point, we do not use a classical
optimizer to find the optimal beta and gamma, but instead, we divide the beta
and gamma region in a 50x50 grid landscape and take the probability of the
optimal solution at the minimum energy. It is expected that using an
optimization method for this step improves the CoP while keeping the same
tendency. Fig. 16 (a) shows the CoP for the TSP problem with 2, 3, 4, 5, 6,
and 7 cities (2, 6, 12, 20, 30, and 42 qubits). For the cases with 30 and 42
qubits, we use JUQCS [33, 34] in its GPU-accelerated version [35] to obtain
the results. Note that the same TSP cases using the slack variables approach
will require 2, 7, 17, 36, 73, 148, and 304 qubits. We thus explore problem
sizes (5, 6, and 7 cities) with the unbalanced penalization approach that are
unfeasible using the slack approach on quantum computer simulators.
Interestingly, the CoP is increasing exponentially with the problem size,
which means that the probability of finding the optimal solution exponentially
increases when we compare it with random guessing in the set of possible
configurations (cf. [45]).

Fig. 16 (b) shows the CoP for the BPP for 2, 3, 4, 5, 6, and 7 items (2, 6,
13, 21, 32, and 43 qubits). For this case, the CoP is improving even more
rapidly compared to TSP and also here for the BPP with 6 and 7 items the
solutions using slack variables (61 and 78 qubits) is beyond range of what can
be simulated with quantum computer simulators. The largest simulations up to
43 qubits were performed on JUWELS Booster [46] as they required more than
$1/8\,\mathrm{PiB}$ of distributed memory and took more than one million core
hours. Finally, Fig. 16 (c) shows the CoP for the KP with 5 to 41 items (5 to
41 qubits) with increments of two (same number of qubits). In this case, the
CoP increases linearly with the number of qubits. This is a poor performance
if we compare it with the TSP and BPP. We suspect the improved CoP is related
to the number of constraints of the problem.

##  V Conclusions

We have presented unbalanced penalization, a new method to encode inequality
constraints of combinatorial optimization problems into QUBO penalizations.
The method does not require extra slack variables to encode the inequality
constraints. This is extremely beneficial because there is no increase in the
number of variables or qubits needed to encode the problems, especially for
the TSP and the BPP.

The method is suitable for VQA on gate-based quantum computers and QA on
quantum annealers, which both require expressing the combinatorial
optimization problems in terms of QUBOs. We have tested the method using QA
for the BPP and QAOA with 1 layer for the TSP, BPP, and KP. We have shown that
in 10 random cases for different problem sizes the optimal solution of the
combinatorial optimization problem is located in the vicinity of the ground
state energy of the cost Hamiltonian. The method is highly generalizable in
the sense that the tuned parameters $\lambda_{0,1,2}$ for small problems
generalize well to larger problems. For QA, we have demonstrated that the
unbalanced penalization method outperforms the slack variables approach,
finding better solutions in terms of quality and quantity. Additionally, the
probabilities to find the optimal solutions among the local minima of the
unbalanced penalization (Figs. 13 to 15) are much better than using slack
variables. Qualitatively, this can be understood because the addition of slack
variables increases the search space and therefore reduces the probability of
finding optimal and nearly-optimal solutions.

The proposed method was evaluated on both the D-Wave Advantage and D-Wave
Hybrid solvers. Our findings demonstrate that the method utilizing unbalanced
penalization consistently outperforms the slack variables approach in solving
the BPP. Specifically, the unbalanced penalization method successfully finds
solutions for problem instances with up to 31 items using the Hybrid solver,
while the slack variables approach achieves solutions for only up to 11 items.
We remark that additional evaluations of this method on D-Wave systems for the
TSP support the same conclusions [47].

We have tested the unbalanced penalization method close to the limits of what
is possible to be simulated on supercomputers using JUQCS. For the case of the
TSP, we went up to 42 qubits, for the BPP up to 43 qubits, and for the KP up
to 41 qubits. Indeed, we are exploring regions unfeasible using the slack
variables method. For example, for the TSP with 6 and 7 cities using the slack
variables method it requires 73 and 148 qubits, and for the BPP with 6 and 7
items it requires 61 and 78 qubits. In terms of probability, the unbalanced
penalization method shows large advantages for the small problems we made a
comparison for. For the TSP with 4 cities, the probability of finding the
optimal solution at the local minimum energy of the QAOA with $p=1$ using the
unbalanced penalization method is more than 10 times larger than the one using
the slack variables approach. For the BPP with 3 items, the probability for
finding the optimal solution with the unbalanced penalization method is more
than 10000 times larger, and for the KP with 10 items it is almost 346 higher.

Finally, we have presented a new metric, the CoP, which compares the
probability of obtaining a specific solution using a quantum algorithm against
random guessing over the set of possible solutions. The CoP is useful to
characterize how an algorithm performs with increasing problem size or to
compare how different combinatorial optimization problems perform for a
specific algorithm. We have shown the CoP for the TSP, BPP, and KP using QAOA
with $p=1$ in the region $\gamma\in(-\pi/\max\\{q_{ij}\\}$, 0) and
$\beta\in(-\pi/2$, 0). Interestingly, for all of them, the CoP increases with
the problem size. For the KP the increase is linear while for the TSP and BPP
the increase is exponential. Further research is worthwhile to see if this
provides an implicit advantage of using QAOA to solve combinatorial
optimization problems, what the characteristics of the TSP and the BPP are
that make the CoP increase exponentially, and to see how increasing the number
of QAOA layers improves the CoP.

###### Acknowledgements.

We would like to thank Leonardo Disilvestro and Vishal Sharma from Entropica
Labs for helpful comments and suggestions about the results of the present
work. J. A. Montañez-Barrera acknowledges support by the German Federal
Ministry of Education and Research (BMBF), funding program Quantum
technologies - from basic research to market, project QSolid (Grant No.
13N16149). D. Willsch acknowledges support from the project Jülich UNified
Infrastructure for Quantum computing (JUNIQ) that has received funding from
the German Federal Ministry of Education and Research (BMBF) and the Ministry
of Culture and Science of the State of North Rhine-Westphalia. The authors
gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-
centre.eu) for funding this project by providing computing time on the GCS
Supercomputer JUWELS at Jülich Supercomputing Centre (JSC).

## References

  * Lucas [2014] A. Lucas, Frontiers in Physics 2, 1 (2014), arXiv:1302.5843 .
  * Kochenberger  _et al._ [2014] G. Kochenberger, J. K. Hao, F. Glover, M. Lewis, Z. Lü, H. Wang, and Y. Wang, Journal of Combinatorial Optimization 28, 58 (2014).
  * Ohzeki [2020] M. Ohzeki, Scientific Reports 10, 1 (2020), arXiv:2002.05298 .
  * Harrigan  _et al._ [2021] M. P. Harrigan, K. J. Sung, M. Neeley, K. J. Satzinger, F. Arute, K. Arya, J. Atalaya, J. C. Bardin, R. Barends, S. Boixo, M. Broughton, B. B. Buckley, D. A. Buell, B. Burkett, N. Bushnell, Y. Chen, Z. Chen, Ben Chiaro, R. Collins, W. Courtney, S. Demura, A. Dunsworth, D. Eppens, A. Fowler, B. Foxen, C. Gidney, M. Giustina, R. Graff, S. Habegger, A. Ho, S. Hong, T. Huang, L. B. Ioffe, S. V. Isakov, E. Jeffrey, Z. Jiang, C. Jones, D. Kafri, K. Kechedzhi, J. Kelly, S. Kim, P. V. Klimov, A. N. Korotkov, F. Kostritsa, D. Landhuis, P. Laptev, M. Lindmark, M. Leib, O. Martin, J. M. Martinis, J. R. McClean, M. McEwen, A. Megrant, X. Mi, M. Mohseni, W. Mruczkiewicz, J. Mutus, O. Naaman, C. Neill, F. Neukart, M. Y. Niu, T. E. O’Brien, B. O’Gorman, E. Ostby, A. Petukhov, H. Putterman, C. Quintana, P. Roushan, N. C. Rubin, D. Sank, A. Skolik, V. Smelyanskiy, D. Strain, M. Streif, M. Szalay, A. Vainsencher, T. White, Z. J. Yao, P. Yeh, A. Zalcman, L. Zhou, H. Neven, D. Bacon, E. Lucero, E. Farhi, and R. Babbush, Nature Physics 17, 332 (2021), 2004.04197 .
  * Niroula  _et al._ [2022] P. Niroula, R. Shaydulin, R. Yalovetzky, P. Minssen, D. Herman, S. Hu, and M. Pistoia,  (2022), arXiv:2206.06290 .
  * Preskill [2018] J. Preskill, Quantum 2, 1 (2018), arXiv:1801.00862 .
  * Shor [1994] P. W. Shor, in  _Proceedings 35th Annual Symposium on Foundations of Computer Science_ (1994) pp. 124–134.
  * Shor [1997] P. W. Shor, SIAM J. Comput. 26, 1484 (1997).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. Benjamin, S. Endo, K. Fujii, J. Mcclean, K. Mitarai, X. Yuan, L. Cincio, and P. Coles, Nature Reviews Physics 3, 1 (2021).
  * Khairy  _et al._ [2020] S. Khairy, R. Shaydulin, L. Cincio, Y. Alexeev, and P. Balaprakash, Proceedings of the AAAI Conference on Artificial Intelligence 34, 2367 (2020).
  * Apolloni  _et al._ [1989] B. Apolloni, C. Carvalho, and D. de Falco, Stoch. Process. Their Appl. 33, 233 (1989).
  * Finnila  _et al._ [1994] A. B. Finnila, M. A. Gomez, C. Sebenik, C. Stenson, and J. D. Doll, Chem. Phys. Lett. 219, 343 (1994).
  * Kadowaki and Nishimori [1998a] T. Kadowaki and H. Nishimori, Phys. Rev. E 58, 5355 (1998a).
  * de Falco and Tamascelli [2011] D. de Falco and D. Tamascelli, RAIRO - Theoretical Informatics and Applications 45 (2011), 10.1051/ita/2011013.
  * Ayanzadeh  _et al._ [2021] R. Ayanzadeh, J. Dorband, M. Halem, and T. Finin, Scientific Reports 11, 1 (2021), arXiv:2010.00115 .
  * Willsch  _et al._ [2022a] D. Willsch, M. Willsch, C. D. Gonzalez Calaza, F. Jin, H. De Raedt, M. Svensson, and K. Michielsen, Quantum Inf. Process. 21, 141 (2022a).
  * Heim  _et al._ [2015] B. Heim, T. F. Rønnow, S. V. Isakov, and M. Troyer, Science 348, 215 (2015), arXiv:1411.5693 .
  * Yan and Sinitsyn [2022] B. Yan and N. A. Sinitsyn, Nature Communications 13, 1 (2022), arXiv:2110.12354 .
  * Tasseff  _et al._ [2022] B. Tasseff, T. Albash, Z. Morrell, M. Vuffray, A. Y. Lokhov, S. Misra, and C. Coffrin,  , 1 (2022), arXiv:2210.04291 .
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, (2014).
  * Willsch  _et al._ [2020] M. Willsch, D. Willsch, F. Jin, H. De Raedt, and K. Michielsen, Quantum Inf. Process. 19, 197 (2020).
  * Jiang  _et al._ [2022] H. Jiang, Z.-J. M. Shen, and J. Liu,  (2022), arXiv:2209.08246 .
  * Orús  _et al._ [2019] R. Orús, S. Mugel, and E. Lizaso, Reviews in Physics 4, 1 (2019), arXiv:1807.03890 .
  * Souza  _et al._ [2022] A. M. Souza, E. O. Martins, I. Roditi, N. Sá, R. S. Sarthour, and I. S. Oliveira, Frontiers in Physics 9, 1 (2022), arXiv:2005.02846 .
  * Mugel  _et al._ [2022] S. Mugel, C. Kuchkovsky, E. Sánchez, S. Fernández-Lorenzo, J. Luis-Hita, E. Lizaso, and R. Orús, Physical Review Research 4, 1 (2022), arXiv:2007.00017 .
  * Sharabiani  _et al._ [2021] M. T. A. Sharabiani, V. B. Jakobsen, M. Jeppesen, and A. S. Mahani,  , 1 (2021), arXiv:2105.11322 .
  * Urgelles  _et al._ [2022] H. Urgelles, P. Picazo-martinez, D. Garcia-roger, and J. F. Monserrat, (2022).
  * Luckow  _et al._ [2021] A. Luckow, J. Klepsch, and J. Pichlmeier, Digitale Welt 5, 38 (2021), arXiv:2103.07433 .
  * Inoue  _et al._ [2021] D. Inoue, A. Okada, T. Matsumori, K. Aihara, and H. Yoshida, Scientific Reports 11, 1 (2021), arXiv:2003.07527 .
  * IBM Q team [2021] IBM Q team, “Qiskit: An open-source framework for quantum computing,”  (2021), https://doi.org/10.5281/zenodo.2573505, release 0.39.4.
  * D-Wave Systems [2022] D-Wave Systems, “D-Wave Ocean SDK,”  (2022), https://github.com/dwavesystems/dwave-ocean-sdk, release 6.2.0.
  * Sharma  _et al._ [2022] V. Sharma, N. S. B. Saharan, S.-H. Chiew, E. I. R. Chiacchio, L. Disilvestro, T. F. Demarie, and E. Munro,  , 1 (2022), arXiv:2210.08695 .
  * De Raedt  _et al._ [2007] K. De Raedt, K. Michielsen, H. De Raedt, B. Trieu, G. Arnold, M. Richter, Th. Lippert, H. Watanabe, and N. Ito, Comput. Phys. Commun. 176, 121 (2007).
  * De Raedt  _et al._ [2019] H. De Raedt, F. Jin, D. Willsch, M. Willsch, N. Yoshioka, N. Ito, S. Yuan, and K. Michielsen, Comput. Phys. Commun. 237, 47 (2019).
  * Willsch  _et al._ [2022b] D. Willsch, M. Willsch, F. Jin, K. Michielsen, and H. De Raedt, Comput. Phys. Commun. 278, 108411 (2022b).
  * Glover  _et al._ [2019] F. Glover, G. Kochenberger, and Y. Du, 4or 17, 335 (2019).
  * Ver [2022] Discrete Optimization 44, 100594 (2022).
  * Grötschel and Nemhauser [2008] M. Grötschel and G. L. Nemhauser, Discrete Optimization 5, 168 (2008).
  * Heßler  _et al._ [2020] K. Heßler, S. Irnich, T. Kreiter, and U. Pferschy,  (2020).
  * Yan  _et al._ [2022] J. Yan, Y. Lu, L. Chen, S. Qin, Y. Fang, Q. Lin, T. Moscibroda, S. Rajmohan, and D. Zhang,  _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, 2022, Washington, DC, USA_, Vol. 1 (Association for Computing Machinery, 2022) arXiv:2207.11122 .
  * Kroes  _et al._ [2020] M. Kroes, L. Petrica, S. Cotofana, and M. Blott, GECCO 2020 - Proceedings of the 2020 Genetic and Evolutionary Computation Conference , 1125 (2020), arXiv:2003.12449 .
  * Kadowaki and Nishimori [1998b] T. Kadowaki and H. Nishimori, Physical Review E 58 (1998b).
  * Johnson  _et al._ [2011] M. Johnson, M. Amin, S. Gildert, and et al., Nature 473, 194 (2011).
  * Wang  _et al._ [2019] Z. Wang, N. C. Rubin, J. M. Dominy, and E. G. Rieffel,  (2019), 10.1103/PhysRevA.101.012320, arXiv:1904.09314 .
  * [45] D. Willsch, M. Jattana, M. Willsch, S. Schulz, F. Jin, H. De Raedt, and K. Michielsen, in  _NIC Symposium 2022_, Publication Series of the John von Neumann Institute for Computing (NIC) NIC Series, Vol. 51, edited by M. Müller, C. Peter, and A. Trautmann (Forschungszentrum Jülich GmbH Zentralbibliothek, Verlag, Jülich) pp. 31–44.
  * Kesselheim  _et al._ [2021] S. Kesselheim, A. Herten, K. Krajsek, J. Ebert, J. Jitsev, M. Cherti, M. Langguth, B. Gong, S. Stadtler, A. Mozaffari, G. Cavallaro, R. Sedona, A. Schug, A. Strube, R. Kamath, M. G. Schultz, M. Riedel, and Th. Lippert, “JUWELS Booster – A Supercomputer for Large-Scale AI Research,” (2021), arXiv:2108.11976 [cs.DC] .
  * Montanez-Barrera  _et al._ [2023] J. A. Montanez-Barrera, P. van den Heuvel, D. Willsch, and K. Michielsen,  (2023), arXiv:2305.18757 .