A Combination of Genetic Algorithm and Particle Swarm Optimization for Vehicle Routing Problem with Time Windows

A combination of genetic algorithm and particle swarm optimization (PSO) for vehicle routing problems with time windows (VRPTW) is proposed in this paper. The improvements of the proposed algorithm include: using the particle real number encoding method to decode the route to alleviate the computation burden, applying a linear decreasing function based on the number of the iterations to provide balance between global and local exploration abilities, and integrating with the crossover operator of genetic algorithm to avoid the premature convergence and the local minimum. The experimental results show that the proposed algorithm is not only more efficient and competitive with other published results but can also obtain more optimal solutions for solving the VRPTW issue. One new well-known solution for this benchmark problem is also outlined in the following.


Introduction
The vehicle routing problem (VRP) is a combinatorial optimization and integer programming problem seeking to service a number of customers with a fleet of vehicles. Proposed by Dantzig and Ramser in 1959, VRP is important to the fields of transportation, scheduling, distribution and logistics [1,2]. The problem involves many real-world considerations, such as time-window OPEN ACCESS constraints, driver waiting costs, backhauls, layovers, etc. The vehicle routing problem with time windows (VRPTW) has been extensively studied by many researchers from the fields of operational research, applied mathematics, network analysis, graph theory, computer applications, traffic transportation, etc. Firstly, VRPTW is still one of the most difficult problems in combinatorial optimization, consequently presenting a great challenge. Secondly, in a more practical aspect, study of this problem could provide a real opportunity to reduce the costs in the area of logistics [3].
The VRPTW is a generalization of the VRP where the service for a customer starts within a given time interval, and it has been the subject of intensive research efforts for both heuristic and exact algorithms. The actual solutions of VRP can be generally classified into two main categories: the exact algorithms and the heuristic algorithms. The main approaches for solving VRPs are shown in Table 1. Table 1. Main approaches for solving VRPs.

Algorithms Remarks
The exact algorithms Branch and bound method [4,5] The Efficiency depends on the depth of the branch and bound tree.
Set segmentation method [6,7] Hard to determine the minimum cost for each solutions.
Dynamic programming method [8,9] Effective to limited-size problems, hard to consider the concrete demands such as time windows.
The heuristic algorithms The traditional heuristic algorithms Savings algorithm [12,13] Computes rapidly, hard to get the optimal solution.
Sweep algorithm [14,15] Suitable to the same number of customers for each route with few routes.
Two-phase algorithm [16,17] Hard to get the optimal solution.
The meta-heuristic algorithms Tabu search algorithm [18][19][20] Has the good ability of local search, but is time consuming, and depends on the initial solution.
Genetic algorithm [13,21] Has the good ability of global search, computes rapidly, hard to obtain the global optimal solution.
Iterated local search [22,23] Has the strength of fast convergence rate and low computational complexity.

Variable neighborhood
Search [26,27] Is suitable for large and complex optimization problems with constraints.
Ant colony algorithm [28][29][30] Has good positive feedback mechanism, but is time consuming and prone to stagnation.
Neural network algorithm [31,32] Computes rapidly, has slow convergence and can easily be trapped in a local optimum Artificial bee colony algorithm [30,33] Achieves a fast convergence speed, is associated with the piecewise linear cost approximation.
The exact algorithms can obtain the exact solution, but the computational effort required to solve this problem increases exponentially with the problem size. The traditional heuristic algorithms often only get the approximate solution close to the optimal solution, and are limited to the smaller problems. When the size of the problems increases, the solution precision of the traditional heuristic algorithms is often poor. The traditional heuristic algorithms adapt to local optimization and combined with the meta-heuristic algorithms to improve the solutions [39]. For large, complex problems, only the meta-heuristic algorithms can be used due to their strong performance of global search [18,40,41].
The VRPTW is a non-deterministic polynomial-time hard (NP-hard) problem. Due to the complexity of the VRPTW, it is not easy to obtain an exact solution for a large problem in real time. For such problems, optimal solutions are found quickly and are sufficiently accurate. Usually this task is accomplished by using various meta-heuristic algorithms, which rely on some insights into the nature of the problem. Particle swarm optimization (PSO) is not only superior in terms of high accuracy speed calculation, as well as its simple program, but it is also robust. Nie and Yue integrated the concept of evolving individuals originally modeled by GA with the concept of self-improvement of PSO [42]. Hao et al. proposed a modified particle swarm optimization which took the crossover between each particle's individual best position [43]. In [44], dynamic parameterized mutation and crossover operators were combined with a PSO implementation individually and in combination to test the effectiveness of these additions. In [45], the proposed method used the concept of particles' best solution and social best solution in the PSO algorithm, followed by combining it with crossover and mutation of GA. Considering the particle real encoding, linear decreasing inertia weight function and crossover operator of genetic algorithms, a combination of genetic algorithm and PSO for VRPTW is proposed that can improve the performance when computing speed to obtain the optimal solutions.

