# An Adaptive Large Neighborhood Search for the Larger-Scale Instances of Green Vehicle Routing Problem with Time Windows

Authors: Zixuan Yu; Ping Zhang; Yang Yu; Wei Sun; Min Huang
Journal: Complexity
 (2020)

---

## Abstract

Due to huge amount of greenhouse gases emission (such as CO2), freight has been adversely affecting the global environment in facilitating the global economy. Therefore, green vehicle routing problem (GVRP), aiming to minimize the total carbon emissions in the transportation, has become a hot issue. In this paper, an adaptive large neighborhood search (ALNS) algorithm is proposed to solve large-scale instances of GVRP. The core of ALNS algorithm is destroy operators and repair operators. In the destroy operators, a new removal heuristic applying to the characteristics of GVRP is proposed. The heuristic can quickly remove customers who bring a large amount of carbon emissions with pertinence, and these customers may be arranged more properly in future repair operators. In the repair operators, a fast insertion method is developed. In the fast insertion method, the feasibility of a new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers. Thus, the computational time of the ALNS algorithm is greatly saved. Computational experiments were performed on Solomon benchmark with 100 customers and Homberger benchmark instances with up to 1000 customers. Given the same computational time, the proposed ALNS improves the average accuracy by 8.49% compared with the classic ALNS. In the optimal situation, the improvement can achieve 33.61%.


---

## Body

## 1. Introduction

The burning of fossil fuels is the main source of greenhouse gases, and the transportation industry is the largest consumer of petroleum products and the main culprit of global air pollution [1]. Therefore, in recent years, more and more attention has been paid to the green issue in the field of VRP [2].Green VRP focuses on energy consumption and environmental issues. In recent years, GVRP has been widely studied [2, 3] and the exact solution of GVRP can be obtained [4–6], but it is far from satisfying the requirements of practical application. Therefore, a lot of metaheuristic algorithms are developed. Among many metaheuristic methods, ALNS is an effective algorithm with outstanding performance in VRP [7], especially for the large-scale instances. Up to now, there are no published papers, in which the proposed algorithm can solve GVRP instances with more than 500 customer points in a reasonable computational time in a practical application. The existing ALNS for VRP/PDP can be used to solve GVRP by changing the objective function to the total carbon emission but did not consider the characteristics of GVRP. Thus, an ALNS algorithm is developed to solve the large-scale instances of GVRP according the characteristics of GVRP.


### 1.1. Related Work

The current research on green routing mainly includes three types of problems: pollution routing problem (PRP), energy minimizing vehicle routing problem (EMVRP), and green vehicle routing problem (GVRP).PRP consists of routing a number of vehicles to serve a set of customers within preset time windows and determining their speed on each route segment, so as to minimize a function comprising fuel, emission, and driver costs [8]. The two mixed-integer convex optimization models for PRP were proposed by Fukasawa et al. [9]. The nonlinear mixed-integer programming model for PRP was presented by Majidi et al. [10]. A bi-objective model for pickup and delivery PRP with minimizing total system cost and fuel consumption to reduce carbon emissions was proposed by Abad et al. [11]. For green transportation scheduling with pickup time and transport mode selections, an Evolution Strategy-Based Memetic Pareto Optimization was developed by Guo et al. [12].EMVRP aims to minimize the total energy used in the routing [13]. In EMVRP, the energy consumption over each arc is defined as the product of the arc length and the weight of the vehicle. For EMVRP, Fukasawa et al. [14] proposed two new mixed-integer linear programming formulations: an arc-load formulation and a set partitioning formulation, and proposed a branch-and-cut algorithm for the arc-load formulation and a branch-cut-and-price algorithm for the set partitioning formulation strengthened by additional constraints.Most studies of GVRP considered vehicle travel time and weight to be key factors in greenhouse gas emissions [3, 6, 15]. Bektaş et al. [3] believe that green vehicle routing is a VRP that considers explicitly carbon dioxide equivalent emissions and reviewed some of the latest developments in this field. Sawik et al. [16] focus on multiobjective GVRP while taking into account the minimization of noise, pollution, and fuel consumption.For the research on the exact solution of green routing, Bektaş and Laporte [4] developed a branch and price (BP) for PRP. A BP algorithm for time-dependent PRP was developed by Franceschetti et al. [17]. Yu et al. [18] used integrated scheduling for reducing carbon emissions of pickup and delivery problem (PDP). Yu et al. [6] developed a branch-and-price algorithm for heterogeneous fleet green vehicle routing problem (HFGVRP). The algorithm can solve the instances with 100 customers. Yu et al. [6] developed an exact algorithm for bi-objective green ride-sharing problem (BGRSP).Due to the complexity of GVRP, the exact algorithm can solve up to 100 customer instances [6], which is far from meeting the requirements of practical applications. Therefore, many studies have focused on heuristic or metaheuristic algorithms. Xiao et al. [19] developed a simulated annealing (SA) algorithm for CVRP with fuel consumption rate (FCVRP), which can solve up to 483 customers. Úbeda et al. [20] developed a tabu search (TS) algorithm for GVRP. Zhang et al. [21] developed a hybrid artificial bee colony algorithm for environmental VRP. Montoya et al. [22] developed a multispace sampling heuristic for GVRP. Ene et al. [23] developed a hybrid metaheuristic algorithm for HFGVRP with up to 100 customers. Shui and Szeto [24] developed a hybrid rolling horizon artificial bee colony algorithm for dynamic green bike repositioning problem. Wang and Szeto [25] developed a large neighborhood search in a bike-sharing network. Based on a multigraph reformulation, Andelmin and Bartolini [26] developed a multistart local search heuristic for GVRP with up to 470 customers. Li et al. [27] developed an ant colony optimization algorithm for multidepot GVRP with multiple objectives. Rastani et al. [28] developed an ALNS for VRP to reduce emissions stemming from logistics operations. Rezaei et al. [29] developed a GA algorithm and a SA algorithm for HFGVRP. Zhu and Hu [30] developed a hybrid algorithm (HA) for time-dependent green vehicle routing problem (TDGVRP) for up to 200 customers.Adaptive large neighborhood search (ALNS) is an effective neighborhood search heuristic for large-scale instances of NP-hard problems. ALNS is composed of a number of competing subheuristics that are used with a frequency corresponding to their historic performance [7]. ALNS greatly improves the common shortcoming of traditional neighborhood search (NS), i.e., the difficulty in moving out of the current search space.Due to the outstanding performance of ALNS, it is widely used in many variants of VRP. Laporte et al. [31] developed an ALNS for CARPSD (capacitated arc-routing problem with stochastic demands). Azi et al. [32] developed an ALNS algorithm for VRP with multiple routes. Keskin and Çatay [33] developed an ALNS for electric VRPTW. Mancini [34] developed a large neighborhood search (LNS) for hybrid VRP. Bach et al. [35] developed an ALNS for distance-constrained capacitated vehicle routing problem (DCVRP). Gu et al. [36] developed an ALNS for the VRP with commodity-constrained split delivery. Hof and Schneider [37] developed an ALNS for simultaneous pickup and delivery VRP with up to 1000 customers. Lahyani et al. [38] developed a hybrid ALNS for multidepot open VRP with up to 288 customers. Sacramento et al. [39] developed an ALNS for drones VRP with up to 375 customers.Up to now, there are no published papers, in which the proposed algorithm can solve GVRP instances with more than 500 customer points in a reasonable computational time in a practical application. Although the ALNS for VRP can be used to solve GVRP by changing the objective function, the ALNS did not consider the characteristics of GVRP. Thus, an ALNS according to the characteristics of GVRP is developed to solve the large-scale instances with 1000 customers of GVRP.


### 1.2. Our Contributions

The main contributions are as follows:(i)
An ALNS for GVRP is developed to solve the large-scale instances with 1000 customers. In the same time, the proposed ALNS improves the average accuracy by 8.49% compared with the classic ALNS. In the optimal situation, the improvement can achieve 33.61%.(ii)
A new type of removal heuristic in destroy operators was proposed according to the characteristics of GVRP.(iii)
To save computational time, a fast insertion method was proposed. In the fast insertion method, the feasibility of the new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers. Thus, the total computational time is saved by 78.9%.The rest of the paper is organized as follows: in Section2, GVRPTW is described and formulated; in Section 3, an ALNS algorithm for the large-scale instances with 1000 customers of GVRPTW is proposed. A new type of removal heuristic according to the characteristics of GVRP is proposed. To save computational time, a fast insertion method was proposed. In the fast insertion method, the feasibility of the new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers. Section 4 presents computational results, including the experiments of Solomon instances with 100 customers and Homberger instances with up to 1000 customers. In the same time, the proposed ALNS improves the average accuracy by 8.49% compared with the classic ALNS. In the optimal situation, the improvement can achieve 33.61%. Section 5 gives conclusions.


## 1.1. Related Work

The current research on green routing mainly includes three types of problems: pollution routing problem (PRP), energy minimizing vehicle routing problem (EMVRP), and green vehicle routing problem (GVRP).PRP consists of routing a number of vehicles to serve a set of customers within preset time windows and determining their speed on each route segment, so as to minimize a function comprising fuel, emission, and driver costs [8]. The two mixed-integer convex optimization models for PRP were proposed by Fukasawa et al. [9]. The nonlinear mixed-integer programming model for PRP was presented by Majidi et al. [10]. A bi-objective model for pickup and delivery PRP with minimizing total system cost and fuel consumption to reduce carbon emissions was proposed by Abad et al. [11]. For green transportation scheduling with pickup time and transport mode selections, an Evolution Strategy-Based Memetic Pareto Optimization was developed by Guo et al. [12].EMVRP aims to minimize the total energy used in the routing [13]. In EMVRP, the energy consumption over each arc is defined as the product of the arc length and the weight of the vehicle. For EMVRP, Fukasawa et al. [14] proposed two new mixed-integer linear programming formulations: an arc-load formulation and a set partitioning formulation, and proposed a branch-and-cut algorithm for the arc-load formulation and a branch-cut-and-price algorithm for the set partitioning formulation strengthened by additional constraints.Most studies of GVRP considered vehicle travel time and weight to be key factors in greenhouse gas emissions [3, 6, 15]. Bektaş et al. [3] believe that green vehicle routing is a VRP that considers explicitly carbon dioxide equivalent emissions and reviewed some of the latest developments in this field. Sawik et al. [16] focus on multiobjective GVRP while taking into account the minimization of noise, pollution, and fuel consumption.For the research on the exact solution of green routing, Bektaş and Laporte [4] developed a branch and price (BP) for PRP. A BP algorithm for time-dependent PRP was developed by Franceschetti et al. [17]. Yu et al. [18] used integrated scheduling for reducing carbon emissions of pickup and delivery problem (PDP). Yu et al. [6] developed a branch-and-price algorithm for heterogeneous fleet green vehicle routing problem (HFGVRP). The algorithm can solve the instances with 100 customers. Yu et al. [6] developed an exact algorithm for bi-objective green ride-sharing problem (BGRSP).Due to the complexity of GVRP, the exact algorithm can solve up to 100 customer instances [6], which is far from meeting the requirements of practical applications. Therefore, many studies have focused on heuristic or metaheuristic algorithms. Xiao et al. [19] developed a simulated annealing (SA) algorithm for CVRP with fuel consumption rate (FCVRP), which can solve up to 483 customers. Úbeda et al. [20] developed a tabu search (TS) algorithm for GVRP. Zhang et al. [21] developed a hybrid artificial bee colony algorithm for environmental VRP. Montoya et al. [22] developed a multispace sampling heuristic for GVRP. Ene et al. [23] developed a hybrid metaheuristic algorithm for HFGVRP with up to 100 customers. Shui and Szeto [24] developed a hybrid rolling horizon artificial bee colony algorithm for dynamic green bike repositioning problem. Wang and Szeto [25] developed a large neighborhood search in a bike-sharing network. Based on a multigraph reformulation, Andelmin and Bartolini [26] developed a multistart local search heuristic for GVRP with up to 470 customers. Li et al. [27] developed an ant colony optimization algorithm for multidepot GVRP with multiple objectives. Rastani et al. [28] developed an ALNS for VRP to reduce emissions stemming from logistics operations. Rezaei et al. [29] developed a GA algorithm and a SA algorithm for HFGVRP. Zhu and Hu [30] developed a hybrid algorithm (HA) for time-dependent green vehicle routing problem (TDGVRP) for up to 200 customers.Adaptive large neighborhood search (ALNS) is an effective neighborhood search heuristic for large-scale instances of NP-hard problems. ALNS is composed of a number of competing subheuristics that are used with a frequency corresponding to their historic performance [7]. ALNS greatly improves the common shortcoming of traditional neighborhood search (NS), i.e., the difficulty in moving out of the current search space.Due to the outstanding performance of ALNS, it is widely used in many variants of VRP. Laporte et al. [31] developed an ALNS for CARPSD (capacitated arc-routing problem with stochastic demands). Azi et al. [32] developed an ALNS algorithm for VRP with multiple routes. Keskin and Çatay [33] developed an ALNS for electric VRPTW. Mancini [34] developed a large neighborhood search (LNS) for hybrid VRP. Bach et al. [35] developed an ALNS for distance-constrained capacitated vehicle routing problem (DCVRP). Gu et al. [36] developed an ALNS for the VRP with commodity-constrained split delivery. Hof and Schneider [37] developed an ALNS for simultaneous pickup and delivery VRP with up to 1000 customers. Lahyani et al. [38] developed a hybrid ALNS for multidepot open VRP with up to 288 customers. Sacramento et al. [39] developed an ALNS for drones VRP with up to 375 customers.Up to now, there are no published papers, in which the proposed algorithm can solve GVRP instances with more than 500 customer points in a reasonable computational time in a practical application. Although the ALNS for VRP can be used to solve GVRP by changing the objective function, the ALNS did not consider the characteristics of GVRP. Thus, an ALNS according to the characteristics of GVRP is developed to solve the large-scale instances with 1000 customers of GVRP.


