Bilayer Local Search Enhanced Particle Swarm Optimization for the Capacitated Vehicle Routing Problem

: The classical capacitated vehicle routing problem (CVRP) is a very popular combinatorial optimization problem in the ﬁeld of logistics and supply chain management. Although CVRP has drawn interests of many researchers, no standard way has been established yet to obtain best known solutions for all the different problem sets. We propose an efﬁcient algorithm Bilayer Local Search-based Particle Swarm Optimization (BLS-PSO) along with a novel decoding method to solve CVRP. Decoding method is important to relate the encoded particle position to a feasible CVRP solution. In bilayer local search, one layer of local search is for the whole population in any iteration whereas another one is applied only on the pool of the best particles generated in different generations. Such searching strategies help the BLS-PSO to perform better than the existing proposals by obtaining best known solutions for most of the existing benchmark problems within very reasonable computational time. Computational results also show that the performance achieved by the proposed algorithm outperforms other PSO-based approaches.


Introduction
The huge importance of solving capacitated vehicle routing problem (CVRP) in companies supporting service and logistics for supplying commodities is a strong motivation for many researchers worldwide to solve the CVRP with innovative algorithms.Capacitated Vehicle Routing Problem in general deals with some customer nodes having some sort of pre-specified demands to be served from a hub with a set of vehicles.Although there are many practical constraints in the real world like stochastic customer demand, resource restrictions and so on.There are many pieces of research that introduce practical variations of CVRP [1][2][3][4][5].However, a large portion of practical applications are simple in nature which can be fit to CVRP.The definition of CVRP has been outlined in [6].A vehicle picks some shipment from the depot to serve some of the customers according to their demands until its capacity permits, makes a trip to serve them, and returns to the depot.The vehicles also have additional restrictions on their maximum travelling distance which is the total of distance travelled by the vehicle and the service time it requires to handle the customers it serves.A set of nodes is given on a graph and a hub location h is specified among them.Then a vehicle routing problem can be portrayed as Figure 1, where the nodes represent customers, cities or so on.The target of solving a CVRP is to find minimum total cost, specifically minimum travelling distance, incurred by the vehicles for visiting the customer nodes to satisfy their demands.CVRP is closely related to the classical Travelling Salesperson's Problem (TSP) with a difference that vehicle capacity is limited in CVRP, but not in TSP.This makes CVRP a more complex problem than TSP.As TSP is a well-known NP hard problem, where no algorithm is available to solve TSP by polynomial time when the number of nodes increases, CVRP too is obviously an NP hard problem [7].
Even though an exact algorithm is unlikely to exist to solve CVRP, there have been a few attempts from researchers based on branch and bound technique, set partitioning technique and so on, which can handle CVRP with only few nodes [8].Therefore, many approximate algorithms are proposed to solve CVRP effectively.Although these algorithms never find the optimal, they can find "good" or near-optimal solutions within a reasonable computational time.
Among heuristics, Clark and Wright's saving algorithm is a very early introduction in the area of CVRP [9].This technique is basically a constructive heuristic which builds tours gradually by inserting a node at each time with respect to some saving conditions till no more cluster can be constructed.There are many improvements and variations available in the literature from many researchers [10][11][12]."Cluster-First Route-Second" [13,14] and "Route-First Cluster-Second" [15] are other two well-known algorithms for handling CVRP.The former one builds some clusters of customers then a TSP tour is obtained with those customers to solve CVRP, and the latter one is the technique to make a TSP tour of all customers then cluster them into routes with respect to the given constraints.
Recently many nature-inspired metaheuristics such as Genetic Algorithm (GA) [16,17], Tabu Search (TS) [18][19][20] and Simulated Annealing (SA) [21,22] have become attractive in the field of optimization.Nevertheless, there are many research articles that report their performances in dealing with CVRP [23][24][25][26][27].It is worth mentioning that swarm optimization techniques are getting popular due to their simplicity and performance to handle optimization problems nowadays [28][29][30].Swarmbased metaheuristics are mainly inspired from unique characteristics of some swarms such as birds (bird flocking), fish (fish schooling), ants (foraging), etc.These creatures are not very intelligent individually whereas they show intelligent and complex behavior while working collectively; see, for example, the foraging behavior of bees for honey collection.Swarm-based metaheuristics are also getting remarkable attention from many researchers in dealing optimization problems like CVRP and its variations [31][32][33][34][35][36][37].
In this work, to explain the applicability of the proposed scheme, we adopt particle swarm optimization (PSO) due to its growing applications in the field of optimization [38,39].PSO was first introduced by Kennedy and Eberhart [40] for handling continuous optimization problems.Afterwards, numerous applications are proposed to trick the discrete optimization problems till date.Among them, articles on CVRP are also available in the literature.Quantum individual-based discrete PSO (DPSO) is proposed in [41] where, the quantum particle swarm optimization is introduced in [42].In DPSO, PSO is hybridized with simulated annealing (SA) to obtain a complete CVRP solution.In their approach they cluster the customers with DPSO, and then make a sequence CVRP is closely related to the classical Travelling Salesperson's Problem (TSP) with a difference that vehicle capacity is limited in CVRP, but not in TSP.This makes CVRP a more complex problem than TSP.As TSP is a well-known NP hard problem, where no algorithm is available to solve TSP by polynomial time when the number of nodes increases, CVRP too is obviously an NP hard problem [7].
Even though an exact algorithm is unlikely to exist to solve CVRP, there have been a few attempts from researchers based on branch and bound technique, set partitioning technique and so on, which can handle CVRP with only few nodes [8].Therefore, many approximate algorithms are proposed to solve CVRP effectively.Although these algorithms never find the optimal, they can find "good" or near-optimal solutions within a reasonable computational time.
Among heuristics, Clark and Wright's saving algorithm is a very early introduction in the area of CVRP [9].This technique is basically a constructive heuristic which builds tours gradually by inserting a node at each time with respect to some saving conditions till no more cluster can be constructed.There are many improvements and variations available in the literature from many researchers [10][11][12]."Cluster-First Route-Second" [13,14] and "Route-First Cluster-Second" [15] are other two well-known algorithms for handling CVRP.The former one builds some clusters of customers then a TSP tour is obtained with those customers to solve CVRP, and the latter one is the technique to make a TSP tour of all customers then cluster them into routes with respect to the given constraints.
Recently many nature-inspired metaheuristics such as Genetic Algorithm (GA) [16,17], Tabu Search (TS) [18][19][20] and Simulated Annealing (SA) [21,22] have become attractive in the field of optimization.Nevertheless, there are many research articles that report their performances in dealing with CVRP [23][24][25][26][27].It is worth mentioning that swarm optimization techniques are getting popular due to their simplicity and performance to handle optimization problems nowadays [28][29][30].Swarm-based metaheuristics are mainly inspired from unique characteristics of some swarms such as birds (bird flocking), fish (fish schooling), ants (foraging), etc.These creatures are not very intelligent individually whereas they show intelligent and complex behavior while working collectively; see, for example, the foraging behavior of bees for honey collection.Swarm-based metaheuristics are also getting remarkable attention from many researchers in dealing optimization problems like CVRP and its variations [31][32][33][34][35][36][37].
In this work, to explain the applicability of the proposed scheme, we adopt particle swarm optimization (PSO) due to its growing applications in the field of optimization [38,39].PSO was first introduced by Kennedy and Eberhart [40] for handling continuous optimization problems.Afterwards, numerous applications are proposed to trick the discrete optimization problems till date.Among them, articles on CVRP are also available in the literature.Quantum individual-based discrete PSO (DPSO) is proposed in [41] where, the quantum particle swarm optimization is introduced in [42].In DPSO, PSO is hybridized with simulated annealing (SA) to obtain a complete CVRP solution.In their approach they cluster the customers with DPSO, and then make a sequence of customers utilizing the SA algorithm.
Though they show quality solutions for some benchmark instances, and the approach requires large computational time to reach the best-known solutions (BKS).Kao [43] proposed a combinatorial PSO mainly based on the work of Chen [41].A technique utilizing an advanced version of PSO, namely GLN-PSO in [44], was introduced by Pongchairerks and Kachitvichyanukul [45].They utilize two decoding techniques namely "Solution Representation-1 (SR-1)" and "Solution Representation-2 (SR-2)" to solve CVRP where the first decoding technique (SR-1) is an improved version of their previous work [46].Mainly a local search is added to their previous solution representation which was a decoding technique based on relative position of vehicles.The second decoding technique (SR-2) is an enhanced version of SR-1 where they consider a vehicle coverage radius in addition to other two reference points for the vehicles in their representation.Even though Ai achieves better results [44] than their previous work [46], the computational time is still noticeable.A hybrid PSO is presented in [47] where they obtained initial solutions to run their PSO utilizing MPNS-GRASP [48].To reduce the computational time, they utilize expanding neighborhood search [49,50] as well.They utilize path relinking technique [51] to guide the PSO.There is another hybrid approach in [52] which hybridizes two swarm optimization techniques, namely ACO and PSO.ACO mainly make clusters of customer and builds the routes, and then PSO is utilized in short term memory to speed up the convergence by laying pheromone on the global and personal best solutions only.However, these algorithms have shown some good results at the expense of considerable computational time as all the approaches are mainly devoted to obtaining optimal solutions to the problems rather than concentrating on the main idea of approximation algorithms which is to obtain good solutions for large problems within reasonable computational time.
A candidate solution represented by a particle in PSO is usually in an encoded form, which needs a decoding procedure to obtain the comprehensive meaning of it.In this article, we represent a solution in a simple vector form to easily apply different operations of PSO, and at the same time we propose an effective decoding algorithm to interpret the solution.Decoded results are then exploited with a unique bilayer local search scheme.One layer of local search is for each particle in every iteration.The second layer of local search is applied on the pool of the best particles generated in different generations.The local search scheme improves the solution very effectively.Addition of these simple steps makes the algorithm very fast as compared to other available techniques and shows remarkable results in handling benchmark problems.
The rest of the article is organized as follows.A clear formulation of the problem statement is presented in Section 2. The proposed bilayer local search-enhanced PSO (BLS-PSO) is introduced describing the methodology in detail in Section 3. Section 4 presents the parameter optimization, the discussion of computational results and the comparison with other best known PSO-based algorithms.Finally, conclusion and future remarks are made in Section 5.

