Quantum column generation for battery-constrained vehicle routing problems

Solving vehicle routing problems to optimality is particularly challenging, not only due to their NP-hardness, but also because naive model formulations often yield weak relaxations, an issue that leads standard solvers, like CPLEX or Gurobi, to perform particularly poorly. In turn, tight formulations of vehicle routing use variables that represent feasible routes. This causes the model to have a number of variables that grows exponentially with problem size, which makes even the relaxed (continuous) problem intractable.

Column generation is a method for solving Linear Programs (LPs) with a large number of variables, such as the relaxed tight formulation of vehicle routing. For this, a restricted version of the LP is solved, where only a subset of variables (routes) is considered, while the other are fixed to zero. An auxiliary problem, known as pricing problem, is solved to find new routes that must be added to the master problem in order to improve the objective function. In short, the master problem selects the best routes out of a restricted pool of possibilities, while the pricing problem introduces new routes into the pool. The algorithm alternates between solving each of the two problems, until no new improving routes are found by the pricing problem. Crucially, the effectiveness of the column generation algorithm relies on the ability to solve the pricing problem efficiently.

For the application of battery-constrained routing, the pricing problem is an NP-hard binary problem and can be cast into QUBO form. Particularly, the pricing problem takes the form

where the variables x_ij represent whether the route visits location j after location i, 0 is the index of the depot, Q is the battery capacity, N is the number of locations, p_i is the prize for visiting location i, and { λ_k^* }_k=0^N comes from the master problem.

The formulation above is known as the arc formulation. In (Kanai, 2024), a different model is used with variables q_it representing whether location i is visited at time step t. This results in a quadratic objective function. Using the model of (Kanai, 2024) in our battery-constrained setting would further result in a quadratic constraint. This poses the following problem:

• Classical solvers deal better with linear objectives and constraints.
• In transforming such a model to QUBO form, the quadratic constraint cannot be efficiently translated into a penalty term.

For these reasons, we choose to use the arc formulation instead. This comes with a caveat: this formulation allows solutions with subtours, i.e., disconnected routes, as opposed to a single route. Subtours need to be eliminated by adding constraints that forbid them. Since the number of possible subtours is exponential, we add these constraints (or penalties in the case of QUBO) on-the-fly.

We implemented the complete column generation algorithm and compared different methods for solving the pricing problem. Specifically, we compare

1. Pricing formulated as a QUBO. The constraints are encoded as penalty terms, which are initially estimated using crude estimates, and subsequently updated using the Augmented Lagrangian Method (ALM) (Djidjev, 2023). The resulting QUBO is solved using simulated annealing (SA) and quantum annealing (QA) in D-Wave's Advantage 4.1 system. Note that both methods are heuristic, e.g., do not guarantee finding optimal solutions.
2. Pricing formulated as a Binary Linear Program. Due to the exponential number of constraints needed to eliminate all possible subtours, we implement a framework that solves a restricted version of the problem via Gurobi and adds constraints on the fly as needed. This method is still guaranteed to find the optimal solution.

The experiments are run on synthetic instances with number of locations 5 ≤ N ≤ 15 and a central depot. For each value of N, ten instances are generated. The locations and prizes are chosen randomly. The vehicle capacity is set to 20 and the number of vehicles to √N.

Figure 2 shows a comparison of Gurobi versus SA, the latter with two variants of penalty tuning methods (naive and ALM). In terms of total collected prize, Gurobi demonstrates a modest improvement over SA as problem size grows. In computational time however, Gurobi is orders of magnitude faster than SA. More sophisticated penalty tuning does not appear to help in bridging this gap.

Figure 3 compares SA to QA on a subset of instances. This is due to both limited hardware access, as well as hardware limitations. Because of the high density of the solved QUBOs we estimate that we can solve problems of up to N = 10 locations. As shown in Figure 3, QA achieves worse solution quality with longer end-to-end execution time. It must be noted that the QA runtime includes the connection overhead and embedding time. The connection overhead cannot be reduced using parallelization techniques, due to the sequential nature of column generation, where each subproblem depends on the solution to the previous subproblem.

Figure 1: Maximum size of the QUBO as a function of the number of clients. These sizes were obtained experimentally. Note that the QUBO sizes do not increase monotonically, due to the varying number of constraints.

Figure 2: Comparison between Gurobi and SA with two variants of penalty tuning methods, naive and ALM.

Figure 3: Comparison between QA and SA.

Our experimental study indicates that the proposed methodology offers no computational advantage for the problem under consideration. The reason in inherent to the structure of the underlying pricing problem.

Most notably, the QUBO derived from the model in Eq. (1) is almost fully connected (see Figure 4). This is also the case for the formulation used by (Kanai, 2024). This issue is intrinsic to any routing problem, because of the fact that a vehicle can travel between any two locations. The density of the resulting QUBO severely limits the scalability of the whole approach.

Independently of the density issue, modelling routes inherently incurs in a large number of constraints, which substantially increases the difficulty of the problem, as appropriate penalties must be determined for each constraint. Moreover, the arc formulation which we use for battery-constrained routing, introduces an additional limitation: constraints must be added on-the-fly, resulting in a QUBO that grows in size and an energy landscape that becomes progressively more restricted. In particular, the common QA technique of reusing embeddings can no longer be applied.

These two main reasons prevent the approach to scaling up to problem instances large enough to challenge established solvers like Gurobi.

From a broader perspective, the quantum-assisted column generation (CG) approach should not be viewed as inherently deficient. Rather, our findings demonstrate that the scalability of this method is critically dependent on the structural properties of the underlying pricing problem. We identify an intrinsic and unsalvable limitation in the context of routing problems. However, this limitation does not necessarily generalize to other classes of combinatorial optimization problems, for which quantum-assisted CG may still offer meaningful advantages.

More generally, it must be emphasized that CG solves linear programs, and, as such, must be embedded in a larger framework that enforces integrality. Not all combinatorial problems benefit from exponentially large formulations like the one used in this study, and therefore the applicability of CG must be studied on a case-by-case basis.

Figure 4: Visual representation of the QUBOs encoding a typical pricing problem for the PCDRP (N = 15). The image on the left represents the pricing problem with the initial estimates for the penalty weights, and the image on the right represents the pricing problem after the last iteration of penalty updating.