Annealing-Assisted Column Generation
for Inequality-Constrained Combinatorial Optimization Problems

We define $\Omega$ as the route set that includes all feasible routes in the specific CVRP instance. $\omega_{k}$ denotes the total distance of the $k$th route $r_{k}$. Let $A=[a_{i,k}]_{N\times|\Omega|}$ be a matrix whose element $a_{i,k}$ is $1$ if route $r_{k}$ contains the customer $v_{i}$ and $0$ if route $r_{k}$ does not contain the customer $v_{i}$. The set cover problem for the CVRP can then be represented by the following equation:

$\theta_{k}$ is the decision variable that takes $1$ if route $r_{k}\in\Omega$ is chosen and $0$ if $r_{k}\in\Omega$ is not selected for the solution. The optimization problem (12)–(15) is referred to as the master problem. Constraint (13) guarantees that every customer is visited at least once. Constraint (13) can be formulated as an equality constraint, but this typically slows down the convergence of the column generation. Consequently, an inequality form is typically used for this purpose, particularly during column generation. If the constraint (13) is expressed as an inequality, this problem is called the a set partition problem. The number of vehicles operated is determined by the constraint $\eqref{eq:mp_c2}$. Finally, the constraint (15) requires the decision variables to be integers. This is because the additional inequality constraint $\theta_{k}\leq 1$ can be eliminated in a restricted master problem, that underwent linear programming relaxation, as discussed later. Because the optimal $\theta_{k}$ is $0$ or $1$ in this model, $\theta_{k}\geq 2$ is not optimal.

The set cover problem is known to be $\mathcal{NP}$-complete, and the greedy algorithm is a well-known approximation algorithm for this problem. The number of feasible routes in $\Omega$ increases exponentially, and it is virtually impossible to consider all the feasible routes. Therefore, column generation focuses on a subset of $\Omega$, that comprises only promising routes. Let $\Omega^{\text{RST}}$ be such a restricted route set. Converting $\Omega$ in the set cover problem (12)–(15) into $\Omega^{\text{RST}}$ yields the following restricted master problem:

The restricted master problem (16)–(19) is differentiated from the original set cover problem (12)–(15) because the integer variable $\theta_{k}$ is relaxed to a real number. During this procedure, we aim to optimize the restricted master problem. This relaxation is effective because the restricted master problem provides a tight lower bound for the original set cover problem. More precisely, the optimal objective value $E^{\text{LP}}$ of the restricted master problem asymptotically reaches the optimal objective value $E^{\star}$ of the original set cover problem when the number of customers $N$ approaches to infinity [41, 42].

The following subsection introduces the dual problem of the restricted master problem, which is crucial to determine the existence of potential routes.

The restricted master problem (12)–(15) can be rewritten in the following standard inequality form:

Multiplying both sides of constraints (22)–(24) by the nonnegative dual variables $y_{i}$, $y_{0}$, $y_{0}^{\prime}$, and then adding the resulting terms, we obtain:

Substituting $(y_{0}-y_{0}^{\prime})$ with $y_{0}\in\mathbb{R}$, the relationship above becomes equivalent to

Assuming the coefficients on the left side of (26) are smaller than $\omega_{k}$,

we obtain the following relationship:

According to the inequality (28), the lower bound of the objective value of the restricted master problem is given by the following term:

Therefore, the dual problem of the restricted master problem (16)–(19) is

The symbol $y_{i}({}^{\forall}{v_{i}}\in V\setminus\{v_{0}\})$ denotes a nonnegative dual variable corresponding to each customer $v_{i}$. Variable $y_{0}$ is a real dual variable corresponding to the depot.

Prior to introducing the pricing subproblem, we describe the fundamental principles underlying the column generation scheme. The feasible region in the restricted master problem (16)–(19) expands proportionally with the number of generated columns. Conversely, the feasible region in the dual problem (30)–(33) contracts in proportion to the increase in the number of potential routes generated. This fact is of great importance because it implies that the dual solution $\bm{y}_{1}^{\star}$, which is optimal in the feasible region $\Omega_{1}^{\text{RST}}$, is a global optimum in $\Omega$ if $\bm{y}_{1}^{\star}$ is feasible in $\Omega$ [43]. Consequently, if the global optimum has not been determined from the current restricted master problem, there should be a concealed new route $r_{k}^{\text{new}}$ whose inequality constraints are violated by the current dual solution. This new route can be identified by solving the following pricing subproblem, which reflects the violation of the inequality constraints in the dual problem:

