A Bee Evolutionary Algorithm for Multiobjective Vehicle Routing Problem with Simultaneous Pickup and Delivery

A new closed-loop supply chain logistics network of vehicle routing problem with simultaneous pickups and deliveries (VRPSPD) dominated by remanufacturer is constructed, in which the customers are originally divided into three types: distributors, recyclers, and suppliers. Furthermore, the fuel consumption is originally added to the optimization objectives of the proposed VRPSPD. In addition, a bee evolutionary algorithm guiding nondominated sorting genetic algorithm II (BEG-NSGA-II) with a two-stage optimizationmechanism is originally designed to solve the proposedVRPSPDmodel with three optimization objectives: minimum fuel consumption, minimum waiting time, and the shortest delivery distance. The proposed BEG-NSGA-II could conquer the disadvantages of traditional nondominated sorting genetic algorithm II (NSGA-II) and algorithms with a two-stage optimization mechanism. Finally, the validity and feasibility of the proposed model and algorithm are verified by simulating an engineering machinery remanufacturing company’s reverse logistics and another three test examples.


Introduction
Recycling and remanufacturing is an integral system by which the old or discarded products are recycled and then processed, recovered, and sold as like-new ones in the markets [1], and the environmental impacts are also reduced by maintaining the geometrical form of the products, thus preserving its environmental and economic values [2,3].It is one of the important ways to realize the green sustainable development [4,5].For example, in recent years, China offers subsidies on a range of vehicles such as small-and mediumsized old vehicles, the rural bus, and yellow-sticker vehicles subsidy in its "automotive replacement" policy, to encourage vehicle owners to submit vehicles to officially recognized endof-life vehicle (ELV) dismantlers [6].Reverse logistics is a central problem of recycling and remanufacturing of waste products [7][8][9].Its path planning has a great influence on distribution efficiency and environmental protection efficiency, which usually includes three main models: vehicle routing problem that delivers the goods before pickups, vehicle routing problem that delivers and pickups the goods in a mixed way, and vehicle routing problem with simultaneous pickups and deliveries (VRPSPD) of the goods, etc. [10].The model of VRPSPD is better than other ones in environmental protection and saving cost, for it can make full use of the vehicle load space [11].Therefore, the VRPSPD problem attracts more and more researchers in recent years.
The VRPSPD problem was first introduced by Min [12], who developed a solution procedure to solve a public library distribution with two vehicles.Based on a set of covering formulations, Klose et al. [13] implemented a Branch and Price approach for the VRPSPD.Li and Lim [14] proposed a tabu-embedded simulated annealing algorithm which restarted a search procedure from the current best solution after several nonimproving search iterations.Crispim and Brandão [15] proposed a hybrid heuristic algorithm which combined the tabu search with the falling variable field method for solving the vehicle routing problem with backhauls.Ropke and Pisinger [16] designed an adaptive large neighborhood search heuristic composed of a number of competing subheuristics, which improved many of the best known solutions.Montané and Galvao [17] used tabu search algorithm and local optimization method for the VRPSPD problem with maximum travel constraint.Based on GLNPSO, a PSO algorithm with multiple social structures, Ai and Kachitvichyanukul [18] developed a new PSO algorithm for VRPSPD.C ¸atay [19] proposed an ant colony algorithm by employing a new saving-based visibility function and pheromone updating procedure.Mingyong and Erbao [20] solved the VRPSPD problem with time windows and travel constraints by improving the differential evolution algorithm.Wu et al. [21] proposed a hybrid chaotic quantum evolutionary algorithm which successfully reduced the number of vehicles and the distance of the distribution according to the data in the literature [15].Subramanian et al. [22] proposed a branch-cut-and-price approach for the VRPSPD problem.By combining a local search named variable neighborhood descent algorithm into PSO, Goksal et al. [23] presented a heuristic approach for VRPSPD and improved several best known solutions.By analysing the vehicle load fluctuation characteristics, Zhang et al. [24] designed a heuristic factor and solved the VRPSPD problem with vehicle travel constraints based on an improved ant colony algorithm.Then Zhang et al. [25] studied the VRPSPD with time dependent vehicle routing problems and developed a hybrid algorithm that integrated both ant colony system and tabu search algorithms for solving it.Polat et al. [26] proposed a perturbation based variable neighborhood search heuristic for solving the VRPSPD problem with time limit.Avci and Topaloglu [27] developed a hybrid local search algorithm in which a nonmonotone threshold adjusting strategy is integrated with tabu search and effectively solved the VRPSPD problem with different types of vehicles.With the deterioration of the environment and the enhancement of environmental protection consciousness, the scholars have begun to take the vehicle fuel consumption as one of the optimization objectives in the past two years.Jin and Pei-Hua [28] used multistarting tabu search algorithm to study the heterogeneous vehicle low-carbon routing problem based on energy consumption and carbon emissions, then proposed a variable domain random search method, and solved the low energy vehicle routing problem about time windows [29].By analysing the solution space and complexity of fuelconsumption-minimizing capacitated vehicle routing problem (FCM-CVRP) and capacitated vehicle routing problem (CVRP), Wu et al. [30] found that the FCM-CVRP was more difficult to solve, and they proposed a two-phase algorithm to solve it.Considering the factors of road grade, Rao et al. [31] proposed a low-carbon vehicle routing problem model which has capacity constraints and objective of minimum fuel consumption.Ombuki et al. [32] used a multiobjective evolutionary genetic algorithm based on the Pareto sort for the Vehicle Routing Problem with Time Windows (VRPTW), and its two optimization objectives are minimum number of vehicles and shortest delivery distance.Lau et al. [33] took the shortest delivery distance and the shortest travel time as the optimization objectives and used NSGA-II algorithm guided by fuzzy logic for the routing problem with multiple warehouses, multiple customers, and a variety of products.Baños et al. [34] put forward a model with two objectives of shortest distance and minimum unbalance load, which considered the workload of the vehicle, and proposed a hybrid heuristic algorithm composed of genetic algorithm and simulated annealing algorithm for it.Hsueh [35] studied the vehicle routing problem with the objectives of minimizing the sum of the fixed costs and the expected fuel consumption costs.Jemai, Zekri, and Mellouli [36] studied the Multiobjective GA for green vehicle routing problem with different metrics: computation time, traveled distance, emissions volume, generational distance, spacing, entropy, and contribution.Soleiman Chaharlang and Ghaderi [37] proposed a multiobjective nonlinear programming model for the green vehicle routing problem (GVRP), including original and remanufactured products distribution (both delivery and pickup) of end-of-life (EOL) products.
Although there is a stream of research on the VRPSPD problem and related algorithms, we find that there are still some important problems which have been overlooked and needed to be solved.It is indicated that the traditional VRPSPD model only considers a single type of customer.However, in a real recycling and remanufacturing system, customers consist of various types, such as distributors, recyclers, and suppliers.In addition, although the optimization objective of fuel consumption has been studied in some models such as FCM-CVRP and CVRP, it has not been considered in the existing literatures with the traditional VRPSPD.However, the fuel consumption is actually a key factor affecting the resource consumption and environmental destruction in VRPSPD.With regard to the overlooked problems mentioned above, in this paper the customers are originally divided into three types which are distributors, recyclers, and suppliers, and a new closed-loop supply chain logistics network of VRPSPD dominated by remanufacturer is constructed, which is more consistent with the actual situation of recycling and remanufacturing.Furthermore, for the purpose of resource saving and environment protection, the fuel consumption is originally added to the optimization objectives of VRPSPD.Besides, inspired by the literature [38], a bee evolutionary algorithm guiding nondominated sorting genetic algorithm II (BEG-NSGA-II) with a two-stage optimization mechanism is originally designed to solve the proposed VRPSPD model with three optimization objectives: minimum fuel consumption, minimum waiting time, and the shortest delivery distance.The optimal solution should meet the travel constraints as well as load constraints.The proposed BEG-NSGA-II could conquer the disadvantages (premature convergence to local solution and low searching efficiency) of the traditional NSGA-II and the disadvantage (inability to gain stable and high quality initial population in the first stage) of the existing algorithms with two-stage optimization mechanism when it is used to solve the problems with rigor constraints.Finally, the validity and feasibility of the proposed model and algorithm are verified by simulating an engineering machinery remanufacturing company's reverse logistics and another three test examples.
The rest of this paper is organized as follows.The proposed VRPSPD model is constructed in the next section.In Section 2, the multiobjective optimization model of the proposed VRPSPD is proposed.In Section 3, the implementation details of the proposed BEG-NSGA-II are presented.Afterwards, experimental studies and discussions are made in Section 4. Finally, the conclusions are described in Section 5.