## 1.2. Our Contributions

The main contributions are as follows:(i)
An ALNS for GVRP is developed to solve the large-scale instances with 1000 customers. In the same time, the proposed ALNS improves the average accuracy by 8.49% compared with the classic ALNS. In the optimal situation, the improvement can achieve 33.61%.(ii)
A new type of removal heuristic in destroy operators was proposed according to the characteristics of GVRP.(iii)
To save computational time, a fast insertion method was proposed. In the fast insertion method, the feasibility of the new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers. Thus, the total computational time is saved by 78.9%.The rest of the paper is organized as follows: in Section2, GVRPTW is described and formulated; in Section 3, an ALNS algorithm for the large-scale instances with 1000 customers of GVRPTW is proposed. A new type of removal heuristic according to the characteristics of GVRP is proposed. To save computational time, a fast insertion method was proposed. In the fast insertion method, the feasibility of the new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers. Section 4 presents computational results, including the experiments of Solomon instances with 100 customers and Homberger instances with up to 1000 customers. In the same time, the proposed ALNS improves the average accuracy by 8.49% compared with the classic ALNS. In the optimal situation, the improvement can achieve 33.61%. Section 5 gives conclusions.


## 2. Description and Formulation of GVRPTW

### 2.1. Problem Description of GVRPTW

GVRPTW focuses on routing vehicles to serve a group of customers within their time windows to minimize the total carbon emissions.G=V,E, defined in Euclidean space, is a complete directed graph with V=0,1,…,n+1. Nodes 0 and n+1 denote the start and the end depot, respectively. Let V0=V\0,n+1 be the set of customers. For customer i, the quantity qi of products is to be served by the vehicle within the time window ai,bi. In addition, the service time si is required. E represents the set of arcs. For each arc i,j∈E, let dij be the distance from node i to node j. An unlimited fleet of vehicles is located at the depot. The objective of GVRPTW is to minimize the carbon emissions, while considering the following assumptions: (1) all customers must be served; (2) each customer must be visited only once with only one vehicle; (3) each feasible route starts and ends at the depot; (4) the service at each customer i must begin within the customer time window ai,bi, and the vehicle will wait if it arrives before the ready time ai; and (5) each vehicle needs to return to the depot before the due time bn+1.


### 2.2. Mathematical Model of GVRPTW

#### 2.2.1. Notations

The following notations are used:0: start depot.n + 1: end depot.v: vehicle speed.V: set of nodes V=0,1,…,n+1.V0: set of customers V0 = V\{0, n + 1}.i, j: nodes i, j ∈ V.dij: distance between node i and node j.fij: vehicle’s payload through arc (i, j).qi: quantity of products at the customer i.[ai, bi]: time window when service may start at node i, where ai and bi represent the earliest and latest time, respectively.si: service time at node i, and the service time at depot (i.e., node 0 and n + 1) is set to 0.eij: carbon emissions between node i and node j.ui: time of vehicle arriving at customer i.di: load of vehicle leaving from customer i.p: a feasible route visited by a vehicle, p=0,1,…,l,n+1, where l is number of customers in the route p.Q: capacity of vehicle.w: curb weight of vehicles.xij: decision variables; if the vehicle is driving from node i to node j, then xij = 1; otherwise, xij = 0.tij: travel time between node i and node j.


#### 2.2.2. Mathematical Formulation

The carbon emissions calculation formula is used [40]:(1)eij=rαl×dijv+α2×w+fij×dijv,where r represents carbon emissions index parameter, by which the fuel consumption on arc i,j can be transformed into carbon emissions; α1 represents the emissions unrelated to the vehicle mass; α2 represents the emissions linear to the vehicle mass; w represents curb weight of vehicles; and v represents the vehicle speed. It can be seen from equation (1) that the carbon emissions from node i and node j are related to the following factors: distance, vehicle speed, and the total weight of the vehicle between the two nodes.The mathematical programming model of GVRP is formulated as follows:(2)min∑i,j∈Eeijxij,(3)s.t.∑i,j∈Exij=1,∀i∈V0,(4)∑0,i∈Ex0i=∑i,n+1∈Exi,n+1,∀i∈V0,(5)∑j,i∈Exji=∑i,j∈Exij,∀i∈V0,(6)uj≥ui+si+tijxij,∀i,j∈E,(7)ai≤ui≤bi,∀i∈V0,(8)dj≥di+qjxij,∀i,j∈E,(9)0≤di≤Q,∀i∈V0,(10)ui≥0,∀i∈V,(11)xij∈0,1⋅∀i,j∈E.Equation (2) is the objective function, where eij is related to the vehicle’s current load fij and the distance dij between node i and node j. Therefore, equation (2) is nonlinear with decision variables xij. Constraint (3) means that each customer needs to be served once. Constraints (4) and (5) represent the flow balance of the depot and the customers, respectively. Constraints (6) and (8) represent arrival time update constraints and load update constraints, respectively. Constraint (7) is the time window constraints. Constraint (9) is the load constraints. Constraint (11) defines the domains of decision variables.Since the model is nonlinear, it cannot be solved by the solver for linear programming problem directly, such as CPLEX. To obtain the optimal solution of GVRP, some researchers developed the branch-and-price (BP) algorithm [6] and branch-and-cut and price (BCP) algorithm [14]. Although both BP and BCP are based on the set partitioning model transformed by the mathematical programming model, the BP and BCP algorithms can solve up to 100 customers. Therefore, some metaheuristic algorithms have been proposed. Up to now, there are no published papers, in which the proposed algorithm can solve GVRP instances with more than 500 customers in a reasonable computational time in a practical application. Although the ALNS for VRP can be used to solve GVRP by changing the objective function, the ALNS did not consider the characteristics of GVRP. Thus, an ALNS algorithm according to the characteristics of GVRP is developed. The proposed ALNS efficiently solves the large-scale instances with more than 500 customers even 1000 customers of GVRP.


## 2.1. Problem Description of GVRPTW

GVRPTW focuses on routing vehicles to serve a group of customers within their time windows to minimize the total carbon emissions.G=V,E, defined in Euclidean space, is a complete directed graph with V=0,1,…,n+1. Nodes 0 and n+1 denote the start and the end depot, respectively. Let V0=V\0,n+1 be the set of customers. For customer i, the quantity qi of products is to be served by the vehicle within the time window ai,bi. In addition, the service time si is required. E represents the set of arcs. For each arc i,j∈E, let dij be the distance from node i to node j. An unlimited fleet of vehicles is located at the depot. The objective of GVRPTW is to minimize the carbon emissions, while considering the following assumptions: (1) all customers must be served; (2) each customer must be visited only once with only one vehicle; (3) each feasible route starts and ends at the depot; (4) the service at each customer i must begin within the customer time window ai,bi, and the vehicle will wait if it arrives before the ready time ai; and (5) each vehicle needs to return to the depot before the due time bn+1.


## 2.2. Mathematical Model of GVRPTW

### 2.2.1. Notations

The following notations are used:0: start depot.n + 1: end depot.v: vehicle speed.V: set of nodes V=0,1,…,n+1.V0: set of customers V0 = V\{0, n + 1}.i, j: nodes i, j ∈ V.dij: distance between node i and node j.fij: vehicle’s payload through arc (i, j).qi: quantity of products at the customer i.[ai, bi]: time window when service may start at node i, where ai and bi represent the earliest and latest time, respectively.si: service time at node i, and the service time at depot (i.e., node 0 and n + 1) is set to 0.eij: carbon emissions between node i and node j.ui: time of vehicle arriving at customer i.di: load of vehicle leaving from customer i.p: a feasible route visited by a vehicle, p=0,1,…,l,n+1, where l is number of customers in the route p.Q: capacity of vehicle.w: curb weight of vehicles.xij: decision variables; if the vehicle is driving from node i to node j, then xij = 1; otherwise, xij = 0.tij: travel time between node i and node j.


### 2.2.2. Mathematical Formulation

The carbon emissions calculation formula is used [40]:(1)eij=rαl×dijv+α2×w+fij×dijv,where r represents carbon emissions index parameter, by which the fuel consumption on arc i,j can be transformed into carbon emissions; α1 represents the emissions unrelated to the vehicle mass; α2 represents the emissions linear to the vehicle mass; w represents curb weight of vehicles; and v represents the vehicle speed. It can be seen from equation (1) that the carbon emissions from node i and node j are related to the following factors: distance, vehicle speed, and the total weight of the vehicle between the two nodes.The mathematical programming model of GVRP is formulated as follows:(2)min∑i,j∈Eeijxij,(3)s.t.∑i,j∈Exij=1,∀i∈V0,(4)∑0,i∈Ex0i=∑i,n+1∈Exi,n+1,∀i∈V0,(5)∑j,i∈Exji=∑i,j∈Exij,∀i∈V0,(6)uj≥ui+si+tijxij,∀i,j∈E,(7)ai≤ui≤bi,∀i∈V0,(8)dj≥di+qjxij,∀i,j∈E,(9)0≤di≤Q,∀i∈V0,(10)ui≥0,∀i∈V,(11)xij∈0,1⋅∀i,j∈E.Equation (2) is the objective function, where eij is related to the vehicle’s current load fij and the distance dij between node i and node j. Therefore, equation (2) is nonlinear with decision variables xij. Constraint (3) means that each customer needs to be served once. Constraints (4) and (5) represent the flow balance of the depot and the customers, respectively. Constraints (6) and (8) represent arrival time update constraints and load update constraints, respectively. Constraint (7) is the time window constraints. Constraint (9) is the load constraints. Constraint (11) defines the domains of decision variables.Since the model is nonlinear, it cannot be solved by the solver for linear programming problem directly, such as CPLEX. To obtain the optimal solution of GVRP, some researchers developed the branch-and-price (BP) algorithm [6] and branch-and-cut and price (BCP) algorithm [14]. Although both BP and BCP are based on the set partitioning model transformed by the mathematical programming model, the BP and BCP algorithms can solve up to 100 customers. Therefore, some metaheuristic algorithms have been proposed. Up to now, there are no published papers, in which the proposed algorithm can solve GVRP instances with more than 500 customers in a reasonable computational time in a practical application. Although the ALNS for VRP can be used to solve GVRP by changing the objective function, the ALNS did not consider the characteristics of GVRP. Thus, an ALNS algorithm according to the characteristics of GVRP is developed. The proposed ALNS efficiently solves the large-scale instances with more than 500 customers even 1000 customers of GVRP.


