A Robust VRPHTW Model with Travel Time Uncertainty

1 Introduction

Vehicle routing problem (VRP) is one of the classical problems in operational research. There are many uncertain parameters in VRP model such as the accident, weather conditions, vehicle number on road, road characteristic factors and so on. These uncertain parameters have a large extent effect on the optimal solution of VRP, especially on the optimal solution of vehicle routing problem with hard time window (VRPHTW) that is extremely strict in travel time. It is an important issue to take travel time uncertainty into consideration in VRPHTW. Robust optimization is one of methods that can solve such kind of problem.

Robust optimization is an optimization method based on robust control theory. Robustness means that control system is insensitive to small changes. The solution obtained by robust optimization can adapt to the change of parameters, that is, robust optimization ensures the feasibility and effectiveness of the solution. Robust optimization is first proposed by Soyster^[1] in the seventies last century. He considers the uncertainty of parameters and presents completely robust optimization that would cause a waste of resources. Mulvey et al.^[2] developed a general model formulation of robust optimization and applied robust optimization into mathematical model containing uncertain information. Bertsimas et al.^[3] applied robust optimization to obtain solution of objective programming containing uncertain value and proposed a robust method with adjustable robustness. They pointed out that when the number of uncertain variables is less than or equal to the degree of robustness that is set by decision maker, the model can be solved easily through robust dual transformation, and the optimal solution of the original model has better robustness; When the number of uncertain variables is greater than the degree of robustness, the optimal solution of original model can still be effective and optimal with high probability.

Robust optimization is concerned by various sectors scholars. Here are some of the endless examples to show the robust optimization’s widely use. Tian et al.^[4] built robust supply chain network optimization model and adjusted the degree of robustness by setting the regret value. Li et al.^[5] used robust optimization in the assembly line production where operating time of production line is an uncertain quantity. Wang et al.^[6] altered robust optimization method instead of traditional optimization method to solve fleet operating problem.

Robust optimization is not applied widely on VRP but until recent years, the application of this method on solving the uncertain vehicle routing problem is gradually increasing. In the study of VRP, robust method is mainly used in the customer demand uncertainty and travel time uncertainty. Sungur et al.^[7] used robust optimization method to minimize the cost of VRP where the demand is a dynamic parameter. Calculation results show that although optimal solution increased a certain cost, the solution can satisfy the uncertain demand. Gounaris et al.^[8] studied the resource constraints vehicle routing problem (CVRP) with uncertain demand. Their robust CVRP model can adapt to the changes in demand and can reduce distribution costs. Han et al.^[9] solved the vehicle routing problem with soft time windows and uncertain travel time. Their method does not need to determine the probability distribution of random variables and still can obtain optimal path with minimum expected cost.

Bertsimas robust method is a classical robust optimization method. Many scholars follow their steps. Zhao et al.^[10] applied the robust method to ATIS route guidance. In general vehicle routing problem the objective function is mainly to minimize the transport cost. However, their objective function of route guidance is to minimize the time. They used Bertsimas robust method to handle travel time uncertainty in view of the path risk. The optimal solution is more stable under the condition of the uncertain path information. Lee et al.^[11] considered the travel time uncertainty and demand uncertainty in VRPTW model and propose label method to solve shortest path problem with resource limits. Solomon numerical example results show that their method can obtain robust solution of model with uncertain factors more accurately. But their model only considers the deadline and travel costs except the earliest time window and other vehicle cost. Agra et al.^[12] built a maritime transport model and use Bertsimas robustness method to solve the travel time uncertain problem. Hierarchical method of graph theory and Xpress Optimizer are also applied into the model to solve the robust vehicle routing problem. Then on the basis of the previous work Agra et al.^[13] put forward two kinds of formulation method to process parameter uncertainties. Numerical results show that their method solving uncertain problem is as easy as solving a certainty problem by general approach. But Bertsimas robustness method is not obviously used in their article since their key point is labeling method.

From the above literature review, it is found that robust method is one of methods suitable for the optimization of problems with uncertainty and can make the VRPTW model solution more stable. But the current study based on robust optimization only considered the deadline or soft time window. It is an important issue to explore the VRP model with hard time window based on robust optimization. In addition, study should perfect their objective function in model to give thorough consideration of cost such as transport cost, fixed vehicle cost, and waiting cost. Motivated by the above observation, this paper formulate a VRPTW model with uncertain travel time using Bertsimas robustness method. The objective of this model is to minimize the total cost of vehicle traveling cost, fixed cost, and waiting cost. We also present a modified Max-Min Ant System Algorithm (MMAS) algorithm proposed by Stützle and Hoos^[14] to solve a numerical example.