The Multiobjective Optimization Model of the Proposed VRPSPD
2.1.The Proposed VRPSPD Model.The proposed VRPSPD model in the paper is shown in Figure 1.This model is described as follows: a group of vehicles start from the recycling and remanufacturing plant with new products; then the vehicles distribute the new products to each distributor, take waste products back from each recycling business, or take new components from supply business during the set time.Some requirements must be met: the load of vehicles with waste products backed from the service points must meet the needs of all distributor points; the load of every service point cannot exceed the allowable values; each service point is serviced only once in a logistics delivery cycle.
The following assumptions are made for the network structure: (i) Due to the uncertainty of quantity and quality of recycled products, each service point (except waste treatment plan) in the model should establish an information-sharing platform at the beginning of each logistics cycle; then, according to the information, recycling and remanufacturing factory arranges distribution vehicles.
(ii) Distributor points are sale points of remanufacturing products as well as recycling points of waste products.
(iii) Recycling process center and recycling and remanufacturing plant are located in the same location.

Fundamental Assumptions and Parameter Description.
The VRPSPD model proposed in this paper is based on the following assumptions: (i) Every vehicle moves at a constant speed between any two customers; namely, the acceleration speed is zero.
(ii) Other energy consumption of distribution vehicles (such as air conditioning) is zero, which means Pacc = 0 in (2).
(iii) Every vehicle departs from the recycling company and returns to it in the end.
(iv) All vehicles' maximum payload is identical and they have the same weight as well.
(v) The fuel consumption for each vehicle is only affected by three factors, namely, travel distance, payload, and travel speed.
(vi) The road slope () between any two customers is set to zero, i.e., without considering the road slope of any two customers.
By analysing the operating process of recycling and remanufacturing reverse logistics network, recycling points and suppliers can be regarded as special distribution points.So the vehicle routing problem of the proposed model can be changed into VRPSPD.V = {v 0 , v 1 , v 2 , . .., v n , r 1 , r 2 , . ..,r m , p 1 ,p 2 , . .., p s } is used to express the set of all service points, in which v 0 is the recycling and remanufacturing plant, {v 1 , v 2 , ... v n }is the set of distributors, {r 1 , r 2 , ..., r m } is the set of recycling points, and {p 1 , p 2 , ..., p s } is the set of suppliers.Use the corresponding sequence number of each element of V to form a new set   ; i.e.,   = {0, 1, . .., n, n +1,. ..n + m,. .., n + m + s}.
The parameters used in the definition of multiobjective VRPSPD are as follows:  If the real delivery time for point i is later than L i , then there exists waiting time which can be formulated as follows: min Subject to: Equations ( 5), (6), and ( 7) indicate the three optimization objectives which are minimum waiting time, fuel consumption, and delivery distance.Equation (8) ensures that each service point must be served once and only once.Equation (9) ensures that each service point must be served by one and only one vehicle.
Equation (10) ensures that each vehicle can only be used at most once.Equation (11) ensures that the vehicle's load starting from V 0 is equal to the demand of all the customers that it serves.Equation ( 12) is the vehicle's load change equation that indicates the load after vehicle k severed i before accessing j.Equation ( 13) ensures that the load of vehicle k equals the pickup quantity of all customers which are severed by vehicle k.Equation ( 14) ensures that each vehicle's load cannot exceed its maximum load, and (15) ensures that each vehicle's travel distance is not more than the maximum distance.

