TRSP

    _transport robot scheduling problem_
QC

    quantum computing
QUBO

    quadratic unconstrained binary optimization
GPGPU

    general purpose computation on graphics processing unit
MIP

    mixed integer program
LP

    linear programming
AGV

    automated guided vehicles
JSSP

    job shop scheduling problem
QAOA

    quantum approximate optimization algorithm
MIS

    maximum independent set
BQM

    binary quadratic model
CQM

    constrained quadratic model
LBQM

    D-Wave’s hybrid _Leap_ framework
FDA

    Fujitsu’s digital annealer
FDAh

    Fujitsu’s digital annealer hybrid framework
TSB

    Toshiba’s Simulated Bifurcation Machine
API

    application programming interface
NISQ

    noisy intermediate-scale quantum
SE-GU

    Gurobi with the sequence model
TI-GU

    Gurobi with the time-indexed model
QU-LBQM

     D-Wave’s hybrid _Leap_ framework (LBQM) with the QUBO model
QU-FDA

     Fujitsu’s digital annealer (FDA) with the QUBO model
QU-FDAh

     Fujitsu’s digital annealer hybrid framework (FDAh) with the QUBO model
TI-GU@QU-LBQM

    the relative runtime of TI-GU w.r.t. QU-LBQM
SE-GU@QU-LBQM

    the relative runtime of SE-GU w.r.t. QU-LBQM
TI-GU@QU-FDA

    the relative runtime of TI-GU w.r.t. QU-FDA
SE-GU@QU-FDA

    the relative runtime of SE-GU w.r.t. QU-FDA
SE-GU@QU-FDAh

    the relative runtime of SE-GU w.r.t. QU-FDAh

# A Quantum Optimization Case Study for a Transport Robot Scheduling Problem

Dominik Leib1, Tobias Seidel1, Sven Jäger1, Raoul Heese1, Caitlin Jones2,
Abhishek Awasthi2, Astrid Niederle3, Michael Bortz1 D. Leib:
<EMAIL_ADDRESS>Awasthi<EMAIL_ADDRESS>
( 1 Fraunhofer ITWM, Kaiserslautern, Germany
2 BASF Digital Solutions GmbH, Ludwigshafen am Rhein, Germany
3 BASF SE, Ludwigshafen am Rhein, Germany
)

###### Abstract

We present a comprehensive case study comparing the performance of D-Waves’
quantum-classical hybrid framework, Fujitsu’s quantum-inspired digital
annealer, and Gurobi’s state-of-the-art classical solver in solving a
transport robot scheduling problem. This problem originates from an
industrially relevant real-world scenario. We provide three different models
for our problem following different design philosophies. In our benchmark, we
focus on the solution quality and end-to-end runtime of different model and
solver combinations. We find promising results for the digital annealer and
some opportunities for the hybrid quantum annealer in direct comparison with
Gurobi. Our study provides insights into the workflow for solving an
application-oriented optimization problem with different strategies, and can
be useful for evaluating the strengths and weaknesses of different approaches.

Keywords— Optimization, Quantum Computing, Robot Scheduling

##  1 Introduction

QC is a field that has witnessed a rapid increase in interest and development
over the past few decades since it was theoretically shown that quantum
computers can provide an exponential speedup for certain tasks [19, 26, 52].
Translating this potential into a practically relevant quantum advantage,
however, has proven to be a very challenging endeavor. Nevertheless, the
emerging field is considered to have a highly disruptive potential for many
domains, for example in machine learning [50], chemical simulations [13] and
optimization [35], the domain of this work. Due to the fact that optimization
problems are of utmost importance also for industrial applications, we
investigated a potential advantage of quantum and quantum-inspired technology
for the so-called _transport robot scheduling problem_ (TRSP), a real-world
use-case in optimization that is derived from an industrial application of an
automatized robot in a high-throughput laboratory. The optimization task is to
plan a time-efficient schedule for the robot’s movements as it transports
chemical samples between a rack and multiple machines to conduct experiments.
This is an NP-hard problem which for certain instances can be challenging to
solve using classical computing techniques, and hence is an attractive
candidate to search for an advantage with non-classical techniques.

In our study, we compared the solution quality and runtime of different
solvers on a large set of instances of the problem. As solvers, we considered
LBQM that makes use of the D-Wave quantum annealer [17], FDA[41], FDAh, as
well as the industry-grade Gurobi solver [27]. As a key element of this work,
we provide three different models for the TRSP that follow different design
philosophies. This is justified by the different ways in which the problem
task can be modelled and the inherent differences in the problem formulations
that the solvers addressed can accept. LBQM, FDA and FDAh are restricted to a
formulation as a quadratic unconstrained binary optimization (QUBO), whereas a
mixed integer program (MIP) with integer and float variables can be used by
Gurobi, which makes a comparison of multiple formulations meaningful.

The TRSP considered in this paper is a special combination of different
scheduling problems that, to our knowledge, has not been considered before.
Scheduling problems have been studied intensively for several decades and
classical algorithms exist for numerous variants [11, 45]. Since most of the
industry-relevant scheduling problems are NP-hard, these classical algorithms
mainly consist of meta-heuristics or use general-purpose MIP solvers, which
basically solve the problem using a branch and bound approach with several
additional improvements like cutting planes. In addition to classical
algorithmic developments, a considerable amount of research has also been done
in hardware-based parallel computing, especially in general purpose
computation on graphics processing unit (GPGPU) parallelization [16, 7]. The
problem discussed in this work is an extension of the typical job shop
scheduling problem (JSSP), where the inclusion of a robot adds additional
restrictions. More specifically, the studied scheduling problem falls into the
category of robotic cell scheduling and automated guided vehicles (AGV)
scheduling problems. Most work on robotic cell scheduling deals with infinite
cyclic schedules [18]. This comprises polynomial-time algorithms and hardness
results [53], MIP techniques [44, 12, 22] and heuristic approaches [36]. Many
efficiently solvable and hard special cases have been identified [51] and
heuristics have been proposed for some of the hard cases [54]. Those problems
differ from our use case in one way or another. The problems considered by the
above mentioned papers allow, unlike our use case, that the jobs can wait at a
machine after their completion before being picked up by the robot. Robotic
cell scheduling problems without this possibility have been studied by [1, 2],
whose problems differ from our, among others, in the considered objective
function. Our objective function, the total job completion time, has been
extensively studied for flow shop scheduling problems without a robot [45, 28,
5, 48], the latter of which shows that the no-wait variant is strongly NP-hard
on two machines. Apart from the no-wait constraint, the problem considered in
our work is characterized by the fact that jobs have to go to the last machine
several times. Such settings are known as a re-entrant flow shops, for which
[31] developed a heuristic algorithm.

We are mainly interested in the performance of non-standard solution
approaches using quantum or quantum-inspired solvers in this study. Because
these solvers rely on heuristics, benchmarks for real-world applications are a
highly relevant research topic. Most quantum optimization approaches fall into
two major groups, one for gate-based hardware and one for annealing-based
hardware [4]. The majority of gate-based approaches to optimization use
parameterized gates to find the ground state of a Hamiltonian related to the
cost function of the optimization problem in a quantum-classical hybrid
fashion, for example via the quantum approximate optimization algorithm (QAOA)
[21, 9].

Approaches based on quantum annealing also seek to find the ground state of a
Hamiltonian, but by aiming for an adiabatic change from an initial state that
can be easily prepared. In contrast to actual quantum computing devices, other
classical software and hardware components are merely _inspired_ by quantum
computing, for example FDA [6] and Toshiba’s Simulated Bifurcation Machine
(TSB) [56]. Typically, optimization tasks for quantum solvers and the
aforementioned quantum-inspired technologies are modeled as QUBO problems
[33]. An in-depth analysis of pure QUBO comparison on four quantum and
quantum-inspired solvers can be found in [43]. In their work, the authors
compare the solutions of a library of quadratic benchmark problems on the
D-Wave quantum annealer, FDA, and TSB against each other.

QC has already been successfully used for optimization in various fields. For
example, in [40], chemical reaction networks are optimized with quantum
computing. In [55], it is shown that using the QAOA, it is possible to beat
some classical heuristic algorithms on the binary paint shop problem. However,
some work has shown that the current circuit model algorithms are not always
adequate enough to reach significant convergence required for a good solution
[8]. Quantum annealing has proven to offer some advantage against the
classical simulated annealing algorithm for a spin-glass problem, using D-Wave
hardware [47], but this is no conclusive evidence. In one of the more recent
works on quantum annealing [49], the authors suggest a nature inspired hybrid
quantum algorithm for robot trajectory optimization for PVC sealing in a real
industrial setting. In [20], the authors present a solution to the maximum
independent set (MIS) problem using a Rydberg atom device, along with a claim
of a possible super-linear quantum speed-up against classical simulated
annealing. Other classical algorithms might still be superior to a quantum
approach on current devices [3]. Several works consider scheduling problems
[61, 14, 57]. In [23], an AGV transportation problem using different classical
and quantum approaches is studied and [30] investigates a nurse scheduling
problem with the usage of a quantum annealer.

The remaining manuscript is structured as follows. We provide a detailed
description of the TRSP and its mathematical modeling in Section 2. In Section
3, we describe the design of our numerical study and list the problem
instances and solvers that we use. The results of this study are presented in
Section 4. Finally, we conclude our study in Section 5.

##  2 Transport Robot Scheduling Problem

In this section, we present a detailed explanation of the TRSP, which is a
real-world use case derived from one of BASF’s high-throughput laboratories.
This optimization problem is about finding the most time-efficient route of a
transport robot tasked with moving chemical samples from one processing
machine to another, as illustrated in Fig. 1. In the following, we first
provide a general description of the problem setup and then present different
modeling approaches. These models build the foundation of the subsequent
benchmarks.

Figure 1: Sketch of the high-throughput laboratory setup for the TRSP. The
laboratory consists of a sample rack and different processing machines (water
mixer, sample shaker and photo booth). The arrows indicate how the samples may
be transported by the robot between the machines and the rack, one at a time.

###  2.1 Problem Description

The laboratory we are modeling here consists of a _sample rack_ and three
different processing machines: a _water mixer_ , a _sample shaker_ and a
_photo booth_. And, finally, the _robot_ itself that is tasked with carrying
chemical samples from one place to another with the goal to conduct chemical
experiments. Only the experimental plan (i. e., how each sample has to be
processed in the laboratory) is predefined in advance, but not the specific
order of the experiments. Initially, a certain number of samples is stored on
the rack. Each of these samples needs to be first taken to the water mixer,
then to the sample shaker. Once the sample shaking is completed, one or more
photos have to be taken of each sample at the photo booth. Consecutive photos
need to be taken after specific (i. e., predefined) time intervals, where the
first photo of each sample has to be taken immediately after the shaking
process. Finally, each sample has to be brought back to the rack. The
processing times for different samples on the same machine can be different as
specified by the experimental plan. We assume that each machine can only hold
(and process) one single sample at any given time, or remain idle, and the
processing steps cannot be interrupted before their completion. It is required
that a machine starts processing a sample as soon as the sample is brought by
the robot. Moreover, we assume that a sample has to be moved by the robot in-
between two processing steps. Hence, a sample has to be lifted from a machine
(and the machine is made available) as soon as it finishes processing.

By definition, the robot requires exactly one time unit to move from any place
to any other, with or without a sample, and picking up or dropping a sample
does not require extra time. Like the machines, the robot can transport only a
single sample at any given time or drive empty or remain idle. In particular,
it is not possible that the robot places a sample at a machine and picks up
another at the same time.

The objective of this scheduling task is to minimize the _sum of sample
completion times_ , i. e., the sum of the times when the samples arrive at the
the rack after their last photo has been taken. The solution of this
optimization problem is a sequence of tasks for the robot that yields an
efficient laboratory operation.

###  2.2 Mathematical Modeling

