Next Article in Journal
Effect of Feed Concentrate Intake on the Environmental Impact of Dairy Cows in an Alpine Mountain Region Including Soil Carbon Sequestration and Effect on Biodiversity
Previous Article in Journal
Self-Organizing Map Network-Based Soil and Water Conservation Partitioning for Small Watersheds: Case Study Conducted in Xiaoyang Watershed, China
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
School of Business Administration, Southwestern University of Finance and Economics, Chengdu 610074, China
2
College of Management, Sichuan Agricultural University, Chengdu 610074, China
3
Department of Finance, Wenzhou Business College, Wenzhou 325035, China
*
Author to whom correspondence should be addressed.
Sustainability 2020, 12(5), 2127; https://0-doi-org.brum.beds.ac.uk/10.3390/su12052127
Submission received: 31 January 2020 / Revised: 1 March 2020 / Accepted: 5 March 2020 / Published: 9 March 2020
(This article belongs to the Section Sustainable Transportation)

Abstract

:
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 distance- and 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.

1. 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 CO2 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.

1.1. 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 CO2 gas emissions. Their computation results showed that, compared with the distance minimization model and the duration minimization model, the proposed CO2 minimization model reduced the CO2 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 CO2 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 CO2 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 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.

1.2. 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.

2. 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 CO2 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:
t o t a l   c o s t = v = 1 V i = 1 m + n j = 1 m + n C v a r l i j x i j v + v = 1 V i = 1 m j = m + 1 m + n C f i x x i j v + v = 1 V i = 1 m + n j = 1 m + n l i j P c o 2 w c o 2 V e f f w i j v w p + v = 1 V i = 1 m + n j = 1 m + n x i j v l i j P c o 2 w c o 2 V e f f w c u r b / k
Model constraints:
v = 1 V i = 1 m x i j v = 0   , j = 1 , , m
v = 1 V i = 1 m + n x i j v = 0   , j = 1 , , m + n   and   j = i
v = 1 V i = 1 m + n x i j v = 1   ,   j = m + 1 , , m + n
v = 1 V i = 1 m + n x i j v = v = 1 V i = 1 m + n x j i v   ,   j = m + 1 , , m + n
j = m + 1 m + n x i j v = j = m + 1 m + n x j i v   ,   i = 1 , , m , v = 1 , , V
i = 1 m j = m + 1 m + n x i j v 1   ,   v = 1 , , V
w i j v 0 ,   i = 1 , , m + n ,   j = m + 1 , , m + n ,   v = 1 , , V
v = 1 V d j x i j v v = 1 V w i j v   ,   i = 1 , , m + n ,   j = m + 1 , , m + n
w i j v ( V c a p d i ) x i j v   ,   i = 1 , , m + n ,   j = m + 1 , , m + n ,   v = 1 , , V
i = 1 m + n w i j v i = 1 m + n w j i v = d j   ,   j = m + 1 , , m + n ,   v = 1 , , V
v = 1 V j = m + 1 m + n x i j v D c a p i   ,   i = 1 , , m
i = m + 1 | S | + m j = m + 1 | S | + m x i j v | S | 1   ,   2 | S | V c a p / m i n { d e m a n d } ;   v = 1 , , V
x i j v { 0 , 1 }
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.

3. 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.

3.1. 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.

3.2. 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.
    Algorithm 1. The framework of the proposed HEA for MDGVRP.
    1: Input: Problem instance I , size p of population P
    2: Output: Best found solution S b e s t
    3: S 1 , …, S 2 P o p u l a t i o n _ I n i t i a l i z a t i o n ( I , p )
    4: for i = {1, …, p } do
    5: S i V a r i a b l e _ N e i g h b o r h o o d _ S e a r c h ( S i )
    6: end for
    7: while the maximum running time T m a x t is not met do
    8: Choose two random individuals S j and S k from population P
    9: S 0 R o u t e _ B a s e d _ C r o s s o v e r ( S j ,   S k )
    10:   S 0 V a r i a b l e _ N e i g h b o r h o o d _ S e a r c h ( S 0 )
    11: if f   ( S 0 )   <   f   ( S b e s t ) then
    12: S b e s t S 0   , f ( S b e s t )   f ( S 0 )
    13: end if
    14: P P o p u l a t i o n _ U p d a t i n g ( P , S 0 )
    15: end while
    16: return S b e s t
  • 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.

