Improved branch-cut-and-price for capacitated vehicle routing

Abstract

The best performing exact algorithms for the capacitated vehicle routing problem developed in the last 10 years are based in the combination of cut and column generation. Some authors only used cuts expressed over the variables of the original formulation, in order to keep the pricing subproblem relatively easy. Other authors could reduce the duality gaps by also using a restricted number of cuts over the master LP variables, stopping when the pricing becomes prohibitively hard. A particularly effective family of such cuts are the subset row cuts. This work introduces a technique for greatly reducing the impact on the pricing of these cuts, thus allowing much more cuts to be added. The newly proposed branch-cut-and-price algorithm also incorporates and combines for the first time (often in an improved way) several elements found in previous works, like route enumeration and strong branching. All the instances used for benchmarking exact algorithms, with up to 199 customers, were solved to optimality. Moreover, some larger instances with up to 360 customers, only considered before by heuristic methods, were solved too.

Appendix 1

We report detailed computational results for a large set of CVRP instances, comparing the proposed BCP with some of the recent exact methods described in Sect. 1. Additional computational results include a comparison of five variants of the pricing presented in Sect. 4.1. Moreover, some techniques used in the BCP to better handle special instances with a high degree of symmetry are presented.

1.1 Detailed computational results

The BCP was run over the standard classes of instances (A, B, E, F, M, and P, available at www.branchandcut.org) used for testing exact methods for the CVRP. The name of an instance from this benchmark follows the general format X-nY-kZ, where X denotes the class it belongs to, Y corresponds to the number of vertices and Z refers to the fixed number of routes (K) required for any feasible solution. As usual in the literature of exact methods, the cost matrix for these instances is calculated from depot and customers coordinates, following the TSPLIB convention of rounding the euclidean distances to the nearest integer. Moreover, in order to better compare our results with those presented in [3], our BCP uses for each of these instances the same initial upper bound they used.

We also report results over instances traditionally used in the literature on heuristic methods, proposed in [5, 6] (instances with initial letter “C”) and [14] (instances with initial letter “G”). In those cases, the convention is not rounding the Euclidean distances and not fixing the number of routes. The initial upper bound used for each of these instances corresponds to the best solution value available in the literature.

Table 3 Detailed results on CVRP instances of class A

Table 4 Detailed results on CVRP instances of class B

Table 5 Detailed results on CVRP instances of classes E and M

Table 6 Detailed results on CVRP instances of class F

Table 7 Detailed results on CVRP instances of class P

Tables 3, 4, 5, 6 and 7 present detailed results for all instances from the datasets A, B, E-M, F and P, respectively. For each instance, column Ins gives the name of the instance, UB presents the initial upper bound used and OPT shows the value of the optimal solution found (values in bold indicate proven optimality for the first time). The following 10 columns are root node information. The first of these columns indicates the lower bound obtained using only robust cuts (RLB) and the next four columns report information immediately before enumeration: the improved lower bound due the addition of non-robust lm-SRCs (NRLB), the number of active lm-SRCs (i.e., those with non-zero dual variables) in the master LP (nS), the percentage of arcs $(i,j)^q \in A_Q$ fixed to zero (Fix) and the accumulated time up to this point ($\texttt {T1}$). Then, the number of routes enumerated (the word “limit” indicates that the enumeration was unsuccessful), (|R|i), the improved root node lower bound after route enumeration, obtained by adding additional non-robust SRCs and clique cuts (NLB), the number of remaining enumerated routes after that (|R|f), the time spent with route enumeration plus the time to go from lower bound NRLB to lower bound NLB (or only the time spent with enumeration, if it is aborted) (T2) and, finally, the time spent solving the resulting MIP (TMIP) complete the information of the root. Columns |R|i to TMIP contain values only for the instances where a successful enumeration is performed at the root node. Column Nds gives the number of nodes in the branch-and-bound tree and TT gives the total time in seconds. Additionally, there are four other columns: BMR11, Con12, Rop12 and CM14. They are, respectively, the time in seconds reported by [3, 8, 32], and [9]. For each method, the number of solved instances, average time and geometric mean time are also reported.