In our benchmark, we test three modeling approaches against each other. On the
quantum and quantum-inspired side we consider a QUBO formulation, whereas on
the classical side we use two MIP formulations. First, a so-called _sequence
model_ and second, a so-called _time-indexed model_. In the following, we
first introduce the common terminology for all modeling approaches. Next, we
shortly sketch the main features of each model. For a more detailed
description, we refer to Appendix A. The motivation for the development of
multiple models is to carry out a comparison between the solutions obtained by
the most suitable problem encoding for quantum and classical solvers. This
ensures that we are comparing the best of both worlds (classical and quantum),
and do not restrict ourselves to a model which is more suitable for quantum
over classical computing.

####  2.2.1 Common Terminology

The processing machines are addressed by $M_{1}$ for the water mixer, $M_{2}$
for the sample shaker and $M_{3}$ for the photo booth. The scheduling time is
discretized into time slots which all have length of one time unit. The
transport robot takes one time unit for each operation that is either
transportation or empty traversal between the machines and the rack. In this
way, each transport robot scheduling problem is uniquely determined by the
following parameters:

  * •

The number of samples to be scheduled $N\geq 1$.

  * •

The number of photos $K\geq 1$, which agrees for each sample
$j\in\\{1,\ldots,N\\}$.

  * •

The processing times $p_{j,1},p_{j,2},p_{j,3}\in\mathbb{N}_{>0}$ for machines
$M_{1},M_{2}$ and $M_{3}$, which can vary for each sample
$j\in\\{1,\ldots,N\\}$.

  * •

The time gaps $g_{j,k}\in\mathbb{N}_{\geq 2}$ to be kept between consecutive
photos $k$ and $k+1$ for $k\in\\{1,\ldots,K-1\\}$, which also can vary for
each sample $j\in\\{1,\ldots,N\\}$.

As an example, Fig. 2 provides a feasible schedule in form of a Gantt chart to
visualize these parameters.

Figure 2: An example Gantt chart of a robot transport scheduling problem (as
sketched in Fig. 1) with $N=2$ samples and $K=2$ photos. When a sample is
processed on one of the machines or carried by the robot in the time-frame
$[t,t^{\prime}]$, a bar is drawn from $t$ to $t^{\prime}$ in the respective
row in a corresponding color. Empty movements of the robot are not drawn
explicitly. For example, at time $t=13$ the robot is at the rack as sample $1$
has been brought to the sample rack from $t=12$ to $13$. It takes one unit of
time for the robot to travel from the rack to the water mixer to pick up
sample $2$ at $t=14$. From $t=22$ to $t=23$, the sample is brought from the
photo booth to the rack and back, which is a consequence of the assumption
that a sample has to be moved by the robot in-between two processing steps.
The objective value of the depicted schedule is $19+26=45$.

####  2.2.2 QUBO Model

A general QUBO reads

$\displaystyle\min_{x}$ $\displaystyle x^{\top}\cdot Q\cdot x$ (1) s.t.
$\displaystyle x\in\\{0,1\\}^{n}$

for some matrix $Q\in\mathbb{R}^{n\times n}$, where $x$ represents a vector of
$n$ binary optimization variables. Two challenging properties of QUBOs must be
taken into account in the modeling. Since only binary variables are allowed,
this implies that other types of variables must be avoided, i. e. a
reformulation into a binary form is necessary. Second, the problem is
unconstrained. This restriction can be overcome by using _penalty terms_ ,
which are quadratic functions in the model variables that evaluate to a
positive value when the current assignment of values to the variables leads to
an infeasible solution. Typically, the penalty terms are designed to yield $0$
if the corresponding solution is feasible, so that they do not contribute to
the objective values of feasible solutions. More general information about
QUBOs and their properties can be found, e. g., in Ref.[33, 38, 24].

Our proposed QUBO model for the TRSP is based on the well-known starting time
formulation (see e. g. Ref.[14]) and can be written as

$\displaystyle\min_{x}$
$\displaystyle\rho_{0}F(x)+\sum_{i=1}^{7}\rho_{i}P_{i}(x)$ (2) s.t.
$\displaystyle x\in\\{0,1\\}^{n},$

where $F$ is the objective function and $P_{1},\ldots,P_{7}$ denote the
penalty functions and $\rho_{0},\ldots,\rho_{7}\in\mathbb{R}_{>0}$ are tunable
parameters that have to be chosen such that the objective and penalty terms
are suitably balanced. As in equation 1, $n$ represents the total number of
binary optimization variables. These have a distinct meaning that can be
identified with three indices. Specifically,

$x_{j,m,t}:=\begin{cases}1,&\text{if sample $j$ starts processing on machine
$M_{m}$ at time $t$},\\\ 0,&\text{otherwise}\end{cases}$ (3)

for all $j\in\\{1,\ldots,N\\}$, $m\in\\{1,2,3\\}$ and
$t\in\\{1,\ldots,T-1\\}$. Here, $T$ denotes the time horizon, which is chosen
in such a way that there is enough time to schedule all samples sequentially,
implying that there is at least one feasible solution. It can be explicitly
computed for each instance as described in Section A.1. In terms of Fig. 2,
one has, for example, $x_{1,1,1}=1$ and $x_{1,2,8}=1$.

The penalty terms for the QUBO model have to be formulated using the binary
optimization variables. This section only provides an example for such a term,
a complete description can be found in Section A.1. Specifically, we consider
here the constraint that each sample must access the machines $M_{1}$ and
$M_{2}$ exactly once, which can be achieved by

$P_{1}:=\sum_{j=1}^{N}\sum_{m=1}^{2}\left[\left(\sum_{t=1}^{T-1}x_{j,m,t}\right)-1\right]^{2}\;.$
(4)

This term evaluates to zero if and only if for each pair of sample $j$ and
machine $M_{m}$, the variable $x_{j,m,t}$ is $1$ for precisely one time slot
$t$. Since $P_{1}$ is bounded below by $0$ due to its quadratic nature, each
local minimum of $P_{1}$ is a feasible solution w.r.t. the rule of machine
access to $M_{1}$ and $M_{2}$. The other penalty terms can be formulated
similarly.

Finally, the objective function $F$ sums up for each sample the time when the
sample arrives at the rack after the entire scheduling process (”sum of sample
completion times“). For example, the objective function in the case of Fig. 2
evaluates to $45$ time units.

####  2.2.3 MIP Models

A general MIP reads

$\displaystyle\min_{x,y}$ $\displaystyle c_{1}^{\top}x+c_{2}^{\top}y$ (5) s.t.
$\displaystyle A_{1}x+A_{2}y\leq b$ $\displaystyle C_{1}x+C_{2}y=d$
$\displaystyle x\in\mathbb{Z}^{n},y\in\mathbb{R}^{m},$

where $A_{1}\in\mathbb{Z}^{k\times n}$, $C_{1}\in\mathbb{Z}^{l\times n}$,
$A_{2}\in\mathbb{Z}^{k\times m}$, $C_{2}\in\mathbb{Z}^{l\times m}$,
$b\in\mathbb{Z}^{k}$, $d\in\mathbb{Z}^{l}$, $c_{1}\in\mathbb{Z}^{n}$ and
$c_{2}\in\mathbb{Z}^{m}$ for $n,m,k,l\in\mathbb{N}_{>0}$. MIP have been used
since the late 1950s as a tool for solving scheduling problems. It is not
possible to model the disjunctive constraints resulting from the discrete
ordering decisions only by means of starting time variables. Different types
of binary variables have been proposed to achieve this. The main types are
position variables $x_{ijk}$ indicating if job $j$ is the $k$th job on machine
$i$ [58], linear ordering variables $\delta_{ijk}$ deciding if job $j$ is
processed before job $k$ on machine $i$ [39] and time-indexed variables
$x_{ijt}$ specifying that job $j$ is started (or processed or completed) on
machine $i$ at time $t$ [10, 46]. Ref.[34] compared these three approaches
experimentally for a job shop scheduling problem.

Due to the powerful nature of (mixed) integer programming in contrast to the
restrictive nature of the QUBO models, we provide two MIP models to be solved
using Gurobi. The first one, in the following named _sequence model_ , makes
use of integer start time and binary linear ordering variables. The second
model, called the _time-indexed model_ , is restricted to a binary formulation
comparable to the QUBO model, where we make use of time-indexed variables.

####  2.2.4 MIP: Sequence Model

In the sequence model, we model sequences of _events_ that affect the behavior
of the transport robot with respect to the machines and the photos of a
sample. We define the _set of events_ as

$E:=\bigl{\\{}(j,i,a)\mid
j\in\\{1,\ldots,N\\},\,i\in\\{1,\ldots,2+K\\},\,a\in\\{0,1\\}\bigr{\\}}.$ (6)

An event $e=(j,i,0)$ represents either that a sample $j$ is placed on machine
$M_{i}$ for $i\in\\{1,2\\}$ or to the $(i-2)$th photo shoot for $i>2$, an
event $(j,i,1)$ corresponds to picking it up again. For each event $e\in E$ we
define an optimization variable $\tau_{e}\in\mathbb{R}_{\geq 0}$ to model the
time for event $e$ to happen. In terms of Fig. 2, we have, for example,
$\tau_{(1,1,0)}=1$ and $\tau_{(1,1,1)}=4$. A simple formulation can be
achieved by additionally introducing a binary variable for each pair $e,f\in
E$, $e\neq f$ of events that indicates if $e$ occurs before $f$. We reduce the
size of the model by exploiting the fact that the ordering of some events is
fixed or coupled. For example, we do not need a variable that specifies the
order in which a given sample is brought to the water mixer and to the sample
shaker. This leads to three sets of linear ordering variables that can be
found in Section A.2, as well as the various constraints to ensure
feasibility. The objective function (i. e., the sum of the sample completion
times) can be easily expressed using the variables $\tau_{e}$ corresponding to
events when a sample is picked up from the last photo.

####  2.2.5 MIP: Time-Indexed Model

The second constrained model makes use of discrete time-indexed variables
similar to the QUBO model from Section A.1. In this formulation, we model the
behavior of the transport robot by defining certain routes a sample can be
transported along, which include those from the rack to all machines and back
or movements between subsequent machines. The numbering of the moves is shown
in Fig. 3.

Figure 3: Numbering of the eight different routes the transport robot can move
along between the machines and the rack as sketched in Fig. 1.

As the model name implies, we have, given a discrete time horizon
$T\in\mathbb{N}_{>0}$, binary variables to model when each sample takes which
route as

$y_{j,r,t}:=\begin{cases}1,\text{ if sample $j$ is transported by the robot on
route $r$ during the time $(t,t+1)$}\;,\\\ 0,\text{ otherwise}\;\end{cases}$
(7)

for all $j\in\\{1,\ldots,N\\}$, $r\in\\{1,\ldots,8\\}$ and
$t\in\\{0,\ldots,T-1\\}$. In terms of the Gantt chart from Fig. 2, this would
imply $y_{1,1,0}=1$, $y_{1,2,4}=1$, $y_{2,1,5}=1$ and so on. The time horizon
$T$ is defined as for the QUBO model, see equation 9.

The constraints of the model are similar to the penalty terms of the QUBO
Model and are listed in Section A.3. The objective function (i. e., the sum of
the sample completion times) is defined in terms of the ancilla optimization
variables $z_{j}$ for $j\in\\{1,\dots,N\\}$, that are bounded below by the
arrival time of sample $j$ at the rack after the schedule has finished.

##  3 Benchmark Setup

In the present section, we describe the design of the benchmark. We start with
an outline of the considered problem instances that are listed in more detail
in Appendix B. Subsequently, we describe the three different commercial
technologies that we use.

###  3.1 Instances

To set the stage for our benchmark, we specify $260$ test instances of our
optimization problem of interest, each defined by a different set of
parameters. Specifically, each instance is uniquely determined by the number
of samples $N$, the number of photos $K$, the gaps $g_{j,k}$ between
subsequent photos $k$ and $k+1$ for $k\in\\{1,\ldots,K-1\\}$ and
$j\in\\{1,\ldots,N\\}$, and, finally, the processing times
$p_{j,1},p_{j,2},p_{j,3}$ of the water mixer, sample shaker and photo booth,
respectively, as explained in Section 2.2. For the sake of simplicity, the
processing time of the photo booth agrees for all samples of the same
instance, that is $p_{j,3}:=p_{3}$ for all $j\in\\{1,\ldots,N\\}$.