The objective function (34) is called the reduced cost $R_{k}$. The prior discussion asserts that route $r_{k}$ with a negative reduced cost, i.e., $R_{k}<0$ has the potential to lower the current objective value. When $R_{k}<0$, the iterative process of solving the restricted master problem, the dual problem, and the pricing subproblem continues. Figure 3 illustrates the evolution of the relevant objective values throughout the column generation process. Initially, a random set of routes is generated, the total distance of which is typically far from the optimal solution. It can be observed that the corresponding reduced cost is negative at the initial iteration. As the iteration progresses, the objective value of the restricted master problem improves, whereas the reduced cost approaches zero. If a new route is not acquired and the reduced cost is negative, the algorithm is terminated. After this iterative process, we solved the original set cover problem to obtain the optimal routes. However, these tentatively good routes may still be problematic owing to the set cover formulation; that is, some of the chosen routes visit the same customers. To resolve this undesirable situation, we finally solve the set partition problem, where the inequality constraint (13) is modified to an equality constraint. Note that the dual solutions $y_{i}^{\star}$ in the pricing subproblem determine the order of priority for each customer. This is because customers whose dual-solution value is high can lower their reduced cost if the route contains such customers.

A new route $r_{k}^{\text{new}}\in\Omega\setminus\Omega^{\text{RST}}$ is obtained as an output equivalent to the column of the matrix $A=[a_{i,k}]_{N\times|\Omega|}$. Therefore, the sequence of optimization procedures described here is called column generation. Although the pricing subproblem is the shortest path problem with negative weight edges, it entails the capacity constraints (35) and is known to be $\mathcal{NP}$-hard. One solution to this difficulty is to allow each customer to be visited more than once by modifying the matrix element $a_{i,k}$ alleviates this difficulty. Consequently, this modification enables the execution of the column generation in pseudopolynomial time using dynamic programming [42]. With respect to the restricted master problem and the dual problem, we used the COIN-OR branch and cut (CBC) solver offered by an open-source library PuLP to optimize them. An overall flowchart of the column generation method is shown in Figure 4.

In this study, we utilized a GPU-based Ising machine to solve this pricing subproblem, without modifying the matrix element $a_{i,k}$. The next section is devoted to explaining how the pricing subproblem is solved using an Ising machine.

The previous section introduced the pricing subproblem expressed in the Equations  (34) and (35). Let $q_{t,i}$ be a binary variable that takes $1$ if the vehicle visits the customer $v_{i}$ at the $t$th step, and $0$ otherwise. Using this binary variable $q_{t,i}$, the QUBO model of the pricing subproblem is as follows:

where $p$ is a penalty coefficient and defined as $p:=\max_{i\neq j}{[c_{i,j}-y^{\star}_{i}]}$ to ensure that the reduced cost term is of the same order of magnitude as for the penalty term. $T$ is the maximum number of steps that each vehicle can take and is determined to satisfy $T\geq N/U$. $z_{i}\in\{0,1\}$ in Eq. (39) is a slack variable. Equations (37) and (38) correspond to the objective function of the pricing subproblem. The terms $\mathcal{H}_{\text{c1}}$, $\mathcal{H}_{\text{c2}}$, and $\mathcal{H}_{\text{c3}}$ are equivalent to the penalty terms. $\mathcal{H}_{\text{c1}}$ stipulates that each customer should be visited at most once. At each time step $t$, the vehicle is assumed to visit only a single customer because of the term $\mathcal{H}_{\text{c2}}$. Because only one variable with a value of $1$ exists in the same row, these types of constraints are referred to as the one-hot constraints. The capacity constraints of CVRP were converted into equality constraints to adapt the penalty function method in $\mathcal{H}_{\text{c3}}$. Here, $h$ in (41) is a polynomial used to implement the conversion. Fixstars Amplify Annealing Engine selects $h$ with the least number of slack variables from the following three options:

where $h_{k}\in\{0,1\}$ denotes a slack variable. All these polynomials are used to encode integers, that is, multiple discrete values, using binary variables. Hence, these encoding methods are known as binary-integer encoding [44, 45].

The first polynomial (43) is referred to as unary encoding, and can represent integers through the summation of slack variables that take $1$. The number of slack variables required for this encoding method is equal to $(Q-1)$. Second, binary encoding, expressed as a polynomial (44), employs binary notation to encode integers. In contrast to unary encoding, binary encoding may include superfluous integers that have a detrimental effect on the inequality constraints. To circumvent this issue, the constant terms are subtracted from the polynomial (44). The number of slack variables required for binary encoding is estimated as $\mathcal{O}(\log{Q})$. Finally, the polynomial (45) is called linear encoding. Linear encoding utilizes a linear combination of slack variables. The number of slack variables required for the linear encoding is reduced to $\mathcal{O}(\sqrt{Q-1})$. Figure 5 shows the number of slack variables required for the three binary-integer encoding techniques. As shown later, we prepared five different problem settings with capacities of $26$, $60$, $122$, $169$, and $355$. Consequently, binary encoding is typically selected for the numerical simulations because it is likely to require the fewest number of slack variables.