The rest of the paper is organized as follows. Section 2 outlines our VRPHTW model with uncertain travel time. In Section 3, we analyze the VRPHTW model. Algorithm and example are presented in Section 4. Finally, Section 5 provides conclusions and points out directions for future research. All the proofs of the theoretical results are given in Appendix A.

2 VRPHTW model with uncertain travel time

In this section, we develop the VRPHTW model with uncertain travel time.

Suppose there is one distribution center (hereafter referred as DC) and n demand points denoted by the vertex set V of one undirected graph G = {V, E}, where V = {0,1,2, · · · , n}, 0 denotes the DC, and the others denote the demand points. In addition, E = {e_ij}is the set of arcs of G, where e_ij denote the path between i and j, i ≠ j. Demand in point i is deterministic and constant denoted by q_i, i = 1, 2, · · · ,n. DC has m vehicles denoted by K = {k|k = 1,2, ··· ,m}. The DC dispatches them to serve every demand point (hereafter referred as customer). Each vehicle has a maximum load denoted by M_k and an average travel speed denoted by v_k. Every vehicle starts off from the DC and finally returns to the DC. Once a vehicle is dispatched, a fixed cost f_k is incurred. Every customer is served by one and only one vehicle. The travel time is the sum of uncertain travel time and the average travel time. The average travel time is determined by the average travel speed and the path length. For the uncertain travel time, there are a lower bound and an upper bound both of which can be obtained from the historical observations. The uncertain travel time distributes in such an interval randomly. Vehicles should arrive at the demand point within the customer time window. If the vehicle arrives late, the demand of this customer would disappear. If the vehicle arrives in advance, waiting cost will be incurred until the customer accepts the service. For the sake of conciseness, service time is not involved in this model as it can be transferred as the upper bound of the time window. The objective of the DC is to minimize the total cost including transport cost, waiting cost and fixed cost by deciding the distribution quantity and the distribution roote for every vehicle.

The following is the other notation used in the model.

Among the above variables, x_ijk is decision variable and the others are exogenous variables. Then we can obtain the VRPHTW model as follows.

The objective function is to minimize the total cost of transport cost, fixed cost and waiting cost. Constraint (1) is the vehicle capacity. Constraint (2) ensures that every demand point is served by one and only one vehicle. Constraint (3) is on the balance of the vehicle stream. Constraint (4) and (5) ensure that vehicles should start from DC and end at DC. Constraint (6) guarantees that there is no loop between demand points. Constraint (7) are the travel time of vehicle k on arc e_ij that contains the uncertain variables η_ij. Constraint (8) ensures that each customer is served before the upper bound of the time window.

3 Robust VRPHTW model analysis

In robust linear programming, a small disturbance of some variables will not affect the optimal solution. The main idea of Bertsimas robust optimization is to obtain an optimal solution that when the number of uncertain variables changes within the predefined scope, the optimal solution is still a feasible or an optimal one; when the uncertain variables change exceeding predefined scope, the optimal solution can be a feasible or optimal solution under a computable bounded probability. Define Γ as a parameter that stands for the uncertainty set by decision maker, i.e., Γ is the degree of robustness or the degree of conservatism. It should be noted that Γ may be not an integer and 0 ≤ Γ ≤ n(n + 1), where n(n + 1) is the total number of travel path between every two points in V considering direction. Since the travel time of some paths may be deterministic in reality, it can be assumed that the travel time of ⌊Γ⌋ paths in η_ij are uncertain and the travel time of only one path in Γ−Γηi′j is uncertainty. Thus, there are ⌊Γ⌋ + 1 paths at most the travel time of which are uncertainty. Based the above idea of robustness, the time window constraints of (P1) in Section 2 can be formulated as (9), (10) and (11).

where S is an arc set contains ⌊Γ⌋ arcs.

Expression (9) can be interpreted as that first rank all possible travel paths e_ij ∈ E according to their biggest uncertain time i.e. η̅_ij at the same time and then choose the first ⌊Γ⌋ + 1 paths to take their travel time uncertainty into calculation in the form of δ_ij. The rest uncertain time variations are smaller than the first ⌊Γ⌋ + 1 ones and the solution of robust VRPHTW model can be feasible when Γ fetches values from an appropriate scope. The original VRPHTW model is transformed into a robust VRPHTW model as following.

