Article A Hybrid Chaos-Particle Swarm Optimization Algorithm for the Vehicle Routing Problem with Time Window

State-of-the-art heuristic algorithms to solve the vehicle routing problem with time windows (VRPTW) usually present slow speeds during the early iterations and easily fall into local optimal solutions. Focusing on solving the above problems, this paper analyzes the particle encoding and decoding strategy of the particle swarm optimization algorithm, the construction of the vehicle route and the judgment of the local optimal solution. Based on these, a hybrid chaos-particle swarm optimization algorithm (HPSO) is proposed to solve VRPTW. The chaos algorithm is employed to re-initialize the particle swarm. An efficient insertion heuristic algorithm is also proposed to build the valid vehicle route in the particle decoding process. A particle swarm premature convergence judgment mechanism is formulated and combined with the chaos algorithm and Gaussian mutation into HPSO when the particle swarm falls into the local convergence. Extensive experiments are carried out to test the parameter settings in the insertion heuristic algorithm and to evaluate that they are corresponding to the data’s real-distribution in the concrete problem. It is also revealed that the HPSO achieves a better performance than the other state-of-the-art algorithms on solving VRPTW.


Introduction
Logistics management plays an important role in reducing costs and improving the competitiveness of enterprise.The logistics companies have to make strategic and operational decisions to schedule the business processes more efficiently.Many companies use the method of optimizing vehicle routes so as to reduce costs and to improve logistic service quality effectively.Among them, the vehicle routing problem (VRP) draws great interest, which presents promising performance in chain management and many other areas.The VRP can be defined as determining an optimal set of routes for vehicles, which would save time and cost, so as to effectively strengthen the logistic companies' competitiveness.One variants of VRP is the vehicle routing problem with time windows (VRPTW).VRPTW imposes the appropriate time window and time constraints to each task node (including dispatch node and receiving node) on the basic VRP.In this way, VRPTW is able to improve service quality.VRPTW has a wide range of applications, such as distribution requirements planning [1,2], material handling systems [3], bank deliveries [4], postal deliveries [5], school bus routing [6] and so on.Actually, VRPTW has been a research field of great interest, and many different solving algorithms have been proposed.
Despite its great potential in solving the vehicle routing problem, VRPTW is still a difficulty, since it belongs to NP-hard problems.Many researchers put forward a variety of accurate methods to solve the problem, such as Multi-Objective Programming, Linear Programming and Integer Programming.These polynomial methods with constraints can perform effectively and get the optimal solution in the small-scale VRP.However, with the increase of the problem scale, it needs much more time to match the explosive growth.Therefore, much effort has been done on the heuristic intelligent algorithm to solve VRPTW.The heuristic algorithm is simple to implement and does not present an obvious time increase with the problem scale increasing.The main problem is that the iteration speed is slow in the early stage, and it is easy to fall into local optimal solutions.In a word, both the accurate solution and heuristic solution have weakness on solving VRPTW.
Nowadays, the problem scale of real-world VRPTW increases rapidly.The accurate algorithms need too much time to model and search for a proper solution, while single heuristic algorithms often fall into local optimal solutions and need a large number of iteration times.This paper gives a new insight into VRPTW and proposes an effective method by combining different VRPTW solutions.Particle swarm optimization (PSO) is a promising method, since its iteration speed is fast.However, it often falls into local optimal solution.The chaos optimization algorithm (COA) and the Gaussian strategy can upgrade solution diversity and avoid falling into the local optimal.The insertion heuristic algorithm can construct the initial solution quickly and upgrade iteration speed.So, this paper combines these algorithms and proposes a hybrid chaos-particle swarm optimization algorithm (HPSO) for VRPTW.The contributions of this paper are as follows: (1) A PSO based on COA is proposed, so as to avoid the problem that the number of the feasible solution particles is too small in the early state and the speed of the convergence is low.When the feasible solution particle first appears, COA is able to re-initialize the particle swarm, which can increase the number of the feasible solution particles and accelerate the convergence.(2) The insertion heuristic algorithm is introduced into our method to construct the route of VRPTW in the decoding process, which is much easier and more effective than conventional methods.
(3) COA and Gaussian strategy are used to judge the premature convergence and to improve the particle swarm when the solving process falls into local optimal solutions.
The remainder of this paper is organized as follows.Section 2 introduces the related work; Section 3 describes the VRPTW model; Section 4 details the proposed HPSO; Section 5 presents various experiments to compare HPSO with other state-of-the-art methods on solving VRPTW.Finally is a summary of the paper and future outlooks.