In Appendix B, we describe the algorithm that was used to generate parameter
sets for the benchmark instances. Since the resulting instances span a wide
range of complexity, we divide the resulting benchmark library into two parts,
where each part is defined by the number of binary variables in the
corresponding QUBO formulation from Section 2.2.2 as explained in Section A.1
in more detail. The first part, which we call _library of minor instances_ ,
contains all $161$ instances that have at least $2071$ and at most $8080$
binary variables. The second part, which we call _library of major instances_
, contains the remaining $99$ instances with at least $10\,822$ and at most
$22\,692$ binary variables. The reason for that specific division is that
$8192$ is the maximal amount of variables that can be solved directly on
Fujitsu’s digital annealer.

We collect groups of instances $(N,K)$ that have the same number of samples
and photos as shown in the boxplots in Fig. 4, i. e., within those groups the
leftover parameters $p_{j,m}$ and $g_{j,k}$ for
$j\in\\{1,\ldots,N\\},m\in\\{1,2,3\\}$ and $k\in\\{1,\ldots,K-1\\}$ may vary.
These groups can be understood as a collection of similar TRSPs. This approach
allows us to consider statistical metrics over several instances when we
compare models and solvers. In Appendix B, we list how many instances each
group contains.

Figure 4: The number of variables of each group $(N,K)$ of instances. The
minor instances are left to the red vertical line, the major instances to the
right.

###  3.2 Quantum and classical solvers

In our benchmark, we solve the generated instances with a selection of model
and solver combinations with the main goal to assess the performance of
quantum and quantum-inspired technology. Specifically, we consider three
solver candidates:

  1. 1.

Gurobi: As a baseline, we use the branch and bound algorithm of Gurobi, which
is a state-of-the-art mathematical programming solver running on classical
hardware [27]. In summary, it relies on an implicit enumeration that allows
the original problem to be split into smaller sub-problems using a decision
tree. The use of lower bounds derived from linear programming (LP) relaxations
allows for a reduction of the search space. In this work we utilized the cloud
based service of Gurobi solver, which ran on a Intel(R) Xeon(R) Platinum
8275CL CPU (3.00 GHz with 8 physical cores).

  2. 2.

D-Wave’s hybrid _Leap_ framework (LBQM): D-Wave provides cloud-based access to
their adiabatic quantum computers with over $5000$ qubits [17]. By design,
their hardware is specifically tailored to solve QUBOs. To this end, the QUBO
is encoded in an Hamiltonian such that each optimization variable is
represented by one qubit [62] and the ground state corresponds to the optimal
solution. The quantum annealing mechanism aims to find the ground state by
performing a suitable time evolution of the quantum system with a subsequent
measurement of all qubits to reveal the optimal solution.

The D-Wave hardware has only limited connectivity, which means that each qubit
can only interact with a certain number of other qubits. This limitation
restricts the correlations between optimization variables that can be
represented by the Hamiltonian. Finding a suitable representation with these
constraints is an NP-hard problem [37] that has to be solved classically to
configure the quantum annealer for a certain problem. In practice, the quantum
annealer can typically only be used for QUBOs with much less than $5000$
optimization variables.

For this reason, D-Wave also provides a hybrid software framework LBQM, which
is a black-box algorithm for binary quadratic models that runs on both
classical and quantum annealing hardware. It allows to handle larger
optimization problems that are too big for the quantum hardware by presenting
only parts of the original problem to the quantum annealer. However, the exact
mode of operation of LBQM is not publicly available. In this study, we use
only the quantum annealer in a hybrid fashion via LBQM. The quantum machine
used in the hybrid framework is the _D-Wave Advantage System 4.1_ and the
region _na-west-1_. We choose to use a constant number of $1000$ samples (or
readouts) for all evaluations and use default settings for all parameters.

  3. 3.

Fujitsu’s digital annealer (FDA) and Fujitsu’s digital annealer hybrid
framework (FDAh): The digital annealer from Fujitsu can be considered as a
quantum-inspired algorithm that runs on dedicated (classical) hardware [6] and
can be accessed using a cloud service. It is based on simulated annealing [32,
15] with two major differences. Firstly, the utilization of an efficient
parallel-trial scheme to exploit the parallelization capabilities of the
hardware and, secondly, a dynamic escape mechanism to avoid locally optimal
solutions. The detailed hardware specifications are confidential. The solver
supports QUBOs with up to $8192$ variables.

In addition, the hybrid solver FDAh is provided to solver bigger problem
instances by utilizing both dedicated and classical hardware [41] similar to
D-Wave’s LBQM. In this study, we use both FDA and FDAh. Both solvers require a
set of parameters that specify how the annealing is done, which also include
the number of repetitions and parallel runs on the chip. The specific
parameters we used for FDA and FDAh are provided in Appendix C.

In a small pre-study, we excluded a few other solvers; see Appendix D.

Each instance can be modelled with each of the three modeling approaches from
Section 2. However, not all solvers are applicable to all problem formulations
and all instances. The MIP sequence model is solved with Gurobi for all
instances. The time-indexed model is solved with Gurobi only for the minor
instances. The QUBO model is solved with LBQM and FDA for minor instances. For
major instances, the QUBO model is only solved with FDAh.

We call each valid model and solver combination an _approach_ and use a unique
name to refer to it. Summarized, we consider Gurobi with the sequence model
(SE-GU), Gurobi with the time-indexed model (TI-GU),  with the  model (QU-
LBQM),  with the  model (QU-FDA) and  with the  model (QU-FDAh). An overview
over all approaches is shown in Fig. 5.

Figure 5: Summary of model (see Section 2.2) and solver (see Section 3.2)
combinations for the benchmarks.

For all problems, we prescribe a runtime limit of $3600$ seconds for Gurobi.
Both LBQM and FDA also require a time limit for each run, which scales with
the problem size in the QUBO formulation as follows. The time limit for LBQM
is set to be $\min\\{100,1.5\cdot\frac{n}{100}\\}$ seconds, where $n$ is the
number of variables in the QUBO formulation for the minor instances. The
runtime of the digital annealer is implicitly set with the _steps_ parameter,
where each step taken in the annealing process takes a constant amount of
time. We set the number of steps to be $1\text{$\cdot$}{10}^{7}$ for the
instances with $$2071$\leq n\leq$4096$$, $5\text{$\cdot$}{10}^{7}$ for the
ones with $$4096$<n\leq$6000$$ and $1\text{$\cdot$}{10}^{8}$ for the instances
with $$6000$<n\leq$8080$$ variables in the QUBO formulation. Lastly the major
instances computed with the hybrid framework FDAh based on the digital
annealer require a time limit as well. For this we distributed the available
time of $5$ hours to the instances, correspondingly to their number of
variables. This computes approximatively as $n\cdot$0.0117$$ seconds where $n$
is the number of variables in the QUBO formulation.

The benchmark setup is summarized in Table 1, where we recall the approaches
from Fig. 5. The table also contains the values of the QUBO parameters
$\rho_{0},\ldots,\rho_{7}$ from equation 23 that were chosen for LBQM, FDA and
FDAh, respectively. The choice was made according to previous experiments with
smaller problem instances. For this purpose, a typical strategy is to
iteratively increase the parameter $\rho_{i}$ if the corresponding penalty
term $P_{i}$ is non-vanishing. Additionally, one needs to make sure that the
parameter $\rho_{0}$ for the target function is set such that it is not in
favor to violate penalty terms and a good optimization is achieved.

Table 1: Benchmark setup: Summary of problem instances from Section 3.1 and solvers from Section 3.2 for the optimization problems (or models) from Section 2. Property | Minor instances | Major instances
---|---|---
Number of instances | $161$ | $99$
Number of variables $(n)$ |  $2071$ to $8080$ |  $10\,822$ to $22\,692$
Approach (cf. Fig. 5) | Used for minor instances | Used for major instances
SE-GU | ✓ | ✓
TI-GU | ✓ | ✗
QU-LBQM | ✓ | ✗
QU-FDA | ✓ | ✗
QU-FDAh | ✗ | ✓
Approach (cf. Fig. 5) | Minor instance limit | Major instance limit
SE-GU | $3600\text{\,}\mathrm{s}$ | $3600\text{\,}\mathrm{s}$
TI-GU | $3600\text{\,}\mathrm{s}$ | —
QU-LBQM | $\min\\{$100$,$1.5$\cdot\frac{n}{$100$}\\}\cdot$1\text{\,}\mathrm{s}$$ | —
QU-FDA | $\begin{cases}$1\text{$\cdot$}{10}^{7}$\text{ iterations,}&$2071$\leq n\leq$4096$\\\ $5\text{$\cdot$}{10}^{7}$\text{ iterations,}&$4096$<n\leq$6000$\\\ $1\text{$\cdot$}{10}^{8}$\text{ iterations,}&$6000$<n\leq$8080$\end{cases}$ | —
QU-FDAh | — | $n\cdot$0.0117\text{\,}\mathrm{s}$$
Solver |  QUBO parameters from equation 23
LBQM | $\rho_{0}=$1$,\rho_{1}=$30\,000$,\rho_{2}=\rho_{3}=\rho_{4}=\rho_{5}=\rho_{7}=$10\,000$,\rho_{6}=$15\,000$$
FDA | $\rho_{0}=$1000$,\rho_{1}=$4000$,\rho_{2}=\rho_{3}=$1000$,\rho_{4}=\rho_{5}=\rho_{6}=\rho_{7}=$1500$$
FDAh | $\rho_{0}=$1000$,\rho_{1}=$2000$,\rho_{2}=\rho_{3}=$500$,\rho_{4}=\rho_{5}=\rho_{6}=\rho_{7}=$750$$

##  4 Benchmark Results

In the current section, we present the results of our previously described
benchmark, which is summarized in Table 1. For this purpose, we first show the
results for the minor instances and subsequently the results for the major
instances.

###  4.1 Results for Minor Instances

In Fig. 6, we show the objective values and runtimes of several approaches as
boxplots. All runtimes are end-to-end runtimes, that is, we consider the
entire evaluation pipeline, beginning with the submission of the problem to
the solver and ending with the return of a solution, including potential
network delays. The programmatic construction of the optimization problem for
the application programming interface (API) of the solver based on the
instance data is not part of the runtime.

From Fig. 6(a), we can observe that both the SE-GU and TI-GU solutions reach a
better objective value than the solutions from QU-LBQM and QU-FDA. When
comparing objective values, it has to be taken into account that the QUBO
model objective, equation 2, also includes penalty terms, which become
positive for infeasible solutions and therefore increase the objective value
accordingly. Specifically, we find that only QU-LBQM yields infeasible
solutions for some instances, whereas all other approaches yield feasible
solutions (SE-GU and TI-GU solutions are by definition always feasible). For
our analysis, we include both feasible and infeasible solutions. By performing
a Welch t-test [59], we find that the means of the results from both SE-GU and
TI-GU are lower than the means of the QU-FDA and QU-LBQM results with a
statistical significance of over $99\%$, respectively. The same holds for the
QU-FDA objective values in comparison to QU-LBQM.

(a) Objective values

(b) Runtimes

Figure 6: Benchmark results for minor instances as boxplots. The results are
grouped into sets of instances $(N,K)$ with the same number of samples $N$ and
photos $K$. A horizontal line marks the upper time limit of
$3600\text{\,}\mathrm{s}$ for Gurobi in Fig. 6(b). Some instances have not
been solved to feasibility by QU-LBQM, as indicated by the peaks above
$10^{4}$ in Fig. 6(a). Abbrevations according to Fig. 5.

