Two-stage Particle Swarm Optimization Algorithm for the Time Dependent Alternative Vehicle Routing Problem

Consider a delivery problem faced by a distribution center or cargo company in urban environments. Traveling conditions change throughout the day, and the problems associated with congestion are most obvious during peak-hours. When this is the case of the delivery of the goods, allowing the drivers to take alternative routes will let them easier to meet the requirement of the delivery time by the customers.


Introduction
According to a report published by the International Energy Agency in 2009, approximately 23% of global carbon dioxide emissions can be traced back to transportation [1]. Reducing the total time spent on traveling and the number of vehicles assigned at distribution centers would be one of the methods to decrease carbon dioxide emissions caused by transportation.
Consider a delivery problem faced by a distribution center or cargo company in urban environments. Traveling conditions change throughout the day, and the problems associated with congestion are most obvious during peak-hours. When this is the case of the delivery of the goods, allowing the drivers to take alternative routes will let them easier to meet the requirement of the delivery time by the customers.
Several researchers have studied the Time Dependent Vehicle Routing Problem (TDVRP) which takes the issues of time-dependent travel time into account. Different from the other VRP variants, the TDVRP is to determine a set of vehicle routes originating and terminating at a single depot when the travel time of vehicles depends on the time of day in order to minimize the total travel cost. All customers are visited within their claimed time windows exactly once, and the total demand of customers assigned to each route does not violate the capacity of the vehicle. TDVRP was first proposed by Malandraki et al., [2] who employed discrete step functions to describe travel time along different time periods of a day; and formulate the problem into a Mixed Integer Programming model. Although this kind of approximation is easier to be integrated into a mathematical programming model; it fails to satisfy the principle of First In-First Out (FIFO). Then, Ichoua et al. [3] defined the continuous travel time functions in the TDVRP and showed their appropriateness and ability to satisfy FIFO. This paved a way to the subsequent studies [4][5][6][7]. Most of these studies derived the travel time functions from the travel speed functions through integration or various mathematical transformations, which in turn, were often represented as step functions.
Although TDVRPs have been discussed intensively, few studies have reported on the existence of alternative routes in the TDVRP. Due to its realistic and useful applications, this study intends to investigate the properties of such problem entitled by Time Dependent Alternative VRP (TDAVRP). Due to the NP-hard nature of the problem, the concept of Particle Swarm Optimization (PSO) is adopted in order to develop an efficient algorithm for solutions. Sensitivity analysis will also be conducted.
The rest of the paper is organized as follows: Section 4 introduces the definitions and mathematical model of the proposed Time Dependent Alternative Vehicle Routing Problem (TDAVRP). Section 5 describes the proposed PSO algorithm for solving the TDAVRP. Section 6 discusses the computational experiments using the proposed PSO on Solomon's benchmark [8]. After the summary of this work in Section 7, Section 8 concludes the results of this research and suggests the directions for future studies.

The Proposed Model
The Time Dependent Alternative Vehicle Routing Problem (TDAVRP) can be represented by a complete undirected multi-graph denoted by G(N, E), where N is the set of nodes, Ν={1, 2, …, n} with additional node 0 as a start depot and node n+1 as a terminal depot; E is a set of edges connecting pairs of nodes by E= E 1 ∪ E 2 with E 1 , the edge set of designated edges and E 2 , the edge set of alternative edges.
V is the set of vehicles, V ={1, 2, …, v max } where v max is the maximum number of the available vehicles. M is the number of intervals of travel time functions considered for each designated edge, M = {1, 2, …, m max } and m max is the maximum index of time interval.
Before modeling, some basic assumptions of the TDAVRP are stated as below: (1) One and only one vehicle serves each customer at the allowed time window to satisfy the demand; any violation will cause additional cost as a penalty.
(2) All routes must originate and end at the depot, that is node 0=node n+1.
(3) Only consider two possible edges between two nodes, one is called the designated edge, the other one is the alternative edge as shown in Figure 1a.
(4) The travel speed distribution of the designated edge is time dependent and is defined by a step function; that of the alternative edge is a constant as presented in Figure 1b. Based on the stated problem, a Mixed Integer Programming model is proposed as below of which the notations are defined in Appendix A: The objective function (1) shows that this model minimizes the total route cost, consisting of fixed transportation cost (Z 1 ), time cost (Z 2 ) and distance cost (Z 3 ). Decision makers can assign different weights ( σ α ) to decide the importance of each cost in the planning of routes. Constraints from (2) to (6) are path constraints. Constraint (2) shows that no more than one designated edge or one alternative edge are selected for the transportation plan, while allowing that neither edge be selected. Specifically, we allow fewer than v max vehicles to be used. Constraints (3.1) and (3.2) ensure that each customer is served by exactly one vehicle. Constraints (4) and (5) are the out-degree and flow balance constraints ensuring that the solutions are consistent with the set of routes. Constraint (6) eliminates routes including edges from the depot directly back to itself. Constraint (7) is the capacity constraint limiting the total demand on a route to the capacity of the vehicle. Constraints (8.1) and (8.2) are time window constraints, indicating that the departure time from any customer node should fall within the departure time windows; otherwise penalty costs would be incurred. Constraints (9) to (17) are continuity constraints on travel time. Constraint (9) is the relation inequality of m ijv t and iv l . Constraints

