Improved artificial bee colony algorithm for vehicle routing problem with time windows

This paper investigates a well-known complex combinatorial problem known as the vehicle routing problem with time windows (VRPTW). Unlike the standard vehicle routing problem, each customer in the VRPTW is served within a given time constraint. This paper solves the VRPTW using an improved artificial bee colony (IABC) algorithm. The performance of this algorithm is improved by a local optimization based on a crossover operation and a scanning strategy. Finally, the effectiveness of the IABC is evaluated on some well-known benchmarks. The results demonstrate the power of IABC algorithm in solving the VRPTW.


Introduction
The Vehicle Routing Problem (VRP) was first described by Dantzig and Ramser in 1959 [1]. The problem involves determining the delivery routes, which start and end at the depot and also serve all the customers. Each customer is visited once and the total demand of all customers in a specific route can't exceed the capability of the vehicle. The goal of the VRP is to minimize the total cost, where cost can be defined as distance or time [2]. In 1981, Lenstra and Rinnooy Kan [3] proved that the VRP is NP-hard. Its main drawback is the well-known curse of dimensionality. Constraint programming is considered as an efficient method for giving the optimal solution [4] but the time to find the optimal is prohibitive in large problems. Some studies improved constraint programming to solve the VRP: Backer et al. [5] developed a method using iterative improvement techniques within a Constraint Programming framework to get a good solution under a relatively short computing time; Ozfirat and Ozkarahan [6] introduced an algorithm which decomposing the problem into smaller scale ones firstly and then solve it by constraint programming; Guimarans et al. [7] proposed a method of combining Constraint Programming, Probabilistic Algorithms and Lagrangian Relaxation.
The vehicle routing problem with time windows (VRPTW) is a variant of the VRP that has an additional time window constraint. In the VRPTW, a fleet of vehicles set off from a depot to serve a number of customers, which have various demands and specific time windows. The objective of the VRPTW includes maximizing the sum of the on-time delivery probabilities to customers, minimizing the expected total cost, and some others [8]. The VRPTW, as an extension of the VRP, has also been proved to be an NP-hard problem [9]. The time window require the time constraint and the traffic speed prediction is very important in the delivery routes in paper. Finally, a list of the notation used in the proposed algorithm is presented in the S1 Table. Problem description A typical VRPTW specifies: the service time constraint, demand of each customer and the capacity of each vehicle. The routes must be designed so that each point is visited once by one vehicle and all routes start and end at the depot. The total demand of all points on one particular route cannot exceed the capacity of the vehicle. The service time of each customer must be in the time constraints of the customers.
The VRPTW model in this paper is represented as a weighted graph, G = (V, E), where V = {1,2,3,Á Á Á,N} represents the vertex set and N is the number of nodes in the graph, E represents the edge set; c ij (c ij > 0, c ii = 1, i, j 2 V) represents the distance between a pair of vertices. The notation used in this model are listed as follows in Table 1.
For the sake of convenience, several assumptions are made in this paper: (1) Customer demands can be split and the maximum possible demand of each customer is smaller than the vehicle capacity; (2) a customer's actual demand is only re-vealed when the vehicle arrives at that customer's location and the customer's earliest acceptable service time (i.e. the lower bound of the time window) begins; (3) service times at customer locations are neglected for clarity. (4) We assume t ij = c ij .
The model is described as follows: The objective is to design a network that can minimize the total travel length under the conditions of all constraints. The constraints are defined as follows: (2) The maximum number of routes constraint: to specify the maximum P routes going out of the depot.
(3) Travel constraint: ensures every route starts and ends at the central depot.
(4)-(5) Service constraints: to assume that every customer node can be visited exactly once by one vehicle.
(6) The capacity constraint: to assume the quantity of a vehicle cannot exceed its capacity.
(7) The maximum travel time constraint: to limit the travel time to less than the maximum. As the ABC algorithm is easily trapped in local optima, we present an improved artificial bee colony (IABC) algorithm to avoid this limitation in our investigation of the vehicle routing problem.

Improved artificial bee colony algorithm
Artificial bee colony (ABC) algorithm is a relatively new optimization method, which simulates the acts of the real bees looking for nectar in nature. It has been successfully applied to solve a series of complex combinational optimization problems. In order to overcome some limitations of the standard ABC algorithm, an improved ABC algorithm is presented in this paper. And in this section, we firstly introduce the fundamental rules of ABC algorithm, then two local optimization methods that prevent the algorithm trapped in local optimal are presented.