On the other hand, according to Fig. 6(b), the computation time for TI-GU and
for some instances of SE-GU exceed the computation time of QU-LBQM and QU-FDA.
Since MIP solvers typically spend a lot of time proving that a solution is
optimal, we are also interested in the time taken by Gurobi (for both SE-GU
and TI-GU) to find solutions of the same quality as those obtained from QU-
LBQM or QU-FDA. Hence, we perform an additional analysis of the iterative
solver progress of each Gurobi run and look for the earliest computation time
at which Gurobi has reached an objective value that is less than or equal to
the corresponding objective value returned by the competing solvers for the
same instance. We call this earliest computation time the _relative runtime_.
Specifically, we consider the relative runtime of  w.r.t.  (@), the relative
runtime of  w.r.t.  (@), the relative runtime of  w.r.t.  (@) and the relative
runtime of  w.r.t.  (@). In the special case that Gurobi is not able to find
an objective value of the desired quality within its limit of $3600$ seconds
(which only occurs for some major instances), this time limit is used in place
of the earliest computation time. Exemplarily, we consider a specific instance
to visualize @ and @ in Fig. 7.

Figure 7: Visualization of the relative runtime of TI-GU w.r.t. QU-LBQM and
QU-FDA, denoted by @ and @, respectively. Here, we consider the example
instance $(7,4,3)(3)$; see supplementary material. The orange dots (connected
by lines for better visualization) mark the resulting objective values of TI-
GU at the corresponding time steps. The horizontal upper, blue and lower,
green line mark the final objective value of QU-FDA and QU-LBQM, respectively,
on the same instance. The blue and green lines intersect with the orange lines
at some point. The time coordinate of the next lower TI-GU objective value
after this intersection represents the relative runtime of TI-GU w.r.t. the
solver, which is marked as a vertical line in the corresponding color. In
other words, the relative runtime represents how long TI-GU has to run until
it reaches an objective value that is at least as good as the result from QU-
LBQM or QU-FDA, respectively.

The results of this analysis are shown in Fig. 8. This plot shows that QU-LBQM
is not able to compete with SE-GU. All problems from the first $4$ out of $9$
instance groups have been solved with SE-GU in under $1$ second while the
remaining instances in less than $10$ seconds, whereas the QU-LBQM runtimes
range between $50$ and $100$ seconds. However, LBQM finds a comparable
solution faster than TI-GU for most problems with $6$ or more samples and
remains competitive for smaller problems. A Welch-t test confirms that the
mean of TI-GU runtime is larger than the one of QU-LBQM runtime with a
significance over $99\%$.

Furthermore, Fig. 8(b) shows that QU-FDA is outperformed by SE-GU as well.
Analogous to Fig. 8(a), the instances in groups $(4,4),(5,3),(5,4)$ and
$(6,3)$ have been solved by SE-GU in $1$ second or less. But in contrast to
Fig. 8(a), the other groups have their median between $1$ second and $10$
seconds, i. e., which reflects that the target objectives from QU-FDA are
lower than those from QU-LBQM (see Fig. 6(a)). Nonetheless, the time taken for
SE-GU to reach the solution quality of QU-FDA is $10$ to $100$ times smaller.
Regarding TI-GU, QU-FDA finds a comparable solution almost always faster with
a few exceptions.

(a) Runtime comparison for @, @ and QU-LBQM from Fig. 6(b).

(b) Runtime comparison for @, @ and QU-FDA from Fig. 6(b).

Figure 8: Benchmark results for minor instances as boxplots. We show the
relative runtimes of TI-GU and SE-GU w.r.t. QU-LBQM and QU-FDA, denoted by @,
@, @ and @, respectively. The results are grouped into sets of instances
$(N,K)$ in analogy to Fig. 6. See Fig. 7 for an example of the relative
runtime computation. Abbrevations according to Fig. 5.

###  4.2 Results for Major Instances

The results for major instances are presented in analogy to the results for
minor instances from the previous section. In Fig. 9, we show the runtime and
the target value of the solvers on the corresponding models as boxplots.

(a) Runtimes

(b) Objective values

Figure 9: Benchmark results for major instances as boxplots. The results are
grouped into sets of instances $(N,K)$ as for previous plots. Abbrevations
according to Fig. 5.

The objective values of QU-FDAh are worse than the ones of SE-GU with a
significance of over $97\%$, but Fig. 9(a) shows that the runtime of SE-GU
increases strictly until it reaches the upper bound for the computation time
of $3600$ seconds, which happens for ca. $15$ samples. On the other hand, the
computation time of QU-FDAh ranges between $120$ and $300$ seconds, where only
a slight increase can be seen.

Analogously to Fig. 8(b), we evaluate the earliest computation times of SE-GU
model to reach objective values equal to or lower than the objective values
obtained from QU-FDAh, denoted by the relative runtime of  w.r.t.  (@). The
results are shown in Fig. 10.

Figure 10: Benchmark results for major instances as boxplots. We show the
relative runtime of SE-GU w.r.t. QU-FDAh, denoted by @, in analogy to Fig. 8.
We also show the runtime of QU-FDAh from Fig. 9(a). The results are grouped
into sets of instances $(N,K)$ as for previous plots. Abbrevations according
to Fig. 5.

In Fig. 10, a strictly increasing computation time can be seen for SE-GU,
whereas the QU-FDAh runtime remains almost constant. For the biggest instances
with $N=20$ samples, QU-FDAh has a clear advantage with respect to the
computation time, whereas it is competitive to SE-GU for the instances with
$15$ samples. In this sense QU-FDAh finds a solution of comparable quality
much faster for problems with $20$ samples than SE-GU and the latter was not
able to prove optimality for some of the instances with $20$ samples. A Welch
t-test confirms with a significance of over $99\%$ that the QU-FDAh mean is
lower than the @ mean.

##  5 Conclusion and Outlook

This paper presents a thorough benchmarking of an industrially relevant use
case of combinatorial optimization, the TRSP with the goal to achieve a time-
optimal robot schedule, as motivated by a BASF high-throughput laboratory, see
Fig. 1. We solve a large set of instances for this optimization problem with
varying difficulty using three commercially available solvers: (i) the D-Wave
quantum annealer, (ii) the quantum-inspired Fujitsu digital annealer and (iii)
the classical state-of-the-art solver Gurobi. To this end, we develop several
mathematical models: a QUBO model for the quantum and digital annealer and two
different MIP models for Gurobi, which we call time-indexed and sequence
model, respectively. Modeling the same problem in different, solver-specific
forms helps us to optimally assess the capabilities of each solver. In total,
we compare five different approaches (i. e., model and solver combinations as
sketched in Fig. 5): (i) TI-GU, (ii) SE-GU, (iii) QU-LBQM, (iv) QU-FDA and (v)
QU-FDAh. For our performance study, we separated all problem instances into
two groups. First, the _minor instances_ with problems less than $10\,000$
binary variables in the QUBO formulation and, second, the _major instances_
with problems with more than $10\,000$ and up to $22\,000$ variables. For
practical reasons, we only solve the minor instances with SE-GU, TI-GU, QU-
LBQM and QU-FDA, whereas the major instances are only solved with SE-GU and
QU-FDAh, respectively.

Our benchmark reveals insights both regarding the objective values of the
optimization problem (i. e., the sum of sample completion times) as well as
the end-to-end runtimes for the considered approaches. Regarding the objective
values, we observe for minor instances that SE-GU and TI-GU give similar
results, outperforming QU-FDA, which in turn outperforms QU-LBQM (cf. Fig.
6(a)). For major instances, SE-GU outperforms QU-FDAh (cf. Fig. 9(b)).
Regarding the runtime, we find that for smaller instances TI-GU takes the
highest time and SE-GU takes mostly the lowest. Between these two extremes,
QU-FDA and QU-LBQM take about the same amount of time (cf. Fig. 6(b)).
However, the runtime of SE-GU significantly increases with increasing instance
complexity. This same observation continues for the large instances, for which
the runtime of SE-GU is mostly larger than that of QU-FDAh (cf. Fig. 9(a)).

To get further insights into the relationship between objective value and
runtime, we also studied the relative runtime of Gurobi, that is the time that
Gurobi took to find an objective value that is at least as good as the final
result from another approach. For minor instances, we find that the relative
runtimes of SE-GU w.r.t. QU-LBQM and QU-FDA, respectively, are strictly lower
than the runtimes of QU-LBQM and QU-FDA, i. e., Gurobi found solutions of
comparable quality faster than the quantum and quantum-inspired approaches
(cf. Figs. 8(a) and 8(b)). This is not surprising since SE-GU tended to find
better objectives in shorter time. For major instances, the relative runtimes
of SE-GU w.r.t. QU-FDAh increase significantly with increasing instance
complexity and clearly exceed the runtime of QU-FDAh for the biggest instances
(cf. Fig. 10). Thus, QU-FDAh shows an advantage on some bigger instances.
Although the resulting objective values of QU-FDAh were not optimal, the
approach shows a clear advantage on some bigger instances when compared to SE-
GU on a similar time scale.

Summarized, no general advantage of the quantum and quantum-inspired solvers
was found. However, for certain instances the quantum-inspired hybrid usage of
the Fujitsu digital annealer turned out to be a very promising alternative to
Gurobi and was clearly superior to a hybrid usage of D-Wave’s quantum
annealer. Our study is not a conclusive result but rather an application-
oriented case study that provides a snapshot of the current technology and
leaves room for performance improvements on the modeling as well as the solver
side. For example, an improvement of the quantum annealer might be possible
with additional problem-specific fine-tuning of the annealing schedule or
other hardware-related parameters. Moreover, the recently released constrained
quadratic model (CQM) solver from D-Wave also promises to provide much better
performance compared to the solver used in this work. Especially in an agile
field such as quantum computing, a technology snapshot such as ours can hardly
provide any forecasts about future developments. Therefore, in order to
preserve an up-to-date assessment, further practical evaluations for real-
world use cases will be necessary. The methods and results from this project
can serve as a blueprint or at least point of reference for this kind of
ongoing research.

##  6 Acknowledgements

We would like to thank Behrang Shafei, Jens Meissner and Horst Weiss for their
invaluable input and support throughout the research process. Without their
ongoing contributions, the work would not have been accomplished. We would
also like to acknowledgement Martin Schuetz, for making us aware of a similar
work [49].

##  Appendix A Models

All models shown in this section use the following problem parameters:

  * •

machines: $M_{1}$ (water mixer), $M_{2}$ (sample shaker), $M_{3}$ (photo
booth),

  * •

samples: $j\in\\{1,2,3\dotsc,N\\}$,

  * •

machine processing times: $p_{j,m}$ for $j\in\\{1,\ldots,N\\}$,
$m\in\\{1,2,3\\}$,

  * •

time (slots): $t\in\\{0,\dots,T-1\\}$,

  * •

photos: $k\in\\{1,\ldots,K\\}$,

  * •

gap between photos $k$ and $k+1$ of sample $j\in\\{1,\dotsc,N\\}$: $g_{j,k}$
for $k\in\\{1,\dotsc,K-1\\}$.

###  A.1 QUBO Model

####  A.1.1 Variables

The QUBO model uses the variables

$x_{j,m,t}:=\begin{cases}1,&\text{if sample $j$ starts processing on machine
$M_{m}$ at time $t$},\\\ 0,&\text{otherwise}\end{cases}$ (8)

with $x_{j,m,t}\in\\{0,1\\}$ for all $j\in\\{1,\ldots,N\\}$, $m\in\\{1,2,3\\}$
and $t\in\\{1,\ldots,T-1\\}$. Here,

$T:=T(N,K,p,g):=\left(\sum_{j=1}^{N}p_{j,1}+p_{j,2}+K\cdot
p_{j,3}\right)+\left(\sum_{j=1}^{N}\sum_{k=1}^{K-1}g_{j,k}\right)+N\cdot(3+K)$
(9)

is the time horizon for an instance with parameters $N$, $K$,
$p:=\\{p_{j,m}\mid j\in\\{1,\ldots,N\\},\ m\in\\{1,2,3\\}\\}$,
$g:=\\{g_{j,k}\mid j\in\\{1,\ldots,N\\},k\in\\{1,\ldots,K-1\\}\\}$.