## 2.2.1. Notations

The following notations are used:0: start depot.n + 1: end depot.v: vehicle speed.V: set of nodes V=0,1,…,n+1.V0: set of customers V0 = V\{0, n + 1}.i, j: nodes i, j ∈ V.dij: distance between node i and node j.fij: vehicle’s payload through arc (i, j).qi: quantity of products at the customer i.[ai, bi]: time window when service may start at node i, where ai and bi represent the earliest and latest time, respectively.si: service time at node i, and the service time at depot (i.e., node 0 and n + 1) is set to 0.eij: carbon emissions between node i and node j.ui: time of vehicle arriving at customer i.di: load of vehicle leaving from customer i.p: a feasible route visited by a vehicle, p=0,1,…,l,n+1, where l is number of customers in the route p.Q: capacity of vehicle.w: curb weight of vehicles.xij: decision variables; if the vehicle is driving from node i to node j, then xij = 1; otherwise, xij = 0.tij: travel time between node i and node j.


## 2.2.2. Mathematical Formulation

The carbon emissions calculation formula is used [40]:(1)eij=rαl×dijv+α2×w+fij×dijv,where r represents carbon emissions index parameter, by which the fuel consumption on arc i,j can be transformed into carbon emissions; α1 represents the emissions unrelated to the vehicle mass; α2 represents the emissions linear to the vehicle mass; w represents curb weight of vehicles; and v represents the vehicle speed. It can be seen from equation (1) that the carbon emissions from node i and node j are related to the following factors: distance, vehicle speed, and the total weight of the vehicle between the two nodes.The mathematical programming model of GVRP is formulated as follows:(2)min∑i,j∈Eeijxij,(3)s.t.∑i,j∈Exij=1,∀i∈V0,(4)∑0,i∈Ex0i=∑i,n+1∈Exi,n+1,∀i∈V0,(5)∑j,i∈Exji=∑i,j∈Exij,∀i∈V0,(6)uj≥ui+si+tijxij,∀i,j∈E,(7)ai≤ui≤bi,∀i∈V0,(8)dj≥di+qjxij,∀i,j∈E,(9)0≤di≤Q,∀i∈V0,(10)ui≥0,∀i∈V,(11)xij∈0,1⋅∀i,j∈E.Equation (2) is the objective function, where eij is related to the vehicle’s current load fij and the distance dij between node i and node j. Therefore, equation (2) is nonlinear with decision variables xij. Constraint (3) means that each customer needs to be served once. Constraints (4) and (5) represent the flow balance of the depot and the customers, respectively. Constraints (6) and (8) represent arrival time update constraints and load update constraints, respectively. Constraint (7) is the time window constraints. Constraint (9) is the load constraints. Constraint (11) defines the domains of decision variables.Since the model is nonlinear, it cannot be solved by the solver for linear programming problem directly, such as CPLEX. To obtain the optimal solution of GVRP, some researchers developed the branch-and-price (BP) algorithm [6] and branch-and-cut and price (BCP) algorithm [14]. Although both BP and BCP are based on the set partitioning model transformed by the mathematical programming model, the BP and BCP algorithms can solve up to 100 customers. Therefore, some metaheuristic algorithms have been proposed. Up to now, there are no published papers, in which the proposed algorithm can solve GVRP instances with more than 500 customers in a reasonable computational time in a practical application. Although the ALNS for VRP can be used to solve GVRP by changing the objective function, the ALNS did not consider the characteristics of GVRP. Thus, an ALNS algorithm according to the characteristics of GVRP is developed. The proposed ALNS efficiently solves the large-scale instances with more than 500 customers even 1000 customers of GVRP.


## 3. Adaptive Large Neighborhood Search for GVRPTW

### 3.1. Framework of Adaptive Large Neighborhood Search

ALNS, proposed by Ropke and Pisinger [7], adaptively chooses among a number of insertion and removal heuristics to intensify and diversify the search and the algorithm is robust and to some extent self-calibrating [41].To describe the framework of ALNS in brief, the pseudocode of ALNS algorithm is shown in Algorithm1.Algorithm 1:ALNS ().
Input: an initial solution (x), the number of removed customers z∈0,…,n.Output: best solution (xbest).xbest⟵x.repeatx′⟵x.(1) Perform destroy operators onx′. See Section 3.2.(2) Perform repair operators onx′. See Section 3.3.if (fx′<fx) or (simulated annealing criteria is met) thenx⟵x′.iffx′<fxbestthenxbest⟵x′.end ifend ifPerform adaptive method. See Section3.4.until termination condition is met.Outputxbest.The initial solution (x) is obtained by a construction heuristic, where each customer i generates route p=0,i,n+1, until all customers are dealt with.When judging whether to update the current solution, the simulated annealing acceptance criterion [8] is used to prevent getting trapped in a local minimum [7].Implicitly defining the solution space through destroy operators (i.e., step (1)) and repair operators (i.e., step (2)), rather than explicitly defining the solution space, is a significant difference between ALNS and traditional neighborhood search (NS). Therefore, ALNS greatly improves the common shortcoming of traditional NS, i.e., the difficulty in moving out of the current search space.


### 3.2. Destroy Operators

In the destroy operators, there are several removal heuristics. In each removal heuristic,x′ is updated by removing z customers from x′. The removal heuristics are selected dynamically according to their past performance [8]. The classic ALNS used random removal heuristic, worst removal heuristic, and Shaw removal heuristic in the destroy operations. In addition, forward load removal heuristic (Section 3.2.4) is proposed according to the characteristics of GVRP. For more removal heuristics, refer to Hemmelmayr et al. [42] and Alinaghian and Shokouhi [43].


#### 3.2.1. Random Removal Heuristic

Random removal heuristic was proposed by Ropke and Pisinger [7]. Random removal heuristic selects z customers randomly and removes them from the solution.


#### 3.2.2. Worst Removal Heuristic

Worst removal heuristic was proposed by Ropke and Pisinger [7]. Worst removal heuristic removes the z customers with the biggest difference in the carbon emissions. The difference in the carbon emissions is the difference in the carbon emissions between the original solution and the solution produced by removing the z customers.


#### 3.2.3. Shaw Removal Heuristic

Shaw removal heuristic was proposed by Shaw [44]. The idea of Shaw removal heuristic is to remove customers that are similar, in order to obtain a promising solution in repair operators.


#### 3.2.4. Forward Load Removal Heuristic

In GVRP, the carbon emissions between two customers are related to the distance and the load. Therefore, in the pickup/delivery problem, if the customer with a high demand but is served earlier/later, then the carbon emissions will be more than when the customer is served later/earlier. Thus, the forward load removal heuristic is proposed according to the characteristics of GVRP, to remove customers with high demand and served early.The pseudocode for forward load removal heuristic is shown in Algorithm2.Algorithm 2:Forward load removal heuristic ().
Input: a feasible solution (x), the number of removed customers (z), parameter c∈0,1.Output: set of removed customers.D⟵∅. D is the set of removed customers.H⟵∅. H is used to record the array of customers sorted by descending qi.M⟵0. M records the order of selected customer.whilez > 0 doRandomly choose a feasible routep=0,1,…,l,n+1∈x with more than 1 customers i.e., l > 1.H ⟵ sorting p by descending qi.M⟵⌊c∗l⌋.(1)forj=1toM−1do(2)ifqj>HM.demanddo(3) Removej from p.D⟵D∪j.z⟵z−1.end ifend forend whileOutputD.In step (1),M is the order of selected customer.In step (2), the customer beforeM in the route p is the customer served earlier, i.e., j=1toM−1.In step (3), if the customer’s demand is higher than the demand ofHM, i.e., qj>HM.demand, then the customer is the high-demand customer. Thus, the customer is removed from p.Steps (2) and (3) guarantee that the customers with high-demand served earlier are removed.Forward load removal heuristic can be seen as a compromise between random removal heuristic and worst removal heuristic. Compared to random removal heuristic, forward load removal heuristic not only increases the calculation time a little but also removes the customers with huge carbon emissions. Compared with worst removal heuristic, forward load removal heuristic does not calculate the carbon emissions of each customer but uses the characteristics of the problem to simplify the judgment of customers with huge carbon emissions, which greatly reduces the calculation time of a single heuristic.


### 3.3. Repair Operators

In the repair operators, there are several insertion heuristics. In each insertion heuristics,x′ is updated by reinserting z customers (removed in destroy operators) into x′. The insertion heuristics are selected dynamically according to their past performance [8]. The proposed ALNS used greedy insertion heuristic and regret heuristic in repair operators. In addition, to save computational time, a fast insertion method was proposed. In the fast insertion method, the feasibility of the new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers. Thus, the judgment time of the feasibility of the new route is significantly saved.


#### 3.3.1. Greedy Insertion Heuristic

Greedy insertion heuristic is proposed by Ropke and Pisinger [7]. In greedy insertion heuristic, for all the feasible positions corresponding to all customers that have been removed, find a customer with the smallest carbon emissions increase to insert.


#### 3.3.2. Regret Insertion Heuristic

Regret insertion heuristic is proposed by Ropke and Pisinger [7]. In the regret heuristic, customers are treated in the order of their regret value. Informally speaking, the insertion, which will bring regret most if it is not done now, will be chosen. Regret-2 and Regret-3 insertion heuristics [7] were used.


#### 3.3.3. Fast Insertion Method