Mathematical Formulation
CVRP is usually stated as a graph of a set C = {c i |i = 1, 2, . . ., n} of n nodes, resembling customers, along with associated demands D = {d i |i = 1, 2, . . ., n} and required service times ST = {st i |i = 1, 2, . . ., n}.If no service time is specified for the customers in a problem definition, then it is set to 0. Among the nodes, a special node, say c h , is called the depot, from where the demands of other nodes are met using some vehicles.For the depot c h , the service time st h = 0.A vehicle starts from the depot, serves the demands of some nodes, and comes back to the depot.The objective is to serve all the customers using a set of vehicles, keeping the cost (usually calculated as the total distance travelled by the vehicles) as low as possible.A capacity constraint is associated with every vehicle, which specifies the maximum amount of total demands of the nodes that can be served by it.Maximum travel distance is another constraint that, if specified, puts a limit on the total distance traveled along with service time spent by a vehicle [43].
The required variables for the CVRP formulation are presented in Table 1.
fulfilling the following constraints Equations ( 2)-( 8) are the constraints to ensure the validity of a solution to Equation (1).Equations ( 2) and (3) ensure that a customer is served by exactly one vehicle.Equation (4) ensures the direction of a vehicle to serve a node i.e., if a vehicle travels from node c j to node c o then it never returns to the node c o from any other node c k or depot c h .Equation (5) restricts the customers to be assigned to a vehicle exceeding the maximum load capacity of that vehicle.Equation (6) confirms the fulfillment of the maximum travel distance constraint of a vehicle.Equation (7) forces a vehicle to start the journey from and end at the depot c h and never visit the depot any more.Equation (8) restricts sub-tour generated within a route.Equation ( 9) is the integrality constraint for the formulation.

