Time-Dependent Multi-Depot Heterogeneous Vehicle Routing Problem Considering Temporal–Spatial Distance

: Aiming at the multi-depot heterogeneous vehicle routing problem under the time-dependent road network and soft time window, considering vehicle ﬁxed cost, time window penalty cost and vehicle transportation cost, an optimization model of time-dependent multi-depot heterogeneous vehicle routing problem is established with the objective of minimizing distribution cost. According to the characteristics of the problem, a hybrid genetic algorithm with variable neighborhood search considering the temporal–spatial distance is designed. Customers are clustered according to the temporal–spatial distance to generate initial solutions, which improves the quality of the algorithm. The depth search capability of the variable neighborhood search algorithm is applied to the local search strategy of the genetic algorithm to enhance the local search capability of the algorithm. An adaptive neighborhood search number strategy and a new acceptance mechanism of simulated annealing are proposed to balance the breadth and depth required for population evolution. The validity of the model and algorithm is veriﬁed by several sets of examples of different scales. The research results not only deepen and expand the relevant research on vehicle routing problem, but also provide theoretical basis for logistics enterprises to optimize distribution scheme.


Introduction
Time-dependent multi-depot heterogeneous vehicle routing problem is in the condition of road traffic changing and multiple depots joint distribution, at the same time, considering the customer's service time, spatial location, heterogeneous vehicles, fuel consumption and other factors. In reality, tableware suppliers according to customers' demand sent a variety of vehicles to deliver. Express delivery industry salesman arranges vehicles reasonably according to regional express quantity. At the same time, the speed of vehicles is affected by the real-time traffic conditions of the route. The speed of distribution vehicles is slow in peak hours, but fast in off-peak hours. This paper covers two hot issues in current vehicle routing problem (VRP) research: time-dependent VRP (TDVRP) and multi-depot VRP (MDVRP). Firstly, TDVRP is proposed by Beasley [1]. Its research focused on vehicle travel time affected by vehicle departure time; the follow-up researches of TDVRP mostly adopted the time-dependent function based on travel speed proposed by Ichoua et al. [2], the vehicle traveling speed is expressed as a piecewise function [3,4]. The Federal Highway Administration [5] mentioned that the speed of a vehicle changes gently, rather than having a step change at a certain point. Moreover, the multi-depot joint distribution mode can optimize the distribution scheme through global information, it can effectively solve the problem of high transportation cost and low service level gradually appearing in the independent distribution mode of single distribution depot [6], so the MDVRP has received more and more attention.