In repair operators, the new route generated by inserting customer needs to be judged whether or not feasible. In traditional repair operators, the constraints of all the customers after the inserted customer are checked to judge whether the new route is feasible. It takes too much time. In some situations (i.e., Situations1 and 2), the feasibility of the new route can be judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all the customers. It is worth noting that this method is a general method and can be used in all variants of VRPTW.Notations used in the section are as follows:j: customer jr: inserted customer rp=0,1,…,j,j+1,…,l,n+1: a feasible route;pins=0,1,…,j,r,j+1,…,l,n+1: a route produced by inserting customer r into p;dtj: the departure time of customer jDT=dt0,dt1,…,dtj,dtj+1,…,dtl,dtn+1: departure time array of the customers in route p;DTins=dt0,dt1,…,dtj,dtr,dtj+1,…,dtl,dtn+1: departure time array of the customers in route pins.Situation 1.
If there exists a customerk after the inserted customer r in the new route (i.e., pins=0,1,…,j,r,j+1,…,k,,l,n+1), with dtkinDTins equal to dtkinDT, then the new route pins is a feasible route.Example 1.
In a feasible routep=0,1,2,3,4,n+1, DT=0,1,4,5,6,7 is the departure time array of the customers in p. For customer 9, t19=t92=1, a9,b9=2,3. Customer 9 can be inserted into customers 1 and 2 in p without violating the time window constraint of customer 2. In the new route pins=0,1,9,2,3,4,n+1, if dt2inDTins=4 is equal to dt2inDT, then pins must be a feasible route. Thus, it is unnecessary to check the constraints of customers 3 and 4.
Situation1 means that the feasibility of the new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers.Situation 2.
Ifr cannot be inserted into the route p, then r cannot be inserted into any route generated by the route p.Example 2.
If customer 9 cannot be inserted in the routep=0,1,2,n+1, then customer 9 cannot be inserted in pins=0,1,2,3,n+1 generated by p.
According to Situation2, the feasibility of the route can be prejudged before inserting r.
The program of the fast insertion method is shown in Algorithm3.
Step (1) initializes the parameters.
Step (2) checks whether the routep is generated by a route in which customer r cannot be inserted into (see situation 1). If yes, output flag = false; otherwise, continue to check the feasibility of the new route.
Step (3) checks the capacity constraint of the new route, i.e.,qp+qr≤Q, and the time constraint of customer r, i.e., dtj+tjr≤br are satisfied. If yes, continue to check the feasibility of the new route; otherwise, output flag = false.
Step (3-1) calculates the departure time of the inserted customerr.
Step (4) checks if there exists a customerk after the inserted customer r in the new route pins with dtkinDTins equal to dtkinDT (see situation 1).
Step (4-1) checks the time constraint of customerk, i.e., dtk−1+tk−1,k≤bk, is satisfied. If yes, continue to check the feasibility of the new route; otherwise, break and output flag = false.
Step (4-1-1) calculates the departure time of the customerk.
Step (4-1-2) judges whetherdtkinDTins is equal to dtkinDTins. If yes, the new route is feasible; therefore, flag is set as true and break. Thus, the feasibility of the new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers.
Step (4-2) means that the feasibility of the new route is judged by checking the constraints of all customers after the inserted customerr.
Experiments confirmed that the total computational time through the fast insertion method is saved by 78.9% compared with the traditional insertion method.
In addition, this method is a general method, not only applicable to GVRP but also to all VRPTW variants.Algorithm 3:Feasible solution judgment ().
Input: route p=0,1,…,j,j+1,…,l,n+1 and customer r.Output: feasibility (flag) of the insertionqp is the total load of the route p.flag⟵false.(1)ifp is generated by a route which r cannot be inserted into then(2) Prejudge the feasibility of the insertion (i.e., situation2).Output flag.elseif (qp+qr≤Q) and (dtj+tjr≤br) then(3)dtr⟵maxdtj+tjr,ar+sr(3-1)fork=j+1toldo(4)ifk=j+1thendtk−1⟵dtr.tk−1,k⟵tr,j+1.end ififdtk−1+tk−1,k≤bkthen(4-1)dtk⟵maxdtk−1+tk−1,k,ak+sk.(4-1-1)ifdtkinDTins=dtkinDTthen(4-1-2)flag⟵true. Judge whether route is feasible using partial customers after r (i.e.,situation 1).breakend ifelsebreakend ififk=lthen(4-2)flag⟵true.end ifend forend ifend ifOutput flag.


### 3.4. Adaptive Method

Adaptive method is used to change the probability of heuristics in destroy operators and repair operators according to the performance of the heuristics [7]. The proposed ALNS algorithm uses a probability update formula as follows:(12)ait+1=ρait+fxt−fxt+1fxt,where ait is the probability of the algorithm i being called in the tth iteration, ρ∈0,1 is the evaporation coefficient, and fxt is the objective value of the solution xt in iteration t.Equation (12) can be seen as a combination of the evaporation method and reward method. To prevent the algorithm from prematurely converging, the probability is given within a fix range amin,amax.


### 3.5. Flowchart of ALNS for GVRP

The flowchart of ALNS for GVRP is shown in Figure1.Figure 1
The flowchart of ALNS for GVRP.As shown in Figure1, the data of a test instance are read first. For more information about experiment instances, refer to Section 4.2. Then, an intimal solution (x) is obtained by construction heuristic.In the destroy operators (Section3.2), one of the three removal heuristics is selected according to their past performance, and the initial probability of each heuristic is equal to 1/3.In the repair operators (Section3.3), one of the three insertion heuristics is selected according to their past performance, and the initial probability of each heuristic is equal to 1/3. The fast insertion method is embedded in all insertion heuristics in the repair operators. When judging the feasibility of each newly generated route, it is used to save judging time.If the solution (x′) generated in this iteration is better than (carbon emission is smaller than) the best solution (xbest), update the best solution, i.e., xbest←x′.If the solution (x′) generated in this iteration is better than the current solution (x) or meets the simulated annealing acceptance criteria [8], update the current solution, i.e., x←x′.The probability of the heuristics being called in the destroy operators and repair operators keeps changing, according to its performance in the adaptive method (Section3.4).When the given running time is reached, the best solution (xbest) is output and the algorithm ends.


## 3.1. Framework of Adaptive Large Neighborhood Search

ALNS, proposed by Ropke and Pisinger [7], adaptively chooses among a number of insertion and removal heuristics to intensify and diversify the search and the algorithm is robust and to some extent self-calibrating [41].To describe the framework of ALNS in brief, the pseudocode of ALNS algorithm is shown in Algorithm1.Algorithm 1:ALNS ().
Input: an initial solution (x), the number of removed customers z∈0,…,n.Output: best solution (xbest).xbest⟵x.repeatx′⟵x.(1) Perform destroy operators onx′. See Section 3.2.(2) Perform repair operators onx′. See Section 3.3.if (fx′<fx) or (simulated annealing criteria is met) thenx⟵x′.iffx′<fxbestthenxbest⟵x′.end ifend ifPerform adaptive method. See Section3.4.until termination condition is met.Outputxbest.The initial solution (x) is obtained by a construction heuristic, where each customer i generates route p=0,i,n+1, until all customers are dealt with.When judging whether to update the current solution, the simulated annealing acceptance criterion [8] is used to prevent getting trapped in a local minimum [7].Implicitly defining the solution space through destroy operators (i.e., step (1)) and repair operators (i.e., step (2)), rather than explicitly defining the solution space, is a significant difference between ALNS and traditional neighborhood search (NS). Therefore, ALNS greatly improves the common shortcoming of traditional NS, i.e., the difficulty in moving out of the current search space.


## 3.2. Destroy Operators

In the destroy operators, there are several removal heuristics. In each removal heuristic,x′ is updated by removing z customers from x′. The removal heuristics are selected dynamically according to their past performance [8]. The classic ALNS used random removal heuristic, worst removal heuristic, and Shaw removal heuristic in the destroy operations. In addition, forward load removal heuristic (Section 3.2.4) is proposed according to the characteristics of GVRP. For more removal heuristics, refer to Hemmelmayr et al. [42] and Alinaghian and Shokouhi [43].


### 3.2.1. Random Removal Heuristic

Random removal heuristic was proposed by Ropke and Pisinger [7]. Random removal heuristic selects z customers randomly and removes them from the solution.


### 3.2.2. Worst Removal Heuristic

Worst removal heuristic was proposed by Ropke and Pisinger [7]. Worst removal heuristic removes the z customers with the biggest difference in the carbon emissions. The difference in the carbon emissions is the difference in the carbon emissions between the original solution and the solution produced by removing the z customers.


### 3.2.3. Shaw Removal Heuristic

Shaw removal heuristic was proposed by Shaw [44]. The idea of Shaw removal heuristic is to remove customers that are similar, in order to obtain a promising solution in repair operators.


### 3.2.4. Forward Load Removal Heuristic

In GVRP, the carbon emissions between two customers are related to the distance and the load. Therefore, in the pickup/delivery problem, if the customer with a high demand but is served earlier/later, then the carbon emissions will be more than when the customer is served later/earlier. Thus, the forward load removal heuristic is proposed according to the characteristics of GVRP, to remove customers with high demand and served early.The pseudocode for forward load removal heuristic is shown in Algorithm2.Algorithm 2:Forward load removal heuristic ().
Input: a feasible solution (x), the number of removed customers (z), parameter c∈0,1.Output: set of removed customers.D⟵∅. D is the set of removed customers.H⟵∅. H is used to record the array of customers sorted by descending qi.M⟵0. M records the order of selected customer.whilez > 0 doRandomly choose a feasible routep=0,1,…,l,n+1∈x with more than 1 customers i.e., l > 1.H ⟵ sorting p by descending qi.M⟵⌊c∗l⌋.(1)forj=1toM−1do(2)ifqj>HM.demanddo(3) Removej from p.D⟵D∪j.z⟵z−1.end ifend forend whileOutputD.In step (1),M is the order of selected customer.In step (2), the customer beforeM in the route p is the customer served earlier, i.e., j=1toM−1.In step (3), if the customer’s demand is higher than the demand ofHM, i.e., qj>HM.demand, then the customer is the high-demand customer. Thus, the customer is removed from p.Steps (2) and (3) guarantee that the customers with high-demand served earlier are removed.Forward load removal heuristic can be seen as a compromise between random removal heuristic and worst removal heuristic. Compared to random removal heuristic, forward load removal heuristic not only increases the calculation time a little but also removes the customers with huge carbon emissions. Compared with worst removal heuristic, forward load removal heuristic does not calculate the carbon emissions of each customer but uses the characteristics of the problem to simplify the judgment of customers with huge carbon emissions, which greatly reduces the calculation time of a single heuristic.


## 3.2.1. Random Removal Heuristic

Random removal heuristic was proposed by Ropke and Pisinger [7]. Random removal heuristic selects z customers randomly and removes them from the solution.


## 3.2.2. Worst Removal Heuristic

Worst removal heuristic was proposed by Ropke and Pisinger [7]. Worst removal heuristic removes the z customers with the biggest difference in the carbon emissions. The difference in the carbon emissions is the difference in the carbon emissions between the original solution and the solution produced by removing the z customers.


## 3.2.3. Shaw Removal Heuristic

Shaw removal heuristic was proposed by Shaw [44]. The idea of Shaw removal heuristic is to remove customers that are similar, in order to obtain a promising solution in repair operators.


## 3.2.4. Forward Load Removal Heuristic

In GVRP, the carbon emissions between two customers are related to the distance and the load. Therefore, in the pickup/delivery problem, if the customer with a high demand but is served earlier/later, then the carbon emissions will be more than when the customer is served later/earlier. Thus, the forward load removal heuristic is proposed according to the characteristics of GVRP, to remove customers with high demand and served early.The pseudocode for forward load removal heuristic is shown in Algorithm2.Algorithm 2:Forward load removal heuristic ().
Input: a feasible solution (x), the number of removed customers (z), parameter c∈0,1.Output: set of removed customers.D⟵∅. D is the set of removed customers.H⟵∅. H is used to record the array of customers sorted by descending qi.M⟵0. M records the order of selected customer.whilez > 0 doRandomly choose a feasible routep=0,1,…,l,n+1∈x with more than 1 customers i.e., l > 1.H ⟵ sorting p by descending qi.M⟵⌊c∗l⌋.(1)forj=1toM−1do(2)ifqj>HM.demanddo(3) Removej from p.D⟵D∪j.z⟵z−1.end ifend forend whileOutputD.In step (1),M is the order of selected customer.In step (2), the customer beforeM in the route p is the customer served earlier, i.e., j=1toM−1.In step (3), if the customer’s demand is higher than the demand ofHM, i.e., qj>HM.demand, then the customer is the high-demand customer. Thus, the customer is removed from p.Steps (2) and (3) guarantee that the customers with high-demand served earlier are removed.Forward load removal heuristic can be seen as a compromise between random removal heuristic and worst removal heuristic. Compared to random removal heuristic, forward load removal heuristic not only increases the calculation time a little but also removes the customers with huge carbon emissions. Compared with worst removal heuristic, forward load removal heuristic does not calculate the carbon emissions of each customer but uses the characteristics of the problem to simplify the judgment of customers with huge carbon emissions, which greatly reduces the calculation time of a single heuristic.