Related Work
Many researchers used different heuristic algorithms for VRPTW.Instead of obtaining the optimal solution of VRPTW, like the accurate algorithms, heuristic algorithms can quickly get the approximate optimal solution.When the problem scale of VRPTW is large, the advantage of the heuristic algorithm is more obvious.Thus, early studies are based on single heuristics and metaheuristic, such as local search [7], minimum K-trees [8], adaptive memory [9], simulated annealing (SA) [10], tabu search (TS) [11,12], genetic algorithm (GA) [13] and ant colony optimization (ACO) [14].When the problem scale of VRPTW is small, these methods mentioned above can get an optimal result with limited time.However, when the problem scale of VRPTW increases, they will easily fall into the local optimal, and the iteration time cost will increase explosively.So, nowadays, many researchers focus on combing two or three heuristics and metaheuristic methods to solve VRPTW.
(1) A typical hybrid method is TS combined with other algorithms [15][16][17].Tabu search keeps a list of forbidden transformations to transfer the solution step that deteriorates the objective function value; while a hybrid algorithm based on TS and variable neighborhood descent uses a sweep method to obtain an initial solution [18].(2) Local search combined with other algorithms was also proposed to solve VRPTW.
Bianchessi et al. citeR19 proposed a local search heuristic method based on TS using a variable neighborhood structure, in which the node-exchange-based and arc-exchange-based movements are combined.Gendreau et al. [20] proposed a hybrid neighborhood search heuristic to optimize the planned route of vehicle in a context, where new requests with a pick-up and a delivery location would occur in real-time.Zachariadis et al. [21] proposed a hybrid local search approach based on TS, which can efficiently examine the rich solution neighborhood by statically encoding tentative moves into the special data of a diversification component, which is based on the aspiration criteria of TS.Hong [22] proposed a hybrid and improved the large neighborhood search (I-LNS) algorithm to solve the VRPTW.Utilizing the remove-reinsert process of the LNS, the latest request nodes are regarded as a part of the removed nodes; these nodes can be inserted into the original or current solution in good time in the reinsertion process; meanwhile, its computing speed is high and effective, for it does not need to resolve the primal problem each time.
(3) SA combined with other algorithms was used to solve VRPTW.Tavakkoli-Moghaddam et al. [23] proposed a liner-integer model based on SA, which constructs VRPTW with the independent route length in order to minimize the heterogeneous fleet cost and the capacity utilization.TTavakkoli-Moghaddam et al. [24] also proposed a new mathematical model based on SA to find the short routes with the minimum fleet travel cost, the minimum travel distances and the maximum sale.
(4) ACO combined with other algorithms provided another way to solve VRPTW [25][26][27].The weakness of ACO for VRPTW is that ants tend to create infeasible solutions with non-routed clients.Most implementations of ACO combined with other algorithms return feasibility to the solutions by applying a simple ad hoc post-insertion procedure, and when the post-insertion fails to route all the clients, the partial solution is discarded.( 5) GA can also be combined with other algorithms to solve VRPTW [28,29].GA mimics the evolution process in solving problems.The basic operation is the mating of two tours in order to form a new tour.Moreover, GA combined with other algorithms uses algorithmic anagoges to mutation and selection.The weakness of GA for VRPTW is the encoding of the chromosome.
With the problem scale increase, the difficultly of encoding the chromosome increases.
From the above brief review of different hybrid methods, it is concluded that many attempts have been made to solve VRPTW by combining different approaches.Two main properties of all the above methods in the literature can be summarized as follows: (1) Firstly, use a heuristic method to obtain a feasible initial solution; (2) Then, hybridize the metaheuristic method with another heuristic method to improve the solution during the iteration process.Particle swarm optimization (PSO) is a recent and promising evolution metaheuristic and has been successfully applied to solve VRPTW.Many researchers conduct research on solving VRPTW by combining PSO with other heuristic methods.Chen et al. [30] proposed a discrete PSO algorithm combined with SA heuristics.He got the high precision vehicle routes for some vehicle routing problems with less than 134 customers.However, the algorithm would cost a lot of computing time.In the worst case, it would cost more than half an hour to get the solution.Kachitvichyanukul [31] proposed two new PSO algorithms (SR-1 and SR-2) for encoding and decoding particles into vehicle routes.By comparing SR-1 with SR-2, he got the conclusion that SR-2 is better than SR-1.Zhao et al. [32] proposed a new encoding method based on PSO.They conducted a 2L-dimensional position vector of particles.Among them, L represented the number of service customers.Then, 2L-dimensional space corresponding to each particle was divided into two L-dimensional vectors: one vector meant the vehicle number for each customer, and the other vector represented the service order for each customer in the vehicle route.As to the particle decoding process, Zhang et al. [33] proposed a simulated annealing algorithm to construct the vehicle route.Firstly, they used PSO to solve the sub-problem of task assignment.Then, they used the simulated annealing algorithm to solve the sub-problem of route optimization.In order to get the infeasible solutions for particle swarm, Wu et al. [34] put forward a local neighborhood mechanism to increase the diversity of the particles and avoid falling into the local optimal state.Gong et al. [35] proposed a set-based PSO to solve the discrete combinatorial optimization VRPTW and achieved better effectiveness and efficiency.All the above PSO based methods can get a satisfied solution, but when the problem scale increases, it will cause problems, such as the speed of the iteration becoming slow in the early stage and it is easy to fall into local optimal solutions.
Aiming at alleviating the above problems, this paper puts forward a new hybrid algorithm, HPSO.It uses the advantage of the COA and Gaussian strategy to upgrade basic PSO solution diversity to avoid falling into the local optimal.COA is also beneficial for accelerating the iteration in the early stage.This paper also proposes an insertion heuristic algorithm to construct the initial solution and upgrade the iteration speed, as basic PSO is difficult for constructing the initial solution.

