Introduction

Quantum annealing [1] has been applied to many combinatorial optimization problems, including NP-hard problems which have numerous applications in logistics, scheduling, and finance. The First In - First Out (FIFO) Stack-Up problem is shown to be NP-hard in [2], therefore it is very unlikely that this problem can be solved in polynomial time. Quantum annealing (QA) is an approach to solve NP-hard problems, formulated as optimization problems, approximately. In general quantum annealing is a hardware based approach for solving optimization problems. Some examples of practical problems analyzed with quantum annealers are in the areas of traffic flow optimization [3], vehicle routing [4, 5], protein folding [6], job scheduling [7], traffic signal optimization [8], air traffic management [9], training of deep neural networks [10], item listing in the web [11], and ionic diffusion [12]. The aim of our work is to investigate whether the FIFO Stack-Up Problem is suitable for quantum annealers. According to [13] it is difficult to know a priori if any given real-world problem will benefit from QA. Its authors claim that only some problems are good candidates for quantum annealing. A summary of error sources can be found e.g. in [14]. QA often gives good results when the number and connectivity of the qubits required for the problem do not scale strongly with the size of the problem. This is not the case for the FIFO Stack-Up problem.

Central to a quantum annealer are qubits. Qubits are different to traditional bits in that they can take states in-between 0 and 1. The idea of quantum annealing can be represented by the following Hamiltonian

$$\begin{aligned} H = -A \sum _i {\hat{\sigma}}^{(i)}_x + B \left(\sum _{i} h_{i} {{\hat{\sigma}}}^{(i)}_z + \sum _{i>j} J_{i,j} {\hat{\sigma}}^{(i)}_z {\hat{\sigma}}^{(j)}_z\right) \end{aligned}$$
(1)

where \({\hat{\sigma }}^{(i)}_{x}\) and \({\hat{\sigma }}^{(i)}_{z}\) are Pauli matrices operating on a qubit \(q_i\), and \(h_i\) and \(J_{i,j}\) are the qubit biases and coupling strengths of the couplers of two qubits. The process starts in the ground state of the first term of the Hamiltonian which is a superposition state of 0 and 1. A and B are variables which are changed by decreasing A and increasing B in such a way that the final state is a ground state of the problem Hamiltonian, which is the second term of Eq. (1). After this process is complete, measuring the qubit states yields a proposed solution to the problem.Footnote 1 The problems are formulated as finding a minimum of an objective function, which describes a binary quadratic model (BQM).

According to [15] the minima of objective functions can be formulated in two forms, QUBO-formulation

$$\begin{aligned} \min _{\textbf{X}} E({\textbf{X}}) = \min _{\textbf{X}} \left( \sum _{i=1}^n Q_{i,i} x_i + \sum _{i<j} Q_{i,j} x_i x_j \right)&\quad&x_i \in \{0, 1\} \end{aligned}$$
(2)

and Ising-formulation

$$\begin{aligned} \min _{\textbf{S}} E({\textbf{S}}) = \min _{\textbf{S}} \left( \sum _{i=1}^n h_i s_i + \sum _{i<j} J_{i,j} s_i s_j \right)&\quad&s_i \in \{-1, +1\} \end{aligned}$$
(3)

The two formulations are equivalent via bijective relations \(h_i = Q_{i,i}/2 + \sum _{i<j} Q_{i,j}/4 + \sum _{j<i} Q_{j,i}/4\), and \(J_{i,j} = {Q_{i,j}}/4\).

BQMs are naturally represented by a graph \(P=(V_P, E_P)\), where in \(V_P\) each variable \(s_i\) is represented as a node \(s_i\) with weight \(h_i\), and in \(E_P\) we have for every pair \(i < j\) with nonzero \(J_{i,j}\) an edge \(e = \{s_i, s_j\}\) with edge weight \(J_{i,j}\). QUBOs are a class of NP-hard optimization problems; as we can use QUBOs to optimize the number of satisfied constraints in an instance of 0/1 Integer Programming – one of Karp’s original 21 NP-complete problems [16].

Consequently, many hard to compute problems can be converted into QUBO-models in polynomial time. Therefore, a way to approximate the lowest energy state of a QUBO-model quickly and precisely would provide a way to approximate the solutions to those problems as well.

A logical variable of a QUBO model can be represented by one or more qubits. The process of transforming the logical variables is called embedding. It is important to note that due to physical limitations, couplers only exist between some qubits. Therefore, the structure of the physical model is limited by the device’s topology. For this report, the device used was D-Wave’s Advantage System, which uses the Pegasus Topology, for details see [17, 18]. Due to these physical constraints, the Pegasus Topology has bounded degree, i.e. a limited number of neighbor vertices, since each qubit can have at most a constant number of couplers and all neighbor qubits have to be in bounded distance. Since the graph is far from fully connected, most QUBO-models do not map directly onto the graph. This means that solving a problem on a quantum annealer generally involves 2 steps: First the problem is reformulated into a QUBO-model, the logical model. In this step, it is important to use as few variables as possible. To this end, we have applied two optimization steps compared to mechanical conversion for the FIFO Stack-Up problem. Secondly, the model is transformed into a model that fits the target hardware, the physical model. The process of transforming a logical model into a physical model is called embedding. The physical model is also sometimes referred to as an embedding. An example for an embedding for a tiny configuration of a FIFO Stack-Up problem is given in section “Pallet Stack-Up Solution”. For this report embeddings were generated using the D-Wave implementation of the heuristic algorithm presented in [19].

If one logical variable is mapped to a chain of multiple qubits all qubits belonging to this chain should have the same value. The coupling strength between the qubits of a chain is set to a negative value in order to make them likely to take an equal value. In this case the physical model has more variables than the logical model and it is less likely for the annealing process to return optimal solutions, because there is no way to guarantee that the qubits in a chain are equal for a reasonable coupling strength. If this is not the case, a chain breaks, which becomes more likely when the chains get longer. There are different ways to resolve chain breaks. For this report, chain breaks were resolved using majority voting, meaning that the value of the logical variable in the result will be the value most common among its physical counterparts, because this method gave the best results.

In this work the definition of the FIFO Stack-Up problem presented in [20] was used. The elements relevant to this application will be reiterated in this section. The problem arises in delivery industry, where bins have to be stacked-up from conveyor belts onto pallets with respect to customer orders. The bins reach the palletizer on the main conveyor of an order-picking system. A distribution conveyor pushes the bins out to several buffer conveyors. Robotic arms are placed at the end of these buffer conveyors where each arm picks-up the first bin of one of the buffer conveyors and moves it onto a pallet located at stack-up places. Full pallets are carried away by automated guided vehicles, or by another conveyor system, while new empty pallets are placed at free stack-up places.