Vehicle Routing Problem with Time Windows (VRPTW)
In the VRPTW, a fleet of K identical vehicles supplies goods to n customers. All the vehicles have the same capacity Q . For each customeri , 1, 2, , Each vehicle serves a subset of customers on the route. The vehicles have the same capacity Q . The total sum of demands of customers served by a vehicle on a route cannot exceed Q . The additional constraints are that the service beginning time at each node ( 0,1, , ) i c i n =  must be greater than or equal to i b , the beginning of the time window, and the arrival time at each node i c must be lower than or equal to i e , i.e., the end of the time window. In case the arrival time is less than i b , the vehicle has to wait until the beginning of the time window before starting servicing the customer. The goal is to find a set of routes which can guarantee each customer to be served by one vehicle within a given time interval and then satisfy the vehicle capacity constraints. Furthermore, the size of the set should be less than the number of vehicles needed and the total travel distance should be minimized. Moreover, the mathematical formulation of the VRPTW is presented as follows [4,46]: Subject to: (1 ) ( , [1, ]; where ijk x is 1 if vehicle k travels from customer i to customer j , and 0 otherwise. i t denotes the time a vehicle starts services at customer i ; M is an arbitrary large constant. Objective function states that the total cost is to be minimized. Constraint (2) specifies that there are no more than K routes going out of the depot. Constraint (3) ensures that one vehicle goes into and out of a customer exactly. Equation (4) is the capacity constraint. The time window is assured in Equation (5), Equation (6) and Equation (7). Equation (8) is the flow conservation constraints that describe the vehicle path.

Particle Swarm Optimization (PSO)
PSO is a population-based stochastic optimization technique developed by Eberhart and Kennedy in 1995, inspired by social behavior of bird flocking or fish schooling [34] and was first intended for simulating these organisms' social behavior. It is best to imagine a swarm of birds that are searching for food. When one of them finds the food, some of them will follow the first bird, while others will find other food. Initially, the birds do not know where the food is, but they know at each time how far the food is. Each bird will follow the one that is nearest to the food. Throughout the course of preying, a bird will use its own experiences and collective information to search for food.
In particle swarm optimization, the particles are moved around in the search space according to a few simple formulae. The position of a particle represents a candidate solution to the optimization problem at hand. Each particle searches for better positions in the search space by changing its velocity according to rules originally inspired by behavioral models of flocking. The movements of the particles are guided by their own best-known position in the search space as well as the entire swarm's best-known position. When improved positions are discovered, these will then come to guide the movements of the swarm. The process is repeated and by doing so it is hoped, but not guaranteed, that a satisfactory solution will eventually be discovered. Compared with other intelligence optimization algorithms such as ant colony optimization, genetic algorithm, simulated annealing algorithm, neural network algorithm and artificial immune algorithm, PSO retains the global search strategy based on the swarm and has no individuals as in crossover and mutation. In PSO, through adjusting the velocities and positions of the particles which fly through the problem space by following the current optimum particles, the optimal solution can be obtained. Due to its simplicity, strong robustness, and fast optimization speed, PSO is suitable for very large and complex optimization problems with constraints. At first, PSO was applied to solve continuous optimization problems; however, several applications were proposed during these years in the area of combinatorial optimization problems including shop scheduling [47,48], project scheduling [49], travelling sales force [50], partitional clustering [51], and vehicle routing [34,52].
PSO is one of the evolution algorithms with the characteristics of evolutionary computing and swarm intelligence. Similar to other evolutionary algorithms, PSO searches for optima by evaluating individual fitness based on cooperation and competition between individuals. In PSO, each individual is considered as a particle without weight and volume in n -dimensional search space and flies through the space with a certain speed. The speed is adjusted dynamically by the individual's experience and the entire swarm's experience.
PSO is initialized with a population of random solutions and searches for the optimal solution by updating the particle's position. Each particle is the feasible solution and is designated a fitness value by the objective function. Each particle keeps the track of its coordinates in the problem space which are associated with the best solution (fitness) it has achieved so far. The fitness value is called best p .
When a particle takes all the population as its topological neighbors, the best position is a global best and is called best g . Through best p and best g , particles update themselves to produce the next generation of swarms. The selection of fitness function depends on the research goals. The fitness function to evaluate the individuals is always related to the objective function. For a VRPTW, the total cost can be viewed as the fitness value. The inverse of the total cost is used to represent the fitness of the individuals, and then the fitness function is defined as follows: In a PSO algorithm the particles represent potential solutions to the problem, and the swarm consists of P particles. Each particle p can be represented through n -dimensional vectors: the first one is defined by 1 2 ( , , , ) that indicates the position of the particle p in the searching space at the iteration t . The second one is 1 2 ( , , , ) that represents the velocity with which the particle p moves. The third one is the best position of the pth particle and the last one is  that represents the global best position in the swarm until tth iteration. The swarm is updated by the following equations: where 1 k and 2 k are acceleration coefficients, which are respectively called cognitive and social parameter; 1() rand and 2() rand are two random numbers uniformly distributed in [0,1] . Acceleration coefficients 1 k and 2 k are positive constants to control how far a particle will move in a single iteration. Low values allow particles to roam far from target regions before being tugged back, while high values result in abrupt movement towards, or past, target regions. Typically, these are both set to a value of 2.0, although assigning different values to 1 k and 2 k sometimes leads to an improved performance.
A constant, max v , is used to arbitrarily limit the velocities of the particles t pn v and improve the resolution of the search. When max v is large ( 5) ≥ , the velocity of the particle is large, too; it is conducive to a global search, though it may fly through the optimal solution [53][54][55]. When max v is , the velocity of the particle is also small; it leads to a fine search in a specific region, but it is easy to fall to local optimum [53][54][55]. In a word, the search efficiency depends on max v .
Each particle moves in the search space with a velocity according to its own previous best solution and its group's previous best solution. In Equation (10), the velocity 1 t pn v + consists of three parts-momentum, cognitive and social parts-respectively each term of the right side of Equation (10). The momentum part denotes the previous velocity of the particle, which improves the ability of the global search. The cognitive part denotes the process of learning from an individual's experience. The social part denotes the process of learning from others' experience, which represents the information sharing and social cooperation between particles. The balance among these parts determines the performance of a PSO algorithm. Without the momentum part, the particle then moves in each step without knowledge of the past velocity. Without the cognitive part, the convergence speed is fast, but can easily fall into a local optimum for a large problem size. Without the social part, it is hard to get the optimal solution due to the lack of the communications among individuals.

Particle Encoding
Encoding is a bridge connecting a problem with an algorithm. Encoding method and initial solution have a great impact on the VRP problem. It is the key step to finding the appropriate encoding method for the particles and the corresponding solutions. In brief, the encoding methods of PSO include real encoding, binary-encoding and integer encoding. Integer encoding is easy to decode and convenient for the fitness function calculation, but it requires a lot of computing resources and tends toward premature convergence. In addition, binary encoding is hard to decode. Therefore, real encoding is adopted in the proposed method.

Inertia Weight Function
In order to better control the searching and exploring abilities and improve the convergence of particle swarm algorithm, the inertia weight ω is introduced into the velocity function. Therefore, 1 t pn v + is changed into the following form: The inertia weight ω is employed to control the impact of the previous history of velocities on the current one. Accordingly, the parameter ω regulates the trade-off between the global and the local exploration abilities. If ω is high, particles can hardly change their direction and turn around, which consequently implies a larger area of exploration as well as a reluctance against convergence towards optimum. On the other hand, If ω is small, only little momentum is presented from the previous time-step, thereby leading to quick changes of direction. A suitable value for the inertia weight usually contributes the balance between global and local exploration abilities and consequently results in a reduction of the number of iterations required to obtain the optimum solution.
Considering whether the inertia weight ω is changed or not, the calculation methods of the inertia weight ω include three categories: fixed-weight method, time-varying weight method and adaptive inertia weight method. The fixed-weight method selects a constant value as the inertia weight ω and is kept unchanged. The fixed-weight method was originally introduced by Shi and Eberhart in [56]. The time-weight method selects an iterative relationship as the inertia weight ω according to a selected range, which is defined as a function of time or iteration number. In [57][58][59][60], time-varying inertia weight strategies were introduced and shown to be effective in improving the fine-tuning characteristic of the PSO. In [61][62][63], the adaptive inertia weight methods used a feedback parameter to monitor the state of the algorithm and adjusted the value of the inertia weight.
In order to better balance the global search ability and local search capabilities, kinds of inertia weight are often used. At first, a higher ω is selected to expand the search space and converge to a region. Then a smaller ω is selected to explore the local region to obtain high accuracy solutions. In this paper, the definition of inertia weight is a linear decreasing function of the number of iterations. ω is given as follows: where 0 ω is the initial value of inertia weight, e ω is the final value of inertia weight, max η is the maximum number of iterations, η is the current number of iterations. From Equation (12), by the linear function decreasing inertia weight, it is easy to obtain better global search ability and make particles enter the area around the optimal value early in the iteration process, and it is easy to obtain better local search ability and the solutions close to the optimum value late in the iteration process.
Since the value of the inertia weight is mainly determined based on the iteration number, this strategy is relatively simple and has fast convergence compared to other methods.

Crossover Operator of Genetic Algorithm
In order to expand the search space to obtain the optimum solution and enhance the diversity of particles, the idea of genetic algorithm is integrated into the PSO to avoid the premature convergence and local minimum value. Crossover operator of genetic algorithm is introduced. The crossover operators of the particle's velocity and position are given as follows: where child p represents the offspring of the particle, parent p represents the parent of the particle, c p represents the crossover probability. c p ( [0,1] c p ∈ ) is a random number. Based on the trial and error, 0.2 was found to be a suitable value for c p . The two parents are combined to produce two new offspring. If the fitness of offspring is more than the fitness of parents, the offspring will be selected to replace the parents; otherwise, the offspring are discarded.

The Proposed Algorithm Flow
The flow of the proposed algorithm is shown in Figure 1. In this flow, the parameters are set in step (1), and the particles are initialized in step (2). Its position is decoded in step (3), its corresponding fitness value is evaluated in step (4), its cognitive and social information is updated in step (5), and then moved by step (6). Step (7) is the controlling step for repeating or stopping the iteration.
Step (8) generates a new set of the population. Note that the improvements of this flow from the original PSO take place in step (3), which uses the real encoding method to decode the route in addition to step (6), which introduces a linear decreasing function of the number of iterations, and step (8), which is combined with the crossover operator of genetic algorithm. The proposed algorithm steps are given as follows:  (3) Particle encoding. According to the particle encoding rules, for (6) Particles updating. Update the velocity and the position of each pth particle according to Equation (10). (7) Termination judgment. If the stopping criterion is met, go to step (9). Otherwise, 1 η η = + and go to step (8). The stopping criterion is that max η η ≥ or finding a better solution. A better solution means that the hierarchical cost objective value is better than that of the best solution found so far.
(8) Crossover operator. Generate random number c p . According to Equations (13)- (16), generate a new set of population. Then return to step (3). (9) Outputting the optimal solution. Decode best g as the best set of vehicle route R and output the optimal solution R .

Experimental Results
In order to evaluate the performance of the proposed algorithm, we implemented the algorithm in Visual C# under Microsoft Windows-XP on a PC with Intel P4 4 GHz CPU and 2 GB RAM. The VRPTW benchmark problems of Solomon, which have been the most commonly chosen to evaluate and compare all algorithms, are tested in this paper.

Solomon Benchmark Problems
It is well known that Solomon Benchmark problems are the widely used standard test set for VRPTW. Solomon Benchmark derives from the website [64]. The Solomon test set consists of 56 problem instances for each dimension category problem, i.e., 25, 50 and 100 customers. Each of these instances comprises 100 customers. The location of the depot and the customers are given as integer values from the range 0,1, ,100  in a Cartesian coordinate system. The distance between two customers is the simple Euclidean distance. It is assumed that the travel times are equal to the corresponding Euclidean distances between the customer locations. One unit of time is necessary to run one unit of distance by any vehicle. Different capacity constraints are considered for the vehicle in each class of instance, as well as the demands from the customers. The test problems are grouped into six problem types: R1, R2, RC1, RC2, C1, and C2, each containing 8 to 12 instances. In R1 and R2, the customer locations are generated randomly in a given area according to a uniform distribution. C1 and C2 have customers located in clusters. RC1 and RC2 contain a mix of randomly distributed and clustered customers. R1, C1 and RC1 have narrow time windows and the vehicles have only small capacities. Therefore, each vehicle serves only a few customers. In contrast, R2, C2 and RC2-have wider windows and the vehicles have higher capacities. Each vehicle supplies more customers and therefore, compared to the type 1 problems, fewer vehicles are needed. Because the Solomon's test set represents relatively well different kinds of scenarios, it has often been chosen to evaluate many solution proposals in the literature.

Evaluation of the Results
In order to compare the solutions, four evaluation indexes are presented: the total traveled distance, the average total traveled distance, the CPU runtime and the quality of the solutions. The total traveled distance and the average total traveled distance are the main criteria in the VRPTW and determines the algorithm merits. The CPU runtime describes the algorithm efficiency. The quality of the solutions is the comprehensive index, which denotes the percentage of deviation from the best-known solution. The formula for calculating the percentage of deviation is as follows: where dev z is the percentage of deviation from the best-known solution, z is the cost of the current solution, best z is the cost of the best-known solution.

Analysis of the Results
Each of 56 Solomon Benchmark problems is performed by some of the best-reported methods for VRPTW, namely, the genetic algorithm, ant colony optimization [29], PSO [65] and the proposed algorithm. Each problem involves 100 customers, randomly distributed over a geographical area. For each problem, 10 replications of the algorithm are attempted. According to Equation (1), we have adopted a formulation in which the total traveled distance is the optimized objective. Results displayed in the following tables contain the total traveled distance and the number of vehicles used in order to facilitate the comparison of our results with those obtained by hierarchical approaches in which the total traveled distance is the primary objective and, for the same total traveled distance, the secondary objective is number of vehicles, and also with multi-objective formulations that simultaneously consider both objectives [66]. It compares the total traveled distance (TD), the average total traveled distance, the CPU runtime and the quality of the solutions for each of the problems. Both best and average results are presented.

The Total Traveled Distance
The results of TD and the number of vehicles (NV) are shown in Table 2. The published best-known solutions are not obtained by one or a particular class of methods [15,[67][68][69][70][71]. Results emphasized in bold in Table 2 represent the new best solutions reached by the proposed algorithm: 1 out of 56 (1.79%). Results emphasized in bold and italic in Table 2 represent the previously best-known solutions that cannot be reached by this algorithm. The results show that many previously best-known solutions from the proposed algorithm have been reached: 31 out of 56 (55.4%). Globally, the TD average results are close to the optimal solutions known in the proposed algorithm. This comparison shows that the results from the proposed algorithm are competitive with other published results. In the genetic algorithm, 27 out of 56 (48.2%) have reached the previously best-known solutions. In PSO, 28 out of 56 (50%) have reached the previously best-known solutions. In ACO, 26 out of 56 (46.4%) have reached the previously best-known solutions. It can be seen that the proposed algorithm is better than the genetic algorithms, ACO and PSO. The efficiency is improved and the results are close to the best-known solutions. It is possible to see that the proposed algorithm continues to be very competitive in terms of total TD.  TD  NV  TD  NV  TD  NV  TD  NV  TD  NV  TD  NV  TD  NV  TD  NV  TD TD  NV  TD  NV  TD  NV  TD  NV  TD  NV  TD  NV  TD  NV  TD  NV  TD  NV  27 R110 1080. 36 11 1080. 36

The Average Total Traveled Distance
The results of the average total traveled distance are shown in Table 3. Each entry refers to the best performance obtained with a specific technique over a particular data set. Rows C1, C2, R1, R2, RC1 and RC2 present the average number of vehicles and average total distance with respect to the six problem groups. The last row refers to the cumulative number of routes and traveled distance over all problem instances. The first column describes the various data sets and corresponding measures of performance defined by the average number of routes (or vehicles), and the total traveled distance. The following columns refer to particular problem-solving methods. The performance of the proposed algorithm is depicted in the last column. The results indicate that the proposed algorithm can match the best-known solutions, and it performs better than other algorithms. In addition, the proposed algorithm is the only method that found the minimum number of tours and the comparable traveled distance consistently for all problem data sets. The results of the average CPU runtime are listed in Table 4. Rows C1, C2, R1, R2, RC1 and RC2-present the average CPU runtime with respect to the six problem groups. From Table 4, for the same instance, ACO is almost the most time consuming; the genetic algorithm takes second place; the PSO takes third place; and the proposed algorithm is last. The result shows that the proposed algorithm is highly efficient and suitable to be used in real-life problems to generate good solutions for further improvement. The average CPU runtime required to solve each type of problem set is between 62 and 149 milliseconds by the proposed algorithm, which is much faster than other algorithms. It can be seen that the proposed algorithm is very effective. This happened because crossover in GA is done between random chromosomes, whereas in the proposed algorithm, crossover is done between a particle's chromosome and best chromosome. Therefore, the proposed algorithm does not only give a better result, but it also reaches convergence results faster than other methods.  (17), the quality of the solutions is shown in Table 5. Here z is the average total traveled distance. Rows C1, C2, R1, R2, RC1 and RC2 present the average dev z with respect to the six problem groups. The last row refers to the cumulative dev z over all problem instances. It can be seen that, with the proposed method, the quality of the results is between 1.32% and 8.17% with average quality equal to 4.07%; in the ACO, the quality of the solutions is between 5.07% and 9.75% with average quality equal to 6.95%; in PSO, the quality of the solutions is between 2.32% and 8.53% with average quality equal to 5.63%; in the genetic algorithm, the quality of the solutions is between 3.39% and 8.96% with average quality equal to 6.29%. The improvement in the quality of the solutions was achieved with the addition of the crossover operator of genetic algorithm. The reason is that, now, the particles moved in a more fast and efficient way to their local optimum or to the global optimum solution (to the best particle in the swarm). It is proved that the addition of the evolution of the population phase before the individuals used in the next generation improves the results of the algorithm, especially in the large-scale vehicle routing instances which are more difficult and time consuming. From the results shown in Tables 2-5, it is observed that the proposed algorithm is fast and produces good solutions, which are only slightly less accurate than the best-known results. This comparison also shows that the results from the proposed algorithm are competitive with other published results. The experimental results reveal that the proposed algorithm is highly capable at minimizing the total travel distance. In particular, it is obvious that the proposed algorithm has great advantages for the solved problem size over the genetic algorithm, PSO and ACO. The reason is that, now, the particles moved in a more fast and efficient way to their local optimum or to the global optimum solution (to the best particle in the swarm). After each environment change (peak movement), the particles are placed on a point far from the optimum through the inertia weight function. The particles take few iterations to reach this point. At the same time, the particle's coding can be ensured that a particle is always able to generate a new feasible solution.

Conclusions
This paper proposes a genetic and PSO algorithm for VRPTW. Moreover, this algorithm was applied to the Solomon Benchmark problems and produced satisfactory results. Compared to other approaches, experimental results indicate that the proposed algorithm is fast and possesses a relatively accurate results. The proposed algorithm has proved to be an effective and competitive algorithm for the optimization problems due to its easy implementation, inexpensive computation and low memory requirements. Three major contributions are as follows: (1) The real encoding method avoids the complex encoding and decoding computation burden.
(2) A linear decreasing function of the number of iterations in PSO have a flexible and well-balanced mechanism to enhance and adapt to the global and local exploration abilities, which can help find the optimal solution with the least number of iterations. (3) The crossover operator of the genetic algorithm is introduced to generate a new population guaranteeing that the offspring inherits good qualities from this parent. The crossover operator avoids premature convergence and local minimum value and increases the diversity of particles.
Although the combination of the genetic algorithm and PSO for VRPTW obtains satisfactory achievements, there is still some room for improvement, such as how to effectively construct the initial solution and how to precisely judge the local optimal solution. Further research applying the proposed algorithm to other VRPTW variants will be carried out in our future work.