The Proposed Algorithm BEG-NSGA-II
3.1.Framework.In this paper, a two-stage optimization mechanism is used to solve the proposed VRPSPD with strict constraints.In the first stage, we use an improved bee evolutionary genetic algorithm to optimize the initial population.The optimization efficiency of the first stage is improved by optimizing the selection operators, selecting the different crossover operators according to the similarity of individuals' parent chromosomes, and selecting the different mutation operators according to the performance of individual's parent chromosome.In the second stage, an improved NSGA-II is used to optimize the proposed VRPSPD model.Based on the improved crossover and mutation operators, we construct the methods of deleting duplicate individuals, introducing new individuals, and using elite population instead of the parent population to improve the population diversity.To deal with the strict constraints, we introduce an external auxiliary population.During the optimizing process, if the constraint violation degree of the infeasible solutions is smaller than a given value, these infeasible solutions are copied to the external auxiliary population and used to evolve in the next generation.By using this external auxiliary population, the infeasible solutions gradually evolve to the boundary of feasible solutions, and the convergence speed is accelerated.The framework of the proposed algorithm is shown in Figure 2.More details of steps are described as follows.
The first stage includes the following.
Step 1. Randomly generate N different chromosomes to generate the initial population, and set t = 1, where N is the number of chromosomes and t is the current generation.
Step 2. Calculate the constraint violation value of each chromosome with (17).If the number of individuals that meets the constraint conditions is greater than the specified value or t > T (T is the maximal iteration time in the first stage), jump over the first phase; otherwise go to step 3.
where Drift pi denotes the ith constraint condition value of chromosome P after decoding; Cap i denotes the upper bound of the ith constraint condition.
Step 3. Calculate each chromosome's fitness based on the degree of constraint violation using (18); then select the best fitness individual as the queen; if t > 1, compare the fitness value of the queen with the parent queen, and select the better one as the new queen; if the number of feasible solutions is more than one, use nondominated sorting and congestion degree computing to deal with the feasible solutions to ensure that the queen selected from the parent population is the best one, then compare it with the parent queen, and select the better one as the new queen.
Step 4. Select  individuals from the population using roulette wheel method.The number of  is calculated by (19), where round [ ] means the operation of round; for example, round Step 5. Randomly generate  individuals different from each other, and combine them with the individuals generated from step 4. The number of  is calculated by Step 6.The new queen takes crossover and mutation with other individuals generated from step 5 to produce the offspring population.
Step 7. Delete the redundant repeating individuals and randomly generate M individuals different from each other into the population, where M equals the number of redundant repeating individuals that were deleted.
Step 8.Return to step 2, and set t = t + 1.
The second stage includes the following.
Step 9. Take fast nondominated sorting and crowding degree calculation on population.
Step 10.Carry out the following three kinds of operation for the population generated by step 9: select the top S (S=0.3N)individuals as the elite population; copy the individuals whose constraint violation values are less than a given value to the external auxiliary population  * ; use a binary tournament selection method to select crossover individuals and they were taken crossover according to the corresponding selected crossover methods to generate offspring.Step 11.The individuals take mutation with the probability of 0.1 after crossover.If the offspring rank = 1, then select the conventional mutation; otherwise the two-binding mutation (TBM) or reverse mutation is selected randomly.
Step 12. Randomly generate  individuals different from each other.The number of  is calculated by where t denotes the current iterations and GEN denotes the total iterations.Round [ ] denotes rounding operating; if  < 5, set  = 5.
Step 13.Combine the individuals of S, S * , and  and the individuals taken from step 11.
Step 14. Use step 7 and step 9 to deal with the individuals from step 13, and select top N individuals as the offspring population.Step 15.If t > GEN, output the result; otherwise, t = t + 1, and return step 10.