Table 8 Detailed results on instances from [5] and [6]

Table 9 Detailed results on instances from [14]

We remark that the Ref. [32] only presents the aggregated results shown in Table 1. The detailed results for Røpke’s method shown in Tables 3, 4, 5, 6 and 7 were kindly provided to us in 2012 by personal communication. In fact, as one of the first tasks of the project that culminated with BCP proposed in this article, we replicated his method and obtained similar results in each instance. Therefore, we can attest that those results are trustworthy. This is an important concern, since they did not appeared yet in a peer reviewed paper.

Table 8 shows detailed results for selected instances proposed by [5] and [6] and Table 9 presents detailed results for instances with up to 360 customers proposed by [14]. Those instances are the most relevant for comparing recent heuristic approaches. For each instance, two new columns indicate the number of customers and the number of routes in the optimal solution found. These tables do not include columns of competing exact methods because no such results are reported in the literature.

Fig. 3

Optimal solution of M-n200-k16, value 1274

Fig. 4

Optimal solution of G14, value 1080.55

Fig. 5

Optimal solution of G19, value 1365.60

Figures 3, 4 and 5 illustrate the optimal solution found for instances M-n200-k16, G14 and G19.

Table 10 shows detailed results over a number of selected instances of five variants of the pricing: 1—forward labeling (Algorithm 3); 2—bidirectional labeling with concatenation at Q / 2; 3—bidirectional labeling with dynamic determination of the concatenation point (Algorithm 5); 4—the later algorithm with mild use of completion bounds ($n_S'$ is small); 5—the same algorithm with aggressive use of completion bounds ($n_S'$ is larger). For each instance, there are five rows, each giving the results for a pricing variant. The name of the columns are: Ins—the name of the instance; UB—the initial upper bound; LB—the root lower bound obtained; Fix—the percentage of arcs $(i,j)^q \in A_Q$ fixed at the end of the root node; Alg—the pricing variant; nS—the number of lm-SRCs with non-zero duals; nFL—the number of forward labels; FT—the time to run the forward labeling; nBL—the number of backward labels; BT—the time to run the backward labeling; CT—the time to run the concatenation; CBT—the time to calculate completion bounds; TT - the total time of the pricing; and qf—the concatenation point. The number of customers is indicated in parentheses for the instances proposed in [14]. Note that for variant 1, $\mathbf{qf} = Q$.

The results (columns from nS to qf) are an average of the rounds of the exact pricing performed in the end of the root node, after the last separation of lm-SRCs. In order to ensure a fair comparison, in this particular experiment, all the pricing variants are run “in parallel” with exactly the same input, the reduced costs. We only use columns found by variant 5, the columns from the other variants are discarded.

The results in Table 10 show that variant 5 is the best pricing algorithm, at least for very hard instances, those with more than 199 customers. Therefore, it was the variant chosen for the proposed BCP.

1.2 Handling symmetrical instances

A high level of symmetry is observed in most instances from [14]. This impacts negatively the BCP performance, making the cutting and branching operations less efficient than usual. In fact, whenever a fractional solution is cut, it is likely that a symmetric solution with the same cost will appear in its place. The lower bounds can only move after all symmetric solutions are cut. We implemented some special techniques in our BCP in order to mitigate this negative impact.

The proposed BCP has a procedure that automatically detects, for instances that use non-rounded Euclidean distances calculated from the customer and depot coordinates, angles of rotation and reflection with respect to the depot that result in equivalent instances. For example, the instance G14 (see Fig. 4) remains the same if it is rotated around the depot in an angle of 90 degrees or reflected over a horizontal line passing through the depot. As a result, for every feasible solution for this instance (fractional or integral), there are eight symmetric solutions with the same cost. For the instance G19, shown in Fig. 5, the number of symmetric solutions is 12 because the rotation angle is 60 degrees. More precisely, the procedure uses the rotation and reflection operations mentioned before to partition the set of arcs into symmetry groups. In instances G9–G12, each group has two members, in instances G13–G19, each group has eight members, in instances G17–G20, the groups have size 12. Instances G1–G8 are also very symmetric, but they are not considered in this CVRP work because they have additional distance constraints.