VRPTW Model Formulation
VRPTW as an extension of the basic vehicle routing problem is generally described as: one yard has K cars.The deadweight of each car is q k (k = 1, 2, . . ., K).There are L customers who need service.Customer i demands the g i (i = 1, 2, . . ., L) goods (M axg i M axq k ).It assumes that the time of service for customer i is T i , and the start time of service is a certain range of time, [ET i , LT i ](ET i means the earliest start time and LT i means the latest start time ).If the vehicle arrives earlier than ET i , the vehicle need to wait until ET i .If the vehicle arrives later than LT j , customer i needs to be delayed.At last, it finds the vehicle route that costs least.The mathematical model is described as follows: The yard number is defined as 0, the customers are defined as 1, 2, . . ., L, and and so point, i(i = 0, 1, . . ., L), can represent the yard and the customers.
The symbols in the model are defined as follows: c ij is the transportation costs of i to j, which refer to the costs of distance and time.enumerates j is the time of vehicle arrival for customer j. c 1 is the cost per unit of time for the vehicle to wait for the customer.c 2 is the fine per unit of time for a later vehicle arrival.c 1 (ET j − S j ) is the cost of waiting when the vehicle arrives earlier than ET j .c 2 (s j − LT j ) is the fine of vehicle arrival later than LT j .K is the number of all vehicles.g i is the demand quantity of customer i. q k is the capacity of vehicle k. o i is service time.t ij is the travel time between customer i and customer j. w i is the waiting time of a vehicle at customer i.
Based on the symbols, we can obtain the VRPTW mathematical model as Equation (1): The constraints are defined as follows: (1) On each route, the sum of each customer's demand cannot exceed the capacity of the vehicle, which is defined as Equation ( 2): (2) Each customer must get the delivery service, which is defined as Equations ( 3) and ( 4): (3) The service for each customer can only be completed by a specific vehicle, which is defined as Equations ( 5) and ( 6): (4) Variables in the range are defined as Equations ( 7) and ( 8): (5) Time window constraints are defined as Equations ( 9) and ( 10): The objective function of the mathematical model consists of three parts, including toll costs, waiting time costs and delay fine costs.The goal is to get the minimum value of the total costs.In addition, set the parameter c 1 and c 2 , we can get the model of the corresponding vehicle routing problem with hard time windows and the model with soft time windows.If c 2 is infinity, we can get the model of the vehicle routing problem with hard time windows.At the same time, the model also ensures that the vehicle route does not have the local sub-loop.Besides, the model set the fine costs when the vehicle arrives earlier or later.To get the least total costs, the solution has to avoid or minimize the fine costs.