If Γ equals to zero, it is a deterministic model without considering travel time uncertainty. If Γ equals to n(n + 1), it means that the travel time uncertainty of all the pathes are taken into account and thus the model is equivalent to the complete robust model of Soyster^[1], i.e. the solution remains feasible for every possible realization of travel time. That is also the reason Γ is regarded as the degree of conservatism. The bigger Γ is, the more conservative the decision maker is. As Γ fetches other values, solution of this model can be adapted to travel time changes in different effects along with the degree of robustness Γ. There is a trend that the robustness of solution will increase as the value of Γ increases. However it is not necessary for Γ to tend too close to n(n + 1) since it can reach in consideration of resource saving as complete robust solution is always a waste of resource.

Then how can we handle the realization of travel times are out of initial prediction since there are more than ⌊Γ⌋ + 1 paths with uncertain travel time? As the arriving time can be expressed as a robust form, we can amplify of upper bound of the probability of time window constraint violation. Theorem 1^[3] shows that when the number of uncertain time variations is larger than ⌊Γ⌋ + 1, the optimal solution of robust VRPHTW model remains feasible at a certain probability.

Theorem 1

Suppose that the probability of time window constraint violation is p(t̃_jk > b_j), that is, the real arriving time of vehicle k arriving at j is later than the latest time window of j. There is an upper bound of the probability of time window constraint violation as following.

p(t~jk>bj)≤exp⁡(−Γ22n(n+1))(12)

It can be found that the upper bound is independent of the solution and only depends on the number of customers and the degree of robustness. Denote exp⁡(−Γ22n(n+1)) by y(Γ). It is obvious that y(Γ)is decreasing in Γ. For n =100, there are 10100 pathes at most and the curve of y(Γ)is shown in Figure 1. It can be found that the probability of time window constraint violation decreases quickly in Γ. As Γ equals to 422, y(422)≈ 0.001 which is a very considerable data compared with 10100. Thus the above model (P2) can be used to solve the original VRPHTW model (P1) since we can obtain a wonderful upper bound of the probability of time window constraint violation.

Figure 1

The curve of y(Γ) when n = 100

4 The modified ant colony algorithm and numerical analysis

In this section a modified ant colony algorithm is proposed to solve model (P2) and a numerical example is conducted to illustrate the theory analysis and the algorithm.

Ant colony algorithm is suitable for this kind of problem. In ant colony algorithm, ants cooperate by exchanging pheromone similarly as max-min ant system algorithm (MMAS). The main idea of MMAS algorithm is to set upper and lower bound of pheromone to avoid premature convergence. Here we use a modified ant colony algorithm to solve such model. The key part of our algorithm is to identify the set Λ that is a contradictory objective with cost objective. The uncertainty of travel time is reflected by random variation of distance between every two demand points. Different random case has different vehicle routing plan. But Λ should be independent of the random case since the decision maker don’t know time variation accurately in advance. In our algorithm, the busy road sections that has high utilization frequency are used to identify Λ initially. This is a robust idea to protect path solution against the worst case. Busy road sections could be obtained from the historic data.

Algorithm 1:

1. Initialization.

   Maximum iteration s = Nmax, the number of the optimal route num = Num, demand points number n, the feasible solution number of each iteration is W, the degree of conservatism is Γ. Input cost function, demand coordinates, demand quantity, time windows, vehicle speed, maximum load of each vehicle, pheromone matrix and other initial value of ant colony algorithm, etc.

2. Find set Λ.

3. Generate random demand point distance to represent the uncertainty of travel time distance = distance(0.5+rand(1)). Generate initial route solution by simple exclusive method algorithm and calculate initial route cost, record it as historical best solution. s = 0, num = 0.

4. s = s + 1, the feasible solution number of each iteration is w = 0.

5. w = w + 1, put m ants on the distribution center and the serial number of each ant is j = 0.

6. Look for feasible demand point where the vehicle load constraints and time windows are all satisfied. If there are some feasible demand points to move on, calculate transition probabilities and use roulette wheel selection method to decide which is the next point to move. then n = n – 1 and keep a record of visited demand point. Turn to step 2.5. If there is no feasible demand point, j = j + 1 and start step 2.4 again.

7. Determines whether moved ants have traversed all demand points. If n ≠ ∅, turn to step 2.4. If n = ∅, output one feasible solution and if w < W, turn to step 2.3, else, turn to step 2.6.

