Solving the Vehicle Routing Problem using Cut-and-Solve

(1)

Solving the Vehicle Routing Problem using

Cut-and-Solve

December 18, 2009

Author: T.C. van Dongen, s1455478

(2)

1 Introduction

In 2006, Climer and Zhang [2] introduced a new technique, named Cut-and-Solve, to solve large Asymmetric Traveling Salesman Problem (ATSP) instances to optimality. Their unoptimized implementation outperformed state-of-the-art solvers for five out of seven real-world problem classes of the ATSP. For four of these classes, Cut-and-Solve was able to solve (substantially) larger problems. In this thesis we apply Cut-and-Solve to solve Asymmetric Vehicle Routing Problem (AVRP) instances to optimality. Moreover, we im-prove upon Cut-and-Solve by introducing a stopping rule. We show that Cut-and-Solve generally improves upon a regular ILP solver (Xpress MP) without but especially with the stopping rule.

In section 2 we provide some background information about the Traveling Salesman Prob-lem (TSP) and an extension: the Vehicle Routing ProbProb-lem (VRP). Furthermore, we elab-orate on the various mathematical formulations of the VRP and briefly explain some tech-niques that are usually used to solve the VRP and other Integer Linear Programs. Then, in section 3, we introduce the alternative to these techniques: Cut-and-Solve. We explain how it works and how we have implemented it in order to solve instances of the VRP. Moreover, we explain how we improve upon Cut-and-Solve by introducing a stopping rule. In section 4 we create 180 random instances in order to compare the computing times of the ILP solver and the computing times of Cut-and-Solve. Moreover, we test Cut-and-Solve on some well-known benchmark instances in the literature [5]. Besides this we try to find the right configuration of Cut-and-Solve by adjusting the parameters in such a way that we can find optimal solutions quicker. Finally, in section 5 we draw our conclusions based on the results.

(4)

2 Background

2.1 Traveling Salesman Problem and Vehicle Routing Problem

The Traveling Salesman Problem (TSP) is that of finding the shortest possible tour among a set of cities. Each city should be visited exactly once and all pairwise distances between the cities are known. For the symmetric case the distance from city x to city y is equal to the distance from city y to city x, but for the asymmetric case this is not necessarily true. Already in the 1800s, mathematical problems related to the TSP were treated by the Irish mathematician Sir William Rowan Hamilton and by the British mathematician Thomas Penyngton Kirkman. However, the general form of the TSP has first been studied in the 1930s by Karl Menger in Vienna and Harvard. The difficulty of this problem lies in the fact that the number of possible tours increases (super-) exponentially with the number of cities n. Assuming that all n cities are connected, there are 1₂n! possible tours for the symmetric case and n! possible tours for the asymetric case. In 1972 Richard M. Karp showed that both the asymmetric and the symmetric TSP are NP-complete problems [10].

A natural generalization of the TSP is the Vehicle Routing Problem (VRP), which was first introduced by Dantzig and Ramser in 1959 [4]. In this problem a set of K identical vehicles with capacity C are based at a depot. The vehicles should be used to satisfy demand at n collection points, such that: (1) each route starts and ends at the depot, (2) each customer is visited exactly once by exactly one vehicle, (3) the total demand of each route does not exceed C and (4) the total distance traveled by the vehicles is minimized. Note that if we set K = 1 and C = +∞, we are dealing with the TSP. Again we may deal with a symmetric or asymmetric VRP, depending on the fact whether the distances between the cities are symmetric or asymmetric. Since the VRP is a generalization of the the TSP it is also NP-complete.

(5)

the introduction of the VRP by Dantzig and Ramser, several variants and lots of methods to solve the VRP exact or approximately were proposed. For the interested reader I refer to [14] and [7].

2.2 Mathematical formulations of the VRP

2.2.1 Two-index flow formulation

In the literature several mathematical formulations for the VRP have been proposed. A popular formulation is the so-called two-index vehicle flow formulation, where O(n2_{) binary}

variables xij are used to indicate whether an arc (i, j) is crossed in the optimal solution.

Here xij will be 1 if the arc is crossed and 0 otherwise. Let K be the number of vehicles,

let cij be the cost of traversing arc (i, j), let di be the demand of city i, let V be the set of

all vertices and let r(S) be the minimum number of vehicles needed to serve set S. Then the mathematical formulation reads:

minX i∈V X j∈V cijxij s.t. (1) X i∈V xij = 1 j ∈ V \{0} (2) X j∈V xij = 1 j ∈ V \{0} (3) X i∈V xi0 = K (4) X j∈V x0j = K (5) X i /∈S X j∈S xij ≥ r(S) S ⊆ V \{0}, S 6= ∅ (6) xij ∈ {0, 1} i, j ∈ V (7)