Proposed Methodology
The traditional PSO algorithm applies the velocities on the particles to explore new solutions.We propose an enhanced PSO by incorporating a bilayer local search approach with it.The first layer of local search is applied on the whole population, i.e., on each particle in the population in every iteration.Thus, it improves the quality of the solutions obtained by traditional PSO by searching their neighborhoods.The second layer of local search is implemented on the pool of the best particles found in different iterations.This layer is to give additional chances to the best particles to look for better solutions in their neighborhoods.Such local searches allow the particles to have a better exploration around the current position as otherwise a particle may move from a position missing out a nearby potential solution.
We represent a particle in PSO in an encoded vector form to apply PSO operations easily on them.Hence, we also propose a decoding method to build a real comprehensive solution out of a particle.The overall workflow of proposed bilayer local search-enhanced particle swarm optimization (BLS-PSO) is sketched in Figure 2.
The proposed BLS-PSO technique to solve a capacitated vehicle routing problem (CVRP) can broadly be thought as having three major parts: the application of the particle swarm optimization (PSO), the decoding step to obtain the detail routes encoded in a solution and the bilayer local search approach to search for a better solution in the neighborhood.We discuss these steps in detail in the following subsections.
of local search is applied on the whole population, i.e., on each particle in the population in every iteration.Thus, it improves the quality of the solutions obtained by traditional PSO by searching their neighborhoods.The second layer of local search is implemented on the pool of the best particles found in different iterations.This layer is to give additional chances to the best particles to look for better solutions in their neighborhoods.Such local searches allow the particles to have a better exploration around the current position as otherwise a particle may move from a position missing out a nearby potential solution.
We represent a particle in PSO in an encoded vector form to apply PSO operations easily on them.Hence, we also propose a decoding method to build a real comprehensive solution out of a particle.The overall workflow of proposed bilayer local search-enhanced particle swarm optimization (BLS-PSO) is sketched in Figure 2.
The proposed BLS-PSO technique to solve a capacitated vehicle routing problem (CVRP) can broadly be thought as having three major parts: the application of the particle swarm optimization (PSO), the decoding step to obtain the detail routes encoded in a solution and the bilayer local search approach to search for a better solution in the neighborhood.We discuss these steps in detail in the following subsections.

PSO for CVRP
PSO is a population-based metaheuristic which consists of particles or chromosomes consisting of candidate solutions.In this work, a chromosome is a "Customer Vector", a permutation of all the customers specified in a problem statement.We can also think such a customer vector as a position vector in the solution space.
At first, we generate a predetermined number of random particles for PSO.We also assign a random velocity for each particle to guide the particles to search for better solutions utilizing the PSO scheme.If the position of the l-th particle is   −1 on iteration (t − 1), then the particle can move to its new position   () influenced by its velocity    , as done in the traditional PSO, using