The fundamental rules of ABC
An ABC consists of three main populations: employed bees (or called leaders), onlookers and scouts. Employed bees are responsible for the mining process and gathering required information. Scouts undertake the exploitation of new food sources. Whereas onlookers are in charge of sharing information with employed bees and scouts by communicating in the dance area. The three roles are interchangeable among the individuals. For example, in Fig 1, a bee initially has two choices; it can scout for food sources (role S in Fig 1), or follow the leader as an onlooker after witnessing the leader's dance (role R in Fig 1). When a scout has found the food source and begins collecting nectar, it becomes a leader. This bee then returns to the honeycomb and unloads the food. At this time, it has three choices; abandon the food source and become a scout or onlooker (blue line in Fig 1), continue the leader role and dance to recruit more bees for food sourcing (red line in Fig 1), or collect nectar from the food source without recruiting other bees (green line in Fig 1) (Liu et al. [34]).
In the ABC algorithm, each alternative solution is called a nectar source. The nectar collecting process is similar to the process of searching for the optimum solution. The quality of each nectar source is determined by its corresponding fitness value. The number of solutions (S) equals the number of leaders or onlookers. The position of food source i is represented by a D- dimensional coordinate vector u i (u i1 ,u i2 ,Á Á Á,u id ,Á Á Á,u iD ) T 2 {u i |1 i S,i is integer}. So the ABC algorithm first obtains an original population with S solutions, where each solution u i (i = 1,2,Á Á Á,S) is a D-dimensional vector. Then during each iteration, each leader finds a new food source by conducting a cyclic neighborhood search on its originally old food source. The cycle number is denoted by C(whereC = 1,2,Á Á Á,M). After the new food position is generated, its corresponding fitness value need to be evaluated.
If the nectar found at the present food source (solution) is of higher quality (fitness) than the former nectar, the old food source is replaced by the new source; otherwise, the position of the old food source is unchanged. As the search approaches the best solution, the neighborhood range becomes smaller. Once all leaders have completed their search work, they return to the dance area and deliver the nectar information of the food source to onlookers. Then, each onlooker selected a food source according to the roulette wheel selection method. The higher the nectar content of a food source is, the higher the probability that that food source will be chosen by onlookers.
Having chosen a food source, onlookers either select a neighborhood search around the food source, or keep the current solution. Furthermore, the ABC algorithm imposes a parameter limit on the number of updates. If a solution is not improved after consecutive cycles, it has sunk into a local optimum and should be abandoned. The corresponding leader then converts into a scout and searches for new food source randomly. After the scout finds a new food source, the scout becomes a leader again. After each leader is assigned to a food source, another iteration of the ABC algorithm begins. The whole process is repeated until a stopping condition is met. For more details about the ABC algorithm, interested readers should refer to Karaboga [20,35] and Szeto [29].
Based on the preceding analysis, four selection processes are specified in the ABC algorithm; first is the global selection process, in which onlookers seek a good food source; second is the local selection process, by which leaders and onlookers search the neighborhood; third is a greedy selection process in which all artificial bees judge and compare the old and new food sources, then retrieve the best solution; and fourth is a stochastic selection process that identifies a food source.
In the application of vehicle routing problem, the bees construct paths in different ways depending on their roles. Leaders only retrace the last path of the iterative procedure without changing the situation. The scouts and followers choose the next nectar source as follows: where p k ij ðtÞ represents the probability that bee k travels from node i to j in iteration t. A k i ¼ F À Tabu k denotes the available nodes for bee k (i.e., the points that satisfy all constraints but have not been previously visited by the bee). Tabu k is the tabu list of bee k, which prevents the nodes already visited by bee k from being repeatedly visited. In some applications (i.e., router optimization), it is also convenient that the bee can return through the original path. s ij (t) and f ij (t) represent the bee colony and heuristic information, respectively, while α and β are the weight coefficients.
The control information differs between scouts and followers. For a scout, the control information is given by where l is the number of unvisited nodes. For the follower, the information is given by Abandon the path of scout Optional paths contain the scout path 1 l Abandon the path of scout where r represents the guiding strength of the leader, and l represents the number of unvisited nodes, as before. However, the above ABC easily becomes trapped in local optima. Moreover, a local optimum is difficult to escape because of the limited diversity of solutions, and the current solution may contain bad information that is passed to the next solution. To overcome these problems, we introduce the crossover operation and a scanning strategy for local optimization.

