A matheuristic approach for the Pollution-Routing Problem

This paper deals with the Pollution-Routing Problem (PRP), a Vehicle Routing Problem (VRP) with environmental considerations, recently introduced in the literature by [Bektas and Laporte (2011), Transport. Res. B-Meth. 45 (8), 1232-1250]. The objective is to minimize operational and environmental costs while respecting capacity constraints and service time windows. Costs are based on driver wages and fuel consumption, which depends on many factors, such as travel distance and vehicle load. The vehicle speeds are considered as decision variables. They complement routing decisions, impacting the total cost, the travel time between locations, and thus the set of feasible routes. We propose a method which combines a local search-based metaheuristic with an integer programming approach over a set covering formulation and a recursive speed-optimization algorithm. This hybridization enables to integrate more tightly route and speed decisions. Moreover, two other"green"VRP variants, the Fuel Consumption VRP (FCVRP) and the Energy Minimizing VRP (EMVRP), are addressed. The proposed method compares very favorably with previous algorithms from the literature and many new improved solutions are reported.


Introduction
The Vehicle Routing Problem (VRP) and its variants have been the subject of considerable research efforts in the past year, mostly due to the growing number of additional constraints and objectives arising from real-world problems (Vidal et al., 2013a). Given the growing global concern about environmental issues, VRPs have recently started to incorporate "green" aspects such as pollution and alternative fuels, among others (see, e.g., the recent review of Lin et al. 2014 on VRPs with environmental issues).
Environmental costs due to greenhouse gas (GHG) emissions are not usually paid directly by the companies. Nevertheless, some countries are developing emissions trading schemes (e.g., the European Union Emissions Trading System) to make the companies responsible for product between load and distance for each arc. Similar approaches, i.e, those that make use of the vehicle load to minimize the fuel consumption or CO 2 emissions, were presented by Peng e Wang (2009), Scott et al. (2010), Ubeda et al. (2011) and Xiao et al. (2012). The latter introduced the Fuel Consumption Vehicle Routing Problem (FCVRP). Kara et al. (2008) presented a VRP with cumulative costs, which generalizes the EMVRP as well as the Minimum Latency Problem and the m-Traveling Repairman Problem. Kopfer et al. (2013) took the load into account for CO 2 emissions evaluations in the presence of heterogeneous vehicle types.
Minimizing the fuel consumption considering only the load and distance can be insufficient since the travel speed plays a major role. This speed is directly affected by the road congestion. In view of this, Kuo (2010) proposed a model for minimizing the total fuel consumption where the speeds are time-dependent and the load is used to estimate the cost. A Simulated Annealing algorithm was implemented to solve a set of instances from the Solomon benchmark (Solomon, 1987). Later, Kuo e Wang (2011) devised a Tabu Search algorithm for the same problem. Other time-dependent problems with emission minimization can be found in Figliozzi (2011), Saberi e Verbas (2012) and Jabali et al. (2012), and many other references relevant to green logistics are mentioned in the surveys of Dekker et al. (2012), Salimifard et al. (2012), Lin et al. (2014) and Demir et al. (2014b). Bektaş e Laporte (2011) proposed the Pollution-Routing Problem (PRP), which seeks to minimize both operational and environmental costs, taking into account the customers timewindow constraints. The total travel distance, the amount of load carried per distance unit, the vehicle speeds and the duration of the routes are the main cost components. Three different variants were also presented, considering either distance, weighted load and energy minimization. The authors performed an extensive experimental analysisto capture the trade-off between each variation, as well as the effect of the speed and time-window constraints on the distance, energy and costs.
The PRP was addressed with a two-phase heuristic in Demir et al. (2012). In the first phase, the VRPTW is solved by means of an Adaptive Large Neighborhood Search (ALNS), including five insertion operators and twelve removal operators. In a second phase, vehicle speeds are optimized using a recursive algorithm. Computational experiments were carried out for instances with up to 200 customers. Other recent developments consider generalizations of the problem. A bi-objective variant considering fuel and driving time minimization is presented in Demir et al. (2014a), andFranceschetti et al. (2013) consider the time-dependent PRP.
It should be noted that several aspects of the PRP are conflicting. Higher speeds imply routes with shorter durations, but at the same time result in a larger amount of emissions and vice-versa. Hence, to reduce pollution, speed on arcs may be decreased to be closer to the speed which minimizes emissions. Yet with lower speed, the set of feasible VRP routes may be empty or drastically smaller. As a consequence, an optimal solution of the reduced speed VRP can have longer distance and even more emissions in some cases. One main goal of our study was to better integrate route and speed decisions in order to find quickly a suitable balance between these antagonist aspects.