PSO for CVRP
PSO is a population-based metaheuristic which consists of particles or chromosomes consisting of candidate solutions.In this work, a chromosome is a "Customer Vector", a permutation of all the customers specified in a problem statement.We can also think such a customer vector as a position vector in the solution space.
At first, we generate a predetermined number of random particles for PSO.We also assign a random velocity for each particle to guide the particles to search for better solutions utilizing the PSO scheme.If the position of the l-th particle is A t−1 l on iteration (t − 1), then the particle can move to its new position A (t) l influenced by its velocity V t l , as done in the traditional PSO, using The velocity also gets updated at each iteration, according to Shi and Eberhart [53], as where V t l is the velocity of the lth particle in the t-th iteration, B l is the best solution attained by the l-th particle till (t − 1) iterations, and G is best solution attained by any particle in the population in (t − 1) iterations.k 1 , k 2 and k 3 are parameters to guide the velocity to move a particle to a certain direction specifying the effect of its own velocity, attraction to the best position ever reached by itself and the best position ever reached by any particle in the population.Effect of the history of velocity is controlled by the parameter k 1 , influence of the personal best position is determined by k 2 and k 3 maintains the balance with the global best solution at any stage.These influences are basically the strength of PSO algorithm which leads a particle to a better position.Thus, we have utilized the proposal of Shi and Eberhart [53] instead of using the straightforward PSO velocity update equation from Eberhart and Kennedy [40] to try to keep the preeminent balance among the influential portions in Equation (11).
PSO was first proposed for continuous problem domains where ⊕ can be considered as vector addition and → as vector subtraction.However, CVRP is defined in discrete domain.Hence, to apply PSO in discrete domain, we adopt the swap sequence-based PSO [54].We need to map the meaning of Equations ( 10) and (11) so that PSO can deal a discrete type problem.In this purpose [54], authors proposed a new definition to the characteristics of two symbols, ⊕ and →.X ⊕ SS means applying a sequence of swap sequence on the elements of X to reach Y, a new position in the solution space, while X → W gives a basic swap sequence BSS, the applying of which W can be converted to X.A detailed description is presented here to understand the mechanism of the operators.
Swap sequence is basically moving a particle from one solution or position to another solution or position in the solution space.Suppose an encoded solution or particle in the solution space is X = (1, 3, 5, 2, 4) assuming the problem has five nodes in the problem statement.Then SS can be a swap sequence consisting of one or more swap operators that move X to its new position Y.For example SS = SO(2,3) may be a random initial swap sequence having a single swap operator or velocity for particle X.Thus, X ⊕ SS = (1, 3, 5, 2, 4) ⊕ SO(2, 3) = (1, 5, 3, 2, 4) would be Y the new particle position.However, SS may be a swap sequence having more than one swap sequence or basic swap sequence or swap operators as well.Let another SS = SO(2, 3) ⊕ SO(4, 5) having two swap operators.Applying a swap sequence with more than one swap operators on a particle is simply applying swap operators one after another on the particle.Thus, X ⊕ SS = (1, 3, 5, 2, 4) ⊕ SO(2, 3) ⊕ SO(4, 5) = (1, 5, 3, 2, 4) ⊕ SO(4, 5) = (1,5,3,4,2).This definition helps to make a swap sequence or velocity from Equation (11) and then swap sequence will apply on a particle according to Equation (10) to move the particle from one position to other.

Decoding the Customer Vector
The customer vectors in the particles presented in the previous section are in an encoded form.It is not only difficult to comprehend them but also not suitable to apply the proposed local search operations on them.We propose a novel distance-based decoding method presented in Algorithm 1 to make a complete CVRP solution from a chromosome or customer vector.
In the proposed decoding method, the required number of vehicles m is taken as input parameter from the problem statement.We first assign the first m customers in the customer vector to the m vehicles (line 2 to line 4 in Algorithm 1).Afterwards we assign each of the remaining (n − m) customers to one of the m routes according to Algorithm 1.For a customer c i , where m < i < n, we first try to add it to the end of a suitable route (line 9 to line 16 in Algorithm 1).For this, we rank the routes according the distances of their last assigned nodes from c i , and try to assign it to the end of the route having c i 's closest possible customer at its end.According to the rank, we keep on trying the routes one after another to find a suitable route where c i can be assigned without exceeding the corresponding vehicle's capacity.
If c i cannot be assigned to any route following this strategy, i.e., its addition to any route exceeds the limit of the vehicle, we try to replace a customer having a lower demand than c i .This procedure is depicted from line 17 to 31 in Algorithm 1.An example scenario presents the major steps of the mechanism of Algorithm 1 followed by the pseudocode.route(r) = x(r); // put the first m nodes in x as the starting nodes of the m routes 3. end for 4. i = m; 5. while i <= |x| do 6.
i = i + 1; // index of the next node in x to insert in any route 7.
sortedRoutes = indices of routes sorted in non-increasing order based on distances of their last nodes to x(i) 8. assigned = 0; 9.
if addition of x(i) to the end of route(r) satisfies Equations ( 5) and ( 6) then 12.
end for 17.
if assigned == 0 // x(i) was not assigned to any route in the loop from line 10 to 17 then 18.
maxRouteSize = maximum number of nodes assigned to a route 19.
for k = maxRouteSize down to 1 do 20.sortedRoutes = indices of routes sorted in non-increasing order based on the increases of the distances if the kth nodes of the routes are replaced by x(i) 21.
add kth node to the end of x; // to be assigned again 25.
replace kth node of route(r) by x(i); Note that a random customer vector (2, 8, 6, 9, 11, 7, 3, 4, 5, 10) is found in an iteration of PSO.Considering the vehicle capacity constraint to be 40, the number of vehicles required is assumed to be given as 4.
Now the first 4 customer nodes from the customer vector will be assigned to the 4 routes, say, R1 = (2), R2 = (8), R3 = (6) and R4 = (9).Afterwards the algorithm will try to insert node 11 to any of the routes.For this, the routes are sorted in non-increasing order based on distances of node 11 from their last nodes, e.g., 2, 8, 6 and 9. Suppose the sorted sequence is (R4, R2, R1, R3).Since adding node 11 to route 4 does not violate the capacity constraint, R4 becomes (9,11).The next node in the customer vector is node 7, for which suppose that sorted route sequence is (R1, R2, R3, R4).R1 can successfully accommodate node 7, and it becomes (2,7).Following this manner, suppose that nodes 3 and 4 are assigned to R2 and R3, making them (8, 3) and (6, 4), respectively.Now based on node 5, suppose that the sorted route sequence is (R2, R1, R3, R4).However, note that adding node 5 exceeds the vehicle capacity for each of the routes.So, it cannot be assigned in this step.Now according to the second strategy, the algorithm tries to the replace a node with 5 from any of the routes.For this the routes are again sorted in non-increasing order based on the increases of the distances if the nodes 7, 3, 4 and 11 are replaced by node 5. Suppose that the sorted sequence is (R4, R1, R3, R2).So, the algorithm tries R4 first, and finds that replacing node 11 with node 5 does not violate the capacity constraint.Hence, R4 becomes (9,5), and node 11 is added to the end of the customer vector, making a temporary customer vector as (2,8,6,9,11,7,3,4,5,10,11).The node 10 is then assigned to R2, making it (8,3,10), according to the first strategy.Now, the algorithm tries to insert node 11 again to any of the routes.According the first strategy, suppose the sorted route sequence is (R1, R4, R3, R2).Hence, node 11 can be accommodated in R1, making it (2,7,11).After decoding steps, now we have route vectors at hand (R1 (2,7,11), R2 (8, 3, 10), R3 (6, 4), R4 (9, 5)) to apply the bilayer local search operations to find if any better solution is available in the neighborhood of the particle.

