A Hybrid Discrete Differential Evolution Algorithm to Solve the Split Delivery Vehicle Routing Problem

For the split delivery vehicle routing problem in the distribution link of logistics where the optimization objective is to minimize the total distance of the transportation path, a hybrid discrete differential evolution (HDDE) algorithm is adopted in this study. The algorithm applies the encoding and decoding method based on the remaining vehicle capacity to split the demand of customer point and adjust the splitting node to the optimal position of its path. In order to avoid the prematurity of the algorithm, the perturbation mechanism is added, i.e., the three-point insertion method is used in the crossover stage. To ensure the higher quality of the solution, the idea of tabu search is introduced in three cases: initial population generation, selection operation and local optimal situation. Finally, through the numerical experiments based on benchmark data, the HDDE algorithm is compared with other competitive algorithms to demonstrate its domination.


I. INTRODUCTION
With a mass outbreak of COVID-19, medical resources are seriously deficient, and hospitals have higher requirements for the transportation speed of medical supplies. More and more hospitals introduce intelligent logistics equipment, and intelligent transport vehicles such as unmanned delivery vehicle and unmanned warehouse transport vehicle have begun to enter people's vision. In order to minimize human contact and reduce the possibility of human-to-human transmission of the virus, unmanned delivery vehicles transport medical supplies to many hospitals in urban areas. The demand of each hospital may be of great quality or volume, even exceed the transport capacity of unmanned vehicles. To reduce the empty load rate and improve the utilization rate of transport The associate editor coordinating the review of this manuscript and approving it for publication was Huaqing Li . resources, the supplies of hospitals in great demand can be split and distributed to multiple unmanned delivery vehicles to transport. Such transport problem can be classified as the split delivery vehicle routing problem (SDVRP). Dror and Trudeau [1] first proposed the split delivery vehicle routing problem, a relaxation of the classical capacitated vehicle routing problem. In its seminal paper, Dror and Trudeau found two characteristics of the optimal solution if the travel cost of any pair of customer points satisfies the triangle inequality: (1) at most one common customer exists in any two paths; and (2) k-split (k ≥ 2) cycle does not exist, where k-split cycle is defined as follows: for k demand points 1, 2,. . . , k and k routes, the points 1 and 2 are included in route 1, the points 2 and 3 are included in route 2,. . . , the points k-1 and k are included in route k − 1, and the points k and 1 are included in route k. This subset of k demand points is called a k-split cycle. Then in 1990, the two authors [2] first presented the mixed integer linear programming model of SDVRP and demonstrated that the number of vehicles and the travel distance can be significantly reduced by splitting if demand is over 10 percent of vehicle capacity. Belenguer et al. [3] indicated that SDVRP has a strict lower bound if the customer demand is less than the vehicle loading capacity and the actual delivery volume is integer. Archetti and Mansini et al. [4] pointed out that when the distance values satisfy the triangle inequality and the vehicle loading capacity Q = 2, the solution of SDVRP can be obtained in polynomial time under the general cost condition. When Q > 2, SDVRP is NPhard. Archetti and Savelsbergh et al. [5] studied the situation where the demand of customers exceeds the loading capacity of vehicle and analyzed the cost savings brought by SDVRP from the aspects of vehicle number and driving distance. They proved that the relationship between average customer demand and vehicle loading capacity was the biggest factor affecting cost savings in SDVRP. The cost savings brought by SDVRP are not related to the distribution of customer points but are affected by the variance of customer demands. Maximum benefit can be achieved in the case that the average demand is slightly more than 50% of the vehicle capacity and customer demands have a smaller variance [6].
Since SDVRP is one of the NP-hard problems, its solution mainly devoted to two aspects, i.e., exact algorithms and heuristic ones. For the former aspect, Dror et al. [7] added sub-tour elimination constraints and designed a branch and bound algorithm with proper constraint combinations, and considerably reduced the differences between the lower and upper bounds at the root of the branching tree. Then, Lee et al. [8] considered SDVRP as a dynamic programming problem and proposed a shortest path search method, greatly reducing the solution space. Jin et al. [9] proposed a two-phase method equipped with valid inequality constraints including customer point clustering stage and TSP path arrangement stage, and their algorithm significantly reduce iterations and computational time required to reach an optimal solution. Moreno et al. [10] designed a column generation and tangent plane method which obtain a more effective lower bound. Archetti and Bianchessi et al. [11] classified and defined the split delivery vehicle routing problem with unlimited fleet (SDVRP-UF) and the split delivery vehicle routing problem with limited fleet (SDVRP-LF), presenting a branchand-price-and-cut method decomposing the problem and a heuristic whose initial solution results from the solution of subproblem (column) to get the upper bound. They presented two branch-and-cut algorithms based on two relaxed formulas that find the lower bound [12].
However, due to the computational limitation, exact algorithms can only solve small instance of SDVRP. Therefore, heuristic ones, especially meta-heuristic algorithms, have become the main research direction of SDVRP. Archetti et al. [13] designed a tabu search algorithm based on the insertion neighborhood construction method which is proved to be more efficient than the previous heuristic algorithm proposed by Dror et al. Aleman and Hill et al [14] presented the tabu search with vocabulary building approach (TSVBA) which is a learning process, i.e., new solutions are generated from a set of initial solutions according to common properties and novel features that can generate superior solutions. Derigs [15] used a variety of neighborhood operators and compares several meta-heuristic algorithms, proving that the attributebased hill climber heuristic has the best performance and the relocate neighborhood operator has the greatest contribution. Wilck and Cavalier [16] designed two kinds of hybrid genetic algorithms, one of which is based on the shortest path algorithm, while the other is combined with the ratio relationship of unit customer demand and unit distance. Berbotto [17] introduced granularity computing technique into tabu search algorithm, presenting the idea of exploring neighborhoods of infeasible solutions that does not satisfy the capacity constraint. Silva et al. [18] proposed an iterative local search heuristic with multiple initial positions, including a novel perturbation mechanism for both two cases of SDVRP-UF and SDVRP-LF. Zhang et al. [19] designed a forest-based tabu search Algorithm with three neighborhood operators, applying a novel solution representation to combine a routes and forest. Rajappa et al. [20] presented an ant colony optimization (ACO) and a hybridization of ACO, genetic algorithm and a heuristic. Shi et al. [21] designed a particle swarm optimization approach including two sets of local searches for optimum solutions, which can determine whether they are performed on a specific solution.
Compared with above heuristic algorithms, Differential Evolution (DE) algorithm has four advantages, which are simple structure, easy implementation, quick convergence, as well as robustness. Given that the idea of hybridization of DE and tabu search have been successfully applied to solve the job-shop scheduling problem by Ponsich and Coello [22], it may be also beneficial to solve SDVRP. To our best knowledge, there is no known research applied DE algorithm to solve the SDVRP, therefore, we propose a hybrid discrete differential evolutionary (HDDE) algorithm in this paper.
The rest of this paper is organized as follows. The mathematical formulation of the SDVRP is provided in Section II, and our HDDE algorithm is presented in Section III. Section IV is devoted to the numerical experiments. Finally, Section V concludes this paper.

II. PROBLEM STATEMENT
The SDVRP studied in this paper is set as single depot, single vehicle type, no time window limit, and pure delivery, which can be described as follows. There are n customer points, which would be served by multiple vehicles with the maximum load of Q. All vehicles start from a certain depot and return it after finishing the delivery task. Each customer demand can be delivered by one or more vehicles, and the objective is to gain the path of each vehicle that the total travel distance is shortest. The following notations in table 1 are used to describe the problem.
The problem is based on the following assumptions. (1) The distance between the two points is symmetric: c ij = c ji .
(2) Distance satisfies the triangle inequality: c ij +c jk ≥ c ik . The mathematical programming model of the problem is provided as follows: Subject to: i∈S j∈S The objective is represented by Equation (1), which minimizes the total travel distance. Equation (2) represents the minimum number of vehicles required. Constraint (3) restricts that each customer and the depot must be visited at least once. Constraint (4) ensures flow conservation con-straint that the number of visiting a point is the same as the number of leaving the point in each path. Constraint (5) limits that customer point is serviced if and only if a vehicle visited it, and its demand assigned to a vehicle is less than or equal to its total demand. Constraint (6)- (7) guarantees that all customer demand should be satisfied, and the actual load of any vehicle will not exceed the maximum load. Constraint (8) means subtour elimination. Constraints (9)-(10) provides variable restrictions.

III. A HYBRID DISCRETE DIFFERNETIAL EVOLUTION ALGORITHM FOR SDVRP
The differential evolution algorithm presented by Storn and Price [23] in 1997 is one of the global population evolutionary optimization algorithms to solve continuous optimization problems. The basic idea of this algorithm is to start from the initial population and generate new individuals through crossover and mutation operations. New individuals can only be selected as the next generation if their fitness is superior to that of their parents. For applying in a discrete problem, the individual in discrete differential evolutionary (DDE) algorithm is expressed as a customer-based sequence. In order to overcome pre-maturity and obtain a better solution, a hybrid discrete differential evolutionary (HDDE) algorithm is presented, combined with tabu search operation in initial population generation, selection operation and the phase at which the solution cannot be improved, and three-point insertion method is adopted in the crossover stage. The procedure of the HDDE algorithm is as follows.

A. ENCODING AND DECODING
The customer point sequence called the target individual in HDDE is a potential solution. Each individual is encoded as a full array of customer points, and then decoded according to its order by the following three phases: preliminary path generation phase according to capacity constraint, spliting point position adjustment stage and distance calculation stage. Some related variables are showed in table 2. (1) Decoding phase 1: generate preliminary paths through capacity constraint relation.
Since distances conform to the relation of triangle inequality, the demand of some customer points should not be split to shorten the total travel distance. In the decoding phase, put customer point one by one into the path according to the order of customer points in an individual, and judge when the path capacity is exceeded: if Q c ≤ Q r , customer point is put into next path; otherwise customer point's demand is split.
In paths obtained from phase 1, the splitting point can only be located at the first or the last in customer sequence not including depot due to the encoding method. In this phase, the splitting point's location will be adjusted to the best location in its path to try to reduce distance and two cases are discussed as follows. The best location is where c is smallest.
Case 1: |SPc| = 1, i.e., a path which is being adjusted exists one splitting point at the first or the last location. For If the change value in the first position is the smallest, the splitting point When x 1 is front of x n * , put x 1 into candidate insertion location λ 1 (1≤ λ 1 ≤ n * − 1) and put x n * into candidate insertion location λ n * (λ 1 ≤ λ n * ≤ n * − 1), then calculate distance change value. If . Analogously, when x n * is front of x 1 , put x n * into candidate insertion location λ n * (1≤ λ n * ≤ n * -1) and put x 1 into candidate insertion location λ 1 (λ n * ≤ λ 1 ≤ n * -1), . Finally, find the best position and put x 1 and x n * into respective best position.
(3) Decoding phase 3: calculate total distance of the paths obtained by phase 2.
For an individual X = [x 1 , x 2 , . . . , x n ], the detailed process of decoding is as follows: Step 1: initialization: Step 3; otherwise: go to Step 4.
Step 3: set i = i + 1, and if i > n: go to Step 6; otherwise: set d = d x i then return to Step 2.
Step 4: if Q c ≤ Q r : set j = j + 1, Q r = Q r − Q c , Q c = Q, then return to Step 2; otherwise: go to Step 5; Step 5: customer point's demand is split and use up the remaining capacity of current path. Then set d = d−Q c , j = j + 1, Q c = Q and return Step 2.
Step 6: start adjusting locations of splitting points and set j = 1.
Step 7: if |SPc| = 0 for current path j, go to Step 8; otherwise go to Step 9.
Step 9: for the splitting point in SPc, find the best adjusting position, then return to Step 8.
Step 10: calculate total travel distance then stop.   neighborhood search and a simulation of human memory function. Its basic idea is to use a tabu list to record the local optimal solution information found recently and to query tabu list in subsequent search to avoid roundabout search, to jump out of the local optimality and ensure a diversified and effective search [24].
In order to enhance the global optimization capability, Part of the idea of tabu search algorithm is introduced as an operation of HDDE. Before the formal description of HDDE, a brief introduction is presented on the tabu search operation (TSO). TSO consists of four following elements.

1) INITIAL INDIVIDUAL
Different from the previous neighborhood search methods that directly operate paths in the existing literatures, because the designed decoding method has met the capacity constraint, tabu search is performed directly for an individual, and one individual in the process of HDDE is designated as initial solution.

2) TABU OBJECTIVE AND TABU LENGTH
For convenience, the objective function value of the individual is tabued in a tabu list. Tabu length is measured experimentally in appendix.

3) NEIGHBORHOOD OPERATOR
Three neighbor operators including relocate operator, exchange operator and 2-opt operator are adopted to increase the diversity of neighborhood solutions (neighborhood individuals). Relocate operator is that a customer point is removed and inserted into another location and an example is presented  in figure 3: remove 3 and insert it between 7 and 1. Exchange operator is a switch of any two customer points' locations and figure 4 shows an example of switching the locations of 5 and 4. 2-opt operator is that the individual breaks at two points and forms three segments, with the middle subsequence reversed, and figure 5 illustrates an example: reverse subsequence 5-7-1-4. Since the capacity relation is judged when decoding, the three neighborhood operators do not need to judge the capacity relation, which simplifies the operation.
To save running time, a series of candidate solutions not tabued are generated by one operator selected randomly of three. After that, the object value of the selected candidate solution is placed into the tabu list.

4) THE NUMBER OF NEIGHBORHOOD SOLUTIONS SELECTED (N S )
Unlike the classical tabu search algorithm, TSO doesn't carry out multiple iterations. After neighborhood operator searching about the initial solution once, a certain number of candidate neighborhood solutions are selected depending on when HDDE is applied.
In subsequent sections, the symbol TSO −X i − n s − XS is applied to indicate the use of TSO, representing that TSO start from the initial individual X i to find the n s best neighborhood individuals to update the individual set XS, where |XS| = n s .

C. HYBRID DIFFERENTIAL DISCRETE EVOLUTION ALGORITHM
Two heuristic methods are introduced to generate the initial population. During each iteration, individuals in the current population are successively mutated, crossed and selected. The DE algorithm may stall in the later stage, so the distur-bance mechanism is designed. The detailed processes of these operations are shown in the following subsections.

1) INITIAL SOLUTION GENERATION
The target population consists of target individuals in which the hth target individual is X τ h = [x τ h,1 , x τ h,2 , . . . ,x τ h,n ], where h = 1, 2, . . . , , τ is the current iteration and is the population size. Firstly, two initial individuals of the initial population are generated by sweep algorithm and nearest neighborhood algorithm.
a: X 0 1 BASED ON SWEEP ALGORITHM Given coordinates of depot and all customer points, a coordinate system is established with the depot as the origin, and an arc is sent out from the depot. The arc starts at an angle of zero and rotates counterclockwise, then the sequence that the customer points are scanned is taken as X 0 1 [25].
b: X 0 2 BASED ON NEAREST NEIGHBORHOOD ALGORITHM First, put the customer point closest to the depot into X 0 2 , and then put the one closest to the current last customer point into X 0 2 in turn, until all the customer points are put into X 0 2 [26]. Initial individuals X 0 1 and X 0 2 have been generated, and the remaining individuals are multiple neighborhood individuals of TSO based on the two initial individuals, where is assumed as the even number. Namely, the following two operations are performed: ( . . , . The purpose of interlacing generation is to evenly distribute the two groups of neighborhood solutions. Therefore, the objective value of individuals in the initial population can be calculated by the decoding method.

2) MUTATION OPERATION
Near-optimal solution is updated as the population iterates. The mechanism of mutation operation adopted in this study is DE/best/2. The mutant individual obtained by mutation

For each target individual, a mutant individual V τ
h is obtained by the following formula: where X τ −1 best is the best individual, X τ −1 α 1 , X τ −1 β 1 , X τ −1 α 2 and X τ −1 β 2 are four different target individuals at iteration (τ − 1), α 1 , β 1 , α 2 and β 2 are random integers in [1, ], and Z ∈ [0,1] is a mutation probability. The operator ⊗ is calculated as follows: . . ,g τ h,n ] is a temporary individual. The operator ⊕ is calculated as follows: where operator mod is modulus. Individual V τ h maybe not valid because some encoding elements may repeat or not exist. An instance is provided in Table 4-6 to explain the mutation process.

3) CROSSOVER OPERATION
Mutant individual can increase the perturbation of target individual and jump the local optimum [27]. The crossover operation is performed after mutation operation to obtain a trial individual U τ h = [u τ h,1 , u τ h,2 , . . . ,u τ h,n ] that is valid. In the crossover operation phase, selective parts of mutant individual V τ h are put into a target individual X τ −1 h to produce VOLUME 8, 2020 a trial individual U τ h . Three-point insertion method is used to increase the disturbance. The insertion points ρ1, ρ2 and ρ3 are three different positions of target individual, where three parts of mutant individual are inserted, and they are different random integers in [1, n]. For each encoding element of V τ h , if a random number rand() < Y , it can be placed into , where Y ∈ [0,1] is a crossover probability; otherwise, it is deleted. The encoding element appear more than once in U τ h are omitted from right to left, which controls that each customer number in U τ h exists once only. For V τ h in Table 6, the crossover process of X τ −1 h is provided in Table 7. The crossover probability Y = 0.7 and the insertion points are ρ1 = 1, ρ2 = 3 and ρ3 = 6.

4) SELECTION OPERATION
The new generation of individuals comes from trial individuals or the last generation's neighborhood individual. Trial individuals with better objective values can be selected as the next generation of individuals. To make full use of the neighborhood solutions generated in the iterative process, if a trial individual U τ h can't update the individual as the initial individual, perform TSO and reserve one best neighborhood individual to replace {X τ h }. The best individual of the current population may need to be updated.

5) DISTURBANCE OPERATION
The selection operation of the DDE algorithm adopts the greedy selection strategy, which can accelerate the convergence, but easily result in the local optimum. To enhance the global optimization ability, DDE algorithm is improved as follow: (1) If there is no improvement on best individual of current population within two consecutive iterations, perform TSO , with X τ best as the initial individual, perform TSO and reserve ( − 1) best neighborhood individuals to replace the current population except X τ best . It can enhance the local optimization ability and improve the quality of the whole population once.
(2) If there is no improvement on best individual of current population within ten consecutive iterations, perform TSO − , with a random individual X τ rand as the initial individual, perform TSO and reserve ( − 1) best neighborhood individuals to replace the current population except X τ rand . Therefore, the local optimal value can be jumped out by disturbing the population.
In the above two cases, the best individual of the current population may need to be updated.

D. ALGORITHM FLOWCHART
The concrete procedures of HDDE are expressed as follows [28].
Step 1: Initialize parameters including mutation probability Z , crossover probability Y , population size , tabu list length TL and maximum iteration τ max .
Step 2.1: Execute two heuristic method to obtain two sequences including all costumers.
Step 3: Set h = 1 and start iteration.
Step 5: Perform the mutation and crossover to obtain the trial individual U τ h .
Step 8: Update X τ best . If there is no improvement of current population in two consecutive iterations, perform TSO − X τ best -( − 1) − ({X τ h } h=1 − X τ best ) and update X τ best , otherwise go to Step 10.
Step 9: Update X τ best . If there is no improvement of current population in ten consecutive iterations, perform TSO − X τ rand -( − 1) − ({X τ h } h=1 − X τ rand ) and update X τ best .
Step 11: Output the historical optimal sequence as the final solution. Stop.
The flowchart of HDDE is expressed in figure 6.