## 3.3. Repair Operators

In the repair operators, there are several insertion heuristics. In each insertion heuristics,x′ is updated by reinserting z customers (removed in destroy operators) into x′. The insertion heuristics are selected dynamically according to their past performance [8]. The proposed ALNS used greedy insertion heuristic and regret heuristic in repair operators. In addition, to save computational time, a fast insertion method was proposed. In the fast insertion method, the feasibility of the new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers. Thus, the judgment time of the feasibility of the new route is significantly saved.


### 3.3.1. Greedy Insertion Heuristic

Greedy insertion heuristic is proposed by Ropke and Pisinger [7]. In greedy insertion heuristic, for all the feasible positions corresponding to all customers that have been removed, find a customer with the smallest carbon emissions increase to insert.


### 3.3.2. Regret Insertion Heuristic

Regret insertion heuristic is proposed by Ropke and Pisinger [7]. In the regret heuristic, customers are treated in the order of their regret value. Informally speaking, the insertion, which will bring regret most if it is not done now, will be chosen. Regret-2 and Regret-3 insertion heuristics [7] were used.


### 3.3.3. Fast Insertion Method

In repair operators, the new route generated by inserting customer needs to be judged whether or not feasible. In traditional repair operators, the constraints of all the customers after the inserted customer are checked to judge whether the new route is feasible. It takes too much time. In some situations (i.e., Situations1 and 2), the feasibility of the new route can be judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all the customers. It is worth noting that this method is a general method and can be used in all variants of VRPTW.Notations used in the section are as follows:j: customer jr: inserted customer rp=0,1,…,j,j+1,…,l,n+1: a feasible route;pins=0,1,…,j,r,j+1,…,l,n+1: a route produced by inserting customer r into p;dtj: the departure time of customer jDT=dt0,dt1,…,dtj,dtj+1,…,dtl,dtn+1: departure time array of the customers in route p;DTins=dt0,dt1,…,dtj,dtr,dtj+1,…,dtl,dtn+1: departure time array of the customers in route pins.Situation 1.
If there exists a customerk after the inserted customer r in the new route (i.e., pins=0,1,…,j,r,j+1,…,k,,l,n+1), with dtkinDTins equal to dtkinDT, then the new route pins is a feasible route.Example 1.
In a feasible routep=0,1,2,3,4,n+1, DT=0,1,4,5,6,7 is the departure time array of the customers in p. For customer 9, t19=t92=1, a9,b9=2,3. Customer 9 can be inserted into customers 1 and 2 in p without violating the time window constraint of customer 2. In the new route pins=0,1,9,2,3,4,n+1, if dt2inDTins=4 is equal to dt2inDT, then pins must be a feasible route. Thus, it is unnecessary to check the constraints of customers 3 and 4.
Situation1 means that the feasibility of the new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers.Situation 2.
Ifr cannot be inserted into the route p, then r cannot be inserted into any route generated by the route p.Example 2.
If customer 9 cannot be inserted in the routep=0,1,2,n+1, then customer 9 cannot be inserted in pins=0,1,2,3,n+1 generated by p.
According to Situation2, the feasibility of the route can be prejudged before inserting r.
The program of the fast insertion method is shown in Algorithm3.
Step (1) initializes the parameters.
Step (2) checks whether the routep is generated by a route in which customer r cannot be inserted into (see situation 1). If yes, output flag = false; otherwise, continue to check the feasibility of the new route.
Step (3) checks the capacity constraint of the new route, i.e.,qp+qr≤Q, and the time constraint of customer r, i.e., dtj+tjr≤br are satisfied. If yes, continue to check the feasibility of the new route; otherwise, output flag = false.
Step (3-1) calculates the departure time of the inserted customerr.
Step (4) checks if there exists a customerk after the inserted customer r in the new route pins with dtkinDTins equal to dtkinDT (see situation 1).
Step (4-1) checks the time constraint of customerk, i.e., dtk−1+tk−1,k≤bk, is satisfied. If yes, continue to check the feasibility of the new route; otherwise, break and output flag = false.
Step (4-1-1) calculates the departure time of the customerk.
Step (4-1-2) judges whetherdtkinDTins is equal to dtkinDTins. If yes, the new route is feasible; therefore, flag is set as true and break. Thus, the feasibility of the new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers.
Step (4-2) means that the feasibility of the new route is judged by checking the constraints of all customers after the inserted customerr.
Experiments confirmed that the total computational time through the fast insertion method is saved by 78.9% compared with the traditional insertion method.
In addition, this method is a general method, not only applicable to GVRP but also to all VRPTW variants.Algorithm 3:Feasible solution judgment ().
Input: route p=0,1,…,j,j+1,…,l,n+1 and customer r.Output: feasibility (flag) of the insertionqp is the total load of the route p.flag⟵false.(1)ifp is generated by a route which r cannot be inserted into then(2) Prejudge the feasibility of the insertion (i.e., situation2).Output flag.elseif (qp+qr≤Q) and (dtj+tjr≤br) then(3)dtr⟵maxdtj+tjr,ar+sr(3-1)fork=j+1toldo(4)ifk=j+1thendtk−1⟵dtr.tk−1,k⟵tr,j+1.end ififdtk−1+tk−1,k≤bkthen(4-1)dtk⟵maxdtk−1+tk−1,k,ak+sk.(4-1-1)ifdtkinDTins=dtkinDTthen(4-1-2)flag⟵true. Judge whether route is feasible using partial customers after r (i.e.,situation 1).breakend ifelsebreakend ififk=lthen(4-2)flag⟵true.end ifend forend ifend ifOutput flag.


## 3.3.1. Greedy Insertion Heuristic

Greedy insertion heuristic is proposed by Ropke and Pisinger [7]. In greedy insertion heuristic, for all the feasible positions corresponding to all customers that have been removed, find a customer with the smallest carbon emissions increase to insert.


## 3.3.2. Regret Insertion Heuristic

Regret insertion heuristic is proposed by Ropke and Pisinger [7]. In the regret heuristic, customers are treated in the order of their regret value. Informally speaking, the insertion, which will bring regret most if it is not done now, will be chosen. Regret-2 and Regret-3 insertion heuristics [7] were used.


## 3.3.3. Fast Insertion Method

In repair operators, the new route generated by inserting customer needs to be judged whether or not feasible. In traditional repair operators, the constraints of all the customers after the inserted customer are checked to judge whether the new route is feasible. It takes too much time. In some situations (i.e., Situations1 and 2), the feasibility of the new route can be judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all the customers. It is worth noting that this method is a general method and can be used in all variants of VRPTW.Notations used in the section are as follows:j: customer jr: inserted customer rp=0,1,…,j,j+1,…,l,n+1: a feasible route;pins=0,1,…,j,r,j+1,…,l,n+1: a route produced by inserting customer r into p;dtj: the departure time of customer jDT=dt0,dt1,…,dtj,dtj+1,…,dtl,dtn+1: departure time array of the customers in route p;DTins=dt0,dt1,…,dtj,dtr,dtj+1,…,dtl,dtn+1: departure time array of the customers in route pins.Situation 1.
If there exists a customerk after the inserted customer r in the new route (i.e., pins=0,1,…,j,r,j+1,…,k,,l,n+1), with dtkinDTins equal to dtkinDT, then the new route pins is a feasible route.Example 1.
In a feasible routep=0,1,2,3,4,n+1, DT=0,1,4,5,6,7 is the departure time array of the customers in p. For customer 9, t19=t92=1, a9,b9=2,3. Customer 9 can be inserted into customers 1 and 2 in p without violating the time window constraint of customer 2. In the new route pins=0,1,9,2,3,4,n+1, if dt2inDTins=4 is equal to dt2inDT, then pins must be a feasible route. Thus, it is unnecessary to check the constraints of customers 3 and 4.
Situation1 means that the feasibility of the new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers.Situation 2.
Ifr cannot be inserted into the route p, then r cannot be inserted into any route generated by the route p.Example 2.
If customer 9 cannot be inserted in the routep=0,1,2,n+1, then customer 9 cannot be inserted in pins=0,1,2,3,n+1 generated by p.
According to Situation2, the feasibility of the route can be prejudged before inserting r.
The program of the fast insertion method is shown in Algorithm3.
Step (1) initializes the parameters.
Step (2) checks whether the routep is generated by a route in which customer r cannot be inserted into (see situation 1). If yes, output flag = false; otherwise, continue to check the feasibility of the new route.
Step (3) checks the capacity constraint of the new route, i.e.,qp+qr≤Q, and the time constraint of customer r, i.e., dtj+tjr≤br are satisfied. If yes, continue to check the feasibility of the new route; otherwise, output flag = false.
Step (3-1) calculates the departure time of the inserted customerr.
Step (4) checks if there exists a customerk after the inserted customer r in the new route pins with dtkinDTins equal to dtkinDT (see situation 1).
Step (4-1) checks the time constraint of customerk, i.e., dtk−1+tk−1,k≤bk, is satisfied. If yes, continue to check the feasibility of the new route; otherwise, break and output flag = false.
Step (4-1-1) calculates the departure time of the customerk.
Step (4-1-2) judges whetherdtkinDTins is equal to dtkinDTins. If yes, the new route is feasible; therefore, flag is set as true and break. Thus, the feasibility of the new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers.
Step (4-2) means that the feasibility of the new route is judged by checking the constraints of all customers after the inserted customerr.
Experiments confirmed that the total computational time through the fast insertion method is saved by 78.9% compared with the traditional insertion method.
In addition, this method is a general method, not only applicable to GVRP but also to all VRPTW variants.Algorithm 3:Feasible solution judgment ().
Input: route p=0,1,…,j,j+1,…,l,n+1 and customer r.Output: feasibility (flag) of the insertionqp is the total load of the route p.flag⟵false.(1)ifp is generated by a route which r cannot be inserted into then(2) Prejudge the feasibility of the insertion (i.e., situation2).Output flag.elseif (qp+qr≤Q) and (dtj+tjr≤br) then(3)dtr⟵maxdtj+tjr,ar+sr(3-1)fork=j+1toldo(4)ifk=j+1thendtk−1⟵dtr.tk−1,k⟵tr,j+1.end ififdtk−1+tk−1,k≤bkthen(4-1)dtk⟵maxdtk−1+tk−1,k,ak+sk.(4-1-1)ifdtkinDTins=dtkinDTthen(4-1-2)flag⟵true. Judge whether route is feasible using partial customers after r (i.e.,situation 1).breakend ifelsebreakend ififk=lthen(4-2)flag⟵true.end ifend forend ifend ifOutput flag.


## 3.4. Adaptive Method

Adaptive method is used to change the probability of heuristics in destroy operators and repair operators according to the performance of the heuristics [7]. The proposed ALNS algorithm uses a probability update formula as follows:(12)ait+1=ρait+fxt−fxt+1fxt,where ait is the probability of the algorithm i being called in the tth iteration, ρ∈0,1 is the evaporation coefficient, and fxt is the objective value of the solution xt in iteration t.Equation (12) can be seen as a combination of the evaporation method and reward method. To prevent the algorithm from prematurely converging, the probability is given within a fix range amin,amax.


## 3.5. Flowchart of ALNS for GVRP