Chromosome Representation.
According to the characteristics of logistics path planning, a string with different integers (the length equals L, where L is the number of customers) is used to express a chromosome.Then n k -1 zeros are inserted into this string, where n k is the number of distribution vehicles and zeros express the distribution centers.For example, if there are 6 customers and three vehicles to serve them, the string of 25603104 expresses the distribution paths which are path 0-2-5-6-0, path 0-3-1-0, and path 0-4-0.

Crossover Operators.
In this paper, in order to enhance the searching space and avoid prematurity of local optimal solutions, we use two types of crossover operators (singlepoint crossover and two-point crossover operator) in this proposed algorithm.The single-point crossover operator helps to play excellent high genetic characteristics and to improve the convergence speed.Nevertheless, when the two parent generations become identical or similar, the crossover offspring almost get no change especially in the middle and later periods of the evolution.Hence, a two-point crossover operator is proposed for its ability to improve the mentioned problems.Equation ( 22) is a difference degree function.The crossover operator is chosen by (23), where L is the length of the chromosome and a il denotes the first gene of ith chromosome.
(, ) ≤ , Single-point crossover  (, ) > , Two-point crossover (23) The procedure of single-point crossover is described as follows (P 1 and P 2 are used to denote two parents; O 1 and O 2 are used to denote two offspring).
Step 1.A random parameter k that meets the inequality 0<k< L (the length of the chromosome) is generated to determine the position of the crossover.
Step 2. The elements from 1 to k in P 1 are duplicated to O 1 in the same positions; and the elements from 1 to k in P 2 are duplicated to O 2 in the same positions.
Step 3. The elements in P 2 but not in O 1 are duplicated to the remaining empty positions in O 1 from left to right.
Step 4. The elements in P 1 but not in O 2 are duplicated to the remaining empty positions in O 2 from left to right.
The procedure of two-point crossover is described as follows (P 1 and P 2 are used to denote two parents; O 1 and O 2 are used to denote two offspring).
Step 1. Two random parameters k 1 and k 2 that meet the inequality 0 < k 1 < k 2 < L as well as k 1 ̸ = k 2 are generated to determine the positions of crossover.
Step 2. The elements from k 1 to k 2 in P 1 are appended to the leftmost positions of O 1 ; and the elements from k 1 to k 2 in P 2 are appended to the leftmost positions of O 2 .
Step 3. The elements in P 2 but not in O 1 are duplicated to the remaining empty positions in O 1 from left to right.
Step 4. The elements in P 1 but not in O 2 are duplicated to the remaining empty positions in O 2 from left to right.
The examples of singe-point crossover and two-point crossover are, respectively, shown in Figures 3 and 4.