Two-stage Particle Swarm Optimization
Particle Swarm Optimization (PSO) is easy to implement, and the parameters which are needed to be adjusted are less than Genetic Algorithms (GA) [9]. PSO's fast convergence [7] and good performances on the continuous search space have drawn our attention, and is adopted in this study.
In this section, the concept of the design of a heuristic algorithm based on PSO is first described in Section 5.1. Then, the procedure of implementing PSO is presented in Section 5.2.

Design Concept of the Proposed PSO
In the solution to the TDAVRP, there are two main variables that must first be established: the service order for each vehicle; the departure time at each node and the selection of edge between each pair of nodes. It is difficult to consider these variables simultaneously, because the departure time of each node is based on the sequence of customers.
To facilitate a solution, this study proposes a two-stage PSO, including the Primary and Secondary PSOs. The Primary PSO is used to determine the service order for each vehicle and the Secondary PSO provides information relates to the travel strategies such as departure time and edge selection. Local improvement is also included to improve service sequences derived from the Primary PSO.
Primary PSO-Solution of Service Sequence: This subsection describes the basic concepts of the Primary PSO, such as the solution representation method, the construction of an initial population, and the decoding method.

st Solution representation
For the Primary PSO, the solution representation method proposed by Ai et al. [10] is adopted to describe particles.
The encoding is divided into two parts. Assume n be the number of customers to be served and v max is the number of available vehicles. The total dimension of a particle representation is n+2v max where n  represents the customer priorities and 2v max denotes the vehicle route orientation respectively.

nd Initial population
The initial population takes the advantage of the information related to the input data of the TDAVRP by : (1) Ranking the time windows of customers in ascending order If two customers have the same starting time windows, we compared the spans of their time windows, and set the shorter window (higher priority) ahead of the longer. If two customers have the same starting time windows and spans, we arrange the order randomly.
To quantize the concept of customer priority for PSO operations, a set of n random numbers is generated from the range of [0, 1]. Then, the smaller random number is assigned to the customers with higher priority.
(2) Randomly generating two real numbers for each vehicle One number is generated from the range of [0, px max ] and the other is from the range of [0, py max ]. The values of px max and py max depend on the maximum values of customer coordinates in the second part of the solution representation in the Primary PSO.

rd Decoding method
To determine the dimension of the particles in the Secondary PSO, we must first decode the particles in the Primary PSO. In accordance with the decoding method developed by Ai et al. (2009), particle positions are decoded into v max or less than v max vehicle routes using the following procedure: (Ai et al., 2009) Step 1 Construct the priority list of customers Step 2 Construct the vehicle priority matrix Step 3 Route construction Local Improvement-Improvement of Solution Quality: Based on the procedure of the Primary PSO, a set of vehicle routes is established. Because time factors are not taken into account during the decoding procedure in the Primary PSO, situations involving waiting or lateness at customer nodes occur frequently in the Secondary PSO. To avoid the serious impact of time-window factor to the final solution, a local improvement method is proposed.
In this phase, the remaining information for each particle that has the same dimension as in the Primary PSO is randomly generated, before entering the Secondary PSO. The departure time for each vehicle and edge type between each pair of nodes are decoded from the generated values.

st Solution representation
The solution representation for the Secondary PSO was based on the routes derived from a local improvement as presented in Figure  3. The total dimension is v max +(n+v used ) where v used is the number of vehicles that have actually been used in local improvement. The second part is generated randomly in the range of [0,1] with (n+v used ) dimension as initial solutions.

rd Decoding method
The decoding process is as follows: Step 1 Determine the departure time of each customer Refer to the first v max positions, position v represents the departure time of vehicle v. Ignore the value of position v if vehicle v is not used.
Step 2 Determine the type of edges in each route From the value of position v max +1 in the solution representation and the existing routes, the position is decoded as the type of the first edge in the first route (the first assigned vehicle during the decoding procedure in the Primary PSO and local improvement) and so on. In each position, if the value is less than 0.5, it is decoded as a designated edge; otherwise, it is an alternative edge. The decoding process of this part is shown in Figure 4.

