# Multiobjective Green Time-Dependent Location-Routing Problem and Algorithms

Authors: Hua-Xin Zhang; Chun-Miao Zhang
Journal: Advances in Operations Research
 (2022)

---

## Abstract

To reduce the logistic cost and carbon emission and improve customer satisfaction, this study proposes a multiobjective green time-dependent location routing problem (MOGTDLRP) model in which the objectives are to minimize the distribution total cost, delivery time, and fuel consumption. This model will be solved by several hyperheuristic algorithms which include the high-level heuristics and the low-level heuristics. There are three acceptance criterions for the solution: improving and equal, all moves and accept all solutions, and dynamic acceptance criteria. Through the case, the performance of the algorithm and the influence of various factors on the solution are analyzed in this study. The experimental results show that the proposed model can effectively reduce logistic costs, carbon emissions, and vehicle travel time.


---

## Body

## 1. Introduction

Urban logistics plays an important role in the economic vitality of the cities. Green logistics is defined as planning and execution of logistics activities in a more environmentally friendly way by considering external factors such as waste, noise, energy usage, and greenhouse gas (GHG) emissions [1, 2]. At the same time, the urban freight transport is known to be responsible for congestion. The location-routing problem (LRP) is one of the most important combinational optimization problems in supply chain management and logistics system planning and conducts joint decision-making regarding the locations of arbitrary types of facilities and the routing of vehicles [3]. The LRP with environmental issues was called the green LRP (GLRP) which aims to reduce fuel consumption and vehicle exhaust emissions [4]. Koç [5] proposed a model which aims to minimize the total costs including the depot costs, vehicle and routing costs, and emissions costs. Yu [6] investigated the location-routing problem with time-dependent demands as an extension of the location-routing problem by considering the time-dependent demand characteristic. Alamatsaz [7] developed a mixed-integer programming model considering time windows for customers and drivers. There are many other researchers [8–11] who investigated the green LRP model with the objective of minimizing the total costs including depot costs, vehicle and routing costs, and emission costs. However, in addition to reducing costs, other factors such as travel time, fuel consumption, and customer satisfaction should be taken seriously. Therefore, multiobjective GLRP (MOGLRP) problems have gradually become a research hotspot since 2019. Toro [12] proposed a multiobjective model for the GLRP to minimize the operational cost and the environmental effect. Rabbani [13] introduced a new variant of the MOGLRP to minimize the total travel distance and the total costs including vehicle fixed cost and CO2 emissions. Tang [14] established a biobjective model to reduce costs and carbon emissions from the perspective of a sustainable supply chain network.For further analysis, the papers on GLRP in Table1 are compared with five factors: (i) year, (ii) the number of objective functions, (iii) time window, (iv) vehicle speed, and (v) algorithm.Table 1
The papers on GLRP.
ReferencesYearObjectiveTime windowSpeedAlgorithmReferencesYearObjectiveTime windowSpeedAlgorithm[5]2016Single√Heuristic[9]2020Single√Heuristic[14]2016MultiHeuristic[17]2020SingleHeuristic[10]2017Multi√Heuristic[23]2020Multi√Heuristic[12]2017Multi√Heuristic[6]2021Single√√Heuristic[13]2017Single√√Exact[7]2021Multi√Heuristic[15]2018Multi√√Heuristic[18]2021MultiHeuristic[16]2018Single√√Heuristic[19]2021MultiHeuristic[4]2019Single√√Heuristic[20]2021Multi√Heuristic[24]2019Single√√Heuristic[21]2021Single√Heuristic[11]2020SingleHeuristic[23]2021SingleHeuristic[8]2020Single√√HeuristicThere are 11 papers about single GLRP and 10 papers about multiobjective GLRP in Table1. In the above models, some do not consider the speed of the vehicle, and the rest consider the speed to be a constant. However, speed is the critical parameter in calculating fuel consumption and emissions [25]. The fuel consumption is accepted which was affected by many factors in which the vehicle speed is the most critical one, and the vehicle emission is proportional to the fuel consumption [26]. In the real road network, the vehicle speed is the time-dependent function because of the traffic flow, weather, accidents, and other factors. The model based on time-dependent road network is more instructive to the real logistics distribution. So, a three-objective model named the multiobjective green time-dependent location routing problem (MOGTDLRP, as shown in Figure 1. Its goal is to reduce the total costs (including open depot costs, vehicle costs, and distribution costs), the vehicle fuel consumption, and the delivery time.Figure 1
The diagram of the MOGTDLRP.Since the LRP is a NP-hard problem, MOGTDLRP is also a NP-hard problem. In the MOGTDLRP, vehicles travel on different types of roads at different speeds at different times, and deliver goods to customers within the time window. Since exact algorithms can only solve small-scale optimization problems in a reasonable amount of time, metaheuristics are selected to solve the large-scale problems. These metaheuristics include genetic algorithm [27–29], the ant colony optimization [30], and A Tabu search [17]. Different from the metaheuristics, several hyperheuristics (HHs) are proposed in this study. The hyperheuristic system, attempting to find the right sequence of heuristics in a given situation rather than trying to solve a problem directly, was defined as a “heuristic selection heuristic” algorithm [32, 33]. In the HH, the domain barrier isolates the high-level heuristic (HLH), including the selection strategy of the operator and the receiving mechanism of the solution and the low-level heuristic (LLH) which contains a series of low-level operators, problem definitions, objective functions, and other information. The role of the selection strategy is to monitor the performance information of the LLH to select a good and suitable operator, and the receiving mechanism determines whether to replace the parent solution according to the quality of the child solution and controls the search direction of the algorithm and the convergence speed. In this study, we design HHs with heuristic algorithms (artificial bee colony algorithm (ABC), ant colony optimization (ACO), and Tabu search algorithm (TS)) as HLH to solve MOGTDLRP. And we design three acceptance criterions (ACs) for the solution, Improving and Equal (IE), All Moves (AM) and Accept all solutions, and Dynamic Acceptance Criteria (DA) and a series of low-level operators that conform.In summary, the main contributions and innovations of this study are as follows:(1)
The effects of multiple objectives, time-dependent network, and fuel consumption are considered, and the MOGTDLRP mathematical model is constructed(2)
Three HLHs and three ACs combine to form nine HH to analyze the performance of the algorithms(3)
We propose a set of benchmarks of the MOGTDLRP and analyze the simulation resultsThe construction of this study is as follows. Section2 is about the MOGTDLRP model. Section 3 is an introduction to the HHs. Section 4 is the results and analysis. Section 5 is the conclusion and outlook.


## 2. The Model of MOGTDLRP

The model of MOGTDLRP with three optimization objectives will be introduced in this section. The MOGTDLRP can be defined on an asymmetric directed graphG = (V, E), where V is composed of the set of the candidate depots D (D=1,2,…,M) and the set of customers C (C=1,2,…,N) and E is a set of edges (E=i,j:i,j∈V,i≠j/i,j:i,j∈D). The demand qi, service time STi, coordinate (xi, yi), and hard time window [ai, ei] of customer iEC are known. The capacity Pj, rent cost FDj, and coordinates (xj, yj) of the candidate depot j (jED) are known. Each vehicle h in the homogeneous fleet has the same capacity Qh and rent cost FVh. The symbol dij is the distance between the two nodes i and j (i≠j∈V). The delivery time of a vehicle depends on the distance, where the speed changes according to the departure time and the arc being traversed. Vht is the travel speed of the vehicle h traversing different arcs. TTij is the travel time between nodes i and j. Kijh is the load vehicle h leaving i and traveling to j. STj is the service time of customer j. ATjh is the time when vehicle h arrives at node j. The decision variable yr equals 1 if depot r is open and equals 0 otherwise. zim equals 1 if customer i is served by depot m and equals 0 otherwise. xijh equals 1 if the vehicle h travels from i to j and equals 0 otherwise. Fijh, which can be calculated by the CMEM model [8], is the fuel consumption of the vehicle h which travels from node i to node j. Table 2 gives the parameters used in the MOGTDLRP model.Table 2
Parameters and values.
SymbolsDescriptionSymbolsDescriptionDSet of the candidate depotsVhtTravel speed of the vehiclehCSet of the customersKijhLoad of vehicleh leaving i and traveling to jqiThe demands of the customeriSTjService time of customerjSTjService time of customerjATjhThe time when vehicleh arrives at node jxi, yiCoordinate of the nodeiTTijTravel time between nodesi and j[ai, ei]Time window of the nodeiyr1 if depotr is open, and 0 otherwisePjCapacity of the depotjFDjCost of the depotjzim1 if customeri is served by depot m, and 0 otherwiseQhCapacity of the vehiclehFVhCost of the vehiclehxijh1 if customeri is served by depot m, and 0 otherwisedijDistance between the two nodesi and jFijhFuel consumption of the vehicleh from node i to node jSome assumptions need to be satisfied: (1) each customer must be served only once, (2) each vehicle must return to the original depot, (3) every vehicle shall not be overloaded, and (4) the load of each depot must not exceed its capacity.Based on the above assumptions and definitions, the model can be represented as follows:(1)minf1=∑j∈DFDjyj+∑i∈V∑j∈V∑h∈HFVhxijh+∑i∈V∑j∈V∑h∈HFijhxijh,(2)minf2=∑i∈G∑j∈G∑k∈K∑m=1MxijkmTTijm+maxei−ATi,0+STi,(3)minf3=∑i∈V∑j∈V∑h∈HFijhxijh.The following constraints must be satisfied:(4)∑i∈V∑h∈Hxijh=1,∀j∈C,(5)∑h∈H∑i∈Vxijh=∑h∈H∑i∈Vxjih,∀j∈C,(6)∑j∈Dzij=1,∀i∈C,(7)xijh+∑k∈H,k≠h∑r∈V,r≠jxjrk≤1,i∈V,j∈V,i≠j,h∈H,(8)∑h∈Hxijh≤zij,∀i∈C,j∈D,(9)∑h∈Hxjih≤zij,∀i∈C,j∈D,(10)∑h∈Hxijh+zik+∑m∈D,m≠kzjm≤2,∀i,j∈C,k∈D,i≠j,(11)∑i∈D∑j∈CKijh=∑i∈C∑j∈Vqixijh,∀h∈H,(12)∑i∈C∑j∈DKijh=0,∀h∈H,(13)∑i∈Cqizik≤Pkyk,∀k∈D,(14)∑i∈V∑h∈HKijh−qjxijh=∑i∈V∑h∈HKjihxjih,∀j∈C,(15)Kijh≤Qhxijh,∀i,j∈V,i≠j,h∈H,(16)Kijh≥qjxijh,∀i∈V,j∈C,h∈H,(17)ATjh=xijh⋅maxei,ATih+STi+TTijh,∀i,j∈V,h∈H,(18)0≤ATjh≤lj,∀j∈C,h∈H,(19)xijh∈0,1,∀i,j∈V,h∈H,(20)yj∈0,1,∀j∈D,(21)zij∈0,1,∀i∈C,∀j∈D.The description and explanation are as follows. Equation (1) is the costs’ objective. Equation (2) is travel time objective. Equation (3) is the third objective which is the total fuel consumption. Constraints (4) and (5) make sure that each customer is served exactly once. Constraints (6) and (7) impose that each customer is served by one depot and one vehicle, respectively. Constraints (8)–(10) forbid illegal routes that the corresponding vehicles do not return to the departure depots. Constraints (11) and (12) make sure that the demand of each customer is satisfied. Constraint (13) guarantees that the load of each selected depot must be less than its capacity. Constraint (14) is the dynamic equilibrium of the load of each vehicle after visiting customer j. Constraints (15) and (16) guarantee that overloading is not allowed. Constraints (17) and (18) are the bounds on the arrival time at each node: a vehicle must arrive at a customer before the closing time windows. Finally, the last three constraints are decision variables.


## 3. Algorithms

Metaheuristics are wide-ranging. However, choosing the best algorithm and configuring the parameters and operator to solve the problem is very difficult [34]. Numerical experiments have shown that it is impossible to develop a single metaheuristics algorithm that is always efficient for a diverse set of optimization problems. Hyperheuristic algorithm (HH) is a super heuristic algorithm, which provides a high-level heuristic (HLH) method. It can solve various combinatorial optimization problems by managing or manipulating a series of low-level heuristics (LLH). HHs can help researchers to reduce the inherent time and effort required to set up a new domain, as it is a difficult task for testers without deep prior domain-specific knowledge [35]. There are two main tasks in the design of HH: one is the high-level selection strategy and the other is the design of problem domain operator. Using the HH to solve the MOGTDGLRP problem can solve the selection problem of the local optimization operators. If iterating all operators, a lot of time will be wasted. The HH can decide whether to enable the operator through a reasonable selection mechanism according to the history of the operator, which not only saves time but also improves the quality of the solution. In the HLH, there are two different target strategies: operator selection strategy and reconciliation receiving mechanism. In addition, there is an information transmitter between the HLH and the LLH, which transmit information, such as the information of the selection operator, the judgment information of the receiving mechanism, the running time and frequency of the operators, and the number of consecutive unimproved times of the current solution. Therefore, the HLH plays a vital role in improving the performance of the HH. In this study, we have designed nine operators and three HLHs including Tabu search algorithm (TS), ant colony optimization (ACO), and artificial bee colony algorithm (ABC). Next, we will introduce HH according to the following steps: (1) coding method and initial population generation; (2) HLHs; (3) acceptance criteria; (4) the LLHs pool of operators, ξ.


### 3.1. Coding Method and Initial Population Generation

In the MOGTDLRP problem domain, a complete solution is a collection of all routesR=r1,r2,…,rn. Each route ri has the same structure: the first and the end node is the opened depot, and the middle part is the customer number. For example, the route {D2-C1-C3-C10-D2} means that a certain vehicle departs from the depot D2 to deliver goods to customers C1, C3, and C10 in order, and then, returns to D2. The initial solution is critical to the quality of the final solution. In the traditional LRP, when generating an initial route, only the capacity constraints of the depots and vehicles need to be satisfied, and the optimization objective is the distribution distance. However, in the MOGTDLRP, it is also necessary to meet the time window of the depots and customers. In order to produce high-quality initial solutions, it is necessary to design more efficient algorithms. There are two classical algorithms for vehicle-routing problem (VRP) with time windows for reference: Solomon’s Insertion algorithm [36] and IMPACT algorithm [37]. However, these two algorithms need to set a lot of parameters. Zhang [8], in her paper, designed one algorithm named the insertion method based on travel time (IMTT) to generate the initial solution. This method has good solution quality and does not need to set many parameters. Therefore, this study used the IMTT to generate the initial solutions. Its steps are as follows:First step: select the opened depots according to the barycentric methodSecond step: randomly select a customer node to generate a new routeThird step: calculate the total travel time of this routeFourth step: calculate the impact values of the remaining customers inserted into the routeFifth step: insert new customers into the route according to the impact values until one or more constraints are not satisfiedFigure2 shows the generation process of the initial solution. Depot0 is the selected depot; Insi represents the insertable position; Ci is the customer in the route.Figure 2
The generation process of the initial solution.


### 3.2. High-Level Heuristic

#### 3.2.1. The Artificial Bee Colony (ABC) Algorithm

The ABC includes three elements: foods, employed bees, and unemployed bees which have the onlooker bees and the scouts make up. The foods are the solutions in the MOGTDLRP. Employed bees are responsible for finding foods and carrying information about the foods. Onlooker bees follow the employed bees with a certain probability to select the foods. Scouts are responsible for searching for foods randomly to enhance the algorithm’s ability to jump out of local optimal solutions. In the initialization phase of the ABC algorithm, SN solutions, that are the location of the food sources, are generated by IMTT. Then, these initial solutions are evaluated to obtain their fitness values.In the ABCHH, the food sources are the bottom operators, and the optimization process is the process of the bee colony searching for foods. The high-level strategy idea of the ABC is as follows.First, 20% of the bees in the colony are set as scout bees whose roles are to find excellent sources of food. Each scout bee calculates 20 times at each operator. The count table score (i) of the operator i will be added by 1 point when the fitness improves. The operators are sorted in the descending order according to their scores when excellent food has been found, all individuals become employed bees to enjoy the food, and individual is updated by the following equation:(22)Xit+1=Xit⋅pi⋅Sa+randS−Sa−Sbest+Sbest,i=1,2,…,SN,a,b=1,…,S,a≠b,pi=SCaSorpi=1,where i represents the ith individual in the population, t is the number of iterations, Sa represents the ath operator in the operator set, SCa is the score ranking of operator Sa, and Pi·Sa means that the Sa operator is selected with the probability of Pi. When the t generation solution Xti is better than the t-1 generation solution, pi = 1; otherwise, it is calculated according to equation. rand (S-Sa-Sbest) means that an operator is randomly selected after the operator set S eliminates the Sa and Sbest operators, and Sbest is the operator with the highest score in the S set.


#### 3.2.2. Related Definitions

Definition 1.
Status minus.
For the problem of real number coding, let the state of beei be xid. According to equation (23), bee i randomly approaches the position xkd of another bee for information transmission:(23)xid′=xid+ϕidxid−xkd,i=1,2,…,SN,d=1,2,…,D,which means an operation bit d is randomly selected, and the element corresponding to the d bit of xid and xkd forms a field [a b], where a is a smaller value and b is a larger value. Then, λ (xid-xkd) is the reverse order of xid from a to b.
For example, letxid=123456, xkd=352164, and random number β = 4; then, xid (4) = 4, xkd (4) = 1, and ab=14. According to the definition, xid−xkd means that the 1–4 bits of xid are arranged in the reverse order to produce a new solution xid′=432156.Definition 2.
State number multiplication.
C × xid represents the number multiplication of states that means taking the first.
⌈C×D⌉ elements of xid, and the symbol ⌈⌉indicates rounding up.
For example,xid=123456 and C = 0.5; then, ⌈C×D⌉=3 and xid′=123.Definition 3.
Status sum.
Status sum are 2-opt switching operations.
For example, letxid=123456 and XO=2314. First, we convert XO to XO = ((2, 3) (1, 4)). Then, xid + XO indicates that the second bit and the third bit of xid exchange, xid1=132456, and then, the first bit and the fourth bit of xid exchange are the result which is xid2=432156. If the XO element is odd, the last bit is omitted.


#### 3.2.3. The Ant Colony Optimization (ACO) Algorithm

The high-level strategy is based on ACO. There arem ants distributed among the bottom operators. Each ant determines the next vertex according to the pheromone in each route, and the quality of the operator is evaluated by visibility function ηi:(24)ηjt=γηjt−1+∑kmIkjTkjt,where Ikj (t) is the improved value (which can be negative) after ant k applies the bottom operator j in generation t, Tkj (t) is the CPU running time occupied by the calculation, and γ is a value between 0 and 1.The equation ofτij (t), the pheromone on each arc, is as follows:(25)τijt=1−ρτijt−1+∑km#ijPkt⋅IPktTPkt,where ρ (0 < ρ < 1) is the evaporation coefficient of pheromone, 1−ρ represents the persistence coefficient of pheromone, Pkt represents the order in which ant k selects the bottom operator, #ijPkt represents the number of times the arc (i, j) appears in the route of ant k, IPkt is the amount of improvement of ant k in the route, and TPkt is the CPU time.The ant colony selects operators based on pheromone and visibility; the equation is as follows:(26)Vijt=αηjt+βτijt.There is a variable PV to improve the selection probability of some operators with poor performance:(27)PVijt=maxVijt,QσVijt,Q=∑h∈Hmax0,Viht+ε10n,where H is the set of bottom operators, n = |H| is the number of the bottom operators, and ε and σ can ensure that the underlying heuristic operator with poor performance still has a small but nonzero probability of being chosen. If arc (i, j) is an illegal arc, PVij(t) is 0. If the performance of all bottom operators is not good, q = 0, and set all PV values to 1 to ensure that all bottom operators are selected with equal probability. The probability that the arc (i, j) is selected can be expressed as(28)probabilityijkt=PVijt∑h∈HPViht.


#### 3.2.4. Tabu Search (TS) Algorithm

The high-level strategy idea of the TS is as follows. Score evaluates the performance of operators. Then, select an operator with a larger score to improve the current solution in the iterative process. Each operator has the same initial score. When the operator improves the current solution, the score of this operator will add 1; otherwise, it will subtract 1. Detailed process is as follows.First, an initial scorerk = 0 is set for each operator k. The score interval of each operator is [rmin, rmax]; rmin and rmax represent the lowest and highest scores, respectively. The operator k that is not in the Tabu list and has the highest score will be selected. If it can improve the objective function, the score of operator k will add 1. Otherwise, the score will subtract 1 and be put into the Tabu list with fixed Tabu length. Operators in the Tabu list are exempted according to the first in first-out principle.


### 3.3. Acceptance Criteria

The acceptance criterion is the process of whether the new solution is obtained after processing by the operator. If the new solution always unconditionally replaces the parent solution, the optimization time will be too long. However, the HH will fall into a local search if it just saves the improved solution and omits the inferior solution. Other acceptance criteria include the accepting improved solutions completely and accepting inferior solutions in a certain proportion.In this study, three acceptance criteria are proposed:(1)
Improving and equal (IE): accept all improved solutions and reject nonimproved solutions.(2)
All moves (AM): accept all solutions.(3)
Dynamic acceptance criteria (DA): if the current optimal solution is not updated too many times, which means that the solution falls into local search, the acceptance probability of inferior solution should be increased to ensure global search. The acceptance probability is calculated as follows:(29)p=ft−1−ftft+ftbest/2+M1M,whereft−1 is the parent solution, ft is the child solution, ftbest is the optimal solution of generation t, M is a constant, and M1 is the number of times that ftbest has not been updated.


### 3.4. Pool of LLH

This section introduces the nine LLH operators including four operators operating within one route, three operators operating between two routes, and two operators operating the selection of the depots:(1)
2-opt inside one route(2)
Or-opt inside one route(3)
Rev inside one route(4)
Move inside one route(5)
Or-opt between two routes(6)
Interchange between two routes(7)
Crossover between two routes(8)
Depot replace(9)
Depot interchange


### 3.5. Pseudocode of the HH

Table3 lists the pseudocode of the HH, and Figure 3 shows the flowchart.Table 3
Pseudocode of the HH.
(1) //Initialization(2) Set the parameters(3) Initialize Solution (Pop)(4) Current Solution = Pop (r)(5) Best Solution = Current Solution(6) Fitness (Best Solution) = Fitness (Current Solution)(7)  //Main loop(8)   whilet < ITER do(9)   //High-level Selection Strategy(10)   operator = Select (ξ)(11)   //low-level heuristics(12)   [Child Solution] = Implement(Current Solution, operator)(13)   //High-level Acceptance Criterion(14)  if Fitness (Child Solution) <Fitness (Current Solution) then(15)   Current Solution = Child Solution(16)  else(17)   Current Solution = Accept or Not (Child Solution, Current Solution)(18)  end if(19)  //Save the global best solution(20)  if Fitness (Best Solution) >Fitness (Child Solution) then(21)   Best Solution = Child Solution(22)  end if(23)  Update related data(23) end whileFigure 3
The flowchart of the HH.


## 3.1. Coding Method and Initial Population Generation

In the MOGTDLRP problem domain, a complete solution is a collection of all routesR=r1,r2,…,rn. Each route ri has the same structure: the first and the end node is the opened depot, and the middle part is the customer number. For example, the route {D2-C1-C3-C10-D2} means that a certain vehicle departs from the depot D2 to deliver goods to customers C1, C3, and C10 in order, and then, returns to D2. The initial solution is critical to the quality of the final solution. In the traditional LRP, when generating an initial route, only the capacity constraints of the depots and vehicles need to be satisfied, and the optimization objective is the distribution distance. However, in the MOGTDLRP, it is also necessary to meet the time window of the depots and customers. In order to produce high-quality initial solutions, it is necessary to design more efficient algorithms. There are two classical algorithms for vehicle-routing problem (VRP) with time windows for reference: Solomon’s Insertion algorithm [36] and IMPACT algorithm [37]. However, these two algorithms need to set a lot of parameters. Zhang [8], in her paper, designed one algorithm named the insertion method based on travel time (IMTT) to generate the initial solution. This method has good solution quality and does not need to set many parameters. Therefore, this study used the IMTT to generate the initial solutions. Its steps are as follows:First step: select the opened depots according to the barycentric methodSecond step: randomly select a customer node to generate a new routeThird step: calculate the total travel time of this routeFourth step: calculate the impact values of the remaining customers inserted into the routeFifth step: insert new customers into the route according to the impact values until one or more constraints are not satisfiedFigure2 shows the generation process of the initial solution. Depot0 is the selected depot; Insi represents the insertable position; Ci is the customer in the route.Figure 2
The generation process of the initial solution.


## 3.2. High-Level Heuristic

### 3.2.1. The Artificial Bee Colony (ABC) Algorithm

The ABC includes three elements: foods, employed bees, and unemployed bees which have the onlooker bees and the scouts make up. The foods are the solutions in the MOGTDLRP. Employed bees are responsible for finding foods and carrying information about the foods. Onlooker bees follow the employed bees with a certain probability to select the foods. Scouts are responsible for searching for foods randomly to enhance the algorithm’s ability to jump out of local optimal solutions. In the initialization phase of the ABC algorithm, SN solutions, that are the location of the food sources, are generated by IMTT. Then, these initial solutions are evaluated to obtain their fitness values.In the ABCHH, the food sources are the bottom operators, and the optimization process is the process of the bee colony searching for foods. The high-level strategy idea of the ABC is as follows.First, 20% of the bees in the colony are set as scout bees whose roles are to find excellent sources of food. Each scout bee calculates 20 times at each operator. The count table score (i) of the operator i will be added by 1 point when the fitness improves. The operators are sorted in the descending order according to their scores when excellent food has been found, all individuals become employed bees to enjoy the food, and individual is updated by the following equation:(22)Xit+1=Xit⋅pi⋅Sa+randS−Sa−Sbest+Sbest,i=1,2,…,SN,a,b=1,…,S,a≠b,pi=SCaSorpi=1,where i represents the ith individual in the population, t is the number of iterations, Sa represents the ath operator in the operator set, SCa is the score ranking of operator Sa, and Pi·Sa means that the Sa operator is selected with the probability of Pi. When the t generation solution Xti is better than the t-1 generation solution, pi = 1; otherwise, it is calculated according to equation. rand (S-Sa-Sbest) means that an operator is randomly selected after the operator set S eliminates the Sa and Sbest operators, and Sbest is the operator with the highest score in the S set.


### 3.2.2. Related Definitions

Definition 1.
Status minus.
For the problem of real number coding, let the state of beei be xid. According to equation (23), bee i randomly approaches the position xkd of another bee for information transmission:(23)xid′=xid+ϕidxid−xkd,i=1,2,…,SN,d=1,2,…,D,which means an operation bit d is randomly selected, and the element corresponding to the d bit of xid and xkd forms a field [a b], where a is a smaller value and b is a larger value. Then, λ (xid-xkd) is the reverse order of xid from a to b.
For example, letxid=123456, xkd=352164, and random number β = 4; then, xid (4) = 4, xkd (4) = 1, and ab=14. According to the definition, xid−xkd means that the 1–4 bits of xid are arranged in the reverse order to produce a new solution xid′=432156.Definition 2.
State number multiplication.
C × xid represents the number multiplication of states that means taking the first.
⌈C×D⌉ elements of xid, and the symbol ⌈⌉indicates rounding up.
For example,xid=123456 and C = 0.5; then, ⌈C×D⌉=3 and xid′=123.Definition 3.
Status sum.
Status sum are 2-opt switching operations.
For example, letxid=123456 and XO=2314. First, we convert XO to XO = ((2, 3) (1, 4)). Then, xid + XO indicates that the second bit and the third bit of xid exchange, xid1=132456, and then, the first bit and the fourth bit of xid exchange are the result which is xid2=432156. If the XO element is odd, the last bit is omitted.


### 3.2.3. The Ant Colony Optimization (ACO) Algorithm

The high-level strategy is based on ACO. There arem ants distributed among the bottom operators. Each ant determines the next vertex according to the pheromone in each route, and the quality of the operator is evaluated by visibility function ηi:(24)ηjt=γηjt−1+∑kmIkjTkjt,where Ikj (t) is the improved value (which can be negative) after ant k applies the bottom operator j in generation t, Tkj (t) is the CPU running time occupied by the calculation, and γ is a value between 0 and 1.The equation ofτij (t), the pheromone on each arc, is as follows:(25)τijt=1−ρτijt−1+∑km#ijPkt⋅IPktTPkt,where ρ (0 < ρ < 1) is the evaporation coefficient of pheromone, 1−ρ represents the persistence coefficient of pheromone, Pkt represents the order in which ant k selects the bottom operator, #ijPkt represents the number of times the arc (i, j) appears in the route of ant k, IPkt is the amount of improvement of ant k in the route, and TPkt is the CPU time.The ant colony selects operators based on pheromone and visibility; the equation is as follows:(26)Vijt=αηjt+βτijt.There is a variable PV to improve the selection probability of some operators with poor performance:(27)PVijt=maxVijt,QσVijt,Q=∑h∈Hmax0,Viht+ε10n,where H is the set of bottom operators, n = |H| is the number of the bottom operators, and ε and σ can ensure that the underlying heuristic operator with poor performance still has a small but nonzero probability of being chosen. If arc (i, j) is an illegal arc, PVij(t) is 0. If the performance of all bottom operators is not good, q = 0, and set all PV values to 1 to ensure that all bottom operators are selected with equal probability. The probability that the arc (i, j) is selected can be expressed as(28)probabilityijkt=PVijt∑h∈HPViht.


### 3.2.4. Tabu Search (TS) Algorithm

The high-level strategy idea of the TS is as follows. Score evaluates the performance of operators. Then, select an operator with a larger score to improve the current solution in the iterative process. Each operator has the same initial score. When the operator improves the current solution, the score of this operator will add 1; otherwise, it will subtract 1. Detailed process is as follows.First, an initial scorerk = 0 is set for each operator k. The score interval of each operator is [rmin, rmax]; rmin and rmax represent the lowest and highest scores, respectively. The operator k that is not in the Tabu list and has the highest score will be selected. If it can improve the objective function, the score of operator k will add 1. Otherwise, the score will subtract 1 and be put into the Tabu list with fixed Tabu length. Operators in the Tabu list are exempted according to the first in first-out principle.


## 3.2.1. The Artificial Bee Colony (ABC) Algorithm

The ABC includes three elements: foods, employed bees, and unemployed bees which have the onlooker bees and the scouts make up. The foods are the solutions in the MOGTDLRP. Employed bees are responsible for finding foods and carrying information about the foods. Onlooker bees follow the employed bees with a certain probability to select the foods. Scouts are responsible for searching for foods randomly to enhance the algorithm’s ability to jump out of local optimal solutions. In the initialization phase of the ABC algorithm, SN solutions, that are the location of the food sources, are generated by IMTT. Then, these initial solutions are evaluated to obtain their fitness values.In the ABCHH, the food sources are the bottom operators, and the optimization process is the process of the bee colony searching for foods. The high-level strategy idea of the ABC is as follows.First, 20% of the bees in the colony are set as scout bees whose roles are to find excellent sources of food. Each scout bee calculates 20 times at each operator. The count table score (i) of the operator i will be added by 1 point when the fitness improves. The operators are sorted in the descending order according to their scores when excellent food has been found, all individuals become employed bees to enjoy the food, and individual is updated by the following equation:(22)Xit+1=Xit⋅pi⋅Sa+randS−Sa−Sbest+Sbest,i=1,2,…,SN,a,b=1,…,S,a≠b,pi=SCaSorpi=1,where i represents the ith individual in the population, t is the number of iterations, Sa represents the ath operator in the operator set, SCa is the score ranking of operator Sa, and Pi·Sa means that the Sa operator is selected with the probability of Pi. When the t generation solution Xti is better than the t-1 generation solution, pi = 1; otherwise, it is calculated according to equation. rand (S-Sa-Sbest) means that an operator is randomly selected after the operator set S eliminates the Sa and Sbest operators, and Sbest is the operator with the highest score in the S set.


## 3.2.2. Related Definitions

Definition 1.
Status minus.
For the problem of real number coding, let the state of beei be xid. According to equation (23), bee i randomly approaches the position xkd of another bee for information transmission:(23)xid′=xid+ϕidxid−xkd,i=1,2,…,SN,d=1,2,…,D,which means an operation bit d is randomly selected, and the element corresponding to the d bit of xid and xkd forms a field [a b], where a is a smaller value and b is a larger value. Then, λ (xid-xkd) is the reverse order of xid from a to b.
For example, letxid=123456, xkd=352164, and random number β = 4; then, xid (4) = 4, xkd (4) = 1, and ab=14. According to the definition, xid−xkd means that the 1–4 bits of xid are arranged in the reverse order to produce a new solution xid′=432156.Definition 2.
State number multiplication.
C × xid represents the number multiplication of states that means taking the first.
⌈C×D⌉ elements of xid, and the symbol ⌈⌉indicates rounding up.
For example,xid=123456 and C = 0.5; then, ⌈C×D⌉=3 and xid′=123.Definition 3.
Status sum.
Status sum are 2-opt switching operations.
For example, letxid=123456 and XO=2314. First, we convert XO to XO = ((2, 3) (1, 4)). Then, xid + XO indicates that the second bit and the third bit of xid exchange, xid1=132456, and then, the first bit and the fourth bit of xid exchange are the result which is xid2=432156. If the XO element is odd, the last bit is omitted.


## 3.2.3. The Ant Colony Optimization (ACO) Algorithm

The high-level strategy is based on ACO. There arem ants distributed among the bottom operators. Each ant determines the next vertex according to the pheromone in each route, and the quality of the operator is evaluated by visibility function ηi:(24)ηjt=γηjt−1+∑kmIkjTkjt,where Ikj (t) is the improved value (which can be negative) after ant k applies the bottom operator j in generation t, Tkj (t) is the CPU running time occupied by the calculation, and γ is a value between 0 and 1.The equation ofτij (t), the pheromone on each arc, is as follows:(25)τijt=1−ρτijt−1+∑km#ijPkt⋅IPktTPkt,where ρ (0 < ρ < 1) is the evaporation coefficient of pheromone, 1−ρ represents the persistence coefficient of pheromone, Pkt represents the order in which ant k selects the bottom operator, #ijPkt represents the number of times the arc (i, j) appears in the route of ant k, IPkt is the amount of improvement of ant k in the route, and TPkt is the CPU time.The ant colony selects operators based on pheromone and visibility; the equation is as follows:(26)Vijt=αηjt+βτijt.There is a variable PV to improve the selection probability of some operators with poor performance:(27)PVijt=maxVijt,QσVijt,Q=∑h∈Hmax0,Viht+ε10n,where H is the set of bottom operators, n = |H| is the number of the bottom operators, and ε and σ can ensure that the underlying heuristic operator with poor performance still has a small but nonzero probability of being chosen. If arc (i, j) is an illegal arc, PVij(t) is 0. If the performance of all bottom operators is not good, q = 0, and set all PV values to 1 to ensure that all bottom operators are selected with equal probability. The probability that the arc (i, j) is selected can be expressed as(28)probabilityijkt=PVijt∑h∈HPViht.


## 3.2.4. Tabu Search (TS) Algorithm

The high-level strategy idea of the TS is as follows. Score evaluates the performance of operators. Then, select an operator with a larger score to improve the current solution in the iterative process. Each operator has the same initial score. When the operator improves the current solution, the score of this operator will add 1; otherwise, it will subtract 1. Detailed process is as follows.First, an initial scorerk = 0 is set for each operator k. The score interval of each operator is [rmin, rmax]; rmin and rmax represent the lowest and highest scores, respectively. The operator k that is not in the Tabu list and has the highest score will be selected. If it can improve the objective function, the score of operator k will add 1. Otherwise, the score will subtract 1 and be put into the Tabu list with fixed Tabu length. Operators in the Tabu list are exempted according to the first in first-out principle.


## 3.3. Acceptance Criteria

The acceptance criterion is the process of whether the new solution is obtained after processing by the operator. If the new solution always unconditionally replaces the parent solution, the optimization time will be too long. However, the HH will fall into a local search if it just saves the improved solution and omits the inferior solution. Other acceptance criteria include the accepting improved solutions completely and accepting inferior solutions in a certain proportion.In this study, three acceptance criteria are proposed:(1)
Improving and equal (IE): accept all improved solutions and reject nonimproved solutions.(2)
All moves (AM): accept all solutions.(3)
Dynamic acceptance criteria (DA): if the current optimal solution is not updated too many times, which means that the solution falls into local search, the acceptance probability of inferior solution should be increased to ensure global search. The acceptance probability is calculated as follows:(29)p=ft−1−ftft+ftbest/2+M1M,whereft−1 is the parent solution, ft is the child solution, ftbest is the optimal solution of generation t, M is a constant, and M1 is the number of times that ftbest has not been updated.


## 3.4. Pool of LLH

This section introduces the nine LLH operators including four operators operating within one route, three operators operating between two routes, and two operators operating the selection of the depots:(1)
2-opt inside one route(2)
Or-opt inside one route(3)
Rev inside one route(4)
Move inside one route(5)
Or-opt between two routes(6)
Interchange between two routes(7)
Crossover between two routes(8)
Depot replace(9)
Depot interchange


## 3.5. Pseudocode of the HH

Table3 lists the pseudocode of the HH, and Figure 3 shows the flowchart.Table 3
Pseudocode of the HH.
(1) //Initialization(2) Set the parameters(3) Initialize Solution (Pop)(4) Current Solution = Pop (r)(5) Best Solution = Current Solution(6) Fitness (Best Solution) = Fitness (Current Solution)(7)  //Main loop(8)   whilet < ITER do(9)   //High-level Selection Strategy(10)   operator = Select (ξ)(11)   //low-level heuristics(12)   [Child Solution] = Implement(Current Solution, operator)(13)   //High-level Acceptance Criterion(14)  if Fitness (Child Solution) <Fitness (Current Solution) then(15)   Current Solution = Child Solution(16)  else(17)   Current Solution = Accept or Not (Child Solution, Current Solution)(18)  end if(19)  //Save the global best solution(20)  if Fitness (Best Solution) >Fitness (Child Solution) then(21)   Best Solution = Child Solution(22)  end if(23)  Update related data(23) end whileFigure 3
The flowchart of the HH.


## 4. Simulation Experiment and Analysis

All programs are coded in MATLAB R2018a and executed on a computer with an Intel (R) Core (TM) i5-5200U CPU @2.20GHZ, 4 GB of RAM, and Windows 7 operating system.


### 4.1. Benchmarks

The experimental data come from [8]. Table 4 shows the detailed depot data, and Table 5 includes the time windows of all depots. In the tables, “Ben” is benchmark’s name, “Para” is the abbreviation of parameter, “cap” is the capacity, “coor” is the coordinate, “costs” is the fixed open costs, and “TW” notes the time window. Table 6 shows the data of the vehicles. The data about customers come from the Solomon benchmarks.Table 4
Data of the depots.
BenParaDepot1Depot2Depot3Depot4Depot5Depot6Depot7Depot8Depot9Depot10C1cap9908009008508409701000910930780coorx40643544291863851137y5013795740829386317costs40000450004200041000480005000038000490004700046000C2cap95080010109709209901030930870890coorx4089135205129609896y509546436910028439742costs90000100000120000950001050009700011500011200099000117000R1cap9607509108207207901000800790890coorx35105246811194788872y35545660245940776649costs20000190002200021000180002300024000170002500024000R2cap102081072079089010707407001100790coorx3592822564100103864y35771782178772519960costs850009400094000890001000009200097000870009900096000RC1cap1050900109085079094097011809001020coorx40154024506279108087y50523937660100951175costs18000190001700021000260002400023000190002400025000RC2cap130012009008001080780109012409001100coorx40862355286872342661y503794100925219618844costs860009100087000990009600010000085000940009300097000Table 5
Time windows of the depots.
BenTWBenTWC1[0.1236]R2[0.1000]C2[0.3390]RC1[0.240]R1[0.230]RC2[0.960]Table 6
The data of the vehicles.
BenC1C2R1R2RC1RC2Capacity20070020010002001000Vehicle800270050025004502500


### 4.2. The Time-Dependent Speed Function

Ioannou [37] first proposed the vehicle travel time-dependent speed function in his study about vehicle-routing problem and calculated the travel time by the length of the road segment and the travel speed, that is, the travel speed is a time-dependent function. This time-dependent speed function can avoid jumps in travel time and ensure the FIFO characteristics of the time-varying network; that is, the vehicle that departs first arrives at the destination first. For the road network, it can be considered to conform to the FIFO characteristic without overtaking.In this study, to simulate urban road conditions and rush hours, we divide the roads into five types and divide the total service time into four equal time periods. Table7 shows the values.Table 7
The time-dependent speed function.
Road typeS1S2S3S411.802.201.602.4021.602.401.802.2031.402.601.003.0041.202.801.402.6051.003.001.202.80The road type is obtained by equation (30), where Grade (i, j) is the type of road between nodes i and j and Mod (a, b) is the remainder of the number a divided by the number b:(30)Moda,b+1.


### 4.3. Parameter Settings

Population size is 100, the algorithm iteration number is 200, ant path lengthLP = 11, weight parameter α = (0.6,0.7,0.8), β = (0.6,0.7,0.8), and γ = (0.6,0.7,0.8), and enhancement coefficient ɛ = 0.001 and σ = 1.001. Tabu search algorithm parameter maximum score is rankingMax = 5, minimum score is rankingMin = 0, and Tabu length is tabuLength = 4.


### 4.4. Comparative Analysis of the Nine Combined HHs

Three high-level algorithms (ABC, ACO, and TS) and three acceptance criterion (AM, IE, and DA) can form nine HHs named ABC + AM, ABC + IE, ABC + DA, ACO + AM, ACO + IE, ACO + DA, TS + AM, TS + IE, and TS + DA These nine HHs solve the benchmarks of the MOGTDLRP. The results can reflect the performance of the algorithm. The relevant data are as follows:Objective function:f1 minimum total cost.Instances: C1 (C101–C109).Comparison parameters: deviation percentageRD obtained as(31)RD=fA−fbestfbest×100,wherefA is the best result obtained by algorithm A and fbest is the average of the results obtained by the nine HHs.The flow of the nine HHs:Step 1: generating the initial population by the IMTT algorithm and calculating the objective functions to get the initial father solutions S11, S21, and S31Step 2: setting relevant algorithm parametersStep 3: selecting the operator of the LLH according to HLH and calculating the objective functions to get the child solutions S12, S22, and S32Step 4: determining whether to accept the current child solutions according to the acceptance criteriaStep 5: updating related parametersStep 6: if the algorithm termination condition is satisfied, the algorithm ends; otherwise, go to step 3The results of C101–C109 obtained by nine HHs are in Table8.Table 8
Results of C101–C109 obtained by 9 HHs.
NumberHHRDScoresC101C102C103C104C105C106C107C108C1091ABC + AM1.111.201.001.500.981.741.151.201.43302ABC + IE1.051.101.040.960.781.120.670.001.09483ABC + DA0.000.000.000.150.00.000.000.000.00804ACO + AM1.951.501.121.430.970.891.201.521.05275ACO + IE1.451.320.800.790.850.001.131.350.76456ACO + DA0.950.750.580.000.180.000.991.130.00697TS + AM1.131.811.280.780.350.931.211.450.79348TS + IE1.121.020.950.340.320.870.851.080.57599TS + DA0.650.750.000.000.200.570.000.890.0067The maximum value of RD in Table8 is 1.95, which indicates that the solutions obtained by these nine algorithms are not very different. We continue to analyze the performance of the algorithm with Boda algorithm. Each HH has 1 to 9 scores according to the RD value. The algorithm with the smallest RD gets 9 points, the algorithm with the second-smallest RD gets 8 points, and so on. Figure 4 shows the score of the nine HHs. The abscissas 1–9 in the figure correspond to 9 HHs. Figure 4 is a pie chart of the number of optimal solutions obtained by the 9 combination algorithms on C1 examples; that is, the number of times that RD is equal to 0.Figure 4
The bar chart of HHs scored on C1.Through the data of Table8 and Figures 4 and 5, it is known that the ABC + DA combination algorithm scores the highest with 80 points, followed by ACO + DA with an equal score of 69, and the third place is TS + DA with a score of 67. The pie chart shows that the number of times ABC + DA algorithm obtains the optimal solution accounts for 47.06%; ACO + DA accounted for 17.65%; TS + DA accounted for 23.53%; ABC + IE and ACO + IE accounted for 5.88%, respectively, and other combinations do not get the optimal solution. Under ABC high-level strategy, the total score is 30 + 48 + 80 = 158; ACO senior strategy score is 27 + 45 + 69 = 141 points; TS high-level strategy score is 34 + 59 + 67 = 160 points. The score of the AM acceptance criteria is 30 + 27 + 34 + 91, IE is 48 + 45 + 59 = 152, and DA is 80 + 69 + 67 = 16 points. The scores of the three strategies have little difference, but the difference of the solution acceptance criteria is large. AE accepts all solutions, which is easy to cause too divergent search range and difficult to find the optimal solution. IE only accepts noninferior solutions, which is easy to fall into the trap of local optimization. DA accepts all noninferior solutions and conditionally selects inferior solutions. This method gets a higher score because of which it takes into account the advantages of AE and IE and avoid their shortcomings. In general, the ABC + DA is the best combined HH algorithm in the C1 case.Figure 5
The pie chart of HHs scored on C1.


### 4.5. Analysis of the Single Objective Optimization Problem

MOGTDLTP is a multiobjective model, which includes three objectives: total cost, travel time, and fuel consumption. Different optimization objectives lead to different distribution routes. In this section, the difference of the solutions of the different single objective optimization problem is discussed.Instances: RC2 (RC201–RC208).Algorithm: ABC-DASingle objective:f1, minimum total cost, f2, minimum delivery time, and f3, minimum fuel consumption.The results are in Table9.Table 9
The result of the three optimal objectives.
BenchmarksOBJTCRD1TIRD2FCRD3RC201f11.8011e + 0502458.95.291611.84.06f21.8925e + 055.072335.30175113.05f31.8905e + 054.962428.23.981548.90RC202f11.7982e + 0502748.712.091322.65.22f21.8629e + 053.62452.101788.842.30f31.9176e + 056.292730.711.3612570RC203f11.7964e + 0502066.16.351144.18.49f21.8868e + 055.031942.801181.512.01f31.8855e + 054.962025.24.241054.80RC204f11.7936e + 0502648.841.45859.118.08f21.8862e + 055.161872.501121.141.04f31.8829e + 054.982433.641.46794.860RC205f11.798e + 0502667.611.061303.95.47f21.8922e + 055.24240201717.538.92f31.8874e + 054.972684.511.761236.30RC206f11.7971e + 050232118.591214.95.08f21.8881e + 0511.231957.101313.413.6f31.8866e + 0511.052170.610.911156.20RC207f11.7975e + 050203811.591247.67.94f21.8877e + 055.021826.301266.59.58f31.8866e + 054.961952.16.891155.80RC208f11.7948e + 0501678.64.64978.53.53f21.8845e + 0551604.10954.911.03f31.8845e + 054.61795.711.94945.140In Table9, OBJ is the objective, TC is total cost, TI is travel time, and FC is fuel consumption; RDi is the difference percentage calculated by(32)RD1=TCi−TCminTCmin×100,RD2=TIi−TIminTImin×100,RD3=FCi−FCminFCmin×100.The results show that the optimal solutions are different under different optimization objectives. Compared with thef2 and f3 objectives, the total cost obtained with the f1 objective is the minimum, but the delivery time and fuel consumption are not the minimum. The values of RD indicate the deviation percentage between the current value and the minimum value. The larger the value, the greater the deviation. When taking travel time as the optimization objective, the RD values get the maximum seven times, indicating that when taking travel time as the optimization objective, although the total travel time gets the optimization, the total cost and fuel consumption are high. Therefore, it is not appropriate to take travel time as a separate optimization objective. When f1, the total cost is the optimization objective; the average value of RD is 6.62, and when f3, the fuel consumption is the optimization objective; the average value of RD is 5.30. Therefore, the solution obtained under the f3 objective is better.


### 4.6. Comparison between Single Objective and Multiobjective

Instances: RC201 and RC208Optimization objectives:f1, minimum total cost, f2, minimum travel time, f3, minimum fuel consumption, and MO, multiobjective.Algorithm: ABC + DAResults are shown in Table10.Table 10
The result of the four optimal objectives.
InstancesOBJTCTIFCInstancesOBJTCTIFCRC201f11801162458.91611.8RC208f11794851678.6978.5f21892532335.31751f21884561604.1954.91f31890552428.21548.9f31884581795.7945.14MO-TC18015524591655MO-TC1794391663939MO-TI18923323321733MO-TI1884461599946MO-FC18814123501641MO-FC1884151775915In Table10, MO-TC is the solution with the minimum total cost, MO-TI is the solution with the minimum delivery time, and MO-FC is the solution with the minimum fuel consumption in the Pareto solution set.In the single objective function of RC201, the minimum total cost is 180116, the minimum travel time is 2335.3, and the minimum fuel consumption is 1548.9. In the multiobjective function of RC201, the minimum total cost is 180155 > 180116, the deviation isRD = 0.02, the minimum travel time is 2332 < 2335.3, and the deviation is RD = −0.14; the minimum fuel consumption is 1641 > 1548.9, and the deviation is RD = 5.95. In the single objective function of RC208, the minimum total cost is 179485, the smallest travel time is 1604.1, and the smallest fuel consumption is 945.14. In the multiobjective function of RC208, the minimum total cost is 179439 < 179485, the deviation is RD = −0.02, the minimum travel time is 1599 < 1604.1, and the deviation RD = −0.32; the minimum fuel consumption is 915 < 945.14, and the deviation RD = −3.2. Four of the six RD values are less than 0, indicating that the combination of cost, time, and fuel consumption under multiobjective function is better than that under single objective function. Figure 6 is the curves of the combined values of RC201 and RC208, respectively. In the figure, the solid lines are the combined curve obtained by single objective and the dotted lines are the combined curve under multiobjective. The dotted lines are mostly below the solid lines, which also indicate that the combined values obtained by the multiobjective function are better.Figure 6
The curves of the combined values of RC201 and RC208.


### 4.7. Results and Analysis of All Benchmarks

In this section, we use ABC + DA algorithm to solve all benchmarks and analyze the results. The number of iterations is set to 200. All results are in Table11 in which “MC” is the minimum total cost, “MT” is the minimum travel time, and “MF” is the minimum fuel consumption. “AC,” “AF,” and “AT,” respectively, represent the average cost, average fuel consumption, and average time of the whole Pareto solution set. “AV” means the average number of activated vehicles.Table 11
The results of C1∼C2.
BenNPMCMTMFACATAFAVBenNPMCMTMFACATAFAVC101487009955595188980957098010R11224646015999604646615999669C1028898811020615779264710301173611R20111180941224013408747327614273C10369013899461138901449973114411R202218269221101192182704215612043C1047903831003013839039010119139010R203719277419671274192847204813473C10568806096971185919469833128911R20441803401949840180549213710493C10648795496431154879619844116111R205151804782132978180677225811773C10768714196471141881009712115810R206318573119141231185749194612493C10858793695581003895349597106811R20781803751811875180537201410373C10968730796571111898479693118110R2083180357161685718039617478963C2011819637511983575196409120596094R209418583618231219189408187213083C202141964169568616196524114637244R2101318054219391042180700223812003C20311196551119667511968081260310084R2114180283152078318028515207853C204819667610549876196789110579894RC10164206020281660420722048167212C2059196348943554819641898166184RC10264184119081441418511933145112C2069196338937953819637896245784RC10354201519951515420482037154810C207111963619727561196442100096424RC10484267518181175427051850120510C208619640310606603196733110309334RC10574226718971417422971983144713R10174814522761603490212290166415RC10634154318861397420051957155511R10254960321421603497192151171914RC10774080917881309408321792133210R103104833521261335484192167141912RC10844078417791284407861779128610R10474769418201194477321925123211RC201418015523321641185120238417043R10574476219601388463362068162212RC202217983425431334179867256013673R10664823618051236482571819125712RC2031317963019771130179674205111743R10764761117891111476451937114511RC2041218753519241035187964194011313R1086465441736104446582178110829RC2051018879724741297189007253515073R10944568817071188456941709119411RC206817969420961172183089216612143R110104426816821048465191730111910RC207718575418291224186313193112803R11174377216881092464301750121611RC2084179439158591518609015949233The data in Table11 indicate the following:(1)
In general, due to the high quality of solutions, the number of solutions of every benchmark is not a huge amount.(2)
Because of the cost of the depots, the values of MCs are too large which indicate that the ratio of depots’ cost to logistic cost is greater than the ratio of fuel consumption and travel time cost. However, as the number of customers increases, the fuel consumption and travel time cost increase, leading the change in selection depots.(3)
More vehicles need to be enabled to participate in distribution to meet customer time window needs. Taking the RC201 case as an example, the total demand of all customer is 1724 and the capacity of vehicles is 1000. Theoretically, only two vehicles need to be enabled, but in the MOGTDLRP, three vehicles need to be enabled. Due to the use of large capacity vehicles, the number of multiobjective optimal solutions of C2, R2, and RC2 is more than that of corresponding C1, R1, and RC1. The average number of Pareto optimal solutions for C1 is 5.78, while C2 is 10.75, which is about twice that of C1. The average number of Pareto optimal solutions for R1 is 6.42; while R2 is 6.72, which is greater than 6.42. The average number of Pareto optimal solutions for RC1 is 5.75; while RC2 is 7.5, which is greater than 5.75. From a practical point of view, logistic companies want to use fewer vehicles, so when designing the algorithm, the number of vehicles is the priority condition for the solution. On the contrary, using fleets of different capacities is also an effective way for logistic networks to achieve a low-carbon economy.(4)
Figures7(a)–7(c) are the two-dimensional projections of C201 total cost and travel time, travel time and fuel consumption, and total cost and fuel consumption. Figure 7(d) is the three-dimensional diagrams of the solutions of C201. The relationship between cost and fuel consumption is linear because the fuel consumption is part of the total cost. Figures 7(a), 7(b), and 7(d) solution distribution is relatively uniform, indicating that the quality of the solution is higher.Figure 7
The solutions of C201.
(a)(b)(c)(d)


## 4.1. Benchmarks

The experimental data come from [8]. Table 4 shows the detailed depot data, and Table 5 includes the time windows of all depots. In the tables, “Ben” is benchmark’s name, “Para” is the abbreviation of parameter, “cap” is the capacity, “coor” is the coordinate, “costs” is the fixed open costs, and “TW” notes the time window. Table 6 shows the data of the vehicles. The data about customers come from the Solomon benchmarks.Table 4
Data of the depots.
BenParaDepot1Depot2Depot3Depot4Depot5Depot6Depot7Depot8Depot9Depot10C1cap9908009008508409701000910930780coorx40643544291863851137y5013795740829386317costs40000450004200041000480005000038000490004700046000C2cap95080010109709209901030930870890coorx4089135205129609896y509546436910028439742costs90000100000120000950001050009700011500011200099000117000R1cap9607509108207207901000800790890coorx35105246811194788872y35545660245940776649costs20000190002200021000180002300024000170002500024000R2cap102081072079089010707407001100790coorx3592822564100103864y35771782178772519960costs850009400094000890001000009200097000870009900096000RC1cap1050900109085079094097011809001020coorx40154024506279108087y50523937660100951175costs18000190001700021000260002400023000190002400025000RC2cap130012009008001080780109012409001100coorx40862355286872342661y503794100925219618844costs860009100087000990009600010000085000940009300097000Table 5
Time windows of the depots.
BenTWBenTWC1[0.1236]R2[0.1000]C2[0.3390]RC1[0.240]R1[0.230]RC2[0.960]Table 6
The data of the vehicles.
BenC1C2R1R2RC1RC2Capacity20070020010002001000Vehicle800270050025004502500


## 4.2. The Time-Dependent Speed Function

Ioannou [37] first proposed the vehicle travel time-dependent speed function in his study about vehicle-routing problem and calculated the travel time by the length of the road segment and the travel speed, that is, the travel speed is a time-dependent function. This time-dependent speed function can avoid jumps in travel time and ensure the FIFO characteristics of the time-varying network; that is, the vehicle that departs first arrives at the destination first. For the road network, it can be considered to conform to the FIFO characteristic without overtaking.In this study, to simulate urban road conditions and rush hours, we divide the roads into five types and divide the total service time into four equal time periods. Table7 shows the values.Table 7
The time-dependent speed function.
Road typeS1S2S3S411.802.201.602.4021.602.401.802.2031.402.601.003.0041.202.801.402.6051.003.001.202.80The road type is obtained by equation (30), where Grade (i, j) is the type of road between nodes i and j and Mod (a, b) is the remainder of the number a divided by the number b:(30)Moda,b+1.


## 4.3. Parameter Settings

Population size is 100, the algorithm iteration number is 200, ant path lengthLP = 11, weight parameter α = (0.6,0.7,0.8), β = (0.6,0.7,0.8), and γ = (0.6,0.7,0.8), and enhancement coefficient ɛ = 0.001 and σ = 1.001. Tabu search algorithm parameter maximum score is rankingMax = 5, minimum score is rankingMin = 0, and Tabu length is tabuLength = 4.


## 4.4. Comparative Analysis of the Nine Combined HHs

Three high-level algorithms (ABC, ACO, and TS) and three acceptance criterion (AM, IE, and DA) can form nine HHs named ABC + AM, ABC + IE, ABC + DA, ACO + AM, ACO + IE, ACO + DA, TS + AM, TS + IE, and TS + DA These nine HHs solve the benchmarks of the MOGTDLRP. The results can reflect the performance of the algorithm. The relevant data are as follows:Objective function:f1 minimum total cost.Instances: C1 (C101–C109).Comparison parameters: deviation percentageRD obtained as(31)RD=fA−fbestfbest×100,wherefA is the best result obtained by algorithm A and fbest is the average of the results obtained by the nine HHs.The flow of the nine HHs:Step 1: generating the initial population by the IMTT algorithm and calculating the objective functions to get the initial father solutions S11, S21, and S31Step 2: setting relevant algorithm parametersStep 3: selecting the operator of the LLH according to HLH and calculating the objective functions to get the child solutions S12, S22, and S32Step 4: determining whether to accept the current child solutions according to the acceptance criteriaStep 5: updating related parametersStep 6: if the algorithm termination condition is satisfied, the algorithm ends; otherwise, go to step 3The results of C101–C109 obtained by nine HHs are in Table8.Table 8
Results of C101–C109 obtained by 9 HHs.
NumberHHRDScoresC101C102C103C104C105C106C107C108C1091ABC + AM1.111.201.001.500.981.741.151.201.43302ABC + IE1.051.101.040.960.781.120.670.001.09483ABC + DA0.000.000.000.150.00.000.000.000.00804ACO + AM1.951.501.121.430.970.891.201.521.05275ACO + IE1.451.320.800.790.850.001.131.350.76456ACO + DA0.950.750.580.000.180.000.991.130.00697TS + AM1.131.811.280.780.350.931.211.450.79348TS + IE1.121.020.950.340.320.870.851.080.57599TS + DA0.650.750.000.000.200.570.000.890.0067The maximum value of RD in Table8 is 1.95, which indicates that the solutions obtained by these nine algorithms are not very different. We continue to analyze the performance of the algorithm with Boda algorithm. Each HH has 1 to 9 scores according to the RD value. The algorithm with the smallest RD gets 9 points, the algorithm with the second-smallest RD gets 8 points, and so on. Figure 4 shows the score of the nine HHs. The abscissas 1–9 in the figure correspond to 9 HHs. Figure 4 is a pie chart of the number of optimal solutions obtained by the 9 combination algorithms on C1 examples; that is, the number of times that RD is equal to 0.Figure 4
The bar chart of HHs scored on C1.Through the data of Table8 and Figures 4 and 5, it is known that the ABC + DA combination algorithm scores the highest with 80 points, followed by ACO + DA with an equal score of 69, and the third place is TS + DA with a score of 67. The pie chart shows that the number of times ABC + DA algorithm obtains the optimal solution accounts for 47.06%; ACO + DA accounted for 17.65%; TS + DA accounted for 23.53%; ABC + IE and ACO + IE accounted for 5.88%, respectively, and other combinations do not get the optimal solution. Under ABC high-level strategy, the total score is 30 + 48 + 80 = 158; ACO senior strategy score is 27 + 45 + 69 = 141 points; TS high-level strategy score is 34 + 59 + 67 = 160 points. The score of the AM acceptance criteria is 30 + 27 + 34 + 91, IE is 48 + 45 + 59 = 152, and DA is 80 + 69 + 67 = 16 points. The scores of the three strategies have little difference, but the difference of the solution acceptance criteria is large. AE accepts all solutions, which is easy to cause too divergent search range and difficult to find the optimal solution. IE only accepts noninferior solutions, which is easy to fall into the trap of local optimization. DA accepts all noninferior solutions and conditionally selects inferior solutions. This method gets a higher score because of which it takes into account the advantages of AE and IE and avoid their shortcomings. In general, the ABC + DA is the best combined HH algorithm in the C1 case.Figure 5
The pie chart of HHs scored on C1.


## 4.5. Analysis of the Single Objective Optimization Problem

MOGTDLTP is a multiobjective model, which includes three objectives: total cost, travel time, and fuel consumption. Different optimization objectives lead to different distribution routes. In this section, the difference of the solutions of the different single objective optimization problem is discussed.Instances: RC2 (RC201–RC208).Algorithm: ABC-DASingle objective:f1, minimum total cost, f2, minimum delivery time, and f3, minimum fuel consumption.The results are in Table9.Table 9
The result of the three optimal objectives.
BenchmarksOBJTCRD1TIRD2FCRD3RC201f11.8011e + 0502458.95.291611.84.06f21.8925e + 055.072335.30175113.05f31.8905e + 054.962428.23.981548.90RC202f11.7982e + 0502748.712.091322.65.22f21.8629e + 053.62452.101788.842.30f31.9176e + 056.292730.711.3612570RC203f11.7964e + 0502066.16.351144.18.49f21.8868e + 055.031942.801181.512.01f31.8855e + 054.962025.24.241054.80RC204f11.7936e + 0502648.841.45859.118.08f21.8862e + 055.161872.501121.141.04f31.8829e + 054.982433.641.46794.860RC205f11.798e + 0502667.611.061303.95.47f21.8922e + 055.24240201717.538.92f31.8874e + 054.972684.511.761236.30RC206f11.7971e + 050232118.591214.95.08f21.8881e + 0511.231957.101313.413.6f31.8866e + 0511.052170.610.911156.20RC207f11.7975e + 050203811.591247.67.94f21.8877e + 055.021826.301266.59.58f31.8866e + 054.961952.16.891155.80RC208f11.7948e + 0501678.64.64978.53.53f21.8845e + 0551604.10954.911.03f31.8845e + 054.61795.711.94945.140In Table9, OBJ is the objective, TC is total cost, TI is travel time, and FC is fuel consumption; RDi is the difference percentage calculated by(32)RD1=TCi−TCminTCmin×100,RD2=TIi−TIminTImin×100,RD3=FCi−FCminFCmin×100.The results show that the optimal solutions are different under different optimization objectives. Compared with thef2 and f3 objectives, the total cost obtained with the f1 objective is the minimum, but the delivery time and fuel consumption are not the minimum. The values of RD indicate the deviation percentage between the current value and the minimum value. The larger the value, the greater the deviation. When taking travel time as the optimization objective, the RD values get the maximum seven times, indicating that when taking travel time as the optimization objective, although the total travel time gets the optimization, the total cost and fuel consumption are high. Therefore, it is not appropriate to take travel time as a separate optimization objective. When f1, the total cost is the optimization objective; the average value of RD is 6.62, and when f3, the fuel consumption is the optimization objective; the average value of RD is 5.30. Therefore, the solution obtained under the f3 objective is better.


## 4.6. Comparison between Single Objective and Multiobjective

Instances: RC201 and RC208Optimization objectives:f1, minimum total cost, f2, minimum travel time, f3, minimum fuel consumption, and MO, multiobjective.Algorithm: ABC + DAResults are shown in Table10.Table 10
The result of the four optimal objectives.
InstancesOBJTCTIFCInstancesOBJTCTIFCRC201f11801162458.91611.8RC208f11794851678.6978.5f21892532335.31751f21884561604.1954.91f31890552428.21548.9f31884581795.7945.14MO-TC18015524591655MO-TC1794391663939MO-TI18923323321733MO-TI1884461599946MO-FC18814123501641MO-FC1884151775915In Table10, MO-TC is the solution with the minimum total cost, MO-TI is the solution with the minimum delivery time, and MO-FC is the solution with the minimum fuel consumption in the Pareto solution set.In the single objective function of RC201, the minimum total cost is 180116, the minimum travel time is 2335.3, and the minimum fuel consumption is 1548.9. In the multiobjective function of RC201, the minimum total cost is 180155 > 180116, the deviation isRD = 0.02, the minimum travel time is 2332 < 2335.3, and the deviation is RD = −0.14; the minimum fuel consumption is 1641 > 1548.9, and the deviation is RD = 5.95. In the single objective function of RC208, the minimum total cost is 179485, the smallest travel time is 1604.1, and the smallest fuel consumption is 945.14. In the multiobjective function of RC208, the minimum total cost is 179439 < 179485, the deviation is RD = −0.02, the minimum travel time is 1599 < 1604.1, and the deviation RD = −0.32; the minimum fuel consumption is 915 < 945.14, and the deviation RD = −3.2. Four of the six RD values are less than 0, indicating that the combination of cost, time, and fuel consumption under multiobjective function is better than that under single objective function. Figure 6 is the curves of the combined values of RC201 and RC208, respectively. In the figure, the solid lines are the combined curve obtained by single objective and the dotted lines are the combined curve under multiobjective. The dotted lines are mostly below the solid lines, which also indicate that the combined values obtained by the multiobjective function are better.Figure 6
The curves of the combined values of RC201 and RC208.


## 4.7. Results and Analysis of All Benchmarks

In this section, we use ABC + DA algorithm to solve all benchmarks and analyze the results. The number of iterations is set to 200. All results are in Table11 in which “MC” is the minimum total cost, “MT” is the minimum travel time, and “MF” is the minimum fuel consumption. “AC,” “AF,” and “AT,” respectively, represent the average cost, average fuel consumption, and average time of the whole Pareto solution set. “AV” means the average number of activated vehicles.Table 11
The results of C1∼C2.
BenNPMCMTMFACATAFAVBenNPMCMTMFACATAFAVC101487009955595188980957098010R11224646015999604646615999669C1028898811020615779264710301173611R20111180941224013408747327614273C10369013899461138901449973114411R202218269221101192182704215612043C1047903831003013839039010119139010R203719277419671274192847204813473C10568806096971185919469833128911R20441803401949840180549213710493C10648795496431154879619844116111R205151804782132978180677225811773C10768714196471141881009712115810R206318573119141231185749194612493C10858793695581003895349597106811R20781803751811875180537201410373C10968730796571111898479693118110R2083180357161685718039617478963C2011819637511983575196409120596094R209418583618231219189408187213083C202141964169568616196524114637244R2101318054219391042180700223812003C20311196551119667511968081260310084R2114180283152078318028515207853C204819667610549876196789110579894RC10164206020281660420722048167212C2059196348943554819641898166184RC10264184119081441418511933145112C2069196338937953819637896245784RC10354201519951515420482037154810C207111963619727561196442100096424RC10484267518181175427051850120510C208619640310606603196733110309334RC10574226718971417422971983144713R10174814522761603490212290166415RC10634154318861397420051957155511R10254960321421603497192151171914RC10774080917881309408321792133210R103104833521261335484192167141912RC10844078417791284407861779128610R10474769418201194477321925123211RC201418015523321641185120238417043R10574476219601388463362068162212RC202217983425431334179867256013673R10664823618051236482571819125712RC2031317963019771130179674205111743R10764761117891111476451937114511RC2041218753519241035187964194011313R1086465441736104446582178110829RC2051018879724741297189007253515073R10944568817071188456941709119411RC206817969420961172183089216612143R110104426816821048465191730111910RC207718575418291224186313193112803R11174377216881092464301750121611RC2084179439158591518609015949233The data in Table11 indicate the following:(1)
In general, due to the high quality of solutions, the number of solutions of every benchmark is not a huge amount.(2)
Because of the cost of the depots, the values of MCs are too large which indicate that the ratio of depots’ cost to logistic cost is greater than the ratio of fuel consumption and travel time cost. However, as the number of customers increases, the fuel consumption and travel time cost increase, leading the change in selection depots.(3)
More vehicles need to be enabled to participate in distribution to meet customer time window needs. Taking the RC201 case as an example, the total demand of all customer is 1724 and the capacity of vehicles is 1000. Theoretically, only two vehicles need to be enabled, but in the MOGTDLRP, three vehicles need to be enabled. Due to the use of large capacity vehicles, the number of multiobjective optimal solutions of C2, R2, and RC2 is more than that of corresponding C1, R1, and RC1. The average number of Pareto optimal solutions for C1 is 5.78, while C2 is 10.75, which is about twice that of C1. The average number of Pareto optimal solutions for R1 is 6.42; while R2 is 6.72, which is greater than 6.42. The average number of Pareto optimal solutions for RC1 is 5.75; while RC2 is 7.5, which is greater than 5.75. From a practical point of view, logistic companies want to use fewer vehicles, so when designing the algorithm, the number of vehicles is the priority condition for the solution. On the contrary, using fleets of different capacities is also an effective way for logistic networks to achieve a low-carbon economy.(4)
Figures7(a)–7(c) are the two-dimensional projections of C201 total cost and travel time, travel time and fuel consumption, and total cost and fuel consumption. Figure 7(d) is the three-dimensional diagrams of the solutions of C201. The relationship between cost and fuel consumption is linear because the fuel consumption is part of the total cost. Figures 7(a), 7(b), and 7(d) solution distribution is relatively uniform, indicating that the quality of the solution is higher.Figure 7
The solutions of C201.
(a)(b)(c)(d)


## 5. Conclusion and Outlook

In this study, we proposed the MOGTDLRP model with three objectives and a time-dependent speed function with dividing one day into four periods and definition five road types. Fifty six instances were obtained based on the Solomon benchmarks by referring to the relevant literature. In every benchmark, there are 100 customers whose coordinates, requirements, and time windows are known and 10 depots whose coordinates, capacity, and time window are also known. The benchmarks are solved by nine efficient HHs with the high-quality initial solution obtained by IMTT. The results indicated that the MOGTDLRP could effectively reduce the fuel consumption and the total costs.According to the analysis of experimental results, this study can get the following conclusions and suggestions:(1)
The proposed HHs have high accuracy in solving the MOGTDLRP model and can obtain high-quality solution in reasonable computing time. The MOGTDLRP model can effectively reduce logistic costs, fuel consumption, and travel times.(2)
Factors such as customer distribution and customer time windows should be the same issues that logistic companies need to pay attention to as depot costs.(3)
Heterogeneous fleet can reduce the distribution cost.(4)
The traffic speed limit affects the delivery time and vehicle fuel consumption, so the nearest delivery scheme does not necessarily have the least delivery time or fuel consumption.In summary, the time-dependent network is closer to the realistic transportation network than the static network. However, the route from node to node is unique in the time-dependent network. However, there are multiple routes between nodes in real road networks. Furthermore, the actual urban road network is complex, and the speed of vehicles on the same road at the different time periods is random and fluctuating. At the same time, when a vehicle moves from one node to another, it often passes through multiple nodes, and the routes between the nodes are not unique. Therefore, how to solve the optimization problem of green location-path of random road network and multinode road network is the next research direction.



---

*Source: 1811689-2022-05-21.xml*