Mutation Operators.
In this paper, we have adopted three different mutation operators for the purpose of expending the solution space as well as maintaining the good solutions, which are conventional mutation operator, reverse mutation operator, and TBM operator.The examples of these mutation operators are shown in Figures 5, 6, and 7.At this stage, every offspring inherits the excellent traits of the queen for it is generated by the mutation of individuals and queen.So we randomly chose a mutation operator with P m probability at this stage.At the second stage, if the rank of offspring equals 1, we select the conventional mutation; otherwise we randomly select reverse or TBM mutation.
The main procedure of conventional mutation operator is described as follows (P 1 and O 1 are used to denote a parent and offspring, respectively).
Step 1. Randomly select two positions in P 1 .
Step 2. Swap the elements in the selected positions to generate O 1 .
The main procedure of reverse mutation operator is described as follows (P 1 and O 1 are used to denote a parent and offspring, respectively).
Step 1. Randomly select two positions in P 1 .
Step 2. Reverse the numbers between the two selected positions to generate O 1 .
The main procedure of TBM operator is described as follows (P 1 and O 1 are used to denote a parent and offspring, respectively).
Step 1.A random parameter m (m <L -3, where L is the length of the chromosome) is generated in P 1 .
Step 2. Exchange the elements m with m+3 and m+1 with m+4 in P 1 to generate O 1 .