HPSO
HPSO is based on PSO.To improve the quality of the particle in the early stage, we use the global optimal particle by COA to re-initialize the particle swarm.At last, we judge the premature convergence and use COA and Gaussian mutation to avoid falling into the local optimal state.Algorithm 1 is shown as follows: Algorithm 1 HPSO Input: maximum number of iterations T Output: the optimal solution 1: Randomly initialize the swarm particle encoding, then decode the particle and get the fitness value of the particle' note F lag = true.(Algorithm 2) 2: When getting the feasible particle solution, use the global optimal particle by COA to re-initialize the particle swarm and get the new particle swarm including n particles; note F lag = true.Go to Step 3, otherwise go to Step 4. (Algorithm 3) 3: Judge the degree of the premature convergence: assume that the average optimal fitness value of the particle in history has not changed for N rounds; if N 1 N N 0 , optimize the global optimal particle by COA; if N > N 1 and the probability of particle randomly generated is P > P 0 , convert the particle by Gaussian Mutation.N 1 , N 0 and P 0 are predetermined constants.6: Stop updating particles; quit.

The Encoding and Decoding of the Particles
This paper proposes an encoding method.It assumes that there are K vehicles, and L customers.It constructs particles in (L + 2 * K)-dimensional space, in which each dimension is described by a real number.The front L-dimensions are the total number of the customers, and the later 2 * K-dimensions are the position of the K vehicles in the Descartes figure.2: Construct V , the matrix of vehicle priority service for the customer: 2.1 Get the position of the vehicle from the position of the particle.For k vehicles, Based on particles encoded, we decode the position of the particles to "the list of customer priority service" and "the matrix of vehicle priority service for the customer".According to them, we construct the vehicle routes.The process is shown as Figure 1 and described as follows: (1) Use the front L-dimensions of the position of the particle to construct the list of customer priority service, which represents the order of the service for customers.The sorting rules of the list are to sort the values of the front L-dimensions in ascending order and then regard the corresponding subscript value as the service order number of each customer.(2) Decode the later 2 * K-dimensions to the position of the vehicles, and construct the matrix of vehicle priority service for customer.This matrix is based on the distance between the position of vehicles and the position of customers.It refers to the priority of each customer receiving the service.The closer the distance, the higher the priority.(3) Build the vehicle routes: according to the list and the matrix, the customers are inserted into the corresponding vehicle route.When it is completed, the set of vehicle routes will be constructed.
The whole process is shown as Algorithm (2).

Insertion Heuristic Algorithm
In the process of particle decoding to construct the vehicle route, we need the operation of insertion.This paper proposes the Insertion Heuristic Algorithm as follows: If after inserting customer u, it meets Equation ( 11), we call p as the best position of u in R. C(i, u, j) is the new route that meets the constraints after inserting u.So, it is necessary to change the start service time of customer (i p , . . ., i m ) and to check whether it meets the constraints.
According to the characteristics of VRPTWs goal, specific heuristic conditions are shown in Equations ( 12) and (13): In these formulas, d iu is the distance between customer i and customer u. µ is a constant.b j is the start time for serving customer j before inserting customer u. b ju is the start time after insertion.C 11 (i, u, j) is the increased traveling distance from the original route after insertion.C 12 (i, u, j) is the delay time serving customer j after insertion.
Considering VRPTW, it is impractical to use the insert method of shortest distance or least time cost.This is because VRPTW is limited by the constraints of time and route.After all customers are inserted by this method, the route will not be the shortest.Even the route cannot meet the constraints after several rounds of insertion.Taken as a compromise, this paper proposes the heuristic insertion conditions that are shown in Equation ( 14) and (15): In these formulas, α 2 , α 2 , λ are the weighting parameters.C 1 (i, u, j) is the whole cost and C 2 (i, u, j) is the whole income.The heuristic insertion conditions try to construct part of the route instead of the entire route to maximize the costs of the new route.For example, when u = α 1 = λ = 1 and α 2 = 0, c 2 (i, u, j) is the reduced distance after inserting customer u.

