Journal of Cleaner Production 201 (2018) 896e908

Contents lists available at ScienceDirect

Journal of Cleaner Production

journal homepage: www.elsevier.com/locate/jclepro

Heterogeneous fixed fleet vehicle routing problem based on fuel and

carbon emissions
Jin Li a, Danping Wang b, Jianghua Zhang c, *
a
School of Management and E-Business, Key Research Institute-Modern Business Research Center, Zhejiang Gongshang University, Hangzhou, 310018,
China
b
School of Business, Shandong Normal University, Ji'nan, Shandong 250014, China
c
School of Management, Shandong University, Ji'nan, Shandong 250100, China

a r t i c l e i n f o a b s t r a c t

Article history: In this paper, we study an emission-based heterogeneous fixed fleet vehicle routing problem (E-HFFVRP)
Received 1 February 2018 with considerations of fuel and carbon emissions. This problem involves routing a fleet of a fixed number
Received in revised form of vehicles with various capacities to serve a set of customers. It seeks to minimize the objective function,
6 July 2018
which incorporates the fixed expenses and variable costs consisting of fuel consumptions and carbon
Accepted 7 August 2018
Available online 13 August 2018
emissions. It is a new variant of the heterogeneous fixed fleet vehicle routing problem (HFFVRP), in
which a fleet consists of a fixed number of vehicles with different capacities, fixed costs and variable
costs. We formulate this problem with a mixed integer programming model by introducing an approach
Keywords:
Vehicle routing
to calculate fuel and carbon emissions. Moreover, a split-based adaptive tabu search (SATS) algorithm
Freight transportation using an optimal split scheme and an adaptive tabu search algorithm is proposed. Its key features and
Heterogeneous fleet components are designed accordingly. Results of numerical experimentations on two sets of generated
Fuel instances confirm the efficiency and effectiveness of the algorithm.
Carbon emissions © 2018 Elsevier Ltd. All rights reserved.

1. Introduction about 20%e30% through improved routing decision on when and

how to visit its distributed stores. Therefore, it is necessary to
In recent years, the energy crisis and global warming have incorporate fuel and carbon emissions into classical vehicle routing
aroused great concern in the world. To deal with this challenge, problem (VRP) to study a new variant and formulation.
most countries in the world have enacted related carbon polices The classical VRP (Dantzig and Ramser, 1959; Dantzig and
such as carbon cap and trade, carbon tax, and carbon quota. EU has Ramser, 1959) focuses on how to minimize the total travel dis-
implemented carbon cap and trade since 2005 and imposed carbon tance or total costs of all vehicles, which generally has a linear
tax from 2012. Thus, how to effectively save energy and reduce relationship with distance. However, the routing arrangement of
carbon emissions has become the urgent task to address the minimizing the travel distance does not necessarily derive the
world's sustainable development. Transportation has been a major optimal solution of minimizing fuel and carbon emissions, because
source of carbon emissions, given the rapid industrial and eco- the fuel and carbon emissions of vehicles are not only related to
nomic development. In China, the carbon emissions of freight distance, but also to various factors such as vehicle weight, speed
transportation have already accounted for 30% of the total carbon and load Turkensteen (2017).
emissions and the fuel consumption of freight companies is at least In our paper, we study an emission-based heterogeneous fixed
40% of their total costs (Li and Zhang, 2014; Li and Zhang, 2014). fleet vehicle routing problem (E-HFFVRP), which extends the
With strict carbon polices and increasing transportation cost, many traditional VRP by taking fuel and carbon emissions into account. In
companies aim to reduce the carbon emissions and fuel con- particular, there is a heterogeneous fleet of vehicles to conduct
sumptions by optimizing vehicle routing decisions. For example, delivery from one depot to all customers. After finishing the de-
Walmart reduced its transport energies and carbon footprint by livery, the vehicles will return to the depot. Moreover, the number
of vehicles is limited and each type of vehicles available has
different load capacity, fuel consumption, carbon emission and
* Corresponding author. fixed cost. The costs of carbon emissions involve the costs of carbon
E-mail address: [email protected] (J. Zhang).

https://doi.org/10.1016/j.jclepro.2018.08.075
0959-6526/© 2018 Elsevier Ltd. All rights reserved.
 J. Li et al. / Journal of Cleaner Production 201 (2018) 896e908 897