Experimental Results
The proposed BEG-NSGA-II algorithm was coded in MAT-LAB R2014a and implemented on a computer configured with Intel Core i3 CPU with 2.67 GHz frequency and 8GB RAM.Three sets of examples are used to illustrate the performance of the proposed algorithm, in which three of the same trucks were used to deliver and pick up the goods.The first one is the VRPSPD problem of an actual engineering machinery remanufacturing company (with 12 customers).To further verify the effectiveness of the proposed algorithm, another two instances with 22 customers and 42 customers are used to simulate the problem of VRPSPD.
Because the problem of this paper has no benchmarks, we generate the test instances based on the SCA3-0 instance from Dethloff [40].The generating standards are described as follows: n customers are randomly generated that are distributed in a square area whose length is 100 km; d i is randomly generated within [0.5ton, 5ton]; p i = 0.5(d i * r+1), where r is randomly generated within [0, 1]; the number of distributors, recyclers, and suppliers is randomly generated within [d i , 2 *  i ], and the total number of distributors, recyclers, and suppliers equals the total customers that we set (such as 12 customers).
The values of parameters are shown in Table 2 [30].
In addition, G m =25ton, G 0 =10ton, and =1, and the best value of v ij is 17m/s.The adopted parameters of the proposed algorithm are set as follows: N = 100, p c = 1, p m = 0.1, T = 100, and GEN = 1000.In order to demonstrate the effectiveness and performance of the proposed algorithm, the traditional algorithm and the proposed algorithm are run 10 times, respectively.3 shows an example of VRPSPD problem of an actual engineering machinery remanufacturing company with 12 customers.The first column shows the information of the customers (there are 1 recycling center, 6 distributers, 3 recyclers, and 3 suppliers).In the second column, in order to facilitate encoding and decoding in   the algorithm, all customers were numbered as a series continuous integer (recycling center 0 is numbered as 0; distributors 1, 2, ..., and 6 are numbered as 1, 2, ..., and 6; recyclers 1, 2, and 3 are numbered as 7, 8, and 9; suppliers 1, 2, and 3 are numbered as 10, 11, and 12.).The third and fourth columns show the X/km and Y/km coordinates of the customers, respectively.The amounts of demand/ton and recycle/ton of the customers are shown in columns 5 and 6.What is more, column 7 determines the service time of the customers (for example, 10 am-3 pm means that the service time of the customer is needed from 10 am until 3 pm).

Test Example 1. Table
In order to demonstrate the performance of the proposed algorithm and research the effect of the duplicate individuals, test example 1 was tested by the traditional NSGA-II (T-NSGA-II), BEG-NSGA-II without eliminating the duplicate individuals (W-BEG-NSGA-II), and BEG-NSGA-II, respectively.Table 4 shows the optimization results of test example 1 using these three algorithms.From the optimization results, we can gain the following information: (i) The algorithms of T-NSGA-II, W-BEG-NSGA-II, and BEG-NSGA-II could find 5, 6, and 6 optimal Pareto solutions with average computing time of 1.2, 0.9, and 0.5 minutes correspondingly during 10 run times.
(ii) The algorithm of BEG-NSGA-II could find 2 Pareto solutions (bold and marked as " a " in Table 4) that dominate the results of T-NSGA-II and 6 Pareto solutions (bold and marked as " b " in Table 4) that dominate the results of W-BEG-NSGA-II.
From the information mentioned above, we can draw the following conclusions: (i) The algorithms of T-NSGA-II and W-BEG-NSGA-II may generate amount of duplicate individuals, especially after 100 iterations.These duplicate individuals may cause the problems of inefficiency, slow convergence, and converging to local Pareto optimal solutions.
(ii) The proposed algorithm of BEG-NSGA-II is more effective than the algorithms of T-NSGA-II and W-BEG-NSGA-II to solve the VRPSPD problem with rigor constraints.