3.3. 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:
  • Intra-insertion ( N 1 ): The operator chooses a customer in a route and relocates it in the current route (Figure 1);
  • Intra-swap ( N 2 ): The operator exchanges the positions of two customers in the same route (Figure 2);
  • Inter-insertion ( N 3 ): Unlike the intra-insertion operator ( N 1 ), the operator chooses a customer node from a given route and relocates it in another different route (Figure 3);
  • Inter-swap ( N 4 ): Unlike the intra-swap operator ( N 2 ), the operator exchanges the positions of two customers in different routes (Figure 4);
  • Intra-2opt ( N 5 ): The operator deletes two non-adjacent edges, reverses the intermediate customer nodes in the route, and relinks the route by adding two new edges (Figure 5);
  • 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 6).
Algorithm 2. Variable neighborhood search
1: Input: Initial solution S 0 ; maximum iterations without improvement Θ
2: Output: Local optimum S *
3: k 0
4: S * S 0
5: S S 0
6: while Number of move operator does not exceed the maximum number (i.e., 4) of move operators (i.e., k k m a x ) do
7:   S L o c a l S e a r c h ( S , k ) by using the neighborhood operator N k
8: if f   ( S )   <   f   ( S * ) then
9: S * S
10:   k 1
11: else
12:   k k + 1
13: end if
14: S * a r g m i n { f ( S ) , f ( S * ) }
15: S S h a k e ( S , S * )
16: end while
17: return S *
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 = 1 K ( 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).
Algorithm 3. The shaking procedure
1: Input: Current solution S , local optimum S * , the shaking strength Ω , the interval of threshold Λ
2: Output: Perturbed solution S
3: ω 0
4: while ω < do
5: Randomly pick a neighboring solution S N 5 ( S ) N 6 ( S )
6: λ R a n d o m ( Λ )
7: if f   ( S )   > ( 1 - λ ) ×   f   ( S * ) then
8: S S
9: ω   ω + 1
10: end if
11: end while
12: return S

3.4. 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 S 1 = { R 1 1 , R 2 1 ,…, R K 1 1 } and S 2 = { R 1 2 , R 2 2 ,…, R2K2}, where K 1 and K 2 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 Ri (1iK1 or K2) based on the iterated greedy strategy from two parent solutions S1 and S2 to the offspring solution. Specifically, the route with maximum value of ∆ f (Ri)/| Ri | is obtained from the parents, where ∆ f (Ri) denotes the incremental objective value after inserting the route, and | Ri | represents the number of customer nodes in route Ri.
  • Step 2: Remove the customer nodes in route Ri from both two parent solutions S1 and S2.
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.

3.5. 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, Sc (lines 6–7). Second, if the offspring solution S0 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 PJ by displacing the worst one S w (lines 8–9). Finally, the new population P is returned after the distance- and quality-based population updating mechanism (line 11).
Algorithm 4. The distance- and quality-based population updating mechanism
1: Input: Population P , offspring solution S 0 , distance factor parameter β
2: Output: The new population P after population updating strategy
3: S w = a r g m a x { f ( S i ) : i = 1 ,…, p }
4: S c the closest solution to the offspring solution S 0 according to the Hamming distance
5:   D i s the Hamming distance between the offspring solution S 0 and the closest solution S c
6: if f   ( S 0 )   <   f   ( S c ) and   D i s < β × n then
7: P P { S 0 } \ { S c }
8: else if f   ( S 0 )   <   f   ( S w ) and   D i s > β × n then
9: P P { S 0 } \ { S w }
10: end if
11: return P

4. 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.

4.1. 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 number of depots n , The number of customers m } = {(2,9), (2,10), (2,11), (2,12), (2,13),(2,14)};
    The solver CPLEX and the previous best-performing ACO procedure and the proposed HEA algorithm are used for solving the mathematical model proposed in Section 2;
  • The second instance set—The large-scale instance set with the number of depots ranging from 3 to 5 and the number of customer nodes ranging from 20 to 40 are considered as follows:
    {The number of depots n , The number of customers m } = {(3,20), (3,30), (4,20), (4,30), (5,20),(5,40)};
    The proposed HEA algorithm and the previous best-performing ACO procedure are used for solving the mathematical model presented in Section 2.
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.

4.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:
L B = j = m + 1 m + n C v a r l γ ( j ) j + ( n + 1 ) C f i x + j = 1 m + n l γ ( j ) j P C O 2 w C O 2 V e f f d j w p + j = 1 m + n l γ ( j ) j P C O 2 w C O 2 V e f f W c u r b / k
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.