In total, there are $N\cdot 3\cdot T(N,K,p,g)$ binary variables. However, a
certain subset of these variables can by definition be fixed to $0$ to reduce
the number of effective optimization variables:

  * i)

The processing of any job on machines $M_{1}$, $M_{2}$ and $M_{3}$ can only
start after the preceding processes are finished. Thus, for all
$j\in\\{1,\dotsc,N\\}$ and $m\in\\{1,2,3\\}$ we set

$x_{j,m,t}=0\ \ \forall t:\ 0\leq t\leq
m-1+\sum_{m^{\prime}=1}^{m-1}p_{j,m^{\prime}}.$ (10)

  * ii)

Analogously the latest start for a process on machines $M_{1}$ and $M_{2}$ can
be bounded by $T-1$ minus the total time taken by the subsequent processes.
For all jobs $j=1,\ldots,N$ and $m=1,2$ we set

$x_{j,m,t}=0\ \ \forall t:\ (T-1)-\left(\sum_{k=1}^{K-1}g_{j,k}+K\cdot
p_{j,3}+\sum_{m^{\prime}=m}^{2}(1+p_{j,m^{\prime}})\right)\leq t\leq T-1.$
(11)

Summarized, $c(N,K,p,g)$ variables are required to attain a constant value of
zero to achieve feasibility with $c(N,K,p,g)$ being the number of variables
set to zero according to equations 10 and 11. Consequently, we arrive at a
QUBO with

$n:=N\cdot 3\cdot T(N,K,p,g)-c(N,K,p,g)$ (12)

binary optimization variables $x_{j,m,t}$.

####  A.1.2 Constraints

In the following, we write $x$ for the vector consisting of these $n$
optimization variables in an arbitrary but fixed order.

With this notation, the penalty terms for the QUBO read:

  * 1)

Each sample is processed precisely once by $M_{1}$ and $M_{2}$.

$P_{1}(x):=\sum_{j=1}^{N}\sum_{m=1}^{2}\left[\left(\sum_{t=1}^{T-1}x_{j,m,t}\right)-1\right]^{2}\;.$
(13)

  * 2)

All samples are processed $K$ times by the photo booth $M_{3}$.

$P_{2}(x):=\sum_{j=1}^{N}\left[\left(\sum_{t=1}^{T-1}x_{j,3,t}\right)-K\right]^{2}\;.$
(14)

  * 3)

No two processings can start at the same time (because processings start
immediately when a sample is brought to a machine).

$P_{3}(x):=\sum_{t=0}^{T-1}\left(\sum_{j=1}^{N}\sum_{m=1}^{3}x_{j,m,t}\right)\left[\left(\sum_{j=1}^{N}\sum_{m=1}^{3}x_{j,m,t}\right)-1\right]\;.$
(15)

  * 4)

A machine $M_{m}$ is occupied in the time frame between $t$ and $t+p_{j,m}$ in
the case that any sample $j$ has been put on it at time $t$.

$P_{4}(x):=\sum_{m=1}^{3}\sum_{j=1}^{N}\sum_{t=0}^{T-1}\Biggl{(}x_{j,m,t}\cdot\sum_{j\neq
j^{\prime}=1}^{N}\sum_{t^{\prime}=t+1}^{\min\\{t+p_{j,m},T-1\\}}x_{j^{\prime},m,t^{\prime}}\Biggr{)}$
(16)

  * 5)

Processing of a sample on consecutive machines must not overlap.

$P_{5}(x):=\sum_{m=1}^{2}\sum_{j=1}^{N}\sum_{t=0}^{T-1}\Biggl{(}x_{j,m+1,t}\cdot\sum_{t^{\prime}=\max\\{t-p_{j,m},0\\}}^{T-1}x_{j,m,t^{\prime}}\Biggr{)}$
(17)

  * 6)

After a sample is processed at machine $M_{2}$, it needs to be immediately
brought to the photo booth, and the subsequent photos need to be taken at
specific times later.

Let
$\gamma_{j,k}:=1+\sum_{k^{\prime}=1}^{k-1}g_{j,k^{\prime}}+p_{j,3}\cdot(k-1)$
for $k\in\\{1,\ldots,K\\}$ be the desired gap between the end of the shaking
process on $M_{2}$ and start of the $k$th photo on $M_{3}$. Set

$P_{6}(x):=\sum_{j=1}^{N}\sum_{k=1}^{K}\sum_{t=0}^{T-1-\gamma_{j,K}}x_{j,2,t}\cdot(1-x_{j,3,t+\gamma_{j,k}+p_{j,2}})\;.$
(18)

  * 7)

Certain pairs of tasks require a gap of one time unit between them for the
robot to travel between the machines and the rack empty-handed. We set
$P_{7}:=P_{7.1}+P_{7.2}+P_{7.3}$ where $P_{7.1}$, $P_{7.2}$ and $P_{7.3}$ are
defined as follows:

7.1) If a job $j$ is put on a machine $M_{m}$ at time $t$, then no job can be
put on any machine at time $t+1$.

$P_{7.1}(x):=\sum_{t=0}^{T-2}\Bigg{(}\sum_{j=1}^{N}\sum_{m=1}^{3}x_{j,m,t}\cdot\sum_{j^{\prime}=1}^{N}\sum_{m^{\prime}=1}^{3}x_{j^{\prime},m^{\prime},t+1}\Bigg{)}$
(19)

7.2) If a job $j$ is picked up from a machine $M_{m}$ at time $t$, then no
other job $j^{\prime}$ can be picked up from any other machine at times $t$
and $t+1$.

$P_{7.2}(x):=\sum_{t=0}^{T-2}\sum_{j=1}^{N}\Bigg{(}\sum_{m=1}^{3}x_{j,m,t-p_{j,m}}\cdot\sum_{j^{\prime}=1}^{N}\sum_{\begin{subarray}{c}m^{\prime}=1\\\
m^{\prime}\neq
m\end{subarray}}^{3}(x_{j^{\prime},m^{\prime},t-p_{j^{\prime},m^{\prime}}}+x_{j^{\prime},m^{\prime},t+1-p_{j^{\prime},m^{\prime}}})\Bigg{)}$
(20)

7.3) If a job $j$ is picked up from a machine $M_{m}$ at time $t$, then no
other job can be put on any machine at times $t$ and $t+1$.

$P_{7.3}(x):=\sum_{t=0}^{T-1}\sum_{m=1}^{3}\Bigg{(}\sum_{j=1}^{N}x_{j,m,t-p_{j,m}}\cdot\sum_{\begin{subarray}{c}j^{\prime}=1\\\
j^{\prime}\neq
j\end{subarray}}^{N}\sum_{m^{\prime}=1}^{3}(x_{j^{\prime},m^{\prime},t}+x_{j^{\prime},m^{\prime},t+1})\Bigg{)}$
(21)

####  A.1.3 Objective

The objective function computes as

$F(x):=\sum_{j=1}^{N}\sum_{t=1}^{T-1}\left(t+p_{j,2}+1+K\cdot
p_{j,3}+\sum_{k=1}^{K-1}g_{j,k}+1\right)\cdot x_{j,2,t}\;.$ (22)

####  A.1.4 Summary

Altogether, the _QUBO model_ is defined as:

$\displaystyle\min_{x}$
$\displaystyle\rho_{0}F(x)+\sum_{i=1}^{7}\rho_{i}P_{i}(x)$ (23) s.t.
$\displaystyle x\in\\{0,1\\}^{n}$

###  A.2 Sequence Model

####  A.2.1 Variables

We define the event set
$E\coloneqq\bigl{\\{}(j,i,a)\bigm{|}j\in\\{1,\dotsc,N\\},\
i\in\\{1,\dotsc,2+K\\},\ a\in\\{\text{0 (place)},\text{1 (pick
up)}\\}\bigr{\\}}$, where event

$(j,i,a)\text{ represents action }a\text{ for sample }j\ \begin{cases}\text{on
machine }i&\text{for }i\leq 2;\\\ \text{at the camera for photo }i-2&\text{for
}i\geq 3.\end{cases}$ (24)

The following variables are used in the sequence model:

  * •

$\tau_{e}\in\mathbb{R}_{\geq 0}$: Time of event $e$ for all $e\in E$.

  * •

$\delta_{m,j,j^{\prime}}=\begin{cases}1,&\text{if sample $j$ is processed by
machine\leavevmode\nobreak\ $m$ before sample\leavevmode\nobreak\
$j^{\prime}$}\;,\\\ 0,&\text{otherwise}\;,\end{cases}$
for all $m\in\\{1,2\\}$, $j,j^{\prime}\in\\{1,\ldots,N\\},\ j\neq j^{\prime}$.

  * •

$\varepsilon_{e,f}=\begin{cases}1,&\text{if event $e$ happens before event
$f$}\;,\\\ 0,&\text{otherwise}\;,\end{cases}$
for all $e=(j,i,a),f=(j^{\prime},i^{\prime},a^{\prime})\in E,\ j\neq
j^{\prime},\ i\neq i^{\prime}$, $\min\\{i,i^{\prime}\\}\leq 2$.

  * •

$\zeta_{k,j,k^{\prime},j^{\prime}}=\begin{cases}1,&\text{if the }k\text{th
photo of sample }j\text{ is taken before the }k^{\prime}\text{th photo of
sample }j^{\prime}\;,\\\ 0,&\text{otherwise}\;,\end{cases}$
for all $j,j^{\prime}\in\\{1,\ldots,N\\}$, for all
$k,k^{\prime}\in\\{1,\ldots,K\\}$ with $j\neq j^{\prime}$, $k\neq k^{\prime}$.

####  A.2.2 Constraints

To formulate the constraints we introduce the following notation:

  * •

For $j\in\\{1,\ldots,n\\}$ and $k\in\\{1,\ldots,K\\}$ we denote by
$\gamma_{j,k}:=1+\sum_{k^{\prime}=1}^{k-1}g_{j,k^{\prime}}+(k-1)\cdot p_{j,3}$
the gap between the end of the shaking process of sample $j$ and the beginning
of the process of the $k$th photo.

  * •

$\displaystyle\Gamma\coloneqq\sum_{j=1}^{N}\Biggl{[}\biggl{(}\sum_{m=1}^{3}p_{j,i}+1\biggr{)}+\gamma_{j,K}+1\Biggr{]}$
is a sufficiently large constant that is used for the “Big M” method [25].

Then we define the following list of constraints:

  1. 1)

No two samples are processed on neither the water mixer nor the sample shaker
at the same time.

$\delta_{m,j,j^{\prime}}+\delta_{m,j^{\prime},j}=1\qquad\forall
m\in\\{1,2\\},\ j,j^{\prime}\in\\{1,\ldots,N\\},\ j\neq j^{\prime}.$ ($C_{1}$)

  2. 2)

The robot cannot perform two events of different samples at different places
at the same time.

$\varepsilon_{e,f}+\varepsilon_{f,e}=1\qquad\forall
e=(j,i,a),f=(j^{\prime},i^{\prime},a^{\prime})\in E,\ j\neq j^{\prime},\ i\neq
i^{\prime}.$ ($C_{2}$)

  3. 3)

No two samples can be at the photo booth for different photos at the same
time.

$\zeta_{k,j,k^{\prime},j^{\prime}}+\zeta_{k^{\prime},j^{\prime},k,j}=1\qquad\forall
j,j^{\prime}\in\\{1,\ldots,N\\},\ k,k^{\prime}\in\\{1,\dotsc,K\\},\ j\neq
j^{\prime},\ k\neq k^{\prime}.$ ($C_{3}$)

  4. 4)

The sample must be picked up exactly after its processing time.

$\tau_{(j,i,1)}=\tau_{(j,i,0)}+p_{j,i}\qquad\forall j\in\\{1,\ldots,N\\},\
i\in\\{1,2\\}.$ ($C_{4}$)

  5. 5)

The transportation time of $1$ time unit must be waited for.