The problem how to distribute bins to the buffer conveyors is investigated in [21]. In this paper, we assume the assignment of the bins to the buffer conveyors and the order of bins within each conveyor as given. Such systems can be modelled by several simple queues that works on the first in first out principle. The FIFO Stack-Up problem is to decide for a given list Q of k sequences of bins and a positive integer p whether there is a processing of Q, such that in each situation during the processing of Q at most p pallets are open. A pallet t is called open, if one but not all bins for t have already been removed from the sequences. An example of this problem is illustrated in Fig. 1.

Fig. 1
figure 1

Left: A FIFO stack-up system with an upstream distribution conveyor. Right: An example of the FIFO Stack-Up problem. There are 2 sequences and 3 stacking places currently in use

Full size image

Each of the buffer conveyors is represented as a sequence

$$\begin{aligned} q_1= & {} (b_1,\ldots ,b_{n_1}), ~ \ldots ,\\q_\ell = & (b_{n_{\ell -1}+1}, \ldots , b_{n_\ell }), ~ \ldots , ~\\ q_k =& {} (b_{n_{k-1}+1}, \ldots , b_{n_k}) \end{aligned}$$

of bins. At each step, only the first unremoved bin in a sequence can be removed. In Fig. 1 this would be the bin at the bottom of the conveyor. To make notation easier, it is defined that bin \(b_i\) is to the left of bin \(b_j\) if they are in the same sequence and \(i<j\).

Each bin b is labeled with a symbol plt(b). Bins with the same symbol must be stacked on the same pallet. Conversely, bins with differing symbols cannot be stacked on the same pallet. The concrete values of the labels are not important, since the only important consideration is whether two bins have the same or different labels. In this report positive integers or characters are used as labels. Unless otherwise stated, n denotes the number of bins, m denotes the number of pallet symbols.

The position of the first bin with a label t in a sequence \(q_i\) is denoted by \({ first}(q_i,t)\). The position of the last bin with label t in the sequence \(q_i\) is denoted by \({ last}(q_i,t)\). To simplify later definitions, when there is no bin with label t in sequence \(q_i\) let \({ first}(q_i,t) = |q_i|+ 1\) and \({ last}(q_i,t) = 0\).

Let \(Q=(q_1,...,q_k)\) be a list of sequences and let \(C_Q=(i_1,..., i_k)\), \(0 \le i_j \le \vert q_j \vert\), be some tuple in \({\mathbb {N}}^k\). Then \(C_Q=(i_1,...,i_k)\) is called a situation, and each value \(i_j\) in \(C_Q\) denotes the number of bins removed from sequence \(q_k\). Since bins can only be removed from each sequence in the order defined by the sequence, this completely defines the situation.

A pallet t is called open in situation \(C_Q\) if there is an \(i_j\in C_Q\) with \({ first}(q_j, t)\le i_j\) and an \(i_\ell \in C_Q\) with \({ last}(q_\ell , t)>i_\ell\). If \({ last}(q_j, t)\le i_j\) for each \(q_j\in Q\) the pallet is called closed. Initially, the pallet is called unprocessed. For the purpose of this report a pallet being closed or unprocessed is interchangeable.

The FIFO Stack-Up problem is shown to be NP-complete in [2] even if the number of bins per pallet is bounded. A decision problem asks whether something is true, in our case, whether there is a processing of the bins such that no more than k stack-up places are used. In an optimization problem we have to find among all solutions the one with the best performance according to some metric. In our case, the optimization version of the problem asks for a processing of the sequences in Q with the lowest possible number of pallets p open at any one step.

Example 1

Let \(Q=(q_1,q_2)\) where

$$\begin{aligned} q_1 = (b_{1},\ldots ,b_{4}) = [a,b,a,b] ~~ \text{ and } ~~ q_2 = (b_{5},\ldots ,b_{10}) = [c,d,c,d,a,b]. \end{aligned}$$

For n not necessarily distinct pallets \(t_1, \ldots , t_n\) let \([t_1,\ldots ,t_n]\) denote some sequence of n pairwise distinct bins \((b_1, \ldots , b_n)\), such that \({ plt}(b_i) = t_i\) for \(i=1,\ldots ,n\). Table 1 shows a processing of \(Q=(q_1,q_2)\) with 2 stack-up places. The underlined bin is always the bin that will be removed in the next transformation step. The already removed bins are shown greyed out.

Table 1 A processing of \(Q = (q_1, q_2)\) from Example 1 with 2 stack-up places
Full size table

There are at least two obvious ways to describe a solution of the stack-up problem. Consider a processing of a list Q of sequences. Let \(B = (b_{\pi (1)}, \ldots , b_{\pi (n)})\) be the order in which the bins are removed during the processing of Q, and let \(T = (t_1, \ldots , t_m)\) be the order in which the pallets are opened during the processing of Q. Then B is called a bin stack-up solution of Q, and T is called a pallet stack-up solution of Q. The transformation in Table 1 defines the bin stack-up solution

$$\begin{aligned} B=(b_5, b_6, b_7, b_8, b_1, b_9, b_2, b_{10}, b_3, b_4), \end{aligned}$$

and the pallet stack-up solution \(T = (c,d,a,b)\).

In the next section we convert both solutions into QUBO formulations and discuss possible optimizations to reduce the number of logical variables. The total number of variables as a function of the number of bins and pallets are derived in section “Calculating the number of variables” while in section “Experimental Results” an experimental evaluation for a few small instances of the FIFO Stack-Up problem on the D-Wave Quantum Annealer is given. Finally in section “Conclusion” we draw a conclusion.

QUBO-Formulations of the Linear Programming Solutions

In [20] two linear programming solutions for the stacking problem are presented. A linear program defines a value to be minimized subject to certain constraints that must be satisfied. For an introduction to mixed-integer linear programming see [22]. A linear program can be converted into QUBO-formulation by constructing a term that takes very large values if any of the constraints are violated. Specifically, a solution that violates any constraint must yield more energy than any solution that doesn’t violate any constraint. The easiest way to ensure this is to multiply the constraint terms by a weight A. For a linear program with constraints \({C_1({\textbf{X}}), ..., C_k({\textbf{X}})}\) that minimizes the value v the energy function \(E({\textbf{X}})\) has the following general structure:

$$\begin{aligned} E({\textbf{X}}) = A\cdot (C_1({\textbf{X}})+...+C_k({\textbf{X}})) + v \end{aligned}$$

A must be chosen to be larger than the largest possible value of v. This ensures that a violation of a constraint costs more energy than any possible value of v. In our case v is the number of stack-up places.

