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

Abstract

Due to huge amount of greenhouse gases emission (such as CO₂), 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%.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 Section 2, 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=\mathrm{0,1},\dots ,n+1$. Nodes 0 and $n+1$ denote the start and the end depot, respectively. Let $V_0=V\backslash 0,n+1$ be the set of customers. For customer $i$, the quantity $q_i$ of products is to be served by the vehicle within the time window $a_i,b_i$. In addition, the service time s_i is required. $E$ represents the set of arcs. For each arc $i,j\in E$, let $d_{ij}$ 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 $a_i,b_i$, and the vehicle will wait if it arrives before the ready time $a_i$; and (5) each vehicle needs to return to the depot before the due time $b_{n+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,\dots ,n+1$.

  V₀: set of customers V₀ = V\{0, n + 1}.

  i, j: nodes i, j ∈ V.

  d_ij: distance between node i and node j.

  f_ij: vehicle’s payload through arc (i, j).

  q_i: quantity of products at the customer i.

  [a_i, b_i]: time window when service may start at node i, where a_i and b_i represent the earliest and latest time, respectively.

  s_i: service time at node i, and the service time at depot (i.e., node 0 and n + 1) is set to 0.

  e_ij: carbon emissions between node i and node j.

  u_i: time of vehicle arriving at customer i.

  d_i: load of vehicle leaving from customer i.

  $p$: a feasible route visited by a vehicle, $p=\mathrm{0,1},\dots ,l,n+1$, where l is number of customers in the route $p$.

  Q: capacity of vehicle.

  $w$: curb weight of vehicles.

  x_ij: decision variables; if the vehicle is driving from node i to node j, then x_ij = 1; otherwise, x_ij = 0.

  t_ij: travel time between node i and node j.

2.2.2. Mathematical Formulation

The carbon emissions calculation formula is used [40]:$$\begin{array}{}\text{(1)}& e_{ij}=r\frac{{\alpha }_l\times d_{ij}}{v}+\frac{{\alpha }_2\times w+f_{ij}\times d_{ij}}{v},\end{array}$$where r represents carbon emissions index parameter, by which the fuel consumption on arc $i,j$ can be transformed into carbon emissions; ${\alpha }_1$ represents the emissions unrelated to the vehicle mass; ${\alpha }_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:$$\begin{array}{}\text{(2)}& \min {\displaystyle \sum _{i,j\in E}{e}_{ij}{x}_{ij}},\text{(3)}& \text{s.t.\hspace{0.17em}}{\displaystyle \sum _{i,j\in E}{x}_{ij}}=1,\hspace{1em}\forall i\in V_0,\text{(4)}& {\displaystyle \sum _{0,i\in E}{x}_{0i}}={\displaystyle \sum _{i,n+1\in E}{x}_{i,n+1}},\hspace{1em}\forall i\in V_0,\text{(5)}& {\displaystyle \sum _{j,i\in E}{x}_{ji}}={\displaystyle \sum _{i,j\in E}{x}_{ij}},\hspace{1em}\forall i\in V_0,\text{(6)}& u_j\ge u_i+s_i+t_{ij}{x}_{ij},\hspace{1em}\forall i,j\in E,\text{(7)}& a_i\le u_i\le b_i,\hspace{1em}\forall i\in V_0,\text{(8)}& d_j\ge d_i+q_j{x}_{ij},\hspace{1em}\forall i,j\in E,\text{(9)}& 0\le d_i\le Q,\hspace{1em}\forall i\in V_0,\text{(10)}& u_i\ge 0,\hspace{1em}\forall i\in V,\text{(11)}& x_{ij}\in \mathrm{0,1}\cdot \forall i,j\in E.\end{array}$$

Equation (2) is the objective function, where e_ij is related to the vehicle’s current load f_ij and the distance d_ij between node i and node j. Therefore, equation (2) is nonlinear with decision variables x_ij. 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 Algorithm 1.

Algorithm 1: ALNS ().

  Input: an initial solution (x), the number of removed customers $z\in 0,\dots ,n$.

  Output: best solution (x_best).

  $x_{\text{best}}⟵x$.

  repeat

   $x^{\prime }⟵x$.

   (1) Perform destroy operators on $x^{\prime }$. See Section 3.2.

   (2) Perform repair operators on $x^{\prime }$. See Section 3.3.

   if ($f{x}^{\prime }<fx$) or (simulated annealing criteria is met) then

    $x⟵x^{\prime }$.

    if$f{x}^{\prime }<f{x}_{\text{best}}$then

     $x_{\text{best}}⟵x^{\prime }$.

    end if

   end if

   Perform adaptive method. See Section 3.4.

  until termination condition is met.

  Output x_best.

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^{\prime }$ is updated by removing z customers from $x^{\prime }$. 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 Algorithm 2.

Algorithm 2: Forward load removal heuristic ().

  Input: a feasible solution (x), the number of removed customers (z), parameter $c\in \mathrm{0,1}$.

  Output: set of removed customers.

  $D⟵\varnothing$. D is the set of removed customers.

  $H⟵\varnothing$. H is used to record the array of customers sorted by descending q_i.

  $M⟵0$. M records the order of selected customer.

  while z > 0 do

   Randomly choose a feasible route $p=\mathrm{0,1},\dots ,l,n+1\in x$ with more than 1 customers i.e., l > 1.

   $H$ ⟵ sorting p by descending q_i.

   $M⟵\lfloor c^{\ast }l\rfloor$.

  (1) for$j=1\hspace{0.17em}\text{to}\hspace{0.17em}M-1$do

   (2) if$q_j>HM.\hspace{0.17em}\text{demand}$do

   (3) Remove j from p.

     $D⟵D\cup j$.

     $z⟵z-1$.

    end if

   end for

  end while

  Output D.

In step (1), M is the order of selected customer.

In step (2), the customer before M in the route $p$ is the customer served earlier, i.e., $j=1\hspace{0.17em}\text{to}\hspace{0.17em}M-1$.

In step (3), if the customer’s demand is higher than the demand of $HM$, i.e., $q_j>HM.\hspace{0.17em}\text{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^{\prime }$ is updated by reinserting z customers (removed in destroy operators) into $x^{\prime }$. 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., Situations 1 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 j

  r: inserted customer r

  $p=\mathrm{0,1},\dots ,j,j+1,\dots ,l,n+1$: a feasible route;

  $p_{\text{ins}}=\mathrm{0,1},\dots ,j,r,j+1,\dots ,l,n+1$: a route produced by inserting customer r into p;

  dt_j: the departure time of customer j

  $DT=d{t}_0,d{t}_1,\dots ,d{t}_j,d{t}_{j+1},\dots ,d{t}_l,d{t}_{n+1}$: departure time array of the customers in route p;

  $D{T}_{ins}=d{t}_0,d{t}_1,\dots ,d{t}_j,d{t}_r,d{t}_{j+1},\dots ,d{t}_l,d{t}_{n+1}$: departure time array of the customers in route p_ins.

Situation 1.

If there exists a customer k after the inserted customer r in the new route (i.e., $p_{\text{ins}}=\mathrm{0,1},\dots ,j,r,j+1,\dots ,k,,l,n+1$), with $d{t}_k\hspace{0.17em}\text{in}\hspace{0.17em}D{T}_{\text{ins}}$ equal to $d{t}_k\text{\hspace{0.17em}in\hspace{0.17em}\hspace{0.17em}}DT$, then the new route pins is a feasible route.

Example 1.

In a feasible route $p=\mathrm{0,1,2,3,4},n+1$, $DT=\mathrm{0,1,4,5,6,7}$ is the departure time array of the customers in p. For customer 9, $t_{19}=t_{92}=1$, $a_9,b_9=\mathrm{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 $p_{\text{ins}}=\mathrm{0,1,9,2,3,4},n+1$, if $d{t}_2\hspace{0.17em}\text{in}\hspace{0.17em}D{T}_{\text{ins}}=4$ is equal to $d{t}_2\hspace{0.17em}\text{in}\hspace{0.17em}DT$, then p_ins must be a feasible route. Thus, it is unnecessary to check the constraints of customers 3 and 4.

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

If r 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 route $p=\mathrm{0,1,2},n+1$, then customer 9 cannot be inserted in $p_{\text{ins}}=\mathrm{0,1,2,3},n+1$ generated by p.

According to Situation 2, the feasibility of the route can be prejudged before inserting r.

The program of the fast insertion method is shown in Algorithm 3.

Step (1) initializes the parameters.

Step (2) checks whether the route $p$ 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., $q_p+q_r\le Q$, and the time constraint of customer r, i.e., $d{t}_j+t_{jr}\le b_r$ 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 customer r.

Step (4) checks if there exists a customer k after the inserted customer r in the new route pins with $d{t}_k\hspace{0.17em}\text{in}\hspace{0.17em}D{T}_{\text{ins}}$ equal to $d{t}_k\hspace{0.17em}\text{in}\hspace{0.17em}DT$ (see situation 1).

Step (4-1) checks the time constraint of customer k, i.e., $d{t}_{k-1}+t_{k-1,k}\le b_k$, 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 customer k.

Step (4-1-2) judges whether $d{t}_k\hspace{0.17em}\text{in}\hspace{0.17em}D{T}_{\text{ins}}$ is equal to $d{t}_k\hspace{0.17em}\text{in}\hspace{0.17em}D{T}_{\text{ins}}$. 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 customer r.

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=\mathrm{0,1},\dots ,j,j+1,\dots ,l,n+1$ and customer r.

  Output: feasibility (flag) of the insertion

  q_p is the total load of the route p.

  $\text{flag}⟵\text{false}$.

   (1) if p is generated by a route which r cannot be inserted into then

   (2) Prejudge the feasibility of the insertion (i.e., situation 2).

   Output flag.

   else

    if ($q_p+q_r\le Q$) and ($d{t}_j+t_{jr}\le b_r$) then

   (3) $d{t}_r⟵\max d{t}_j+t_{jr},a_r+s_r$

    (3-1) for$k=j+1\hspace{0.17em}\text{to}\hspace{0.17em}l$do

   (4) if$k=j+1$then

    $d{t}_{k-1}⟵d{t}_r$.

    $t_{k-1,k}⟵t_{r,j+1}$.

   end if

   if$d{t}_{k-1}+t_{k-1,k}\le b_k$then

    (4-1) $d{t}_k⟵\max d{t}_{k-1}+t_{k-1,k},a_k+s_k$.

     (4-1-1) if$d{t}_k\text{in}\hspace{0.17em}D{T}_{\text{ins}}=d{t}_k\text{in}\hspace{0.17em}DT$then

     (4-1-2) $\text{flag}⟵\text{true}$. Judge whether route is feasible using partial customers after r (i.e.,situation 1).

      break

     end if

    else

     break

    end if

    if$k=l$then

    (4-2) $\text{flag}⟵\text{true}$.

     end if

    end for

   end if

  end if

  Output 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:$$\begin{array}{}\text{(12)}& a_i^{t+1}=\rho a_i^t+\frac{f{x}^t-f{x}^{t+1}}{f{x}^t},\end{array}$$where $a_i^t$ is the probability of the algorithm i being called in the t^th iteration, $\rho \in \mathrm{0,1}$ is the evaporation coefficient, and $f{x}^t$ is the objective value of the solution x^t 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 $a_{\min },a_{\max }$.

3.5. Flowchart of ALNS for GVRP

The flowchart of ALNS for GVRP is shown in Figure 1.

[figure omitted; refer to PDF]

As shown in Figure 1, 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 (Section 3.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 (Section 3.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^{\prime }$) generated in this iteration is better than (carbon emission is smaller than) the best solution (x_best), update the best solution, i.e., $x_{\text{best}}\leftarrow x^{\prime }$.

If the solution ($x^{\prime }$) 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\leftarrow x^{\prime }$.

The probability of the heuristics being called in the destroy operators and repair operators keeps changing, according to its performance in the adaptive method (Section 3.4).

When the given running time is reached, the best solution (x_best) 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 distance d_ij 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 $d_{ij}$ was changed to $d_{ij}^{\prime }=C_d{d}_{ij}$ [47] and kept to four decimal places. For the quantity of products of each customer, the quantity was scaled to $q_i^{\prime }=C_q{q}_i$. 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/b_0$. Moreover, the time window $a_i,b_i$ of customer i in Solomon instances and Homberger instances was modified as $a_{i}24/b_0,b_{i}24/b_0$ [47]. Furthermore, the service time was set to 1.5 h for all customers [47] in Solomon instances (100 customers), and the service time s_i was modified as $s_{i}24/b_0$ 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.

Notation	Description	Typical values
R	Carbon emissions index parameter	3.164
${\alpha }_1$	Emissions unrelated to the vehicle mass	1.078086
${\alpha }_2$	Emissions linear to the vehicle mass	0.072253
W1	Represents curb weight of vehicle type 1	2500 kg
Q1	Capacity of vehicle 1	2000 kg
W2	Represents curb weight of vehicle type 2	9000 kg
Q2	Capacity of vehicle 2	10000 kg
V	Vehicle speed	42 km/h
Cd	Distance coefficient	$2N=100$
		$0.5200\le N\le 800$
		$0.1N=1000$
Cq	Load coefficient	10

Table 2 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 (Table 3). 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 Section 3.5.

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

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

Notation	Description	Typical values
Z	Number of removed customers	$8N\le 100$
		$20200\le N\le 1000$
C	Coefficient in forward load removal heuristic	0.5
amin	Minimum value of the called probability	0.01
amax	Maximum value of the called probability	0.99
$\rho$	Evaporation parameter in adaptive method	0.999

Table 3

Given computational time of the heuristics for the instances with different customers.

Number of customers	100	200	400	600	800	1000
GCT (s)	20	40	80	120	160	200

GCT: given computational time; unit is seconds.

Table 4

Computational results of the instances with 100 customers.

Instance	Exact		Classic ALNS		The proposed ALNS
	ECE	CT (s)	MCE	Gap (%)	MCE	Gap (%)
r101	1127.68	171.72	1166.28	3.42	1164.65	3.28
r201	1027.12	22363.17	1095.67	6.67	1057.08	2.92
c101	818.56	44.20	850.54	3.91	835.15	2.03
c105	597.03	32.72	610.56	2.27	610.56	2.27
Average				4.07		2.62

ECE: 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., $\text{Gap}=\text{MCE}-\text{ECE}/\text{ECE}\ast 100\%$).

Table 5

Computational results of the instances with 200 customers.

Instance	Classic ALNS		Our ALNS		RCE (%)	QOI
	MCE	ITER	MCE	ITER
R1_2_1	813.73	984	794.45	2712	2.37	2.76
R1_2_2	721.85	1129	708.63	2290	1.83	2.03
R2_2_1	1096.21	433	1029.41	1474	6.09	3.40
R2_2_2	1004.70	367	924.46	962	7.99	2.62
Average					3.43	2.70
C1_2_1	827.44	945	683.45	6379	17.40	6.75
C1_2_2	727.80	899	587.61	4333	19.26	4.82
C2_2_1	655.25	519	609.00	2765	7.06	5.33
C2_2_2	626.34	479	572.09	1209	8.66	2.53
Average					13.10	4.86
RC1_2_1	648.92	1024	633.69	2263	2.35	2.21
RC1_2_2	620.04	1175	598.64	1991	3.45	1.70
RC2_2_1	880.85	404	765.83	1880	13.06	4.66
RC2_2_2	746.70	447	726.07	956	2.76	2.14
Average					5.41	2.68
Total average					7.31	3.41

ITER: number of iterations; RCEs: reduced carbon emissions (i.e., $\text{RCE}=\text{MCE\hspace{0.17em}of\hspace{0.17em}classic\hspace{0.17em}ALNS}-\text{MCE\hspace{0.17em}of\hspace{0.17em}our\hspace{0.17em}ALNS}/\text{MCE\hspace{0.17em}of\hspace{0.17em}classic\hspace{0.17em}ALNS}\ast 100\%$); QOI: quotient of number of iterations (i.e., $\text{QOI\hspace{0.17em}}=\text{ITER\hspace{0.17em}of\hspace{0.17em}our\hspace{0.17em}ALNS}/\text{ITER\hspace{0.17em}of\hspace{0.17em}classic\hspace{0.17em}ALNS}$).

Table 6

Computational results of the instances with 400 customers.

Instance	Classic ALNS		Our ALNS		RCE (%)	QOI
	MCE	ITER	MCE	ITER
R1_4_1	1952.90	1755	1887.10	5563	3.37	3.17
R1_4_2	1793.61	1945	1743.05	4449	2.82	2.29
R2_4_1	2686.05	475	2319.57	2845	13.64	5.99
R2_4_2	2668.67	384	2220.06	1509	16.81	3.93
Average					9.16	3.84
C1_4_1	2085.39	872	1881.19	6115	9.79	7.02
C1_4_2	1810.15	1029	1680.64	3526	7.15	3.43
C2_4_1	1818.39	574	1388.01	4144	23.67	7.22
C2_4_2	1621.95	482	1341.20	1564	17.31	3.25
Average					14.48	5.23
RC1_4_1	1689.13	1630	1652.51	3539	2.17	2.17
RC1_4_2	1578.49	2065	1561.58	3751	1.07	1.82
RC2_4_1	1918.34	460	1823.66	2389	4.94	5.20
RC2_4_2	1846.80	416	1667.03	1058	9.73	2.54
Average					4.48	2.93
Total average					9.37	4.00

Table 7

Computational results of the instances with 600 customers.

Instance	Classic ALNS		Our ALNS		RCE (%)	QOI
	MCE	ITER	MCE	ITER
R1_6_1	4298.08	1517	4224.55	5451	1.71	3.59
R1_6_2	3851.56	1674	3725.91	4408	3.26	2.63
R2_6_1	5264.92	395	4718.97	3925	10.37	9.93
R2_6_2	5153.55	299	4152.05	1808	19.43	6.04
Average					8.69	5.55
C1_6_1	4328.34	828	3948.57	5877	8.77	7.10
C1_6_2	3950.98	874	3499.99	4421	11.41	5.06
C2_6_1	3483.27	486	2991.08	3852	14.13	7.93
C2_6_2	3305.61	448	2966.92	1364	10.25	3.05
Average					11.14	5.79
RC1_6_1	3483.14	1337	3408.39	4089	2.15	3.06
RC1_6_2	3272.80	1745	3251.05	3680	0.66	2.11
RC2_6_1	3975.13	490	3761.96	2015	5.36	4.11
RC2_6_2	3465.40	444	3366.21	1167	2.86	2.63
Average					2.76	2.98
Total average					7.53	4.77

Table 8

Computational results of the instances with 800 customers.

Instance	Classic ALNS		Our ALNS		RCE (%)	QOI
	MCE	ITER	MCE	ITER
R1_8_1	7417.07	1573	7210.68	6080	2.78	3.86
R1_8_2	6868.49	1837	6627.29	4688	3.51	2.55
R2_8_1	9113.60	532	8206.30	3696	9.96	6.95
R2_8_2	7803.53	460	6995.30	1767	10.36	3.84
Average					6.65	4.30
C1_8_1	9107.37	896	7177.04	5503	21.20	6.14
C1_8_2	6824.46	942	6348.22	4194	6.98	4.45
C2_8_1	5550.20	499	5032.88	3893	9.32	7.81
C2_8_2	5484.32	420	4606.45	2798	16.01	6.66
Average					13.38	6.26
RC1_8_1	6245.84	1396	6155.85	5240	1.44	3.75
RC1_8_2	5954.10	1604	5779.15	3983	2.94	2.48
RC2_8_1	6894.01	433	6575.12	2863	4.63	6.61
RC2_8_2	6451.26	399	6074.73	1835	5.84	4.59
Average					3.71	4.36
Total average					7.91	4.97

Table 9

Computational results of the instances with 1000 customers.

Instance	Classic ALNS		Our ALNS		RCE %	QOI
	MCE	ITER	MCE	ITER
R1_10_1	2217.49	3164	2201.55	5812	0.72	1.84
R1_10_2	2095.07	3360	2042.57	5163	2.51	1.54
R2_10_1	2953.18	317	2605.26	1532	11.78	4.84
R2_10_2	2717.36	310	2352.83	1020	13.41	3.29
Average					7.11	2.88
C1_10_1	2741.77	816	2412.21	4664	12.02	5.72
C1_10_2	2781.21	844	2334.44	4265	16.06	5.05
C2_10_1	1237.29	399	1123.46	2836	9.20	7.11
C2_10_2	1618.34	359	1074.46	1496	33.61	4.16
Average					17.72	5.51
RC1_10_1	1977.34	4227	1948.04	6307	1.48	1.49
RC1_10_2	1831.02	3175	1821.70	4398	0.51	1.39
RC2_10_1	2067.78	326	1883.58	1200	8.91	3.68
RC2_10_2	2002.16	305	1729.22	808	13.63	2.65
Average					6.13	2.30
Total average					10.32	3.56

It can be seen from Table 5 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 Table 6 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 Table 7 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 Table 8 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 Table 9 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 Tables 5–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 Table 10.

Table 10

Performance evaluation of forward load removal heuristic.

Instance	ALNS with the fast insertion method		Our ALNS		RCE (%)	QOI
	MCE	ITER	MCE	ITER
R1_10_1	2204.01	4519	2201.55	5812	0.11	1.29
R1_10_2	2042.00	3821	2042.57	5163	−0.03	1.35
R2_10_1	2644.45	1046	2605.26	1532	1.48	1.46
R2_10_2	2372.16	567	2352.83	1020	0.81	1.80
Average					0.60	1.48
C1_10_1	2459.72	4060	2412.21	4664	1.93	1.15
C1_10_2	2384.07	4297	2334.44	4265	2.08	0.99
C2_10_1	1132.06	1779	1123.46	2836	0.76	1.59
C2_10_2	1096.12	1070	1074.46	1496	1.97	1.40
Average					1.69	1.28
RC1_10_1	1946.03	4823	1948.04	6307	−0.10	1.31
RC1_10_2	1834.86	3318	1821.70	4398	0.72	1.33
RC2_10_1	1900.53	703	1883.58	1200	0.89	1.71
RC2_10_2	1745.88	450	1729.22	808	0.95	0.46
Average					0.62	1.20
Total average					0.97	1.32

As shown in Table 10, 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 Table 11.

Table 11

Performance evaluation of the fast insertion method.

Instance	ALNS with traditional insertion and forward load removal heuristic		Our ALNS		RCE (%)	QOI
	MCE	ITER	MCE	ITER
R1_10_1	2213.39	2983	2201.55	5812	0.54	1.95
R1_10_2	2088.66	3013	2042.57	5163	2.21	1.71
R2_10_1	3013.91	251	2605.26	1532	13.56	6.10
R2_10_2	2459.99	201	2352.83	1020	4.36	5.07
Average					5.16	3.71
C1_10_1	2717.74	749	2412.21	4664	11.24	6.23
C1_10_2	2663.50	724	2334.44	4265	12.35	5.89
C2_10_1	1173.93	265	1123.46	2836	4.30	10.70
C2_10_2	1651.18	229	1074.46	1496	34.93	6.53
Average					15.71	7.34
RC1_10_1	1977.16	4512	1948.04	6307	1.47	1.40
RC1_10_2	1838.48	3151	1821.70	4398	0.91	1.40
RC2_10_1	2125.90	200	1883.58	1200	11.40	6.00
RC2_10_2	1847.64	206	1729.22	808	6.41	3.92
Average					5.05	3.18
Total average					8.64	4.74

As shown in Table 11, 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%.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (71571037, 71831006, 71420107028, 71601089, and 71620107003) and the Fundamental Research Funds for the Central Universities (N170405005).

References

[1] Y. Zhou, J.-B. Sheu, J. Wang, "Robustness assessment of urban road network with consideration of multiple hazard events," Risk Analysis, vol. 37 no. 8, pp. 1477-1494, DOI: 10.1111/risa.12802, 2017.

[2] Y. Park, J. Chae, "A review of the solution approaches used in recent G-VRP (green vehicle routing problem)," International Journal of Advanced Logistics, vol. 3 no. 1-2, pp. 27-37, DOI: 10.1080/2287108x.2014.956976, 2014.

[3] T. Bektaş, E. Demir, G. Laporte, "Green vehicle routing," Green Transportation Logistics, pp. 243-265, 2016.

[4] T. Bektaş, G. Laporte, "The pollution-routing problem," Transportation Research Part B: Methodological, vol. 45 no. 8, pp. 1232-1250, DOI: 10.1016/j.trb.2011.02.004, 2011.

[5] J. Qian, R. Eglese, "Finding least fuel emission paths in a network with time-varying speeds," Networks, vol. 63 no. 1, pp. 96-106, DOI: 10.1002/net.21524, 2014.

[6] Y. Yu, S. Wang, J. Wang, M. Huang, "A branch-and-price algorithm for the heterogeneous fleet green vehicle routing problem with time windows," Transportation Research Part B: Methodological, vol. 122, pp. 511-527, DOI: 10.1016/j.trb.2019.03.009, 2019.

[7] S. Ropke, D. Pisinger, "An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows," Transportation Science, vol. 40 no. 4, pp. 455-472, DOI: 10.1287/trsc.1050.0135, 2006.

[8] E. Demir, T. Bektaş, G. Laporte, "An adaptive large neighborhood search heuristic for the pollution-routing problem," European Journal of Operational Research, vol. 223 no. 2, pp. 346-359, DOI: 10.1016/j.ejor.2012.06.044, 2012.

[9] R. Fukasawa, Q. He, Y. Song, "A branch-cut-and-price algorithm for the energy minimization vehicle routing problem," Transportation Science, vol. 50 no. 1, pp. 23-34, DOI: 10.1287/trsc.2015.0593, 2016.

[10] S. Majidi, S.-M. Hosseini-Motlagh, J. Ignatius, "Adaptive large neighborhood search heuristic for pollution-routing problem with simultaneous pickup and delivery," Soft Computing, vol. 22 no. 9, pp. 2851-2865, DOI: 10.1007/s00500-017-2535-5, 2018.

[11] H. K. E. Abad, B. Vahdani, M. Sharifi, F. Etebari, "A bi-objective model for pickup and delivery pollution-routing problem with integration and consolidation shipments in cross-docking system," Journal of Cleaner Production, vol. 193, pp. 784-801, DOI: 10.1016/j.jclepro.2018.05.046, 2018.

[12] Z. Guo, D. Zhang, H. Liu, Z. He, L. Shi, "Green transportation scheduling with pickup time and transport mode selections using a novel multi-objective memetic optimization approach," Transportation Research Part D: Transport and Environment, vol. 60, pp. 137-152, DOI: 10.1016/j.trd.2016.02.003, 2018.

[13] I. Kara, B. Y. Kara, M. K. Yetis, "Energy minimizing vehicle routing problem," Proceedings of the 1st Conference on Combinatorial Optimization and Applications, .

[14] R. Fukasawa, Q. He, Y. Song, "A disjunctive convex programming approach to the pollution-routing problem," Transportation Research Part B: Methodological, vol. 94, pp. 61-79, DOI: 10.1016/j.trb.2016.09.006, 2016.

[15] C. Lin, K. L. Choy, G. T. S. Ho, S. H. Chung, H. Y. Lam, "Survey of green vehicle routing problem: past and future trends," Expert Systems with Applications, vol. 41 no. 4, pp. 1118-1138, DOI: 10.1016/j.eswa.2013.07.107, 2014.

[16] B. Sawik, J. Faulin, E. Pérez-Bernabeu, "A multicriteria analysis for the green VRP: a case discussion for the distribution problem of a Spanish retailer," Transportation Research Procedia, vol. 22, pp. 305-313, DOI: 10.1016/j.trpro.2017.03.037, 2017.

[17] A. Franceschetti, D. Honhon, T. V. Woensel, T. Bektaş, G. Laporte, "The time-dependent pollution-routing problem," Transportation Research Part B: Methodological, vol. 56, pp. 265-293, DOI: 10.1016/j.trb.2013.08.008, 2013.

[18] Y. Yu, J. Tang, J. Li, W. Sun, J. Wang, "Reducing carbon emission of pickup and delivery using integrated scheduling," Transportation Research Part D: Transport and Environment, vol. 47, pp. 237-250, DOI: 10.1016/j.trd.2016.05.011, 2016.

[19] Y. Xiao, Q. Zhao, I. Kaku, Y. Xu, "Development of a fuel consumption optimization model for the capacitated vehicle routing problem," Computers & Operations Research, vol. 39 no. 7, pp. 1419-1431, DOI: 10.1016/j.cor.2011.08.013, 2012.

[20] S. Úbeda, J. Faulin, A. Serrano, F. J. Arcelus, "Solving the green capacitated vehicle routing problem using a tabu search algorithm," Lecture Notes in Management Science, vol. 6 no. 1, pp. 141-149, 2014.

[21] S. Zhang, C. K. M. Lee, K. L. Choy, W. Ho, W. H. Ip, "Design and development of a hybrid artificial bee colony algorithm for the environmental vehicle routing problem," Transportation Research Part D: Transport and Environment, vol. 31, pp. 85-99, DOI: 10.1016/j.trd.2014.05.015, 2014.

[22] A. Montoya, C. Guéret, J. Mendoza, J. Villegasd, "A multi-space sampling heuristic for the green vehicle routing problem," Transportation Research Part C: Emerging Technologies, vol. 70, pp. 113-128, 2015.

[23] S. Ene, I. Küçükoğlu, A. Aksoy, N. Öztürk, "A hybrid metaheuristic algorithm for the green vehicle routing problem with a heterogeneous fleet," International Journal of Vehicle Design, vol. 71 no. 1–4, pp. 75-102, DOI: 10.1504/ijvd.2016.078771, 2016.

[24] C. S. Shui, W. Y. Szeto, "Dynamic green bike repositioning problem—a hybrid rolling horizon artificial bee colony algorithm approach," Transportation Research Part D: Transport and Environment, vol. 60, pp. 119-136, DOI: 10.1016/j.trd.2017.06.023, 2018.

[25] Y. Wang, W. Y. Szeto, "Static green repositioning in bike sharing systems with broken bikes," Transportation Research Part D: Transport and Environment, vol. 65, pp. 438-457, DOI: 10.1016/j.trd.2018.09.016, 2018.

[26] J. Andelmin, E. Bartolini, "A multi-start local search heuristic for the green vehicle routing problem based on a multigraph reformulation," Computers & Operations Research, vol. 109, pp. 43-63, DOI: 10.1016/j.cor.2019.04.018, 2019.

[27] Y. Li, H. Soleimani, M. Zohal, "An improved ant colony optimization algorithm for the multi-depot green vehicle routing problem with multiple objectives," Journal of Cleaner Production, vol. 227, pp. 1161-1172, DOI: 10.1016/j.jclepro.2019.03.185, 2019.

[28] S. Rastani, T. Yüksel, B. Çatay, "Effects of ambient temperature on the route planning of electric freight vehicles," Transportation Research Part D: Transport and Environment, vol. 74, pp. 124-141, DOI: 10.1016/j.trd.2019.07.025, 2019.

[29] N. Rezaei, S. Ebrahimnejad, A. Moosavi, A. Nikfarjam, "A green vehicle routing problem with time windows considering the heterogeneous fleet of vehicles: two metaheuristic algorithms," European Journal of Industrial Engineering, vol. 13 no. 4, pp. 507-535, DOI: 10.1504/ejie.2019.10022249, 2019.

[30] L. Zhu, D. Hu, "Study on the vehicle routing problem considering congestion and emission factors," International Journal of Production Research, vol. 57 no. 19, pp. 6115-6129, DOI: 10.1080/00207543.2018.1533260, 2019.

[31] G. Laporte, R. Musmanno, F. Vocaturo, "An adaptive large neighbourhood search heuristic for the capacitated arc-routing problem with stochastic demands," Transportation Science, vol. 44 no. 1, pp. 125-135, DOI: 10.1287/trsc.1090.0290, 2010.

[32] N. Azi, M. Gendreau, J.-Y. Potvin, "An adaptive large neighborhood search for a vehicle routing problem with multiple routes," Computers & Operations Research, vol. 41, pp. 167-173, DOI: 10.1016/j.cor.2013.08.016, 2014.

[33] M. Keskin, B. Çatay, "Partial recharge strategies for the electric vehicle routing problem with time windows," Transportation Research Part C: Emerging Technologies, vol. 65, pp. 111-127, DOI: 10.1016/j.trc.2016.01.013, 2016.

[34] S. Mancini, "The hybrid vehicle routing problem," Transportation Research Part C: Emerging Technologies, vol. 78,DOI: 10.1016/j.trc.2017.02.004, 2017.

[35] L. Bach, G. Hasle, C. Schulz, "Adaptive large neighborhood search on the graphics processing unit," European Journal of Operational Research, vol. 275 no. 1, pp. 53-66, DOI: 10.1016/j.ejor.2018.11.035, 2019.

[36] W. Gu, D. Cattaruzza, M. Ogier, F. Semet, "Adaptive large neighborhood search for the commodity constrained split delivery VRP," Computers & Operations Research, vol. 112,DOI: 10.1016/j.cor.2019.07.019, 2019.

[37] J. Hof, M. Schneider, "An adaptive large neighborhood search with path relinking for a class of vehicle-routing problems with simultaneous pickup and delivery," Networks, vol. 74 no. 3, pp. 207-250, DOI: 10.1002/net.21879, 2019.

[38] R. Lahyani, A.-L. Gouguenheim, L. C. Coelho, "A hybrid adaptive large neighbourhood search for multi-depot open vehicle routing problems," International Journal of Production Research, vol. 57 no. 22, pp. 6963-6976, DOI: 10.1080/00207543.2019.1572929, 2019.

[39] D. Sacramento, D. Pisinger, S. Ropke, "An adaptive large neighborhood search metaheuristic for the vehicle routing problem with drones," Transportation Research Part C: Emerging Technologies, vol. 102, pp. 289-315, DOI: 10.1016/j.trc.2019.02.018, 2019.

[40] E. Demir, T. Bektaş, G. Laporte, "A review of recent research on green road freight transportation," European Journal of Operational Research, vol. 237 no. 3, pp. 775-793, DOI: 10.1016/j.ejor.2013.12.033, 2014.

[41] D. Pisinger, S. Ropke, "A general heuristic for vehicle routing problems," Computers & Operations Research, vol. 34 no. 8, pp. 2403-2435, DOI: 10.1016/j.cor.2005.09.012, 2007.

[42] V. C. Hemmelmayr, J.-F. Cordeau, T. G. Crainic, "An adaptive large neighborhood search heuristic for two-echelon vehicle routing problems arising in city logistics," Computers & Operations Research, vol. 39 no. 12, pp. 3215-3228, DOI: 10.1016/j.cor.2012.04.007, 2012.

[43] M. Alinaghian, N. Shokouhi, "Multi-depot multi-compartment vehicle routing problem, solved by a hybrid adaptive large neighborhood search," Omega, vol. 76, pp. 85-99, DOI: 10.1016/j.omega.2017.05.002, 2018.

[44] P. Shaw, A New Local Search Algorithm Providing High Quality Solutions to Vehicle Routing Problems, 1997.

[45] M. M. Solomon, "Algorithms for the vehicle routing and scheduling problems with time window constraints," Operations Research, vol. 35 no. 2, pp. 254-265, DOI: 10.1287/opre.35.2.254, 1987.

[46] H. Gehring, J. Homberger, "A parallel hybrid evolutionary metaheuristic for the vehicle routing problem with time windows," vol. 2, pp. 57-64, .

[47] S. Dabia, E. Demir, T. V. Woensel, "An exact approach for a variant of the pollution-routing problem," Transportation Science, vol. 51 no. 2, pp. 607-628, 2016.

[48] 360che, 2020. Dongfeng Ruiling L 143 horsepower single-row light truck chassis (load: 1.99 tons), https://product.360che.com/m87/21973_index.html

[49] 360che, 2020. CIMC Huajun 11.5m low flatbed semi-trailer (load: 10 tons), https://product.360che.com/m280/70031_index.html

[50] Y. Yu, Y. Wu, J. Wang, "Bi-objective green ride-sharing problem: model and exact method," International Journal of Production Economics, vol. 208, pp. 472-482, DOI: 10.1016/j.ijpe.2018.12.007, 2019.