(6)

Problem lower bound dd(S)/Ce (see [3]). The problem with constraint 6 is that it has a cardinality increasing exponentially with the number of cities n. Note that there will be 2n− 1 of these constraints. An alternative for constraint 6 can be obtained by extending the subtour elimination constraints by Miller et al. [12] to the VRP (see [13] and [11]). Constraints 8 and 9 have a cardinality increasing polynomially with n.

uj − ui+ Cxij ≤ C − dj i, j ∈ V \{0}, i 6= j, s.t. di+ dj ≤ C (8)

di ≤ ui ≤ C i ∈ V \{0} (9)

Here ui, is an additional continuous variable representing the load of a vehicle after visiting

customer i. Note that constraint 8 is not binding if xij = 0, since uj ≥ dj and ui ≤ C.

When xij = 1 the constraint reads uj ≥ ui+ dj, making sure that if an arc (i, j) is traversed

the vehicle load should increase by the demand of city j.

2.2.2 Three-index flow formulation

A natural extension of the two-index vehicle flow formulation is the three-index vehicle flow formulation that requires O(n2K) binary variables xijk, each indicating whether arc (i, j)

is crossed by vehicle k. In this formulation the vehicle capacities may differ and different costs may be associated with the vehicles. The formulation also needs an additional O(nK) variables yik indicating whether customer i is served by vehicle k.

(7)

di ≤ xijk ≤ C i, j ∈ V \{0}, k = 1, ..., K (15)

yik ∈ {0, 1} i ∈ V, k = 1, ..., K (16)

xijk ∈ {0, 1} i, j ∈ V, k = 1, ..., K (17)

Constraint 11 ensures that each city is served and constraint 12 ensures that the depot is visited K times. Here constraint 13 ensures that the in and out degrees of each vertex i both equal yi. Finally, constraints 14 and 15 are the extended versions of constraints 8 and

9 from the two-index vehicle flow formulation.

2.2.3 Two-commodity flow formulation

Recently, the two-commodity flow formulation became a more popular formulation of the VRP. This formulation was first introduced by Garvin et al. [6], who used it to solve an oil delivery problem. In this formulation, the graph G = (V, E) has to be extended with an additional copy of the depot. The depots are located at vertices 0 and n + 1 and tours are now paths from one depot to the other. Let G0 = (V0, E0) be this extended graph and let D be the set of all vertices excluding the depots. Like in the two-index flow formulation O(n2_{) binary variables x}

ij are needed, indicating whether an arc (i, j) is crossed in the

optimal solution. In addition we have two sets of flow variables yij and zji corresponding

to each arc (i, j). Now if an arc (i, j) is crossed, yij and zji are the load and residual

capacity of a vehicle such that yij + zji = C. Constraints 19 and 20 ensures that the flow

variables are adjusted according to the demand at each vertex i ∈ D. Constraint 21, 22 and 23 impose the correct values of the flow variables incident into the depot. Finally, constraint 24 ensures that load and residual capacity add up to C if an arc is traversed and 25 imposes the right vertex degree.

(8)

X j∈D y0j = d(D) (21) X j∈D zj0 = KC − d(D) (22) X j∈D zn+1,j = KC (23) yij + zji = Cxij i, j ∈ V0 (24) X j∈V0 (xij + xji) = 2 i ∈ D (25) yij ≥ 0 i, j ∈ V0 (26) zij ≥ 0 i, j ∈ V0 (27) xij ∈ {0, 1} i, j ∈ V0 (28) 2.2.4 Set-partitioning formulation

The set-partitioning formulation was originally proposed by Balinski and Quandt in 1964 [1] and requires the enumeration of all feasible circuits Hj. Now xj takes value 1 if the

circuit is used in the optimal solution and 0 otherwise. Furthermore, cj are the costs

associated with circuit Hj and aij is a coefficient with value 1 if vertex i is covered by

circuit Hj and 0 otherwise. Now we need to select K of these circuits (one for each vehicle)

(9)

2.3 Branch-and-Bound and Branch-and-Cut

In order to find exact solutions to the VRP, techniques like and-Bound and Branch-and-Cut are commonly used. We briefly discuss these techniques.

In discrete optimization problems, like the VRP, there is usually a very large number of feasible solutions which cannot all be evaluated within reasonable time limits. The chal-lenge lies in finding an optimal solution without needing to evaluate all possible solutions. Branch-and-Bound is such a technique, were the aim is to exclude large parts of the so-lution space. Without loss of generalization let us consider a minimization problem. We start by dividing the solution space into two or more subspaces. This dividing of (sub-)solution spaces is referred to as branching. Then a lower bound is computed for each of these subspaces which is referred to as bounding. This is usually done by relaxing the integrality constraints. If this lower bound is higher than the current upper bound (or current best solution), the node corresponding to this subspaces is pruned. This subspace will no longer be searched, since it will not contain an optimal solution. If the lower bound is a feasible solution this subspace will also no longer be searched since this solution will be optimal within this subspace. This process continues until all nodes are either pruned or feasible. For the interested reader we refer to [9].