Problem description
The PRP can be defined as follows. Let G = (V, A) be a complete and directed graph with a set V = {0, 1, 2, . . . , n} of vertices and a set A = {(i, j) ∈ V 2 , i = j} of arcs. Vertex 0 represents the depot where a fleet of m identical vehicles with capacity Q is based. Vertices V − {0} correspond to customers, characterized by a non-negative demand q i for a single product, a service time τ i and a specified time-window interval [a i , b i ] for service. We assume that q 0 = 0 and τ 0 = 0. Each arc (i, j) ∈ A represents a travel possibility from node i to j for a distance d ij .
A particularity of the PRP is that the speed v ij on each (i, j) is itself a decision variable, valued between [v min , v max ]. Indeed, each vehicle emits on a certain amount of GHG which depends of weight and speed, among other factors. The PRP aims to find a speed matrix (v) ij and a set of routes R (such that |R| ≤ m) to serve all customers while minimizing environmental and operational costs. Each route σ ∈ R, σ = (σ 1 , . . . , σ |σ| ) starts and ends at the depot, i.e., σ 1 = 0 and σ |σ| = 0, the total demand on each route should not exceed the vehicle capacity, and every customer must be visited within its time window.
Equation (1) defines for a route σ the vehicle load f σ i σ i+1 when traveling on arc (σ i , σ i+1 ), and Equation (2) defines recursively the arrival time t i to customers, knowing that each route starts at time zero, that it is allowed to arrive early at a customer location and wait the start of its time window a i , but that a late arrival is not permitted.
The PRP objective is based on the assumption that CO 2 emissions are approximately proportional to fuel consumption. Using the comprehensive emissions model of Barth et al. (2005), Scora e Barth (2006) and Barth e Boriboonsomsin (2008), Bektaş e Laporte (2011) Parameters w 1 , w 2 , w 3 , w 4 are based on fuel properties, vehicle and network characteristics. Their meaning and value are discussed in Appendix 1, based on Bektaş e Laporte (2011) and Demir et al. (2012). Finally, the overall objective of the PRP, given in Equation (4), considers both the fuel consumption with a cost of ω fc pounds (£) per liter, and the driving costs at a rate of ω fd per unit of time. The The minimum v * f on this function, i.e., the speed value that minimizes fuel costs, is given in Equation (5).
Similarly, for any arc (σ i , σ i+1 ), assuming that there is no waiting time in the route after σ i , the travel cost F fd σ i σ i+1 (v σ i σ i+1 ) including driver wages is given in Equation (6). The speed value that minimizes fuel and driver costs is then expressed in Equation (7). Both values, v * f and v * fd , are independent of the arc under consideration.
The two other problems considered in this work, EMVRP and FCVRP, do not take into account time-window constraints and speed decisions. In practice, these problems can be seen as CVRPs with green-oriented objective functions. The objective of the EMVRP is based on a simplified emission model, which takes solely into account the distance and the load×distance factor (Kara et al., 2007). It is expressed in Equation (8), where ω represents the weight of a vehicle without cargo.
Finally, the FCVRP is based on a linear function of fuel consumption per distance, which considers the vehicle's no-load/curb weight, the load carried, and the fixed cost of the vehicle. The objective is expressed in Equation (9), where h is the vehicle fixed cost, ρ * and ρ 0 are the fuel consumption rate of the vehicle without any load and fully loaded, respectively.
Using the same values of Xiao et al. (2012), with h = 0, ω fc = 1, ρ * = 2, and ρ 0 = 1, the objective can be reformulated as in Equation (10). It is noteworthy that, as in the case of the EMVRP, the FCVRP involves a linear combination of the distance and the load×distance factor.
4 The proposed ILS-SP-SOA matheuristic The proposed algorithm, called ILS-SP-SOA, combines an iterated local search with speed optimization procedures and integer programming optimization over a set partitioning formulation. As a blend of metaheuristic and exact procedures, this category of method is usually called matheuristic (Maniezzo et al., 2009). In ILS-SP-SOA, the initial solution construction and local search, as well as three perturbation procedures aim at minimizing the cost by considering routing decisions without changing the speeds. This sub-problem can be seen as a VRPTW with the objective of Equation (4). We demonstrate in Section 4.1 how local-search moves can be evaluated in O(1) time for this setting. The search is complemented by a recursive speed-optimization procedure, described in Section 4.3, which is applied on each local minimum of the "routing" local search. New speed decisions are included in a dynamic speed matrix, which is in turn used in the subsequent VRPTW sub-problems. Finally, the routes associated to local minimums are stored in a pool, and used by an integer optimization procedure over a set partitioning formulation to generate possible better solutions composed of a different recombination of routes. In this process, each route may be associated to a different speed matrix, and thus this exact procedure adds an additional level of integration between route and speed decisions.
The outline of ILS-SP-SOA is presented in Algorithm 1. The method performs n r restarts (Lines 3-23) of a hybrid procedure combining ILS, SOA and SP. At each restart, the speed matrix v is first initialized with the maximum speed, v ij = v max for all (i, j) ∈ A. A solution S is obtained by applying local search and SOA on an initial solution (Line 8). This initial solution is generated using the modified cheapest insertion heuristic of Penna et al. (2013), considering the PRP objective and the time-window relaxation of Section 4.1. The speed of arcs associated to S is then updated in the matrix v.
Algorithm 1 ILS-SP-SOA(n r , n ils , n sp , n pool , T mip , seed) 1: S best-all ← ∅; f (S best-all ) ← ∞; {S best-all is the overall best solution} 2: P perm ← ∅; i r = 0; { P perm is the pool of permanent routes in the SP} 3: while i r < n r do 4: i r ← i r + 1;

5:
S best ← ∅; f (S best ) ← ∞; {S best is the best solution of the restart phase} i ils ← i ils + 1;

21:
if f (S best ) < f (S best-all ) then

24:
if number of consecutive restarts ≥ n pool then 25: The ILS is then run until n ils successive iterations of local search and perturbation are reached without improvement. At each iteration (Lines 10-18) the local optimal solution is modified by one of the three perturbation mechanisms, selected at random with different probabilities as described in Section 4.2. Time-window constraints are relaxed when performing a perturbation move. This modified solution is possibly improved by applying local search and SOA (Line 12), and the speed matrix is updated (Line 13). A temporary pool of routes is updated during the ILS loop by adding routes associated to local optimal solutions (Line 14). If the number of ILS iterations without improvement is equal to n, the perturbation Change Speeds is applied, as explained later in Subsection 4.2. Finally, the speed matrix v is reinitialized when the number of ILS iterations without improvement is greater than n ils /2. (Lines 17-18).
If the size of the instance is smaller than a given parameter n sp , then the SP method is called only after the last restart phase; otherwise, it is called after every restart. The SP method attempts to create a new solution from routes of a temporary pool of routes P temp derived from local optimums of the local search, which is cleared after each n pool restarts (Lines 24-25), and from a permanent pool of routes P perm which contains the routes associated to the best solutions S best of each restart phase (Line 23). The SP problems are solved using a Mixed Integer Programming (MIP) solver which calls the ILS scheme every time a new incumbent solution is found (Line 20). This collaborative approach, described in Subramanian et al. (2013), not only potentially speeds up the solver runtime, but also may contribute to find better solutions than the best possible combination of routes in the pools. If this happens, the pool P temp is also updated. The MIP solver is run several consecutive times until no improvement is found over S best , and a time limit T mip is imposed to each MIP execution to avoid any unpredictable excessive CPU time. Finally, the algorithm returns the best solution found among all restarts (Line 26).