Local optimization
Due to the weaknesses of the standard ABC, it is necessary to expand the solution space and prevent bad information entering the following search. In this section, crossover operation is first introduced to expand the range of the solution search, then, scanning strategy is developed to eliminate the effect of bad information.
Crossover operation. In local optimization, the crossover operation effectively prevents the algorithm from trapping in a local optimum and reaches further solutions in the search space [36]. In this paper, the performance of ABC is improved by a simple crossover operation, which proceeds as follows: Step 1. Randomly select two paths from the solutions. Among the points in each path, randomly select one crossover point. For example, suppose that clients c 8 and c 2 are selected, as shown in Fig 2. And by exchanging c 8 and c 2 , we obtain a new solution (Fig 3).
Step 2. Because the clients are selected randomly, the solution is likely to be inferior. Therefore, the local optimization of the solution must be improved by a local search algorithm. In this paper, the solution is improved by a 2-opt algorithm (Yu et al., 2012).
Similar to the genetic algorithm, we introduce a crossover rate (p m ) that decides whether a solution needs a crossover. If p m is too large, it will slow the convergence speed of the algorithm. In the early stages, the optimization should usually be performed over the largest possible search space.
Over time, the accuracy of the optimization will gradually improve. A large jump in accuracy will impede the convergence. Therefore, we adopt an adaptive algorithm that computes the crossover rate as follows: where p min m and p max m represent the minimum and maximum crossover rates, respectively, p min m ¼ 1=H (H is the number of clients in the net), p max m ¼ 1=N (N is the number of distribution centers), and T and t denote the maximum and current iteration generations, respectively.
Scanning strategy. If two paths in an optimization program intersect at one point, the local optimization can be effectively performed by eliminating that crossover point [37,38]. Most researches find the crossover point by an inefficient enumeration method (that is, by comparing all paths). Especially in large-scale VRPs, this algorithm requires a long computing time. In this paper, the local optimization program also reduces the crossover phenomenon, but quickly finds the crossover phenomenon by a scanning method [39]. This rapid technique greatly improves the efficiency of the local optimization algorithm. The scanning method proceeds by the following steps: Step 1 Initialization. The proposed scanning strategy considers the angle between the customer point and central depot; therefore, the locations of these points must be known. To this end, the scanning strategy establishes a coordinate system with the origin set to the central depot coordinates, then draws a line from the origin to the customer point in two-dimensional space. As an example, Fig 4 shows the spatial distribution of customer points around the central depot. Each customer is initialized to note its location and corresponding paths. For example, if customer c 4 is located 115.3˚from the horizontal and the selected path is 1, the (location, path) of that customer is initialized to (115.3 o , 1).
Step 2 Search for the potential crossover paths based on the scanning strategy. The customer points are sorted by their angular values, and their path numbers are recorded as shown in Fig 5. Starting from the customer with the smallest angle, all paths containing the same customer are scanned. The path number can be changed (e.g., from customer c 3 to c 2 , the path number changes from 1 to 3). However, the new path may have been already scanned (e.g., from customer c 2 to c 4 , the path number changes from 3 to 1), a phenomenon called fallback.  If the customer points in one path are continuous and never encounter the fallback phenomenon, that path does not cross any other path (i.e., path 2 in the present example). The fallback phenomenon indicates a possible intersection. The possibly intersecting path is then checked in the tabu list. If the path resides in the tabu list, the algorithm determines whether there is a crossover point between two paths (Step 3). If the result is negative, the scanning continues. If the scan reaches the final customer point without encountering the fallback phenomenon, the algorithm advances to Step 6.
Step 3 Compare all sections between two paths to judge whether they intersect. The judgement of an intersection between two paths is as follows: Assume that there exist two paths ab and cd, in which path ab has endpoint coordinates (x 1 , y 1 ) and (x 2 , y 2 ) and path cd has endpoint coordinates (v 1 , w 1 ) and (v 2 , w 2 ). 1) If x 1 6 ¼ x 2 and v 1 6 ¼ v 2 First, the slopes θ 1 and θ 2 of paths ab and cd are respectively computed as follows: If θ 1 = θ 2 , the two paths cannot intersect (Fig 6 (A)); otherwise, they may or may not intersect. In the second case, the existence of the intersection is checked by calculating the coordinates of the potential intersection. The abscissa of the potential intersection is: where, b 1 and b 2 are two judgment parameters. If z lies between x 1 , x 2 , and v 1 , v 2 , then the intersection exists (Fig 6 (B)); otherwise, the intersection is absent (Fig 6(C)).
2) If x 1 = x 2 or v 1 = v 2 If x 1 = x 2 and v 1 = v 2 , the two paths cannot intersect (Fig 6(D)). If x 1 = x 2 and v 1 6 ¼ v 2 , the vertical axis of the potential intersection Z is computed by inserting z = x 1 into the linear equation of path cd. If Z lies between y 1 and y 2 , the two paths intersect (Fig 6(E)); otherwise, they do not intersect (Fig 6(F)). The case x 1 6 ¼ x 2 and v 1 6 ¼ v 2 is treated similarly.
If the two paths intersect, the scanning strategy proceeds to Step 4; otherwise, the traversal continues. If no intersection has been found after the traversal, the strategy proceeds to Step 5, which adds these two paths to the tabu list to avoid repeated judgements in subsequent iterations.
Step 4 The generation of new tours is ensured by the 2-opt algorithm. If the two new tours can improve the solution, they replace the old tours; otherwise, the strategy returns to Step 2.
Step 5 Judgments are terminated if no feasible crossover is found or if the algorithm has reached the specified iteration limit; otherwise, the strategy returns to Step 2.

