Solving the Multi-Depot Green Vehicle Routing Problem by a Hybrid Evolutionary Algorithm

The growing concerns about human pollution has motivated practitioners and researchers to focus on the environmental and social impacts of logistics and supply chains. In this paper, we consider the environmental impact of carbon dioxide emission on a vehicle routing problem with multiple depots. We present a hybrid evolutionary algorithm (HEA) to tackle it by combining a variable neighborhood search and an evolutionary algorithm. The proposed hybrid evolutionary algorithm includes several distinct features such as multiple neighborhood operators, a route-based crossover operator, and a distanceand quality-based population updating strategy. The results from our numerical experiments confirm the effectiveness and superiority of the proposed HEA in comparison with the best-performing methods in the literature and the public exact optimization solver CPLEX. Furthermore, an important aspect of the HEA is studied to assess its effect on the performance of the HEA.


Introduction
During the past decades, the acceleration of energy consumption and environmental pollution have become worldwide concerns. In particular, greenhouse gas emissions and the resulting climate change may have caused more and more disasters worldwide. Consequently, many countries and regions are making huge efforts to reduce greenhouse gas emissions [1]. Meanwhile, consumer pressures and environmental regulations have incentivized many companies to incorporate the environmental considerations into their logistics and supply chain management. As a result, many logistics providers have endeavored to integrate environmental management with their logistics. This has increased the complexity in logistics optimization problems partly because of the potential conflicts between economic and ecological concerns [2].
The environmentally sensitive logistics system requires designing sustainable distribution networks with less negative impact on the environment and the ecology. There are many factors related to green transportation, such as alternative fuels, electric vehicles, and environmental protection infrastructure. In addition, an explicit consideration is given to reducing CO 2 levels through better operating plans. Particularly, the uncertainty about the carbon emissions and mitigation design have obtained lots of attention and research in recent years [3,4] since it has become a huge issue in real life. In order to better study the impact of carbon emissions, in this study we consider and address the carbon emissions as an economic factor.
The vehicle routing problem (VRP) and its related variants have become the subject of lots of research in the literature, which is mainly due to the additional constraints and increasing number of targets caused by realistic problems. For the classic vehicle routing problem, the focus of research is based on the economic part of vehicle routing on companies that provide distribution services. Given the increasing global concern over environmental problems, VRP problems have recently begun to include "green" issues such as alternative fuels and pollution. In particular, multiple objectives and more operational constraints have been incorporated in vehicle routing models for green logistics, which lead to more difficult and more complex optimization problems. As a matter of fact, green logistics have been extensively studied in the literature [5,6] focusing on the environmental effects of different distribution strategies, managing waste disposal, reducing the energy consumption, and adopting the delivery service of the unmanned aerial vehicles. In addition, many new models of transport have recently been presented to improve fuel efficiency and reduce adverse environmental impacts, such as the latest survey on green freight. [7,8].
In order to enrich the research and provide effective ideas and competitive methods for solving the VRP variant regarding these sustainable transportation issues, in this study we attempt to consider a VRP variant with a realistic multi-depot scenario by harmonizing the environmental and economic costs and determining the assignment of customers to depots.