$\tau_{(j,i^{\prime},0)}\geq\tau_{(j,i,1)}+1\qquad\forall
j\in\\{1,\ldots,N\\},\ i,i^{\prime}\in\\{1,\dotsc,2+K\\},\ i^{\prime}>i.$
($C_{5}$)

  6. 6)

The time frames at water mixer and sample shaker must match the orders
specified by the $\delta$-variables. The robot must first bring a sample away
and can then get another one.

$\tau_{(j^{\prime},i,0)}\geq\tau_{(j,i,1)}+2-\Gamma\cdot\delta_{i,j^{\prime},j}\qquad\forall
i\in\\{1,2\\},\ j,j^{\prime}\in\\{1,\ldots,N\\},\ j\neq j^{\prime}.$ ($C_{6}$)

  7. 7)

The $k$th photos of different samples must not overlap.

$\tau_{(j^{\prime},2+k,0)}\geq\tau_{(j,2+k,1)}+2-\Gamma\cdot\delta_{2,j^{\prime},j}\qquad\forall
j,j^{\prime}\in\\{1,\ldots,N\\},\ j\neq j^{\prime},\ k\in\\{1,\dotsc,K\\}.$
($C_{7}$)

  8. 8)

Different photos of different samples must not overlap.

$\tau_{(j^{\prime},2+k^{\prime},0)}\geq\tau_{(j,2+k,1)}+2-\Gamma\cdot\zeta_{k^{\prime},j^{\prime},k,j}\\\
\forall j,j^{\prime}\in\\{1,\ldots,N\\},\ k,k^{\prime}\in\\{1,\dotsc,K\\},\
j\neq j^{\prime},\ k\neq k^{\prime}.$ (25)

  9. 9)

The times of the events of different samples at different machines must match
the order and the robot must have enough time in between.

$\displaystyle\tau_{f}$
$\displaystyle\geq\tau_{e}-\Gamma\cdot\varepsilon_{f,e}+2$
$\displaystyle\forall e=(j,i,1),f=(j^{\prime},i^{\prime},a^{\prime})\in E,\
j\neq j^{\prime},\ i\neq i^{\prime},\ \min\\{i,i^{\prime}\\}\leq 2.$
($C_{9.1}$) $\displaystyle\tau_{f}$
$\displaystyle\geq\tau_{e}-\Gamma\cdot\varepsilon_{f,e}+2$
$\displaystyle\forall e=(j,i,a),f=(j^{\prime},i^{\prime},0)\in E,\ j\neq
j^{\prime},\ i\neq i^{\prime},\ \min\\{i,i^{\prime}\\}\leq 2.$ ($C_{9.2}$)
$\displaystyle\tau_{f}$
$\displaystyle\geq\tau_{e}-\Gamma\cdot\varepsilon_{f,e}+1$
$\displaystyle\forall e=(j,i,0),f=(j^{\prime},i^{\prime},1)\in E,\ j\neq
j^{\prime},\ i\neq i^{\prime},\ \min\\{i,i^{\prime}\\}\leq 2.$ ($C_{9.3}$)

  10. 10)

The photo times must have the correct distance to the end of the shaking
process.

$\displaystyle\tau_{(j,2+k,0)}$ $\displaystyle=\tau_{(j,2,1)}+\gamma_{j,k}$
$\displaystyle\forall j\in\\{1,\ldots,N\\},\ k\in\\{1,\dotsc,K\\}.$
($C_{10.1}$) $\displaystyle\tau_{(j,2+k,1)}$
$\displaystyle=\tau_{(j,2,1)}+\gamma_{j,k}+p_{j,3}$ $\displaystyle\forall
j\in\\{1,\ldots,N\\},\ k\in\\{1,\dotsc,K\\}.$ ($C_{10.2}$)

####  A.2.3 Objective

The objective function computes as

$F(\tau,\delta,\varepsilon,\zeta):=\sum_{j=1}^{N}(\tau_{(j,2+K,1)}+2).$ (26)

Note that by definition of the model, it is a valid solution that a sample is
placed at time $0$ on machine $M_{1}$, which implies that the robot transports
the sample at time $-1$ from the rack to $M_{1}$. To overcome this we add an
extra $+1$ to each sample to shift the whole schedule to the right. This does
not affect the optimization process itself, but is rather a linear translation
of the whole schedule w.r.t. the objective function to get a schedule starting
at $0$ in compliance with the other models. Lastly the sample arrives at the
rack one time unit after the $K$’th photo has been shot, thus we add another
$+1$.

####  A.2.4 Summary

Altogether, the _sequence model_ is defined as:

$\displaystyle\min_{\tau,\delta,\varepsilon,\zeta}$ $\displaystyle
F(\tau,\delta,\varepsilon,\zeta)$ (27) s.t. $\displaystyle
C_{1},\ldots,C_{8},C_{9.1},C_{9.2},C_{9.3},C_{10.1},C_{10,2}$
$\displaystyle\tau_{e}\in\mathbb{R}_{\geq 0}\ \forall e\in E$
$\displaystyle\delta_{m,j,j^{\prime}}\in\\{0,1\\}\ \forall
m\in\\{1,2\\},j,j^{\prime}\in\\{1,\ldots,N\\},\ j\neq j^{\prime}$
$\displaystyle\varepsilon_{e,f}\in\\{0,1\\}\ \forall
e=(j,i,a),f=(j^{\prime},i^{\prime},a^{\prime})\in E,\ j\neq j^{\prime},\ i\neq
i^{\prime},\min\\{i,i^{\prime}\\}\leq 2$
$\displaystyle\zeta_{k,j,k^{\prime},j^{\prime}}\in\\{0,1\\}\ \forall
j,j^{\prime}\in\\{1,\ldots,N\\},\ \forall
k,k^{\prime}\in\\{1,\ldots,K\\}\text{ with }j\neq j^{\prime},k\neq k^{\prime}$

###  A.3 Time-Indexed Model

####  A.3.1 Variables

The time-indexed model makes use of the binary variables

$y_{j,r,t}:=\begin{cases}1,\text{ if sample $j$ is transported by the robot on
route $r$ at time $t$}\;,\\\ 0,\text{ otherwise}\;\end{cases}$ (28)

for all $j\in\\{1,\ldots,N\\}$, $r\in\\{1,\ldots,8\\}$ and
$t\in\\{0,\ldots,T-1\\}$, where the routes $r$ are given by Fig. 3. The time
horizon $T$ is defined as in equation 9. Additionally we have variables
$z_{j}\in\mathbb{R}_{\geq 0}$ for $j\in\\{1,\ldots,N\\}$ that are used to
define the objective function (see equation 30).

Setting $M_{0}$ to be the rack we define $r^{-}(i)$ to be the set of incoming
routes and $r^{+}(i)$ the set of outgoing routes of $M_{i}$ for
$i\in\\{0,1,2,3\\}$, i. e. $r^{-}(0)=\\{2,5,8\\}$, $r^{+}(0)=\\{1,4,7\\}$,
$r^{-}(1)=\\{1\\}$ and so on. The set of all routes is $R:=\\{1,\ldots,8\\}$.

####  A.3.2 Constraints

We define the following constraints:

  1. 1)

The transport robot can carry at most one sample at a time.

$\sum_{j=1}^{N}\sum_{r=1}^{8}y_{j,r,t}\leq 1\qquad\forall
t\in\\{0,\dotsc,T-1\\}.$ ($C_{1}$)

  2. 2)

The transport robot cannot teleport from one machine to another. For every
time $t$ and every machine $i$ the robot cannot serve some route to $i$ during
$(t,t+1)$ and a route starting from another machine during $(t+1,t+2)$.

$\sum_{j=1}^{N}\Biggl{[}\biggl{(}\,\sum_{r\in
r^{-}(i)}y_{j,r,t}\biggr{)}+\sum_{r\in R\setminus
r^{+}(i)}y_{j,r,t+1}\Biggr{]}\leq 1\qquad\forall i\in\\{0,\ldots,3\\},\
t\in\\{0,\dotsc,T-2\\}.$ ($C_{2}$)

  3. 3)

Every machine may be occupied by at most one sample.

$\sum_{j=1}^{N}\sum_{t^{\prime}=0}^{t}\Biggl{[}\bigg{(}\sum_{r\in
r^{-}(i)}y_{j,r,t^{\prime}}\bigg{)}-\sum_{r\in
r^{+}(i)}y_{j,r,t^{\prime}}\Biggr{]}\leq 1\qquad\forall i\in\\{1,2,3\\},\
t\in\\{0,\dotsc,T-1\\}.$ ($C_{3}$)

  4. 4)

Every sample is transported to the water mixer $M_{1}$ exactly once.

$\sum_{t=0}^{T-1}y_{j,1,t}=1\qquad\forall j\in\\{1,\ldots,N\\}.$ ($C_{4}$)

  5. 5)

Samples can only be picked up from a machine, when they have been brought to
it and stayed there for the appropriate time.

$\sum_{t^{\prime}=0}^{t}\sum_{r\in
r^{+}(i)}y_{j,r,t^{\prime}}=\sum_{t^{\prime}=0}^{t-p_{i,j}-1}\sum_{r\in
r^{-}(i)}y_{j,r,t^{\prime}}\qquad\forall j\in\\{1,\ldots,N\\},\
i\in\\{1,2,3\\},\ t\in\\{0,\dotsc,T-1\\}.$ ($C_{5}$)

  6. 6)

Routes between rack and the first two machines must be driven along in the
correct order.

$\sum_{t^{\prime}=0}^{t}y_{j,4,t^{\prime}}\leq\sum_{t^{\prime}=0}^{t-1}y_{j,2,t^{\prime}}\qquad\forall
j\in\\{1,\ldots,N\\},\ t\in\\{0,\dotsc,T-1\\}.$ ($C_{6}$)

  7. 7)

At no time has a sample been brought from the rack to the camera more often
than from the camera or the shaker to the rack.

$\sum_{t^{\prime}=0}^{t}y_{j,7,t^{\prime}}\leq\sum_{t^{\prime}=0}^{t-1}(y_{j,5,t^{\prime}}+y_{j,8,t^{\prime}})\qquad\forall
j\in\\{1,\ldots,N\\},\ t\in\\{0,\dotsc,T-1\\}.$ ($C_{7}$)

  8. 8)

A sample is picked up from the camera at time $t$ if and only if it has been
placed there at time $t-p_{j,3}$, i. e., $y_{j,8,t}=0$ for $t\leq p_{j,3}$ and
$y_{j,8,t}=y_{j,6,t-p_{j,3}-1}+y_{j,7,t-p_{j,3}-1}$ for $t>p_{j,3}$. These
constraints are equivalent to

$\sum_{t^{\prime}=0}^{t}y_{j,8,t^{\prime}}=\sum_{t^{\prime}=0}^{t-p_{j,3}-1}(y_{j,6,t^{\prime}}+x_{j,7,t^{\prime}})\qquad\forall
j\in\\{1,\ldots,N\\},\ t\in\\{0,\dotsc,T-1\\}.$ ($C_{8}$)

  9. 9)

Every sample has to be picked up from the shaker early enough for the last
photo.

$\sum_{t=0}^{T-1-p_{j,3}-\gamma_{j,K}}(y_{j,5,t}+y_{j,6,t})=1\qquad\forall
j\in\\{1,\ldots,N\\}.$ ($C_{9}$)

  10. 10)

Every sample has to be picked up from the camera the same number of times as
it is placed there.

$\sum_{t=0}^{T-1}(y_{j,6,t}+y_{j,7,t}-y_{j,8,t})=0\qquad\forall
j\in\\{1,\ldots,N\\}.$ ($C_{10}$)

  11. 11)

Samples must be at the camera at photo dates.

$y_{j,5,t}+y_{j,6,t}\leq\sum_{t^{\prime}=0}^{t+\gamma_{j,k}-1}(y_{j,6,t^{\prime}}+y_{j,7,t^{\prime}})-\sum_{t^{\prime}=0}^{t+\gamma_{j,k}+p_{j,3}-1}y_{j,8,t^{\prime}}\qquad\forall
j\in\\{1,\ldots,N\\},\\\ \ k\in\\{1,\ldots,K\\},\
t\in\\{0,\dotsc,T-\gamma_{j,k}-p_{j,3}\\}.$ (29)

  12. 12)