PSO Algorithm for TDAVRP
The procedure of the two-stage PSO is displayed in Figure 5. The algorithm begins with the Primary PSO, after decoding the particles into different routes and improving the quality of the solutions by local improvement. It then precedes to the Secondary PSO and performed Secondary PSO operations. Finally, the best particles with the lowest objective value are used in the Secondary PSO to perform the Primary PSO operations. The stopping rules in these two PSOs are all the given number of iterations.

Computational Results
The two-stage PSO algorithm is coded by C ++ . These experiments are all carried out by a PC with Intel(R) Core(TM) 2 Duo CPU E7500 processor 2.93 GHz, 2.94 GHz, 4.00GB RAM.
To demonstrate and evaluate the proposed algorithm, the data and parameters we defined for these purposes are listed in Appendix B, of which Table 1 lists the speed distribution adopted from Donati et al. [5]. Table 2 lists the data for the parameters of the proposed model; and Table 3 sets the data used in the proposed algorithm.
Based on these data, analyses are conducted by first, comparing the results of the PSO with ILOG CPLEX to test the accuracy and efficiency of the two-stage PSO in Section 6.1. Then, for the purposes of analyzing the properties of the algorithm and the results, tests with different benchmark problems and sensitivity analysis are carried out in Sections 6.2 and 6.3. In particular, the comparison of the TDAVRP with the TDVRP is conducted with discussion in Section 6.4.

Computability analysis
By using the speed distribution functions defined in Table 1 of Appendix B, Table 4 shows the computational results of CPLEX and the PSO algorithm using Solomon's small scaled test problems R101 with 5 nodes, in which the values are the averages of the results of 10 test problems with different combinations of travel speed functions.
The error rates of the two-stage PSO with respect to optimal solutions obtained from CPLEX were less than 0.50%, and all of the computation times of the PSO are less than 10 seconds as within 0.10% of using CPLEX. There were 2 or 3 vehicles needed to serve customers,  Route v used The type of edge is decided by this dimenstion      which indicates that the results tends to assign more vehicles to avoid violations of time windows due to the penalties.
When the size of the problem increases to 25 nodes, CPLEX is unable to provide a feasible solution. Table 5 presents test results of instances with 25, 50, and 100-nodes, respectively. The results have shown that the algorithm is able to obtain feasible solutions within 13

Sensitivity analysis
In the objective function, there are three different costs of fixed cost, time cost and distance cost. In order to analyze the influences of the corresponding parameters on the transportation planning, three evaluation indexes have been considered to evaluate the performance of each solution. They are the usage rate of vehicles, the violation of time windows, and the percentage of times that alternative edges are selected, respectively. The usage rate of vehicles has been evaluated by the average number of customers served by one vehicle. And the violation of time windows is the sum of the waiting time and lateness time in one solution. The percentage of alternative edges is the number of selected alternative edges compared to the number of all used edges.
Parametric analysis: From the test results in 100 nodes, the usage rate of vehicle is about 100 1.43 70 ≈ . Obviously, this rate is too low to be realistic in this case and has to be improved. To determine a method to enhance this rate, possible factors of unit fixed cost, average speed and tolerances have been taken into account and evaluated.
The test results in the parametric analysis indicated that there is trade-off between number of vehicles and violation of time windows in the same instance, that is, the lower is the number of vehicles, the higher is the violation. Since the usage rate of vehicles is sensitive to the average speed and tolerances, more attention needs to be paid on these parameters in practice. The increase of the average speed would also decrease the percentage of alternative edges due to the reduction of the average travel time of the designated edges. The adjustment of the width of time windows has influences on the usage rate of vehicles and violation of time windows, especially the violations would decrease rapidly while the width of time windows changed from 10 to 20 (Unit time) as shown in Figure 6.