Literature Review
The literature on green vehicle routing problems mainly considered three factors: energy consumption, pollution reduction, and waste management. In the following paragraphs, we will review the related papers in terms of the three aspects.
Palmer [9] was the first to study the environmental issues in the VRPs. Unlike the previous studies that investigated the environmental costs based on the total distance or duration, the authors considered other issues such as road terrain, vehicle speed, and traffic congestion to generate a matrix of CO 2 gas emissions. Their computation results showed that, compared with the distance minimization model and the duration minimization model, the proposed CO 2 minimization model reduced the CO 2 gas emissions.
Later on, Kara et al. [10] proposed a mathematical formula called the energy-minimizing vehicle routing problem (EMVRP), whose purpose was to minimize the sum of the product of the load and distance for each arc. A similar approach was presented to reduce fuel consumption or CO 2 emissions [11,12]. Specifically, Demir et al. [11] made numerical comparisons of several freight vehicle emission models in consideration of their performance in field studies. They demonstrated that reducing greenhouse gas emissions in freight required the use of appropriate emission models in the planning process.
The pollution routing problem (PRP) is a class of green vehicle routing problems. Bektas and Laporte [13] first proposed the PRP that aimed to minimize both environmental and operational costs, taking the greenhouse emissions, fuel, and travel times into account. The authors made a computational analysis to achieve the balance between each variation and speed constraints on the energy, distance, and costs. After that, Demir et al. [14] presented an adaptive large neighborhood search (ALNS) for the PRP problem with the time-window constraints. As a further extension, Demir et al. [15] considered a bi-objective PRP, in which both objectives, namely the minimization of fuel consumption and the minimization of driving time, were conflicting and tackled, respectively. They solved the bi-objective PRP by combining the ALNS method with four a posteriori methods. Their experimental results showed the presented hybrid approach can solve PRP instances with up to 100 nodes. Gajanand and Narendran [16] developed a multiple-route VRP to assess alternative routes between benchmarks and addressed the issue to minimize fuel consumption. Tiwari and Chang [17] employed the distance-based method to calculate the CO 2 emissions, where vehicle load was tackled as an important factor. They generated different clusters for each city they traveled to using different trucks and applied a block reorganization method to the GVRP problem, where each cluster denoted a Sustainability 2020, 12, 2127 3 of 19 block. Suzuki [18] considered the minimization of the emissions from vehicles and fuel consumption. For solving it, they proposed a solution method which employed the mechanism of efficient frontier to exploit the promising area in the solution space.
Considering the importance of multiple depots and maximum capacity constraints in practical scenarios, Jabir et al. [19] integrated emissions and economic costs for a capacitated multi-depot GVRP. An integer linear programming model was formulated to solve a lot of small-scale instances by employing the optimization solver Lingo. In addition, they employed an ant colony optimization (ACO) algorithm to tackle small-scale as well as large-scale problems within a reasonable amount of time.

Motivations and Contributions
The key point to design an effective heuristic method relies on an appropriate trade-off between search intensification and diversification. Hybrid evolutionary algorithms have shown great performance for solving many optimization problems, such as the job shop scheduling problem [20], cutwidth minimization problem [21], and course timetabling problem [22]. To the best of our knowledge, there have been no previous papers using evolutionary algorithms to solve the MDGVRP. Additionally, there are few neighborhood operators in the previous studies in the literature. To help tackle this computationally challenging problem in MDGVRP, in this study we propose a hybrid evolutionary algorithm (HEA) to achieve a better trade-off between the exploitation and exploration in the solution space.
Generally, the HEA algorithm [20,23] is a well-known variant algorithm of the evolutionary algorithm that combines the intensification strength of local optimization and the diversification power of the evolutionary algorithm. The proposed HEA can further be employed to solve the single-depot GVRP, since it can be considered as a special issue of multi-depot GVRP when the number of depots equals 1.
The main contributions in our study are summarized as follows: • We first present a hybrid evolutionary algorithm by incorporating a variable neighborhood search method into the framework of an evolutionary algorithm to solve MDGVRP; • We propose several dedicated neighborhood operators extending the previous search operators in ACO method [19] for search intensification, and the route-based crossover operator for search diversification; • We present a distance-and quality-based replacement strategy to update the population; • Extensive experimental results demonstrate the proposed HEA method can obtain a better trade-off in terms of the computational efficiency and solution quality in comparison with the previous ACO method in literature; • These ideas of the proposed combined evolutionary-based framework and the variable neighborhood search suitable for solving MDGVRP are general and can be employed as a reference for other related GVRPS.
The sections of the paper are organized as follows. Section 2 gives the MDGVRP and its mathematic model. Section 3 proposes the hybrid evolutionary algorithm. Section 4 gives the experimental protocols and parameter setting and reports the computational results and performance comparisons with the best-performing heuristics in the literature and the exact optimization solver CPLEX. The effectiveness of the route-based crossover operator, a key component in the proposed algorithm, is evaluated in Section 5. Section 6 shows the advantages and limitations of the proposed HEA method. At last, conclusion comments and future potential research areas are presented in Section 7.