8. Calculate W routes’ cost, the smallest one is the optimal vehicle route of the s time iteration. Compare it with historical best solution and update historical best solution. Use super excellent ant colony optimization strategy to update the ant pheromone. Control the concentration of pheromone on the base of Max-Min Ant algorithm. If s = Nmax, turn to step 2.7, else, turn to step 2.2.

9. Note the historical best solution as global optimal and num = num + 1. If num = Num, turn to step 2.8, else, turn to step 2.1.

10. Count the total number of each road section that Num optimal routes contain and the statistical result is the busy road section set, that is, the set Λ.

11. Find a robust solution in a random situation.

12. Generate random demand point distance. Find the top Γ busiest road sections from set Λ and reset the demand point distance by roust theory. Generate initial route solution and calculate initial route cost. s = 0.

13. Same with step 2.2.

14. Same with step 2.3.

15. Same with step 2.4.

16. Same with step 2.5.

17. Same with step 2.6, but if s = Nmax, turn to step 3.7, else, turn to step 3.2.

18. Output the historical best solution which is the robust solution of this problem.

Then we conduct the following numerical example transforming from R101 in Solomon benchmark examples to illustrate our model and algorithm. There are 20 demand points and their coordinates are in the first and second column. DEMAND stands for quantity demanded. READY and DUE represent ready time and due time of demand point, respectively. SERVICE is the service time at demand point. Other initial parameters are as follows. Ant colony size is 20. Maximum iteration is 5000. Transport cost, fixed vehicle cost, and waiting cost are 1, 2 and 0.1, respectively. Persistence of pheromone is 0.2. Weight of pheromone and visibility of path are 1 and 5. Weight of pheromone and time window used in calculating transition probability are 0.6 and 0.4. The basic data of this numerical example can be seen in Table 1.

In the numerical example, decision maker can identify a busy road sections set after obtaining 200 path solutions in 200 different situations. Record the frequency of each road section as shown in Table 2. The number 115 in the second row and the first column mean that there are 115 vehicles pass through road section e₁₀ in 200 path solutions. Then Λ can be determined by selecting Γ or ⌊Γ⌋ + 1 busy road sections. For all road sections that belongs to Λ, we record the largest time variation when Γ is an integer. If Γ isn’t an integer, we record (Γ – ⌊Γ⌊) folds of the last road section’s largest time variation. The time variation of road sections that does not belong to Λ is without any consideration. Thus the optimal path solution can be obtained by programming. The constraint violation probability of each solution is computed from ten thousand times random experiment.

In Figure 2, the points are optimal solutions violation probability calculated from random selection of Λ instead of busy road sections shown with solid line. The drop-down dot dash curve is the upper bound of time window constraint violation probability changed with Γ in theory. It can be found that the time window constraint violation probabilities of each optimal solution calculated from busy road section set are all below their upper bound curve. However, the solutions calculated from random Λ set are not when Γ is larger than 64 which means that it is an efficient way to select Λ by busy road section set.

Figure 2

Constraint violation probability calculated by different Λ with n = 20

In Figure 3, the dotted line is the total cost modified by dividing 700 in order to draw with the probability together. If Γ = 0, the robust model degenerate into general VRPHTW model and the solution constraint violation probability is relatively very high from the figure. As there are only 20 points in graph, the route choice is so limited that the downtrend of optimal solution violation probability is not so obvious. It can still be found that the violation probability get smaller when Γ increases, that is, the robustness of solution gets stronger as Γ increasing, while it cost more that can be found from dotted curve. In other words, strong robustness can be achieved with high cost which is coincident with the theoretical analysis.

Figure 3

The constraint violation probability and the total cost with n = 20

The numerical example illustrates that an optimal solution can be obtained as long as the decision maker gives a reasonable conservatism degree Γ even when the travel time is uncertain. Furthermore, our robust model can balance cost and constraint violation probability well.

5 Conclusion

In this paper, the travel time uncertainty is taken into consideration in a VRPHTW model. In addition, transport cost, fixed cost and waiting cost are also included in the total cost. Based on Bertsimas’s robust optimization idea, robust optimization is used in solving such VRPHTW model. The probability of time window constraint violation bounds is derived. From the trend of probability of time window constraint violation, it is argued that if the decision maker has a reasonable estimate of robustness degree Γ, the solution of robust VRPHTW model would not violate time window constraints. We also propose a modified ant colony algorithm and a numerical example is presented. The results of numerical experiment and the theoretical analysis illustrate that the time window constraint violation probabilities of every optimal solution are all smaller than its upper limit.

For the future research, it is worthwhile to extend our method to explore vehicle routing problem with soft time window.