Comparison of Different Benchmark Instances
Because our interest in this study is on the vehicle routing when the alternative edges are considered, particular attention is placed in this section to the behavior of alternative edges while different test problems categorized by R, C and RC in Solomon's benchmark problems are adopted [8]. Table 6 presents the further test results of six problem sets that the values of each set are the averages from the results of eight instances (each instance is run 10 times). Error rate=(PSO Min_cost−(CPLEX optimal solution)) / (CPLEX optimal solution).
Ave_Vehicles=the average number of vehicles (routes) in 5 running times.   In the same category of problems, the percentages of alternative edges are close, but the percentage in the problem with a long scheduling horizon is higher than the problem with a short scheduling horizon. This is possibly because that the earlier is the service, the lower is the total cost. When the designated edges are in off-hour, the designated edges must be selected because the time and distance cost of the designated edges are lower than the alternative edges while the other factors are fixed. However, if the designated edges are congested, the decision would evaluate the weighted sum of the time and distance cost of these two edges, in which the designated edges have lower distance cost but higher time cost. In problems of Category R, the longer is the scheduling horizon, the larger the width is. When the scheduling horizon is long and the designated edges are congested, the possibility of choosing the alternative edges would increase because the reduction of the violations would reduce the effects of higher distance cost of the alternative edges. This leads to the higher percentages of alternative edges in long scheduling horizons.
In addition, we also observed that the percentages of choosing the alternative edges in problems of Category C were higher than other categories. Since the difference of the distance cost between two types of edges is close in the problems with clustered distribution, the alternative edge with lower violations will be chosen more frequently than the designated edge in the peak hour.

Comparison with TDVRP
When compare the results of the TDVRP and TDAVRP, the numbers of vehicles are nearly equivalent but the violation of time windows in the TDVRP is higher than the TDAVRP which caused the differences in the total cost. The violations of time windows of the TDVRP and TDAVRP in different instances are shown in Figure 7, of which the violation of the TDVRP in every instance is higher than that of the TDAVRP and the differences are significant in instance R101 which indicates the high complexity of this problem to serve customers within their time windows, so the alternative edges is able to reduce the violation in the TDAVRP.
In conclusion, the time-dependent property would influence the usage rate of vehicles and violations significantly compared with the VRPTW. If the usage rates in the TDVRP and TDAVRP are close, violation of time windows in the former problem is higher than the latter. This has shown the significance of the existence of alternative edges in reducing the violations.

Summary and Discussion
This study intends to solve a Time Dependent Vehicle Routing Problem (TDVRP) in a multi-graph network of which alternative edges are considered. By taking account of the fluctuation of travel time of the designated edges, the considered Time Dependent Alternative Vehicle Routing Problem (TDAVRP) has shown to provide additional choices of the alternative edges in travelling and thus reduced the transportation cost.
In order to facilitate the applications and analysis, we have proposed a Mixed Integer Programming (MIP) model with weighted cost of number of vehicles, time and distance in the objective function to support decision makers in evaluation of different combinations of costs when the service time window and the vehicle capacity are restricted.
Since the TDAVRP is an NP-hard problem, a PSO-based algorithm has been developed. The solution process of the proposed algorithm is divided into two stages of Primary PSO and Secondary PSO. Besides, the procedure of local improvement is included to improve the solutions obtained from the Primary PSO before the Secondary PSO is carried out.
To evaluate the performance of the proposed algorithm, Solomon's benchmark instances are adopted [8]. From the test results of the small problems in comparison with those from CPLEX, the proposed twostage PSO can reach the optimum solutions with higher computation efficiency. For large scaled problems of which CPLEX cannot obtain feasible solutions, the proposed PSO has provided promising result. For instance, it can solve 100-node problems within 13 minutes.
In addition, from the relations between the usage rate of vehicles and the violation of time windows, we have indicated the factors which influence the percentage usage of alternative edges and have confirmed by our sensitivity analysis as summarized in the following.
The results of sensitivity analysis indicated that the usage rate of vehicles and violation of time windows in the TDAVRP are sensitive to the average speed, tolerances and width of time windows. The properties of the problems, like the distribution of customers and the length of scheduling horizon, also affect the usage rate and violations significantly. On the other hand, the percentage of alternative edges  would decrease when the congestion situation of the designated edges is improved by, for instance, the increase of average speed. The test results have also shown that the existence of the alternative edges is more significant in the problems with clustered customers. Finally, TDAVRP provides lower violations in time windows than that in TDVRP, which also shows that the former has better performance than the latter.

Conclusion
In conclusion, the contributions of this study are concluded as follows: First, the model helps decision makers reduce the violation cost by the consideration of the alternative edges. Second, the proposed heuristic method could obtain feasible solutions in a short time. Finally, sensitivity analysis and comparisons of the results have provided useful suggestions for transportation planning.
There are still some issues that need to be discussed further: for instance, validation of the travel times functions; the operations of local improvement to enhance the connection between the Primary PSO and Secondary PSO; and finally, real-world applications are expected to explore the influences of other parameters on the TDAVRP.