A tight upper bound can sometimes be obtained quickly using a heuristic. In the liter-ature, many lower and upper bounds have been proposed sometimes specifically adapted to a particular kind of problem. Also several branching procedures are proposed in order to find the optimal solution as quickly as possible. The quality of the lower bounds and the upper bounds is crucial as these bounds are used to exclude parts of the solution space.

(10)

preserved. In this way the lower bounds at each node are improved, enlarging the proba-bility that this node is pruned. In fact we sacrifice computing time at each node in order to be able to exclude parts of the solution space in an earlier stage. Again we refer to [9] for the interested reader. An example of such cuts are the Gomory cuts [8]. We borrow an example from Climer and Zhang [2].

Example: Assume we are given three binary decision variables, x1, x2, x3and a constraint

10x1+ 16x2+ 12x3 ≤ 20. Now we relax the binary constraints by substituting the binary

constraints by 0 ≤ xi ≤ 1. It is observed that the following cut can be added to the

problem: x1 + x2 + x3 ≤ 1 without removing any of the solutions to the original

prob-lem. However, solutions would be removed from the relaxed problem (such as x1 = 0.5,

x2 = 0.25, and x3 = 0.5).

In both Branch-and-Bound and Branch-and-Cut a depth-first search strategy is commonly applied. With dept-first search one subtree of the search tree is searched entirely until a first feasible solution is found. Normally, we explore the subspace with the most promising lower bound, i.e. the lowest lower bound, as we hope that a (near) optimal solution will be in this subspace. A drawback of this approach is that we might be exploring a subspace of the Branch-and-Bound tree with no solution or a very poor solution. In this way a lot of search time is lost. Another problem is that both Branch-and-Bound and Branch-and-Cut require a lot of memory as large search trees should be stored. For larger problems this can create a lot of difficulties, especially when it is not possible to prune a lot of nodes. In the past decades many attempts were made to overcome these problems. The Cut-and-Solve algorithm overcomes these problems and will be presented in the next section.

3 The Cut-and-Solve algorithm

(11)

cut is made that cuts out a part of SS0. This part of the solution space should potentially contain the optimal solution to the optimization problem that we are trying to solve. We will refer to this part cut off by the piercing cut by SS_sparse0 . Now the best solution in SS0

sparse of the solution space should be obtained using an exact method. For this

rea-son SS_sparse0 should be small enough, such that it remains solvable within reasonable time limits. In fact we search that part of the solution space where we expect the optimal solution to be. Now the best solution in SS_sparse0 will be our first incumbent (or current best) solution. Since we already found the best solution in this part of the original solution space, there is no need for searching it again. For this reason we remove SS_sparse0 from the original solution space SS0_{. Let us denote the new reduced solution space by SS}1

(= SS0\SS0

sparse). Now we obtain a lower bound on the (new) reduced solution space

SS1_{. It is clear that if our lower bound is already worse than the best solution obtained}

so far, it makes no sense to continue the search and we have found the best solution already.

However, when this is not the case, we repeat the process. We make a second piercing cut removing a new part SS1

sparse from SS1, where we now expect the best solution to be. We

find the best solution in this part of the solution space and update our incumbent (current best) solution if we find a better solution. It is important to note that the incumbent solution should be used as an upper bound, while solving the sparse problem. This process will repeat itself until the lower bound of the reduced solution space becomes worse than the best solution so far.

(12)

is true since at each iteration more than one solution is removed from the solution space and their are finitely many solutions. It is clear that when the algorithm terminates, the current incumbent solution must be an optimal solution. This is true since it is the best solution in all solution spaces cut away by the piercing cuts and the remaining part of the solution space cannot contain a better solution, since the lower bound of this solution space is already higher than the incumbent solution.

We borrow an example from [2], where the following Integer Linear Problem (ILP) is solved with Cut-and-Solve: min y − 4 5x s.t. y +3 5x ≥ 6 5 y +13 6 x ≤ 9 y − 5 13x ≥ 1 14 x ≥ 0 x ∈ I y ≤ 3 y ∈ I

(13)

that the solution of the relaxation of the remaining part is worse than the best solution obtained, so the algorithm terminates and the optimal solution has been found.

Figure 1: Solving an IP using Cut-and-Solve.

An advantage of Cut-and-Solve is that no enormous search trees have to be stored, since only small subproblems have to be solved. After solving such a subproblem we only have to store the best solution and can throw away the search three. Besides this we search parts of the solution space that seem to be most promising. In this way we might find a good solution earlier in the process which makes it possible to speed up the process. Remember that a good solution is crucial as it is used to prune nodes.