Let $Sym(\phi )$ be the set of all elements of a given instance I that correspond to the element $\phi$ in one of the equivalent instances that belong to symmetry group of I, an element being a vertex, arc, route or solution. The symmetry detection procedure works as follows. First, it seeks for rotation symmetries. For an arbitrarily chosen customer vertex i, it tests every other vertex j having $d_j = d_i$ and $c_{0j} = c_{0i}$ as candidates to belong to Sym(i). For each candidate, it checks whether or not the corresponding rotation angle with respect to the depot can be applied to every vertex leading to another vertex with the same demand value. In the positive case, the rotation angle is accepted. After all candidates are checked, the smallest accepted rotation angle is used. Then, the procedure seeks for a straight line passing through the depot that defines a reflection symmetry operation. The used algorithm is similar but, since any vertex can be exactly over the reflection line all pairs of vertices are tested as candidates. For each candidate, the corresponding reflection line is applied to all vertices, being accepted if every vertex that do not belong to the line has a corresponding symmetric. In this case the procedure stops upon the first accepted reflection line. Although the time to detect symmetries can be potentially significant, we observed that it is negligible (less than one hundredth of second) for all instances from [14].

Table 10 Results for 5 variants of the pricing

The BCP benefits from the symmetry detection in two ways.

• 1. For every variable $x_a^q$ fixed to zero by reduced cost, all variables $x_{a'}^q$ where $a' \in Sym(a)$ can also be fixed to zero. This is valid because the same argument that was used to fix a can be turned into an argument to fix $a'$, by only switching indices.

• 2. Consider the simple disjunction $(x_a = 1) \vee (x_a = 0)$ on the branching. The stronger disjunction $(x_a = 1) \vee (\sum _{a \in Sym(a)} x_a = 0)$ is also valid. This can be explained as follows: if all optimal solutions use at least one arc in Sym(a), there is an optimal solution that uses a. It is important to remark that the left node (corresponding to $x_a = 1$) from this branching is not completely symmetrical anymore. This is good, but also means that the symmetry related tricks can not be used in the corresponding subtree. On the other hand, the right node (corresponding to $\sum _{a \in Sym(a)} x_a = 0$) is still symmetrical. Now, consider a disjunction $(\sum _{a \in \delta ^-(S)} x_a = 1) \vee (\sum _{a \in \delta ^-(S)} x_a \ge 2)$ on the branching. Let X be a non-empty proper subset of the arcs in A. We say that the subset $X' \subset A, |X'|=|X|$, belongs to Sym(X) if every arc in X can be mapped into an arc in $X'$ by the same combination of rotation and reflection operations. It can be seen that all possible sets Sym(X) have the same cardinality |Sym|, which is also the cardinality of the groups of symmetrical arcs. We claim that the right side of disjunction can be strengthened to $(\wedge _{X' \in Sym(\delta ^-(S))} \sum _{a \in X'} x_a \ge 2)$. In other words, |Sym| inequalities can be imposed on the right node. This is valid because if in all optimal solutions there is a set $X' \in Sym(\delta ^-(S))$ such that $\sum _{a \in X'} x_a = 1$, there is an optimal solution where $\sum _{a \in \delta ^-(S)} x_a = 1$. Using this improved disjunction, the symmetry is partially broken in the left child but is maintained in the right child, and thus can be exploited again.

The use of those techniques makes a lot a difference in the BCP performance on highly symmetrical instances. For example, instance G17 takes more than 3 days to be solved without them.