Exact algorithms for the chance-constrained vehicle routing problem

Abstract

We study the chance-constrained vehicle routing problem (CCVRP), a version of the vehicle routing problem (VRP) with stochastic demands, where a limit is imposed on the probability that each vehicle’s capacity is exceeded. A distinguishing feature of our proposed methodologies is that they allow correlation between random demands, whereas nearly all existing exact methods for the VRP with stochastic demands require independent demands. We first study an edge-based formulation for the CCVRP, in particular addressing the challenge of how to determine a lower bound on the number of vehicles required to serve a subset of customers. We then investigate the use of a branch-and-cut-and-price (BCP) algorithm. While BCP algorithms have been considered the state of the art in solving the deterministic VRP, few attempts have been made to extend this framework to the VRP with stochastic demands. In contrast to the deterministic VRP, we find that the pricing problem for the CCVRP problem is strongly $\mathcal {NP}$-hard, even when the routes being priced are allowed to have cycles. We therefore propose a further relaxation of the routes that enables pricing via dynamic programming. We also demonstrate how our proposed methodologies can be adapted to solve a distributionally robust CCVRP problem. Numerical results indicate that the proposed methods can solve instances of CCVRP having up to 55 vertices.

Appendices

Appendix 1: Comparison of the CCVRP with recourse models

We conducted one additional experiment to compare the solution obtained by the CCVRP with the one obtained by a recourse model for VRPSD, where the assumed recourse is that a return trip must be made to the depot whenever a vehicle’s capacity is exceeded. This experiment was done for instances having independent normal customer demands, since these are the only instances that can be solved for the recourse model. For such instances, the optimal solution for the recourse-based model and for the CCVRP were obtained, denoted as $x^{r}$ and $x^{cc}$, respectively. Both solutions were then compared according to two quantities: $z^{r}(x)$, the objective value of the recourse-based two-stage stochastic program, and $\eta (x)$, which represents the largest probability that a vehicle will have its capacity violated. Note that by design, $z^{r}(x^{r}) \le z^{r}(x^{cc})$ and that $\eta (x^{cc}) \le 0.05$—though it is possible that $\eta (x^{r}) < \eta (x^{cc})$. The instances selected for this experiment are smaller than those in the previous sections due to difficulty in solving the recourse-based model for larger instances.

Table 4 Comparison of solutions from CCVRP and recourse models for VRPSD

The results are presented in Table 4. We find that when evaluating the CCVRP solution in the recourse model, the value $z^r(x^{cc})$ was, on average, only about 1% more than the optimal value $z^r(x^{r})$, and the largest increase was 3.4%. On the other hand, while $\eta (x^{cc})$ was always (by design) less than 0.05, $\eta (x^{r})$ was greater than 0.15 in three of the instances, and as high as 0.5, meaning that in the solution there was a vehicle whose capacity would be exceeded on average 50% of the time. We thus conclude that the CCVRP model tends to yield solutions that are high quality for the recourse model, whereas the reverse is not true. In addition, the CCVRP model is not dependent on a particular assumption of the recourse taken, and can be solved also when customer demands are not independent.

Appendix 2: Instance generation details

Instances with independent normal random demands were generated by letting the mean customer demand $d_i$ be equal to the demand of the customer i in the deterministic instance. Low variance instances were generated by choosing the standard deviations $\sigma _i$ uniformly at random in the interval $[0.07*d_i, 0.13*d_i]$. For high variance instances, $\sigma _i$ was selected uniformly at random in the range $[0.14*d_i, 0.26*d_i]$.

Instances with joint normal distributions were generated by first choosing the means and standard deviations using the same procedure as for the independent normal instances. To determine correlation between two customers $i,j, i \ne j$, we first let $\gamma _{ij} =1/(\ell _{ij}*U(0.4,1.6))$, where U(0.4, 1.6) is a number chosen uniformly at random in [0.4, 1.6]. We then set the correlation between customers $i \ne j$ as $\rho _{ij} = 0.2 \gamma _{ij}/(\overline{\gamma } + \underline{\gamma })$, where $\overline{\gamma }$ and $\underline{\gamma }$ are the largest and smallest values of $\gamma _{ij}$ over all $i \ne j$. Correlations are determined in this way so that customers that are closer together tend to have higher correlation, and the scaling ensures that all correlations are significantly less than 1.0. The covariance between customers i and j is then set as $\varSigma _{ij} = \rho _{ij}\sigma _i\sigma _j$. This procedure successfully yielded a positive definite matrix in each of our test instances.