pollution and emission reduction under various environmental amount of fuel and carbon emissions. Traditional approaches usu-
policies, such as carbon tax and carbon cap-and-trade mechanism. ally use a linear function with respect to travel distance or vehicle
Vehicle fixed costs are mainly related with the costs of hiring, load. In contrast, our approach takes into account the effects of not
maintenance, depreciation and insurances, etc. The task of the just travel distance but also lots of other factors, such as vehicle
decision-maker is how to arrange the fleet and routing of these weight, load and speed. Moreover, we propose a mixed integer
different vehicles to achieve the minimization of total transport programming model for the E-HFFVRP, and adapt adaptive memory
network costs. programming to develop a split-based adaptive tabu search (SATS)
Our study is related to green vehicle routing problem (GVRP), algorithm. A few promising characteristics of this algorithm are
which is an alternative to regular well-established vehicle routing introduced to improve its computational performance. In partic-
models. Literature reviews on GVRP are presented in the papers ular, the algorithm incorporates an optimal split scheme and an
such as Lin et al. (2014) and Zhou et al. (2016). GVRP has been used adaptive tabu search algorithm. The optimal split scheme is applied
recently for purposes such as making routing decisions to minimize to generate the adaptive memory, and the adaptive tabu search
fuel consumptions Zhu et al. (2014); Hosseini-Nasab and Lotfalian algorithm is used to improve each temporary solution.
(2017) or reduce carbon emissions Niu et al. (2018); Xiao and The remainder of this paper is organized as follows. Section 2
Abdullah Konak (2017), recovery and distribution of remanufac- introduces an approach of calculating the fuel and carbon emis-
tured products considering sustainability and green criteria sions, and then describes the E-HFFVRP model. In Section 3, the
Soleimani et al. (2018) or routing a capacitated multi-depot dis- algorithm named SATS is presented. Its main features and different
tribution system Jabir et al. (2017). components are given in detail. In Section 4, results of numerical
As a typical GVRP, the pollution routing problem (PRP) is also an studies on two sets of test problems are shown to illustrate the
extension of the traditional VRP and has a broad objective function efficiency of the proposed algorithm and analyze the effect of split
that accounts for carbon emissions. It was first introduced by strategy. Finally, Section 5 summarizes the final conclusions.
Bektas and Laporte (2011), and then some effective approaches are
developed to solve this problem, such as adaptive large neighbor- 2. Mixed integer programming model
hood search Demir et al. (2012) and disjunctive convex program-
ming Fukasawa et al. (2016). Some other variants of PRP are paid 2.1. Problem description
more attention in the literature such as the time-dependent PRP
Franceschetti et al. (2017), PRP with multiple conflicting objective The notations of E-HFFVRP are given as follows. We assume
functions Suzuki (2016), pickup and delivery PRP Tajik et al. (2014) that G ¼ ðV; EÞ is a digraph, V ¼ f0; 1; /; ng is the nodes set, and
and PRP with uncertain demand and travel time (Eshtehadi et al. E ¼ fði; jÞji; j2V; isjg is the arcs set between each pair of nodes. The
(2017). However, in practical business activities, outsourcing has depot denoted by node 0 contains different types of vehicles, and
been widely used. Outsourcing carriers may have limited number of other n nodes constitute a customer set C ¼ f1; 2; /; ng.
heterogeneous vehicles, each type with different capacities, fixed x ¼ f1; 2; /; Kg is the set of vehicles and j ¼ f1; 2; /; Mg is the set
costs, fuel consumptions and carbon emissions. This is so called E- of vehicle types. For each m2j, the vehicle capacity is Qm , and the
HFFVRP, which is to study how to arrange the routing of fixed fixed cost is Fm . The fuel cost per unit is Cf , and the emission cost
number of vehicles from the depot to a series of customers to
per unit fuel is Ce . The number of vehicle type m is fixed and equals
minimize a total cost consisting of fixed costs, fuel and carbon P
emissions. to Nm . Therefore, the total number of vehicles is M m¼1 Nm . For each

Another related problem is heterogeneous fixed fleet vehicle customer i2C, its demand to be delivered from the depot is qi . For
routing problem (HFFVRP), which routes a fixed number of vehicles each arcði; jÞ2E, dij is assumed to be a non-negative distance.
from the depot to visit all customers with minimizing the overall Given above conditions, the objective function of the E-HFFVRP
fixed and variable costs. To address HFFVRP, different heuristics and is to minimize the sum costs of variable expenses (fuel consump-
metaheuristics have been developed, such as adaptive memory tions and carbon emissions) and fixed expenses, by determining a
programming metaheuristic Li et al. (2010), tabu search algorithm vehicle scheduling strategy. Unlike the HFFVRP, which only focuses
Branda ~o (2011), and lifted polynomial size formulations Leggieri on the economic costs, the E-HFFVRP considers both the economic
and Haouari (2017). However, different with our problem, and environmental costs (i.e. fuel cost and carbon emission cost) in
HFFVRP did not consider the effect of fuel and carbon emissions on its objective function.
vehicle routing arrangement.
One closely related literature is from Zhang et al. (2015), whose 2.2. Calculation of fuel and carbon emission costs
focus is on a capacitated vehicle routing problem to minimize the
sum of fuel cost, carbon emission cost and the vehicle usage cost. By Carbon emissions of vehicles depend on a number of factors
contrast, our study is centered on an E-HFFVRP to minimize not Sundarakani et al. (2010), such as transport mode, fuel type, fuel
only the variable costs of fuel and carbon emissions but also fixed consumption, and travel distance. For practical reasons, we intro-
costs. Hence, the study of E-HFFVRP in this paper is a new extension duce an emission calculation approach for road transport from
of VRP in (Zhang et al. (2015), and is quite different from them as Hoen et al. (2010). The total carbon emissions on an arc ði; jÞ are
well. To the best of our knowledge, E-HFFVRP has not yet attracted assumed to be EMij , which is related to FCij and the fuel emission
much attention in the extant literature and an appropriate method factor (FE). Thus, it can be calculated as follows:
that addresses this problem remains unexplored.
Existing literature review shows that pure HFFVRP without fuel EMij ¼ FE$FCij (1)
and carbon emissions is well studied. However, low-carbon trans-
portation and logistics are paid more attention recently by gov- where FE is defined as a gram of CO2 emitted per liter of fuel, i.e., the
ernments and customers. To fill this gap, in this paper, we present efficiency of converting the fuel into carbon emissions, which is
the E-HFFVRP that integrating both economic and environmental affected by vehicle type, fuel type and other factors.
impacts. To reduce fuel consumptions and carbon emissions of The traditional calculation of fuel consumption is a linear
vehicle routing, we use a comprehensive approach to calculate the function of travel distance, which is not in line with industry ap-
plications. Instead, we adopt a comprehensive fuel consumption
898 J. Li et al. / Journal of Cleaner Production 201 (2018) 896e908

model as reflected in Bektas and Laporte (2011) that considers more

precise factors, including vehicle load, speed, and distance as X
n X
K
xm
0jk  Nm ; cm (9)
shown in the following.
j¼1 k¼1
h   i
FCij ¼ aij w þ zij þ bv2ij dij (2)
X
M X
K X
M X
K
qj xm
ijk  zij  ðQm  qi Þxm
ijk ; ci; j (10)
where vij is a vehicle's average speed on an arc ði;jÞ; L ¼ w þ zij is the m¼1 k¼1 m¼1 k¼1
total vehicle load, w denotes the empty vehicle weight, and zij is a
 
vehicle's load on this arc; aij represents an arc-dependent constant  
X  
related to acceleration, road angle, etc.; and b is a vehicle-specific xm
ijk  S  1; S3V\f0g; Ssf; cm; k (11)
constant that depends on the vehicle frontal surface area, air den- i;j2SS  
sity, and rolling resistance.
h
In terms of Eq. (2), the total cost TCij that consists of fuel con-
xm m
ijk 2f0; 1g; yik 2f0; 1g; zij 2 0; þ∞Þ; ci; j; m; k (12)
sumptions and carbon emissions is shown as follows:
  The objective function (4) has two parts. The first part measures
TCij ¼ cf þ ce FCij (3) the total variable costs of fuel consumptions and carbon emissions.
The second part measures the total fixed costs for vehicles.
Table 1 summarizes some symbols and values (ranges) used in Constraint (5) ensures that every customer can only be visited by a
the E-HFFVRP model and experiments. To be consistent with vehicle once. Constraint (6) and (7) define that every customer
practical applications, these data are derived from one large logis- enters and leaves one vehicle exactly once. Constraint (8) is the
tics company in China and some related literature Bektas and capacity limitation of vehicles. Constraint (9) imposes a limitation
Laporte (2011); Zhang et al. (2015). on the number of each vehicle type. Constraint (10) guarantees that
the vehicle's load on an arc ði; jÞ cannot exceed its remaining ca-
pacity. Constraint (11) is used to eliminate sub-loops. Constraint
2.3. Model formulation (12) defines xm , ym , and Zij as decision variables; Given a vehicle k of
ijk ik
type m and arc ði; jÞ, xm
ijk
¼ 1 if and only if the vehicle of type m
Given above analysis and integer programming approaches Wu
passes through arc ði;jÞ, otherwise, xm
ijk
¼ 0, Given a vehicle k of type
et al. (2010); Wu and Zhou (2015), we formulate E-HFFVRP as a
mixed integer programming model below: m and node i, ym ik
¼ 1 if and only if this vehicle visits node i, other-
wise it is 0. Given an arc ði; jÞ, the load zij  0 if it is visited by a
X
M X
K X
n n 
X h   i vehicle, otherwise zij ¼ 0.
min Cf þ Ce aij w þ zij þ bv2ij Þdij xm
ijk
m¼1 k¼1 i¼0 j¼0;jsi

X
M X
K X
n 3. Metaheuristic for the E-HFFVRP
þ Fm xm
0jk
m¼1 k¼1 j¼1 Since the E-HFFVRP is a new extension of the classical VRP, a
(4) metaheuristic algorithm should be used to solve this NP-hard
problem. Based on the general framework of adaptive memory
programming procedure Olivera and Viera (2007), a split-based
M X
X K
ym adaptive tabu search (SATS) algorithm that incorporates an
s:t: ik ¼ 1; cis0 (5)
m¼1 k¼1 optimal split scheme and embedded adaptive tabu search algo-
rithm is proposed. To facilitate the presentation, we define the
following notations: a is the adaptive memory that comprises some
X
n
xm m
ijk ¼ yjk ; cj; m; k (6) solutions; b is the set of temporary solutions constructed based on
i¼1 a at each iteration; nb is the number of temporary solutions found
at each iteration; xi is the ith solution in b; x0i is the ith improved
X
n solution for xi ; xcb is the best solution obtained in the current
xm m
ijk ¼ yik ; ci; m; k (7) iteration; xgb is the global-best solution obtained from the start of
j¼1
the algorithm; and nb is the number of best solutions in the current
iteration. In order to help the readers understand the whole picture
X
n
of the proposed algorithm, the basic structure of the SATS is shown
qi ym
ik  Qm ; ck; m (8) by a flow chart in Fig. 1.
i¼1
The steps of SATS are described in the following. First, a is set to
store the E-HFFVRP solutions, which are denoted as sequences
without delimiters and then organized by a split algorithm. To
enhance the diversification strategies, the solutions in a are
Table 1
selected in probability to generate temporary solutions in each
Symbols representation referred in the E-HFFVRP model.
iteration step. Next, an adaptive tabu search (ATS) is used to
Symbols Descriptions Values (ranges) improve each temporary solution. To better utilize the search
cf Fuel cost per liter (CNY) 7.30 memory and history, several best improved solutions are selected
ce Emissions cost per liter (CNY) 0.64 to update a. The above steps are repeated until the maximum
FE Fuel emissions factor (kilogram/liter) 2.32 number of iterations IM is satisfied. Finally, after the main loop, a
vij Average arc speed (kilometers/hour) 40e100
post-optimization procedure (POP) is added to further improve the
aij Arc specific constant 0.09e0.15
best solution. Specifically, the algorithmic framework of SATS
 J. Li et al. / Journal of Cleaner Production 201 (2018) 896e908 899

W contains weights of arcs. The weight wij 2W reflects the trans-

portation cost (condition (15)).

X
j
cði; jÞ2A : qSk  Qm (13)
k¼iþ1

Nam ¼ Nm  Num  1 (14)

cði; jÞ2A : wij

!
  X
j1
¼ cf þ ce FC0;Siþ1 þ FCSk ;Skþ1 þ FCSj ;0 þ Fm
k¼iþ1
(15)
Therefore, an optimal E-HFFVRP solution expressed by S is
actually the minimum cost route r from 0 to n. Let qmin ¼

minfqSk k ¼ 1; 2; /; ng, and the number of customers traversed by
a sub-route is not more thana ¼ bQ m =qmin c. Then, route r can be
derived within OðmnaÞ using OS.
Algorithm 2 in Appendix presents a detail description of OS.
Algorithm 2 computes several labels for each node j2N\f0g: V½j is
the cost of the minimum cost route from node 0 to node j; P½j
denotes the predecessor of node j on this route; T½j is the type of
vehicle used on this route; and Nm ½j, cm2j, is the number of
available vehicles of type m on this route. In step 2, the “do while”
loop enumerates all feasible sub-sequences Si ; /; Sj with different
types of vehicles to find the min-cost route and directly update
these labels. If the label of the last node is updated, the final total
cost is expressed as V½n. However, some un-routed customers may
remain at the end of the iterations.
Algorithm 2 does not generate the E-HFFVRP solution explicitly.
Thus, Procedure 1 in Appendix is used to extract the E-HFFVRP
solution from vectors of labels Pand T. For each route shown as a list
Fig. 1. Flow chart of the SATS. of customers, Procedure 1 starts with putting each customer node
in one route and thus initializes n routes and labels of vehicle type.
If a route is empty at the end of the algorithm, then the route does
incorporates an optimal split algorithm and an adaptive tabu search not belong to the solution. The number of routes used is computed
as shown in Algorithm 1 in Appendix.
by the counter t. The insertion cost cki of route rðkÞ including
customer i is calculated by the following formula:
3.1. Solution denotation and split   
cki ¼ cf þ ce FCi1 ;i þ FCi;i2  FCi1 ;i2 (16)
Most solution denotations for the VRP regard the depot as the
route's delimiters. For example, the tour (0-1-2-3-0-4-5-0) repre- where i1 and i2 are the predecessor node and successor node for a
sents the sub-route (0-1-2-3-0) and (0-4-5-0). Instead, we simply given insertion location, respectively. Each un-routed customer in R
denote a solution tour as a sequence S composed of n customer is put into one of the available routes with minimum insertion cost.
without delimiters, such as the above tour (0-1-2-3-0-4-5-0) be- As the vehicle capacity is not considered, infeasible solutions may
comes (1-2-3-4-5). This denotation approach can transform VRP be generated. As a diversification strategy, SATS allows cross-border
into a simpler traveling salesman problem. On the other hand, it is search between feasible and infeasible spaces to improve the
essential to split S into sub-routes to evaluate the solution. Some quality of the final solutions.
splitting procedures have been proposed to compute the fitness,
such as Christian (2004), but they did not consider the effect of fleet
composition, limited number of vehicles and the computation of 3.2. Solution evaluation
fuel and carbon emissions. Accordingly, we develop an improved
optimal splitting procedure called Optimal Split (OS) to derive the SATS allows exploring the infeasible solutions, which means
best E-HFFVRP solution corresponding to the denotation sequence. that there might exist a solution of violating the vehicle capacity.
The main idea of OS is that for any kind of optimal tour split, the We define a penalty function FðxÞ to evaluate a solutionx.
remaining tour split is also optimal when subsequent several sub-
routes are removed. OS works on an auxiliary graph H ¼ ðN; A; WÞ, X
K
FðxÞ ¼ f ðxÞ þ pli (17)
where N includes n þ 1 nodes indexed from 0 to n; A contains route
i¼1
arcs, for any arc ði;jÞ2A, i < j, if one vehicle of type m traversing Siþ1
to Sj is in existence satisfying vehicle load (condition (13)) and where f ðxÞ represents the original objective value of solution x; K
number limitation (condition (14), where Nam is the number of type denotes the number of vehicles used; li is the load excess in vehicle
m vehicles available, Num is the number of type m utilized vehicles); (or route) i; and p is a penalty term.
900 J. Li et al. / Journal of Cleaner Production 201 (2018) 896e908

3.3. Adaptive memory initialization and temporary solutions one of the Nd ðSi Þ. The double insertion removes any pair of cus-
generation tomers ðSi ; Sj Þ from one tour and inserts them into any different
position of another tour that has at least one of the Nd ðSi Þ or Nd ðSj Þ.
In the initial stage of SATS, a is composed of various solutions. The swap allows a customer node Si from one tour to be exchanged
Let am be the size of a. We generate random permutation of with a customer node Sj from another tour while satisfying that at
customer nodes as initial solutions. Then, these solutions are ar- least one customer in route j falls inNd ðSi Þ. For the 2-opt, non-
ranged in descending order and put in a. First, Each solution x is adjacent arcs ðSi ; Siþ1 Þ and ðSj ;Sjþ1 Þare replaced by ðSi ; Sj Þ and ðSiþ1 ;
assigned with a triple ðfðxÞ; f ðxÞ; FðxÞÞ, where fðxÞ is the total load Sjþ1 Þ in the same tour, and then the direction of the sub-tour ðSiþ1 ;
excess value of the capacity in x. Second, there are three possible
/; Sj Þ which has at least one of the Nd ðSi Þ or Nd ðSjþ1 Þ is reversed.
cases in which these solutions are sorted: (1) When x1 and x2 are
both feasible, x1 has a higher level than x2 if and only if f ðx1 Þ < f ðx2 Þ;
(2) when x1 is feasible and x2 is infeasible, x1 has a strictly higher 3.4.2. Tabu status and tabu tenure
level than x2 ; (3) when x1 and x2 are both infeasible, x1 has a higher To avoid local cycling, we don't allow recent solutions found to
level than x2 if and only if Fðx1 Þ < Fðx2 Þ. Note that if f ðx1 Þ ¼ f ðx2 Þ or be used in the next q iterations and put some relevant attributes of
Fðx1 Þ ¼ Fðx2 Þ, x1 and x2 is chosen randomly to be assigned with a these solutions into a tabu list. For each action that generates a new
higher level. solution, the inverse action will be taboo in subsequent q iterations
A multi-start strategy is used to enhance diversification and thus and is then inserted into the tabu list. The status of an operator is
increase the probability of generating better solutions, that is, nb required that a customer removed from a position of a tour cannot
solutions in a are probabilistically selected to be temporary solu- go back to it during the tabu tenure. The tubu tenure is changed
tions at each iteration. As mentioned above, a solution with a higher systematically as the search evolution. Due to the fact that the tabu
level is more likely to be selected in probability. Hence, the prob- restriction may prohibit promising moves, a common aspiration
ability of selecting the ith solution is defined as pi ¼ ðns  i þ 1Þ=n2s , criterion is used if the move action produces a better solution.
where ns represents the total number of available solutions in a. In
terms of the solution level in a, the first solution with the highest
level is certain to be selected, and the final one with the lowest level 3.4.3. Restart and shaking
has least selection probability. As tabu search algorithm is easy to fall into local optimal, to
address this problem, restart and shaking is often used to obtain
3.4. ATS for solution improvement better solutions. By changing parameters, restart can make a good
tradeoff between intensification and diversification. The restart of
An improvement-based procedure is applied to improve all our algorithm is done from the best-so-far solution xgb if no
solutions in b and push to their local optima. We use an ATS as a improvement is found within iL iterations. Furthermore, to escape
solution improvement procedure in SATS. Some related notations local optimum at the end of the search, it is necessary to shake the
in the ATS are listed as follows: i0 represents the iteration counter of best known solution and reach better solutions. Shaking technique
no improvement in current best feasible or infeasible solution; iL is is often used when the search falls into local optimum within a
the limit of i0 in the current cycle; q is the tabu tenure; d is the given neighborhood structure. In our implementation, a random
number of nearest neighbors of a customer; IA is the total limitation move is performed on the current neighborhood by randomly
number of iterations; FS FD ,FW and FO are the frequencies of the executing one of the three movements: GENI, US Gendreau et al.
neighborhood operators of single, double, swap and 2-opt, (1994) and 2-opt.
respectively. A typical feature of ATS is that it can adaptively adjust
several key parameters with search evolution, such as the tabu
3.5. Adaptive memory update
tenure and the frequencies of the neighborhood operators. In
particular, once the search reaches to the middle of the cycle (iL =2)
After the solution improvement procedure ATS is implemented,
or the end of the cycle (iL ), accordingly, the ATS changes some
several poor solutions are removed from a and replaced by the
parameters (qdFS FD FW FO ) to enhance the diversification. POP and improved solutions. First, we select nb (nb < nb ) best solutions from
restart strategy are also used to achieve the intensification. The
nb improved solutions according to the following conditions: (1)
outline of the ATS is described in Algorithm 3 in Appendix.
The feasible solution should be selected in advance when feasible
and infeasible solutions are all available; (2) When only feasible
3.4.1. Neighborhood structure solutions are derived, we select the one with the least objective
In a tabu search algorithm, the degree and quality of solution value; and (3) When only infeasible solutions are generated, we
search are determined by neighborhood structure. We use three select the one with the least penalty value. Then, these solutions are
intra-route neighborhood operators (single insertion, double assigned and sorted with the triple given in Section 3.3. Finally, the
insertion and swap) and one inter-route neighborhood operator (2- nb worst solutions in a are deleted and replaced by the same
opt). We limit the search of these operators by d  neighborhood number of best solutions with assigned order level.
restriction to reduce the number of potential moves at each itera-
tion, which thus not only saves the computing time but also drives
the search to find better solutions. For any node s2S, we define its 3.6. POP
d  neighborhood Nd ðsÞ as the set of the d nodes on the tour that is
closet to s with respect to the total costs. This concept was proposed At the end of the ATS or SATS, in order to promote intensifica-
and used by Gendreau et al. (1994). However, the main difference tion, we use a POP to further improve the best found solution.
with them is that the d in our algorithm varies with the search Similar to those operators in the shaking stage, GENI, US, and 2-opt
evolution, that is, the initial d is very small and increases gradually are also executed. However, they are not used by random selection,
during the search evolution. instead they are implemented in sequence using the best-accepted
The single insertion removes a customer Si from one tour and criterion, i.e., they search the entire neighborhood of the solution
inserts it into various positions in another tour that contains at least space to find the best improvement.
 J. Li et al. / Journal of Cleaner Production 201 (2018) 896e908 901

4. Numerical experiments dependent on number of customers. Some parameters will be

varied adaptively with the search evolution of the algorithm to
4.1. Benchmark problems induce diversification. The initial q is defined to be n=2 according to
a set of experiments. In the middle of the cycle, the q should be
SATS is programmed in C language and executed on a PC with reduced to allow intensification. However, at the end of the cycle,
Intel Duo CPU 2.26G, 2G of RAM, under Windows 7. This paper is current solutions in the search space cannot be improved, and q
the first attempt to use instance to test the E-HFFVRP. The sets of should be added to drive the search toward more distant spaces.
benchmark problems in existing papers are unable to be used. Therefore, for a cycle, we set q ¼ maxðq=2; 6Þ in the middle, and at
Hence, we plan to test the performance of SATS by designing two the end, we change it to be q ¼ minð2q þ 3; 0:6nÞ.
sets of benchmark problems, whose configurations for the E- Similar with the setting of the tabu tenure, we change d with the
HFFVRP are shown in Tables 1 and 2. Set 1 is generated in Table 2 search evolution as well. A set of experiments showed that the
based on Taillard (1999), which comprises eight problems labeled initial d must be small. Here, it is set as d ¼ 1. In the middle of the
from 13 to 20. Six vehicle types labeled from A to F are given. In Set cycle, it is set as d ¼ minðd þ 1;n  1Þ. And at the end of the cycle, no
1, we generate the number of customers according to the integers better solution could be found. To enhance diversification in the
between 50 and 100. search space, d is desirable to increase, and thus we set it as d ¼
We also generate a new set of benchmark problems, i.e., Set 2 in minðd þ 2;n  1Þ. According to our previous experiments, the value
Table 3, which are composed of eight problems denoted by P1eP8 of p(penalty term) was set as follows: At the initial stage, it equals to
with the number of customers ranging from 150 to 300. All cus- Pn
i¼1 qi =100n. If ten successive iterations all generate infeasible
tomers for each problem are located in a Euclidean plane. The co-
solutions, then p ¼ 2p. If ten successive iterations all generate
ordinates ðxi ; yi Þ of customer i (i ¼ 1; 2; /; n) are randomly
feasible solutions, then p ¼ p=2.
produced in square interval ½0; 1002 in a uniform distribution, Finally, we study how to set appropriate frequencies of four
where (30, 60) represents the depot. We set the demand of a neighborhood operators. According to the degree of conversion to a
customer as any integer randomly selected in interval (Zhou et al., solution, the four neighborhood operators are sorted from low to
2016; Christian, 2004). The fleet contains six types of vehicles: A, B, high as follows: single insertion, swap, double insertion, and 2-opt.
C, D, E, and F. In the headings of Tables 2 and 3, VT and VP represent With multiple iterations, single insertion can result in the same
the vehicle type and property respectively. For vehicles of type i, bi conversion caused by the other two operators (swap and double).
is a vehicle-specific constant, and r0 is the ratio of total demand and Therefore, single insertion should be used most frequently. Since
total capacity, i.e., (total demand/total capacity)  100. According to double insertion and swap can produce similar conversions, they
the characteristics of these two sets of problems, they have should not appear in the same iteration. Through a set of trials, we
important implications in the resolution. have the following frequencies of the operators: At the beginning of
the cycle, FS ¼ 1, FD ¼ 10, and FW ¼ 5; in the middle of the cycle,
FD ¼ 5, FW ¼ 10, and FO ¼ 0; at the end of the cycle, FD ¼ 0, FW ¼ 2,
4.2. Parameter tuning and FO ¼ 6.
According to our experiments, good values used for other pa-
In SATS, how to set some key parameters has a great impact on rameters are set and summarized as follows:
the algorithmic performance. Some characteristics in ATS are in line
~o (2011). Thus, pffiffiffi
with that in an algorithm for solving HFFVRP Branda (1) IA ¼ 100 n;
previous knowledge is used to tune the value of several ATS pa- (2) Initially iL ¼ 1.5n, and at the end of the cycle, iL ¼ iL þ0.5n;
pffiffiffi pffiffiffi pffiffiffi
rameters, such as qd, and the frequencies of the neighborhood op- (3) am ¼ 10n, nb ¼ n, nb ¼ n=2, IM ¼ 3000 n.
erators. For other parameters setting, an F-Race approach is used to
automatically achieve the best configuration of SATS in a set of
experiments. The F-Race can set some evaluation metrics and then
automatically derive the best configuration of an algorithm through 4.3. Numerical results
preliminary experiments (Birattari, 2005) (Birattari, 2005). The
tuning results indicated that the parameter combination settings 4.3.1. Effect of split strategy
work well in our proposed algorithm. In this paper, we use a novel optimal split strategy in the pro-
In order to ensure that the SATS performance is consistent with posed algorithm, SATS. We conduct an experiment based on the
the size of the problem, most of the parameters are set to be given parameters of configuration tuning to assess the quality of

Table 2
Data for the testing problems of Set 1.

VT VP Problem VT VP Problem

13 14 15 16 17 18 19 20 13 14 15 16 17 18 19 20

n 50 50 50 50 75 75 100 100 n 50 50 50 50 75 75 100 100

r0 95.39 88.45 94.76 94.76 95.38 95.38 76.74 95.92 r0 95.39 88.45 94.76 94.76 95.38 95.38 76.74 95.92

A QA 20 120 50 40 50 20 100 60 D QD 70 350 150

FA 20 100 100 100 25 10 500 100 FD 120 320 180
NA 4 4 4 2 4 4 4 6 ND 4 1 2
bA 0.54 2.67 1.59 0.73 1.59 0.54 2.60 1.84 bD 1.91 4.05 3.05
B QB 30 160 100 80 120 50 200 140 E QE 120 250
FB 35 1500 250 200 80 35 1200 300 FE 225 400
NB 2 2 3 4 4 4 3 4 NE 2 1
bB 0.70 3.13 2.60 2.11 2.67 1.59 3.37 3.01 bE 2.67 3.70
C QC 40 300 160 140 200 100 300 200 F QF 200 400
FC 50 3500 450 400 150 100 2100 500 FF 400 800
NC 4 1 2 3 2 2 3 3 NF 1 1
bC 0.73 3.75 3.13 3.01 3.37 2.60 3.75 3.37 bF 3.37 4.14
902 J. Li et al. / Journal of Cleaner Production 201 (2018) 896e908

Table 3
Data for the testing problems of Set 2.

VT VP Problem VT VP Problem

P1 P2 P3 P4 P5 P6 P7 P8 P1 P2 P3 P4 P5 P6 P7 P8

n 150 150 200 200 250 250 300 300 n 150 150 200 200 250 250 300 300

r0 88.53 95.43 96.40 93.90 88.15 95.77 95.52 94.29 r0 88.53 95.43 96.40 93.90 88.15 95.77 95.52 94.29

A QA 80 100 150 100 60 120 200 150 D QD 200 400 350 300 250 300 350 300
FA 300 350 400 320 200 340 600 410 FD 550 2200 1500 1200 750 1100 1800 1200
NA 4 3 4 2 6 5 3 5 ND 3 1 2 2 3 2 2 2
bA 2.11 2.60 3.05 2.60 1.84 2.67 3.37 3.05 bD 3.37 4.14 4.05 3.75 3.70 3.75 4.05 3.75
B QB 120 160 200 160 120 160 250 200 E QE 250 400 350 300 350 400 400
FB 350 500 550 450 350 500 800 600 FE 600 2100 1800 1200 1600 2100 2100
NB 4 3 3 4 4 5 3 3 NE 1 2 1 2 1 1 2
bB 2.67 3.13 3.37 3.13 2.67 3.13 3.70 3.37 bE 3.70 4.14 4.05 3.75 4.05 4.14 4.14
C QC 150 300 300 250 200 200 300 250 F QF 350 350 400
FC 500 1500 1200 800 510 560 1500 850 FF 2200 1600 2100
NC 3 3 3 3 2 4 3 3 NF 1 2 1
bC 3.05 3.75 3.75 3.70 3.37 3.37 3.75 3.70 bF 4.05 4.05 4.14

the optimal split strategy. The experiment is conducted by per- 10.11% for Set 1 instances and from 2.93% to 12.92% for Set 2 in-
forming four versions of SATS with different split strategies as stances. In particular, the average deviation from the random split
follows: (1) SSA-SATS, a SATS with a simple split strategy in which strategy is above 10%. Moreover, when the split strategy is used, the
available vehicles of different types are assigned in ascending order computational time required to obtain the best solutions is not
of capacity; (2) SSD-SATS, another simple split strategy in which significantly increased. These figures suggest that the optimal split
available vehicles of various types are assigned in descending order strategy can be used to obtain good-quality solutions. Therefore,
of capacity; (3) RS-SATS, where a random split strategy is imple- the optimal split strategy should be adopted in the SATS to derive a
mented, i.e., available vehicles of different types are assigned at better balance between computational time and solution quality.
random; and (4) OS-SATS, where the proposed optimal split
strategy is considered in SATS. Table 4 summarizes the comparison
results of testing problems in Set 1 and Set 2. For each instance, we
run four algorithms each 20 trials, and present the average value of 4.3.2. Comparison of SATS with other algorithms
the best solutions in the column labeled “Average” along with the In this subsection, we carry out simulation experiment to show
related CPU average time used in computing the algorithm. For the the performance of SATS by comparing it with some existing al-
last three columns, the column entitled DevA (%), DevD (%) and gorithms. For SATS, we also use the best parameter settings
DevR (%) shows the percentage deviations of OS-SATS from SSA- selected in the above sections. The selected algorithms include a
SATS, SSD-SATS, and RS-SATS, respectively. Furthermore, suppose pure tabu search and two methods used in HFFVRP as follows:
SðOSÞ is the optimal solution found with OS-SATS, and SðADRÞ
represents the solution value produced by one of the other three (1) TS: a tabu search algorithm, which is similar to ATS, except
split-based algorithms. Then, the percentage deviation is calculated that it use saving algorithm rather than the split-based
as 100ðSðADRÞ  SðOSÞÞ=SðADRÞ. procedure to construct initialization solution. Its termina-
As shown in Table 4, from both solution quality and computa- tion condition is also the maximum number of iterations,
pffiffiffi
tional efforts, the optimal split strategy outperforms the other split 3000 n.
strategies. Moreover, the optimal split strategy is able to discover (2) BATA: a back-tracking adaptive threshold accepting meta-
the best solution in nearly all instances. The average deviations heuristic Tarantilis et al. (2004), which was implemented on
from the other three split strategies are significant from 4.98% to MS Visual Cþþ, Pentium II, 400 MHz, 128 MB RAM.

Table 4
The effect of split strategy in SATS.

Instances SSA-SATS SSD-SATS RS-SATS OS-SATS DevA (%) DevD (%) DevR (%)

Average CPUs Average CPUs Average CPUs Average CPUs

13 2081.12 100 2085.96 109 2211.47 113 2065.96 95 0.73 0.96 6.58
14 11121.32 107 10192.37 105 11444.35 118 10034.51 98 9.77 1.55 12.32
15 2604.07 113 2591.58 108 2690.05 110 2520.39 90 3.21 2.75 6.31
16 2675.01 98 2658.66 104 2652.51 107 2631.22 89 1.64 1.03 0.80
17 2050.64 149 1991.98 153 2192.88 161 1825.06 143 11.00 8.38 16.77
18 2921.95 153 2859.11 160 3178.75 154 2723.14 151 6.80 4.76 14.33
19 11248.84 236 11455.49 227 11778.69 231 10215.39 225 9.19 10.83 13.27
20 4474.95 215 4493.63 208 4538.69 209 4060.80 206 9.25 9.63 10.53
Average 6.45 4.98 10.11
P1 12822.67 314 15005.69 290 15296.87 285 12950.24 283 0.99 13.70 15.34
P2 14386.42 326 16320.36 313 16318.78 328 13668.71 307 4.99 16.25 16.24
P3 21527.50 447 23377.88 435 23412.30 427 21273.17 416 1.18 9.00 9.14
P4 21384.31 412 22302.72 405 24760.75 410 20837.67 397 2.56 6.57 15.84
P5 23439.18 491 25400.64 462 24239.18 481 21658.77 446 7.60 14.73 10.65
P6 24509.41 524 26456.88 482 27030.09 497 23053.85 468 5.94 12.86 14.71
P7 26428.75 655 28060.79 624 28991.76 628 25985.84 611 1.68 7.39 10.37
P8 31338.07 673 34714.78 637 35043.96 649 31175.32 633 0.52 10.20 11.04
Average 2.93 11.34 12.92
 J. Li et al. / Journal of Cleaner Production 201 (2018) 896e908 903

(3) MAMP: a multistart adaptive memory programming met-

CPUs/Tflop
aheuristic (Li et al., 2012; Li et al., 2010), which was
implemented on MS Visual Cþþ, Intel CPU 2.2G.

287/14350
200/10000
110/5500

148/7400
119/5950
As the algorithms of BATA and MAMP have been tested on

34/1700
46/2300
99/4950
benchmark problems Set 1, we use Set 1 to conduct an experi-
ment for all algorithms. Note that different algorithms are
executed on different computer environments. To offer a fair
MAMP

comparison, we measure the computational speed in ten millions

Best

of floating-point operations per second (Tflop/s) based on Don-

10194.17
garra tables Dongarra (2006) as follows: Pentium IId8 Tflop/s,
1945.06
9710.34
2362.29
2473.69
1731.02
2467.35

4015.58
Intel 2.2Gd50 Tflop/s, Intel 2.26Gd55 Tflop/s. All benchmark
problems found through 20 runs of SATS and TS are compared.
Table 5 presents detailed results of the experiment. For each pair
of instance and algorithm, this table includes the best solution
(Best) and the computational time (CPUs/Tflop) obtaining the
CPUs/Tflop

1156/9248
843/6744
387/3096
368/2944
341/2728
363/2904
971/7768
428/3424

best solution over 20 trials. Furthermore, for SATS, we also show

the total number of vehicles (NV) available and the number of
vehicles generated in the best solutions (Vehicles). Table 6 shows
the relative deviations between the best solution of each algo-
rithm (Best) and the best solution found by all these four algo-
rithms (BSA). The percentage deviation is computed as 100(Best
10210.61
1990.05
9755.32
2362.29
2510.63
1739.19
2472.82

4037.54
BATA

eBSA)/Best.
Best

As shown in Tables 5 and 6, we can conclude the following

observations. Generally, SATS has a good performance on the
benchmark problems Set 1. Specifically, SATS can improve the
best solution of other methods to most testing problems except
CPUs/Tflop

283/15565
271/14905
120/6600
129/7095
115/6325
122/6710
167/9185
175/9625

problems 16 and 19. MAMP obtains the best solution on problem

16. For problem 19, both SATS and MAMP have the same ability to
find the best solution. In addition, SATS has a stable performance
because the average percentage deviation is the smallest, that is,
0.0109%. These values for TS, BATA and MAMP are 1.5976%,
1.3217% and 0.6287%, separately. Moreover, the vehicles used in
10213.57
1986.33
9757.88
2376.56
2541.86
1738.22
2488.68

4035.25

SATS are fewest in available vehicles for many instances, such as

Best

14, 17 and 18. Finally, on average, the computing time for SATS is
TS

faster in CPU seconds, but a bit lower in the measurement of

floating-point operations. This finding means that SATS can
significantly and efficiently reduce the transportation cost and
CPUs/Tflop

241/13255
238/13090

carbon emissions in industry operations.

103/5665
112/6160

101/5555
158/8690
163/8965
97/5335

Our study provides some managerial and application insights.

First, due to the environmental policies emerged in the recent
decades, such as carbon tax and carbon cap and trade mechanism,
which are costly in practical transportation and distribution, lo-
gistic service providers should bring environmental factors (such
10194.17
1920.31
9670.36
2358.49
2475.85
1683.63
2458.04

4012.80

as fuel and carbon emissions) into the routing arrangement to

SATS
Comparison results of different algorithms on benchmark instances.

carry out an environment-friendly transportation. The proposed

Best

method can optimize their operational decisions and help them

reduce their environmental costs effectively. Second, we still
consider the special cases of our model, which seek to minimize
the economic costs and carbon emissions respectively. We find
Vehicles

that there exists tradeoff relationship between economic costs

17

10
12

13
6
9
9

and carbon emissions. For the logistic service providers, it is

valuable to use the proposed method to reduce the carbon

Table 6
NV

17

11
14
10
13
7
9
9

The percentage deviations of each algorithm from the best solution of all algorithms.

Instances BSA SATS TS BATA MAMP

13 1920.31 0 3.3237 3.5044 1.2725

100
100
50
50
50
50
75
75

14 9670.36 0 0.8969 0.8709 0.4117

n

15 2358.49 0 0.7603 0.1609 0.1609

16 2473.69 0.0872 2.6819 1.4713 0
17 1683.63 0 3.1406 3.1946 2.7377
18 2458.04 0 1.2312 0.5977 0.3773
Instances

19 10194.17 0 0.1899 0.1610 0

Table 5

20 4012.80 0 0.5563 0.6127 0.0692

13
14
15
16
17
18
19
20

Average 0.0109 1.5976 1.3217 0.6287

904 J. Li et al. / Journal of Cleaner Production 201 (2018) 896e908

emissions at the minimal costs. Finally, from this paper, the policy- given instances.
maker can know how to formulate and adopt effective govern- Although the SATS algorithm yields good results in the given test
mental polices to affect the environmental performance of the lo- instances, some limitations should be noticed. Firstly, since the E-
gistic service providers. For example, sufficiently higher carbon HFFVRP is an NP-hard problem, if the data structure of the prob-
price may reduce the carbon emissions at a lower degree, but lems is large dimension and uniform, some additional difficulties
depress the logistic firms by increasing more economical costs. will be imposed on the algorithm. Secondly, in terms of the good
performance of the SATS, we believe that this algorithm can be
adapted to efficiently deal with other VRPs found in real-life, but
5. Conclusion several features with the algorithm may be returned, such as tabu
tenure, the frequencies of the neighborhood operators and shaking
In light of growing concerns about global warming, all modern mechanism. Moreover, other efficient algorithms will be incorpo-
industries are starting to increase their investment in reducing fuel rated in the solution improvement to solve more difficult problems
consumptions and carbon emissions. To achieve this objective, the with other strong constraints. In the future, we shall test the per-
transportation industries should optimize vehicle routings by formance of our algorithm through other large-scale real data.
focusing on minimizing variable costs (fuel and emission) and fixed
costs. Accordingly, we studied the emission-based heterogeneous
fixed fleet vehicle routing problem (E-HFFVRP) with reducing fuel Acknowledgments
and carbon emissions. A mixed integer programming model was
built along with introducing an approach of calculating fuel and This work was partly supported by a grant from the Natural
carbon emissions. A split-based adaptive tabu search algorithm is Science Foundation of China (Grant No. 71571111); Zhejiang Pro-
proposed, which combines an optimal split procedure with an vincial Philosophy and Social Sciences Foundation of China (Grant
adaptive tabu search. A vehicle routing is denoted as a sequence No. 18NDJC180YB); Zhejiang Provincial Natural Science Foundation
without delimiters. The split algorithm is then used to extract an of China (Grant No. LY17G020005), and the National Natural Sci-
optimal solution from the sequence. In addition, an adaptive tabu ence Foundation of China (NSFC) (Grant Nos. 71302035, 71672179).
search (ATS) is embedded in SATS to improve the solution. We The authors also would like to thank the Qilu Young Scholars and
generated different sets of instances to evaluate the effectiveness of Tang Scholars of Shandong University for financial and technical
SATS. Then, some parameters are tuned systematically based on support.
previous knowledge and F-Race procedure. We tested this algo-
rithm on a series of problems. The results showed that optimal split Appendix
strategy could be incorporated to enhance the performance of SATS,
and the SATS can effectively search for high quality solutions on the
J. Li et al. / Journal of Cleaner Production 201 (2018) 896e908 905
906 J. Li et al. / Journal of Cleaner Production 201 (2018) 896e908
 J. Li et al. / Journal of Cleaner Production 201 (2018) 896e908 907

References problem with demand and travel time uncertainty. Transport. Res. Transport
Environ. 51, 351e363.
Franceschetti, A., Demir, E., Honhon, D., Van Woensel, T., et al., 2017. A metaheuristic
Bektas, T., Laporte, G., 2011. The pollution-routing problem. Transport. Res. Part B 45 for the time-dependent pollution-routing problem. Eur. J. Oper. Res. 259 (3),
(8), 1232e1250. 972e991.
Birattari, 2005. The Problem of Tuning Metaheuristics as Seen from a Machine Fukasawa, R., He, Q., Song, Y.J., 2016. A disjunctive convex programming approach to
Learning Perspective. PhD thesis. Universite Libre De Bruxelles. the pollution-routing problem. Transp. Res. Part B Methodol. 94, 61e79.
Branda~o, J.A., 2011. Tabu search algorithm for the heterogeneous fixed fleet vehicle
Gendreau, M., Hertz, A., Laporte, G., 1994. A Tabu search heuristic for the vehicle
routing problem. Comput. Oper. Res. 38 (1), 140e151. routing. Manag. Sci. 40, 1276e1290.
Christian, P., 2004. A simple and effective evolutionary algorithm for the vehicle Hoen, K.M.R., Tan, T., Fransoo, J.C., et al., 2010. Effect of Carbon Emission Regulations
routing problem. Comput. Oper. Res. 31 (12), 1985e2002. on Transport Mode Selection in Supply Chains[R]. Retrieved from the Eind-
Dantzig, G., Ramser, J., 1959. The truck dispathing problem. Manag. Sci. 6 (1), 80e91. hoven University of Technology website: http://alexandria.tue.nl/repository/
Demir, E., Bektas, T., Laporte, G., 2012. An adaptive large neighborhood search books/672727.pdf.
heuristic for the pollution-routing problem. Eur. J. Oper. Res. 223 (2), 346e359. Hosseini-Nasab, H., Lotfalian, P., 2017. Green routing for trucking systems with
Dongarra, J., 2006. Performance of Various Computers Using Standard Linear classification of path types. J. Clean. Prod. 146, 228e233.
Equations Software. University of Tennessee. ReportCS-89-85. Jabir, E., Panicker, V.V., Sridharan, R., 2017. Design and development of a hybrid ant
Eshtehadi, R., Fathian, M., Demir, E., 2017. Robust solutions to the pollution-routing colony-variable neighbourhood search algorithm for a multi-depot green
908 J. Li et al. / Journal of Cleaner Production 201 (2018) 896e908

vehicle routing problem. Transport. Res. Transport Environ. 57, 422e457. pollution routing problem. Int. J. Prod. Econ. 176, 143e153.
Leggieri, V., Haouari, M., 2017. Lifted polynomial size formulations for the homo- Taillard, E., 1999. A heuristic column generation method for the heterogeneous fleet
geneous and heterogeneous vehicle routing problems. Eur. J. Oper. Res. 263 (3), VRP. Oper. Res. 33, 1e34.
755e767. Tajik, N., Moghaddam, R.T., Vahdani, B., Mousavi, S.M., 2014. A robust optimization
Li, J., Zhang, J.H., 2014. Study on the effect of carbon emission trading mechanism on approach for pollution routing problem with pickup and delivery under un-
logistics distribution routing decisions. Syst. Eng.-Theory Pract. China 34 (7), certainty. J. Manuf. Syst. 33 (2), 277e286.
1779e1787. Tarantilis, C., Kiranoudis, C., Vassiliadis, V., 2004. A threshold accepting meta-
Li, X.Y., Tian, P., Aneja, Y.P., 2010. An adaptive memory programming metaheuristic heuristic for the heterogeneous fixed fleet vehicle routing problem. Eur. J. Oper.
for the heterogeneous fixed fleet vehicle routing problem. Transport. Res. Part E Res. 152 (1), 148e158.
46, 1111e1127. Turkensteen, M., 2017. The accuracy of carbon emission and fuel consumption
Li, X.Y., Stephen, C.H.L., Tian, P., 2012. A multistart adaptive memory-based Tabu computations in green vehicle routing. Eur. J. Oper. Res. 262 (2), 647e659.
search algorithm for the heterogeneous fixed fleet open vehicle routing prob- Wu, J., Zhou, Z.X., 2015. A mixed-objective integer DEA model. Ann. Oper. Res. 228
lem. Expert Syst. Appl. 39 (1), 365e374. (1), 81e95.
Lin, C., Choy, K., Ho, G., Chung, S., Lam, H., 2014. Survey of green vehicle routing Wu, J., Zhou, Z.X., Liang, L., 2010. Measuring the performance of nations at Beijing
problem: past and future trends. Expert Syst. Appl. 41, 1118e1138. Summer Olympics using integer-valued DEA model. J. Sports Econ. 11 (5),
Niu, Y.Y., Yang, Z.H., Chen, P., Xiao, J.H., 2018. Optimizing the green open vehicle 549e566.
routing problem with time windows by minimizing comprehensive routing Xiao, Y.Y., Abdullah Konak, A., 2017. A genetic algorithm with exact dynamic pro-
cost. J. Clean. Prod. 171, 962e971. gramming for the green vehicle routing & scheduling problem. J. Clean. Prod.
Olivera, A., Viera, O., 2007. Adaptive memory programming for the vehicle routing 167, 1450e1463.
problem with multiple trips. Comput. Oper. Res. 34 (1), 28e47. Zhang, J., Zhao, Y., Xue, W., Li, J., 2015. Vehicle routing problem with fuel con-
Soleimani, H., Chaharlang, Y., Ghaderi, H., 2018. Collection and distribution of sumption and carbon emission. Int. J. Prod. Econ. 170, 234e242.
returned-remanufactured products in a vehicle routing problem with pickup Zhou, M., Jin, H., Wang, W.S., 2016. A review of vehicle fuel consumption models to
and delivery considering sustainable and green criteria. J. Clean. Prod. 172, evaluate eco-driving and eco-routing. Transport. Res. Transport Environ. 49,
960e970. 203e218.
Sundarakani, B., de Souza, R., Goh, M., et al., 2010. Modeling carbon footprints Zhu, X.Y., Garcia-Diaz, A., Jin, M.Z., Zhang, Y., 2014. Vehicle fuel consumption
across the supply chain [J]. Int. J. Prod. Econ. 128 (1), 43e50. minimization in routing over-dimensioned and overweight trucks in capaci-
Suzuki, Y., 2016. A dual-objective metaheuristic approach to solve practical tated transportation networks. J. Clean. Prod. 85, 331e336.