Lastly we define a constraint that bounds the objective variables $z_{j}$
below by the arrival time of sample $j$ at the rack after the schedule has
finished.

$(t+1)\cdot y_{j,8,t}\leq z_{j}\qquad\forall j\in\\{1,\ldots,N\\},\
t\in\\{0,\ldots,T-1\\},$ ($C_{F}$)

####  A.3.3 Objective

Let $y,z$ be vectors achieved by writing the corresponding variables into a
vector in some arbitrary but fixed order. Using the constraint $C_{F}$,
equation $C_{F}$, the objective function simply computes as

$F(y,z):=\sum_{j=1}^{N}z_{j}.$ (30)

####  A.3.4 Summary

Altogether, the _time-indexed model_ is defined as:

$\displaystyle\min_{y,z}$ $\displaystyle F(y,z)$ (31) s.t. $\displaystyle
C_{1},\ldots,C_{11},C_{F}$ $\displaystyle y_{j,r,t}\in\\{0,1\\}\ \forall
j\in\\{1,\ldots,N\\},r\in\\{1,\ldots,8\\},t\in\\{0,\ldots,T-1\\}$
$\displaystyle z_{j}\in\mathbb{R}_{\geq 0}\ \forall j\in\\{1,\ldots,N\\}$

##  Appendix B Benchmark Library

It is our goal to define a representative benchmark library with meaningful
instances. For this purpose, we specify a set of relevant values for each
parameter as listed in Table 2. We consider all $120$ parameter combinations
$N\in D_{N}$, $K\in D_{K}$ and $p_{3}\in D_{p_{3}}$. For each combination, the
remaining parameters are then drawn uniformly from the corresponding set of
allowed values $10$ times, which leads to $1200$ instances in total.

Table 2: Parameter configurations for our benchmark instances. Parameter | Choice | Domain
---|---|---
$N$ | all | $D_{N}:=\\{1,2,\ldots,10,15,20,25,50,100\\}$
$K$ | all | $D_{K}:=\\{1,2,3,4\\}$
$p_{3}$ | all | $D_{p_{3}}:=\\{1,3\\}$
$p_{j,1}$ | random | $D_{p_{1}}:=\\{4,5,6,7,8\\}$
$p_{j,2}$ | random | $D_{p_{2}}:=\\{1,2,3,4\\}$
$g_{j,1}$ | random | $D_{g_{1}}:=\\{4,5\\}$
$g_{j,2}$ | random | $D_{g_{2}}:=\\{8,\ldots,12\\}$
$g_{j,3}$ | random | $D_{g_{3}}:=\\{16,\ldots,24\\}$

However, we filter these instances such that only a relevant subset remains.
Specifically, we only keep the instances with at least $2071$ and at most
$22\,692$ binary variables in the corresponding QUBO formulation from Section
2.2.2. This excludes the very easy and the very hard cases. A complete list
with all computed instances, their parameters and the benchmark results is
appended as supplementary information.

We collect groups of instances $(N,K)$ that have the same number of samples
$N$ and photos $K$. In Table 3, we list how many instances each group
contains.

Table 3: Size informations about the instance groups: Number of instances as well as minimum and maximum number of variables for the the corresponding QUBO formulation. $(N,K)$ | Num. instances | Min num. vars | Max num. vars
---|---|---|---
$(4,4)$ | $19$ | $2096$ | $2488$
$(5,3)$ | $15$ | $2071$ | $2494$
$(5,4)$ | $20$ | $3384$ | $4118$
$(6,3)$ | $20$ | $2849$ | $3834$
$(6,4)$ | $20$ | $4935$ | $6109$
$(7,3)$ | $20$ | $4022$ | $5146$
$(7,4)$ | $16$ | $6787$ | $8080$
$(8,3)$ | $20$ | $5215$ | $6880$
$(9,3)$ | $11$ | $6750$ | $8042$
$(8,4)$ | $4$ | $10\,822$ | $11\,058$
$(9,4)$ | $20$ | $11\,590$ | $14\,106$
$(10,3)$ | $1$ | $11\,026$ | $11\,026$
$(10,4)$ | $20$ | $14\,307$ | $17\,580$
$(15,2)$ | $20$ | $12\,270$ | $16\,129$
$(15,3)$ | $10$ | $19\,775$ | $20\,581$
$(20,1)$ | $20$ | $14\,541$ | $18\,765$
$(20,2)$ | $4$ | $22\,064$ | $22\,692$

##  Appendix C Fujitsu Solver Settings

Both FDA and FDAh have various parameters that can be used to tune the solver
runs. Table 4 and Table 5 contain the values for each parameter of both
solvers that was used in the benchmark.

Table 4: FDA parameters. Here, $k_{\mathrm{iter}}$ denotes the number of iterations for QU-FDA as declared in the column “Minor instance limit” of Table 1.  Parameter | Value
---|---
“access_profile_file” | “AUTO”
“processor” | “DAv2”
“connection_mode” | “async”
“DAv2_optimization_method” | “annealing”
“DAv2_number_iterations” | $k_{\mathrm{iter}}$
“DAv2_number_runs” | $32$
“DAv2_temperature_start” | $10\,000$
“DAv2_temperature_end” | $10$
“DAv2_temperature_mode” | $0$
“DAv2_temperature_interval” | $k_{\mathrm{iter}}/$1\,000\,000$$
“DAv2_offset_increase_rate” | $100$
“DAv2_solution_mode” | “COMPLETE”
“DAv2_bit_precision” | $16$
Table 5: FDAh parameters. Here, $k_{\mathrm{time}}$ denotes the time limit for QU-FDAh as declared in the column “Major instance limit” of Table 1. Parameter | Value
---|---
“access_profile_file” | “AUTO”
“processor” | “DAv3”
“connection_mode” | “async”
“DAv3_time_limit_sec” | $k_{\mathrm{time}}$
“DAv3__use__target_energy” | “False”
“DAv3_target_energy” | $0$
“DAv3_num_solution” | $16$
“DAv3_num_group” | $4$
“DAv3_num_output_solution” | $5$
“DAv3_gs_num_iteration_factor” | $5$
“DAv3_gs_num_iteration_cl” | $200$
“DAv3_gs_penalty_auto_mode” | $1$
“DAv3_gs_penalty_coef” | $1$
“DAv3_gs_penalty_inc_rate” | $150$
“DAv3_gs_max_penalty_coef” | $0$

##  Appendix D Pre-study

Before we started with the actual benchmarks for this paper, we performed a
small pre-study in which we exemplarily tested our solver candidates on a few
instances to assess their overall behavior. This allowed us to exclude three
potential solver candidates beforehand. In this appendix section, we briefly
summarize our qualitative observations for these excluded solvers:

  * •

First, it turned out that QAOA on the IBM gate-based quantum computers [29] is
very limited in problem size (up to $\sim 10^{2}$ variables) and, in addition,
has an end-to-end runtime that is magnitudes longer than that of all other
solvers due to the queuing system of the cloud service which can lead to major
delays.

  * •

Second, Quera’s Rydberg quantum computer _Aquila_ [60] does not support local
detuning, which is necessary for the QUBO encoding from Ref.[42], and we did
not find an alternative MIP representation of our use case.

  * •

Third, TSB had a very bad performance on all test instances and did not yield
feasible solutions.