Problem Description and Mathematic Model
In line with [19], we assume that the vehicles and depots were capacitated with maximum limit constraint. Additionally, we considered similar vehicles with the same speed, capacity, and emissions characteristics. There were environmental and economic costs related with the distribution plan. To be specific, there were two components, i.e., the variable cost and fixed cost, in the economic part. The fixed cost is usually associated with one-time expenditures, such as the rent, the insurance expenses, and other maintenance costs of vehicles. The variable cost considers the expenditure on fuel consumption. Additionally, the environmental costs in this study evaluated the monetary values of CO 2 emissions in the route of vehicle.
We followed the mathematical model from Jabir et al. [19]. The symbols and their definitions for MDGVRP are given in Table 1. The mathematical model for MDGVRP is given as follows: w ijv ≤ 0, i = 1, . . . , m + n, j = m + 1, . . . , m + n, v = 1, . . . , V The objective function presented in Equation (1) seeks to minimize the total cost. Constraint (2) ensures the transportation from a depot to another depot is not permitted. Constraint (3) avoids that the source and the destination are the same. Constraint (4) ensures that a vehicle only can serve one customer each time. Constraint (5) ensures that a single link arrives at and departs from the customer routine. Constraint (6) ensures that the number of outbound links from a depot should be equal to the number of inbound links. Constraint (7) ensures that each vehicle should depart from a depot. Constraint (8) ensures that an empty vehicle returns to a depot. Constraint (9) ensures that the lower bound of the vehicle load in which the vehicle should satisfy the requirement of the destination node as the minimum value. Constraint (10) ensures the upper bound of vehicle load, where the difference between the maximum capacity of vehicle and the requirement of customer node j when a vehicle departs from vertex i to vertex j. Constraint (11) ensures the vehicle load flow balance. Constraint (12) ensures that the depot capacity cannot exceed the maximum capacity upper bound. Constraint (13) illustrates the sub-tour elimination constraints. Constraint (14) guarantees the binary integrality.

A Hybrid Evolutionary Algorithm for MDGVRP
The hybrid evolutionary algorithm proposed in our study was a population-based approach that included several important phases: initial population procedure, variable neighborhood search phase, crossover operator, and population updating mechanism. The initial population procedure aimed to generate diversified solutions with good quality. The variable neighborhood search phase was able to effectively obtain the local optima for search intensification. The route-based crossover operator and distance-and-quality population updating strategy can enhance the search diversification of the proposed method. We present the general framework of HEA and its ingredients in the following subsection.

The Main Framework of HEA for MDGVRP
The proposed HEA algorithm follows the general scheme of evolutionary algorithms and consists of three major components: the population initialization phase to produce a random initial population, the variable neighborhood search (VNS) procedure to improve the incumbent solution, and the route-based crossover operator to generate an offspring. The algorithm progresses through a number of iterative cycles involving VNS and crossover operations. A pictorial representation of one cycle of the algorithmic framework is given in Algorithm 1.

Initial Population Phase
The initial population procedure aims to generate diversified solutions with good quality. Specifically, in order to obtain an initial solution S, which consists of several routes, we present a random greedy construction method inspired from the iterated greedy construction procedure in [24] given as follows: • Step 1: Each route of S is initialized by generating a randomly selected depot Rc as the first node and a randomly customer node as the second node.
Step 2: As for each route, a nearest node to the last node of the current route Rc is iteratively selected and inserted into the current route Rc until the current vehicle cannot serve the remaining customers.
The initial construction procedure is iteratively executed p times to generate several random and initial solutions. Additionally, each solution can be further optimized with a dedicated variable neighborhood search phase. The initialization phase generates an initial population consisting of several solutions with relatively good quality.

Variable Neighborhood Search Procedure
The local refinement procedure is a very important component in the method to obtain the local optima. In this study, the proposed HEA employs variable neighborhood search (VNS) as the local refinement method to improve the initial solution and offspring solutions in population. The key ingredients of a VNS-based approach are several neighborhood operators and shaking phase.