Literature Review
As mentioned previously, a number of scholars have studied the TDVRP. For TDVRP, Sabar et al. [7] established a mathematical model aiming at minimizing total distribution cost with different traffic conditions in different time periods, and proposed an adaptive evolutionary algorithm using variable parameters to solve it effectively. Liao et al. [8] considered that the travel time of vehicles in the network may change, and proposed a two-stage method in which the frequency sweep method was used to allocate vehicles in the first stage and tabu search algorithm was used to improve the route according to real-time information in the second stage. Zhang et al. [9] considered the different travel times of different sections of the distribution network in different time periods, and established a vehicle routing optimization model considering road dynamic congestion with the objective of minimizing the distribution cost and designed an improved genetic algorithm to solve it. Haghani et al. [10] considered the change of travel time between two nodes and the demand of new customers, and studied two cases of whether to consider new demand points. He established a route optimization model aiming at minimizing delivery time, and designed a genetic algorithm to solve it. Duan et al. [11] considered the time-varying and stochastic traffic conditions of the road network at the same time. He established a stochastic time-varying vehicle routing optimization model considering hard time windows based on the principle of minimum and maximum, and solve it by a nondominated sorting ant colony algorithm. Cai et al. [12] considered the time-varying road network environment and the correlation between fuel consumption and load capacity, and established an optimization model of vehicle routing problem under the time-varying road network. He adopted an adaptive ant colony algorithm to solve it. Ge et al. [13] researched the route optimization problem of real-time traffic information, established a mathematical model of open pollution route problem with load and working time constraints, and designed an improved adaptive genetic algorithm to solve it. Wu et al. [14] considered the time-varying characteristics of road network traffic in the actual distribution process, established an optimization model of perishable food integrated production and distribution problem with the lowest total cost, and designed a hybrid genetic algorithm to solve it. Taniguchi et al. [15] researched the advanced intelligent transportation system used in urban distribution, and established a vehicle routing problem solving model considering real-time changes in vehicle travel time, and dynamically updated vehicle travel time using traffic simulation data. Mancini et al. [16] took repetitive congestion into consideration, and proposed to use multiple functional relations to express the speed change of roads in one day, and took Torino as an example for analysis. Liu et al. [17] considered the impact of vehicle speed on carbon emissions, and proposed a cross-time domain calculation method and a method to avoid traffic congestion during rush hours, and designed an improved ant colony algorithm to solve the model.
For MDVRP, Karakatič et al. [18] summarized various genetic algorithms that are designed for solving MDVRP, and evaluated the efficiency of different existing genetic algorithms, and compared the solutions based on genetic algorithms with other existing algorithms. Bezerra et al. [19] proposed a general variable neighborhood search algorithm to solve the MDVRP model aiming at minimizing the total cost. Oliveira et al. [20] established an optimization model with the goal of minimizing the total distribution cost, and proposed a decomposition method for MDVRP, transforming the problem into a single depot VRP, and designed a cooperative coevolutionary algorithm to solve it. Ray et al. [21] established a new optimization model for the MDVRP including depot selection and shared commodity delivery, and proposed an efficient heuristic algorithm to solve it. Afshar-Nadjafi et al. [22] considered the influence of departure time on travel time between nodes, and established a mixed integer programming model with the goal of minimizing total costs, and designed a heuristic algorithm to solve it. Soto et al. [23] established an optimization model aiming at minimizing travel distance for the multi-depot open vehicle routing problem, and proposed a general multiple neighborhood search hybridized with a tabu search strategy to solve the problem. Liu et al. [24] established a mixed integer programming model aiming at minimizing the cost for the MDVRP with time windows based on vehicle leasing and sharing, transforming the problem into a single depot vehicle routing problem by introducing a virtual distribution depot, and proposed a hybrid genetic algorithm. Li et al. [25] considered the constraints of time window, vehicle capacity, and travel time establishing an integer programming model with the minimum total travel cost, and proposed a hybrid genetic algorithm with adaptive local search to solve it. Fan et al. [6] proposed a half-open MDVRP based on joint distribution mode of fresh food, considering the timeliness requirements of fresh products transportation, and constructed an optimization model aiming at minimizing the total distribution cost, and designed an ant colony algorithm to solve the optimization model. Aiming at solving the MDVRP with vehicle capacity and route length constraints, Contardo et al. [26] established an optimization model with the objective of minimizing the total travel time and proposed a new exact algorithm to solve it. Salhi et al. [27] established an optimization model aiming at minimizing the total cost and designed a variable neighborhood search implementation to solve the MDVRP with heterogeneous vehicle fleet.

Problem Description
The TDMDHVRP-TSD studied in this paper can be described as: suppose that there is a complete directed graph G = (V, E), there are different types of roads, and the driving speed of each type of road v = {v 1 , v 2 , · · · , v l } is continuously changing; V = D ∪ C is the all nodes set, D = {1, 2, · · · , m} is the distribution depot set, C = {m + 1, m + 2, · · · , m + n} is the customer points set, and the working time window of each distribution depot is [T s , T f ]; E = {(i, j)|i, j ∈ V } is edge set, l ij is the distance between nodes i and j; the distribution depot has R vehicle types, r(r ∈ R) is the vehicle type, K r is the r type of vehicles set, K r is the r type of k vehicle, its capacity is Q r , vehicle fixed cost is c r1 , unit distance transportation cost is c r2 ; the demand of customer i is d i , T ik r is the time when the vehicle k r arrives at the node i, t ik r is the service time of the vehicle k r at the node i; [ET i , LT i ] is the customer's service time window, ET i is the earliest acceptable service time for customer i, LT i is the latest acceptable service time for customer i, the vehicle arrives earlier than ET i or later than LT i will incur penalty cost, unit time waiting cost is c 3 , unit time delay cost is c 4 . The decision variable x ijk r represents the vehicle k r arrives from point i to point j; the decision variable y ij indicates the customer j is served by the distribution depot i.

Determination of Time-Dependent Function of Vehicle Speed
Most of the existing researches on TDVRP use piecewise function to express the vehicle speed, and the speed is abrupt at a certain point in a time. However, in reality, the vehicle speed varies continuously, such as the trigonometric relation v(t) = ϕ sin(γt) + δ between the vehicle speed (v) and the time (t) proposed in [28] (which ϕ, γ, δ are related to the road conditions). In this paper, the continuous variation of the road speed in a day is approximately expressed by a plurality of trigonometric relations. The trigonometric expression between speed (v) and the time (t) can be expressed as follows: The parameters a β , b β , c β , d β , β ∈ {1, 2, . . . , n} are related to road conditions. According to the data in [29], the change trend of vehicle speed throughout the day can be obtained, as shown in Figure 1. According to the change trend of vehicle speed, the whole day is divided into multiple time periods, and the functional relationship between vehicle speed (v) and the time (t) is different in each time period. Assuming that the time T i when the vehicle leaves node i is within T β , T β+1 , there are two possibilities for the vehicle to travel from node i to node j, i.e., cross the time period and not cross. If l ij ≤ T β+1 T i v(t)dt, the vehicle arrives at the node j within T i , T β+1 without across the time period, the traveling time t ij can be obtained by calculating the upper limit of integration according to the speed function relation of the time period; if l ij > T β+1 T i v(t)dt, the vehicle needs to across time periods, assuming that the vehicle travels from node i to node j, the travel time is t ij = T β+M−1 − T i + t

Build Mathematical Model
Based on the above, the TDMDHVRP-TSD model established in this paper is as follows: Constraint (2) represents the objective function, which is to minimize the sum of vehicle fixed cost, vehicle transportation cost and time window penalty cost. Constraint (3) represents that the customer is only served by a vehicle sent by a distribution depot once, and it is constrained by the balance of entry and exit. Constraint (4) indicates that the number of vehicles starting from any distribution depot is equal to the number of vehicles returning to the distribution depot; Constraint (5) represents the capacity constraint of the vehicle; Constraint (6) indicates that the number of vehicles used by the distribution depot does not exceed the limit of the number of vehicles; Constraint (7) represents the time when the vehicle arrives i from the customer j; Constraint (8) indicates that there is no access between different distribution depots; Constraint (9) is to eliminate subloop constraints; Constraint (10) ensures that the return time of any vehicle to the distribution depot does not exceed the closing time of the distribution depot; Constraint (11) indicates that each customer can only be served by one distribution depot; Constraints (12) and (13) are attributes of decision variables.

Solution Methods
VRP is a classical NP-hard problem. The heuristic algorithm has obvious advantages in solving such a problem. Genetic algorithm has the characteristics of robustness, high parallelism and strong searching ability, but its convergence speed is slow and it is easy to fall into local optimization, which cannot guarantee the overall optimization. The variable neighborhood search algorithm uses multiple different neighborhood structures for systematic search and has strong local search capability. Therefore, this paper combines genetic algorithm and variable neighborhood search algorithm, and on this basis, firstly clusters customers according to the temporal-spatial distance between depots and customers, thus generating a better initial population, and secondly designs a hybrid genetic algorithm with variable neighborhood search considering the temporal-spatial distance (HGAVNS). The variable neighborhood search algorithm is applied to the local search strategy of the genetic algorithm to enhance the optimization ability of the algorithm. The algorithm designed in this paper combines the advantages of genetic algorithm and variable neighborhood search algorithm, and introduces the concept of temporal-spatial distance according to the customer's time window requirements. The proposed algorithm accords with the characteristics of the problem and is very suitable for solving the problem. The algorithm flow is shown in Figure 2.

Customer Clustering and Initial Population Generation
Customers with similar temporal-spatial distance are clustered to generate clustering clusters. Customers within the same cluster have close temporal-spatial distance, while customers between different clusters have far temporal-spatial distance. The clustering method is shown in constraint (14), which aims to minimize the sum of the temporalspatial distance from other customers in each cluster to the center of the cluster, where u is the number of clusters and l ST ij is the temporal-spatial distance between customer i and customer j, and the calculation method reference [30].
Finally, the customers in the cluster are sorted according to their distance from the depot, and the appropriate vehicle types are flexibly selected for route division. The route can be vividly depicted in the three-dimensional map, as shown in Figure 3. The ordinate represents time, the planar two-dimensional coordinate represents spatial position, 0 represents depot, A, B, C, D, E represent customers, and the cylinder represents time window. The temporal-spatial route 1 is 0 3 5 , the time when the vehicle arrives at the customer A and E is earlier than the earliest expected time, resulting in a certain waiting cost, while the time when the vehicle arrives at the customer D is later than the latest expected time, resulting in a corresponding delay cost.

Encoding and Decoding
In this paper, the form of integer coding is adopted. For example, if there are 9 customers and 1 depot, the numbers 1-9 represent customers, and the customers are randomly arranged, the number 0 represents depot. When decoding, customers are divided into vehicles according to the initial arrangement sequence according to the time when the vehicles return to the depot and the load constraints. When the next customer is inspected and found that the current vehicle cannot meet the requirements, the customer can flexibly select other vehicle to serve the customer. It inserts 0 before the first customer, that is to say, the vehicle comes from the depot. This paper will illustrate with a simple example. As shown in Figure 4, assuming that the full order of customers is 1, 5, 8, 3, 4, 6, 9, 2, 7, the first vehicle is sent to serve customer 1. If the constraint is not met when customer 4 is served, then the first route is 1 − 5 − 8 − 3 and vehicle type 1 is selected. By analogy, the second route is 4 − 6 vehicle type 2 is selected; the third route is 9 − 2 − 7, and vehicle type 1 is selected.

Fitness Evaluation
The fitness function of chromosomes can be constructed according to Equation (2) in the model. The fitness function of chromosome K can be expressed as follows: The objective function of this paper is to minimize the total cost, and the objective of optimization is to select chromosomes with larger fitness function values. Where z K is the objective function value of the chromosome K.

Selection
The selection operation adopts a combination of roulette and elite reservation. By roulette, the probability of each chromosome being selected is proportional to the fitness function value, i.e., the higher the fitness function value, the higher probability of the chromosome being selected and, on the contrary, the lower the probability of being selected. After the selection operation is completed, the elite retention strategy is adopted to retain the optimal chromosome of each generation, i.e., the individual with the highest fitness value of the previous generation replaces the individual with the lowest fitness value of the offspring. The strategy of combining roulette and elite reservation is adopted to ensure that the population size remains unchanged and the algorithm converges quickly.

Evolution
The evolutionary operation in this paper selects sequential crossover operator. As shown in Figure 5, when the parent Pop1 is crossed sequentially, the parent Pop2 is randomly selected from the population. Firstly, the nodes i 11 i 12 i 21 "and i 22 are randomly generated. The part between the random nodes i 11 and i 12 of the parent Pop1 is taken as the first segment of the child newPop1. The subsequent nodes of the child newPop1 are related to the parent Pop2, the customer nodes between the random points i 11 and i 12 of the parent Pop2 are firstly eliminated. In the elimination process, the position order of the customer nodes in the parent Pop2 is not changed, and the eliminated customer nodes are arranged as the second segment of the child newPop1 to form a new child newPop1, and the same is true for the child newPop2.
The pseudo code of the crossover operator is as follows:

Local Search Strategy
(1) Neighborhood structure Firstly, neighborhood structure sets N k = {N 1 , N 2 , · · · N l } are constructed, the individual x in the population starts to disturb from the first domain structure N 1 , and if no improved solution is found within the preset neighborhood search times S n , the next neighborhood structure is executed; otherwise, if an improved solution x is obtained in a certain neighborhood structure, then x = x is made, and the iteration is restarted by returning to the first neighborhood structure, until the iteration is cycled to the last neighborhood structure, and when the improved solution is not found, the search is terminated; when the number of variable neighborhood search cycles reaches the preset MaxS n , the search is terminated and the algorithm enters the next stage. In this paper, five neighborhood structures are used to enhance the local search capability of the algorithm: (1) Insert: Randomly select two customer points i and j, insert i after customer j. As shown in Figure 6a, customer 3 and customer 6 are randomly selected and customer 3 is inserted behind customer 6. (2) Exchange: Randomly select two customer nodes i and j to exchange the positions of the two customer nodes. As shown in Figure 6b, the positions of customer 3 and customer 6 are exchanged. (3) 2-insert: In the original scheme, two consecutive customer nodes are randomly selected and inserted after the randomly selected customer point j. As shown in Figure 6c, customers 3 and 4 are inserted behind customer 6. (4) 2-opt: Randomly select two customer nodes i and j, and exchange the order between customer nodes. As shown in Figure 6d, the position of customer 3 is kept unchanged, and customers 4, 5, 7 and 6 are in reverse order. (5) or-opt: In the original scheme, two consecutive customer nodes are randomly selected and inserted into the back of randomly selected customer point j in reverse order. As shown in Figure 6e, customers 3 and 4 are inserted behind customer 6 in reverse order. (2) Adaptive mechanism and new solution acceptance In this paper, an adaptive neighborhood search times strategy and a new solution acceptance mechanism are designed to enhance the breadth and depth of the algorithm so that the variable neighborhood search algorithm can jump out of the local optimization.
Neighborhood search times have great influence on the search ability of the algorithm, which directly leads to the performance of the algorithm. In the iterative process of the algorithm, the disturbance intensity required by the population is different. At the beginning of the iteration, the number of neighborhood searches should be small, to make the population converge quickly. However, with the continuous iteration of the population, the number of neighborhood searches is increased to enhance the search capability of the algorithm. The strategy of adaptive neighborhood search times in this paper is as follows: (1) Setting the initial neighborhood search number S n = 1 and the number of times the optimal solution is continuously unchanged con num; (2) If the optimal solution of the population after this iteration is not improved, let con num = con num + 1, S n = S n + 1; If the perturbed solution is improved, let con num = 0, S n = 1; (3) When the number of times con num that the optimal solution has not changed continuously reaches the preset value stop num, the algorithm terminates and outputs the optimal solution.
In order to further improve the perturbation of the population and expand the search space, this paper uses the acceptance rule of the solution in simulated annealing algorithm to accept the poor solution with a certain probability, thus effectively avoiding the algorithm from falling into local optimization prematurely, improving the optimization ability of the algorithm, and realizing the diversity of the population at the same time. The calculation method of the new solution acceptance probability is shown in Equation (16).
The pseudo code of the HGAVNS is as follows:

16.
Individual P i (t + 1) continues to be disturbed by the next neighborhood structure N l , iter ← 1 ; 17. repeat 18.
Best solution

Algorithm Test
In order to verify the effectiveness of the algorithm, this paper tests the algorithm by selecting instances with different customer sizes. This algorithm is programmed with MATLAB R2018b, the operating system is windows 10, the computer memory is 8G, the CPU is Intel i7-7700M, and the main frequency is 3.60 GHz. After repeated tests, the parameters of this algorithm are set as follows: population size Pop size = 30 ∼ 150, maximum iteration number MAXGEN = 800, initial variable neighborhood search time S n = 1, maximum neighborhood cycle times MaxS n = 1000, stop num = 20 ∼ 50. The setting value of the parameter is related to the customer scale n in the corresponding instances, when n ≤ 50, Pop size = 30, stop num = 20 when 50 <n ≤ 100, Pop size = 100, stop num = 30. When n >100, Pop size= 150, stop num= 50.

