Heterogeneous fixed fleet vehicle routing problem based on fuel and carbon emissions

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: Received 1 February 2018 Received in revised form 6 July 2018 Accepted 7 August 2018 Available online 13 August 2018

In this paper, we study an emission-based heterogeneous fixed fleet vehicle routing problem (E-HFFVRP) with considerations of fuel and carbon emissions. This problem involves routing a fleet of a fixed number of vehicles with various capacities to serve a set of customers. It seeks to minimize the objective function, which incorporates the fixed expenses and variable costs consisting of fuel consumptions and carbon 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 to calculate fuel and carbon emissions. Moreover, a split-based adaptive tabu search (SATS) algorithm using an optimal split scheme and an adaptive tabu search algorithm is proposed. Its key features and components are designed accordingly. Results of numerical experimentations on two sets of generated instances confirm the efficiency and effectiveness of the algorithm. © 2018 Elsevier Ltd. All rights reserved.

Keywords: Vehicle routing Freight transportation Heterogeneous fleet Fuel Carbon emissions

1. Introduction In recent years, the energy crisis and global warming have aroused great concern in the world. To deal with this challenge, most countries in the world have enacted related carbon polices such as carbon cap and trade, carbon tax, and carbon quota. EU has implemented carbon cap and trade since 2005 and imposed carbon tax from 2012. Thus, how to effectively save energy and reduce carbon emissions has become the urgent task to address the world's sustainable development. Transportation has been a major source of carbon emissions, given the rapid industrial and economic development. In China, the carbon emissions of freight transportation have already accounted for 30% of the total carbon emissions and the fuel consumption of freight companies is at least 40% of their total costs (Li and Zhang, 2014; Li and Zhang, 2014). With strict carbon polices and increasing transportation cost, many companies aim to reduce the carbon emissions and fuel consumptions by optimizing vehicle routing decisions. For example, Walmart reduced its transport energies and carbon footprint by

* Corresponding author. 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.

about 20%e30% through improved routing decision on when and how to visit its distributed stores. Therefore, it is necessary to incorporate fuel and carbon emissions into classical vehicle routing problem (VRP) to study a new variant and formulation. The classical VRP (Dantzig and Ramser, 1959; Dantzig and Ramser, 1959) focuses on how to minimize the total travel distance or total costs of all vehicles, which generally has a linear relationship with distance. However, the routing arrangement of minimizing the travel distance does not necessarily derive the optimal solution of minimizing fuel and carbon emissions, because the fuel and carbon emissions of vehicles are not only related to distance, but also to various factors such as vehicle weight, speed and load Turkensteen (2017). In our paper, we study an emission-based heterogeneous fixed fleet vehicle routing problem (E-HFFVRP), which extends the traditional VRP by taking fuel and carbon emissions into account. In particular, there is a heterogeneous fleet of vehicles to conduct delivery from one depot to all customers. After finishing the delivery, 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 fixed cost. The costs of carbon emissions involve the costs of carbon

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