4.3. 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 C P L E X solver C P L E X L B , the upper bound C P L E X U B , and the gap in percentage (Gap (%)) between the C P L E X L B and C P L E X U B , 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 b e s t and f a v g , the average running time T a v g , 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.

4.4. Experimental Results on the Second 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 b e s t (6490.5 vs. 6034.1), the average total cost f a v g (6514.3 vs. 7033.1), the average computing time T a v g (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.

5. 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.

6. 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.

7. 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].

Author Contributions

Conceptualization, B.P. and L.W.; methodology, B.P.; software, B.P.; validation, L.W.,; formal analysis, B.P.; investigation, Y.Y.; resources, Y.Y.; data curation, Y.Y.; writing—original draft preparation, B.P.; writing—review and editing, B.P., L.W. and X.C.; visualization, B.P.; supervision, X.C.; project administration, X.C.; funding acquisition, B.P. and X.C. All authors have read and agree to the published version of the manuscript.

Funding

This research was funded by Fundamental Research Funds for the Central Universities of China (Grant Number: JBK2001013 and JBK190504), Key Projects of the National Social Science Fund (Grant Number: 16AJL004), Natural Science Foundation of Zhejiang (Grant Number: LQ19G020002), the 13th Five-Year Plan Teaching Reform Project of Higher Education of Zhejiang (Grant Number: jg20180428), the Project of Education Department of Zhejiang (Grant Number: Y201940902), Basic Scientific Research Project of Wenzhou(Grant Number: R20180002), and Wenzhou Social Sciences Planning Project (Grant Number: 19wsk225).

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Zhou, Y.; Lee, G. A Lagrangian relaxation-based solution method for a green vehicle routing problem to minimize greenhouse gas emissions. Sustainability 2017, 9, 776. [Google Scholar] [CrossRef] [Green Version]
  2. Ubeda, S.; Arcelus, F.J.; Faulin, J. Green logistics at Eroski: A case study. Int. J. Prod. Econ. 2011, 131, 44–51. [Google Scholar] [CrossRef]
  3. Guo, J.X.; Tan, X.; Gu, B. The impacts of uncertainties on the carbon mitigation design: Perspective from abatement cost and emission rate. J. Clean. Prod. 2019, 232, 213–223. [Google Scholar] [CrossRef]
  4. Nocera, S.; Galati, O.I.; Cavallaro, F. On the uncertainty in the economic valuation of carbon emissions from transport. J. Transp. Econ. Policy 2018, 52, 68–94. [Google Scholar]
  5. Li, J.; Guo, H.; Zhou, Q.; Yang, B. Vehicle Routing and Scheduling Optimization of Ship Steel Distribution Center under Green Shipbuilding Mode. Sustainability 2019, 11, 4248. [Google Scholar] [CrossRef] [Green Version]
  6. Eun, J.; Song, B.D.; Lee, S.; Lim, D.E. Mathematical Investigation on the Sustainability of UAV Logistics. Sustainability 2019, 11, 5932. [Google Scholar] [CrossRef] [Green Version]
  7. Lin, C.; Choy, K.L.; Ho, G.T.; Chung, S.H.; Lam, H. Survey of green vehicle routing problem: Past and future trends. Expert Syst. Appl. 2014, 41, 1118–1138. [Google Scholar] [CrossRef]
  8. Braekers, K.; Ramaekers, K.; Van Nieuwenhuyse, I. The vehicle routing problem: State of the art classification and review. Comput. Ind. Eng. 2016, 99, 300–313. [Google Scholar] [CrossRef]
  9. Palmer, A. The Development of an Integrated Routing and Carbon Dioxide Emissions Model for Goods Vehicles. Available online: https://dspace.lib.cranfield.ac.uk/handle/1826/2547 (accessed on 7 March 2020).
  10. Kara, I.; Kara, B.Y.; Yetis, M.K. Energy minimizing vehicle routing problem. In International Conference on Combinatorial Optimization and Applications; Springer: Berlin/Heidelberg, Germany, 2007. [Google Scholar]
  11. Demir, E.; Bektas, T.; Laporte, G. A comparative analysis of several vehicle emission models for road freight transportation. Transp. Res. Part D Transp. Environ. 2011, 16, 347–357. [Google Scholar] [CrossRef]
  12. Kopfer, H.W.; Schönberger, J.; Kopfer, H. Reducing greenhouse gas emissions of a heterogeneous vehicle fleet. Flex. Serv. Manuf. J. 2014, 26, 221–248. [Google Scholar] [CrossRef]
  13. Bektas, T.; Laporte, G. The pollution-routing problem. Transp. Res. Part B Methodol. 2011, 45, 1232–1250. [Google Scholar] [CrossRef]
  14. Demir, E.; Bektas, T.; Laporte, G. An adaptive large neighborhood search heuristic for the pollution-routing problem. Eur. J. Oper. Res. 2012, 223, 346–359. [Google Scholar] [CrossRef]
  15. Demir, E.; Bektas¸, T.; Laporte, G. The bi-objective pollution-routing problem. Eur. J. Oper. Res. 2014, 232, 464–478. [Google Scholar] [CrossRef]
  16. Gajanand, M.; Narendran, T. Green route planning to reduce the environmental impact of distribution. Int. J. Logist. Res. Appl. 2013, 16, 410–432. [Google Scholar] [CrossRef]
  17. Tiwari, A.; Chang, P.C. A block recombination approach to solve green vehicle routing problem. Int. J. Prod. Econ. 2015, 164, 379–387. [Google Scholar] [CrossRef]
  18. Suzuki, Y. A dual-objective metaheuristic approach to solve practical pollution routing problem. Int. J. Prod. Econ. 2016, 176, 143–153. [Google Scholar] [CrossRef]
  19. Jabir, E.; Panicker, V.V.; Sridharan, R. Design and development of a hybrid ant colony-variable neighbourhood search algorithm for a multi-depot green vehicle routing problem. Transp. Res. Part D Transp. Environ. 2017, 57, 422–457. [Google Scholar] [CrossRef]
  20. Cheng, T.; Peng, B.; Lü, Z. A hybrid evolutionary algorithm to solve the job shop scheduling problem. Ann. Oper. Res. 2016, 242, 223–237. [Google Scholar] [CrossRef]
  21. Bansal, R.; Srivastava, K.; Srivastava, S. A hybrid evolutionary algorithm for the cutwidth minimization problem. In Proceedings of the 2012 IEEE Congress on Evolutionary Computation, Brisbane, Australia, 10–15 June 2012. [Google Scholar]
  22. Abdullah, S.; Burke, E.K.; McCollum, B. A hybrid evolutionary approach to the university course timetabling problem. In Proceedings of the 2007 IEEE Congress on Evolutionary Computation, Singapore, 25–28 September 2007. [Google Scholar]
  23. Wang, Y.; Cai, Z.; Guo, G.; Zhou, Y. Multiobjective optimization and hybrid evolutionary algorithm to solve constrained optimization problems. IEEE Trans. Syst. Man Cybern. Part B 2007, 37, 560–575. [Google Scholar] [CrossRef]
  24. Karabulut, K.; Tasgetiren, M.F. A variable iterated greedy algorithm for the traveling salesman problem with time windows. Inf. Sci. 2014, 279, 383–395. [Google Scholar] [CrossRef]
  25. Birattari, M.; Yuan, Z.; Balaprakash, P.; Stützle, T. F-Race and Iterated F-Race: An Overview. In Experimental Methods for the Analysis of Optimization Algorithms; Springer: Berlin/Heidelberg, Germany, 2010; pp. 311–336. [Google Scholar]
  26. Kellegöz, T.; Toklu, B.; Wilson, J. Comparing efficiencies of genetic crossover operators for one machine total weighted tardiness problem. Appl. Math. Comput. 2008, 199, 590–598. [Google Scholar] [CrossRef]
  27. Lu, Y.; Hao, J.K.; Wu, Q. Hybrid evolutionary search for the traveling repairman problem with profits. Inf. Sci. 2019, 502, 91–108. [Google Scholar] [CrossRef]
  28. Zhou, R.; Liao, Y.; Shen, W. Channel selection and fulfillment service contracts in the presence of asymmetric service information. Int. J. Prod. Econ. 2019, 107504. [Google Scholar] [CrossRef]