The Performance of Proposed Algorithm
To verify the effectiveness of the algorithm in this paper, the standard MDVRP instances (source: http://neo.lcc.uma.es/vrp/vrp-instances/multiple-depot-vrp-instances/, accessed on 15 February 2021) is analyzed. Table 1 shows the results of Greedy randomized adaptive search algorithm (GRASP/VND) [31], general variable neighborhood search algorithm (GVNS) [19], cooperative evolutionary algorithm (CCA) [20] and HGAVNS algorithm. The 'n' column represents the number of customers; the 'd' column represents the number of depots; the 'BKS' column represents the known optimal solution; the 'Best' column represents the best solution found by the algorithm; the '%Dev' column represents the deviation between the best solution found by the algorithm and the known optimal solution. It can be seen from Table 1 that GRASP/VND algorithm has the worst solving quality among the 10 groups of instances, with an average deviation of 5.04% and a maximum deviation of 13.68%. The solution quality of GVNS algorithm is general, with an average deviation of 0.66% and a maximum deviation of 2.52%. The solution quality of CCA algorithm is better, with an average deviation of 0.41% and a maximum deviation of 1.85%. In this paper, the HGAVNS solution quality is the best, the average deviation is 0.34%, and the maximum deviation is 2.09%. Therefore, it can be concluded that the algorithm presented in this paper has a good optimization ability for different scale instances, which verifies the effectiveness and applicability of the algorithm presented in this paper.

Instance Verification
Since there are no instances of TDMDHVRP-TSD in the existing research, this paper randomly generates TDMDHVRP-TSD instances. The calculation example includes 3 distribution depot and 50 customers. The customer's delivery and pickup volume are generated by the demand separation rule proposed by [32]: x i and y i are the coordinate of the customer,q i is the customer's original demand, calculate the ratio r i = min{x i /y i ; y i /x i } of each customer and the delivery volume of the customer is obtained by d i = q i r i . This paper assumes that the working time window of the distribution center is [06 : 00, 19 : 00], the waiting cost per unit time for vehicles is 20, the delay cost per unit time is 30, and the standard time for unloading of per unit goods is 1 min. The cost parameters of different vehicles used by the distribution center for delivery and pickup are shown in Table 2. As shown in Figure 7, the experiment classifies the roads in the distribution network into three types: main roads, secondary roads and branch roads, in which red lines represent main roads with a speed of v 1 , black lines represent secondary roads with a speed of v 2 , and lines not shown in the figure represent branch roads with a speed of v 3 , the all-day variation of the vehicle speed of the three types of roads is shown in Figure 8.  The HGAVNS designed in this paper is used to solve the problem, and the optimal distribution planning is obtained. The specific distribution planning is shown in Table 3.  Table 3 shows that the distribution center sends 5 vehicles to serve the 50 customers, including 1 vehicle with capacity of 60, 3 vehicles with capacity of 100 and 1 vehicle with capacity of 120. The total cost of delivery and pickup is 1599.97 after adding the costs.
In order to verify the validity of the temporal-spatial distance in this paper, this experiment compares the results of the algorithm whether or not considering the temporalspatial distance under the same other conditions. Without considering the temporal-spatial distance, the optimal solution of the algorithm is shown in Table 4. The distribution center dispatched a total of six vehicles for delivery, with a cost of 1828.71. Compared with the optimal scheme above, the total cost increases by 14.30%. Therefore, considering the temporal-spatial distance in the algorithm, the optimal solution is improved, which further proves that the algorithm has strong optimization ability.

Comparison with Different Types of Vehicles
In this section, the different types of vehicles are compared. We compare the vehicles with the payload capacity 60, 100 and 120. Table 5 shows the results that use different types of vehicles. From Table 5, we can see that the cost that use vehicles with payload capacity 120 is highest, as the fixed cost and operational cost of vehicles with payload capacity 120 is highest. The cost that using vehicles with payload capacity 60 is 1950.23, due to the limit of payload capacity, the number of vehicles in use has increased. Use heterogeneous vehicles can reduce the cost. The cost has been reduced by as much as 21.40% and by as little as 16.04%.  1599.97

Conclusions
In this paper, the following conclusions can be drawn from the research on the timedependent multi-depot heterogeneous vehicle routing problem considering temporalspatial distance.
(1) In order to minimize the total cost, an optimization model is established that comprehensively considers customers' service time and spatial location, which can more objectively and accurately reflect the actual operation of the logistics system. The obtained vehicle scheduling optimization scheme can effectively reduce the operation cost of distribution depots.
(2) The HGAVNS algorithm considering the temporal-spatial distance is designed. Clustering customers according to the temporal-spatial distance between customers avoids the situation that the spatial distance of customers is small (large) but the time distance is large (small), and improves the quality of the initial solution.
(3) Regarding selecting operation, the algorithm adopts the strategy of combining elite reservation and roulette to ensure the effective convergence of the algorithm. Moreover, evolutionary operation, adaptive search and new solution acceptance mechanism are used to improve the local search capability and the solution quality, and the obtained results have good stability. The effectiveness of the algorithm is verified by testing instances of different scales and comparing with other literatures.
(4) The TDMDHVRP-TSD model established can solve the multi-depot heterogeneous vehicle routing problem (MDHVRP) with constantly changing road speed and various road types in reality. This is a deepening and expansion of the research on MDHVRP.