## References

  * [1] A. Agnetis “Scheduling no-wait robotic cells with two and three machines” In _European Journal of Operational Research_ 123.2, 2000, pp. 303–314 DOI: 10.1016/S0377-2217(99)00258-1
  * [2] A. Agnetis and D. Pacciarelli “Part sequencing in three-machine no-wait robotic cells” In _Operations Research Letters_ 27.4, 2000, pp. 185–192 DOI: 10.1016/S0167-6377(00)00046-8
  * [3] Tameem Albash and Daniel A. Lidar “Demonstration of a Scaling Advantage for a Quantum Annealer over Simulated Annealing” In _Physical Review X_ 8.3 American Physical Society, 2018, pp. 031016 DOI: 10.1103/PhysRevX.8.031016
  * [4] Yuri Alexeev et al. “Quantum Computer Systems for Scientific Discovery” In _PRX Quantum_ 2.1 American Physical Society, 2021, pp. 017001 DOI: 10.1103/PRXQuantum.2.017001
  * [5] Ali Allahverdi “A survey of scheduling problems with no-wait in process” In _European Journal of Operational Research_ 255.3, 2016, pp. 665–686 DOI: 10.1016/j.ejor.2016.05.036
  * [6] Maliheh Aramon et al. “Physics-Inspired Optimization for Quadratic Unconstrained Problems Using a Digital Annealer” In _Frontiers in Physics_ 7 Frontiers Media SA, 2019 DOI: 10.3389/fphy.2019.00048
  * [7] Abhishek Awasthi, Jöerg Läessig, Jens Leuschner and Thomas Weise “GPGPU-Based Parallel Algorithms for Scheduling Against Due Date” In _2016 IEEE International Parallel and Distributed Processing Symposium Workshops_ , 2016, pp. 766–775 DOI: 10.1109/IPDPSW.2016.66
  * [8] Abhishek Awasthi et al. “Quantum Computing Techniques for Multi-Knapsack Problems” arXiv, 2023 DOI: 10.48550/ARXIV.2301.05750
  * [9] Kostas Blekos et al. “A Review on Quantum Approximate Optimization Algorithm and its Variants”, 2023 DOI: 10.48550/ARXIV.2306.09198
  * [10] Edward H. Bowman “The Schedule-Sequencing Problem” In _Operations Research_ 7.5, 1959, pp. 621–624 DOI: 10.1287/opre.7.5.621
  * [11] Peter Brucker “Scheduling Algorithms” BerlinHeidelberg: Springer-Verlag, 2007 DOI: 10.1007/978-3-540-69516-5
  * [12] Peter Brucker, Edmund K. Burke and Sven Groenemeyer “A mixed integer programming model for the cyclic job-shop problem with transportation” In _Discrete Applied Mathematics_ 160.13, 2012, pp. 1924–1935 DOI: 10.1016/j.dam.2012.04.001
  * [13] Yudong Cao et al. “Quantum chemistry in the age of quantum computing” In _Chemical reviews_ 119.19 ACS Publications, 2019, pp. 10856–10915 DOI: 10.1021/acs.chemrev.8b00803
  * [14] Costantino Carugno, Maurizio Ferrari Dacrema and Paolo Cremonesi “Evaluating the job shop scheduling problem on a D-wave quantum annealer” In _Scientific Reports_ 12.1 Nature Publishing Group UK London, 2022, pp. 6539 DOI: 10.1038/s41598-022-10169-0
  * [15] V. Černý “Thermodynamical Approach to the Traveling Salesman Problem: An Efficient Simulation Algorithm” In _Journal of Optimization Theory and Applications_ 45.1, 1985, pp. 41–51 DOI: 10.1007/BF00940812
  * [16] I. Chakroun, N. Melab, M. Mezmaz and D. Tuyttens “Combining multi-core and GPU computing for solving combinatorial optimization problems” In _Journal of Parallel and Distributed Computing_ 73.12, 2013, pp. 1563–1577 DOI: 10.1016/j.jpdc.2013.07.023
  * [17] D-Wave Systems Inc. “D-Wave Hybrid Solver Service: An Overview” Last accessed 2022-11-08, https://www.dwavesys.com/media/4bnpi53x/14-1039a-b_d-wave_hybrid_solver_service_an_overview.pdf, 2020
  * [18] Milind W Dawande, H Neil Geismar, Suresh P Sethi and Chelliah Sriskandarajah “Throughput optimization in robotic cells”, ISOR 101 Springer Science & Business Media, 2007 DOI: 10.1007/0-387-70988-6
  * [19] David Deutsch and Richard Jozsa “Rapid solution of problems by quantum computation” In _Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences_ 439.1907 The Royal Society London, 1992, pp. 553–558 DOI: 10.1098/rspa.1992.0167
  * [20] S. Ebadi et al. “Quantum optimization of maximum independent set using Rydberg atom arrays” In _Science_ 376.6598 American Association for the Advancement of Science (AAAS), 2022, pp. 1209–1215 DOI: 10.1126/science.abo6587
  * [21] Edward Farhi, Jeffrey Goldstone and Sam Gutmann “A Quantum Approximate Optimization Algorithm”, 2014 DOI: 10.48550/ARXIV.1411.4028
  * [22] Jianguang Feng and Ada Che “Robotic cell cyclic scheduling considering cell layout” In _Proceedings of the 32nd Chinese Control Conference_ , 2013, pp. 2622–2626 URL: https://ieeexplore.ieee.org/abstract/document/6639868
  * [23] Marc Geitz et al. “Solving The Extended Job Shop Scheduling Problem With AGVs – Classical And Quantum Approaches” In _Integration of Constraint Programming, Artificial Intelligence, and Operations Research: 19th International Conference, (CPAIOR 2022), Proceedings_ , LNCS 13292 Los Angeles, CA, USA: Springer-Verlag, 2022, pp. 120–137 DOI: 10.1007/978-3-031-08011-1˙10
  * [24] Fred Glover, Gary Kochenberger and Yu Du “A Tutorial on Formulating and Using QUBO Models” arXiv, 2018 DOI: 10.48550/ARXIV.1811.11538
  * [25] I. Griva, S.G. Nash and A. Sofer “Linear and Nonlinear Optimization: Second Edition”, Other Titles in Applied Mathematics Society for IndustrialApplied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104), 2009
  * [26] Lov K Grover “A fast quantum mechanical algorithm for database search” In _Proceedings of the twenty-eighth annual ACM symposium on Theory of computing_ , 1996, pp. 212–219 DOI: 10.1145/237814.237866
  * [27] Gurobi Optimization, LLC “Gurobi Optimizer Reference Manual”, 2023 URL: https://www.gurobi.com
  * [28] Nicholas G. Hall and Chelliah Sriskandarajah “A Survey of Machine Scheduling Problems with Blocking and No-Wait in Process” In _Operations Research_ 44.3, 1996, pp. 510–525 DOI: 10.1287/opre.44.3.510
  * [29] IBM Quantum “IBM Quantum services” Last accessed 2023-01-31, https://quantum-computing.ibm.com, 2023
  * [30] Kazuki Ikeda, Yuma Nakamura and Travis S. Humble “Application of Quantum Annealing to Nurse Scheduling Problem” In _Scientific Reports_ 9, 2019 DOI: 10.1038/s41598-019-49172-3
  * [31] Caixia Jing, Wanzhen Huang and Guochun Tang “Minimizing total completion time for re-entrant flow shop scheduling problems” In _Theoretical Computer Science_ 412.48, 2011, pp. 6712–6719 DOI: 10.1016/j.tcs.2011.08.030
  * [32] S. Kirkpatrick, C.. Gelatt and M.. Vecchi “Optimization by Simulated Annealing” In _Science_ 220.4598, New Series American Association for the Advancement of Science, 1983, pp. 671–680 DOI: 10.1126/science.220.4598.671
  * [33] Gary Kochenberger et al. “The unconstrained binary quadratic programming problem: a survey” In _Journal of Combinatorial Optimization_ 28, 2014, pp. 58–81 DOI: 10.1007/s10878-014-9734-0
  * [34] Wen-Yang Ku and J. Beck “Mixed Integer Programming models for job shop scheduling: A computational analysis” In _Computers & Operations Research_ 73, 2016, pp. 165–173 DOI: 10.1016/j.cor.2016.04.006
  * [35] Yangyang Li et al. “Quantum Optimization and Quantum Learning: A Survey” In _IEEE Access_ 8, 2020, pp. 23568–23593 DOI: 10.1109/ACCESS.2020.2970105
  * [36] Shi Qiang Liu and Erhan Kozan “A hybrid metaheuristic algorithm to optimise a real-world robotic cell” In _Computers & Operations Research_ 84, 2017, pp. 188–194 DOI: 10.1016/j.cor.2016.09.011
  * [37] Elisabeth Lobe and Annette Lutz “Minor Embedding in Broken Chimera and Pegasus Graphs is NP-complete”, 2022 DOI: 10.48550/ARXIV.2110.08325
  * [38] Andrew Lucas “Ising formulations of many NP problems” In _Frontiers in Physics_ 2 Frontiers Media SA, 2014 DOI: 10.3389/fphy.2014.00005
  * [39] Alan S. Manne “On the Job-Shop Scheduling Problem” In _Operations Research_ 8.2, 1960, pp. 219–223 DOI: 10.1287/opre.8.2.219
  * [40] Yuta Mizuno and Tamiki Komatsuzaki “Finding Optimal Pathways in Chemical Reaction Networks Using Ising Machines”, 2023 DOI: 10.48550/arXiv.2308.04544
  * [41] Hiroshi Nakayama, Junpei Koyama, Noboru Yoneoka and Toshiyuki Miyazawa “Description: Third Generation Digital Annealer Technology” Last accessed 2023-06-13, https://www.fujitsu.com/global/documents/about/research/techintro/3rd-g-da_en.pdf, 2021
  * [42] Minh-Thi Nguyen et al. “Quantum Optimization with Arbitrary Connectivity Using Rydberg Atom Arrays” In _PRX Quantum_ 4.1 American Physical Society, 2023, pp. 010316 DOI: 10.1103/PRXQuantum.4.010316
  * [43] Hiroki Oshiyama and Masayuki Ohzeki “Benchmark of quantum-inspired heuristic solvers for quadratic unconstrained binary optimization” In _Scientific Reports_ 12.1 Springer ScienceBusiness Media LLC, 2022 DOI: 10.1038/s41598-022-06070-5
  * [44] L.. Phillips and P.. Unger “Mathematical Programming Solution of a Hoist Scheduling Program” In _AIIE Transactions_ 8.2 Taylor & Francis, 1976, pp. 219–225 DOI: 10.1080/05695557608975070
  * [45] Michael L. Pinedo “Scheduling” Springer International, 2016 DOI: 10.1007/978-3-319-26580-3
  * [46] A. B. Pritsker, Lawrence J. Waiters and Philip M. Wolfe “Multiproject Scheduling with Limited Resources: A Zero-One Programming Approach” In _Management Science_ 16.1, 1969, pp. 93–108 DOI: 10.1287/mnsc.16.1.93
  * [47] Jack Raymond et al. “Hybrid quantum annealing for larger-than-QPU lattice-structured problems” In _ACM Transactions on Quantum Computing_ 4.3 Association for Computing Machinery (ACM), 2023, pp. 1–30 DOI: 10.1145/3579368
  * [48] Hans Röck “Some New Results in Flow Shop Scheduling” In _Zeitschrift für Operations Research_ 28.1, 1984, pp. 1–16 DOI: 10.1007/BF01919082
  * [49] Martin J.A. Schuetz et al. “Optimization of Robot-Trajectory Planning with Nature-Inspired and Hybrid Quantum Algorithms” In _Phys. Rev. Appl._ 18.5 American Physical Society, 2022, pp. 054045 DOI: 10.1103/PhysRevApplied.18.054045
  * [50] Maria Schuld, Ilya Sinayskiy and Francesco Petruccione “An introduction to quantum machine learning” In _Contemporary Physics_ 56.2 Taylor & Francis, 2015, pp. 172–185 DOI: 10.48550/arXiv.1409.3097
  * [51] Dvir Shabtay and Kfir Arviv “Optimal robot scheduling to minimize the makespan in a three-machine flow-shop environment with job-independent processing times” In _Applied Mathematical Modelling_ 40.5, 2016, pp. 4231–4247 DOI: 10.1016/j.apm.2015.11.021
  * [52] Peter W Shor “Algorithms for quantum computation: discrete logarithms and factoring” In _Proceedings 35th annual symposium on foundations of computer science_ , 1994, pp. 124–134 IEEE DOI: 10.1109/SFCS.1994.365700
  * [53] George Steiner and Zhihui Xue “Scheduling in Reentrant Robotic Cells: Algorithms and Complexity” In _Journal of Scheduling_ 8.1, 2005, pp. 25–48 DOI: 10.1007/s10951-005-5314-6
  * [54] Helman I. Stern and Gad Vitner “Scheduling Parts in a Combined Production-transportation Work Cell” In _Journal of the Operational Research Society_ 41.7 Taylor & Francis, 1990, pp. 625–632 DOI: 10.1057/jors.1990.86
  * [55] Michael Streif et al. “Beating classical heuristics for the binary paint shop problem with the quantum approximate optimization algorithm” In _Physical Review A_ 104.1 American Physical Society (APS), 2021 DOI: 10.1103/physreva.104.012403
  * [56] Kosuke Tatsumura, Alexander R. Dixon and Hayato Goto “FPGA-Based Simulated Bifurcation Machine” In _2019 29th International Conference on Field Programmable Logic and Applications (FPL)_ , 2019, pp. 59–66 DOI: 10.1109/FPL.2019.00019
  * [57] Dawid Tomasiewicz, Maciej Pawlik, Maciej Malawski and Katarzyna Rycerz “Foundations for Workflow Application Scheduling on D-Wave System” In _Computational Science–ICCS 2020_ Springer, 2020, pp. 516–530 DOI: 10.1007/978-3-030-50433-5˙40
  * [58] Harvey M. Wagner “An integer linear-programming model for machine scheduling” In _Naval Research Logistics Quarterly_ 6.2, 1959, pp. 131–140 DOI: https://doi.org/10.1002/nav.3800060205
  * [59] B.. Welch “THE GENERALIZATION OF ‘Student’s’ PROBLEM WHEN SEVERAL DIFFERENT POPULATION VARLANCES ARE INVOLVED” In _Biometrika_ 34.1–2, 1947, pp. 28–35 DOI: 10.1093/biomet/34.1-2.28
  * [60] Jonathan Wurtz et al. “Aquila: QuEra’s 256-qubit neutral-atom quantum computer”, 2023 DOI: 10.48550/arXiv.2306.11727
  * [61] Sheir Yarkoni et al. “Multi-car paint shop optimization with quantum annealing” In _2021 IEEE International Conference on Quantum Computing and Engineering (QCE)_ , 2021, pp. 35–41 DOI: 10.1109/QCE52317.2021.00019
  * [62] Stefanie Zbinden, Andreas Bärtschi, Hristo Djidjev and Stephan Eidenbenz “Embedding Algorithms for Quantum Annealers with Chimera and Pegasus Connection Topologies” In _High Performance Computing_ Cham: Springer International Publishing, 2020, pp. 187–206 DOI: 10.1007/978-3-030-50743-5˙10

  *[LBQM]: D-Wave’s hybrid Leap framework
  *[QUBO]: quadratic unconstrained binary optimization
  *[FDA]: Fujitsu’s digital annealer
  *[FDAh]: Fujitsu’s digital annealer hybrid framework
  *[TI-GU]: Gurobi with the time-indexed model
  *[QU-LBQM]:   () with the  model
  *[SE-GU]: Gurobi with the sequence model
  *[QU-FDA]:   () with the  model
  *[QU-FDAh]:   () with the  model
  *[TRSP]: transport robot scheduling problem
  *[MIP]: mixed integer program
  *[GPGPU]: general purpose computation on graphics processing unit
  *[JSSP]: job shop scheduling problem
  *[AGV]: automated guided vehicles
  *[QAOA]: quantum approximate optimization algorithm
  *[TSB]: Toshiba’s Simulated Bifurcation Machine
  *[MIS]: maximum independent set
  *[QUBOs]: quadratic unconstrained binary optimization
  *[TRSPs]: transport robot scheduling problem
  *[LP]: linear programming
  *[]:
  *[API]: application programming interface
  *[CQM]: constrained quadratic model