Table I illustrates an example solution when the pricing subproblem is solved. Each row represents the $t$th vehicle movement, and each column corresponds to a depot or the customer. Note that this QUBO table satisfies the aforementioned one-hot constraints. In this case, the route is $[v_{0}\rightarrow v_{2}\rightarrow v_{3}\rightarrow v_{0}\rightarrow v_{4}\rightarrow v_{7}\rightarrow v_{6}\rightarrow v_{0}]$. This route is illustrated in Figure 6(a).

The present objective is to minimize the total distance traveled by each vehicle; thus, we require routes to comply with the following two rules:

1. 1.

   The distance matrix $[c_{i,j}]_{N\times N}$ satisfies the triangle inequality:

   $$c_{i,j}+c_{j,k}\geq c_{i,k}\quad^{\forall}v_{i},v_{j},v_{k}\in V.$$

   (46)

2. 2.

   Once the vehicle returns to the depot $v_{0}$, it cannot depart again.

With regard to these two rules, we have incorporated two postprocessing algorithms into the column generation: the $2$-opt neighborhood search and the triangular modification, which we will discuss in greater detail later. In this study, we refer to this annealing-assisted column generation as CG.

For comparison, we also applied Ising machines to CVRP naively. The general QUBO model for CVRP is relatively straightforward because we only need to extend the QUBO model (36)–(41) so that it can handle the multiple vehicles. To achieve this, the decision variable $q_{t,i}$ must be modified to $q_{t,i}^{(u)}$. Appended index $u$ denotes the $u$th vehicle. Using this newly created decision variable, we obtain:

where $p$ is the penalty coefficient defined as $p=\max_{i\neq j}{[c_{i,j}]}$. The forms of the objective function and constraints are nearly identical to those of the QUBO model written for the pricing subproblem. However, note that $\mathcal{H}_{\text{c1}}$ does not require slack variables because this QUBO model considers the routes of all the vehicles simultaneously. Hereafter in this paper, we refer to this method as AE.

We adopted the $2$-opt neighborhood search for the obtained route in order to take advantage of the rule $1)$. As the distance matrix satisfies the triangle inequality, reordering some crossing routes results in a shorter total distance for that route.

Algorithm 1 presents the pseudocode for the $2$-opt neighborhood search. The procedure for the $2$-opt neighborhood search is as follows: First, we select two edges, $edge_{1}=(i,j)$ and $edge_{2}=(k,l)$ from the obtained route. Subsequently, the selected edges are cut, the vertices are swapped, and new edges are created, namely $edge_{3}=(i,k)$ and $edge_{4}=(j,l)$. If the lengths of the newly generated edges $edge_{3}$ and $edge_{4}$ are shorter than those of the former edge, $edge_{1}$ and $edge_{2}$, we reconnect $edge_{3}$ and $edge_{4}$ to the route. This sequential process constitutes one iteration of the $2$-opt neighborhood search. This process was repeated until no further improvement was observed. The time complexity of an iteration of this neighborhood search is $\mathcal{O}(T^{2})$.

The $2$-opt neighborhood search was applied to CG when the initial set of feasible routes was generated and each time the pricing subproblems were solved. On the other hand, we applied this local search method to AE when the final solution was obtained.

The route obtained from the Table I infringes on the rule $2)$. Accordingly, we implemented a triangular modification of the solution as soon as it was generated. This postprocessing method is followed by a $2$-opt neighborhood search. The example route in Figure 6(a) includes subtour $[v_{3}\rightarrow v_{0}\rightarrow v_{4}]$ where the vehicle returns to the depot and departs from it once again. Because the distance matrix satisfies the triangle inequality, we amended this subpath to $[v_{3}\rightarrow v_{4}]$ which is certainly shorter than the previous one. The time complexity of the postprocessing method is estimated to be $\mathcal{O}(T)$. Figure 6(b) shows the modified route. Note that a $2$-opt neighborhood search was also executed.

Just as we apply the $2$-opt neighborhood search to CG, this triangular modification is called when the initial routes are generated and every time the pricing subproblems are solved. In contrast, prior to the $2$-opt neighborhood search, a triangular modification was employed in the event of the final solution acquisition regarding AE.