Test Example 2. Table
In order to demonstrate the performance of the proposed algorithm and research the effect of the duplicate individuals, test example 2 was also tested by T-NSGA-II, W-BEG-NSGA-II, and BEG-NSGA-II, respectively.Table 6 shows the optimization results of test example 2 using these three algorithms.From the optimization results, we can gain the following information: (i) The algorithms of T-NSGA-II, W-BEG-NSGA-II, and BEG-NSGA-II could find 10, 12, and 10 optimal Pareto solutions with average computing time of 2.5, 1.6, and 0.9 minutes correspondingly during 10 run times.
(ii) All of the results (bold and marked as " a " in Table 6) gained by T-NSGA-II dominate the ones gained by W-BEG-NSGA-II.
(iii) The algorithm of BEG-NSGA-II could find 6 Pareto (bold and marked as " a " or " b " in Table 6) solutions that dominate the results of T-NSGA-II as well as W-BEG-NSGA-II.
From the information mentioned above, we can draw the following conclusions: (i) The algorithms of T-NSGA-II and W-BEG-NSGA-II may generate amount of duplicate individuals, especially after 100 iterations.These duplicate individuals may cause the problems of inefficiency, slow convergence, and converging to local Pareto optimal solutions.
(ii) The proposed algorithm of BEG-NSGA-II is more effective than the algorithms of T-NSGA-II and W-BEG-NSGA-II to solve the VRPSPD problem with rigor constraints.

Test Example 3. Table
(ii) The algorithm of BEG-NSGA-II could find 9 Pareto solutions (bold and marked as " a " in Table 8) that dominate the results of T-BEG-NSGA-II.
From the information mentioned above, we can draw the following conclusions:  to test example 1. Figure 16 shows the number of duplicate individuals which changes with the number of iterations in test example 3.

Test Example 4.
Dethloff 's benchmark [40] which has been widely used in these years is adopted to illustrate the pervasiveness and superiority of the proposed BEG-NSGA-II.To the best of our knowledge, the best known upper bounds for Dethloff 's benchmark [40] were found by the following algorithms: PILS: parallel iterated local search [41] VLBR: variable length bone route [42] HPSO: hybrid particle swarm optimization [23] ACS: ant colony system [43] ACSEVNS: ant colony system empowered variable neighborhood search algorithm [44].
Therefore, we compare the performance of the proposed BEG-NSGA-II against these algorithms reported in the literature.Table 9 shows the optimization results of test example 4 using these algorithms and our algorithm.The optimization objective of PILS, VLBR, HPSO, and ACS is vehicle distance, by which the best known solution (BKS) for Dethloff 's benchmark [40] is found.The optimization objective of ACSEVNS is also vehicle distance, which shows the best solutions (Best), average solution over 10 replications (Avg.), and average computation time (T).Because the waiting time of customers is not considered in Dethloff 's benchmark [40], our proposed BEG-NSGA-II only optimizes the objectives of vehicle distance and fuel consuming (Fuel).From the optimization results, we can gain the following information: (i) The algorithm of BEG-NSGA-II could find 9 BKSs in 40 instance.(ii) The objective of vehicle distance and computation time obtained from algorithm of BEG-NSGA-II is little worse than the ones obtained from ACSEVNS.

Algorithms
(iii) In addition to the objective of vehicle distance and computation time, the optimization objective of fuel consumption that is not considered in the compared algorithm is also obtained in the proposed algorithm of BEG-NSGA-II.
From the information mentioned above, we can draw the following conclusions: (i) The proposed algorithm of BEG-NSGA-II can obtain the BKS or near BKS in reasonable computation time considering the objectives of vehicle distance and fuel consumption simultaneously.
(ii) The reason why some results gained from proposed BEG-NSGA-II are little worse than the compared algorithms is that we consider the optimization objectives of vehicle distance and computation simultaneously, while the compared ones only consider the single objective of vehicle distance.
(iii) From the simulations of examples 1-4, we can see that our proposed algorithm of BEG-NSGA-II has high performance in multiobjective optimization as well as meeting the single objective optimization.