Local Search with efficient move evaluation
The local search procedure is based on the Randomized Variable Neighborhood Descent (RVND) of Subramanian et al. (2010). It relies on five inter-route neighborhoods, namely: Shift(1,0), Shift(2,0), Swap(1,1), Swap(2,2), 2-opt * and five intra-route neighborhoods, namely: Reinsertion, Or-opt2, Or-opt3, Exchange and 2-opt. The neighborhood structures based on Shift operators move one or more consecutive customers from a route to another one. Those based on Swap consists of interchanging one or more customers from one route with one or more customers from another. Reinsertion and Or-opt follow the idea of Shift, but involving customers of a single route. Exchange is the intra-route version of Swap(1,1). Finally, 2-opt is the classical Traveling Salesman Problem (TSP) intra-route operator, whereas 2-opt * is the interroute version of 2-opt. A detailed description of these neighborhoods, as well as the Auxiliary Data Structures (ADSs) used to enhance the performance the local search can be found in , Penna et al. (2013) andVidal et al. (2013b).
A local search for the VRPTW may have difficulties to generate feasible solutions, possibly compromising the convergence towards good solutions. To circumvent this issue, penalized infeasible solutions w.r.t. time windows are considered, and the objective includes "time-warp" penalties as in (Vidal et al., 2013b). To compute the cost of new routes generated by the local search in amortized constant time, other ADSs based on subsequences concatenation were implemented. For any subsequence σ the algorithm stores and maintains: • minimum duration T (σ), • minimum time-warp use T W (σ), • earliest E(σ) and latest visit L(σ) to the first vertex allowing a schedule with minimum duration and minimum time-warp use, For a subsequenceσ involving a single customer, i, these ADSs are computed as follows: If the first node in the sequence is a depot, e.g., ifσ = 0, then E(σ) = L(σ) = 0. This prevents a delayed departure. The following equations enable then to derive ADSs for larger subsequences obtained by concatenation ⊕: The penalized cost of a route σ can be derived from these structures as shown in Equation (23) where ω tw is the penalty for one unit of time warp.

Perturbation Mechanisms
The three mechanisms applied during the perturbation phase are described as follows.
• Shift to End -Shifts one customer from one route to the end of other route, randomly with uniform distribution.
• Change Speeds -This perturbation modifies the speed matrix by changing the speeds associated to arcs of one random route in the current best solution S best with one random speed from the set {v * f , v * fd , v max }, with uniform distribution.
During the perturbation phase, Shift to End and Merge routes are randomly selected with different probabilities, 90% and 10%, respectively. Any perturbation leading to an infeasible solution w.r.t. capacity constraints is undone and a new perturbation is attempted. Finally, the additional Change Speeds perturbation is applied after n consecutive ILS iterations without improvement, e.g., when i ils = n.

Speed Optimization Algorithm
Speeds have a direct impact on customers' arrival time, i.e., in meeting their time windows. When routes are fixed, the PRP leads to a speed optimization problem which consists of finding the optimal speeds for each arc while respecting customers' time windows. This problem is solved by a recursive algorithm (Algorithm 2). This method is an adaptation of the RSA procedure of Norstad et al. (2011) and Hvattum et al. (2013), also considering drivers wages, possible waiting times, and a non-fixed arrival time at the last customer.
Algorithm 2 is applied on the complete route, by setting s = 1 and e = |σ|. Let t σ i be the time to start the service of customer σ i . Early and late arrivals are not considered in the course of the RSA procedure such that t σ |σ| is also the arrival time at customer σ i . The algorithm first computes the arrival time t σ |σ| at the last customer when traveling at speed v * fd (minimizing fuel consumption plus driver costs). If this time t σ |σ| is greater or lower than the time windows bounds, it is updated to the nearest bound (Line 6). Then, the necessary speed v ref to arrive at time t σe is computed (Line 7), and the algorithm finds the customer σ p with greatest timewindow violation when using this speed (Lines 10-13). If no violation is found, the solution is returned. Otherwise, the arrival time t σp is updated to the nearest time window bound (Line 15) and SOA is called recursively on two subproblems: from s to p (Line 16), and from p to e (Line 17). Once times t σ i , ∀i = 1, . . . , |σ|, are computed, the associated speeds are revised in such a way that any speed below the optimal speed v * f which minimizes fuel consumption is replaced by a speed v * f and a waiting time (Lines 18-22). The final arrival dates t i at route nodes are obtained (Line 22).
Algorithm 2 is illustrated in Figure 1 on a problem with eight stops. The horizontal axis represents the time, while the vertical axis indicates sequence of customers from bottom to top. The brackets correspond to the time windows and the dots denote the customers' arrival times. Round black dots are related to feasible arrival times; diamond gray dots indicate an early or late arrival time; and square dots indicate the new times after correction. For simplicity, the service times are assumed to be zero. Figure 1a illustrates the arrival times at each customer when the vehicle travels with the optimal speed (v * fd ) starting from σ 1 = 0 at time t σ 1 = a σ 1 . There are two violations and σ 7 Algorithm 2 Speed Optimization Algorithm -SOA if violation > maxV iolation then is the customer with the greatest violation (p = 7). The arrival time at σ 7 is thus adjusted to t σ 7 = b σ 7 (Figure 1b) and two subproblems, P1 and P2, are solved, where P1 considers the subsequence from σ 1 until σ 7 , whereas P2 considers the subsequence from σ 7 until σ |σ| . In P1, the v ref is computed to arrive at σ 7 at time t σ 7 . A new violation is observed at customer σ 4 (p = 4). This violation is corrected ( Figure 1c) and two new subproblems, from σ 1 to σ 4 (P1.1), and from σ 4 to σ 7 (P1.2) are solved. When solving P1.1 and P1.2, no violations are identified (Figure 1c-1d), and thus the recursion is not applied further. This leads to a feasible solution to P1. Subproblem P2 is solved using almost the same rationale of P1 (see Figure1e), except that the customers' arrival times are computed using the optimal speed v * fd (see Line 6 of Algorithm 2). Again, no violation is found, thus completing the whole solution for the problem. At the end, speeds are revised to insert waiting time whenever a speed v < v * f is encountered (Lines 18-22) of Algorithm 2).