Figure 1. Insert customer node 5 between customer nodes 4 and 6 in the same route.
Figure 1. Insert customer node 5 between customer nodes 4 and 6 in the same route.
Sustainability 12 02127 g001
Figure 2. Exchange the positions of two customer nodes 2 and 5 in the same route.
Figure 2. Exchange the positions of two customer nodes 2 and 5 in the same route.
Sustainability 12 02127 g002
Figure 3. Remove customer 7 from the current route and relocate it to a different route.
Figure 3. Remove customer 7 from the current route and relocate it to a different route.
Sustainability 12 02127 g003
Figure 4. Exchange the positions of two different customer nodes 7 and 9 from two different routes.
Figure 4. Exchange the positions of two different customer nodes 7 and 9 from two different routes.
Sustainability 12 02127 g004
Figure 5. Delete two non-adjacent edges (3 4 and 7 D), reverse the intermediate customer nodes (4 5 6 7) to (7 8 5 4), and add two new edges (4 D and 3 7) in the same route.
Figure 5. Delete two non-adjacent edges (3 4 and 7 D), reverse the intermediate customer nodes (4 5 6 7) to (7 8 5 4), and add two new edges (4 D and 3 7) in the same route.
Sustainability 12 02127 g005
Figure 6. Remove two edges (3 4) and (7 D 1 ) in the first route and two edges (8 9) and (10 D 2 ) in the second route. Then, reverse the intermediate customer nodes (9 10) to (10 9) in second route, and add two new edges (3 10) and (9 D 1 ) in the current route. Reverse the intermediate customer nodes (4 5 6 7) to (7 6 5 4) in first route, and add two new edges (4 D 2 ) and (8 7) in the current route.
Figure 6. Remove two edges (3 4) and (7 D 1 ) in the first route and two edges (8 9) and (10 D 2 ) in the second route. Then, reverse the intermediate customer nodes (9 10) to (10 9) in second route, and add two new edges (3 10) and (9 D 1 ) in the current route. Reverse the intermediate customer nodes (4 5 6 7) to (7 6 5 4) in first route, and add two new edges (4 D 2 ) and (8 7) in the current route.
Sustainability 12 02127 g006
Figure 7. Comparisons between HEA and its variants that use two crossover operators.
Figure 7. Comparisons between HEA and its variants that use two crossover operators.
Sustainability 12 02127 g007
Table 1. Symbols and the corresponding definitions for the MDGVRP.
Table 1. Symbols and the corresponding definitions for the MDGVRP.
SymbolsDefinitions
n Maximum number of customer nodes
m Maximum number of depots
V Maximum number of vehicles
x i j v 1, if vehicle v travels from customer node i to customer node j ,
i , j , v Customer index and vehicle index, respectively
I i j Driving distance from customer node i to customer node j
C f i x Fixed vehicle cost
C v a r Variable vehicle cost per unit distance
V e f f Volume of fuel consumption per unit vehicle weight per unit distance
F C O 2 CO2 emissions cost per vehicle weight per unit distance
w p Weight of every transportation product
P C O 2 Average cost per unit weight of CO2
w C O 2 CO2 emission weight per liter fuel consumption w c u r b
w c u r b Average gross weight per vehicle on each route d j
d j Requirement of customer node j
k Ratio of vehicle volume to curb weight
D c a p i Depot i capacity
V c a p Vehicle capacity
Table 2. Settings of the parameters used in HEA.
Table 2. Settings of the parameters used in HEA.
ParametersDescriptionCandidate ValuesFinal Value
p The number of individual solutions in the population{5, 10, 20}5
Ω The threshold of shaking strength in the shaking procedure3, 6, 93
Λ Interval of threshold ratio values in the shaking procedure(0,1)[0.1, 0.2]
β Distance factor parameter in the population updating strategy0.05,0.1,0.20.1
Table 3. The performance of the HEA compared with the ACO algorithm and the general solver CPLEX for the first instance set.
Table 3. The performance of the HEA compared with the ACO algorithm and the general solver CPLEX for the first instance set.
InstancemnCPLEXLBCPLEXUBGap (%)ACOHEA
f b e s t f a v g T a v g Gap % f b e s t f a v g T a v g Gap %
1291246.71246.701246.71246.71.101246.71246.70.10
2291481.61481.601481.61481.61.301481.61481.60.90
3291434.91434.901434.91434.91.501434.91434.90.20
4291222.61222.601222.61222.60.901222.61222.60.60
5291285.71285.701285.71285.72.301285.71285.71.10
62101688.31688.301688.31688.32.801688.31688.31.10
72101553.81553.801553.81553.81.701553.81553.80.80
82101341.81341.801341.81341.82.101341.81341.80.90
92101464.61464.601464.61464.62.801464.61464.60.50
102101416.21416.201416.21416.22.201416.21416.20.90
112111581.81581.801581.81581.81.201581.81581.81.20
122111769.31769.301769.31769.32.901769.31769.31.90
132111792.31792.301792.31792.31.501792.31792.320
142111667.61667.601667.61667.61.701667.61667.60.90
152111807.91807.901807.91807.91.901807.91807.90.60
162121807.71846.22.11807.718921.801807.71849.40.30
1721219952053.72.919952013.12.1019952013.30.50
182121839.61940.25.51839.618611.601839.61881.61.90
19212177918514.117791838.82.3017791780.21.70
202121799.71853.32.91799.71851.13.501799.71838.71.60
212131710.81765.53.21722.81765.31.30.71717.61752.210.4
222131876.91928.72.81884.419763.40.41878.81931.40.70.2
232131809.51916.65.91814.91850.33.90.31816.71837.91.40.5
2421317161802.15.01738.31804.721.31721.11782.21.40.3
252131996.52097.15.02008.52080.32.10.62000.52040.51.10.2
262141902.22101.910.51921.21988.72.31.21907.91943.61.30.4
272141718.41908.611.117391823.33.91.21720.11764.42.50.2
282141759.32035.515.81778.717943.81.11764.618071.10.4
292141824.3208114.11837.11917.32.90.81827.918401.50.3
3021417842099.417.71803.61833.62.81.21785.81812.62.20.2
#Avg 1669.11734.57.21674.21701.52.30.81670.61687.71.10.2
#Best 20 21 29
Table 4. The performance of HEA compared with the ACO algorithm and the general solver CPLEX for the second instance set.
Table 4. The performance of HEA compared with the ACO algorithm and the general solver CPLEX for the second instance set.
InstancemnLBACOHEA
f b e s t f a v g T a v g Gap % f b e s t f a v g T a v g Gap %
313204069.94688.54719.76.315.243674372.35.77.3
323204793.25564.95638616.15267.75269.74.99.9
333203203.63780.23832.25.3183463.13490.42.58.1
343203972.64457.34532.46.912.24302.34302.75.88.3
353203279.23656.33686.1711.53495.63507.62.56.6
363309361.611,215.211,31526.119.810,110.510,1399.68
373307422.48172.18336.924.210.17904.97924.27.56.5
383307160.48470.88622.729.218.37547.17581.8105.4
393307702.98935.4899724.6168326.88336.913.68.1
403308056.49506.69703.627.6188507.68514.918.75.6
414203990.34417.34499.56.810.74201.84221.84.55.3
424204082.44535.54642.9511.14355.94382.83.16.7
434206181.66972.87063.98.112.86645.26650.13.67.5
444205476.26127.96156.8711.95930.75979.65.18.3
454206633.67641.97754.95.715.27151.37157.927.8
464306239.57075.67202.844.513.46564.26597.919.25.2
474306195.46969.87067.940.512.56548.56566.924.25.7
484305792.36881.37026.638.918.86284.66346.723.48.5
494306537.87505.47675.839.814.87047.77099.242.97.8
504306761.877497821.141.614.67370.47376.536.89
515205007.46003.96147.412.419.95503.15522.96.49.9
5252067797680.67696.39.813.37429.87479.14.29.6
535203871.44363.14442.68.412.74243.14274.49.19.6
545204035.146004623.79.2144382.143843.58.6
555203057.93458.53534.112.313.13265.83273.112.56.8
565406774.67865.37931.832.116.17201.47201.527.66.3
5754061406809.36856.15710.96606.66650.433.67.6
5854010,628.111,999.112,258.45012.911,467.71148238.77.9
595409201.710,858.410,926.155.7189910.29987.830.97.7
605408614.410,268.410,280.140.519.29312.29356.142.78.1
#Avg 6034.16805.97033.12314.76490.56514.315.17.6
#Best 0 30

Share and Cite

MDPI and ACS Style

Peng, B.; Wu, L.; Yi, Y.; Chen, X. Solving the Multi-Depot Green Vehicle Routing Problem by a Hybrid Evolutionary Algorithm. Sustainability 2020, 12, 2127. https://0-doi-org.brum.beds.ac.uk/10.3390/su12052127

AMA Style

Peng B, Wu L, Yi Y, Chen X. Solving the Multi-Depot Green Vehicle Routing Problem by a Hybrid Evolutionary Algorithm. Sustainability. 2020; 12(5):2127. https://0-doi-org.brum.beds.ac.uk/10.3390/su12052127

Chicago/Turabian Style

Peng, Bo, Lifan Wu, Yuxin Yi, and Xiding Chen. 2020. "Solving the Multi-Depot Green Vehicle Routing Problem by a Hybrid Evolutionary Algorithm" Sustainability 12, no. 5: 2127. https://0-doi-org.brum.beds.ac.uk/10.3390/su12052127

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop