Abstract
This paper presents a flexible solution methodology for the capacitated vehicle routing problem with stochastic travel times (CVRPSTT). One of the basic ideas of the methodology is to consider a vehicle working time lower than the actual maximum vehicle working time when designing CVRPSTT solutions. In this way, the working time surplus can be used to cope with unexpected congestions when necessary. Another important idea is to transform the CVRPSTT instance to a limited set of capacitated vehicle routing problems (CVRP), each of which is defined by a given percentage of the maximum vehicle working time. Thus, our approach can take advantage of any efficient heuristic that already exists for the CVRP. Based on the two key ideas, this paper presents a simulation-based algorithm, in which Monte Carlo simulation is used to obtain estimates of the cost and the reliability of each solution, and the Clarke and Wright heuristic is improved to generate more reliable solutions. Finally, a number of numerical experiments are done in the paper with the purpose of analyzing the efficiency of the described methodology under different uncertainty scenarios.
1. Introduction
Congestion is a common phenomenon in most urban areas of the world, and it can be caused by a variety of factors, such as accidents, traffic conditions, and weather conditions. In urban logistics, congestion creates substantial fluctuations in travel times and hence affects carriers’ cost structure and the relative weight of wages and overtime expenses [1]. Ignoring congestion and the travel time fluctuations when developing route plans for pick-up and/or delivery vehicles can result in inefficient and suboptimal solutions. Therefore, the travel time among customers and depot is found to be a crucial factor of vehicle routing in urban logistics.
Despite a great amount of research on the vehicle routing problem (VRP), considerable research has been devoted to the general problem by using constant values to represent the travel times. Research on the problem with variant travel times is comparatively sparse. To our knowledge, there are mainly two types of research that are related to the problem with variant travel times in the existing literature. These research work are discussed in the following.
Considering time-dependent travel times in solving vehicle routing problems can reduce the costs of ignoring the changing environment to some extent [2]. The time-dependent VRP (TDVRP) was first formulated by Malandraki and Daskin [3, 4], who proposed dividing the time horizon in slices and assigning a constant travel time for each arc and slice . Therefore, the travel time of each arc is a piecewise step function [5]. Most studies on TDVRP have paid attention to the FIFO property when developing models and algorithms for the problem. And it is only in recent years that a number of heuristics had emerged [6–9] that are enabled to deal with large-scale problems.
However, the TDVRP is not completely correspondent with the actual congestions as a result of the following two points despite the fact that the planning gap (difference between plan and actual) has been reduced by the research of TDVRP. First, a number of researches in the transportation science have revealed that many factors can lead to congestions, such as length and width of road, number of crossings, tunnels, and bridges, accident possibilities, weather conditions, and even seasonal shopping [10]; congestion levels are never the same from day to day on the same route because the variety of traffic-influencing events that influence congestion are never the same [10], and hence, travel times on the same route are usually different even at the same period of a day [10–12]. Second, there are usually some alternative routes between two places in the actual transportation network, each of which has its own distribution of travel times. Therefore, the travel time distribution between two places should be a composite function of the travel time distribution of each alternative route. A recent study by Garaix et al. [13] showed that not considering alternative paths can be disadvantageous in many situations.
The above findings undermine the hypothesis of the studies of TDVRP because the travel time for each arc and slice should be a variable, not a constant, even at the same time of a day. The study of Kenyon and Morton [14], from their experiments of small stochastic vehicle routing problems (SVRP) with only nine nodes, indicates that solutions to the stochastic model, in which the travel time , for each arc and slice , is considered as a stochastic variable, can be significantly better than solutions obtained by solving the associated mean-value model (i.e., the deterministic vehicle routing problem in which all random parameters are replaced with their population means). This can be explained from a simple fact that given a variable for the travel time of edge , there is little probability for a solution to contain the edge as a result of a large variance of . However, if is regarded as a constant, the characteristic of , such as its variance and its probability density, will be certainly neglected. Furthermore, a great deal of research in the transportation science field has founded their theory on a stochastic time-dependent transportation network [15], which also indicates the importance of considering the stochasticity of travel times.
Studies of the vehicle routing problems with stochastic travel times (VRPSTT), in which the travel time between every node pair is a stochastic variable, can make up the above deficiencies of TDVRP. However, the VRPSTT is arguably one of the most challenging and practical variants of the VRP [16], and in the literature, only a few studies addressed the VRPSTT. Laporte et al. [17] proposed a variant of VRPSTT chance-constrained programming model and compensation model, where a penalty function was introduced to prorate the delay time. Park and Song [18] subsequently constructed three new heuristic algorithms based on an extension of VRP algorithms. Kenyon and Morton [14] developed two stochastic programming models. Lecluyse et al. [19] introduced the variability in traffic flow into the model, which was used to evaluate the routes based on the uncertainty involved. Connors and Sumalee [20] and Chen and Zhou [21] studied the stochasticity of travel times from the view of travelers’ equilibrium. Zhang et al. [16] studied the stochastic travel-time vehicle routing problem with simultaneous pick-ups and deliveries and developed a new scatter search approach for the problem by incorporating a new chance-constrained programming method. Furthermore, an efficient way to generate multivariate random variables that are interrelated is presented in [22], which gives a better way to look at the randomness of the travel times from a global level, rather than redistricting it to each local route. And papers [23–25] show a novel approach, that is, the joint diagonalization strategy, for solving stochastic algebraic equations, which is an alternative to solve the VRPSTT.
Although there have been a number of studies on VRPSTT in the literature, practical algorithms and solution approaches are still needed for such problem which is a very challenging and difficult combinatorial optimization problem due to the stochasticity of travel times.
To solve the CVRPSTT, we present a novel simulation-based algorithm. Different from other methods on VRPSTT, this algorithm solves the problem by “reducing” a complex CVRPSTT, where no efficient metaheuristics have been developed yet, to a limited set of more tractable CVRPs where excellent, fast, and extensively tested metaheuristics exist. Furthermore, this algorithm is valid for any travel time statistical distribution with a known mean, and cannot only solve CVRPSTT instances with hundreds of nodes efficiently but also provide decision makers with various solutions to practical problems.
The remainder of this paper is arranged as follows. First, the CVRPSTT is described and analyzed in the next section. Section 3 describes the main framework and the details of the solution procedures. The experimental settings and the results are presented in Section 4, and the paper finishes with the conclusions.
2. Problem Description
Let be a graph, as exemplified in Figure 1, where is the set of customer nodes and depot node and is the set of edges. Each node is associated with a positive demand except . Each edge is associated with a positive random travel time . Define as a set of identical vehicles available to deliver goods from the depot to the customers. Each vehicle has the same capacity and the same maximum vehicle working time which can be also defined as the depot closing time in a day given that the depot opening time is 0. Some additional constraints associated with the problem are as follows:
The route of vehicle , , which is defined as where and represent the depot and the other vertices represent customers, is a cycle starting and ending at the depot and serving a given subset of sequential customers. Let if the edge (, ) is included in and 0 otherwise. Define as the completion time that vehicle completes all its delivery tasks and returns to the depot given its route in. It is called “deficient route” for the route of vehicle if . In a deficient route, a penalty cost has to be incurred.
The CVRPSTT studied in this paper is a single-objective problem in which the total cost, including the travel cost of vehicles and the overtime wages, should be minimised. The objective function for CVRPSTT can be defined as formula (1), where and are used to transform the travel time and the vehicle’s lag time of returning to depot to travel cost and overtime wage, respectively. Generally, /> 1, implying a penalty of overtime wage, and hence, formula (1) can be naturally simplified to formula (2), where / indicating a penalty coefficient when vehicles arrive at depot later than their closing time. And the value of can be much larger than 1, for example, 10, 20. Consider
Formula (2) can be also understood as the total cost when . Thus, the total cost contains two components, the travel cost defined as and the penalty cost defined as . Let if and 0 otherwise. Then the reliability of a solution can be defined as , indicating the percentage of the total number of deficient routes in a solution, and it is an important index in analyzing a solution.
3. The Simulation-Based Algorithm
There are two points that make it difficult to solve the CVRPSTT. First, it is challenging to integrate the stochasticity of travel times into an algorithm. Note that in formula (2) is a random value which can be different from day to day even in the same time of a day. And it is not reasonable to represent it with a single value in solving the problem. Second, how to balance the travel cost and the penalty cost is another important factor in constructing a solution. In general, a solution with less travel cost usually has more penalty cost. This is mainly because a route in a solution with less travel cost is often fully engaged in delivery tasks, and hence, there is often not much flexible time left to cope with unexpected transportation congestions. And once the unexpected event occurs, the penalty cost will probably be incurred. Therefore, an algorithm that is devoted to minimizing the total cost should be able to integrate the stochasticity of travel times and generate a solution that can balance well the travel cost and the penalty cost and hence minimize the total cost.
To solve the CVRPSTT, we present a new simulation-based algorithm (SBA), which can cope well with the above difficulties. The key idea behind our algorithm is to transform the issue of solving CVRPSTT into a new issue which consists of solving multiple “conservative” CVRPs, each characterized by a specific risk (probability) of suffering deficient routes. The term “conservative” refers here to the fact that only a certain percentage of the maximum vehicle working time will be considered when solving CVRPs. In other words, part of the maximum vehicle working time will be reserved for attending possible “emergencies” caused by underestimated random congestions during the actual distribution (routing execution) process. The motivation of transforming CVRPSTT to CVRP comes from the facts that: (a) reducing the value of generally means that more time can be reserved for a vehicle to cope with the possible congestion, and thus, the risk of returning late to the depot will simultaneously be reduced; (b) but on the other hand, if too much time is reserved for vehicles, the travel cost will increase because long routes that have exceeded the specified percentage of the maximum vehicle working time have to be split; (c) there will be a certain percentage of the maximum vehicle working time that can minimize the total cost defined by formula (1); and (d) our algorithm is devoted to finding the optimal percentage by transforming CVRPSTT to multiple CVRPs and solving every CVRP.
The main framework of SBA is described in Algorithm 1.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
The specific details of our algorithm are explained in the following.(1)The random travel time for edge can follow any known statistical distribution, either theoretical or empirical, as long as its expected value exists. Our algorithm does not assume that all travel times must necessarily follow one designated distribution so that different edges can follow different distributions. Furthermore, since our algorithm uses simulations, it could also consider dependences among different travel times.(2)The loop variable represents the percentage of that is considered for solving CVRP. In other words, when solving CVRP(r), the maximum vehicle working time is set to , and represents the reserved time for a vehicle to cope with travel time delays due to congestions.(3)The CVRP(r) can be solved by using any efficient CVRP methodology, for example, the classic Clarke and Wright heuristic [26] which is used in our experiments because the heuristic is able to generate a very good solution in short time and is thus suitable for large-scale simulation.(4)The total cost of solution OS obtained by solving CVRP(r) in step 5 does not contain any penalty cost, and this is an assumed situation that the travel time is equal to the expected travel time for edge , and no congestions or unexpected events exist during the distribution process. However, congestions are inevitable in actual distributions. Therefore, it is necessary to simulate the actual travel times in order to determine the expected cost for a solution.(5)Steps 7–11 use Monte Carlo simulation to generate random travel time for each edge according to its distribution function and calculate the total cost, the travel cost, the penalty cost, and the reliability of OS based on the disrupted travel time value of any edge in OS. After iterating this process for some hundred/thousand times, a random sample of observations regarding these variable travel time are obtained, and an estimate for the expected total cost of OS can be calculated in step 10 by the following expression: where represents the total cost of th time simulation. And estimates for the expected travel cost, the expected penalty cost, and reliability of OS can be simultaneously calculated in step 10 by similar expressions.(6)In step 8, is generated as the travel time for an edge in a simulation and is used to substitute its expected travel time which is employed in solving its corresponding CVRP in step 5. Given a reasonable travel time distribution function, it is obvious that should not be negative or less than or more than values that violate common sense.(7)The solution CS obtained in step 14 can be improved because we have found from a large number of observations that some routes in the obtained solution CS are fully engaged in deliveries while others are relatively idle and able to undertake more tasks in their spare time. The reason that produces the result of the imbalance is that the solution CS is generated by an existing heuristic for CVRP which aims at minimizing the travel cost but not penalty cost of a solution. In other words, CVRP does not consider penalty cost at all. As a result, the fully engaged route is prone to return late to depot and hence results in a higher penalty cost once any unexpected congestion occurs.
To solve the problem, we develop a method for step 15 to balance the workloads among multiple vehicle routes so that the penalty cost and the associated total cost of solution CS could be reduced simultaneously. The key idea behind the method is to iterate the process of deleting a customer node from a relatively engaged route, inserting it into a relatively idle route and evaluating the new solution by using Monte Carlo simulation until the total cost cannot be reduced. The method in step 15 is described in Algorithm 2.
| ||||||||||||||||||||||||||||||||||||||||
Note that step 14 in Algorithm 2 is the most time-consuming step whose efficiency should be improved. We implement the improvement by storing the expected total cost for each route of solution IS, calculating the expected total costs for only two changed routes er and lr in step 15, and adding up the expected total costs of all the changed routes and unchanged routes as the total cost of solution IS. This way is proved from our experiments to greatly accelerate the simulation process of step 14.
4. Computational Results
In this section, we present some computational results of the SBA, which has been coded in Java and run on a laptop with an Intel Core 2 Duo CPU T7100 at 1.8 GHz, 2 GB RAM, and the Microsoft Windows XP operating system.
4.1. Test Instances and Parameter Setting
As mentioned by Li et al. [27], there are no commonly used benchmarks in the SVRP literature, and therefore, each paper presents a different set of randomly generated instances for the different SVRP variants that they studied [14, 16, 17, 27]. This situation makes it difficult to compare the performance of different approaches. As a result, we generalize a well-known set of 59 classical CVRP instances which includes a diversity of clustered and disperse problems of different sizes by using random travel times instead of constant ones and compare the results of our algorithm with the best-known CVRP solutions reported in [28]. The details of these instances can also be found here [28], and the distribution and its parameters used to randomly generate travel times are given below so that other authors can use the same data sets and distribution for verifying and benchmarking purposes.
Since the travel times in VRPSTT can follow any distribution and our heuristic does not have any limitation as long as its mean value is known, we choose the Log-Normal distribution to model the travel time for each edge as a result of its nonnegative values. The two parameters of the distribution, the location parameter and the scale parameter , are formulated by the following two expressions, respectively:
For each instance, we changed the travel time of an edge (, , ) in CVRP instances to a random value by setting to the Euclidean distance between the two nodes according to their coordinates. The original CVRP instances can be particularly defined with Var[] . But for CVRPSTT instances, we have the inequation Var[] > 0. A larger variance of an edge indicates a higher uncertainty when travelling the edge. Referencing the idea from Juan et al. [29], our heuristic is tested under three different variance levels, the relatively low variances, the medium variances, and the relatively high variances. And these variance levels are further illustrated in Table 1, in which the relationship between the variance and the mean value for each variance level is given.
In order to decide the parameter values of our algorithm, preliminary experiments were performed on some instances with different parameter configurations to determine which configuration would be effective in solving the CVRPSTT. The preliminary experiments showed a satisfactory performance when the parameter values were set to , , , , , and (the average route length of each instance’s best-known solution).
4.2. Detailed Results for One Instance
In this section, we present the numerical results for only one of the 59 instances in order to better illustrate our algorithm. Given the classical A-n54-k7 instance, we generalized it by considering the travel time for edge (, and ) as a Log-Normal random variable with E[] = the Euclidean distance and Var[] = 0.5 (medium-variance scenario). Next, we ran our algorithm for the A-n54-k7 instance 21 () times in which different runs have different and values. But in one run we set in order to make the algorithm generate a solution for only one percentage of the maximum working time. Table 2 shows alternative solutions (for the A-n54-k7 instance with medium variances) that we obtained for different values of and . For each of the obtained solutions, the following additional information is provided: number of routes employed, travel costs, penalty costs, total costs, and reliabilities.
In order to better demonstrate the character of alternative solutions to the A-n54-k7 instance, graphical representations are given in Figure 2 to illustrate the changing trends of travel costs, penalty costs, total costs, and reliabilities for all values of . Definitions and solving procedures of these four indexes are given in Sections 2 and 3. From Figure 2, the following five points can be observed and discussed.
4.2.1. Discussions on Travel Costs
Overall, the index of travel costs shows a first decreased and then steady trend. The overall decreased trend when the value is between 0.8 and 0.9 is because with the growth of the value, a vehicle is able to complete more tasks, and hence, some good routes with more tasks can be adopted for a solution. And one of the main reasons to explain the overall steady trend when the value is between 0.9 and 1 is that the major factor that limits the task number of a route becomes the vehicle capacity, not the maximum working time, when the value is bigger than 0.9. In other words, vehicles that satisfy the constraint of the maximum working time do not necessarily meet the capacity requirement, but vehicles that meet the capacity requirement usually satisfy the constraint of the maximum working time. And the situation is maintained for this instance when the value is between 0.8 and 1.
However, there are some exceptions in the overall trend, such as a small increase when the value is between 0.8 and 0.9. This is consistent with common sense under the situation of random travel times.
4.2.2. Discussions on Penalty Costs
There are two factors that influence the index of penalty costs: the total travel distance of all vehicles and the time span reserved to cope with congestions. The longer the travel distance, the greater the possibility of returning late to depot and hence the greater the penalty cost. The more the reserved time, the smaller the possibility of returning late to depot and hence the smaller the penalty cost. If the travel distance is long and simultaneously much time has been reserved to cope with congestions, whether the penalty cost increases or decreases depends on the factor which plays a greater role. When the value is less than 0.87, the penalty cost is gradually reduced with the growth of the value, which indicates that the travel distance factor plays a greater role. When the value is greater than 0.87, the penalty cost gradually increases with the growth of the value, which indicates that the reserved time factor plays a greater role.
4.2.3. Discussions on Total Costs
The figure of the total costs shows a first decreased and then increased trend. This is an index of the weighted sum of travel cost and penalty cost according to formula (2) in Section 2. Note that when is 0.86, neither the travel cost of 1244.77 nor the penalty cost of 20.54 is the minimum, but the total cost reaches the minimum value of 1450.13.
4.2.4. Discussions on Reliabilities
The two indexes of both penalty costs and reliabilities can illustrate the level of a solution’s probability of suffering deficient routes. However, the index of reliabilities shows a first increased and then decreased trend, which is just the opposite to the overall trend of penalty costs. This is because the index of reliabilities denotes the percentage of the total number of deficient routes while the penalty cost is the sum of penalty costs of deficient routes. But note that the maximum reliability (0.912) and the minimum penalty cost (20.35) do not simultaneously appear in the same value as a result of the difference in their expressions.
Relative to the penalty cost index, a solution’s reliability provides decision makers with another index to evaluate the reliability of the solution. If a decision maker cares about the length of (), he can consider the index of penalty cost; otherwise, if a decision maker wants to minimize the number of deficient routes, the index of reliabilities may come in handy.
4.2.5. The Running Time of Our Algorithm
As for the running time of our algorithm, the process of generating the results for one alternative solution to the A-n54-k7 instance can be completed in just three seconds (average 3062 milliseconds in our laptop whose configuration is relatively low). Even for an instance with 135 nodes (instance F-n135-k7), the process of generating the results for one alternative solution can be completed in about 19 seconds. Therefore, it is completely affordable to implement the process for most real cases.
4.3. Results for 59 Instances
Table 3 shows the average results obtained for all 59 classical instances that we generalized and tested under three different scenarios. The meaning of each column is as follows.(i)BKS expected total cost: the expected total cost of the best-known solution (BKS), whose computing process is similar to SBA (Algorithm 1) except the following differences: the BKS computing process does not include step 15; meanwhile, steps 4 and 5 are replaced by directly reading the best-known solution instead of generating a solution for CVRP(r), and and are set equal to ensure the loop (steps 2–13) can be run only once.(ii)SBAS expected total cost: the expected total cost of the solution generated by the SBA (SBAS) described in Algorithm 1.(iii)# Routes: number of routes.(iv)Gap: the gaps in expected total costs between BKS and SBAS, calculated as , where BKS and SBAS stand for the values of their expected total costs.
The table corresponds to three uncertainty scenarios (low-variance, medium-variance, and high-variance) which were previously described. As discussed in this section, similar ideas and conclusions to the ones reached for the A-n54-k7 instance with medium variance also apply for other instances and uncertainty levels.
From the above tables, the following four points can be observed and discussed.(1)In a low-variability scenario defined by , results show that the expected total costs of SBAS tend to provide a better average reliability level of 0.976 when applied to the CVRPSTT. As is discussed before, lower reliability levels imply larger number of deficient routes and thus higher expected variable costs. This explains the average gap of 1.25% between expected total costs of BKS and SBAS, the former containing the total expected costs associated with the BKS when demands are stochastic instead of deterministic. Notice also that the mean number of routes for the chosen solutions has slightly increased from 7.44 in BKS to 8.25 in SBAS, reserving some extra time tends to increase the number of necessary routes, but the resulting solution tends to be more reliable or robust. Finally, in this low-variability scenario the average value for the recommended reserving time level is ; that is, each vehicle will reserve up to 6.6% of the maximum vehicle working time to attend unexpected congestions.(2)In the case of the medium-variability scenario with Var[] 0.5[]2, the BKS is less reliable than before (the average reliability level is about 0.855), and therefore, their associated expected total costs tend to be higher in this new scenario characterized by a higher degree of uncertainty. And the average gap in expected total costs between BKS and SBAS is about 4.19%. Meanwhile, our methodology is able to provide more reliable solutions (approximated average reliability of 0.887). Other interesting results in this second scenario are the average number of routes (8.44), which is slightly higher than the case of the low-variability (8.25), and the average value (0.915), which is lower than the case of the low-variability (0.934).(3)The high-variability scenario reports an estimated average reliability of just 0.692 for the BKS and a corresponding average gap of 3.73% in expected total costs between BKS and SBAS. At the same time, our methodology generates solutions with an estimated reliability index of 0.740. Also, the average number of routes in our solutions is 8.46, which is slightly higher than in previous scenarios, and the average value is 0.902, which is slightly smaller than in previous scenarios.(4)It is interesting that all average reliabilities of SBAS in Table 3 are higher than the average reliabilities of BKS. Figure 3 graphically summarizes the average reliability indices for each uncertainty level and solution (BKS versus SBAS).
5. Conclusions
In this paper, we present a simulation-based algorithm for the capacitated vehicle routing problem with stochastic travel times (CVRPSTT). One of the basic ideas of our methodology is to consider a vehicle working time lower than the actual maximum vehicle working time when designing CVRPSTT solutions. In this way, the working time surplus can be used to cope with unexpected congestions when necessary. Another important idea is to transform the CVRPSTT instance to a limited set of CVRPs, each of which is defined by a given percentage of the maximum vehicle working time. Finally, a number of numerical experiments are done in the paper with the purpose of analyzing the efficiency of the described methodology under different uncertainty scenarios.
In our view, the main contributions of this paper include the following.(1)The methodology is virtually valid for any statistical distribution with a known mean, either theoretical, for example, Normal, Log-Normal, Weibull, Gamma, and so forth, or experimental, in which historical data could be used to model each edge’s travel time. Also, the methodology can be naturally extended to consider different distributions for different edges and possible dependences among these travel times. The methodology is thus developed because in the practical situations travel time between any pair of nodes can be different from day to day even at the same time of a day and there may be some correlations among travel times of multiple routes, which, to our knowledge, cannot be easily dealt with by the existing methodologies.(2)Our methodology can be used to solve CVRPSTT instances with hundreds of nodes in a reasonable time by “reducing” a complex CVRPSTT, where no efficient metaheuristics have been developed yet, to a limited set of more tractable CVRPs where excellent, fast, and extensively tested metaheuristics exist.(3)Moreover, our methodology can provide decision makers with various solutions to practical problems, each of which considers a different percentage of the actual maximum vehicle working time, indicating different levels of risk or reliability. This makes the methodology more flexible in order to satisfy different decision makers.
Conflict of Interests
The authors declared that they have no conflict of interests.
Acknowledgments
This work is partially supported by the Grants from the National Natural Science Foundation of China (nos. 71271037, 71171029, and 71272093), the Foundation of Liaoning Educational Committee (no. L2012019), and the Fundamental Research Funds for the Central Universities (no. DUT12JR09). Those supports are gratefully acknowledged.