Neighborhood Structure
The proposed six neighborhood operators employed in the VNS procedure can be given as follows: Inter-2opt (N 6 ): The operator eliminates two edges in different routes. Then, each route is cut into two parts, respectively. After that, it reverses two middle customer nodes in two routes and relinks each first part of a route with the last part of the other route to obtain two new routes ( Figure 4).                  The proposed variable neighborhood search algorithm, based on these four neighborhood operators ( N 1 − N 4 ), is applied to obtain local optima in Algorithm 2. For each phase, the algorithm starts by initializing a neighborhood generated by a random neighborhood operator. In the local search procedure, we only considered the neighborhood operators that could generate a feasible solution, i.e., one satisfying the maximum capacity of the depot and vehicle. The size and complexity of the neighborhood structure significantly affected the algorithm's performance. The intra-insertion neighborhood operator ( N 1 ) has n k candidate customers to be removed, and n k − 1 possible candidate positions to consider for each route, where n k denotes the number of vehicles in route R k . Hence the size of N 1 neighborhood is K k=1 (n k × (n k − 1)) where k denotes the number of routes. For intra-swap neighborhood operator ( N 2 ), the size of the neighborhood is the same as that of intra-insertion neighborhood operator. The inter-insertion neighborhood operator ( N 3 ) has n candidate customers to be removed, and ( n − n k ) possible candidate positions to consider for each route, where n k denotes the number of vehicles in route R k . Hence the size of N 3 neighborhood is n × (n − n k ). As for inter-swap neighborhood operator ( N 4 ), the size of the neighborhood is the same as that of the inter-insertion neighborhood operator. Therefore, the complexity of the four neighborhoods operators employed in the proposed method is bounded by O(n 2 ). To escape from the local optima trap, we employed a shaking procedure that used the intra-2opt operator ( N 5 ) and the inter-2opt operator ( N 6 ) for search diversification. Since both move operators are able to cause a significant improvement in terms of solution quality, we employed a threshold function to regulate its degree of convergence. Specifically, a randomly chosen neighborhood solution constructed by N 5 ∪ N 6 would replace the current solution S if the neighboring solution matches the threshold (i.e., f (S ) > (1 − λ) × f (S * ) ) (Algorithm 3 line 7).

Route-Based Crossover Operator
For each iteration of HEA, a crossover operator is employed to generate a new offspring by recombining two randomly chosen parent solutions from the population. The key to design an HEA depends on the crossover operator, which should not only produce different solutions but also pass on valuable parts from parent solutions to offspring solutions. The route-based crossover operator was proposed for MDGVRP by composing several routes. Given two parent solutions and K2 denote the number of routes in parent solutions, respectively, the crossover operator is composed of two steps shown as follows:

•
Step 1: Copy one route R i (1≤i≤K1 or K2) based on the iterated greedy strategy from two parent solutions S 1 and S 2 to the offspring solution. Specifically, the route with maximum value of ∆ f (R i )/| R i | is obtained from the parents, where ∆ f (R i ) denotes the incremental objective value after inserting the route, and | R i | represents the number of customer nodes in route R i .

•
Step 2: Remove the customer nodes in route R i from both two parent solutions S 1 and S 2 .
The route-based crossover operator iteratively alternates with these two steps until all the customer nodes and routes are copied to the offspring solution. The crossover operator can not only inherit some route with good solution quality from two parent solutions, but it also modifies other existing routes by deleting duplicated customer nodes and inserting customer nodes. The crossover operator is able to result in an offspring one that is significantly different from the parents.