Computational Results
In this section we report the computational experiments for the three variants considered in this work, but first we describe the benchmark instances used for evaluating the performance of our algorithm.

Benchmark Instances
The first set of instances considered in this work (Set A) is the one from the PRPLIB, suggested by Demir et al. (2012) and available at http://www.apollo.management.soton.ac. uk/prplib.htm. This set consists of nine different groups, ranging from 10 to 200 customers, each one containing 20 instances. For these instances, the same objective parameter values as However, these instances have a large time windows width, such that it is possible to visit many customers within their respective time windows when traveling at optimal speed as further discussed in Section 5.2. In view of this, we created two additional sets of instances with tighter time windows by modifiying those of the PRPLIB. The time horizon of these new sets is 32400. The time-window width of each customer in Set B is randomly selected between the interval 2000 and 5000 with uniform probability, whereas in Set C it is randomly selected between 2000 and 15000. Once the time windows width are defined, the time window lower bound is randomly selected with uniform probability in the time interval I i which allows to feasibly reach the customer and return to the depot before the end of the time horizon. For a customer i, let W i be the chosen width, and let δ min 0i = d 0i /v max and δ min i0 = d i0 /v max be the minimum driving time from and to the depot, respectively. Then, For the FCVRP and the EMVRP, the well-known benchmark instances of Christofides et al. (1979) and Golden et al. (1998) are used. As in Kara et al. (2007), we use ω = 0.15Q for the EMVRP.