The flowchart of ALNS for GVRP is shown in Figure1.Figure 1
The flowchart of ALNS for GVRP.As shown in Figure1, the data of a test instance are read first. For more information about experiment instances, refer to Section 4.2. Then, an intimal solution (x) is obtained by construction heuristic.In the destroy operators (Section3.2), one of the three removal heuristics is selected according to their past performance, and the initial probability of each heuristic is equal to 1/3.In the repair operators (Section3.3), one of the three insertion heuristics is selected according to their past performance, and the initial probability of each heuristic is equal to 1/3. The fast insertion method is embedded in all insertion heuristics in the repair operators. When judging the feasibility of each newly generated route, it is used to save judging time.If the solution (x′) generated in this iteration is better than (carbon emission is smaller than) the best solution (xbest), update the best solution, i.e., xbest←x′.If the solution (x′) generated in this iteration is better than the current solution (x) or meets the simulated annealing acceptance criteria [8], update the current solution, i.e., x←x′.The probability of the heuristics being called in the destroy operators and repair operators keeps changing, according to its performance in the adaptive method (Section3.4).When the given running time is reached, the best solution (xbest) is output and the algorithm ends.


## 4. Computational Results

### 4.1. Software and Hardware Specifications

The BP algorithm was coded in C#, whose RMP and integer branch method were solved by ILOG CPLEX 12.6. Both classic ALNS and the proposed ALNS algorithms are coded in C#. All tests were executed on a laptop with a 2.80-GHz Intel Core TM i5-4200 processor using the Microsoft Windows 10 operating system using 12.00 GB of RAM. The classic ALNS and the proposed ALNS were run 10 times in each instance.


### 4.2. Experiment Instances

#### 4.2.1. Introduction of the Benchmark Instance

The used GVRPTW instances are derived from Solomon VRPTW benchmark instances [45] with 100 customers and Homberger VRPTW benchmark instances [46] with 200, 400, 600, 800, and 1000 customers. Benchmark instances are divided in three classes according to the geographical distribution of the customers: R (random), C (clustered), or RC (semiclustered) [45]. For more details on the benchmark instances, refer to [45, 46].


#### 4.2.2. Construction of Data Acquisition System

The standard test sets from Solomon benchmark and Homberger benchmark are used in the experiment, the coordinate, and the quantity of products, and time windows and service time of the depot and customers can be obtained.Based on the coordinate of the depot and customers, the Euclidean distancedij from node i to node j is calculated. Since the vehicle speed is fixed, the calculated distance needs to be scaled to meet the time window constraints, and the distance dij was changed to dij′=Cddij [47] and kept to four decimal places. For the quantity of products of each customer, the quantity was scaled to qi′=Cqqi. For simplicity of computation, the time windows were set as 24 h, and time windows of all customers were scaled accordingly using the coefficient 24/b0. Moreover, the time window ai,bi of customer i in Solomon instances and Homberger instances was modified as ai24/b0,bi24/b0 [47]. Furthermore, the service time was set to 1.5 h for all customers [47] in Solomon instances (100 customers), and the service time si was modified as si24/b0 in the Homberger instances (200, 400, 600, 800, and 1000 customers).According to the investigation in a logistics company in Shenyang, China, two types of vehicles for the computational tests of Homberger instance [46, 48, 49] were selected by the authors. Table 1 describes the value of parameters used in experiments of GVRP.Table 1
Value of used parameters of GVRP.
NotationDescriptionTypical valuesRCarbon emissions index parameter3.164α1Emissions unrelated to the vehicle mass1.078086α2Emissions linear to the vehicle mass0.072253W1Represents curb weight of vehicle type 12500 kgQ1Capacity of vehicle 12000 kgW2Represents curb weight of vehicle type 29000 kgQ2Capacity of vehicle 210000 kgVVehicle speed42 km/hCdDistance coefficient2N=1000.5200≤N≤8000.1N=1000CqLoad coefficient10Table2 describes the value of parameters used in experiments of ALNS.In the experiment, the running time of the heuristic algorithm is given according to the scale of the problem (Table3). Given the same running time, the performance of the algorithm is evaluated by its minimal carbon emissions.


### 4.3. Computational Results of Improved ALNS Algorithm

Since the branch-and-price algorithm developed by Yu et al. [50] can solve up to 100 customer instances, for 100-customer instances, the proposed algorithm is compared with the exact solution to evaluate its performance. For instances with more than 200 customers, there is no exact solution until now, and the proposed algorithm is compared with the classic ALNS algorithm proposed by Ropke and Pisinger [7].The destroy operators of classic ALNS consist of random removal heuristic, worst removal heuristic, and Shaw removal heuristic, and the repair operators of classic ALNS consist of greedy insertion heuristic and regret insertion heuristic (regret 2 and regret 3). In addition, the objective function is changed to the total carbon emission in the routing.The proposed ALNS in the paper and its details can refer to Section3.5.Table4 shows the results of Solomon benchmark instances with 100 customers using the vehicles with Q = 2000.For the instances with 100 customers, only 4 instances can be exactly solved by the branch-and-price for GVRPTW [6]. As shown in Table 4, in the same given computational time (20 s), the proposed ALNS improves the accuracy significantly compared with the classic ALNS.Tables5–9 show the results of Homberger benchmark instances with 200, 400, 600, 800, and 1000 customers, respectively.Table 2
Value of used parameters of ALNS.
NotationDescriptionTypical valuesZNumber of removed customers8N≤10020200≤N≤1000CCoefficient in forward load removal heuristic0.5aminMinimum value of the called probability0.01amaxMaximum value of the called probability0.99ρEvaporation parameter in adaptive method0.999Table 3
Given computational time of the heuristics for the instances with different customers.
Number of customers1002004006008001000GCT (s)204080120160200GCT: given computational time; unit is seconds.Table 4
Computational results of the instances with 100 customers.
InstanceExactClassic ALNSThe proposed ALNSECECT (s)MCEGap (%)MCEGap (%)r1011127.68171.721166.283.421164.653.28r2011027.1222363.171095.676.671057.082.92c101818.5644.20850.543.91835.152.03c105597.0332.72610.562.27610.562.27Average4.072.62ECE: exact carbon emissions solved by the BP algorithm proposed by Yu et al. [50]; MCE: minimal carbon emission; gap: the relative difference between the found solution and the exact solution (i.e., Gap=MCE−ECE/ECE∗100%).Table 5
Computational results of the instances with 200 customers.
InstanceClassic ALNSOur ALNSRCE (%)QOIMCEITERMCEITERR1_2_1813.73984794.4527122.372.76R1_2_2721.851129708.6322901.832.03R2_2_11096.214331029.4114746.093.40R2_2_21004.70367924.469627.992.62Average3.432.70C1_2_1827.44945683.45637917.406.75C1_2_2727.80899587.61433319.264.82C2_2_1655.25519609.0027657.065.33C2_2_2626.34479572.0912098.662.53Average13.104.86RC1_2_1648.921024633.6922632.352.21RC1_2_2620.041175598.6419913.451.70RC2_2_1880.85404765.83188013.064.66RC2_2_2746.70447726.079562.762.14Average5.412.68Total average7.313.41ITER: number of iterations; RCEs: reduced carbon emissions (i.e.,RCE=MCE of classic ALNS−MCE of our ALNS/MCE of classic ALNS∗100%); QOI: quotient of number of iterations (i.e., QOI=ITER of our ALNS/ITER of classic ALNS).Table 6
Computational results of the instances with 400 customers.
InstanceClassic ALNSOur ALNSRCE (%)QOIMCEITERMCEITERR1_4_11952.9017551887.1055633.373.17R1_4_21793.6119451743.0544492.822.29R2_4_12686.054752319.57284513.645.99R2_4_22668.673842220.06150916.813.93Average9.163.84C1_4_12085.398721881.1961159.797.02C1_4_21810.1510291680.6435267.153.43C2_4_11818.395741388.01414423.677.22C2_4_21621.954821341.20156417.313.25Average14.485.23RC1_4_11689.1316301652.5135392.172.17RC1_4_21578.4920651561.5837511.071.82RC2_4_11918.344601823.6623894.945.20RC2_4_21846.804161667.0310589.732.54Average4.482.93Total average9.374.00Table 7
Computational results of the instances with 600 customers.
InstanceClassic ALNSOur ALNSRCE (%)QOIMCEITERMCEITERR1_6_14298.0815174224.5554511.713.59R1_6_23851.5616743725.9144083.262.63R2_6_15264.923954718.97392510.379.93R2_6_25153.552994152.05180819.436.04Average8.695.55C1_6_14328.348283948.5758778.777.10C1_6_23950.988743499.99442111.415.06C2_6_13483.274862991.08385214.137.93C2_6_23305.614482966.92136410.253.05Average11.145.79RC1_6_13483.1413373408.3940892.153.06RC1_6_23272.8017453251.0536800.662.11RC2_6_13975.134903761.9620155.364.11RC2_6_23465.404443366.2111672.862.63Average2.762.98Total average7.534.77Table 8
Computational results of the instances with 800 customers.
InstanceClassic ALNSOur ALNSRCE (%)QOIMCEITERMCEITERR1_8_17417.0715737210.6860802.783.86R1_8_26868.4918376627.2946883.512.55R2_8_19113.605328206.3036969.966.95R2_8_27803.534606995.30176710.363.84Average6.654.30C1_8_19107.378967177.04550321.206.14C1_8_26824.469426348.2241946.984.45C2_8_15550.204995032.8838939.327.81C2_8_25484.324204606.45279816.016.66Average13.386.26RC1_8_16245.8413966155.8552401.443.75RC1_8_25954.1016045779.1539832.942.48RC2_8_16894.014336575.1228634.636.61RC2_8_26451.263996074.7318355.844.59Average3.714.36Total average7.914.97Table 9
Computational results of the instances with 1000 customers.
InstanceClassic ALNSOur ALNSRCE %QOIMCEITERMCEITERR1_10_12217.4931642201.5558120.721.84R1_10_22095.0733602042.5751632.511.54R2_10_12953.183172605.26153211.784.84R2_10_22717.363102352.83102013.413.29Average7.112.88C1_10_12741.778162412.21466412.025.72C1_10_22781.218442334.44426516.065.05C2_10_11237.293991123.4628369.207.11C2_10_21618.343591074.46149633.614.16Average17.725.51RC1_10_11977.3442271948.0463071.481.49RC1_10_21831.0231751821.7043980.511.39RC2_10_12067.783261883.5812008.913.68RC2_10_22002.163051729.2280813.632.65Average6.132.30Total average10.323.56It can be seen from Table5 that, in the instances of 200 customers, the proposed algorithm is better than the classic algorithm. Compared with the classic ALNS, our proposed ALNS reduced carbon emissions by 7.31% on average. In addition, in the same given computational time, the number of iterations or our ALNS is 3.41 times that of classic ALNS. The performance of the proposed ALNS in class C (clustered) is particularly outstanding: carbon emissions are reduced by 13.10% on average and, in the instance C1_2_2, reduced by 19.26%.It can be seen from Table6 that, in the instances of 400 customers, the proposed algorithm is better than the classic algorithm. Compared with the classic ALNS, our proposed ALNS reduced carbon emissions by 9.37% on average. In addition, in the same given computational time, the number of iterations or our ALNS is 4 times that of classic ALNS. The performance of the proposed ALNS in class C (clustered) is particularly outstanding: carbon emissions are reduced by 14.48% on average and, in the instance C1_2_2, reduced by 23.67%. In addition, our algorithm performs more prominently in the instances with vehicle type 2 (e.g., C2_4_1, C2_4_2, and R2_4_2).It can be seen from Table7 that, in the instances of 600 customers, the proposed algorithm is better than the classic algorithm. Compared with the classic ALNS, our proposed ALNS reduced carbon emissions by 7.53% on average. In addition, in the same given computational time, the number of iterations or our ALNS is 4.77 times that of classic ALNS. The performance of the proposed ALNS in class C (clustered) is particularly outstanding: carbon emissions are reduced by 11.14% on average and, in the instance C2_6_1, reduced by 14.13%. In addition, our algorithm performs more prominently in the instances with vehicle type 2.It can be seen from Table8 that, in the instance of 800 customers, the proposed algorithm is better than the classic algorithm. Compared with the classic ALNS, our proposed ALNS reduced carbon emissions by 7.91% on average. In addition, in the same given computational time, the number of iterations or our ALNS is 4.97 times that of classic ALNS. The performance of the proposed ALNS in class C (clustered) is particularly outstanding: carbon emissions are reduced by 13.38% on average and, in the instance C1_8_1, reduced by 21.20%.It can be seen from Table9 that, in the instance of 1000 customers, the proposed algorithm is better than the classic algorithm. Compared with the classic ALNS, our proposed ALNS reduced carbon emissions by 10.32% on average. In addition, in the same given computational time, the number of iterations or our ALNS is 3.56 times that of classic ALNS. The performance of the proposed ALNS in class C (clustered) is particularly outstanding: carbon emissions are reduced by 17.72% on average and, in the instance C2_10_2, reduced by 33.61%. In addition, our algorithm performs more prominently in the instances with vehicle type 2.It can be seen from Tables5–9 that the proposed ALNS algorithm is superior to the traditional ALNS algorithm in all instances. From the perspective of the geographic distribution of the customers, the proposed algorithm is particularly effective in the instances of class C (clustered) among the three classes (R, C, and RC). Take the instances of class C with 1000 customers as an example, and the algorithm improves the average accuracy by 15.09% compared with the classic ALNS algorithm. In the optimal situation, the improvement can achieve 33.61%.