Each of the scenario distribution consists of 200 equally likely scenarios. The low variance instances were created by generating a sample of 200 scenarios from the joint normal distribution used in the low variance joint normal instances. The high variance instances were generated similarly, except that each customer i also had a probability, $p_i$, of having zero demand. These probability values were first generated randomly in such a way that about half of the customers have $p_i=0$, and the rest have $p_i$ between 0 and 0.4. Thus, to generate each scenario, the demands for all customers were first generated according to the joint normal distribution. Then, for each customer i, its demand in that scenario was set to zero with probability $p_i$. This distribution was used to provide a test with a distribution in which it is difficult to exactly calculate $\mathbb {P}\{ D(S) \le b \}$ for a subset of customers S, motivating the use of the scenario approximation. For both of the high and low variance instances, the sampled demands were rounded to the nearest integer.

Appendix 3: Computational results of primal heuristic variants

We tested the following four variants of the Clarke–Wright heuristic.

• CW: Basic Clarke-Wright heuristic.

• CW$+$LA: Clarke-Wright heuristic adapted to use look-ahead procedure to select routes to merge at each iteration.

• CW$+$LP: Clarke-Wright heuristic initialized with routes based on LP relaxation solution.

• CW$+$LP$+$LA: Heuristic initialized with routes based on LP relaxation solution and using the look-ahead procedure.

Table 5 Summary results for heuristic variants on 20 independent normal instances

Table 6 Summary results for heuristic variants on 20 joint normal instances

Table 7 Summary results for heuristic variants on 20 scenario instances

In these experiments, in the versions that initialize the heuristic based on the solution of the LP relaxation, the LP relaxation solution that is used to form the initial routes is the solution obtained after solving the root node of the BCP formulation (i.e., after the initial column and cut generation phase is complete). For the independent normal instances, the formulation based on the improved relaxation of the pricing subproblem proposed in Sect. 3.3 is used. When using the CW$+$LP$+$LA heuristic to obtain an initial feasible solution in the exact algorithms, the LP relaxation used to initialize the heuristic is the one used in the algorithm, and hence may have different results than those reported here.

Tables 5, 6, and 7 present summary results of the heuristics for the independent normal, joint normal, and scenario instances, respectively. Each table reports summary statistics over the 20 instances in that class, including the average time to solve the LP relaxation (relevant for the two versions that initialize based on the LP relaxation), the average time to run the heuristic (‘Avg. heuristic time’), the number of instances for which the heuristic found a feasible solution (‘Num success’), and the average percent by which the solution found by the heuristic exceeded the best known objective value for that instance (‘% Above best UB’), where this average is taken only over instances for which the heuristic found a feasible solution. For these experiments a time limit of 2000 s was imposed, and if a heuristic hit the time limit the best upper bound obtained is returned, and 2000 s is used as the heuristic time. The time limit was only reached for the CW$+$LA variant on the scenario distribution instances, where the limit was reached in 10 of the 20 instances. These results indicate that the methods based on the look-ahead procedure are most consistent in terms of yielding feasible solutions, and also yield the solutions of best quality. On the other hand, when initialized with just an individual customer per route, the look-ahead procedure is significantly slower than the alternatives. Initializing the look-ahead procedure with partial routes from the LP relaxation solution significantly reduces the time of the heuristic, while yielding solutions of similar quality.

Table 8 Detailed computational results for instances with independent normal distribution

Table 9 Detailed computational results for instances with joint normal distribution

Appendix 4: Detailed computational results

We present the additional computational statistics that were omitted in Sect. 5. Tables 8, 9 and 10 are the corresponding versions of Tables 1, 2 and 3, respectively. For completeness we repeat the original columns described in Sect. 5 and also present the number of branch-and-bound nodes, number of columns, and number of cuts, respectively nBB, nCols and nCuts. Column Alg. represents the variant of the exact algorithm that is being reported. Instances marked with * are ones for which the algorithm failed to find any upper bound and thus the branch-and-bound process is very ineffective, since there is no pruning by bound.

Table 10 Detailed computational results for instances with scenario distribution

Appendix 5: Results on distributionally robust test instances

Table 11 presents the detailed results for the distributionally robust (DR) experiments. We performed experiments by using the scenario-based instances and then computing the sample mean and sample covariance matrix, and using these as d and $\varSigma$ in (16). The number of vehicles was increased on these instances in order to obtain feasible solutions to the DR version of the instance. Column $\varDelta V$ represents the increase in number of vehicles relative to the original number of vehicles used for the non-DR version. With this new number of vehicles, we solved both the original scenario instance (non-DR) and the DR version of the instance, both using $\epsilon = 0.05$. For each instance and each version (DR and non-DR), we report the result from the method that yield the best performance, i.e., the method that gives the best value of feasible solution and, in case of a tie, the one giving the lowest solution time. The best method is reported in column Alg. as BC (the best branch-and-cut version for this type of instance) and/or BCP (the best branch-cut-and-price version for this type of instance). The instances marked with * are ones for which no method found a feasible solution within the time limit. The statistics reported in Table 11 are from the best method of each instance and version. We report the solution time T and the value of the best feasible solution UB and when using a time limit of 7200 s.

Table 11 Computational results on distributionally robust test instances