Determination of scouter and leader proportions
The proportions of scouters and leaders in the colony are denoted as PScouter and Pleader respectively. These proportions largely determine the performance of the ABC algorithm. In particular, PScouter reflects the strength of the randomness in the path searching; the larger the PScouter, the less likely that the colony will select a previously visited path. In this way, the colony can extend the solution space. However, the strengthened search randomness tends to reduce the convergence speed. Meanwhile, Pleader reflects the strength of the convergence factor in the path searching; the larger the Pleader, the more likely that the colony will select the best previously visited path. Thus, a large Pleader improves the capability of maintaining the current optimal solution, but increases the chance of falling into a local optimum. In summary, Pscouter and Pleader exert opposite effects on the convergence rate and local optima trapping. Thus, Pscouter and Pleader are mutually restricted such that one cannot be very much larger or smaller than the other. Improved artificial bee colony algorithm for vehicle routing problem with time windows Properly balancing the PScouter and Pleader values is crucial for increasing the convergence speed and enhancing the ability to find the global optimal solution. Increasing the convergence speed must not unduly compromise the capability of finding the global optimum, which itself requires a good random search method. In this paper, the appropriate parameter combination is determined in an experiment.
The appropriate range of PScouter and Pleader was found in a simulation experiment conducted in R2-01, which can be found in the literature [20]. Assuming 85 bees, the other parameters, whose values were suggested in the literature [40], were employed as follows: maximum iterative time (MIT) = 200; transition rate of onlookers (γ) = 0.9, cycle limit = 50, α = 1, β = 4. PScouter was varied as 0, 0.1, and 0.3, whereas Pleader was varied as 0.25, 0.33, 0.5, 0.75, and 1. The averages and optimal values were obtained from 10 independent computations. The results are listed in Table 2.
From Table 2, we find that (PScouter, Pleader) = (0.3, 0.33) yielded the best performance.  In addition, we randomly created food sources at the initialization step. The number of candidate food sources was increased by a neighborhood moving algorithm. In the experiments, the number of random original food sources m affected the convergence speed and accuracy of the algorithm. This simulation experiment was carried out 20 times, setting m from 1 to 12, and maintaining the other parameters at their previous values (stated above). And Fig 7 plots the convergence when the problem was solved by the IABC algorithm with m = 1,4 and 12. It is concluded that the solution was less accurate at m = 1 or m = 12 than at m = 4. Precise solutions are obtained for m 2 [3,8]. Thus, m = 5 was adopted in the remaining experiments.