### 4.4. Performance Evaluation of the Forward Load Removal Heuristic

In the proposed ALNS algorithm for GVRP, forward load removal heuristic is proposed according to the characteristics of GVRP to replace random removal heuristic in the classic ALNS algorithm. The performance of forward load removal heuristic is evaluated in Table10.Table 10
Performance evaluation of forward load removal heuristic.
InstanceALNS with the fast insertion methodOur ALNSRCE (%)QOIMCEITERMCEITERR1_10_12204.0145192201.5558120.111.29R1_10_22042.0038212042.575163−0.031.35R2_10_12644.4510462605.2615321.481.46R2_10_22372.165672352.8310200.811.80Average0.601.48C1_10_12459.7240602412.2146641.931.15C1_10_22384.0742972334.4442652.080.99C2_10_11132.0617791123.4628360.761.59C2_10_21096.1210701074.4614961.971.40Average1.691.28RC1_10_11946.0348231948.046307−0.101.31RC1_10_21834.8633181821.7043980.721.33RC2_10_11900.537031883.5812000.891.71RC2_10_21745.884501729.228080.950.46Average0.621.20Total average0.971.32As shown in Table10, compared with the classic removal heuristic proposed by [7], the proposed ALNS improved the average accuracy by 0.97%. In addition, in the same time, the number of iterations of the proposed ALNS is 1.3 times that of the classic ALNS.


### 4.5. Performance Evaluation of the Fast Insertion Method

To evaluate the performance of the fast insertion method, the traditional insertion method and the fast insertion method are compared with improved heuristics. The performance evaluation by the accuracy and iteration of the algorithm is shown in Table11.Table 11
Performance evaluation of the fast insertion method.
InstanceALNS with traditional insertion and forward load removal heuristicOur ALNSRCE (%)QOIMCEITERMCEITERR1_10_12213.3929832201.5558120.541.95R1_10_22088.6630132042.5751632.211.71R2_10_13013.912512605.26153213.566.10R2_10_22459.992012352.8310204.365.07Average5.163.71C1_10_12717.747492412.21466411.246.23C1_10_22663.507242334.44426512.355.89C2_10_11173.932651123.4628364.3010.70C2_10_21651.182291074.46149634.936.53Average15.717.34RC1_10_11977.1645121948.0463071.471.40RC1_10_21838.4831511821.7043980.911.40RC2_10_12125.902001883.58120011.406.00RC2_10_21847.642061729.228086.413.92Average5.053.18Total average8.644.74As shown in Table11, in the same time, the number of iterations of the fast insertion method is 4.74 times than that of the traditional insertion method, which means that the fast insertion method significantly saves the judgment time of the feasibility of the new route. In the same time, the fast insertion method improves the average accuracy by 8.64% compared with the traditional insertion method. In the optimal situation, the improvement can achieve 34.93%.


## 4.1. Software and Hardware Specifications

The BP algorithm was coded in C#, whose RMP and integer branch method were solved by ILOG CPLEX 12.6. Both classic ALNS and the proposed ALNS algorithms are coded in C#. All tests were executed on a laptop with a 2.80-GHz Intel Core TM i5-4200 processor using the Microsoft Windows 10 operating system using 12.00 GB of RAM. The classic ALNS and the proposed ALNS were run 10 times in each instance.


## 4.2. Experiment Instances

### 4.2.1. Introduction of the Benchmark Instance

The used GVRPTW instances are derived from Solomon VRPTW benchmark instances [45] with 100 customers and Homberger VRPTW benchmark instances [46] with 200, 400, 600, 800, and 1000 customers. Benchmark instances are divided in three classes according to the geographical distribution of the customers: R (random), C (clustered), or RC (semiclustered) [45]. For more details on the benchmark instances, refer to [45, 46].


### 4.2.2. Construction of Data Acquisition System

The standard test sets from Solomon benchmark and Homberger benchmark are used in the experiment, the coordinate, and the quantity of products, and time windows and service time of the depot and customers can be obtained.Based on the coordinate of the depot and customers, the Euclidean distancedij from node i to node j is calculated. Since the vehicle speed is fixed, the calculated distance needs to be scaled to meet the time window constraints, and the distance dij was changed to dij′=Cddij [47] and kept to four decimal places. For the quantity of products of each customer, the quantity was scaled to qi′=Cqqi. For simplicity of computation, the time windows were set as 24 h, and time windows of all customers were scaled accordingly using the coefficient 24/b0. Moreover, the time window ai,bi of customer i in Solomon instances and Homberger instances was modified as ai24/b0,bi24/b0 [47]. Furthermore, the service time was set to 1.5 h for all customers [47] in Solomon instances (100 customers), and the service time si was modified as si24/b0 in the Homberger instances (200, 400, 600, 800, and 1000 customers).According to the investigation in a logistics company in Shenyang, China, two types of vehicles for the computational tests of Homberger instance [46, 48, 49] were selected by the authors. Table 1 describes the value of parameters used in experiments of GVRP.Table 1
Value of used parameters of GVRP.
NotationDescriptionTypical valuesRCarbon emissions index parameter3.164α1Emissions unrelated to the vehicle mass1.078086α2Emissions linear to the vehicle mass0.072253W1Represents curb weight of vehicle type 12500 kgQ1Capacity of vehicle 12000 kgW2Represents curb weight of vehicle type 29000 kgQ2Capacity of vehicle 210000 kgVVehicle speed42 km/hCdDistance coefficient2N=1000.5200≤N≤8000.1N=1000CqLoad coefficient10Table2 describes the value of parameters used in experiments of ALNS.In the experiment, the running time of the heuristic algorithm is given according to the scale of the problem (Table3). Given the same running time, the performance of the algorithm is evaluated by its minimal carbon emissions.


## 4.2.1. Introduction of the Benchmark Instance

The used GVRPTW instances are derived from Solomon VRPTW benchmark instances [45] with 100 customers and Homberger VRPTW benchmark instances [46] with 200, 400, 600, 800, and 1000 customers. Benchmark instances are divided in three classes according to the geographical distribution of the customers: R (random), C (clustered), or RC (semiclustered) [45]. For more details on the benchmark instances, refer to [45, 46].


## 4.2.2. Construction of Data Acquisition System

The standard test sets from Solomon benchmark and Homberger benchmark are used in the experiment, the coordinate, and the quantity of products, and time windows and service time of the depot and customers can be obtained.Based on the coordinate of the depot and customers, the Euclidean distancedij from node i to node j is calculated. Since the vehicle speed is fixed, the calculated distance needs to be scaled to meet the time window constraints, and the distance dij was changed to dij′=Cddij [47] and kept to four decimal places. For the quantity of products of each customer, the quantity was scaled to qi′=Cqqi. For simplicity of computation, the time windows were set as 24 h, and time windows of all customers were scaled accordingly using the coefficient 24/b0. Moreover, the time window ai,bi of customer i in Solomon instances and Homberger instances was modified as ai24/b0,bi24/b0 [47]. Furthermore, the service time was set to 1.5 h for all customers [47] in Solomon instances (100 customers), and the service time si was modified as si24/b0 in the Homberger instances (200, 400, 600, 800, and 1000 customers).According to the investigation in a logistics company in Shenyang, China, two types of vehicles for the computational tests of Homberger instance [46, 48, 49] were selected by the authors. Table 1 describes the value of parameters used in experiments of GVRP.Table 1
Value of used parameters of GVRP.
NotationDescriptionTypical valuesRCarbon emissions index parameter3.164α1Emissions unrelated to the vehicle mass1.078086α2Emissions linear to the vehicle mass0.072253W1Represents curb weight of vehicle type 12500 kgQ1Capacity of vehicle 12000 kgW2Represents curb weight of vehicle type 29000 kgQ2Capacity of vehicle 210000 kgVVehicle speed42 km/hCdDistance coefficient2N=1000.5200≤N≤8000.1N=1000CqLoad coefficient10Table2 describes the value of parameters used in experiments of ALNS.In the experiment, the running time of the heuristic algorithm is given according to the scale of the problem (Table3). Given the same running time, the performance of the algorithm is evaluated by its minimal carbon emissions.


## 4.3. Computational Results of Improved ALNS Algorithm