Method performance on the PRP
The proposed matheuristic algorithm was coded in C++ and executed on an Intel Core i7 3.40 GHz with 16 GB of RAM, running under Linux Mint 13. CPLEX 12.4 was used to solve the SP problems. Only a single thread was used, and the algorithm was run 10 times for each instance with different random seeds.
The following parameter values were adopted: n r = 20, n ils = n + 5m, n sp = 150, n pool = 2, T mip = 360 s, and ω tw = 10 8 . The first four parameter values directly follow from  and Subramanian et al. (2013), whereas the latter is a large number to prevent infeasible solutions due to late arrivals. A higher value of T mip than Subramanian et al. (2013) is used, as more time seems to be needed by the SP solver to produce improved solutions.
To investigate the interaction between the local search components, the SP solver, and the speed optimization procedure, we performed further experiments with a static version called ILS-SP-SOA-Stat. For a fair comparison, the CPU time limit for each instance has been set to the average CPU time of the dynamic method for the corresponding instance. In ILS-SP-SOA-Stat, v does not change during an iteration of the algorithm, i.e., the "UpdateSpeedMatrix" function and the Change Speeds perturbation is not used in Algorithm (1). Hence, the task of intelligently recombining routes with different speed distributions is exclusively operated by the SP solver.
The performance of each method is reported in Table 1. Each line corresponds to averaged results on a set of 20 instances. The columns report for each method and instance set the average Gap(%) and CPU time. The Gap(%) for each instance is computed as 100(Z − Z bks )/Z bks , where Z is the objective value of the solution and Z bks is the value of the Best Known Solution (BKS) ever found. Detailed results on each problem instance are provided in B.
It appears that both versions of the proposed algorithm largely outperform the ALNS heuristic of Demir et al. (2012) in terms of solution quality. The new methods also appear to be significantly faster than ALNS, with a speed-up ranging from ×2 to ×50. Note that ALNS experiments have been conducted on a different computing environment and a 3.00GHz CPU. The speed difference between computers of the same generation remains moderate. It can be observed that instances of Sets B and C require more CPU time for convergence, and solution improvements are found more continuously as the search progresses. This is due to tighter time windows. The increased interplay between route and speed optimization leads to a longer sequence of solution improvements, and thus a slower convergence. This aspect is further developed in the next section.
Comparing the results obtained with the two versions of ILS-SP-SOA, the average solutions obtained with the dynamic version are better than those from the static version in all sets of instances but 200-B and 200-C. The dynamic version examines a wider diversity of speed choices. For smaller instances -here up to 100 customers -this allows to converge towards better solutions within the allowed termination criteria. For larger instances with tighter time windows, the benefits of enhanced speed optimization and diversity seems to be counterbalanced by a slower convergence, and the static method, more focused on intensification, may perform marginally better. A longer termination criterion may thus benefit more to the dynamic version that to the static version.

Some insights on instance difficulty for the PRP
To better understand what makes a difficult PRP instance, and how these instances can be better addressed, we further analyze instance characteristics and study how often the optimal speeds are used in the best known solutions. As shown in the following, higher occurrences of "exotic" speeds in the BKS, i.e., speeds that are neither v * f nor v * fd , seem to translate into increased problem difficulty. Figure 2 shows the percentage of the number of arcs in which the vehicle travels at a certain speed in the BKS of each group of instances. Three speed categories are discerned: v * f , v * fd , and other speeds v other . It can be seen that v * fd appears to be used more often than the other speeds in Set A, regardless of the size of the instance. As expected, v * f tends to be used much more often in Set B, which is the one with the tightest time windows. Finally, in Set C, where the time windows are larger than those in Set B, v * fd is the most common in the 10-customer instances, but there is an equilibrium in the 50-, 100-and 200-customer instances between v * f and v * fd . Let %Dist be the percentage of the total distance in which the vehicles are traveling with other speeds than v * f and v * fd in the BKS of a given instance. Figure    A similar analysis is performed in Figure 4 considering the CPU time of ILS-SP-SOA-Dyn, and thus the convergence speed of the method. Overall, the CPU time needed to converge towards a final solution tends to increase for problems with higher %Dist. This suggests that the instance difficultly is correlated with the need of a variety of speeds in the solutions.  This analysis shows that, even if instances with large time-windows are reasonably well solved by methods with separate phases such as the ALNS of Demir et al. (2012) or ILS-SP-SOA-Stat, there is still a notable research effort to be done for more intricate problems, when a larger diversity of speeds must be used to fulfill the time constraints. Because of the variety of antagonist interplays between, speed, distance, and feasibility, more integrated methods that consider jointly speed and route optimization within the neighborhoods are likely to be necessary. ILS-SP-SOA-Dyn is already a first step in this direction, and has demonstrated increased performance on some of the most difficult instances.

