{'fraction_non_alphanumeric': 0.02953259083439839, 'fraction_numerical': 0.012004747124338142, 'mean_word_length': 4.550798074486952, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multimodal Container Planning: a QUBO Formulation and Implementation on a Quantum Annealer

Quantum computing is developing fast. Real world applications are within reach in the coming years. One of the most promising areas is combinatorial optimisation, where the Quadratic Unconstrained Binary Optimisation (QUBO) problem formulation is used to get good approximate solutions. Both the universal quantum computer as the quantum annealer can handle this kind of problems well. In this paper, we present an application on multimodal container planning. We show how to map this problem to a QUBO problem formulation and how the practical implementation can be done on the quantum annealer produced by D-Wave Systems.


Introduction
In container logistics various transportation concepts are known [28]. Unimodal transport usually refers to loading cargo onto a single transportation mode, usually a truck, and driving it from origin to destination. Multimodal transport refers to transporting cargo from origin to destination by more than one transportation mode. In intermodal transport, again more than one transportation mode may be used, but the containment unit (container) must always be of standardised size. Co-modal transport looks like multimodal transport, but requires a consortium of shippers and has a focus on exploiting the benefits of each transportation mode in a smart way. Finally, synchromodal transport is a version of intermodal transport that focuses on real-time planning flexibility and coordination between different shippers, both using large amounts of real-time data.
When looking at the planning of container flows, it is often categorised on three levels of problems: on strategic, tactical and operational level [27]. Strategic problems concern long-term investments in the transportation network, for example, where to build new terminals. Tactical problems may concern service design, for example, determining how many times in a month a barge should make a round-trip. Operational problems concern using a current network in an optimal way for problems occurring in the present. Little attention has been on operational problems [22]. Operational problems can be divided into various classes: (a) problems in which containers are assigned to existing barge services, (b) problems in which the routes of barges are determined for a given demand, and (c) problems in which both the assignment of containers to barges and the route of barges are decided upon [29].
In an operational, dynamic, environment with a lot of uncertainty, or in a synchromodal environment where real time system characteristics are used in the planning, short computation times are crucial, which make it possible to repeat the planning process with a high frequency. To this end, efficient algorithms are proposed [28] or techniques are developed to reduce the complexity of the problem [19]. Next to software and algorithm enhancements to speed up the computation, also hardware developments can play a role. A promising direction here is quantum computing. Note that this will not be a solution for the short term, where the current generation quantum computers is far too small to compete with classical computers and their advanced solvers and algorithms. This work is meant to sketch a new direction for medium and long term solutions.
The past decade has seen the rapid development of the two paradigms of quantum computing, quantum annealing and gate-based quantum computing. In 2011 D-Wave Systems announced the release of the world's first commercial quantum annealer 1 operating on a 128-qubit architecture, which has since been continually extended up to the 2048-qubit version, available from 2017 2 . D-Wave announced in 2019 a (more than) 5000 qubit system available mid-2020, using their new Pegasus-technology based chip with 15 connections per qubit 3 . Even more recently, quantum supremacy is claimed to have been achieved by Google's 54-qubit Sycamore gate-based computer [2]. These technological advances have led to a renewed interest in finding classical intractable problems suited for quantum computing.
A particular category of problem that are expected to fit well on quantum devices are (combinatorial) optimisation problems. An important trail in implementing optimisation problems on quantum devices is the Quadratic Unconstrained Binary Optimisation (QUBO) problem. These QUBO problems can be solved by a Universal (gate-based) quantum computer using a QAOA (Quantum Approximate Optimisation Algorithm) implementation [11] or by a quantum annealer, such as D-Wave [18]. For already a large number of combinatorial optimisation problems the QUBO representation is known [15,20]. Well-known examples that have been implemented on one of D-Wave's quantum processors include maximum clique [5], capacitated vehicle routing [13], minimum vertex cover [24], set cover with pairs [4], Multi-Service Location Set Covering Problem [6], traffic flow optimisation [23] and integer factorisation [17]. These studies have shown that although the current generation of D-Wave annealers may not yet have sufficient scale, precision and connectivity to allow faster or higher quality solutions, they have the suitable infrastructure for modelling real-world instances of these problems, effectively decomposing these into smaller sub-problems and solving these on a real Quantum Processing Unit (QPU).The QUBO problems can also be solved by simulated (quantum) annealing on a conventional computer. Then, however, no quantum advantage can be expected.
In this paper we propose two QUBO formulations for a container assignment problem and show the implementation on both simulated annealing as well as on the D-Wave QPU. The size of the current generation of quantum computers limits us to small problems, with a limited number of decision variables. Therefor, we propose first a QUBO-formulation that only uses two potential paths for each container, and expand this to four alternatives per container.

Mathematical framework
In this work a deterministic container-to-mode problem is studied. To give a formal definition the framework of [8] is used. The followingR|D, [D2R]|social(1)problem is meant: to assign freight containers to transportation modes, such that the containers reach their destinations before a deadline against minimum total cost, given that the transport modes have fixed given schedules and all features of the problem are deterministic [16]. These problems are often solved using minimum cost multi-commodity flow problems on time-dependent graphs and space-time networks [1,9]. Finding a non-negative two-commodity flow on a directed graph, however, is proven to be NP-complete [10]. Finding a flow with minimum cost and with at least two commodities must be at least as difficult, however, the LP-relaxation of the studied problems almost always has an integral optimum [16]. To study the potential of quantum computing for this problem we propose a QUBO formulation.

Approach
Given the state of the quantum devices, which are quite small at this moment, we will limit the size of the problem by giving each container a limited number of paths through the transportation network, here either 2 or 4. For this we propose the following approach: 1. Select a set (2 or 4) of alternative routes per container; In the case of 2 alternatives we take the best multimodal path and the uni-modal (trucking) path. In the case of 4 alternative routes, we take again the uni-modal path and the 3 maximum dissimilar paths. 2. Define the QUBO-problem. 3. Solve the QUBO-problem. 4. Evaluate the solution and go back to step 1 if necessary to select different sets of alternatives.

QUBO Formulation
We now give the definition of the QUBO-problem [15] and propose the formulation of the QUBO for the two problems. The QUBO is expressed by the optimisation problem: where x ∈ {0, 1} n , the decision variables and Q is a n × n coefficient matrix. QUBO problems belong to the class of NP-hard problems. Another formulation of the problem, often used, equals or a combination of multiple of these terms where A, B, . . . are weights that can used to tune the problem and include constraints into the QUBO. We will use this representation in the next sections, where we will present the QUBO formulation for the 2 and 4 alternative route problems.
Creating a QUBO For already a large number of combinatorial optimisation problems the QUBO representation is known [15,20]. Many constrained integer programming problems can be transformed easily to a QUBO representation. Assume that we have the problem then we can bring the constraints to the objective value, using a penalty factor λ for a quadratic penalty: Using P = Ic, the matrix with the values of vector c on its diagonal, we get where matrix R and the constant d follow from the matrix multiplication and the constant d can be neglected, where it does not influence the optimisation problem.

2-alternative route QUBO
We look at a situation where containers have to be shipped in an intermodal or synchromodal network. To limit the number of options, we give every container two possible paths through the network, one using a truck only, the other a multimodal or multitrack path through the network. Each track or modality in this network has a specific capacity. We define a QUBO representation for this problem. We assume that all routes are feasible, with regard to due dates etc. The binary integer linear program for the problem is formulated as follows: such that Here for every container i it has to be decided whether it is sent using a truck (x i = 1) or using the multitrack multimodal path (x i = 0). Eq. 7 minimises the costs of the choices, where per container i the choice x i = 0 leads to costs c b i and Eq. 8 makes sure that the capacity constraints of the modalities or tracks are met. Based on the model in Eq. (7)(8), the QUBO formulation of the problem can also be derived. The overall solution for this QUBO is given by: with where A and B denote penalty coefficients to be applied such that the constraints will be satisfied and y ik denotes additional slack variables. These K · m binary slack variables are necessary to remodel Eq. 8 into equality constraints, as they are required in a QUBO-formulation. The parameter K follows from K = max j log 2 (V j ) . When determining the penalty coefficients, we can set A = 1 and look for a good value for B. Rule of thumb is that the gain of violating a constraint must be lower than the costs. This means in this problem that B > max i c b i . Note that if capacity is bounding for all tracks, Eq. 8 changes to an equality, which removes the need for the slack variables in Eq. 11. If the capacity is not bounding for all tracks, i.e.., the capacity is very high, Eq. 8 and Eq. 11 disappear from the problem for certain values of j, instead of allowing the slack variables for very high values. The binary integer linear program for the problem is formulated as follows:

4-alternative route QUBO
such that Here again Eq. 12 minimises the costs of the choices and Eq. 13 makes sure that the capacity constraints of the modalities or tracks are met.
Based on the model in Eq. (12)(13), the QUBO formulation of the problem can also be derived. The overall solution for this QUBO is given by: with where y ik denotes additional slack variables. These binary slack variables are necessary to remodel Eq. 13 into equality constraints, as they are required in a QUBO-formulation. Again, A, B and K denote penalty coefficients and the number of slack variables per track. They will be set conform the rules explained earlier.

Solving the QUBO problem
For solving the QUBO problem we use the D-Wave system. Implementing problems on a quantum device asks for specific Quantum Software Engineering [25]. We will sketch a number of specific implementation issues for this problem.
Quantum Annealing The quantum devices produced by D-Wave Systems are practical implementations of quantum computation by adiabatic evolution [12]. The evolution of a quantum state on D-Wave's QPU is described by a time-dependent Hamiltonian, composed of initial Hamiltonian H 0 , whose ground state is easy to create, and final Hamiltonian H 1 , whose ground state encodes the solution of the problem at hand: The system in Eq. (17) is initialised in the ground state of the initial Hamiltonian, i.e., H(0) = H 0 . The adiabatic theorem states that if the system evolves according to the Schrödinger equation, and the minimum spectral gap of H(t) is not zero, then for time T large enough, H(T ) will converge to the ground state of H 1 , which encodes the solution of the problem. This process is known as quantum annealing. Goal is to find the eigenstates of H, where the eigenvalues are the energy values. The eigenstates with the lowest eigenvalue are the ground states of the system. Although here we are not concerned with the technical details, it is worthwhile to mention that it is not possible to estimate an annealing time T to ensure that the system always evolves to the desired state. Since there is no estimation of the annealing time, there is also no optimality guarantee. The D-Wave quantum annealer can accept a problem formulated as an Ising Hamiltonian, corresponding to the term H 1 in Eq. (17), or rewritten as its binary equivalent, in QUBO formulation. Next, this formulation needs to be embedded on the hardware. In the most developed D-Wave 2000Q version of the system, the 2048 qubits are placed in a Chimera architecture [21].: a 16 × 16 matrix of unit cells consisting of 8 qubits. This allows every qubit to be connected to at most 5 or 6 other qubits. With this limited hardware structure and connectivity, fully embedding a problem on the QPU can sometimes be difficult or simply not possible. In such cases, the D-Wave system employs built-in routines to decompose the problem into smaller sub-problems that are sent to the QPU, and in the end reconstructs the complete solution vector from all sub-sample solutions. The first decomposition algorithm introduced by D-Wave was qbsolv [3], which gave a first possibility to solve larger scale problems on the QPU. Although qbsolv was the main decomposition approach on the D-Wave system, it did not enable customisations, and therefore is not particularly suited for all kinds of problems. The new decomposition approaches D-Wave offers are D-Wave Hybrid and the Hybrid Solver Service, offering more customisations.
Next to the QPU also a CPU can be used within the programming environment of D-Wave, using for example simulated annealing (SA), a conventional meta-heuristic.

Embedding
Some restrictions are introduced by the hardware design of the D-Wave hardware. The current implementation has a qubit connectivity based on a Chimera structure. The QUBO problem at hand has to be transformed to this structure. Because of the limited chip sizes we have currently, a compact formulation of the QUBO is important, but also a compact transformation to the chip design. This problem is known as minor-embedding. Here the vertices correspondent to problem variables and edges exist if Q ij = 0, where Q represents the qubo-matrix. Because of the limited connectivity of the chip, a problem variable has to be duplicated to multiple (connected) qubits. Those qubits should have the same value, meaning the weight of their connection should be such that it holds in the optimisation process. All these qubits representing the same variables are part of a so-called chain, and their edge weights is called the chain strength (λ), which is an important value in the optimisation process. In [7], Coffrin gives a thorough analysis. He indicates that if λ is sufficiently large, optimal solutions will match λ ≥ ij |Q ij |. However, the goal is to find the smallest possible value of λ, to avoid re-scaling of the problem. Coffrin also indicates that finding the smallest possible setting of λ can be NP-hard. Also other research has been performed on selecting an optimal choice for this chain strength [26], but at the moment there is no solid criterion for choosing a value. A rule of thumb that is suggested 4 is λ = max ij Q ij . It may be necessary to use the quantum annealer with multiple values of the chain strength in order to determine which value for λ is optimal for a given problem [14].

2-Alternative Route Case Study
We use a case consisting of 10 containers (n = 10) and 12 tracks (m = 12). The parameters are given in Table 3.  In the optimal solution, containers 4, 7 and 8 are transported by trucks and the others by barges, resulting in a total costs of 85, using the standard Matlab solver. As explained, this size of problems is no problem for classical solvers. The goal of this exercise is to show th potential of quantum computing when the quantum computers have more mature sizes. The QUBO following from Eq. 9 leads to a 46x46 matrix.

Simulated Annealing
The annealing processes (both simulated as quantum) are stochastic processes. This means that if we do a high number of queries we get a probability distribution for the resulting solution. The first important parameter to fix is this number of queries. We chose 500 samples. Next the values for A and B are important. We can set A = 1 and look for a good value for B. Rule of thumb is that the gain of violating a constraint must be higher than the consequent reduction in cost. This means in this problem that B > c b , so B > 10. In Figure 1 the probability distribution of four values for B. We see that choosing B too low (B = 6) gives all solutions that are better than the optimal solution, but they all violate one or more constraints. Choosing B too high gives solutions that are feasible, but they are (far) away of the optimal solution. Choosing B = 12 gives a range of solutions that include the optimal value, which will be found in most runs. The computation time for approach is around one minute, much more than conventional solvers will ask, who aks for less than one second.

Quantum Annealing
To use D-Wave's QPU now two parameters have to be determined. Again the value for A and B are important. By setting A = 1 we can choose a good value for B. Here we expect that the same value as established for simulated annealing will work. Here, we also have to determine a suitable value for the chain strength λ. We start with the rule of thumb indicated in the previous section, λ = max ij Q ij . Note that Eq. 1-3 imply that λ is dependent on the value of B, because B is part of the Q matrix. If we choose the values 3, 6 and 12 for B, then the values 120, 240 and 480 follow for λ. If we use λ ≥ ij |Q ij |, the values 4, 673, 9, 341 and 18, 687 follow. In Table 4 the results from the annealing are depicted. For each combination {λ, B} we present the minimum value, the average value and the percentage of feasible solutions found using 1000 readouts. The minimum value comes from the solution with the lowest energy that does not violate any of the constraints. The average value is the average objective value of all solutions that do not violate any of the constraints.  Note in Table 4 that the best (the optimal) solution is found in the left upper corner, for low values for both λ and B. The suggested value by the rule of thumb does not give good solutions in this example nor do the upper bound values. We also see that the best average values appear at a value λ = 10 for all values for B. This value for λ also comes with the lowest percentage feasible solutions, indicating that the probability distribution, like we see in Figure 1, shifts to the left for these values of λ. For two values (λ = 10 and λ = 240 and B = 6 the distribution is plotted in Figure 2, confirming this idea. Note that for λ = 240 the best solution found, based on minimum energy, has objective value 115. In Figure 2 we see that much better feasible solutions were found (even the optimal solution). These solutions follow from solution that have a chain break(s) and have a higher energy value. Same holds for λ = 10. Here the solution with objective value 90 is found, where 100 is reported as best solution. Calculation time of this approach is around 10 ms (pure QPU time) for each readout, leading to 10 second for 1000 readouts. This is more than conventional solvers ask for this size of problems, however, the idea is that Quantum Annealing will scale better for bigger problems, when the hardware is ready for that size of problems.

Conclusions
The rise of quantum computing and the promising application to combinatorial optimisation ask for re-formulation of many real-world problems to match the form these quantum computer can handle, the QUBO formulation, where QUBO stands for Quadratic Unconstrained Binary Optimisation. In this paper we proposed two QUBO-formulations for multimodal container planning, where sets containing a limited number of paths are used. For each container an assignment is found that minimises the overall costs of transporting all containers to their destination. These formulations are proposed to sketch to potential of quantum computing for this kind of problems when the quantum computer will be a more matured technology.
When implementing these problem formulations on D-Wave's quantum annealer, there are a number of implementation issues such as finding the right embedding, defining the chain strength and finding the right penalty functions. We showed how this is addressed by doing a parameter grid search.
For further research a more general idea for finding these parameters is recommended. Also the extension of the set with alternative paths is recommended and a integrated way to include promising paths to this set. Here we think of introducing a quantum variant of the known classical technique of column generation, Quantum Column Generation.