COA
The basic idea of COA is: Firstly, because the values of the chaotic variable have to be in the range of the limited values, it is necessary to convert the variables.Then, optimize the converted variables according to chaotic rules.Lastly, convert the optimization variables into the original space.Relative to the irregular and aimless search operation in the problem space, chaos has many advantages, which is determined by its characteristics.In addition, chaos can avoid the problem of falling into local optimization that often appears in the evolutionary algorithm and maintains population diversity to improve the global search ability of the algorithm.
In general, there are a lot of methods to generate chaos, such as Lozi's mapping [36], logistic mapping, Ikeda mapping [37], Henon mapping and Tent mapping [38].Among them, Tent mapping (also known as a tent mapping) is a kind of piecewise linear one-dimensional mapping, with uniform power spectral density and probability density and the ideal related characteristic.The mathematical expression is shown in Equation ( 16): From Equation ( 14), we can see that Tent mapping is used to make the chaos process of the float on the interval [0, 1], and the binary number of the fractional part does an unsigned left shift.These kinds of mapping characteristics make full use of the computer's shift operation characteristics (the processor has a shift operation register, so the operation in the computer processing is highly efficient), so it's suitable for computing a large data sequence.The whole process is shown as Algorithm 3: Algorithm 3 Tent mapping algorithm Input: particle position vector Output: chaos transformed particle position vector 1: Do chaos optimization of X(x 1 , x 2 , . . ., x D ) (the position of the particle).According to the mapping relationship, z = (x − x min )/(x max − x min ), map x i (i = 1, 2, . . ., D) into [-1,1] and get z i (i = 1, 2, . . ., D).

Premature Convergence Judgment
In the process of constructing vehicle routes, due to the limits of the vehicles, the time window for serving the customers, even the insertion algorithm, there are particles that cannot be decoded, that is, the position of the particle is invalid and cannot be used as a feasible solution.
Based on the above, assume f i is the optimal fitness value of particle i in history.The optimal position of i in history is a feasible solution.We call this particle the "effective particle".
In Equation ( 17), M is the number of effective particles, and f avg is the average of all optimal fitness values.In the iteration process, if f avg has not changed many times, such as 20 or 50, it is necessary to take measures to increase the diversity of the particle swarm.

Experimental Results and Comparisons
This section is devoted to the performance evaluation of HPSO for VRPTW.The adopted benchmarks and parameter setting will be described in Section 5.1.Section 5.2 describes the parameters analysis of the Insertion Heuristic Algorithm of HPSO.Section 5.3 details the comparison result of the performance of HPSO by different problem sets.Section 5.4 describes the comparison results of HPSO and PSO.Section 5.5 describes the comparison results of HPSO with other state-of-the-art methods.Section 5.6 describes a real-world numerical example of VRPTW.

Benchmark Description and Parameter Setting
In order to compare the performance of our algorithms on VRPTW, the basic data of our testing problems adopt Solomon's benchmark problems [39].Solomon's benchmark problems focus on the factors affecting the route and the scheduling algorithm (including geographical location, the number of customers in each route, the ratio of customers with time windows and the density and flexibility of time windows).They are divided into six types of data sets: R1, C1, RC1, R2, C2 and RC2.There are two goals to solve the data sets: (1) use the least vehicles, that ism the route that is least; (2) the total length of all routes that is shortest.The parameters of HPSO are set as Table 1.We use linear decreasing optimization for inertia factor, w, in the iteration.