IV. COMPUTATIONAL RESULTS
To demonstrate the domination of the proposed HDDE algorithm, experiments with the following two sets of SDVRP benchmark instances were conducted: Set 1: this set contains 25 instances generated by Belenguer et al. [3]. The number of customers ranges from 21 to 100. The coordinates originate in TSPLIB instances. The customer demand is randomly generated through uniform distribu- Set 2: this set contains 42 instances generated by Archetti et al. [11], derived from the capacitated VRP (CVRP) bench- mark instances. The number of customers ranges from 50 to 199, and the generation scenarios of customer demands are the same with set 1.
In this paper, HDDE is compared with two competitive algorithms: (1) the forest-based tabu search (FBTS) algorithm by Zhang et al. [19]; and (2) Genetic algorithm (GA). The proposed HDDE algorithm and GA algorithm were coded in C++ language and run on a PC with an Intel Core i7 CPU (2.6GHz×12) and 8 GB RAM. The parameters of HDDE is tested through orthogonal experimental design: mutation probability: 0.9, crossover probability: 0.6, population size: 50, tabu list length: 300, maximum iterations: 300. The specific process of orthogonal experimental design is shown in appendix.
The experimental results are shown in the tables 8 and 9. The results of FBTS were obtained directly from the literature [19]. The initial population generation method of GA is the same as HDDE. Because the comparation algorithms were carried out on different machines, running times of FBTS and GA are not reported. Column ''Cost'' of HDDE, column ''FBTS'' and column ''GA'' respectively show final objective values (final total travel distance) of three algorithms.   In the experimental results of the two data sets, the percentage of cases that HDDE is superior to FBTS and GA are 61.2% and 100%, respectively. The average GAP 1 between the two data sets is 0.49%, where set 1 is 0.14% and set 2 is 0.70%. Although these gaps are all less than 1%, considering that SDVRP has been extensively studied over the past 30 years, the improvements of the proposed HDDE are considerable. The average GAP 2 between the two data sets is 5.96%, where set 1 is 6.67% and set 2 is 5.53%. It is found that HDDE algorithm is much better than GA. Therefore, HDDE algorithm designed in this study has achieved good results.  The HDDE algorithm that employs the tabu search operator provides a good balance between the global and local optimization capabilities. Moreover, the capacity constraint is only placed in the decoding part, which simplifies the operation of neighborhood operators. However, the running time of this algorithm is relatively long. The time is mainly spent on the path distance calculation during the decoding phase, so path structure and programming skills need to be further improved.

V. CONCLUSION
SDVRP relieves the constraint that one customer can only be accessed by only one vehicle, which is conducive to save vehicle resources and shorten the travel distance. However, this kind of problem is more complicated than the traditional VRP problem, and solution method is quite different. This article investigates SDVRP to minimize total travel distance. For this NP-hard problem, a hybrid discrete differential evolution algorithm is presented to obtain near-optimal solution of high quality. Numerical simulations verify the domination of the HDDE.
Future research could mainly focus on two aspects: (1) Customer point location will be considered to determine whether a customer's demand should be split so that further optimization can be carried out in the decoding stage; and (2) future research will consider the additions of time windows, multi-types vehicle and other practical common factors to design an appropriate algorithm.

Parameters setting for HDDE:
The five factors that affect the performance of HDDE are tested. Factors and levels needed to test are provided in table 10 and the orthogonal experimental result is shown in table 11. Figure 7 indicate mean of experimental results at each level. Each group of experiment is tested randomly ten times: n = 50 and improved percentage of the algorithm is recorded. The column ''Experimental result'' is the average improved percentage of 10 random tests. Improved percentage Z IS −Z FS Z FS ×100% is employed, where Z IS and Z FS denote the objective value of the initial solution and the final solution obtained by HDDE respectively.
According to the experimental results in figure 7, the value of each parameter is determined as: mutation probability: 0.9, crossover probability: 0.6, population size: 50, maximum iteration: 300, tabu list length: 300. YUANYUAN