3.1 Using reduced costs to make piercing cuts

We will use the technique proposed by Climer and Zhang [2] to make the piercing cuts. This technique makes use of the concept of reduced costs.

Reduced cost: Given a minimization problem min{cT_{x : Ax ≤ b, x ≥ 0}, the reduced}

(14)

have value zero because the corresponding coefficient in the objective is too large. The reduced cost indicate how much this coefficient has to decrease before the corresponding variable becomes non-zero in the optimal solution [9]. Alternatively, the reduced cost can be interpretated as the amount by which the objective will increase (deteriorate) if the corresponding variable is forcidbly set to one.

We start by relaxing the integrality constraints of our original problem and solve this relaxation using an LP solver. Climer and Zhang used CPLEX to do this, but we will use Xpress MP (v2.4.1) as we do not have CPLEX at our disposal. From this solution to the relaxation the reduced cost of all of our variables can easily be obtained. Now it is important to note that low reduced costs based on our relaxation make it more likely that the corresponding variable will be nonzero in the optimal solution of our original (unrelaxed) problem. Consequently, it also gives us an indication which variables will probably not be in the optimal solution, i.e. have value zero in the optimal solution. Now let S be the the set of decision variables with reduced cost lower than a predetermined α > 0 and let L be the set of all other decision variables. Now SSsparsewill be the solution

space where all decision variables in L are set equal to zero. This is easily found by adding the following constraint to the problem:

xij = 0 xij ∈ L (33)

Notice that in this way we will be looking for a solution among the decision variables that are more likely to appear in the optimal solution. Remember that we excluded all decision variables that were less likely to be positive in the optimal solution of our problem by setting them equal to zero. Basically we are looking for a solution on a very sparse graph only containing the potentially optimal edges. Moreover, since all basic variables have a reduced cost of zero, the solution to our relaxation will be in SSsparse. This is another

(15)

the reduced solution space of the next iteration will be greater or equal than the current lower bound. Note that it is important to set α neither too small nor too big, because otherwise there will either be no solution in SSsparse or it will take too much time to find

its solution. After searching the sparse solution space we add constraint 34 to the problem to obtain the (reduced) solution space. Note that we use ’≥ 1’ instead of ’> 0’ since we are looking for an integral solution. This will make sure that our current lower bound will improve with at least α at each iteration. Of course we need to drop constraints 33.

X

xij∈L

xij ≥ 1 (34)

After adding the new cut to our original problem we again solve our relaxation on the (reduced) solution space that we just obtained. Now this relaxation helps us in two ways. First of all it is the lower bound that we need to examine if we can already terminate the algorithm, but second of all it will be the source for our next piercing cut in the way we just explained.

It is crucial in our implementation that a good ILP solver is used, since it will be used at each iteration to explore the solution space cut off by the piercing cut. It would be even better to use a solving method specifically adapted for solving the AVRP on sparse graphs. Since this might well be the most time consuming step, this is a straightforward improvement.

3.2 Stopping rules

(16)

already. Does it make sense to keep making piercing cuts and search small parts of the solution space?

With Cut-and-Solve we only add one extra constraint to the original solution space at each iteration. But we do not fix any variables so the number of unfixed variables remains the same. In some sense we throw away an entire search tree at each iteration and have to start over again with the next iteration. When the optimal solution has already been found, it only remains to prove that the remaining solutions are worse than our current best solution. We noticed that especially when the initial integrality gap is large, a lot of iterations are needed to prove that the current best solution is optimal. This takes a lot more time than solving the remaining part of the solution space in one step, i.e. in one last iteration.

Hence, suppose that we have strong believe that the optimal solution has been found. From this point on it does not make sense anymore to search parts of the solution space extensively hoping to find a better solution! We believe that we will improve upon the implementation by Climer and Zhang by introducing stopping rules. If some criteria are met we will search the remaining part of the solution space in one final step, i.e. without making piercing cuts anymore. It would be interesting to see how the algorithm performs if we stop making piercing cuts if one of the following conditions hold:

• Our current best solution is only 5 percent higher than our current lower bound • We did not improve the current best solution during a subsequent iteration

Of course we could also have chosen a different percentage or a different number of non-improving iterations. However, based on our experiences in preliminary runs1 _these

condi-tions seem to be a good indication that a (near-) optimal solution has already been found. It might also be wise to adjust the second condition if we choose different values for α.

(17)

Note that it will be more likely that we will not improve in a succesive iterations if we make smaller piercing cuts.