Parameter value
population size 100 Number of iterations "t" 1,000 Minimum of the position X min 0 Experiments are carried out under the configuration of Windows VistaTM Home Premium, with AMD Turion (tm) 64 X2 Mobile Technology TL-62 2.10GHz and 2GB RAM.We realize the CGAMO on the platform, MyEclipse 6.6.The solution of the main problem invokes the linprog function of Matlab7.0.

Parameters Analysis of Insertion Heuristic Algorithm
Section 4.2 describes an Insertion Heuristic Algorithm of HPSO proposed in this paper.This section will make some experiments to analyze the parameters of the Insertion Heuristic Algorithm.The Insertion Heuristic Algorithm plays a key role in the process of constructing the vehicle route.It involves four parameters: α 1 , α 2 ,µ and λ.Among them,α 1 + α 2 = 1 and α 1 affects the increasing distance after inserting the customer, and α 2 affects the delay time after insertion.The six types of data sets have their own characteristics for all factors from the geographical location to the required time windows.In order to unify the results, the number of the customers is 100, µ = 0.5,λ = 0.5, α 1 is 0.11, 0.21, . . ., 0.91 in turn and α 2 is 0.890.79...0.09 in turn.Select the value and other parameters for each data set; execute in turn.The result is shown as Figure 2. Based on these results of Figure 2, we can make the following conclusion on the parameter analysis of the Insertion Heuristic Algorithm.
(1) The parameters can hardly affect the results of C101 and C102.That is because their customer locations are concentrated and the time windows are [0,1,236] and [0,3,390]; it is relatively large.(2) On the other hand, R101 and R102 are largely affected by the parameters.Because their customer locations are a random distribution and the time windows are [0,230] and [0,1,000], it is relatively small.RC101 and RC201 are between the above.In the process of executing, we must set the weights of the parameters according to the data characteristics so that we get the better solution.(e) (f)

HPSO Performance Comparison by Different Problem Sets
In order to evaluate the proposed HPSO on large-scale problems, we tested on a set of 56 Solomon's test cases.The aim of the experiments is performance comparison of HPSO for different problem sets in different iteration times, which can test HPSO convergence speed and sensibility to the problem scale.The performance comparison results are shown in Tables 2-4.Tables 2-4 compare the performance produced by the proposed HPSO in the set of C1, C2, R1, R2, RC1 and RC2 for 100, 200 and 500 iterations.The rate in Tables 2-4 expresses the ratio of the difference value of the distance between the best known result and the different iteration time results.Based on these results of Tables 2-4, we can make the following conclusion on the HPSO performance of different problem set for solving the VRPTW.
(1) In a set of C, the coverage speed to a good solution is fast.It gets an approximately best result by 100 iteration times, many problems of C coverage to the best known results by 200 iteration times and all problems of C coverage to the best known results by 500 iteration times.We can see that HPSO has a better performance on problem set C than Tavakkoli-Moghaddam et al. [24].(2) In a set of R, the average rate of 100 iteration times is nearly 15% for each problem, the average rate of 200 iteration times is nearly 4% and almost all problems of C coverage to the best known results by 500 iteration times.We can see that HPSO has a better performance on problem set R than Tavakkoli-Moghaddam et al. [24].(3) In a set of RC, the average rate of 100 iteration times is nearly 13% for each problem, the average rate of 200 iteration times is nearly 6% and almost all problems of C coverage to the best known results by 500 iteration times, except RC101, RC102, RC105, RC107, RC201, RC202, RC203 and RC207.We can see that HPSO has a better performance on problem set RC than Tavakkoli-Moghaddam et al. [24].