Conclusions
In this paper, we construct a new closed-loop supply chain logistics network of VRPSPD dominated by remanufacturer, in which the customers are originally divided into three types: distributors, recyclers, and suppliers.The proposed model is more consistent with the actual situation of recycling and remanufacturing.For the purpose of resource saving and environment protection, we originally add the fuel consumption as an optimization objective into VRPSPD.In order to solve the proposed model and conquer the disadvantages of the traditional NSGA-II and the existing algorithms with two-stage optimization mechanism, a bee evolutionary algorithm guiding nondominated sorting genetic algorithm II (BEG-NSGA-II) is proposed with the optimization objectives of minimum fuel consumption, minimum waiting time, and the shortest delivery distance while meeting the travel constraints and load constraints at the same time.To verify the effectiveness of the proposed algorithm, test examples are used to simulate the proposed VRPSPD and the traditional VRPSPD, respectively.From the simulation results of test examples 1-3, we can see that our proposed algorithm BEG-NSGA-II is superior to the traditional NSGA-II (T-NSGA-II) and the proposed BEG-NSGA-II without eliminating the duplicate individuals (W-BEG-NSGA-II).From the simulation results of test example 4, we can see that the proposed algorithm of BEG-NSGA-II can obtain the best known solution (BKS) or near BKS in reasonable computation time considering the objectives of vehicle distance and fuel consumption simultaneously.That is to say, our proposed algorithm of BEG-NSGA-II works better in multiobjective optimization as well as meeting the single objective optimization.It would be prosperous to apply our algorithm to other multiobjective optimization problems, especially in the logistics distribution field.
(i) G m : the vehicle capacity (ii) G 0 : the weight of an empty vehicle (iii) L: the maximum travel distance of a vehicle (iv) d ij : the distance between any two service points i and j (v) v ij : the speed between any two service points i and j (vi) d i : the demand quantity of service point i (vii) p i : the pickup quantity that service point i has (viii) K: the set of vehicles (ix) k: the vehicle of k, k∈K (x) n k : the number of vehicles
Note: (a) " a "meant that the marked result dominates the results of T-NSGA-II.(b) " b " meant that the marked result dominates the results of W-BEG-NSGA-II.(c) " a,b " meant that the marked result dominates the results of T-NSGA-II and W-BEG-NSGA-II simultaneously.

Note:
(a) Cen., Dis., Rec., and Sup. are the abbreviations of recycling center, distributors, recycler, and supplier, respectively.(b) "-" means that these data do not need to be set in advance.

Figure 10 :
Figure 10: The number of duplicate individuals in test example 1 with 12 customers.
" a " meant that the marked result dominates the results of T-NSGA-II.(b) " b " meant that the marked result dominates the results of W-BEG-NSGA-II.(c) " a,b " meant that the marked result dominates the results of T-NSGA-II and W-BEG-NSGA-II simultaneously.

Figure 16 :
Figure 16: The number of duplicate individuals in test example 3 with 42 customers.
Note:(a) " a " means that the marked result dominates the results of T-NSGA-II.

Table 2 :
The values of parameters used in the proposed VRPSPD model.

Table 3 :
The basic information of test example 1 with 12 customers.
Note: (a) Cen., Dis., Rec., and Sup. are the abbreviations of recycling center, distributors, recycler, and supplier, respectively.(b) "-" means that these data do not need to be set in advance.

Table 4 :
Comparisons of the optimization results of test example 1 with 12 customers.

Table 5 :
The basic information of test example 1 with 22 customers.

Table 7 :
The basic information of test example 3 with 42 customers.

Table 8 :
Comparisons of the optimization results of test example 3 with 42 customers.
Note: the best solution of each problem instance is highlighted in bold.BKS: the best known solution.Ref.: the best solution reference.Best: the best solution in 10 replications.Avg.: average solution over 10 replications.T: average computation time.G. Avg.: average of 40 instances.BKS found: number of best solutions found.