Results on other problems
Computational experiments have been conducted on two additional problems with environmental considerations, the FCVRP and EMVRP, in order to assess the performance of ILS-SP-SOA on a wider variety of settings. As previously, a single thread was used and 10 executions were performed for each instance. For a fair comparison with previous methods developed for the FCVRP, we reduced the termination criterion by setting n r = 4, n ils = n/5 + 5m and T mip = 60 s. For the EMVRP we only changed the MIP solver time limit to T mip = 60 s.
For the FCVRP, the results obtained by ILS-SP-SOA were compared with those of Xiao et al. (2012), as can be observed in Table 2. The values highlighted only in boldface indicate that the BKS found by our algorithm equaled the best one reported by Xiao et al. (2012), whereas those in boldface and underline indicate that the BKS found by ILS-SP-SOA improved the one presented by the same authors. To our knowledge, no results on the instances of Christofides et al. (1979) were available for the EMVRP. Hence, Table 3 reports the gap between the average solutions and the BKS found by ILS-SP-SOA.
These results show that ILS-SP-SOA consistently produces solutions of high quality on these two problems. We observe that 22 new BKS solutions were found for the FCVRP, and the average gaps between the average solutions of ILS-SP-SOA and the BKSs is 0.16% for the set of instances of Christofides et al. (1979) and 0.31% for the one of Golden et al. (1998)  do not have a dramatic impact in the solution cost when compared to the associated instances that do not consider these characteristics. Overall, the gaps between the average solutions and the BKS are small, more precisely 0.09% and 0.08%, on the instances of Christofides et al. (1979) and Golden et al. (1998), respectively. When observing such a small variance for VRP problems we can reasonably conjecture that the problem is fairly-well solved.

Concluding remarks
This paper dealt with the Pollution-Routing Problem (PRP), a "green"-oriented variant of the Vehicle Routing Problem (VRP). In order to solve it, we proposed a matheuristic approach, called ILS-SOA-SP, that effectively integrates Iterated Local Search (ILS) with a Set Partitioning (SP) procedure and a Speed Optimization Algorithm (SOA). This approach was also used to solve two other environmental based VRPs, namely the Fuel Consumption Vehicle Routing Problem (FCVRP) and the Energy Minimizing Vehicle Routing Problem (EMVRP). The results of extensive computational experiments demonstrate that ILS-SP-SOA is capable of generating high-quality solutions in a very consistent manner, outperforming previous methods from the literature.
Two new sets of PRP instances were introduced, and our computational experiments show that they are more challenging than the previous one available in the literature. We also showed that the level of difficulty of an instance tends to increase when the vehicles are likely to travel at speeds other than the optimal ones. In other words, good (or even feasible) routing solutions are not easily obtained without using non-optimal speeds, and choosing the correct ones can be a challenging task. The proposed method has the advantage to perform joint route and speed optimization within several local search and integer programming components, and thus performs much better on difficult instances than previously available algorithms.
As a promising line of research, we believe that further efforts are still necessary to incorporate speed decisions more tightly within the local search. This can be done, for example, by developing larger neighborhood structures that are capable of exploring promising regions of the solution space by considering different speeds during the local search. The exploration of these combined neighborhoods is likely to be highly time-consuming, and thus very clever move selection and move evaluation methods may need to be developed.

A Fuel consumption model
The fuel consumption model used in the PRP formulation, given by Eq. (25), is based on the comprehensive emissions model described by Barth et al. (2005), Scora e Barth (2006) and Barth e Boriboonsomsin (2008). It is given by where λ = ξ/κψ and γ = 1/1000η tf η are constants related to the fuel properties, β = 0, 5C d ρA and w are constants associated with the vehicle characteristics and α ij is constant that depends on the road characteristics and on the vehicle acceleration. More precisely, α ij = τ ij + g sin θ ij + gC r cos θ ij , where τ ij is the acceleration and θ ij is the road angle inclination. The meaning and value of all parameters is given in Table 4. Assuming that acceleration and road inclination are null, the fuel consumption can alternatively be written as where w 1 = λkN V , w 2 = λwγgC r , w 3 = λγgC r and w 4 = λβγ. Conversion factor (g/s para l/s) 737 vmin Lower speed limit (m/s) 5.5 (or 20km/h) vmax Upper speed limit (m/s) 25 (or 90km/h) Source:  B Results for the PRP In the tables presented hereafter, Instance denotes the name of the instance, Cost indicates the value of the objective function, |R| is the number of vehicles, CPU T(s) represents the computing time in seconds, Gap(%) is the gap between the solution cost and CPLEX for the 10-customer instances or the best known solution obtained among all experiments for the other instances, given in the column BKS. Note that the approximate optimal solutions for the 10customer instances were computed with the MIP formulation of Demir et al. (2012). It is impractical to consider all possible speed values in the model, and thus a speed-discretization was operated with 500 levels ranging from v * f to v * fd . This leads in general to very high-quality upper bounds.