Proposed Bilayer Local Search
A local search can exploit the neighborhood of a particle with the expectation of getting a better solution.We propose two layers of local searches to be applied during the iterations of the PSO algorithm.Here we explain them in detail.

The First Layer Local Search
The first layer local search is applied on the all "feasible" particles, for which routes could be successfully made by Algorithm 1, in the population in every iteration of PSO.We apply Intra-Route Local Search (Intra-RLS) followed by Inter-Route Local Search (Inter-RLS) to exploit the neighborhood of a particle to look for any better solution before moving to another position in the solution space in the next iteration of PSO.In the Intra-RLS, we apply "swap" and "insertion" strategies on every route constructed for the particle at hand.In the "swap" stage, we interchange the neighboring nodes in a route if it minimizes the total distance of the route.For "insertion", we pick a node, and try to insert it in other possible positions in the route if it decreases the total distance of the route.After all the nodes have been checked, the move with the best gain can be executed.We apply "insertion" operation using the nodes one after another onward from the starting of a route until any node succeeds to move to another position.We do not perform all possible exploration here to keep the complexity of the algorithm tolerable.Moreover, there remains a possibility for the proposed second layer of local search to explore more "insertion" options on the potential solutions that are selected in the "pool" (please see Section 3.3.2for details).
A portrayal of the attempts made by the "insert" and "swap" operations on a route (c d , c e , c f and c g ) are illustrated in Figure 3.The first layer local search is applied on the all "feasible" particles, for which routes could be successfully made by Algorithm 1, in the population in every iteration of PSO.We apply Intra-Route Local Search (Intra-RLS) followed by Inter-Route Local Search (Inter-RLS) to exploit the neighborhood of a particle to look for any better solution before moving to another position in the solution space in the next iteration of PSO.In the Intra-RLS, we apply "swap" and "insertion" strategies on every route constructed for the particle at hand.In the "swap" stage, we interchange the neighboring nodes in a route if it minimizes the total distance of the route.For "insertion", we pick a node, and try to insert it in other possible positions in the route if it decreases the total distance of the route.After all the nodes have been checked, the move with the best gain can be executed.We apply "insertion" operation using the nodes one after another onward from the starting of a route until any node succeeds to move to another position.We do not perform all possible exploration here to keep the complexity of the algorithm tolerable.Moreover, there remains a possibility for the proposed second layer of local search to explore more "insertion" options on the potential solutions that are selected in the "pool" (please see Section 3.3.2for details).
A portrayal of the attempts made by the "insert" and "swap" operations on a route (c d , c e , c f and c g ) are illustrated in Figure 3.The Inter-RLS strategy also applies the similar "insertion" and "swap" operations, but between two different routes.Inter-RLS is portrayed in Figure 4.The arrows show the "insert" options to be tried by node 1, from which the best one is to be performed if it decreases the total distance of the two concerned routes.It is worth mentioning here that the "insertion" operation is not allowed to remove all the nodes of a route into another route thus the required number of vehicles will be constant throughout the program run.

The Second Layer Local Search
The operations of the second layer of local search are similar to those in the first layer.As mentioned earlier, this layer of local search is applied on a pool of the global best solutions attained The Inter-RLS strategy also applies the similar "insertion" and "swap" operations, but between two different routes.Inter-RLS is portrayed in Figure 4.The arrows show the "insert" options to be tried by node 1, from which the best one is to be performed if it decreases the total distance of the two concerned routes.

The First Layer Local Search
The first layer local search is applied on the all "feasible" particles, for which routes could be successfully made by Algorithm 1, in the population in every iteration of PSO.We apply Intra-Route Local Search (Intra-RLS) followed by Inter-Route Local Search (Inter-RLS) to exploit the neighborhood of a particle to look for any better solution before moving to another position in the solution space in the next iteration of PSO.In the Intra-RLS, we apply "swap" and "insertion" strategies on every route constructed for the particle at hand.In the "swap" stage, we interchange the neighboring nodes in a route if it minimizes the total distance of the route.For "insertion", we pick a node, and try to insert it in other possible positions in the route if it decreases the total distance of the route.After all the nodes have been checked, the move with the best gain can be executed.We apply "insertion" operation using the nodes one after another onward from the starting of a route until any node succeeds to move to another position.We do not perform all possible exploration here to keep the complexity of the algorithm tolerable.Moreover, there remains a possibility for the proposed second layer of local search to explore more "insertion" options on the potential solutions that are selected in the "pool" (please see Section 3.3.2for details).
A portrayal of the attempts made by the "insert" and "swap" operations on a route (c d , c e , c f and c g ) are illustrated in Figure 3.The Inter-RLS strategy also applies the similar "insertion" and "swap" operations, but between two different routes.Inter-RLS is portrayed in Figure 4.The arrows show the "insert" options to be tried by node 1, from which the best one is to be performed if it decreases the total distance of the two concerned routes.It is worth mentioning here that the "insertion" operation is not allowed to remove all the nodes of a route into another route thus the required number of vehicles will be constant throughout the program run.