Comparison Results of HPSO and PSO
In HPSO, in order to avoid falling into the local optimum, this paper uses COA and Gaussian mutation.In order to add the number of particles in the early stage and to accelerate the searching speed in solution space, this paper uses COA to reinitialize the particle swarm.To verify results, this section makes experiments to compare the HPSO and basic PSO and tests the effect of COA in HPSO.The results are shown in Table 5.Table 5 compares the performance produced by the proposed HPSO and basic PSO in the set C1, R1 and RC1.The run-time in Table 5 expresses the solving time of problems by HPSO or PSO.The rate in Table 5 expresses the ration of different value of distance between the best known result and HPSO or PSO results.If the rate value is negative, it means that the distance value of HPSO or PSO is shorter than the best known results.If the rate value is a positive number, it means that the distance value of HPSO or PSO is bigger than the best known results.Based on these results of Table 5, we can make the following conclusion on the performance of HPSO and PSO on solving the VRPTW.
(1) HPSO can get all the best known results, and sometimes, it gets better results than the best known result, such as for problems R101, R103, R105, R106 and RC101; while basic PSO can hardly cover the best known results for almost all problems of C1, R1 and RC1.We can see that HPSO has better performance than basic PSO.(2) HPSO covers the best known results at nearly two minutes for R1, one minute for C1 and at most three minutes for RC1.PSO covers its best results at least three minutes for R1, two minutes for C1 and four minutes for RC1.We can see that COA can effectively accelerate the searching speed in solution space, and the HPSO has better performance coverage speed than basic PSO.Because the number of the vehicles and routes in different data sets are different, this paper only chooses R201, RC101 and RC201 to make an experiment.The comparisons of the number and the average optimal fitness value in history of the effective particle are shown as Figure 3. From the data of Figure 3, it is easy to get HPSO better than PSO.Although, in RC101, the solution of HPSO is worse than PSO.The number of effective particles and the average optimal fitness value in history of HPSO are better than PSO.This shows that HPSO is more effective.

Comparison Result of HPSO with Other State-of-the-art Methods
In order to compare the efficiency of different algorithms to VRPTW, The state-of-the-art method of Goal Programming and Genetic Algorithm (GPGA) [28], Improved Large Neighborhood Search (I-LNS) [22] and Ant Colony Optimization Hybridized with Insertion Heuristics (ACO-IH) [25] are used to compare HPSO.The results of comparison experiments of VRPTW are shown in Tables 6-8.The rate in tables 6-8 expresses the ratio of the difference value of the distance between HPSO and other state-of-the-art methods.If the rate value is negative, it means that the distance value of HPSO is shorter than other state-of-the-art methods.If the rate value is a positive number, it means that the distance value of HPSO is bigger than other state-of-the-art methods.Based on these results of Tables 6-8, we can make the following conclusion on the effectiveness of the HPSO approach for solving the VRPTW.

A Real-World Numerical Example
A real-world numerical example is provided here to illustrate the effectiveness of HPSO.The data of Yunda Expressage Company of Wuhan City is used to test the sensitivity of the algorithm to the specimen data.One-day data of the Yunda Expressage Company of Wuhan City of Hongshang district is transformed as a Solomon problem set format, which is shown in Table 9. Yunda Expressage Company of Wuhan City has 10 vehicles.They have 60 delivery nodes in Hongshang district.The delivery nodes can be defined as the customers of Solomon problem data set.We use HPSO to solve this real-world numerical example, and the computing result is shown in Table 10.We use other state-of-the-art methods to solve this example, and the comparison result of HPSO with other state-of-the-art methods is shown in Table 11.Based on these results of Tables 9-11, we can make the following conclusion on the performance of the HPSO and other state-of-the-art methods on solving a real-world numerical example.Table 10.The computing result of the real-word numerical example by HPSO.
Vehicle ID The Amount of delivery node delivery node (3) In conclusion, HPSO achieves better performance than other state-of-the-art methods on solving real-world VRPTW.

Conclusions
VRPTW has drawn great interest in the field of vehicle routing.The research area of the problem also expanded to water carriage, aviation, logistics, industrial management and communication.Most researchers use heuristic evolutionary algorithms to solve it.However, the principles of these heuristic algorithms are mostly complex and need more optimization work.This paper proposes a hybrid PSO algorithm, combining the chaos optimization algorithm and the heuristic algorithm, with the purpose of improving the number of effective particles and avoiding falling into local optimization.
The experiments mainly analyze the influence of the insertion heuristic algorithm parameters and comparison of HPSO with state-of-the-art methods.From the experiments, it concludes that: (1) it is necessary to set the weights for the parameters according to the characteristics of customer locations and time windows.When there are feasible solutions, using the Chaos Optimization Algorithm to re-initialize the particle swarm can speed up the convergence rate.Using CAO and Gaussian mutation can further avoid falling into the local optimum.(2) Using the Solomon benchmark, HPSO achieves much better performance than the other state-ofthe-art methods in converging speed and solutions.(3) Using a real-world numerical example, HPSO covers a better solution than the other state-of-theart methods, and its coverage iteration times and run-time times are also the least.
Although HPSO for VRPTW obtains satisfactory achievements, there are still some spaces for improvement, such as how to effectively construct the initial solution and how to precisely judge the local optimal solution.These would be explored in our future work.gratefully acknowledge the helpful comments and suggestions of the reviewers, which have improved the presentation.