Since the branch-and-price algorithm developed by Yu et al. [50] can solve up to 100 customer instances, for 100-customer instances, the proposed algorithm is compared with the exact solution to evaluate its performance. For instances with more than 200 customers, there is no exact solution until now, and the proposed algorithm is compared with the classic ALNS algorithm proposed by Ropke and Pisinger [7].The destroy operators of classic ALNS consist of random removal heuristic, worst removal heuristic, and Shaw removal heuristic, and the repair operators of classic ALNS consist of greedy insertion heuristic and regret insertion heuristic (regret 2 and regret 3). In addition, the objective function is changed to the total carbon emission in the routing.The proposed ALNS in the paper and its details can refer to Section3.5.Table4 shows the results of Solomon benchmark instances with 100 customers using the vehicles with Q = 2000.For the instances with 100 customers, only 4 instances can be exactly solved by the branch-and-price for GVRPTW [6]. As shown in Table 4, in the same given computational time (20 s), the proposed ALNS improves the accuracy significantly compared with the classic ALNS.Tables5–9 show the results of Homberger benchmark instances with 200, 400, 600, 800, and 1000 customers, respectively.Table 2
Value of used parameters of ALNS.
NotationDescriptionTypical valuesZNumber of removed customers8N≤10020200≤N≤1000CCoefficient in forward load removal heuristic0.5aminMinimum value of the called probability0.01amaxMaximum value of the called probability0.99ρEvaporation parameter in adaptive method0.999Table 3
Given computational time of the heuristics for the instances with different customers.
Number of customers1002004006008001000GCT (s)204080120160200GCT: given computational time; unit is seconds.Table 4
Computational results of the instances with 100 customers.
InstanceExactClassic ALNSThe proposed ALNSECECT (s)MCEGap (%)MCEGap (%)r1011127.68171.721166.283.421164.653.28r2011027.1222363.171095.676.671057.082.92c101818.5644.20850.543.91835.152.03c105597.0332.72610.562.27610.562.27Average4.072.62ECE: exact carbon emissions solved by the BP algorithm proposed by Yu et al. [50]; MCE: minimal carbon emission; gap: the relative difference between the found solution and the exact solution (i.e., Gap=MCE−ECE/ECE∗100%).Table 5
Computational results of the instances with 200 customers.
InstanceClassic ALNSOur ALNSRCE (%)QOIMCEITERMCEITERR1_2_1813.73984794.4527122.372.76R1_2_2721.851129708.6322901.832.03R2_2_11096.214331029.4114746.093.40R2_2_21004.70367924.469627.992.62Average3.432.70C1_2_1827.44945683.45637917.406.75C1_2_2727.80899587.61433319.264.82C2_2_1655.25519609.0027657.065.33C2_2_2626.34479572.0912098.662.53Average13.104.86RC1_2_1648.921024633.6922632.352.21RC1_2_2620.041175598.6419913.451.70RC2_2_1880.85404765.83188013.064.66RC2_2_2746.70447726.079562.762.14Average5.412.68Total average7.313.41ITER: number of iterations; RCEs: reduced carbon emissions (i.e.,RCE=MCE of classic ALNS−MCE of our ALNS/MCE of classic ALNS∗100%); QOI: quotient of number of iterations (i.e., QOI=ITER of our ALNS/ITER of classic ALNS).Table 6
Computational results of the instances with 400 customers.
InstanceClassic ALNSOur ALNSRCE (%)QOIMCEITERMCEITERR1_4_11952.9017551887.1055633.373.17R1_4_21793.6119451743.0544492.822.29R2_4_12686.054752319.57284513.645.99R2_4_22668.673842220.06150916.813.93Average9.163.84C1_4_12085.398721881.1961159.797.02C1_4_21810.1510291680.6435267.153.43C2_4_11818.395741388.01414423.677.22C2_4_21621.954821341.20156417.313.25Average14.485.23RC1_4_11689.1316301652.5135392.172.17RC1_4_21578.4920651561.5837511.071.82RC2_4_11918.344601823.6623894.945.20RC2_4_21846.804161667.0310589.732.54Average4.482.93Total average9.374.00Table 7
Computational results of the instances with 600 customers.
InstanceClassic ALNSOur ALNSRCE (%)QOIMCEITERMCEITERR1_6_14298.0815174224.5554511.713.59R1_6_23851.5616743725.9144083.262.63R2_6_15264.923954718.97392510.379.93R2_6_25153.552994152.05180819.436.04Average8.695.55C1_6_14328.348283948.5758778.777.10C1_6_23950.988743499.99442111.415.06C2_6_13483.274862991.08385214.137.93C2_6_23305.614482966.92136410.253.05Average11.145.79RC1_6_13483.1413373408.3940892.153.06RC1_6_23272.8017453251.0536800.662.11RC2_6_13975.134903761.9620155.364.11RC2_6_23465.404443366.2111672.862.63Average2.762.98Total average7.534.77Table 8
Computational results of the instances with 800 customers.
InstanceClassic ALNSOur ALNSRCE (%)QOIMCEITERMCEITERR1_8_17417.0715737210.6860802.783.86R1_8_26868.4918376627.2946883.512.55R2_8_19113.605328206.3036969.966.95R2_8_27803.534606995.30176710.363.84Average6.654.30C1_8_19107.378967177.04550321.206.14C1_8_26824.469426348.2241946.984.45C2_8_15550.204995032.8838939.327.81C2_8_25484.324204606.45279816.016.66Average13.386.26RC1_8_16245.8413966155.8552401.443.75RC1_8_25954.1016045779.1539832.942.48RC2_8_16894.014336575.1228634.636.61RC2_8_26451.263996074.7318355.844.59Average3.714.36Total average7.914.97Table 9
Computational results of the instances with 1000 customers.
InstanceClassic ALNSOur ALNSRCE %QOIMCEITERMCEITERR1_10_12217.4931642201.5558120.721.84R1_10_22095.0733602042.5751632.511.54R2_10_12953.183172605.26153211.784.84R2_10_22717.363102352.83102013.413.29Average7.112.88C1_10_12741.778162412.21466412.025.72C1_10_22781.218442334.44426516.065.05C2_10_11237.293991123.4628369.207.11C2_10_21618.343591074.46149633.614.16Average17.725.51RC1_10_11977.3442271948.0463071.481.49RC1_10_21831.0231751821.7043980.511.39RC2_10_12067.783261883.5812008.913.68RC2_10_22002.163051729.2280813.632.65Average6.132.30Total average10.323.56It can be seen from Table5 that, in the instances of 200 customers, the proposed algorithm is better than the classic algorithm. Compared with the classic ALNS, our proposed ALNS reduced carbon emissions by 7.31% on average. In addition, in the same given computational time, the number of iterations or our ALNS is 3.41 times that of classic ALNS. The performance of the proposed ALNS in class C (clustered) is particularly outstanding: carbon emissions are reduced by 13.10% on average and, in the instance C1_2_2, reduced by 19.26%.It can be seen from Table6 that, in the instances of 400 customers, the proposed algorithm is better than the classic algorithm. Compared with the classic ALNS, our proposed ALNS reduced carbon emissions by 9.37% on average. In addition, in the same given computational time, the number of iterations or our ALNS is 4 times that of classic ALNS. The performance of the proposed ALNS in class C (clustered) is particularly outstanding: carbon emissions are reduced by 14.48% on average and, in the instance C1_2_2, reduced by 23.67%. In addition, our algorithm performs more prominently in the instances with vehicle type 2 (e.g., C2_4_1, C2_4_2, and R2_4_2).It can be seen from Table7 that, in the instances of 600 customers, the proposed algorithm is better than the classic algorithm. Compared with the classic ALNS, our proposed ALNS reduced carbon emissions by 7.53% on average. In addition, in the same given computational time, the number of iterations or our ALNS is 4.77 times that of classic ALNS. The performance of the proposed ALNS in class C (clustered) is particularly outstanding: carbon emissions are reduced by 11.14% on average and, in the instance C2_6_1, reduced by 14.13%. In addition, our algorithm performs more prominently in the instances with vehicle type 2.It can be seen from Table8 that, in the instance of 800 customers, the proposed algorithm is better than the classic algorithm. Compared with the classic ALNS, our proposed ALNS reduced carbon emissions by 7.91% on average. In addition, in the same given computational time, the number of iterations or our ALNS is 4.97 times that of classic ALNS. The performance of the proposed ALNS in class C (clustered) is particularly outstanding: carbon emissions are reduced by 13.38% on average and, in the instance C1_8_1, reduced by 21.20%.It can be seen from Table9 that, in the instance of 1000 customers, the proposed algorithm is better than the classic algorithm. Compared with the classic ALNS, our proposed ALNS reduced carbon emissions by 10.32% on average. In addition, in the same given computational time, the number of iterations or our ALNS is 3.56 times that of classic ALNS. The performance of the proposed ALNS in class C (clustered) is particularly outstanding: carbon emissions are reduced by 17.72% on average and, in the instance C2_10_2, reduced by 33.61%. In addition, our algorithm performs more prominently in the instances with vehicle type 2.It can be seen from Tables5–9 that the proposed ALNS algorithm is superior to the traditional ALNS algorithm in all instances. From the perspective of the geographic distribution of the customers, the proposed algorithm is particularly effective in the instances of class C (clustered) among the three classes (R, C, and RC). Take the instances of class C with 1000 customers as an example, and the algorithm improves the average accuracy by 15.09% compared with the classic ALNS algorithm. In the optimal situation, the improvement can achieve 33.61%.


## 4.4. Performance Evaluation of the Forward Load Removal Heuristic

In the proposed ALNS algorithm for GVRP, forward load removal heuristic is proposed according to the characteristics of GVRP to replace random removal heuristic in the classic ALNS algorithm. The performance of forward load removal heuristic is evaluated in Table10.Table 10
Performance evaluation of forward load removal heuristic.
InstanceALNS with the fast insertion methodOur ALNSRCE (%)QOIMCEITERMCEITERR1_10_12204.0145192201.5558120.111.29R1_10_22042.0038212042.575163−0.031.35R2_10_12644.4510462605.2615321.481.46R2_10_22372.165672352.8310200.811.80Average0.601.48C1_10_12459.7240602412.2146641.931.15C1_10_22384.0742972334.4442652.080.99C2_10_11132.0617791123.4628360.761.59C2_10_21096.1210701074.4614961.971.40Average1.691.28RC1_10_11946.0348231948.046307−0.101.31RC1_10_21834.8633181821.7043980.721.33RC2_10_11900.537031883.5812000.891.71RC2_10_21745.884501729.228080.950.46Average0.621.20Total average0.971.32As shown in Table10, compared with the classic removal heuristic proposed by [7], the proposed ALNS improved the average accuracy by 0.97%. In addition, in the same time, the number of iterations of the proposed ALNS is 1.3 times that of the classic ALNS.


## 4.5. Performance Evaluation of the Fast Insertion Method

To evaluate the performance of the fast insertion method, the traditional insertion method and the fast insertion method are compared with improved heuristics. The performance evaluation by the accuracy and iteration of the algorithm is shown in Table11.Table 11
Performance evaluation of the fast insertion method.
InstanceALNS with traditional insertion and forward load removal heuristicOur ALNSRCE (%)QOIMCEITERMCEITERR1_10_12213.3929832201.5558120.541.95R1_10_22088.6630132042.5751632.211.71R2_10_13013.912512605.26153213.566.10R2_10_22459.992012352.8310204.365.07Average5.163.71C1_10_12717.747492412.21466411.246.23C1_10_22663.507242334.44426512.355.89C2_10_11173.932651123.4628364.3010.70C2_10_21651.182291074.46149634.936.53Average15.717.34RC1_10_11977.1645121948.0463071.471.40RC1_10_21838.4831511821.7043980.911.40RC2_10_12125.902001883.58120011.406.00RC2_10_21847.642061729.228086.413.92Average5.053.18Total average8.644.74As shown in Table11, in the same time, the number of iterations of the fast insertion method is 4.74 times than that of the traditional insertion method, which means that the fast insertion method significantly saves the judgment time of the feasibility of the new route. In the same time, the fast insertion method improves the average accuracy by 8.64% compared with the traditional insertion method. In the optimal situation, the improvement can achieve 34.93%.


## 5. Conclusion

In this paper, an adaptive large neighborhood search (ALNS) algorithm is proposed to solve large-scale instances of GVRP. In the destroy operators, a new removal heuristic applying to the characteristics of GVRP is proposed. The heuristic can quickly remove customers who bring a large amount of carbon emissions with pertinence, and these customers may be arranged more properly in future repair operators. In the repair operators, a fast insertion method, suitable for all variants of VRPTW, is developed. In the fast insertion method, the feasibility of a new route is judged by checking the constraints of partial customers after the inserted customer, instead of checking the constraints of all customers. Thus, the computational time of ALNS is greatly saved. Computational experiments were performed on Solomon benchmark with 100 customers and Homberger benchmark instances with up to 1000 customers. Given the same computational time, the proposed ALNS improves the average accuracy by 8.49% compared with the classic ALNS. In the optimal situation, the improvement can achieve 33.61%.



---

*Source: 8210630-2020-10-29.xml*