Now suppose that the difference between our relaxation (i.e. our current lower bound) and the incumbent (i.e. our upper bound) became D. Moreover, suppose that we excluded all variables with reduced cost higher than D to make a next piercing cut. Then for sure our lower bound will increase with D at the next iteration after adding equation 34. It is important to note that the reduced cost of a variable can also be interpreted as the amount that the objective function will deteriorate if the variable is forcibly set to one in the optimal solution. Now equation 34 makes sure that at least one or a linear combination of the variables will be forced to one in the solution. Consequently, the objective function (i.e. our lower bound) will increase with at least D. Hence, we know that we can set all variables with reduced cost higher than D equal to zero before solving the remaining part of the solution space.

We would like to stress that with each iteration of Cut-and-Solve we also make a regular cut as with constraint 34 we use ’≥ 1’ instead of ’> 0’. In this way we cut out solutions that were feasible to our relaxation but are infeasible to our original problem. If we stop making iterations we do not have this advantage. Moreover, we lose our memory advantage as we still have to store a possibly enormous tree when searching the remaining solutions in one final step. On the other hand, when we apply the final step the lower and upper bound already improved a few times, which makes it easier to solve the large problem.

4 Results

4.1 A comparison on randomly generated instances

(18)

of vehicles and the excess capacity. These classes are named n-k-s, where n indicates the number of cities, k the number of vehicles and s the percentage of excess capacity. Consider for example the class 20-2-15: here we have 20 cities, 2 vehicles and total capacity is 15 percent higher than total demand.

For each of these classes we randomly created 10 instances. The demand of all cities and the pairwise distances among all cities are randomly generated from a uniform(1,100) dis-tribution. We will refer to Cut-and-Solve with parameter α with CnSα and to the Integer

Linear Problem solver from Xpress MP with ILP. This solver makes use of Branch-and-Cut. With CnSS_α we refer to Cut-and-Solve where the stopping rule is applied.

We compare ILP with CnSα on all of these 180 instances for two different formulations of

the AVRP, namely the two-index flow formulation and the two-commodity flow formula-tion. Preliminary results2 _{showed that the three-index flow formulation yields very poor}

results and the set-partitioning formulation cannot be used due to the enormous number of variables. We choose to solve all problem instances successively with ILP, CnS0.5 and

CnSS_0.5, CnS5 and CnSS5. We do this successively such that the results will not be influenced

by changing computer performance. For all instances we report the integrality gap, since we expect this to be an important factor in explaining differences in performance. In Table 1 we present the average computing times for all classes of instances and the two different formulations.

We see that CnS generally has lower computing times than the ILP for both the two-index flow and the two-commodity flow formulation. Moreover, we see that CnS5 in general

has lower computing times then CnS0.5 and that applying the stopping rule has a positive

effect on the computing times. Also note that although the average integrality gap of the two-index flow formulation is clearly larger, the total computing times for this formulation are smaller. This is mainly caused by the very poor performance in classes 20-6-5, 20-6-10

(19)

Two-index Two-commodity

Class Gap ILP CnS0.5 CnSS0.5 CnS5 CnS5S Gap ILP CnS0.5 CnSS0.5 CnS5 CnSS5

20-2-5 8.14 0.87 0.56 0.56 0.52 0.54 7.31 1.72 2.52 1.56 1.1 1.19 20-2-10 7.75 0.88 0.84 0.78 0.66 0.7 5.15 0.9 1.2 0.8 0.81 0.73 20-2-15 4.63 0.24 0.25 0.26 0.19 0.2 4.29 0.58 0.92 0.75 0.5 0.63 20-4-5 17.53 8.84 13.08 9.73 8.07 8.09 14.34 19.83 30.86 17.49 16.63 16.28 20-4-10 12.21 2.5 2.34 2.31 2.03 2.03 10.88 8.19 9.69 5.51 4.1 3.92 20-4-15 12.77 4.33 4.45 3.96 4.16 4.08 11.3 7.1 10.19 5.05 4.81 4.15 20-6-5 24.04 22.74 73.89 24.55 35.83 28.64 22.05 636.9 1527 601.13 463.5 462.79 20-6-10 17.13 3.84 14.85 8.61 7.54 6.47 15.85 72.96 215.8 105.29 121.8 107.36 20-6-15 20.82 7.71 25.18 9.31 11.06 8.35 18.54 81.35 189.9 67.82 89.22 76.95 30-3-5 9.42 41.48 14.57 14.09 17 16.5 7.57 41.12 18.84 10.9 10.89 10.7 30-3-10 4.49 4.93 2.96 2.98 1.3 1.51 3.62 4.22 3.34 2.11 1.52 1.53 30-3-15 5.33 4.33 3.08 3.15 2.1 2.4 4.57 8.29 5.42 2.54 2.6 2.58 30-5-5 11.82 245.9 124.6 146.35 157 156.92 10.04 820.3 318.1 344.11 190.8 196.86 30-5-10 10.83 87.59 30.25 29.58 51.53 50.91 8.91 107.9 52.11 33.29 29.07 28.38 30-5-15 9.78 46.78 31.93 33.38 42.82 42.79 8.81 117 46.76 30.53 30 30.86 40-4-5 7.21 568.8 169.5 166.75 146.4 141.42 6.39 1348 150.2 130.85 106.9 109.81 40-4-10 5.78 202.6 92.41 92.95 75.65 75.86 4.69 153.6 30.78 17.53 15.51 15.09 40-4-15 5.47 109.9 57.58 50.12 40.92 41.05 4.69 82.43 19.32 10.19 10.56 10.88 Total time - 13644 6624 5994 6047 5884 - 35121 26328 13875 11004 10807

Table 1: Average computing times for ILP, CnS0.5, CnSS0.5, CnS5 and CnSS5

and 20-6-15. It also looks like the two-index flow formulation especially works better on small instances and the two-commodity flow formulation works better on larger instances. An overview of the relative performance of the methods can be found in Table 2.

(20)

Method Two-index Two-commodity

A B Instances Classes Instances Classes

CnS0.5 ILP 118/180 (66%) 12/18 (67%) 85/180 (47%) 9/18 (50%) CnS5 ILP 133/180 (74%) 15/18 (83%) 135/180 (75%) 16/18 (89%) CnS5 CnS0.5 123/180 (68%) 14/18 (78%) 164/180 (91%) 18/18 (100%) CnSs_0.5 CnS0.5 120/180 (67%) 12/18 (67%) 166/180 (92%) 18/18 (100%) CnSs 5 CnS5 97/180 (54%) 9.5/18* (53%) 40/180 (22%) 11/18 (61%)

Table 2: Number of times and percentage of time that method A outperforms method B, * for one class the computing times were equal

computing time is generally the method with the lowest standard deviation. Moreover, applying the stopping rule generally decreases the standard deviation. We expected CnS to be less robust, since CnS sacrifices time to explore an area extensively hoping to find an optimal or near-optimal solution quickly. If we were searching a wrong area extensively, we might have lost a lot of time. When we increase α we expect the performance to become more robust, since the probability of finding a good solution with the (first) cuts becomes greater. Nevertheless, from the results it seems as if CnS also outperforms ILP from a robustness perspective as it has lower standard deviations. Finally, note that CnSα will

become equivalent to ILP if α becomes very large as the first cut will already contain the entire solution space.

4.2 A relation between the performance difference and the

inte-grality gap

(21)

Class Gap ILP CnS0.5 CnSS0.5 CnS5 CnS5S Gap ILP CnS0.5 CnSS0.5 CnS5 CnSS5

20-2-5 3.12 0.48 0.32 0.33 0.25 0.26 2.36 1.22 1.28 0.63 0.49 0.55 20-2-10 4.75 0.86 1 0.89 0.64 0.7 3.85 0.74 1.32 0.69 0.68 0.68 20-2-15 3.71 0.23 0.2 0.21 0.14 0.16 3.53 0.51 1.15 0.83 0.45 0.62 20-4-5 6.19 10.81 12.81 9.75 8.64 8.61 4.65 19.17 20.4 14.12 13.04 12.91 20-4-10 3.26 1.59 1.59 1.63 1.31 1.28 2.67 6.18 3.61 2.84 2.09 1.87 20-4-15 6.09 8.16 7.14 6.42 6.84 6.75 5.03 8.05 8.88 3.47 4.02 3.2 20-6-5 5.25 15.13 70.69 26.84 34.27 30.14 5.37 1040.44 2378.47 900.51 588.4 743.86 20-6-10 4.43 2.41 16.34 9.48 7.84 7.07 4.15 47.1 300.36 141.78 165.71 153.7 20-6-15 5.25 7.48 31.5 8.02 11.31 7.42 4.29 79.32 246.55 66.9 93.78 78.72 30-3-5 3.11 48.19 11.25 11.33 17.82 17.58 2.69 38.52 11.24 9.44 9.75 10.54 30-3-10 2.77 5.63 6.3 6.29 1.77 1.82 1.82 5.88 2.14 1.09 0.75 0.83 30-3-15 2.61 3.06 3.61 3.52 1.8 1.89 2.27 8.94 2.84 1.57 1.18 1.18 30-5-5 4.38 281.88 181.5 244.67 179.64 179.55 2.92 1280.92 422.25 497.45 220.73 227.08 30-5-10 3.25 145.87 34.02 32.67 95.43 93.52 2.4 109.05 48.09 36.64 45.26 44.72 30-5-15 2.5 39.51 28 30.07 45.45 45.18 1.99 79.37 29.26 20.87 26.48 27.36 40-4-5 2.01 687.41 176.8 176.48 154.3 150.19 1.5 1027.88 99.12 92.45 66.43 69.77 40-4-10 2.54 416.27 184.67 180.35 130.37 129.21 1.45 170.67 23.82 13.96 9.35 8.79 40-4-15 2.87 265.29 140.8 118.83 99.21 97.98 2.11 84.88 12.09 6.71 7.99 8.32

Table 3: Standard deviation of the computing time for ILP, CnS0.5, CnSS0.5, CnS5and CnSS5

(22)

Also remember that when the integrality gap is large a lot or iterations are needed to prove optimality of the optimal solution. We observed that this had a negative influence on the computing times. We (partly) overcome this problem by introducing stopping rules. Hence, we expect stopping rules to have effect especially when the integrality gap is large.

In order to see the relation between the integrality gap and the performance of both methods we constructed Table 4. We calculated how much of the instances belong to a certain interval for the integrality gap and what percentage of times CnS outperformed ILP for these intervals. Besides this we divided the total computing time of the ILP by that of CnS for the various intervals and formulations.

Without the stopping rule With the stopping rule Integrality gap (%) 0-7.5 7.5-15 15-22.5 22.5< 0-7.5 7.5-15 15-22.5 22.5<

Percentage of instances

Two-index 35.56 40.56 16.67 7.22 35.56 40.56 16.67 7.22

Two-commodity 45.00 38.89 13.33 2.78 45.00 38.89 13.33 2.78

Percentage of times CnS outperforms ILP

Two-index (α=0.5) 89.1 69.9 33.3 0 87.5 75.3 43.3 30.8

Two-index (α=5) 85.9 83.6 46.7 23.1 85.9 80.8 50 46.2

Two-commodity (α=0.5) 56.8 50 16.7 0 72.8 70 58.3 40

Two-commodity (α=5) 83.9 74.3 54.2 40 83.9 72.9 75 40

Total computing time of ILP divided by that of CnS

Two-index (α=0.5) 3.3 2.5 1.5 0.3 3.3 2.6 1.5 0.8

Two-index (α=5) 3.7 2.6 1.4 0.6 3.6 2.6 1.4 0.8

Two-commodity (α=0.5) 7.5 2.8 0.5 0.4 9.5 3.5 1 1.6

Two-commodity (α=5) 12 4.1 1.5 1.1 11.7 4.1 1.4 1.6

(23)

It is clear from Table 4 that performance decreases drastically with higher integrality gaps. We also see that the stopping rule especially improves the performance on instances with large integrality gaps. This matches with our reasoning as we introduced the stopping rule to overcome this problem.

4.3 The performance of the stopping rule

We clearly see from Tables 1, 2 and 4 that applying the stopping rule has a positive influ-ence on the computing times. This is especially the case if small piercing cuts are made, i.e. α = 0.5 and the integrality gap is large. This is expected as we introduced the stop-ping rule to overcome the problem of making a lot of iterations to prove that a particular solution is optimal. If we make larger piercing cuts, i.e. α = 5, the improvements are less clear. This is mainly because when we decide to solve the remaining part in one step, we are generally almost finished and setting α = 5 in the next step is almost equivalent to solving the remaining solutions in one final step.

(24)

In Table 5 we list some facts for the instances where the stopping rule is applied. With the starting solution we refer to the best solution at the time the stopping rule was applied. Similarly, with the starting gap we refer to the gap at the moment that the stopping rule is applied.

α = 0.5 α = 5 α = 0.5 α = 5

Starting solution was the optimal solution 155/180 163/180 105/180 164/180

Stopping rule is never applied due to a negative gap 20/180 38/180 9/180 18/180

Average improvement of starting solution 0.24% 0.14% 1.6% 0.2%

Average improvement of starting solution if improved 1.7% 1.5% 4.0% 2.1%

Starting time as a percentage of computing time 82.8% 91.6% 29.2% 58.8%

Average starting gap when rule is applied 5.8% 5.2% 8.3% 5.4%

Table 5: Characteristics

(25)

4.4 The choice of α

A drawback of CnS is that the choice of α is arbitrary. We have seen earlier that setting α = 5 gave the best result for both the two-index flow formulation and the two-commodity flow formulation. It would be interesting to see if there is some relation between the value of α and the performance of the algorithm.

Figure 2: Relation between computing time and alpha, black: without stopping rule, red: with stopping rule

(26)

4.5 A comparison on benchmark instances

Finally, we compared ILP and CnS on some well known benchmark problems in the lit-erature (see [5]). The largest objective coefficient is about 300 in these instances so we decided to set α = 15. For all instances we apply the stopping rule. The results can be found in Table 6.

N K % of capacity used Gap (%) ILP CnS Gap (%) ILP CnS

34 2 97.0 14.0 299.9 2578.6 9.4 120.6 56.9 36 3 81.3 8.9 958.8 610.3 8.1 161.8 48.4 39 3 93.2 5.8 38.1 60.3 4.5 38.5 11.9 45 3 89.3 6.2 50.8 34.2 5.0 358.6 60.0 48 3 76.3 6.1 88.9 55.0 5.7 237.3 48.4 56 3 82.1 11.0 M E 9.7 4291.4 916.7 65 3 94.2 7.1 8703.1 2848.2 5.9 5600.3 581.9 71 3 80.2 8.0 277.9 275.5 7.3 4971.7 3810.9

Table 6: Computing times for some well-known benchmark instances [5], M: not possible to compute due to memory problems , E: runs into a fatal error without explanation

(27)

5 Conclusions

CnS generally outperformed ILP on both the randomly created test instances and the benchmark problems. The computing time and the standard deviation of the computing time was generally lower. It became clear that the performance of CnS had a relation with the integrality gap. This is the case for two reasons. First of all, we used the integrality relaxation to do local searches hoping to find a good solution quickly. When the integral-ity gap is large, the probabilintegral-ity that we were doing local searches in wrong areas becomes larger and consequently we may have lost a lot of time. Second of all, we noted that CnS needed a lot of time to prove that a solution was optimal, mainly because it had to do a lot of iterations to close the remaining gap. After each of these iteration we add a cut, but all variables remained unfixed. It appeared that solving the remaining part of the solution space at once is a lot quicker. Because of this we introduced a stopping rule: if we believed that we already found the optimal or a near-optimal solution we solved the remaining part of the solution space at once. This only gave big improvements if the remaining gap was rather large and the piercing cuts were rather small. Otherwise making another piercing cut was almost equivalent to solving the remaining solutions at once. This was confirmed by the computing times that we found when applying CnS with the stopping rule on the randomly generated instances.

(28)

We believe that further research can be done to make piercing cuts more effectively. For example: one could think of ways to gradually increase α as the lower bound and upper bound become better. Note that better bounds make it easier to solve larger problems. In this way it might be possible to solve larger problems without running into memory problems and without needing to make too many unnecessary iterations. One could also think of different stopping rules. For example by starting the final step if the incumbent solution did not improve enough or increasing the necessary starting gap. It would also be interesting to see how CnS performs if a combinatorial relaxation is used. Finally, it would be very interesting to see how CnS performs on different (Binary) Integer Problems. For the ATSP and the AVRP it already showed a promising performance.

Acknowledgements: Special thanks to Prof. Daniele de Vigo for providing the benchmark instances.

References

[1] M. Balinski and R. Quandt. On an integer program for a delivery problem. Operations Research, 12:300–304, 1964.

[2] Sharlee Climer and Weixiong Zhang. Cut-and-solve: an iterative search strategy for combinatorial optimization problems. Artif. Intell., 170(8):714–738, 2006.

[3] G. Cornuejols and F. Harche. Polyhedral study of the capacitated vehicle routing problem. Mathematical Programming, 60:21–52, 1993.

[4] G. B. Dantzig and J. H. Ramser. The truck dispatching problem. Management Science, 6(1):80–91, 1959.

(29)

[6] W.M. Garvin, H.W. Crandall, J.B. John, and R.A. Spellman. Applications of linear programming in the oil industry. Management Science, 3:407–430, 1957.

[7] Bruce Golden, S. Raghavan, and Edward Wasil. The Vehicle Routing Problem: Latest Advances and New Challenges. Operations Research/Computer Science Interfaces Series , Vol. 43. Springer, 2008.

[8] R. E. Gomory. Outline of an algorithm for integer solutions to linear programs. Bulletin of the American Mathematical Society, 64:275–278, 1958.

[9] F. S. Hillier and G. J. Lieberman. Introduction to operations research. McGraw-Hill Science/Engineering/Math, 2005.

[10] R. M. Karp. Reducibility among combinatorial problems. In R. E. Miller and J. W. Thatcher, editors, Complexity of Computer Computations, pages 85–103. Plenum Press, 1972.

[11] R.V. Kulkarni and P.V. Bhave. Integer programming formulations of vehicle routing problems. European Journal of Operations Research, 20:58–67, 1985.

[12] C.E. Miller, A.W. Tucker, and R.A. Zemlin. Integer programming formulations and traveling salesman problems. Journal of the ACM, 7:326–329, 1960.

[13] P. Toth N. Christofides, A. Mingozzi. Combinatorial Optimization, chapter The vehicle routing problem, page 315338. Wiley, Chichester, 1979.

Solving the Vehicle Routing Problem using Cut-and-Solve