Since the variables of the energy function can only take two different values, natural numbers, such as the value minimized by the linear program, must be encoded by multiple variables. This is usually achieved by creating multiple binary variables to represent a number and modifying their weights, so activating different combinations of binary variables will contribute to the appropriate energy. In this report, the log-trick presented in [23] is used. This essentially means encoding natural numbers in QUBO-formulations as binary numbers. For example, the term encoding the natural number v with the binary variables \(v_0,..., v_{B-1} \in \{0,1\}\) is \(\sum _{i=0}^{B-1}2^{i}\cdot v_{i}\). Value B must be chosen so that encoding can represent the largest possible value of v. Since the encoding is essentially just a binary representation, the largest value that can be represented with this encoding is \(2^B-1\).

Before we take a closer look at the QUBO-formulations, let us make a short remark. The FIFO Stack-Up problem does not have to be solved by using an integer program. In [24], a solution is described by some breadth-first search procedure, which can solve instances from practice quite efficiently. On the contrary, it was shown in [25] that a solution by means of integer programs is not efficient. Moreover it is not clear if the FIFO Stack-Up problem in QUBO formulation can be solved at all on a quantum annealer, since other combinatorial problems like the knapsack problem could not be solved on a quantum annealer so far [26]. Furthermore other problems could only be solved for very small instances or the best solutions were not found [27, 28].

In [29, 30] they increased the system size they could solve by splitting instances up into subproblems submitted to D-Wave’s quantum annealer, but still endup with quite small systems. However, integer programs are widely used in industry and the operations research community and sophisticated solvers are provided by many companies such as IBM (cplex), Gurobi (gurobi) or SAP (genios). We want to analyse if, by using a quantum annealer, such stack-up problems can be solved more efficiently than using breadth-first search or conventional solvers. For this analysis it is important to use as few logical variables as possible as well as to reduce the errors,Footnote 2 if possible, by alternating the input QUBO.

Bin Stack-Up Solution

The idea of the bin stack-up solution is to explicitly plan in which order the bins are removed while minimizing the amount of stacking places used concurrently.

To represent a plan, \(n^2\) binary variables \(x_i^j\) can be used, where n is the number of bins of the problem instance. For each bin there are n plan variables, one for each step. \(x_i^j\) is set to one if and only if bin i is removed at step j. The following linear program for calculating the bin stack-up solution is given in [20].

$$\begin{aligned} \text {minimize } p \end{aligned}$$
(4)

subject to

$$\begin{aligned} { PERMUTATION}(n, x_i^j) \end{aligned}$$
(5)

and

$$\begin{aligned} { SEQUENCE\_ORDER}(Q, x_i^j) \end{aligned}$$
(6)

and

$$\begin{aligned} \sum _{t=1}^m f(t,c) \le p \text { for every } c \in [n-1] \end{aligned}$$
(7)

and

$$\begin{aligned} f(t,c)=(\bigvee _{\begin{array}{c} i \in [n] \\ j \le c \text {, } plt(b_i) = t \end{array}} x_i^j)\wedge (\bigvee _{\begin{array}{c} i \in [n] \\ j> c \text {, } plt(b_i) = t \end{array}} x_i^j). \end{aligned}$$
(8)

The constraint \({ PERMUTATION}\) (5) demands that each bin is only removed once and exactly one bin is removed at every step. To realize a bijection \(\pi :[n]\rightarrow [n]\) we define \(n^2\) binary variables \(x_i^j\in \{0,1\}\), \(i,j\in [n]\), such that \(x_i^j\) is equal to 1, if and only if \(\pi (i) = j\). In order to ensure \(\pi\) to be surjective and injective, we use the conditions \(\sum _{i=1}^{n} x_i^j = 1\) for every \(j \in [n]\) and \(\sum _{j=1}^{n} x_i^j = 1\) for every \(i \in [n]\). Further we have to ensure that all variables \(x_i^j\), \(i,j \in [n]\), are in \(\{0,1\}\).

By \(\pi\) the relative order of the bins of each sequence \(q \in Q\) has to be preserved, consider for example the bin solution of the sequences of Example 1. The bins of the corresponding sequence are shown in black, the others are greyed out.