pollution and emission reduction under various environmental policies, such as carbon tax and carbon cap-and-trade mechanism. Vehicle fixed costs are mainly related with the costs of hiring, maintenance, depreciation and insurances, etc. The task of the decision-maker is how to arrange the fleet and routing of these different vehicles to achieve the minimization of total transport network costs. Our study is related to green vehicle routing problem (GVRP), which is an alternative to regular well-established vehicle routing models. Literature reviews on GVRP are presented in the papers such as Lin et al. (2014) and Zhou et al. (2016). GVRP has been used recently for purposes such as making routing decisions to minimize fuel consumptions Zhu et al. (2014); Hosseini-Nasab and Lotfalian (2017) or reduce carbon emissions Niu et al. (2018); Xiao and Abdullah Konak (2017), recovery and distribution of remanufactured products considering sustainability and green criteria Soleimani et al. (2018) or routing a capacitated multi-depot distribution system Jabir et al. (2017). As a typical GVRP, the pollution routing problem (PRP) is also an extension of the traditional VRP and has a broad objective function that accounts for carbon emissions. It was first introduced by Bektas and Laporte (2011), and then some effective approaches are developed to solve this problem, such as adaptive large neighborhood search Demir et al. (2012) and disjunctive convex programming Fukasawa et al. (2016). Some other variants of PRP are paid more attention in the literature such as the time-dependent PRP Franceschetti et al. (2017), PRP with multiple conflicting objective functions Suzuki (2016), pickup and delivery PRP Tajik et al. (2014) and PRP with uncertain demand and travel time (Eshtehadi et al. (2017). However, in practical business activities, outsourcing has been widely used. Outsourcing carriers may have limited number of heterogeneous vehicles, each type with different capacities, fixed costs, fuel consumptions and carbon emissions. This is so called EHFFVRP, which is to study how to arrange the routing of fixed number of vehicles from the depot to a series of customers to minimize a total cost consisting of fixed costs, fuel and carbon emissions. Another related problem is heterogeneous fixed fleet vehicle routing problem (HFFVRP), which routes a fixed number of vehicles from the depot to visit all customers with minimizing the overall fixed and variable costs. To address HFFVRP, different heuristics and metaheuristics have been developed, such as adaptive memory programming metaheuristic Li et al. (2010), tabu search algorithm ~o (2011), and lifted polynomial size formulations Leggieri Branda and Haouari (2017). However, different with our problem, HFFVRP did not consider the effect of fuel and carbon emissions on vehicle routing arrangement. One closely related literature is from Zhang et al. (2015), whose focus is on a capacitated vehicle routing problem to minimize the sum of fuel cost, carbon emission cost and the vehicle usage cost. By contrast, our study is centered on an E-HFFVRP to minimize not only the variable costs of fuel and carbon emissions but also fixed costs. Hence, the study of E-HFFVRP in this paper is a new extension of VRP in (Zhang et al. (2015), and is quite different from them as well. To the best of our knowledge, E-HFFVRP has not yet attracted much attention in the extant literature and an appropriate method that addresses this problem remains unexplored. Existing literature review shows that pure HFFVRP without fuel and carbon emissions is well studied. However, low-carbon transportation and logistics are paid more attention recently by governments and customers. To fill this gap, in this paper, we present the E-HFFVRP that integrating both economic and environmental impacts. To reduce fuel consumptions and carbon emissions of vehicle routing, we use a comprehensive approach to calculate the

897

amount of fuel and carbon emissions. Traditional approaches usually use a linear function with respect to travel distance or vehicle load. In contrast, our approach takes into account the effects of not just travel distance but also lots of other factors, such as vehicle weight, load and speed. Moreover, we propose a mixed integer programming model for the E-HFFVRP, and adapt adaptive memory programming to develop a split-based adaptive tabu search (SATS) algorithm. A few promising characteristics of this algorithm are introduced to improve its computational performance. In particular, the algorithm incorporates an optimal split scheme and an adaptive tabu search algorithm. The optimal split scheme is applied to generate the adaptive memory, and the adaptive tabu search algorithm is used to improve each temporary solution. The remainder of this paper is organized as follows. Section 2 introduces an approach of calculating the fuel and carbon emissions, and then describes the E-HFFVRP model. In Section 3, the algorithm named SATS is presented. Its main features and different components are given in detail. In Section 4, results of numerical studies on two sets of test problems are shown to illustrate the efficiency of the proposed algorithm and analyze the effect of split strategy. Finally, Section 5 summarizes the final conclusions. 2. Mixed integer programming model 2.1. Problem description The notations of E-HFFVRP are given as follows. We assume that G ¼ ðV; EÞ is a digraph, V ¼ f0; 1; /; ng is the nodes set, and E ¼ fði; jÞji; j2V; isjg is the arcs set between each pair of nodes. The depot denoted by node 0 contains different types of vehicles, and other n nodes constitute a customer set C ¼ f1; 2; /; ng. x ¼ f1; 2; /; Kg is the set of vehicles and j ¼ f1; 2; /; Mg is the set of vehicle types. For each m2j, the vehicle capacity is Qm , and the fixed cost is Fm . The fuel cost per unit is Cf , and the emission cost per unit fuel is Ce . The number of vehicle type m is fixed and equals P to Nm . Therefore, the total number of vehicles is M m¼1 Nm . For each customer i2C, its demand to be delivered from the depot is qi . For each arcði; jÞ2E, dij is assumed to be a non-negative distance. Given above conditions, the objective function of the E-HFFVRP is to minimize the sum costs of variable expenses (fuel consumptions and carbon emissions) and fixed expenses, by determining a vehicle scheduling strategy. Unlike the HFFVRP, which only focuses on the economic costs, the E-HFFVRP considers both the economic and environmental costs (i.e. fuel cost and carbon emission cost) in its objective function. 2.2. Calculation of fuel and carbon emission costs Carbon emissions of vehicles depend on a number of factors Sundarakani et al. (2010), such as transport mode, fuel type, fuel consumption, and travel distance. For practical reasons, we introduce an emission calculation approach for road transport from Hoen et al. (2010). The total carbon emissions on an arc ði; jÞ are assumed to be EMij , which is related to FCij and the fuel emission factor (FE). Thus, it can be calculated as follows:

EMij ¼ FE$FCij

(1)

where FE is defined as a gram of CO2 emitted per liter of fuel, i.e., the efficiency of converting the fuel into carbon emissions, which is affected by vehicle type, fuel type and other factors. The traditional calculation of fuel consumption is a linear function of travel distance, which is not in line with industry applications. 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 shown in the following.

i h   FCij ¼ aij w þ zij þ bv2ij dij

(2)

where vij is a vehicle's average speed on an arc ði;jÞ; L ¼ w þ zij is the 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 related to acceleration, road angle, etc.; and b is a vehicle-specific constant that depends on the vehicle frontal surface area, air density, and rolling resistance. In terms of Eq. (2), the total cost TCij that consists of fuel consumptions and carbon emissions is shown as follows:

  TCij ¼ cf þ ce FCij

(3)

Table 1 summarizes some symbols and values (ranges) used in the E-HFFVRP model and experiments. To be consistent with practical applications, these data are derived from one large logistics company in China and some related literature Bektas and Laporte (2011); Zhang et al. (2015).

2.3. Model formulation Given above analysis and integer programming approaches Wu et al. (2010); Wu and Zhou (2015), we formulate E-HFFVRP as a mixed integer programming model below:

min

M X K X n n  X X

Cf þ Ce

h





i

aij w þ zij þ bv2ij Þdij xm ijk

n X K X

xm 0jk  Nm ; cm

(9)

j¼1 k¼1 M X K X

qj xm ijk  zij 

m¼1 k¼1

X

xm ijk

i;j2SS

M X K X ðQm  qi Þxm ijk ; ci; j

(10)

m¼1 k¼1

       S  1; S3V\f0g; Ssf; cm; k  

h m xm ijk 2f0; 1g; yik 2f0; 1g; zij 2 0; þ∞Þ; ci; j; m; k

(11)

(12)

The objective function (4) has two parts. The first part measures the total variable costs of fuel consumptions and carbon emissions. The second part measures the total fixed costs for vehicles. Constraint (5) ensures that every customer can only be visited by a vehicle once. Constraint (6) and (7) define that every customer enters and leaves one vehicle exactly once. Constraint (8) is the capacity limitation of vehicles. Constraint (9) imposes a limitation on the number of each vehicle type. Constraint (10) guarantees that the vehicle's load on an arc ði; jÞ cannot exceed its remaining capacity. Constraint (11) is used to eliminate sub-loops. Constraint (12) defines xm , ym , and Zij as decision variables; Given a vehicle k of ijk ik ¼ 1 if and only if the vehicle of type m type m and arc ði; jÞ, xm ijk ¼ 0, Given a vehicle k of type passes through arc ði;jÞ, otherwise, xm ijk ¼ 1 if and only if this vehicle visits node i, otherm and node i, ym ik wise it is 0. Given an arc ði; jÞ, the load zij  0 if it is visited by a vehicle, otherwise zij ¼ 0.

m¼1 k¼1 i¼0 j¼0;jsi

þ

M X m¼1

Fm

K X n X

3. Metaheuristic for the E-HFFVRP

xm 0jk

k¼1 j¼1

(4)

s:t:

M X K X

ym ik ¼ 1; cis0

(5)

m¼1 k¼1 n X

m xm ijk ¼ yjk ; cj; m; k

(6)

m xm ijk ¼ yik ; ci; m; k

(7)

qi ym ik  Qm ; ck; m

(8)

i¼1 n X j¼1 n X i¼1

Table 1 Symbols representation referred in the E-HFFVRP model. Symbols

Descriptions

Values (ranges)

cf ce FE vij

Fuel cost per liter (CNY) Emissions cost per liter (CNY) Fuel emissions factor (kilogram/liter) Average arc speed (kilometers/hour) Arc specific constant

7.30 0.64 2.32 40e100 0.09e0.15

aij

Since the E-HFFVRP is a new extension of the classical VRP, a 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 adaptive tabu search (SATS) algorithm that incorporates an optimal split scheme and embedded adaptive tabu search algorithm is proposed. To facilitate the presentation, we define the following notations: a is the adaptive memory that comprises some solutions; b is the set of temporary solutions constructed based on 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 solution for xi ; xcb is the best solution obtained in the current iteration; xgb is the global-best solution obtained from the start of the algorithm; and nb is the number of best solutions in the current iteration. In order to help the readers understand the whole picture of the proposed algorithm, the basic structure of the SATS is shown by a flow chart in Fig. 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 selected in probability to generate temporary solutions in each iteration step. Next, an adaptive tabu search (ATS) is used to improve each temporary solution. To better utilize the search memory and history, several best improved solutions are selected to update a. The above steps are repeated until the maximum number of iterations IM is satisfied. Finally, after the main loop, a post-optimization procedure (POP) is added to further improve the 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 transportation cost (condition (15)). j X

cði; jÞ2A :

qSk  Qm

(13)

k¼iþ1

Nam ¼ Nm  Num  1

(14)

cði; jÞ2A : wij 

¼ cf þ ce



FC0;Siþ1 þ

j1 X

! 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

Fig. 1. Flow chart of the SATS.

incorporates an optimal split algorithm and an adaptive tabu search as shown in Algorithm 1 in Appendix. 3.1. Solution denotation and split 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) represents the sub-route (0-1-2-3-0) and (0-4-5-0). Instead, we simply denote a solution tour as a sequence S composed of n customer without delimiters, such as the above tour (0-1-2-3-0-4-5-0) becomes (1-2-3-4-5). This denotation approach can transform VRP into a simpler traveling salesman problem. On the other hand, it is essential to split S into sub-routes to evaluate the solution. Some 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 fuel and carbon emissions. Accordingly, we develop an improved optimal splitting procedure called Optimal Split (OS) to derive the best E-HFFVRP solution corresponding to the denotation sequence. The main idea of OS is that for any kind of optimal tour split, the remaining tour split is also optimal when subsequent several subroutes are removed. OS works on an auxiliary graph H ¼ ðN; A; WÞ, where N includes n þ 1 nodes indexed from 0 to n; A contains route 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 number limitation (condition (14), where Nam is the number of type m vehicles available, Num is the number of type m utilized vehicles);

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 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 not belong to the solution. The number of routes used is computed by the counter t. The insertion cost cki of route rðkÞ including customer i is calculated by the following formula:

   cki ¼ cf þ ce FCi1 ;i þ FCi;i2  FCi1 ;i2

(16)

where i1 and i2 are the predecessor node and successor node for a given insertion location, respectively. Each un-routed customer in R is put into one of the available routes with minimum insertion cost. As the vehicle capacity is not considered, infeasible solutions may be generated. As a diversification strategy, SATS allows cross-border search between feasible and infeasible spaces to improve the quality of the final solutions. 3.2. Solution evaluation SATS allows exploring the infeasible solutions, which means that there might exist a solution of violating the vehicle capacity. We define a penalty function FðxÞ to evaluate a solutionx.

FðxÞ ¼ f ðxÞ þ

K X

pli

(17)

i¼1

where f ðxÞ represents the original objective value of solution x; K denotes the number of vehicles used; li is the load excess in vehicle (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 generation In the initial stage of SATS, a is composed of various solutions. Let am be the size of a. We generate random permutation of customer nodes as initial solutions. Then, these solutions are arranged in descending order and put in a. First, Each solution x is assigned with a triple ðfðxÞ; f ðxÞ; FðxÞÞ, where fðxÞ is the total load excess value of the capacity in x. Second, there are three possible 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 level than x2 ; (3) when x1 and x2 are both infeasible, x1 has a higher level than x2 if and only if Fðx1 Þ < Fðx2 Þ. Note that if f ðx1 Þ ¼ f ðx2 Þ or Fðx1 Þ ¼ Fðx2 Þ, x1 and x2 is chosen randomly to be assigned with a higher level. A multi-start strategy is used to enhance diversification and thus increase the probability of generating better solutions, that is, nb solutions in a are probabilistically selected to be temporary solutions at each iteration. As mentioned above, a solution with a higher level is more likely to be selected in probability. Hence, the probability of selecting the ith solution is defined as pi ¼ ðns  i þ 1Þ=n2s , 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 has least selection probability. 3.4. ATS for solution improvement An improvement-based procedure is applied to improve all solutions in b and push to their local optima. We use an ATS as a solution improvement procedure in SATS. Some related notations in the ATS are listed as follows: i0 represents the iteration counter of no improvement in current best feasible or infeasible solution; iL is the limit of i0 in the current cycle; q is the tabu tenure; d is the number of nearest neighbors of a customer; IA is the total limitation number of iterations; FS FD ,FW and FO are the frequencies of the neighborhood operators of single, double, swap 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 tenure and the frequencies of the neighborhood operators. In particular, once the search reaches to the middle of the cycle (iL =2) or the end of the cycle (iL ), accordingly, the ATS changes some parameters (qdFS FD FW FO ) to enhance the diversification. POP and restart strategy are also used to achieve the intensification. The outline of the ATS is described in Algorithm 3 in Appendix. 3.4.1. Neighborhood structure In a tabu search algorithm, the degree and quality of solution search are determined by neighborhood structure. We use three intra-route neighborhood operators (single insertion, double insertion and swap) and one inter-route neighborhood operator (2opt). We limit the search of these operators by d  neighborhood restriction to reduce the number of potential moves at each iteration, which thus not only saves the computing time but also drives the search to find better solutions. For any node s2S, we define its 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 and used by Gendreau et al. (1994). However, the main difference with them is that the d in our algorithm varies with the search evolution, that is, the initial d is very small and increases gradually during the search evolution. The single insertion removes a customer Si from one tour and inserts it into various positions in another tour that contains at least

one of the Nd ðSi Þ. The double insertion removes any pair of customers ð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 Þ. The swap allows a customer node Si from one tour to be exchanged with a customer node Sj from another tour while satisfying that at least one customer in route j falls inNd ðSi Þ. For the 2-opt, nonadjacent arcs ðSi ; Siþ1 Þ and ðSj ;Sjþ1 Þare replaced by ðSi ; Sj Þ and ðSiþ1 ; Sjþ1 Þ in the same tour, and then the direction of the sub-tour ðSiþ1 ; /; Sj Þ which has at least one of the Nd ðSi Þ or Nd ðSjþ1 Þ is reversed. 3.4.2. Tabu status and tabu tenure To avoid local cycling, we don't allow recent solutions found to be used in the next q iterations and put some relevant attributes of these solutions into a tabu list. For each action that generates a new solution, the inverse action will be taboo in subsequent q iterations and is then inserted into the tabu list. The status of an operator is required that a customer removed from a position of a tour cannot go back to it during the tabu tenure. The tubu tenure is changed systematically as the search evolution. Due to the fact that the tabu restriction may prohibit promising moves, a common aspiration criterion is used if the move action produces a better solution.

3.4.3. Restart and shaking As tabu search algorithm is easy to fall into local optimal, to address this problem, restart and shaking is often used to obtain better solutions. By changing parameters, restart can make a good tradeoff between intensification and diversification. The restart of our algorithm is done from the best-so-far solution xgb if no improvement is found within iL iterations. Furthermore, to escape local optimum at the end of the search, it is necessary to shake the best known solution and reach better solutions. Shaking technique is often used when the search falls into local optimum within a given neighborhood structure. In our implementation, a random move is performed on the current neighborhood by randomly executing one of the three movements: GENI, US Gendreau et al. (1994) and 2-opt.

3.5. Adaptive memory update After the solution improvement procedure ATS is implemented, several poor solutions are removed from a and replaced by the improved solutions. First, we select nb (nb < nb ) best solutions from nb improved solutions according to the following conditions: (1) The feasible solution should be selected in advance when feasible and infeasible solutions are all available; (2) When only feasible solutions are derived, we select the one with the least objective value; and (3) When only infeasible solutions are generated, we select the one with the least penalty value. Then, these solutions are assigned and sorted with the triple given in Section 3.3. Finally, the nb worst solutions in a are deleted and replaced by the same number of best solutions with assigned order level.

3.6. POP At the end of the ATS or SATS, in order to promote intensification, we use a POP to further improve the best found solution. Similar to those operators in the shaking stage, GENI, US, and 2-opt are also executed. However, they are not used by random selection, instead they are implemented in sequence using the best-accepted criterion, i.e., they search the entire neighborhood of the solution space to find the best improvement.

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

901

dependent on number of customers. Some parameters will be varied adaptively with the search evolution of the algorithm to 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 reduced to allow intensification. However, at the end of the cycle, current solutions in the search space cannot be improved, and q should be added to drive the search toward more distant spaces. Therefore, for a cycle, we set q ¼ maxðq=2; 6Þ in the middle, and at the end, we change it to be q ¼ minð2q þ 3; 0:6nÞ. Similar with the setting of the tabu tenure, we change d with the search evolution as well. A set of experiments showed that the initial d must be small. Here, it is set as d ¼ 1. In the middle of the cycle, it is set as d ¼ minðd þ 1;n  1Þ. And at the end of the cycle, no better solution could be found. To enhance diversification in the search space, d is desirable to increase, and thus we set it as d ¼ minðd þ 2;n  1Þ. According to our previous experiments, the value of p(penalty term) was set as follows: At the initial stage, it equals to Pn i¼1 qi =100n. If ten successive iterations all generate infeasible solutions, then p ¼ 2p. If ten successive iterations all generate feasible solutions, then p ¼ p=2. Finally, we study how to set appropriate frequencies of four neighborhood operators. According to the degree of conversion to a solution, the four neighborhood operators are sorted from low to high as follows: single insertion, swap, double insertion, and 2-opt. With multiple iterations, single insertion can result in the same conversion caused by the other two operators (swap and double). Therefore, single insertion should be used most frequently. Since double insertion and swap can produce similar conversions, they should not appear in the same iteration. Through a set of trials, we 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, and FO ¼ 6. According to our experiments, good values used for other parameters are set and summarized as follows:

4. Numerical experiments 4.1. Benchmark problems SATS is programmed in C language and executed on a PC with Intel Duo CPU 2.26G, 2G of RAM, under Windows 7. This paper is the first attempt to use instance to test the E-HFFVRP. The sets of benchmark problems in existing papers are unable to be used. Hence, we plan to test the performance of SATS by designing two sets of benchmark problems, whose configurations for the EHFFVRP are shown in Tables 1 and 2. Set 1 is generated in Table 2 based on Taillard (1999), which comprises eight problems labeled from 13 to 20. Six vehicle types labeled from A to F are given. In Set 1, we generate the number of customers according to the integers between 50 and 100. We also generate a new set of benchmark problems, i.e., Set 2 in Table 3, which are composed of eight problems denoted by P1eP8 with the number of customers ranging from 150 to 300. All customers for each problem are located in a Euclidean plane. The coordinates ðxi ; yi Þ of customer i (i ¼ 1; 2; /; n) are randomly produced in square interval ½0; 1002 in a uniform distribution, where (30, 60) represents the depot. We set the demand of a customer as any integer randomly selected in interval (Zhou et al., 2016; Christian, 2004). The fleet contains six types of vehicles: A, B, C, D, E, and F. In the headings of Tables 2 and 3, VT and VP represent the vehicle type and property respectively. For vehicles of type i, bi is a vehicle-specific constant, and r0 is the ratio of total demand and total capacity, i.e., (total demand/total capacity)  100. According to the characteristics of these two sets of problems, they have important implications in the resolution.

4.2. Parameter tuning In SATS, how to set some key parameters has a great impact on the algorithmic performance. Some characteristics in ATS are in line ~o (2011). Thus, with that in an algorithm for solving HFFVRP Branda previous knowledge is used to tune the value of several ATS parameters, such as qd, and the frequencies of the neighborhood operators. 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 preliminary experiments (Birattari, 2005) (Birattari, 2005). The tuning results indicated that the parameter combination settings work well in our proposed algorithm. In order to ensure that the SATS performance is consistent with the size of the problem, most of the parameters are set to be

pffiffiffi (1) IA ¼ 100 n; (2) Initially iL ¼ 1.5n, and at the end of the cycle, iL ¼ iL þ0.5n; pffiffiffi pffiffiffi pffiffiffi (3) am ¼ 10n, nb ¼ n, nb ¼ n=2, IM ¼ 3000 n.

4.3. Numerical results 4.3.1. Effect of split strategy In this paper, we use a novel optimal split strategy in the proposed algorithm, SATS. We conduct an experiment based on the given parameters of configuration tuning to assess the quality of

Table 2 Data for the testing problems of Set 1. VT

VP

n r0 A

QA FA NA

B

QB FB NB

C

QC FC NC

bA

bB

bC

Problem

VT

13

14

15

16

17

18

19

20

50

50

50

50

75

75

100

100

95.39

88.45

94.76

94.76

95.38

95.38

76.74

95.92

20 20 4 0.54 30 35 2 0.70 40 50 4 0.73

120 100 4 2.67 160 1500 2 3.13 300 3500 1 3.75

50 100 4 1.59 100 250 3 2.60 160 450 2 3.13

40 100 2 0.73 80 200 4 2.11 140 400 3 3.01

50 25 4 1.59 120 80 4 2.67 200 150 2 3.37

20 10 4 0.54 50 35 4 1.59 100 100 2 2.60

100 500 4 2.60 200 1200 3 3.37 300 2100 3 3.75

60 100 6 1.84 140 300 4 3.01 200 500 3 3.37

VP

Problem 13

14

15

16

17

18

19

20

n

50

50

50

50

75

75

100

100

88.45

94.76

94.76

95.38

95.38

76.74

95.92

350 320 1 4.05

150 180 2 3.05 250 400 1 3.70 400 800 1 4.14

r0

95.39

D

QD FD ND

E

QE FE NE

F

QF FF NF

70 120 4 1.91 120 225 2 2.67 200 400 1 3.37

bD

bE

bF

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

P1

P2

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

80 300 4 2.11 120 350 4 2.67 150 500 3 3.05

100 350 3 2.60 160 500 3 3.13 300 1500 3 3.75

150 400 4 3.05 200 550 3 3.37 300 1200 3 3.75

100 320 2 2.60 160 450 4 3.13 250 800 3 3.70

60 200 6 1.84 120 350 4 2.67 200 510 2 3.37

120 340 5 2.67 160 500 5 3.13 200 560 4 3.37

200 600 3 3.37 250 800 3 3.70 300 1500 3 3.75

150 410 5 3.05 200 600 3 3.37 250 850 3 3.70

A

QA FA NA

B

QB FB NB

C

QC FC NC

bA

bB

bC

P3

P4

P5

P6

P7

VP

P8

the optimal split strategy. The experiment is conducted by performing four versions of SATS with different split strategies as follows: (1) SSA-SATS, a SATS with a simple split strategy in which available vehicles of different types are assigned in ascending order of capacity; (2) SSD-SATS, another simple split strategy in which available vehicles of various types are assigned in descending order of capacity; (3) RS-SATS, where a random split strategy is implemented, i.e., available vehicles of different types are assigned at 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 the best solutions in the column labeled “Average” along with the related CPU average time used in computing the algorithm. For the last three columns, the column entitled DevA (%), DevD (%) and DevR (%) shows the percentage deviations of OS-SATS from SSASATS, SSD-SATS, and RS-SATS, respectively. Furthermore, suppose SðOSÞ is the optimal solution found with OS-SATS, and SðADRÞ represents the solution value produced by one of the other three split-based algorithms. Then, the percentage deviation is calculated as 100ðSðADRÞ  SðOSÞÞ=SðADRÞ. As shown in Table 4, from both solution quality and computational efforts, the optimal split strategy outperforms the other split strategies. Moreover, the optimal split strategy is able to discover the best solution in nearly all instances. The average deviations from the other three split strategies are significant from 4.98% to

Problem P1

P2

P3

P4

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

D

QD FD ND

400 2200 1 4.14

E

QE FE NE

350 1500 2 4.05 400 2100 2 4.14

300 1200 2 3.75 350 1800 1 4.05

300 1100 2 3.75 350 1600 1 4.05 400 2100 1 4.14

300 1200 2 3.75 400 2100 2 4.14

QF FF NF

250 750 3 3.70 300 1200 2 3.75 350 1600 2 4.05

350 1800 2 4.05 400 2100 1 4.14

F

200 550 3 3.37 250 600 1 3.70 350 2200 1 4.05

bD

bE

bF

P5

P6

P7

P8

10.11% for Set 1 instances and from 2.93% to 12.92% for Set 2 instances. In particular, the average deviation from the random split strategy is above 10%. Moreover, when the split strategy is used, the computational time required to obtain the best solutions is not significantly increased. These figures suggest that the optimal split strategy can be used to obtain good-quality solutions. Therefore, the optimal split strategy should be adopted in the SATS to derive a better balance between computational time and solution quality.

4.3.2. Comparison of SATS with other algorithms In this subsection, we carry out simulation experiment to show the performance of SATS by comparing it with some existing algorithms. For SATS, we also use the best parameter settings selected in the above sections. The selected algorithms include a pure tabu search and two methods used in HFFVRP as follows: (1) TS: a tabu search algorithm, which is similar to ATS, except that it use saving algorithm rather than the split-based procedure to construct initialization solution. Its termination condition is also the maximum number of iterations, pffiffiffi 3000 n. (2) BATA: a back-tracking adaptive threshold accepting metaheuristic Tarantilis et al. (2004), which was implemented on MS Visual Cþþ, Pentium II, 400 MHz, 128 MB RAM.

Table 4 The effect of split strategy in SATS. Instances

13 14 15 16 17 18 19 20 Average P1 P2 P3 P4 P5 P6 P7 P8 Average

SSA-SATS

SSD-SATS

RS-SATS

OS-SATS

Average

CPUs

Average

CPUs

Average

CPUs

Average

CPUs

2081.12 11121.32 2604.07 2675.01 2050.64 2921.95 11248.84 4474.95

100 107 113 98 149 153 236 215

2085.96 10192.37 2591.58 2658.66 1991.98 2859.11 11455.49 4493.63

109 105 108 104 153 160 227 208

2211.47 11444.35 2690.05 2652.51 2192.88 3178.75 11778.69 4538.69

113 118 110 107 161 154 231 209

2065.96 10034.51 2520.39 2631.22 1825.06 2723.14 10215.39 4060.80

95 98 90 89 143 151 225 206

12822.67 14386.42 21527.50 21384.31 23439.18 24509.41 26428.75 31338.07

314 326 447 412 491 524 655 673

15005.69 16320.36 23377.88 22302.72 25400.64 26456.88 28060.79 34714.78

290 313 435 405 462 482 624 637

15296.87 16318.78 23412.30 24760.75 24239.18 27030.09 28991.76 35043.96

285 328 427 410 481 497 628 649

12950.24 13668.71 21273.17 20837.67 21658.77 23053.85 25985.84 31175.32

283 307 416 397 446 468 611 633

DevA (%)

DevD (%)

DevR (%)

0.73 9.77 3.21 1.64 11.00 6.80 9.19 9.25 6.45 0.99 4.99 1.18 2.56 7.60 5.94 1.68 0.52 2.93

0.96 1.55 2.75 1.03 8.38 4.76 10.83 9.63 4.98 13.70 16.25 9.00 6.57 14.73 12.86 7.39 10.20 11.34

6.58 12.32 6.31 0.80 16.77 14.33 13.27 10.53 10.11 15.34 16.24 9.14 15.84 10.65 14.71 10.37 11.04 12.92

110/5500 34/1700 46/2300 99/4950 148/7400 119/5950 287/14350 200/10000

Best 843/6744 387/3096 368/2944 341/2728 363/2904 971/7768 428/3424 1156/9248 1990.05 9755.32 2362.29 2510.63 1739.19 2472.82 10210.61 4037.54 120/6600 129/7095 115/6325 122/6710 167/9185 175/9625 283/15565 271/14905 1986.33 9757.88 2376.56 2541.86 1738.22 2488.68 10213.57 4035.25 103/5665 112/6160 97/5335 101/5555 158/8690 163/8965 241/13255 238/13090 1920.31 9670.36 2358.49 2475.85 1683.63 2458.04 10194.17 4012.80 17 6 9 9 10 12 8 13 50 50 50 50 75 75 100 100 13 14 15 16 17 18 19 20

17 7 9 9 11 14 10 13

CPUs/Tflop Best

BATA

Best

CPUs/Tflop TS

Best

CPUs/Tflop SATS Vehicles NV n Instances

Table 5 Comparison results of different algorithms on benchmark instances.

903

(3) MAMP: a multistart adaptive memory programming metaheuristic (Li et al., 2012; Li et al., 2010), which was implemented on MS Visual Cþþ, Intel CPU 2.2G.

1945.06 9710.34 2362.29 2473.69 1731.02 2467.35 10194.17 4015.58

MAMP

CPUs/Tflop

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

As the algorithms of BATA and MAMP have been tested on benchmark problems Set 1, we use Set 1 to conduct an experiment for all algorithms. Note that different algorithms are executed on different computer environments. To offer a fair comparison, we measure the computational speed in ten millions of floating-point operations per second (Tflop/s) based on Dongarra tables Dongarra (2006) as follows: Pentium IId8 Tflop/s, 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 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 algorithm (Best) and the best solution found by all these four algorithms (BSA). The percentage deviation is computed as 100(Best eBSA)/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 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 SATS are fewest in available vehicles for many instances, such as 14, 17 and 18. Finally, on average, the computing time for SATS is 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 carbon emissions in industry operations. 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, logistic service providers should bring environmental factors (such as fuel and carbon emissions) into the routing arrangement to carry out an environment-friendly transportation. The proposed 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 that there exists tradeoff relationship between economic costs and carbon emissions. For the logistic service providers, it is valuable to use the proposed method to reduce the carbon Table 6 The percentage deviations of each algorithm from the best solution of all algorithms. Instances

BSA

SATS

TS

BATA

MAMP

13 14 15 16 17 18 19 20 Average

1920.31 9670.36 2358.49 2473.69 1683.63 2458.04 10194.17 4012.80

0 0 0 0.0872 0 0 0 0 0.0109

3.3237 0.8969 0.7603 2.6819 3.1406 1.2312 0.1899 0.5563 1.5976

3.5044 0.8709 0.1609 1.4713 3.1946 0.5977 0.1610 0.6127 1.3217

1.2725 0.4117 0.1609 0 2.7377 0.3773 0 0.0692 0.6287

904

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

emissions at the minimal costs. Finally, from this paper, the policymaker can know how to formulate and adopt effective governmental polices to affect the environmental performance of the logistic service providers. For example, sufficiently higher carbon price may reduce the carbon emissions at a lower degree, but depress the logistic firms by increasing more economical costs. 5. Conclusion In light of growing concerns about global warming, all modern industries are starting to increase their investment in reducing fuel consumptions and carbon emissions. To achieve this objective, the transportation industries should optimize vehicle routings by 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 and carbon emissions. A mixed integer programming model was built along with introducing an approach of calculating fuel and carbon emissions. A split-based adaptive tabu search algorithm is proposed, which combines an optimal split procedure with an adaptive tabu search. A vehicle routing is denoted as a sequence without delimiters. The split algorithm is then used to extract an optimal solution from the sequence. In addition, an adaptive tabu search (ATS) is embedded in SATS to improve the solution. We generated different sets of instances to evaluate the effectiveness of SATS. Then, some parameters are tuned systematically based on previous knowledge and F-Race procedure. We tested this algorithm on a series of problems. The results showed that optimal split strategy could be incorporated to enhance the performance of SATS, and the SATS can effectively search for high quality solutions on the

given instances. Although the SATS algorithm yields good results in the given test instances, some limitations should be noticed. Firstly, since the EHFFVRP is an NP-hard problem, if the data structure of the problems is large dimension and uniform, some additional difficulties 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 several features with the algorithm may be returned, such as tabu tenure, the frequencies of the neighborhood operators and shaking mechanism. Moreover, other efficient algorithms will be incorporated in the solution improvement to solve more difficult problems with other strong constraints. In the future, we shall test the performance of our algorithm through other large-scale real data. Acknowledgments This work was partly supported by a grant from the Natural Science Foundation of China (Grant No. 71571111); Zhejiang Provincial Philosophy and Social Sciences Foundation of China (Grant No. 18NDJC180YB); Zhejiang Provincial Natural Science Foundation of China (Grant No. LY17G020005), and the National Natural Science Foundation of China (NSFC) (Grant Nos. 71302035, 71672179). The authors also would like to thank the Qilu Young Scholars and Tang Scholars of Shandong University for financial and technical support. Appendix

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

References Bektas, T., Laporte, G., 2011. The pollution-routing problem. Transport. Res. Part B 45 (8), 1232e1250. Birattari, 2005. The Problem of Tuning Metaheuristics as Seen from a Machine Learning Perspective. PhD thesis. Universite Libre De Bruxelles. ~o, J.A., 2011. Tabu search algorithm for the heterogeneous fixed fleet vehicle Branda routing problem. Comput. Oper. Res. 38 (1), 140e151. Christian, P., 2004. A simple and effective evolutionary algorithm for the vehicle routing problem. Comput. Oper. Res. 31 (12), 1985e2002. Dantzig, G., Ramser, J., 1959. The truck dispathing problem. Manag. Sci. 6 (1), 80e91. Demir, E., Bektas, T., Laporte, G., 2012. An adaptive large neighborhood search heuristic for the pollution-routing problem. Eur. J. Oper. Res. 223 (2), 346e359. Dongarra, J., 2006. Performance of Various Computers Using Standard Linear Equations Software. University of Tennessee. ReportCS-89-85. Eshtehadi, R., Fathian, M., Demir, E., 2017. Robust solutions to the pollution-routing

907

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 for the time-dependent pollution-routing problem. Eur. J. Oper. Res. 259 (3), 972e991. Fukasawa, R., He, Q., Song, Y.J., 2016. A disjunctive convex programming approach to the pollution-routing problem. Transp. Res. Part B Methodol. 94, 61e79. Gendreau, M., Hertz, A., Laporte, G., 1994. A Tabu search heuristic for the vehicle routing. Manag. Sci. 40, 1276e1290. Hoen, K.M.R., Tan, T., Fransoo, J.C., et al., 2010. Effect of Carbon Emission Regulations on Transport Mode Selection in Supply Chains[R]. Retrieved from the Eindhoven University of Technology website: http://alexandria.tue.nl/repository/ books/672727.pdf. Hosseini-Nasab, H., Lotfalian, P., 2017. Green routing for trucking systems with classification of path types. J. Clean. Prod. 146, 228e233. Jabir, E., Panicker, V.V., Sridharan, R., 2017. Design and development of a hybrid ant 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. Leggieri, V., Haouari, M., 2017. Lifted polynomial size formulations for the homogeneous and heterogeneous vehicle routing problems. Eur. J. Oper. Res. 263 (3), 755e767. Li, J., Zhang, J.H., 2014. Study on the effect of carbon emission trading mechanism on logistics distribution routing decisions. Syst. Eng.-Theory Pract. China 34 (7), 1779e1787. Li, X.Y., Tian, P., Aneja, Y.P., 2010. An adaptive memory programming metaheuristic for the heterogeneous fixed fleet vehicle routing problem. Transport. Res. Part E 46, 1111e1127. Li, X.Y., Stephen, C.H.L., Tian, P., 2012. A multistart adaptive memory-based Tabu search algorithm for the heterogeneous fixed fleet open vehicle routing problem. Expert Syst. Appl. 39 (1), 365e374. Lin, C., Choy, K., Ho, G., Chung, S., Lam, H., 2014. Survey of green vehicle routing problem: past and future trends. Expert Syst. Appl. 41, 1118e1138. Niu, Y.Y., Yang, Z.H., Chen, P., Xiao, J.H., 2018. Optimizing the green open vehicle routing problem with time windows by minimizing comprehensive routing cost. J. Clean. Prod. 171, 962e971. Olivera, A., Viera, O., 2007. Adaptive memory programming for the vehicle routing problem with multiple trips. Comput. Oper. Res. 34 (1), 28e47. Soleimani, H., Chaharlang, Y., Ghaderi, H., 2018. Collection and distribution of returned-remanufactured products in a vehicle routing problem with pickup and delivery considering sustainable and green criteria. J. Clean. Prod. 172, 960e970. Sundarakani, B., de Souza, R., Goh, M., et al., 2010. Modeling carbon footprints across the supply chain [J]. Int. J. Prod. Econ. 128 (1), 43e50. Suzuki, Y., 2016. A dual-objective metaheuristic approach to solve practical

pollution routing problem. Int. J. Prod. Econ. 176, 143e153. Taillard, E., 1999. A heuristic column generation method for the heterogeneous fleet VRP. Oper. Res. 33, 1e34. Tajik, N., Moghaddam, R.T., Vahdani, B., Mousavi, S.M., 2014. A robust optimization approach for pollution routing problem with pickup and delivery under uncertainty. J. Manuf. Syst. 33 (2), 277e286. Tarantilis, C., Kiranoudis, C., Vassiliadis, V., 2004. A threshold accepting metaheuristic for the heterogeneous fixed fleet vehicle routing problem. Eur. J. Oper. Res. 152 (1), 148e158. Turkensteen, M., 2017. The accuracy of carbon emission and fuel consumption computations in green vehicle routing. Eur. J. Oper. Res. 262 (2), 647e659. Wu, J., Zhou, Z.X., 2015. A mixed-objective integer DEA model. Ann. Oper. Res. 228 (1), 81e95. Wu, J., Zhou, Z.X., Liang, L., 2010. Measuring the performance of nations at Beijing Summer Olympics using integer-valued DEA model. J. Sports Econ. 11 (5), 549e566. Xiao, Y.Y., Abdullah Konak, A., 2017. A genetic algorithm with exact dynamic programming for the green vehicle routing & scheduling problem. J. Clean. Prod. 167, 1450e1463. Zhang, J., Zhao, Y., Xue, W., Li, J., 2015. Vehicle routing problem with fuel consumption and carbon emission. Int. J. Prod. Econ. 170, 234e242. Zhou, M., Jin, H., Wang, W.S., 2016. A review of vehicle fuel consumption models to evaluate eco-driving and eco-routing. Transport. Res. Transport Environ. 49, 203e218. Zhu, X.Y., Garcia-Diaz, A., Jin, M.Z., Zhang, Y., 2014. Vehicle fuel consumption minimization in routing over-dimensioned and overweight trucks in capacitated transportation networks. J. Clean. Prod. 85, 331e336.