The Second Layer Local Search
The operations of the second layer of local search are similar to those in the first layer.As mentioned earlier, this layer of local search is applied on a pool of the global best solutions attained It is worth mentioning here that the "insertion" operation is not allowed to remove all the nodes of a route into another route thus the required number of vehicles will be constant throughout the program run.

The Second Layer Local Search
The operations of the second layer of local search are similar to those in the first layer.As mentioned earlier, this layer of local search is applied on a pool of the global best solutions attained in different iterations the PSO algorithm.This is an approach to exploit the best-found solution more deeply as some of them might be very close to the optimal solution.
To limit the growth of the pool, we refine it by deleting some of the solutions in the pool according to Algorithm 2 at every iteration.It keeps the best solution in the pool.It also keeps the solutions that were updated by the local search in the previous iteration for the hope that they may be further updated in the upcoming iterations to attain better solutions.The solution that neither is the best in the pool nor was updated in the previous iteration is deleted from the pool because it can no longer be updated by the proposed local search.

Computational Result Analysis
We have implemented the proposed bilayer local search enhanced particle swarm optimization (BLS-PSO) in MATLAB R2017a environment on a personal computer with a processor Intel ® Pentium ® Dual-Core CPU T4300 @ 2.10 GHz (only single core is used), memory (RAM) of 2.0 GB, and Windows 10 Education version.
We have considered seven sets of capacitated vehicle routing benchmark problems of the well-known CVRPLIB problem set which are also used for the verification of many PSO-based proposals [41], can be downloaded from http://vrp.atd-lab.inf.puc-rio.br/index.php/en/.The benchmark problem sets are called Set A, Set B, Set P [55], Set E, Set F, Set M [56] and Set CMT [57].The problems in the Sets A, B, E, F, M and P specify only capacity constraints for the vehicles, but not any constraint on the distance travelled by a vehicle.Total number of customers varies from 15 to 199 in these sets, and locations of the customers, except Sets B and M, are generated in random scattered manner or, following semi-clustered fashion [52].The problems in Set B and M are generated in a clustered layout [52].Though some researchers [41] pick a subset of the 92 problems to demonstrate the performances of their proposed algorithms, we have considered all the problems to show the real portrayal of the performance of our proposed BLS-PSO on a variety of problems.
Chen et al. [41] considered all the 14 problems of the CMT problem set.CMT1, CMT2, CMT3, CMT4, CMT5, CMT11 and CMT12 consider only the capacity constraint for the vehicles whereas CMT6, CMT7, CMT8, CMT9, CMT10, CMT13 and CMT14 have additional restrictions on the maximum travelling distance of a vehicle.In CMT, problem sets include 50 to 199 customers, whose locations are randomly scattered in CMT1 to CMT10 problem sets while clustered locations are set in CMT11 to CMT14.

Setting the Parameters
A fundamental characteristic of metaheuristic algorithms is that their performance is greatly parameter value intensive.Quality of solution obtained and required computational time to solve a problem critically depend on the parameter sets.As optimizing the parameter sets for any metaheuristic algorithm is itself an optimization problem, we have run all the problems for different parameter sets and selected the best result achieving parameter set for this article.In the parameter set optimization step, there are several parameter values to be decided such as population size or the number of particles PS Particle , maximum number of iterations, T, constants in PSO equation k 1 , k 2 and k 3 , termination criteria and so on.We have performed a set of experiments to fix the parameter sets to obtain the best out of the proposed bilayer local search-based particle swarm optimization (BLS-PSO) algorithm.The experimental designed setup for the experiments are portrayed in Table 3.Then, k 1 value is predicted via Equation ( 12) as used in [58].
The termination criteria used in our approach is that the program run until the maximum number of iterations, T. For BLS-PSO, we have reported the average of the performances of 20 runs for the best combination of k 2 and k 3 for a problem, unless specified otherwise.
The experiments show that performance of the proposed algorithm get improved with increasing both population size, PS Particle and maximum number of iteration, T portrayed in Figure 5.It is observed that even if the PS Particle and T is increased up to certain limit, the proposed algorithm can find results equal to best known solutions at most for 85 benchmark instances.When the setting PS Particle is equal to the number of nodes in respective problem statement and T is 50 or more it can find the results equal to the best-known solutions for maximum 85 benchmark problems.Not finding the outputs equal to the best-known solutions for the remaining problem instances are likely due to having two reasons either the small population size or the number of maximum iteration which is set as low as 100.However, as an approximate algorithm is expected to find a near optimal or, a "good" solution within a small computational time, we have chosen the PS Particle as equal to n and T as 50.
For setting the value of k 2 and k 3 , we have run the program again with the parameter setting of PS Particle and T as mentioned above.In these experimental setups, for a particular k 2 we have run the program for all k 3 values.For any setup, the number of best known results found is marked on the pie chart in Figure 6.Close observation portrays that there is no significant converging point to fix the value of k 2 and k 3 to obtain the best performance of the proposed algorithm.However, the computational time required for one complete run of the proposed algorithm for an instance is not significant.Thus, experiments for all the combinations of k 2 and k 3 for each instance is performed thoroughly and the best result achieved by the proposed BLS-PSO is reported in later subsections.
PSParticle is equal to the number of nodes in respective problem statement and T is 50 or more it can find the results equal to the best-known solutions for maximum 85 benchmark problems.Not finding the outputs equal to the best-known solutions for the remaining problem instances are likely due to having two reasons either the small population size or the number of maximum iteration which is set as low as 100.However, as an approximate algorithm is expected to find a near optimal or, a "good" solution within a small computational time, we have chosen the PSParticle as equal to n and T as 50.For setting the value of k2 and k3, we have run the program again with the parameter setting of PSParticle and T as mentioned above.In these experimental setups, for a particular k2 we have run the program for all k3 values.For any setup, the number of best known results found is marked on the pie chart in Figure 6.Close observation portrays that there is no significant converging point to fix the value of k2 and k3 to obtain the best performance of the proposed algorithm.However, the computational time required for one complete run of the proposed algorithm for an instance is not significant.Thus, experiments for all the combinations of k2 and k3 for each instance is performed thoroughly and the best result achieved by the proposed BLS-PSO is reported in later subsections.