The constraint \({ SEQUENCE\_ORDER}\) (6) demands that the ordering of the sequences is respected, meaning that a bin can only be removed if all of the bins preceding it in its sequence have been removed at a previous point in time. For every sequence \(q_s\), \(1\le s \le k\), and every bin \(b_i\) of this sequence, i.e. \(n_{s-1}+1 \le i \le n_s\), we know that if \(b_i\) is removed at position j, i.e. \(x_{i}^{j}=1\), then

  • every bin \(b_{i'}\), \(i'>i\) of sequence \(q_s\), i.e. \(i < i' \le n_s\), is not removed before \(b_i\), i.e. \(x_{i'}^{j'}=0\) for all \(j'<j\), which is ensured by \(x_{i'}^{j'}\le 1-x_{i}^{j}\) and

  • every bin \(b_{i'}\), \(i'<i\) of sequence \(q_s\), i.e. \(n_{s-1}+1 \le i'< i\), is not removed after \(b_i\), i.e. \(x_{i'}^{j'}=0\) for all \(j'>j\), which is ensured by \(x_{i'}^{j'}\le 1-x_{i}^{j}\).

The constraints given by Eqs. (7) and (8) cause p to be set to at least the number of stacking places used concurrently at a step c. f(t, c) is equal to 1 if and only if a stacking place is required for label t at step c. This is the case if and only if in the considered ordering of the bins there is a bin \(b'\) with \({ plt}(b')=t\) opened at a step \(\le c\) and there is a bin \(b''\) with \({ plt}(b'')=t\) opened at a step \(>c\), if and only if pallet t is open after the c-th bin has been removed. Taking the sum of f(t, c) over all labels gives the number of stacking places used at step c.

To create a QUBO-term for Eq. (5), the \({ PERMUTATION}\) constraint, it is interpreted as exactly-one-true constraints in two directions. Once in the direction of i for each j and once in the direction of j for each i. A simple formulation of an exactly-one-true constraint is given in [23]. Applying that formulation here yields the following term.

$$\begin{aligned} P_{bin}({\textbf{X}}) = \sum _{i=1}^n\left(\left(\sum _{j=1}^n x_i^j\right) -1 \right)^2 + \sum _{j=1}^n\left(\left(\sum _{i=1}^n x_i^j \right) -1 \right)^2. \end{aligned}$$
(9)

It is worth noting that the formulation of this term using a QUBO-matrix adds some constants that cannot be represented in a QUBO-matrix because they are not multiplied with any variables. While omitting these constants changes the value of the minima by \(-2n\), it does not change their location, so the energy function still becomes minimal only for optimal solutions.

The \({ SEQUENCE}\_{ ORDER}\) constraint stated in Eq. (6) can be interpreted as certain plan variables being forbidden from being equal to 1 at the same time. Since \(x_i^j\in \{0,1\}\) this can be formulated as a sum of products

$$\begin{aligned} S_{bin}({\textbf{X}})=\sum _{i'=1}^n\sum _{j=1}^n\sum _{\begin{array}{c} i=i'+1 \\ \exists Q: x_i\in Q\wedge x_{i'}\in Q \end{array}}^n\sum _{j'=j+1}^n x_i^j x_{i'}^{j'}. \end{aligned}$$
(10)

For the sake of this formulation, the indices of the bins are considered to be incremental in each sequence, meaning if \(x_{i'}\) occurs before \(x_i\) in a sequence, then \(i'<i\). By pursuing this approach instead of mechanically converting the inequalities used to implement this constraint in the converted program [20] we can avoid introducing additional slack variables which would make the model more complex and harder to solve.

Table 2 QUBO equivalents of relevant Boolean expression formulations given in [31], where z is usually an auxiliary variable that holds the result of the expression
Full size table

Next, the Boolean expressions (Eq. 8) will be converted, which are used for the inequality (Eq. 7). To represent a Boolean expression in an energy function, an auxiliary variable must be created to hold the result of each Boolean operation. A term is added that takes its minimal value when this auxiliary variable is set correctly. For pairwise Boolean operation such terms are given in Table 2. For longer expressions, the expression can be solved pairwise by substituting in the auxiliary variable that holds the result of the previous junction for either a or b.

This means that converting Eq. (8) into a QUBO-formulation will generate 2 types of auxiliary variables. First there will be intermediate variables, which hold intermediate results to be reused in later pairwise expressions to construct larger expressions. Secondly there will be the final auxiliary variables, which will hold the value that each f(t, c) should take if the stack-up solution is valid. These final auxiliary variables will be called \(aux_{f(t,c)}\). They will be used for the QUBO-formulation of Eq. (7).

The process of generating these auxiliary variables can be realized by a function that takes two arguments, the operator to be applied, either a Boolean OR or a Boolean AND, and a stack of the Boolean variables to apply it to. For example, the expression \(F = a \vee b \vee c \vee d\) would be passed in as the operator OR and a stack containing a, b, c and d. The elements are then extracted from the stack in pairs and the appropriate term to model the Boolean expression is added to the QUBO-model. This creates an auxiliary variable which represents this pairwise expression. The auxiliary variable is then pushed back onto the stack. Using a stack here means that the auxiliary variable will immediately be used in the next pairwise expression. This increases the potential of reusing auxiliary variable from other calls, minimizing the overall amount generated. Once there is only one variable left in the stack, the terms created so far model the entire expression that was given to the function with the remaining variable taking the value of the overall result. For the given example \(F = a \vee b \vee c \vee d\), the algorithm would first combine a and b into an auxiliary variable \(x = a \vee b\) by adding the term representing the condition \(a\vee b\) to the QUBO-model as given in Table 2. The expression becomes \(x \vee c \vee d\). In the next iteration we combine \(y = x \vee c\) and the model becomes \(y \vee d\). Finally \(aux_F = y \vee d\) is created, which represents the entire original expression as one variable.

To generate as few auxiliary variables as possible, some observations about the specific Boolean expressions that arise in Eq. (8) can be made. Firstly, the expressions f(t, c) for different values of t use separate variables since a bin can only have one label. Therefore they can be processed independently from each other. Secondly, whether a specific plan variable \(x_i^j\) appears in the left disjunction or the right disjunction depends only on the value of j, which is compared with the parameter c of f(t, c).

A \(\vee\)-operation connecting all bins of one label and one j value is called a block. Assuming there are q bins of label t, numbered \(t_1,\ldots t_q\), a block is \(\big (x_{t_1}^j\vee \ldots \vee x_{t_q}^j\big )\).

For different values of c, the plan variables will move between the disjunctions in blocks with equal j which are never split up. Therefore it makes sense to first define auxiliary variables which can stand in for these blocks during further processing.

Table 3 Structure of auxiliary variables for one label t
Full size table

The different expressions in these variables can be used to minimize the amount of auxiliary variables generated, which is important to get a larger system running on the D-Wave Annealer.

Table 3 illustrates this problem. The auxiliary variables represented by the underbraces show the ideal distribution with the least auxiliary variables. As can be seen from this representation, the optimal way to generate these auxiliary variables is in the same direction the expression grows, meaning the left disjunction should be read left to right and the right disjunction should be read right to left. This way, combining the same variables in multiple different auxiliary variables unnecessarily can be avoided. In this direction, every pair is combined into an auxiliary variable which now represents that pair, e.g. \(a_1\) and \(a_3\). This process can then be repeated for larger expressions, creating for example \(a_2\) and \(a_4\).

After the Boolean expressions are converted in this way for each label, there is an auxiliary variable \(aux_{f(t,c)}\) for each relevant parameter pair (t, c) that holds the appropriate value that f(t, c) would return in a valid stack-up solution. The sum of the terms generated by this function will be referred to as \(B_{bin}({\textbf{X}})\). These auxiliary variables can be used to model the inequality given in Eq. (7).

To model an inequality, it helps to convert it into an equality first because an equality is easy to represent as a quadratic function. The quadratic term to demand that some value a is equal to some value b is simply

$$\begin{aligned} Eq(a,b) = (a-b)^2. \end{aligned}$$

This term takes its minimal value 0 when a and b are equal. To allow natural numbers as variables they are replaced by an encoding as described at the beginning of this section.

To convert an inequality into an equality a slack variable \(s\in {\mathbb {N}}_0\) can be used. The slack variable makes up the difference between the compared values. For \(a\le b\) the equivalent expression would be \(a+s = b\) for some value \(s\ge 0\). The term in Eq. (11) represents the inequality in Eq. (7) in a QUBO-formulation. A separate slack variable \(s_c\) must be used for each value of c considered. The variables \(s_c\) and p are natural numbers and must be encoded. Let \(s_{c,i}\) denote the i-th bit of the binary representation of number \(s_c\), and let \(p_i\) denote the i-th bit of the binary representation of number p. Since the maximum value of p and the \(s_c\) is the number of pallets m, B is chosen to be \(\lfloor \log _2(m)\rfloor +1\).

$$\begin{aligned} I_{bin}({\textbf{X}}) = \sum _{c=1}^{n-1}\left(\sum _{t=1}^{m}aux_{f(t,c)} + \sum _{i=0}^{B-1} 2^is_{c,i} - \sum _{i=0}^{B-1} 2^ip_i\right)^2. \end{aligned}$$
(11)

Using these formulations of the constraints the final term of the energy function for the bin stack-up solution is:

$$\begin{aligned} \begin{aligned} E_{bin}({\textbf{X}}) =&A\cdot (P_{bin}({\textbf{X}})+S_{bin}({\textbf{X}})\\&+I_{bin}({\textbf{X}})+B_{bin}({\textbf{X}})) + \sum _{i=0}^{B-1} 2^ip_i. \end{aligned} \end{aligned}$$
(12)

A must be chosen so that \(A> 2^B-1\). For \(E_{bin}({\textbf{X}}) < A+ { offset}\) the solutions satisfy all constraints. Otherwise one or more constraints are violated and the solution is considered invalid. Due to the constant components in \(P_{bin}({\textbf{X}})\) the values will be offset by \({ offset} = -2n\cdot A\).

There are two more considerations to be made to reduce the number of variables used by the bin stack-up solution.

Firstly, some inequalities can be neglected. At the first point in time only one stacking place can be in use, since only one bin has been removed. As long as there are at least two bins with each label, this means that adding the first inequality as a constraint does not add any information. The same holds true for the last point in time, since there is only one bin left to be removed. Therefore, the first and last inequalities can be removed. This eliminates two slack variables and 2m values of f(t, c). This logic continues to hold for later steps. At the second and second to last steps at most two stacking places can be in use and so on. Therefore, removing the first and last k inequalities would mean that a value of p smaller than k could mean that the stack-up solution actually requires up to k stacking places. Consequently, the system does not strive towards an optimal solution when below k stacking places. This means that the problem is converted back into a decision problem that asks whether there is a stack-up solution with k stacking places. The answer is yes if \(E_{bin}({\textbf{X}})\le k\) and no otherwise. The larger k becomes, the less precise the solution, but less auxiliary variables are used.

Secondly, there are some plan variables that can never be set to 1 because of the position of the corresponding bin in its sequence. For example the second bin in a sequence cannot be removed at the first step, the third bin cannot be removed at the first or second step and so on. Inversely, the first bin must be removed at the latest when all other sequences are exhausted. The corresponding plan variables can therefore be removed. This saves these plan variables and some auxiliary variables in the Boolean equations for f(t, c).

Pallet Stack-Up Solution

The idea of the pallet stack-up solution is to plan the order, in which the pallets are opened instead of explicitly planning the order, in which the bins are removed. The order of opening the pallets can then be used to make the optimal decision about which bin has to be removed next. This method is also described in [20].

It is based on the directed pathwidth of the Sequence Graph of the problem Instance. The sequence graph is a directed graph \(G_Q=(V,A)\) in which each pallet in the problem instance has a corresponding node in V. There is an arc \((u,v)\in A\) if and only if a bin belonging to the pallet associated with u appears before a bin belonging to the pallet associated with v in one of the sequences of the given problem instance. In other words, the arc exists if the pallet associated with u must be opened before the pallet associated with v can be closed.

An algorithm that computes the sequence graph is given in [20]. For a given digraph G, a directed pathwidth can be calculated. The definition is also stated in [20] but will be restated here for convenience.

A directed path-decomposition of a digraph \(G=(V,A)\) is a sequence

$$\begin{aligned} {{\mathcal {X}}} = (X_1,..., X_r) \end{aligned}$$

of subsets of V, called bags, with the following three properties

  1. 1.

    \(X_1\cup ...\cup X_r = V\)

  2. 2.

    For each \((u,v)\in A\) there is a pair \(i\le j\) such that \(u\in X_i\) and \(v\in X_j\)

  3. 3.

    For all \(i,j,\ell\) with \(1\le i<j<\ell \le r\) it holds \(X_i\cap X_\ell \subseteq X_j\).

The width of a directed path-decomposition \({{\mathcal {X}}}\) is \(\begin{array}{c} \max \\ 1\le i\le r \end{array} |X_i|- 1\). The directed pathwidth of a graph G is the smallest width w over all possible directed path-decompositions of G.

Since the pallet stack-up solution computes the order in which pallets are opened instead of the order bins are removed, the number of plan variables is now dependent on the number of pallets m instead of the number of bins n. There are \(m^2\) plan variables \(x_i^j\), one for each possibility of opening a pallet i at step j. \(x_i^j =1\) if and only if pallet i is opened at step j. Since there are m pallets to be opened there are m steps to completely define the order.

In the following we use the fact that the directed pathwidth equals the directed vertex separation number [32]. For a digraph \(G=(V,A)\) on n vertices, we denote by \(\Pi (G)\) the set of all bijections \(\pi :[n]\rightarrow [n]\) of its vertex set. Given a bijection \(\pi \in \Pi (G)\) we define for \(i \in [n]\) the vertex sets \(L(i,\pi ,G)=\{u\in V \mid \pi (u)\le i \}\) and \(R(i,\pi ,G)=\{u\in V \mid \pi (u) > i \}\). Every position \(i \in [n]\) is called a cut. This allows us to define the directed vertex separation number for digraph G as follows.

$$\begin{aligned} \text {d-vsn}(G)= & {} \min _{\pi \in \Pi (G)} \max _{1 \le i \le \vert V \vert } \vert \{u\\{} & {} \in L(i,\pi ,G) \mid \exists v \in R(i,\pi ,G): (v,u)\in A\} \vert \end{aligned}$$

Given the sequence graph \(G_Q=(V,A)\) with node set V and edge set A, [20] states the following integer linear program to compute the pallet stack-up solution by calculating the directed vertex separation number w of \(G_Q\)

$$\begin{aligned} {\textit{minimize } }w \end{aligned}$$
(13)

subject to

$$\begin{aligned} { PERMUTATION}(m,x_i^j) \end{aligned}$$
(14)

and

$$\begin{aligned} \sum _{j=1}^c{Y(j,c) \le w} \text { for every } c \in [m-1] \end{aligned}$$
(15)

and

$$\begin{aligned} Y(j,c) = \bigvee _{\begin{array}{c} j'\in \{c+1,...,m\} \\ i,i'\in [m],(v_{i'},v_i) \in A \end{array}} (x_i^j \wedge x_{i'}^{j'}). \end{aligned}$$
(16)

The Boolean expression Y(j, c) of Eq. (16) is equal to 1 if and only if there exists an arc from a vertex on a position \(j'>c\) to a vertex on position j and the corresponding plan variables are set to 1. The sum \(\sum _{j=1}^c Y(j,c)\) of Eq. (15) is therefore the number of nodes that are at the end of arcs from positions right of c. The number of stacking places required by a given solution of this optimization problem is \(w+1\).

Just as with the bin stack-up solution, the linear program can be translated into a QUBO-formulation by creating a term for each individual constraint and adding w to the sum of those terms.

The constraint given in Eq. (14) is almost identical to its counterpart in the bin stack-up solution (Eq. 5). The constraint demands that each pallet must be opened exactly once and exactly one pallet is opened at each step. This can be modeled as two exactly-one-true constraints over the plan variables without auxiliary variables

$$\begin{aligned} P_{plt}({\textbf{X}}) = \sum _{j=1}^{m}\left(\left(\sum _{i=1}^mx_i^j\right) -1\right)^2+\sum _{i=1}^{m}\left(\left(\sum _{j=1}^mx_i^j\right) -1\right)^2. \end{aligned}$$
(17)

The term in Eq. (17) is equal to 0 if and only if the constraint given in Eq. (14) is satisfied. In any other case, the value will be greater than 0.

Next, the constraint given in Eq. (16) will be modeled. This constraint mandates, that the values Y(j, c) used in Eq. (15) are set correctly. This constraint can be transferred into QUBO-formulation in the same way as was done for Eq. (8) by creating auxiliary variables for each conjunction and each disjunction.

Some auxiliary variables can be saved by reusing some of the auxiliary variables representing Y(j, c). Because some disjunctions are created for each \(j'\in \{c+1,..,m\}\), the number of disjunctions grows larger the smaller c becomes and includes all disjunctions of the larger c values. Since disjunctions are associative, this means that Eq. (16) can be rewritten as:

$$\begin{aligned} Y(j,c) = \big (\bigvee _{\begin{array}{c} j'=c+1 \\ i,i'\in [m],(v_{i'},v_i) \in A \end{array}} (x_i^j \wedge x_{i'}^{j'})\big ) \vee Y(j,c+1) \end{aligned}$$
(18)
Table 4 The QUBO-matrix generated for the two sequences \(q_1=[0,1]\) and \(q_2=[1,0]\) to calculate the pallet stack-up solution
Full size table

As a break condition, it is stated that \(Y(j,c)=0\) if \(c > m-1\).

By rewriting the expression in this way, repeating equivalent expressions in the constraints for different Y(j, c) can be avoided, reducing the number of auxiliary variables.

Using the formulations of Boolean operators given in [31], an auxiliary variable holding the value of each relevant Y(j, c) can be created. These auxiliary variables will be called \(aux_{Y(j,c)}, 1 \le c<m, 1\le j\le c\).

First, each conjunction can be replaced by an auxiliary variable. The remaining disjunctions can be solved iteratively by substituting auxiliary variables and adding the necessary constraints from left to right (see the description in section “Bin Stack-Up Solution”). With these constraints in place, each \(aux_{Y(j,c)}\) will be set to the correct value. For ease of notation, the sum of the terms that ensure \(aux_{Y(j,c)}\) and the other auxiliary variables are set correctly will be called \(B_{plt}({{\textbf {X}}})\).

Lastly, Eq. (15) will be modeled. Since this is another inequality, slack variables \(s_j\) will again be used to convert the inequality into an equality. Furthermore, w and the slack variables will be represented using the log trick as described in section “QUBO-Formulations of the Linear Programming Solutions”

$$\begin{aligned} I_{plt}({\textbf{X}}) = \sum _{c=1}^{m-1} \big (\sum _{j=1}^{c}aux_{Y(j,c)} - \sum _{i=0}^{B-1}2^is_{j,i} - \sum _{i=0}^{B-1}2^iw_i\big )^2 \end{aligned}$$
(19)

with \(s_{j,i}, w_i\) being the i-th bit of the binary representation of \(s_j\) and w, respectively.

With all the constraints modeled the full energy function can be written as

$$\begin{aligned} E_{plt}({\textbf{X}}) = D(P_{plt}({\textbf{X}}) + B_{plt} ({\textbf{X}}) + I_{plt}({\textbf{X}})) + \sum _{i=0}^{B-1}2^iw_i \end{aligned}$$
(20)

with \(D > 2^B-1\). For a given solution \({\textbf{X}}\), \(E_{plt}({\textbf{X}})\) will be one less than the number of stacking places required if the solution is valid. If the solutions is invalid, meaning it violates at least one of the constraints, the energy will be greater than or equal to D. Because of the constant components in \(P_{plt}({\textbf{X}})\) the values will be offset by \(-2m\cdot D\) as before.

As an example consider a system with two labels 0 and 1, and 2 sequences. The order of the labels of the bins in the first sequence is [0, 1] and in the second sequence is [1, 0]. 11 variables are needed to compute the energy function (20). The corresponding matrix \(Q_{i,j}\) of the energy function (see Eq. (3)) is shown in Table 4. The factor D of Eq. (20) was chosen to be \(D=200\), which gave the best results concerning the number of correct solutions for this matrix.

One possible embedding for this system is shown in Fig. 2. It can be seen that even for only 11 variables it was not possible to find a one-to-one embedding of the QUBO-graph on the Pegasus topology. The variable Y(1, 1) is represented by 2 qubits. For a larger number of variables a lot more chains occur and the chains of qubits get much longer.

Fig. 2
figure 2

Embedding of the \(q_1=[0,1]\) and \(q_2=[1,0]\) system on the Pegasus topology with size parameter 2. A larger version of this topology with size parameter 16 is used in the Advantage System

Full size image

Calculating the Number of Variables

Because the feasibility of solving a problem on a quantum annealer depends largely on the number of variables, next to the structure of the couplings, it is worth considering how many variables are generated by these solutions given the problem parameters.

For these calculations the following variables are used:

  • n, the number of bins

  • m, the number of pallets

  • sl(s), the number of bins in the sequence s

  • lc(l), the number of bins with label l

Bin Stack-Up Solution

To calculate the number of variables generated by the bin stack-up solution the base amount will first be determined, then the number of variables saved by optimizations will be subtracted.

The number of plan variables is \(n^2\), one variable for the possibility of removing each of n bins at each of n points in time. n natural numbers must be encoded to calculate the bin stack-up solution, \(n-1\) slack variables and p. Each number is encoded with \(\lfloor \log _2(m)\rfloor +1\) binary representation variables.

To calculate the number of auxiliary variables generated to calculate Boolean expressions, each label must be considered separately because the amount depends on the number of bins with the label. As discussed in section “Bin Stack-Up Solution”, the first step is generating an auxiliary variable for each block. To generate a single auxiliary variable that represents the value of the disjunction of the lc(l) variables in a block for label l, \(lc(l)-1\) auxiliary variables are necessary. Since there are n blocks, this makes the total number of auxiliary variables to represent the blocks \(n\cdot (lc(l)-1)\) for each label. These auxiliary variables representing time steps are then further used to calculate the \(n-1\) required values of f(t, c). Because auxiliary variables are reused as much as possible, only the largest possible expression on each side of the split has to be considered. These expressions both have \(n-1\) elements. Therefore \(2\cdot (n-2)\) auxiliary variables are used to represent the disjunctions. Additionally \(n-1\) auxiliary variables are used to represent the conjunctions which provide the value of f(t, c). The expression for one label l is

$$\begin{aligned} BV_{bin}(n,l)&= {} (n(lc(l)-1)+2(n-2)+(n-1))\\&= {} n(lc(l)+2)-5. \end{aligned}$$

In total, the number of variables used to calculate the bin stack-up solution before optimizations is

$$\begin{aligned} V_{bin}(n,m)&= {} n^2 + n(\lfloor \log _2(m)\rfloor +1) \nonumber \\{} & {} +\sum _{l=1}^m (n(lc(l)+2)-5). \end{aligned}$$
(21)

Next, the number of variables saved by fixing plan variables that always have to be 0 because of the position of the corresponding bin in the sequence will be calculated. As described in section “Bin Stack-Up Solution”, a bin in the second position of a sequence cannot be removed in the first step, a bin in the third position cannot be removed in the first and second steps etc.. The number of bins that can be removed through this consideration for a sequence s are described by the sum

$$\begin{aligned} \sum _{i=1}^{sl(s)-1}i = \frac{(sl(s)-1)sl(s)}{2} \end{aligned}$$

meaning it is the \((sl(s)-1)th\) triangular number.

For the second consideration, that the first bin in a sequence must be removed at the latest after all other sequences are exhausted, similar logic can be applied. Specifically, it states that in the last \(sl(s)-1\) steps the first bin in the sequence must already have been removed, meaning that many plan variables can be omitted. The same holds true for later bins, reducing the number of plan variables by one less for each position in the sequence. This is described by the \((sl(s)-1)th\) triangular number again, meaning that the total number of plan variables that can be removed for a given sequence s is described by the following function.

$$\begin{aligned} R_{bin}(s) = sl(s)^2-sl(s) \end{aligned}$$
(22)

Because this method reduces the number of plan variables, the number of auxiliary variables necessary to represent the points in time for the calculation of f(t, c) is also likely to be reduced. Unfortunately, the exact number of variables that can be saved depends on the ordering of the sequences. The reason for this is that it takes \(v-1\) variables to represent a disjunction with v variables. This means that if there is more than 1 variable for a point in time, removing one of the variables will also remove an auxiliary variable. However, if a variable is the last variable for a label at a point in time, no auxiliary variable will be saved, because the variable could have represented that point on its own. This means it is impossible to calculate the exact number of auxiliary variables saved without considering the ordering, but an upper bound on the final number of variables can be given by assuming that no auxiliary variables are saved by this step.

A second optimization can be done by assuming that at least k stack-up places are needed. If the optimal solution for the instance does take up less than k stacking places it may be considered equal to a worse solution which uses up to k stacking places. For each label, this removes the inequalities for the k smallest c and the k largest c. Firstly, this saves 2k slack variables. Secondly, for each label the 2k corresponding f(t, c) will not have to be calculated. This means that the k largest disjunctions on the left and right side of the conjunction do not need to be calculated, saving one auxiliary variable each. Additionally, the auxiliary variables representing the conjunctions are also saved for a total of \(4\cdot m\cdot k\) auxiliary variables saved.

The final equation for the number of variables used by the bin stack-up solution considering the optimizations is given below.

$$\begin{aligned} V_{bin}(n,m,k)&= {} n^2 + (n-2k)(\lfloor \log _2(m)\rfloor + 1) \nonumber \\ &\quad + \sum _{l=1}^m(n(lc(l)+2)-5) \nonumber \\ &\quad - \sum _{s\in Q}(sl(s)^2-sl(s)) - 4mk. \end{aligned}$$
(23)

Pallet Stack-Up Solution

The plan for the pallet stack-up solution is represented by \(m^2\) variables, one for each pallet at each point in time.

There are \(m-1\) inequalities in (15) for which natural numbers to represent slack variables are necessary and one natural number to represent w. Each natural number is encoded by \(\lfloor \log _2(m)\rfloor + 1\) variables.

Due to Eq. (15), there are c different Y(j, c) for each value \(c\in [m-1]\). This comes out to the \((m-1)th\) triangular number

$$\begin{aligned} num_{Y(j,c)}(m)= \frac{m^2-m}{2} \end{aligned}.$$

For each Y(j, c) Eq. (18) has to be modeled for one value of \(j'\). This means modeling a conjunction for each edge of the sequence graph $G=(V,A)$, generating \(|A |\) auxiliary variables. Then the recursion needs to be modeled, which means modeling a disjunction with \(|A|+1\) variables. This generates another \(|A|\) auxiliary variables. However, for \(c=m-1\), the recursive part can be omitted since it is always 0, saving \(m-1\) auxiliary variables, one for each j considered in that specific inequality. Overall, the number of auxiliary variables used to model Eq. (18) is given by the following function

$$\begin{aligned} BV_{pal}(m,|A|) = \frac{m^2-m}{2}\cdot 2 \cdot |A|- (m-1). \end{aligned}$$
(24)

Combining and simplifying these components yields the following function for the number of variables used by the pallet stack-up solution

$$\begin{aligned} V_{pal}(m,|A|) = m^2(|A|+1)+m(\lfloor \log _2(m)\rfloor -|A|)+1. \end{aligned}$$
(25)

Experimental Results

Various small instances of the FIFO Stack-Up problem were examined in both QUBO formulations on the D-Wave Quantum Annealer and using the simulated annealing implementation by D-Wave. To simplify notation, the bins in the sequences are given by stating their respective labels \(plt(b_i)\) in place of \(b_i\). A list of all instances is given in Table 5. The smallest system had 2 sequences and 2 labels, \(q_1=[0,1]\) and \(q_2= [1,0]\). For this system 26 variables are needed for the bin stack-up solution and 11 for the pallet stack-up solution.

Table 5 List of instances of the FIFO Stack-Up problem with the optimal number of stacking places opt., examined on the D-Wave Quantum Annealer and using the simulated annealing implementation of D-Wave
Full size table

First the parameter A in the energy equation for the bin stackup solution (12) and the parameter D in the energy equation for the pallet stackup solution (20) were adjusted for a few instances in such a way that the most optimal solutions were found. Then these parameters were kept fixed for all other instances.

Next we modified several parameters to reduce the errors and to increase the number of optimal solutions found. We modified chain_strength, annealing_time and chain_break_method. Furthermore we alternated the normalization of the biases \(h_i\) and the strength of the couplers \(J_{i,j}\) in multiple ways. For all test the results didn’t show any significant improvement compared to the default values in the number of optimal solutions within the statistical errors. Therefore all shown results are gained using the default parameters of the solver.

Table 6 Results for Instances given in Table 5 for 25 runs with 4000 annealing cycles for each problem instance formulated with the bin stack-up solution, returned by the Advantage System 6.1
Full size table

The results for the configurations analyzed using the QUBO formulation based on the bin stack-up solution are given in Table 6, ordered in the number of variables and qubits needed. For the QUBO formulation based on the pallet stack-up solution the results for the configurations are given in Table 7. All simulations were done 25 times with 4000 annealing cycles each.

It can be seen that quite a lot of violations of the constraints for very small instances occurred. For the smallest instance the mean value of the probability of finding a correct solution is only 0.483 and for finding an optimal solution 0.158 with a standard deviation of 0.038 using the pallet stack-up solution with 11 variables. For the bin stack-up solution with 26 variables the probabilities drop down to 0.08 for a correct solution and 0.043 for an optimal solution with the standard deviation of 0.014.

For the bin stack-up solution valid stack-up solutions are only returned for very small problem instances. This can be counteracted somewhat by assuming that at least k stack-up places are needed instead of solving the general optimization problem reducing the number of variables and effectively making the problem smaller. Increasing k in the bin stack-up solution decreases the complexity of this system of terms and the number of variables, which causes correct solutions to become more likely. On the other hand, an optimal solution using less than k pallets can not be found or might be incorrect. The largest system which gave correct results was a system with 3 sequences and 3 labels, but for this system no optimal solution was found. For configurations where more than 62 variables are needed, no optimal solution could be found.

Furthermore, it appears that the constraint that is violated the most is f(t,c), which describes the Boolean expressions. This is to be expected, because the terms that describe the expressions are comparatively complex. As can be seen in the tables the number of qubits increases rapidly with increasing number of variables. This is the case, because the embedding on the Pegasus topology leads to more and more chains with increasing complexity as can be seen in the last rows of Tables 6 and 7, therefore the computation is prone to error. Increasing the chain’s strength to avoid chain breaks leads to less violations of the constraints, but also to less optimal solutions, because the coupling between qubits representing variables decreases.Footnote 3

Table 7 Results for instances given in Table 5 for 25 runs with 4000 annealing cycles for each problem instance formulated with the pallet stack-up solution
Full size table

Results for the pallet stack-up solution show the same behaviour as the results for the bin stack-up solution, but due to the fact, that less variables are needed, a slightly larger system can still be solved on the quantum annealer. Again, for all systems with up to about 60 variables an optimal stack-up solution was found. Instances 4, 1, and 2 of Table 7 as well as instances 8, 6, and 3 show that the success probability of finding the optimal solution are about the same under similar or identical conditions, i.e. comparable number of variables, number of chains and the maximal length of the chains. We will investigate in a future work what influence the embedding has on the quality of the solutions.

A similar behaviour can also be seen in [33] for the MAX-CUT problem, where the probability of success decreases if the number of vertices in the graph increases or the maximum vertex degree increases. For a comparison of the success probability between the pallet stack-up solution and the bin stack-up solution as a function of the number of variables, look at instance 3 with \(k=2\) of the bin stack-up solution and the instance 7 of the pallet stack-up solution. Both problems need 49 logical variables and similar number of qubits, 100 for instance 3 with \(k=2\) and 83 for instance 7. The success probability is in the same order, but a bit better for instance 7 with less qubits. Therefore it seems that the differences in the details of the QUBO matrices can be seen in the number of qubits and therefore also in the success probability, but don’t play an important role. Overall it can be seen that the number of variables is the dominant quantity for the success probability followed by the number of qubits or chains.

In addition to solving the QUBO-formulations using the Quantum Annealing, D-Wave also offers a classical simulated annealing solver. Both methods were compared to put the current capabilities of quantum annealer hardware into perspective.

Table 8 Probability of violated constraints for different problem instances formulated with the bin stack-up solution with \(k=1\)
Full size table
Table 9 Probability of violated constraints for different problem instances formulated with the pallet stack-up solution
Full size table

The results in Tables 8 and 9 show, that simulated annealing is capable of solving systems larger than what quantum annealing is capable of. Even on problem sizes quantum annealing was not able to solve, most anneals return solutions that at least satisfy all constraints. One should keep in mind that the parameters of the energy functions were adjusted such that the number of optimal solutions were optimized. Another possibility would be to optimize the number of correct solutions which satisfy all constraints.

The results show that simulated annealing using the QUBO-formulation still fails to return any valid solution for relatively small problem sizes, when the number of variables increases to a few hundreds. This is of the same order as the number of the physical variables or qubits needed on the quantum annealer, where we didn’t get any optimal solution. However, as shown in [20], the linear programs converted here work for larger problem sizes, though still below realistic sizes.

Conclusion

We developed two QUBO objective functions for solving the FIFO Stack-Up problem on a quantum annealer. We minimized the number of logical variables and discussed the complexity in the number of bins and pallets. The results on a quantum annealer show that it is generally possible to solve the FIFO Stack-Up problem on a quantum annealer, however only for tiny instances. Reducing the variables is not effective enough to solve major problems. Even if the performance of modern hardware increases, the basic problem remains: Large instances cannot be solved. Varying several parameters coudn’t improve the number of optimal solutions found either. For the pallet stack-up solution formulation less variables and therefore also less qubits are needed compared to the bin stack-up solution formulation for the same instance. Therefore for the pallet stack-up solution formulation slightly larger instances could be solved successfully. However this approach does not scale well. On one hand, the combination of using a penalty factor to enforce the constraints of a linear program and large weights generated encoding natural numbers cause some outlier weights that effectively reduce the energy range available for other weights. On the other hand, this is caused by hardware limitations necessitating embedding. The embedding process amplifies the inherent variable growth with problem size. To overcome such limitations, one can try to improve the embedding process, as described in [17] or [34], among others. Another possibility would be to use a direct embedding of Boolean functions as presented in [31]. This would leave the problems of embedding the other elements of the models and connecting the Boolean and non Boolean parts. Furthermore splitting instances up into subproblems could be tried. We will investigate these possibilities in our future work and we will do a detailed analyses of the influence of chains on the number of optimal solutions.