Numerical analysis
To check the performance of the proposed IABC algorithm, we employed some well-known benchmarks (Solomon) in the past literatures [9,[40][41][42][43][44][45][46][47][48][49][50]. The results from IABC and the bestknown solutions are listed in Table 3. Category C represents clustered problem set. Category R is randomly generated geographical distribution of customer nodes. Category RC is a mix of random and clustered structures. The IABC was comparable to the known optimal solution (bold font in Table 3) in 29 experiments, and surpassed the best-known solutions in 4 experiments (C101, C105, C109, R108). Also, in 10 results, the number of vehicles was clearly lower in IABC than in the other algorithms and the optimal solutions were almost equal, indicating that the IABC lowers the vehicle cost. And the average deviation from best known solutions are 1.64%. The solutions calculated in set C, in which customers are clustered in groups, are encouraging. This may be because the introduction of crossover operation diversifies the bee colony, widens the search space and prevents the algorithm from trapped in local optimization. For most heuristics, it is difficult to find the best solution in local and is easily trapped in local optimal because the similarities between customers in instance set C. The crossover operation with bigger crossover rate can expand solution space and get relatively good solution.
Considering both number of vehicles and the total distance, IABC presents as an effective method for solving theoretical and practical VRPTWs. The simulation results reveal three main features of the IABC algorithm: 1) the structure of the algorithm is simple and easily understood. It requires minimal preliminary knowledge and has very strong readability.
2) The solution is highly stable. In an unchanging running environment, multiple runs of the algorithm yield the same outcome.
3) The algorithm runs in parallel; that is, every bee simultaneously begins its working cycle from different missions, so the optimal solution is rapidly found in a large solution space. Therefore, IABC is effective for VRPTW problems. Furthermore, to test the performances of the crossover operation and scanning strategy, we separately introduced these strategies to the standard ABC. The ABCs with the crossover operation and scanning strategy are denoted as ABC-C and ABC-S, respectively. The results are shown in Table 4. Table 4 shows that the ABC-C performs very similarly to IABC, whereas ABC-S obtains the best-known solution in some instances. The quality of the total distance and number of vehicles was comparable in IABC and in ABC-C and ABC-S. However, some of the total distances are better in IABC than in ABC-C and ABC-S, reflecting that the crossover operation and scanning strategy collectively prevent the algorithm from trapping in local optima. The crossover operation improves the solutions, but also increases the computation time. In Table 4, ABC-C consumes more runtime than ABC-S, possibly because the crossover operation requires time to search the crossover nodes. Furthermore, IABC integrates the crossover operation and scanning strategy to provide high-quality solutions in less time than ABC-C and ABC-S.
To analyze the convergences of IABC, the data of C108 is computed by our proposed IABC in ten times and the results can be seen in Fig 8. From Fig 8, it can be attained that the total distance decreases fast before the 60th generation, and then it changes smoothly. The least prediction error appears at about 120th generation, and then it almost remains unchanged. It can be also found that the calculation results of the ten times are almost equal. This suggests that the proposed IABC algorithm has a good converge and this algorithm is effective for the vehicle routing problem with time window.

Conclusion
This paper develops an improved ABC algorithm for solving VRPTW, a complicated combinatorial optimization problem. ABC algorithms, which randomly search for the optimal solution, have been gradually adopted in combinatorial optimization, artificial intelligence and other fields. Other uses of ABCs are robots that search for optimal routes, scheduling and services that seek to transport fresh food in the best order of delivery. The positive feedback and coordination of the ABC algorithm are effective for distribution systems. This paper improves the performance of standard ABC by a crossover operation borrowed from the genetic algorithm, and a scanning strategy. The Comparison with best-known solution in classic VRPTW experiment verifies the capability of the IABC algorithm. The comparison of the results of ABC-C, ABC-S, and IABC shows that the incorporation of crossover operation and scanning strategy can improve the performance of ABC for VRPTW. However, the crossover operation may increase the computation time for the convergence value with the reason that it requires time to search the crossover nodes. In addition, good performance of IABC for small instances cannot ensure the validity in large instances. In future work, we will improve the convergence speed and demonstrate the effectiveness of the algorithm in large instances and practical case studies. And the performance of the IABC in artificial intelligence and other fields are also expected. Besides, information and communication technology advances have encouraged the development of advanced traveler information systems (ATIS) [51]. Applying ATIS to this paper should be studied in the future.
Supporting information S1