Result Analysis
Tables 4-10 present the computational results obtained by BLS-PSO for several benchmark problem sets, where m represents the required number of vehicles for the corresponding problem instance, Cost reports the best fitness obtained by an algorithm along with the corresponding gap in percentage, representing the quality of a solution, calculated as where CostBLS-PSO is the best solution obtained by the proposed method and BKS is the best-known solution for the same instance reported in the literature.AvgT is average of CPU times required by

Result Analysis
Tables 4-10 present the computational results obtained by BLS-PSO for several benchmark problem sets, where m represents the required number of vehicles for the corresponding problem The proposed BLS-PSO algorithm has reached the best-known solutions for 84 instances out of the 106 Euclidean distance-based benchmark instances used in our experiments.BLS-PSO has an average gap of 0.07%, 0.14%, 0.31%, 0.26%, 3.45%, 0.2% and 0.64% to solve instances in Sets A, B, E, F, M, P and CMT, respectively.It requires an average CPU Time of 5.5 s, 6.13 s, 11.51 s, 25.59 s, 68.69 s, 4.72 s and 13.33 s to solve instances in Set A, B, E, F, M, P and CMT, respectively.Thus, BLS-PSO has been found to obtain the results equal to the best-known results for most of the existing benchmark problems within short computational time.
BLS-PSO has not met the best-known results for some instances.The main reason behind this may be the small population size or the lower value of the maximum number of iterations.Increasing them are expected to attain better solutions.

Comparative Performance of PSO-Based Approaches
To portray the comparison of the performance of our proposed algorithm with other PSO-based algorithms in the literature, we have taken the same instances used by other researchers [41,44,47,52].Table 11 contains the results for 16 instances from the Sets A, B, E, F, M and P, and Table 12 lists results for all the 14 instances of Set CMT.Here, PSO-based algorithms DPSO [41], SR1, SR2 [44,46] and PSCO [52] are compared with proposed BLS-PSO.Here, Cost is the result obtained by respective algorithms and a value is in boldface font when it is equal to the best-known result reported in the literature.The superscript "a" indicates that the cost is inferior to what is achieved by the proposed BLS-PSO.BT refers to the computational time in seconds to reach the corresponding result.For the proposed algorithm, we have reported the AvgT as mentioned earlier.As mentioned above there are no convergence of solution for specific k 2 and k 3 value, we have mentioned AvgCost as the average of all the solutions obtained by BLS-PSO for different k 2 and k 3 .The deviation of the AvgCost with the best-known solutions is mentioned as SD in Table 11.It is evident that even the average time taken by the proposed method is far better, let alone the shortest time, than the shortest time of other approaches for each of the problem instances under consideration.A-n33-k5 indicates the problem statement is from Set A, having total 33 nodes and 5 vehicles.a indicates the results compared to those BLS-PSO has obtained better results.
The overall results are also graphically represented in Figure 7 to show the comparative performances more comprehensively.The negative gap in Figure 7a corresponds to the problem instance, for which the BLS-PSO has obtained a new optimal result.Figure 7b demonstrates that time consumed by BLS-PSO is significantly less than the other approaches, where the plot is almost touching the horizontal axis for all the instances.The overall results are also graphically represented in Figure 7 to show the comparative performances more comprehensively.The negative gap in Figure 7a corresponds to the problem instance, for which the BLS-PSO has obtained a new optimal result.Figure 7b demonstrates that time consumed by BLS-PSO is significantly less than the other approaches, where the plot is almost touching the horizontal axis for all the instances.Table 12 portrays the comparison of performance in handling CMT problem set instances.Here we have compared SR1, SR2 [44,46], HybPSO [47] and PSCO [52] with the proposed BLS-PSO.Table 12 portrays the comparison of performance in handling CMT problem set instances.Here we have compared SR1, SR2 [44,46], HybPSO [47] and PSCO [52] with the proposed BLS-PSO.Inf indicates when an algorithm could not achieve any feasible solution.a indicates the results compared to those BLS-PSO has obtained better results.
Figure 8 shows clearer comparative graphical representations of the performances of the PSO-based algorithms dealing instances in Set CMT.Here, BLS-PSO has obtained results equal to the best known results with an average gap of 0.64% whereas SR1, SR2, PACO and HybPSO have gap of 1.788%, 0.874%, 0.913% and 0.084% respectively.However, the average of the reported CPU times required by BLS-PSO is only 13.33 s to solve the instances in Set CMT whereas the average of the best times required by SR1, SR2, PACO and HybPSO are 84.1 s, 163.1 s, 302.58 s and 48.13 s, respectively.Hence, although the gap values obtained by the proposed algorithm are higher than the other methods for some instances, BLS-PSO is found to reach to the best known results for most of the cases in significantly lower computational times as compared to the other algorithms.This advocates for its efficiency.Inf indicates when an algorithm could not achieve any feasible solution.a indicates the results compared to those BLS-PSO has obtained better results.
Figure 8 shows clearer comparative graphical representations of the performances of the PSObased algorithms dealing instances in Set CMT.Here, BLS-PSO has obtained results equal to the best known results with an average gap of 0.64% whereas SR1, SR2, PACO and HybPSO have gap of 1.788%, 0.874%, 0.913% and 0.084% respectively.However, the average of the reported CPU times required by BLS-PSO is only 13.33 s to solve the instances in Set CMT whereas the average of the best times required by SR1, SR2, PACO and HybPSO are 84.1 s, 163.1 s, 302.58 s and 48.13 s, respectively.Hence, although the gap values obtained by the proposed algorithm are higher than the other methods for some instances, BLS-PSO is found to reach to the best known results for most of the cases in significantly lower computational times as compared to the other algorithms.This advocates for its efficiency.