Because the objective function of the pricing subproblem is given by the Equation (34), the customer $v_{i}$, whose dual values are large, tends to be chosen to lower the objective value when using column generation. This potentially entraps the entire algorithm in a local optimum. Therefore, we propose fixing some variables in accordance with the following rule: customers in the $i$th iteration of the obtained route must not be visited in the $(i+1)$th iteration. More specifically speaking, let $\tilde{V}$ be the vertex set of the customers selected at the $i$th iteration. We introduce the variable fixation shown below:

This fixation allows us to explore the solution space globally because the generated columns should have little similarity to each other while simultaneously reducing the number of logical spins. The time complexity of this fixation can be estimated as $\mathcal{O}(|\tilde{V}|T)$. Most importantly, the essence of this method is to remove the overlapping of the columns that are to be generated. Consequently, this fixation method can be applied to other optimization problems with inequality constraints if it can be formulated as a set cover problem. We call this method Limited Column Generation (Limited CG). Note that we performed the $2$-opt neighborhood search and triangular modification in the same manner as for CG.

Ten instances of CVRP were solved for each problem setting, as shown in Table II. Figure 7 shows the dependence of the objective values and the number of iterations on $\sigma_{\text{d}}$. As shown in Figure 7(a), even AE could indeed efficiently solve the problem when $d_{\text{max}}$ is equal to $1$, which is equivalent to a pure VRP in this case. However, CG and Limited CG were observed to be superior to AE when $d_{\text{max}}$ increased to $10$ and $30$. When $d_{\text{max}}$ is $50$ and $100$, instance problems become excessively challenging to the extent that CG and Limited CG are unable to obtain better solutions than AE. Although the objective values for integer variables between CG and Limited CG were not significantly different from each other, we confirmed that the relaxed objective value acquired by Limited CG was superior to that of CG. This is because Limited CG overcame the issue of CG in which a feasible set of routes was rarely obtained by eliminating the overlap of routes. Figure 7(b) shows the number of iterations that both CG and Limited CG took. There was no discernible enhancement in the speed at which each method reached convergence, whereas Figure 7(b) indicates that Limited CG can achieve better results without wasting further computational resources.

Figure 8 illustrates the time evolution of the objective value for a restricted master problem. When $d_{\text{max}}=1$, the two lines almost overlap and we cannot observe the superiority of Limited CG to CG. Nevertheless, as $d_{\text{max}}$ increases, Limited CG can converge to a better lower bound at a faster rate than CG. As previously stated, obtaining a good lower bound is of paramount importance because the objective value gap between the original integer problem and its linear relaxation still exists when $N$ is finite.

When $\sigma_{\text{d}}$ was within $0.156\leq\sigma_{\text{d}}\leq 9.36$, we confirmed the superiority of annealing-assisted column generation methods (CG and Limited CG) compared to the naive method (AE) in the last subsection. Further investigations were conducted to ascertain whether these annealing-assisted column generation methods were effective in resolving these problems. In this context, we introduce another evaluation factor $\bar{N}=N/U$, which is the average number of customers per vehicle. Table III presents the problem settings for the additional simulations. Owing to the difference in the model formulations, the performance of AE should be strongly affected by $\bar{N}$. As shown in Table III, the number of slack variables $n_{\text{AE}}$ increases for a small $\bar{N}$ whereas $n_{\text{sp}}$ remains the same.

Figure 9 illustrates the dependence of the objective values on $\bar{N}$. We revealed that CG and Limited CG are effective in regions in which $\bar{N}$ is small. This result does not contradict the finding that problems with many slack variables are difficult to resolve. Conversely, when $\bar{N}$ increases, the superiority of CG and Limited CG will disappears.

We also found that this relationship also applies to when the problem size scales. Figure 10 shows the scalability of each method for two values of $\bar{N}$. Figure 10(a) shows the case where $\bar{N}=5$, while we can see the case where $\bar{N}$ is increased to $10$ from Figure 10(b). When $\bar{N}=5$ and $N=20$, the three methods were not significantly different because they had similar objective values. This result is corroborated by Equation (20), which indicates that column generation is an effective approach for large $N$. As $N$ increases, the difference in the number of slack variables becomes significant. Figure 10(b) shows the case in which $\bar{N}$ is $10$. In this case, the column generation scheme had almost no supremacy over AE because $\bar{N}$ is relatively large, and the effects of the slack variables are diminished. This outcome is also attributable to the postprocessing algorithms’ capacity to function affirmatively for AE, which has led to a significant improvement in the final path length. However, whether the ultimate result is favorable for CG and Limited CG does not highly depend on the postprocessing algorithms, but rather on the number of feasible combinations acquired until the final iteration. As seen earlier, the objective values of Limited CG are not significantly different from those of CG. However, this result indicates that the column generation scheme is powerful enough to solve larger problems under specific conditions, offering a more effective alternative to the conventional naive approach.