The Distance-and Quality-Based Population Updating Mechanism
For the evolutionary algorithm, the population updating mechanism is employed to determine if the offspring solution should replace the worst one in the population. In this study, we propose a distance-and quality-based population updating mechanism to obtain a better balance between computational efficiency and solution quality.
The population updating mechanism employed in the proposed HEA is presented in Algorithm 4, where we denote the offspring solution by S 0 , the closest solution to the offspring solution by S c , and the worst one in the population by S w . As shown in Algorithm 4, the population updating procedure can be described as follows: First, if offspring solution S 0 is better than the closest solution S c , and the distance between the solutions S 0 and S c according to the Hamming distance is not more than the distance threshold β × n, then the offspring solution S 0 can be added in the new population P by displacing the closest one, S c (lines 6-7). Second, if the offspring solution S 0 is better than the worst one S w and the distance between the offspring solution S 0 and the worst solution S w is not less than the distance threshold β × n, then the offspring solution S 0 is added in the new population P J by displacing the worst one S w (lines [8][9]. Finally, the new population P is returned after the distanceand quality-based population updating mechanism (line 11).

Experiment Results
To evaluate the algorithm performance of HEA, we conducted extensive experiments in the following paragraph. In addition, we compared the proposed algorithm with the best-performing methods in the literature and the exact optimization solver CPLEX.

Benchmark Instances and Experiment Setting
The proposed HEA was tested with two instance sets generated with a similar method employed in [19] for the MDGVRP.

•
The first instance set-The small-scale instance set with two depots and the number of customer nodes that ranged from 9 to 14 are given as follows: The data for the benchmark instances are designed in accordance with [7,19]. We coded the HEA algorithm in C++ and ran it on a PC with a 2.60GHz Intel Core processor with the Windows 10 operating system. To avoid over-fitting of the parameters, we adopted ten representative instances to select the best setting of parameters of the proposed method.
The parameters of ( p, Ω, Λ, and β) were tuned with Iterated F-race (IFR) [25]. The total time budget for IRACE was set to 300 executions for HEA, with the maximum time limit being 60 seconds for each benchmark instance. The parameter settings given by IFR are reported in the final value column in Table 2.

Lower Bound
As the mathematical model presented in Section 2, which can be suitable for solving the small-scale instances, a lower bound is required to estimate the solution quality obtained by HEA. Thus, we provided a simple but effective lower bound for this problem. Precisely, we assumed each customer node was visited from the closest customer or depot node, and all the customer nodes were traveled in a route by disregarding capacity constraints of vehicle and depot. Then, the lower bound can be calculated as follows: (15) where γ( j) denotes the closest customer or depot node to the current customer node j. Although this lower bound is quite crude, it presents an effective and important reference indicator of solution quality for the proposed HEA approach when the CPLEX solver cannot solve the instances.

Experimental Results on the First Instance Set
Our comparisons of HEA with the best-performing heuristic ACO and the general optimization solver CPLEX are given in Table 3. The first three columns give the instance number, the number of depots m, and the number of customers n, respectively. Additionally, the following three columns show the lower bound found by CPLEX solver CPLEX LB , the upper bound CPLEX UB , and the gap in percentage (Gap (%)) between the CPLEX LB and CPLEX UB , respectively. The experimental results obtained by the heuristics (i.e., ACO and HEA) are presented in the following eight columns, including the best and average results f best and f avg , the average running time T avg , and gap Gap (%) in percentage to lower bound. The row #Avg indicates average value of each measure, and the row #Best demonstrates the number of instances where the associated approach finds the best results among all compared methods. Table 3 reports 30 instances in the first instance set consisting of problems with the number of depots m = 2, where the CPLEX solver can find optimum solutions for 15 out of 30 instances when the number of customer nodes n is less than 12, within the allotted 3600 second time limit. As the scale of the instance or the number of customer nodes increases, the results found by CPLEX deteriorate. In particular, the optimal solutions of the following 15 instances cannot be found by CPLEX. Here the HEA algorithm obtains the best results for all 30 instances, obtaining better results in terms of the minimum total cost than the best reference heuristic ACO (1670.6 vs. 1674.2) with significantly better computing time (1.1 s vs. 2.3 s). More significantly, HEA finds the best solutions (comparing the solutions obtained by ACO and CPLEX) for 29 out of 30 instances, while ACO and CPLEX are able to obtain best solutions only for 21 instances and 20 instances, respectively.
In summary, one can observe that our proposed HEA method is more effective than the available best-performing algorithms, including the exact optimization solver CPLEX for the first instance set.  Table 4 reports the second instance set with a larger number of depots (m = 3, 4, and 5). The fourth column shows the lower bound (LB) calculated by Equation (15). In addition, columns 8 and 12 show the gap between LB and best results obtained by ACO and HEA algorithms, respectively. As shown in Table 4, the HEA algorithm again obtained the best results among the compared algorithms in terms of all the measures, which included the best total cost f best (6490.5 vs. 6034.1), the average total cost f avg (6514.3 vs. 7033.1), the average computing time T avg (15.1 vs. 23 s), and the average gap to the lower bound Gap (%) (7.6 vs. 14.7) for all the second instance set. Based on the computational results presented above, we can clearly observe that the HEA approach was more effective than the best-performing algorithms, including the exact optimization solver CPLEX for the second instance set.

Analysis for the Route-Based Crossover Operator
HEA applies the route-based crossover operator presented in Section 3.4 to generate feasible and promising offspring solutions. To study the impact of this crossover operator, we compared HEA with two variants that used only one-point and two-point crossover operators, respectively. Both operators are classic crossover operators used in scheduling and routing problems [26,27], and they can be described as follows. In the one-point crossover operator, the first step is to randomly select one cutting point. The subsequent customer node on one side of the cutting point is copied from a parent solution to the offspring solution. Then, the remaining customers are copied to the offspring solution by keeping the same order unchanged. In the two-point crossover operator, two cutting points are randomly chosen to cut parents. By reference to the two parent solutions S1 and S2, the offspring is generated by firstly copying the subsequent nodes between two cutting points in S1, and then transferring the remaining nodes in the same order from S2. In contrast to the route-based operator, both two reference operators usually generate infeasible solutions, violating the maximum capacity constraint. Hence, we employed the inter-insertion neighborhood operator (N3) to reschedule the position of the conflicting customer node so the infeasible solution could satisfy the maximum capacity constraint.
The experimental condition of this analysis was the same as that of the previous section. The gap to lower bound found by CPLEX and Equation (15) for each test instance is plotted in Figure 7. In terms of the gap indicator, one can observe the better performance of the HEA algorithm when the route-based crossover operator is used, which highlights its importance in HEA.

Figure 7.
Comparisons between HEA and its variants that use two crossover operators.

Discussion
The main advantages of the proposed hybrid evolutionary algorithm for MDGVRP can be summarized as follows: First, the variable neighborhood search procedure uses several dedicated neighborhood moves to quickly obtain the local optima for search intensification, which extends the neighborhood operators used in the previous ACO in the literature. Second, the route-based crossover operator and a distance-and quality-based population updating mechanism can effectively enhance the search diversification. Third, the experimental results mentioned above show the proposed HEA method can obtain a better balance in terms of computational efficiency and solution quality in comparison with the ACO algorithm. The limitations of the proposed HEA is that our study only presents the effectiveness of the proposed HEA on the MDGVRP, and the effectiveness on other vehicle routing problems needs to be explored.

Conclusions and Future Research
The green vehicle routing problem has emerged as an important and practical problem in green logistics. To help provide effective methods to solve MDGVRP and its related variants, in this study we considered a VRP variant with a realistic multi-depot scenario by harmonizing the environmental and economic costs and determining the assignment of customers to depots. Specifically, we proposed a hybrid evolutionary algorithm (HEA) for MDGVRP. In addition, we demonstrated the effectiveness of the HEA and its key features, which included a variable neighborhood search embedded into several neighborhood operators, a route-based crossover operator utilizing the ideas of iterated greedy strategy, and a distance-and-quality population updating strategy.
Computational results on small-scale and large-scale benchmark instances demonstrate that our proposed HEA performs very favorably compared with the existing heuristic ACO in the literature and the exact solver CPLEX. Particularly, HEA can find better results in terms of solution quality as well as computational speed for both instance sets.
Future work can be extended in the following directions. First, in order to further enhance the search intensification, we can employ the local search breakout to replace the variable neighborhood search. Second, other diversification techniques can be applied for use in these problems such as scatter search or path-relinking. Finally, the effectiveness of the proposed ideas for tackling the MDGVRP suggests that in the future we test the performance for solving other variants of the vehicle routing problems or other similar combinatorial optimization problems in logistics and supply management [28].