Conclusions
In this article, we have proposed a bilayer local search-based Particle Swarm Optimization (BLS-PSO) incorporated with a simple chromosome structure to apply PSO operations easily.We have also proposed a decoding technique that depends on the number of vehicles needed to serve the customers which is specified in the problem statement of benchmark instances considered in this article, and then build the routes in a comprehensible solution on which we apply the proposed local

Conclusions
In this article, we have proposed a bilayer local search-based Particle Swarm Optimization (BLS-PSO) incorporated with a simple chromosome structure to apply PSO operations easily.We have also proposed a decoding technique that depends on the number of vehicles needed to serve the customers which is specified in the problem statement of benchmark instances considered in this article, and then build the routes in a comprehensible solution on which we apply the proposed local searches.The first layer of the proposed bilayer local search approach searches the neighborhood of all particles to find better solution whereas the second layer is applied on the pool of the best particle of different generations.The first layer is important for exploring the neighborhoods of the particles before modifying the chromosomes by PSO operations, and the second layer gives more chances to the best selected particles to exploit more of their neighborhoods.Such a contribution helps the search strategy to obtain optimal or near optimal solutions with the expense of very short computational time, which is also evident from the experimental results obtained by applying the method on different benchmark problem sets.
The proposed decoding method can also be easily integrated with other population-based metaheuristics to solve CVRP where a solution is encoded as a permutation of customer nodes.The bilayer search strategy may also enhance the performances of other nature-inspired searching strategies.Moreover, introduction of a probabilistic approach to the decoding algorithm may also contribute to generate different potential solutions.The estimation procedure of the required number of vehicles may also be enhanced by using a stochastic technique other than taking as input from the problem statement to solve generalized vehicle routing problems.We leave all these possibilities to be investigated in our future endeavors.

Figure 1 .
Figure 1.A portrayal of capacitated vehicle routing problem.

Figure 1 .
Figure 1.A portrayal of capacitated vehicle routing problem.

Algorithm 1 .
Distance-based decoding method.Inputsx: customer vector m: number of vehicles given in the problem statement Outputs routesMade: a Boolean value Indicating whether routes could successfully been made route: The routes made Procedure 1. for r = 1 to m do 2.

Figure 4 .
Figure 4.An illustration of Inter-Route Local Search.

Figure 4 .
Figure 4.An illustration of Inter-Route Local Search.

Figure 4 .
Figure 4.An illustration of Inter-Route Local Search.

Figure 5 .
Figure 5. Variation in convergence of BLS-PSO depending on different PSParticle & T.Figure 5. Variation in convergence of BLS-PSO depending on different PS Particle & T.

Figure 5 .
Figure 5. Variation in convergence of BLS-PSO depending on different PSParticle & T.Figure 5. Variation in convergence of BLS-PSO depending on different PS Particle & T.

Figure 6 .
Figure 6.Scattered variation of convergence of the proposed BLS-PSO depending on different values of parameters k2 and k3.

Figure 6 .
Figure 6.Scattered variation of convergence of the proposed BLS-PSO depending on different values of parameters k 2 and k 3 .

Table 1 .
Related notations and variables for CVRP formulation defined in this article.

Table 2 .
Suppose a CVRP example problem scenario with 11 nodes as given in Table2.Hub ID, Customer IDs and associated demands in example problem instance.

Table 3 .
Different values of parameters used in experimental setup.

Table 10 .
Performance of BLS-PSO in solving the benchmark problems in Set CMT.

Table 11 .
Performance comparison of BLS-PSO with respect to other PSO-based approaches for solving some selected instances.

Table 12 .
Performance comparison of BLS-PSO with respect to other PSO-based approaches for solving the instances of Set CMT.

Table 12 .
Performance comparison of BLS-PSO with respect to other PSO-based approaches for solving the instances of Set CMT.