Figure 1 .
Figure 1.Particle decoding and the set of vehicle routes.

Algorithm 2
Decoding the algorithm of HPSO Input: the position of the particle Output: the set of vehicle routing 1: Construct U , the set of customer priority service: 1.1 set S = {1, 2, . . ., L}, each element of S corresponds to the customer's number.And U = ∅; 1.2 From S, choose customer C, which corresponds to the minimum value of the front L-dimensions in the position of the particle, namely x ic = min x ie (c ∈ S); 1.3 Insert c after u (u is the last element of U ); 1.4 Delete c from S 1.5 If S is not empty, go to Step 1.2.

Figure 3 .
Figure 3.The comparisons of the number and the average optimal fitness value in the history of the effective particle.(a) the number of effective particles of R201; (b) the average optimal fitness value in history of R201; (c) the number of effective particles of RC101; (d) the average optimal fitness value in history of RC101; (e) the number of effective particles of RC201; (f) the average optimal fitness value in history of RC201.
Insert v into the end of V i ; 2.2.5 Delete v from G; 2.2.6If G is not empty, go to Step 2.2.3.According to the Insertion Heuristic Algorithm, insert customer c into the route so that the total costs of the new route is least.This route is regarded as the candidate route; 3.2.4Check the constraints of the candidate route; 3.2.5If the candidate route meets the constraints, update R v by the candidate route, optimize R v by 2-opt algorithm, otherwise go to Step 3.3; 3.2.6If j = K, go to Step 3.3; otherwise, if j = j + 1, go to Step 3.2.2;3.3 If h = L, quit; otherwise, if h = h + 1, go to Step 3.2.
2, . . ., L}, j ∈ {1, 2, . . ., K}) ; 2.2 For each customer i: 2.2.1 Calculate the distance from the positions of all vehicles, and θ j is the distance between customer i and vehicle J; 2.2.2 set G = {1, 2, . . ., K}; each element corresponds to the number of vehicles, V i = ∅; 2.2.3 From G, choose vehicle v, which has the shortest distance from customer i, namely: θ v = min θ e , e ∈ G; 2.2.4 3: Construct R, the set of the vehicle routes: 3.1 Assume h = 1; 3.2 Add customers to the route: 3.2.1 Assume c = U h (the current customer to be served for) and variable j = 1; 3.2.2v = V c,j (the vehicle v serves for customer c); 3.2.3

Table 2 .
Comparison of the performance of the proposed HPSO for a set C.

Table 3 .
Comparison of the performance of the proposed HPSO for a set R.

Table 4 .
Comparison of the performance of the proposed HPSO for a set RC.

Table 5 .
Comparison of the performance of HPSO and PSO.

Table 6 .
Comparison results of different methods for vehicle routing problem with time windows (VRPTW) of R1.

Table 7 .
Comparison results of different methods for VRPTW of C1.The solution quality computed by HPSO can approximate the best known results well.If the travel distances and the number of used vehicles are taken into account seriously, the HPSO can almost get the best results of C1, R1 and RC1.(2) Comparing with the results generated by GPGA, I-LNS and ACO-IH, the results of HPSO are almost better than them.The average CPU times of HPSO are less than 3 min.The results are hard to systematically compare with other state-of-the-art methods, because there are big differences in the run environments.

Table 8 .
Comparison results of different methods for VRPTW of RC1.

Table 9 .
One-day data of Yunda expressage company of Wuhan city.

Table 11 .
The result of comparison experiments of the real-world numerical example.