Fuel Emissions Optimization in Vehicle Routing Problems with Time-Varying Speeds

The problem considered in this paper is to produce routes and schedules for a fleet of delivery vehicles that minimize the fuel emissions in a road network where speeds depend on time. In the model, the route for each vehicle must be determined, and also the speeds of the vehicles along each road in their paths are treated as decision variables. The vehicle routes are limited by the capacities of the vehicles and time constraints on the total length of each route. The objective is to minimize the total emissions in terms of the amount of Greenhouse Gas (GHG) produced, measured by the equivalent weight of CO2 (CO2e). A column generation based tabu search algorithm is adapted and presented to solve the problem. The method is tested with real traffic data from a London road network. The results are analyzed to show the potential saving from the speed adjustment process. The analysis shows that most of the fuel emissions reduction is able to be attained in practice by ordering the customers to be visited on the route using a distance-based criterion, determining a suitable path between customers for each vehicle and travelling as fast as is allowed by the traffic conditions up to a preferred speed.


Introduction
Technical developments and the growth in road traffic pose new challenges for research in vehicle routing and scheduling for freight transport. Remote vehicle tracking techniques enable the road traffic data for different times of day and different days of the week to be collected, so as to provide detailed information on transit times for different roads by time of day and day of week. This provides an opportunity to plan vehicle routes and schedules taking time-varying speeds into account. In addition, the growth in road traffic and the use of road freight transport also bring problems of environmental pollution. Concerns about the environmental impact of transport activities have led to new vehicle routing models where the objective is to minimize the harmful effects of transportation on the environment.
An increasing number of papers are being published where fuel emissions are explicitly modeled. However many of these simplify the model by assuming that paths between customers are fixed or that the speeds of the vehicles are time-independent. In the model described in this paper, the speed of the traffic on the underlying road network is time dependent. In addition, the path used by a vehicle between a pair of customers and the speeds on the road segments are decision variables. This paper will describe a column generation based tabu search algorithm, which can work together with a solution method for single paths, in order to minimize fuel emissions for Vehicle Routing Problems (VRPs) with time-varying speeds.
The algorithm is then used for modeling a distribution operation using real traffic data from a road network located in London. The aim of these experiments is to discover how much reduction in CO2e can be obtained by using the algorithm described in this paper, compared with other approaches that are faster to compute. Experiments are also carried out to determine the effect of allowing more waiting time at customers. The paper is organized as follows. Section 2 contains a review of relevant literature. The problem is described in Section 3, which is followed by a set-partitioning model for the VRP. Section 4 introduces the framework of the column generation based tabu search algorithm, which is used to find a prospective sequence for a set of customers, and goes on to discuss details of the algorithm. The computational experiments and their results are then presented in Section 5. Finally, conclusions are drawn in Section 6, and the main findings are highlighted.

Literature Review
In recent years, there has been increasing interest in estimating the environmental effects of vehicle routing policies. A survey of recent work in this area can be found in [9].
Various models have been proposed for estimating the fuel used by vehicles when travelling on roads. Examples include one published by the European Commission in the MEET report described by Hickman et al. [17] and the Comprehensive Modal Emissions Model (CMEM) described by Scora and Barth [27]. The CO 2 emissions are normally calculated as being proportional to the fuel used. The fuel consumption and hence the emissions may relate to factors such as the vehicle type, weight and speed. Demir et al. [3] provide a comparison of a number of such models.
Recent research on minimizing emissions in vehicle routing models can be divided into two main categories: the first is the set of models where time-independence is assumed and the second set includes models where the road conditions are subject to traffic congestion and so the time needed to travel along a road segment depends on the time of day.
Among the time-independent models, Palmer [23] developed a model where vehicle speeds are inputs to the model and the approach is tested on a case study of home deliveries for grocery stores in the UK. He found an average saving of 4.8% in CO 2 emissions was possible compared to using routes that minimize time, but at the expense of a 3.8% increase in the time required. His model does not take vehicle loads explicitly into account, but Suzuki [28] uses a model where load is taken into account and finds that delivering relatively heavy items early in a tour can be worthwhile in reducing the fuel consumption. Several case studies have been reported using time-independent approaches with the objective of minimizing fuel consumption and hence emissions. For example, Ubeda et al. [29] consider emission factors in planning routes for a food delivery operation. They show savings of around 25.5% in CO 2 emissions, but this is mainly due to reducing the number of routes needed compared to the original plan.
Other time-independent models allow the speeds of vehicles to be decision variables. The approach adopted by Bektaş and Laporte [1] in their Pollution-Routing Problem (PRP) uses a CMEM-based model and considers both load and speed in estimating a cost function to be minimized. They propose a non-linear mixed integer mathematical programming formulation and show how it can be linearized. The formulation can only solve small PRP instances, but Demir et al. [4] provide an adaptive large neighbourhood search algorithm for much larger PRP instances. Van Woensel et al. [30] develop a model showing how queuing theory can be used to describe traffic flows and calculate emissions using the model described in the MEET report.
In the set of time-dependent models, Eglese et al. [10] make use of traffic speed information collected at different times on sections of a road network to create a Road Timetable showing the quickest times between origins and destinations starting at different times of the day. In Maden et al. [21] the Road Timetable is used with a tabusearch called LANTIME to minimize the total time required for a distribution operation.
Vehicles are assumed to travel at the speed which minimizes their emissions per unit distance unless the congestion indicates that this is not possible, when the vehicles travel at the average speed of the traffic recorded for that road segment at that time. Results from a case study based on the distribution plans for an electrical goods wholesaler in the UK show that CO 2 emissions can be reduced by around 7% with this approach. This is because routes with high congestion and hence, enforced low speeds and high emissions, are avoided. Figliozzi [13] also takes into account congestion in minimizing emissions using a model based on the MEET report. An integer programming formulation is presented and a solution algorithm is described which is tested on modified Solomon benchmark problems. In contrast to Maden et al. [21], the model allows vehicles to travel faster than their optimum speed that minimizes emissions if the traffic conditions and Vehicle speeds may also be used as decision variables in time-dependent models. Jabali et al. [19] use a similar model to Figliozzi [14] but with speed as an additional decision variable, though without the use of time windows. Their model is based on a complete network where the nodes represent the depot and customers, while the maximum speeds on the arcs linking the nodes are subject to similar profiles. They describe a tabu search heuristic for solving the problem and test it on standard benchmark instances. The results suggest that a reduction of about 11.4% in CO 2 emissions can be achieved, but with a 17.1% increase in travel times. Franceschetti et al. [15] follow a similar approach which also takes costs into account in a similar way to Bektaş and Laporte [1]. A mathematical formulation is produced and provides insights on when it is profitable to wait at customers.
The model presented in this paper is in the last category of models which take into account time-dependent conditions and where vehicle speeds are decision variables. It is designed for use on a road network where information is available on the speed of traffic on individual road segments at different times of the day. The solution provides the path to follow between customers and the speeds to be applied on each road segment. It thus provides a more detailed model than the one used by Jabali et al. [19]; the path used between a pair of customers may change depending on the time of travel in our model. Also, it allows time window constraints for serving customers which are not included in Jabali et al. [19].
If it is assumed that the path used between customers is fixed, then some other recent research on speed optimization is relevant. Fagerholt et al. [12] present the Speed Optimization Problem (SOP) in the context of shipping, provide models to formulate the problem and a solution algorithm. Norstad et al. [22] provide a recursive smoothing algorithm for the SOP that runs fast and has been shown to be optimal by Hvattum et al. [18].
It is often the case that reductions can be made in the emissions resulting from a distribution operation, but at the expense of more time or cost. There are methods explicitly aimed at modelling this issue through a multi-objective approach. One example is provided by Jemai et al. [20] where an evolutionary algorithm is used to solve a biobjective VRP where one objective minimizes total distance, while the other minimizes CO 2 emissions. Demir et al. [5] consider the bi-objective PRP where fuel consumption and driving time are the two relevant objectives.
Finally, there is an emerging strand of research considering vehicle routing problems for alternatively powered vehicles that are designed to be more environmentally friendly.
Such vehicles may have a more limited range before requiring refuelling and there may be a limited availability of refuelling points. An example is given by Erdoğan and Miller-Hooks [11] in which they define a "Green Vehicle Routing Problem" where there are additional constraints on how far the vehicles may travel without refuelling and the refuelling stations are located at specific places. They formulate a mixed integer program to minimise the total distance and develop heuristics for its solution. Tests are carried out based on the location of stations supplying biodiesel fuel in a part of the USA. corresponding emissions, f p , and are determined by a new heuristic approach referred to as NHA and described in [26]. The quantity is a binary variable used to denote whether customer is visited by a vehicle travelling along route p, where if customer is included in route p, and , otherwise. Let be a decision variable to count the number of times route p is used. In a feasible solution, should be 0 or 1 for all . The formulation is shown as follows.

Column Generation Based Tabu Search Algorithm
The proposed solution method uses a column generation algorithm that takes advantage of the power of the branch-and-price technique to solve a set partitioning problem. It is based on a branch-and-price-based large neighborhood search algorithm for the Vehicle Routing Problem with Time Windows (VRPTW) proposed by Prescott-Gagnon et al. [24].
The current solution is destroyed in the destruction step by selecting one of four operators randomly. This leaves a set of partial routes and isolated customers. The large neighbourhood then contains all feasible solutions that are compatible with the partial routes. A heuristic column generation algorithm is used to reconstruct the solution, where tabu search is used to generate columns of negative reduced cost.
The formulation of the restricted master problem (RMP) called within the method is described as follows.
The set contains feasible routes that include the fixed parts of the solution that have not been removed in the destruction step. We wish to find whether the objective value of the RMP could be further improved by adding any route that has not been included in . Any route that improves the solution will have a negative reduced cost when added to the RMP. The reduced cost for any route can be calculated using the values of the dual prices associated with the constraints in (5). If no route exists that was not in with a negative reduced cost, then the solution is optimal. Therefore, the sub-problem of the column generation process is to find the column, which satisfies the time window and capacity constraints, with minimum reduced cost. The sub-problem can be formulated as: , where is the dual price associated with constraint i in (5).
A tabu search algorithm is used to solve the sub-problem heuristically, and the goal is to produce feasible routes with negative reduced costs that have not yet been included in set . The algorithm framework is briefly introduced in the following steps.
Step 1: Generate an initial solution with Clarke and Wright Savings Algorithm (parallel version). Go to Step 2.
Step 2: Apply a destroy operator to determine the neighbourhood. Go to Step 3.
Step 3: Solve the RMP, i.e. a set-partitioning problem formulated as an LP relaxation based on the current set of columns (routes). If the stopping rule for the column generation process is met, go to Step 6; otherwise, go to Step 4.
Step 4: Apply tabu search to generate new columns (routes). Go to Step 5.
Step 5: Find the fuel emissions for the new columns (routes) using NHA. Go to Step 3.
Step 6: If the solution is not integer, apply a branching strategy to get an integer solution.
While the number of iterations for the algorithm is less than a specified maximum number, go to Step 2; otherwise, stop.
Details of the algorithm corresponding to these steps are discussed in the following sections. Additional details can be found in the PhD thesis [25]. if the weighted average speed is greater than the optimal speed. The optimal speed is the speed of the vehicle that minimizes the fuel emissions per kilometer. Once a speed for each arc is determined, the emissions for each arc can be calculated and treated as a cost for each arc. The optimal path between customers and between depot and customer in terms of fuel emissions is calculated using Dijkstra's algorithm [8] and the corresponding travel cost (fuel emissions), travel distance and travel time are used as the initial values of the elements in , and .

The initial solution
The initial customer assignment is computed by the Clarke and Wright Savings Algorithm (parallel version) [2]. Each customer is assigned to an individual route at the beginning. If the arrival time at customer , which is approximately computed by is before the early time window , unlimited waiting is allowed. Thus, the starting service time at customer denoted as is, , and the corresponding waiting time denoted as is, .
Then routes are combined using the estimated travel costs based on fuel emissions.
Capacity and time window constraints are considered when deciding if two routes can be merged.

The Destroy Operators
There are four destroy operators used to modify the current solution in order to diversify the search. The destroy operators are analogous to the ones presented by Prescott-Gagnon et al. [24] and more details are given in [25]. With each destroy operator, only part of a complete route will be removed.
The large neighbourhood search is a process where each iteration starts with the application of a destroy operator to the current solution. Only one destroy operator is applied at each iteration through a roulette-wheel selection procedure. Each of the destroy operators is assigned a value . The destroy operator is chosen randomly at the beginning of each iteration, where each operator i has probability of to be selected. Following the method proposed by Prescott-Gagnon et al. [24], the initial value of is set to 5. Whenever operator i is called and its use improves the current solution, the value of is increased by 1.

Solving the RMP
The column generation process has already been outlined at the beginning of Section 3. In Step 3 of the algorithm framework, the stopping rule is that either an improved integer solution has been found or the value of the objective in (4) has not been improved for a certain number of iterations.
A solution pool is kept to store the routes associated with the basic variables obtained by solving the last RMP at each iteration of the large neighbourhood search. After a destroy operator is selected and applied, the remaining customers in the route define the fixed part.
Then, the routes in the solution pool with the same structure as the fixed part will be added to set ′ . As a result, the reconstruction process will start with a set of prospective columns, some of which may not be obtained by applying tabu search.
Moreover, the fuel emissions for these routes do not need to be recalculated, which saves computational time. However, extra time is required to scan the structure of the routes in the solution pool, so the size of the solution pool is limited in order to control the extra time needed to manage it. Columns are removed by following the 'first in first out' rule once the solution pool is full. The solution pool can also be thought of as a long-term memory.

The tabu search algorithm
The tabu search algorithm is used to solve the sub-problem heuristically, and the goal is to produce feasible routes with negative reduced costs that have not yet been included in set ′. The costs used in the tabu search are based on the approximate costs, but they are adjusted by the dual prices, , corresponding to the constraints in the RMP as follows.
Routes associated with the basic variables after solving the RMP are used as initial solutions, and the tabu search approach is applied individually for each of them. Since the reduced costs of the basic variables are zero, they should be good initial solutions when searching for routes with negative reduced costs. The total number of iterations for the tabu search is at most . As a result, the number of iterations allowed for each column is set to .

The tabu search neighbourhood
Two simple operators are used to reconstruct new routes, which are removing and inserting a customer. Every time a customer node is removed/inserted, the reverse move is tabu for the next iterations. Hence, at one iteration, a move is defined as removing a (non-tabu) customer and inserting another (non-tabu) customer at a possible insertion place. All possible removals and all possible insertion places are tested to find the best one. Only feasible moves are allowed in the search process, so the capacity and time window constraints have to be checked for each customer insertion and only the time window constraints for customer removal.

The challenge of neighbourhood move evaluation
One of the advantages of applying a heuristic method to solve a VRP in a static network is that only the sub-routes changed by a neighbourhood move have to be re-evaluated.
However, in a time-varying road network, the evaluation of a neighbourhood move is not so straightforward. The fuel emissions from any customer to customer may change according to different departure times from customer , so the re-evaluation process is no longer simply an addition and subtraction operation of a few static values relating to the links between customers that have been changed. The time dimension has to be taken into account which leads to the need to re-evaluate significant parts of the affected routes.
Harwood et al. [16] has examined different ways to estimate the cost of a neighbourhood move within a single tour with time-varying traversal speeds, when the objective is to minimize the total time. The tour is divided into three parts according to the nodes being moved: the pre-change part is the tour from the depot until the first node to be changed; the post-change part is the tour from the last node to be changed back to the depot; and the remaining section is called the changed part. To determine whether a move leads to an improvement, the pre-change part of the tour does not need to be recalculated and provided the First-In-First-Out (FIFO) property holds, the post-change part of the tour does not need to be recalculated either. Furthermore, the necessary and sufficient condition for an overall improvement to be achieved is that the tour from the depot to just before the post-change part should be improved. Therefore, only the changed part of the routes has to be revaluated to find out whether a neighbourhood move will improve the cost, if the FIFO property is maintained in a time-varying network. However the postchange part of the route still has to be recalculated in order to obtain the overall improvement, once a neighbourhood move has been identified as leading to an improvement.
When the objective is to minimize emissions or costs that depend on the time of travel, then these results do not apply. Suppose the optimal solution has been found from an origin node to node which passes through . That solution may not contain the optimal solution for a path from to . This is because it may be better to travel faster/slower than the optimal speed with lower fuel efficiency from to , so as to arrive at earlier/later to avoid the congestion from to .
Consequently, any neighbourhood change requires the whole route to be re-evaluated to Under the following two circumstances, NHA is applied to evaluate fuel emissions for the current complete route, so that more accurate values can be used:  The estimated reduced cost of a new column is better than the best solution to the sub-problem associated with the RMP (7) so far. In this case, NHA is used to calculate the actual fuel emissions of the new column, in order to obtain the actual reduced cost. If the reduced cost is smaller than the best solution so far, the aspiration criterion is met, and the tabu restriction is overridden.
 After iterations of tabu search have been applied to a column, the fuel emissions of the column will be calculated by NHA. If the reduced cost based on the re-estimated fuel emission is negative, this column is added into set ′.
If there is an improved integer solution obtained by solving the RMP, the reconstruction process will stop, and a destroy operator will be invoked to determine a new large neighbourhood of the integer solution. Otherwise, the column generation process will stop if the value of (4) cannot be improved in the last iterations of solving the RMP, where represents the maximum number of column generation iterations and is a parameter to be set.

Branching Strategy
When the solution obtained by RMP in the last iteration is not integer, a heuristic branching strategy is applied to derive an integer solution. The branching strategy is to simply fix the decision variable with the largest fractional value at 1, and solve the linear problem again. This process is repeated until the solution is integer or no feasible solution can be found.
If any new integer solution can be obtained by applying the branching strategy, no matter whether it is better or worse than the initial solution at the beginning of the current iteration of large neighbourhood search, the solution will be used for the next iteration.
This can help to diversify the search. If there is no feasible integer solution, the initial solution at the beginning of the current iteration of large neighbourhood search is used again in the next iteration.

Experiments with Real Traffic Data
The VRP algorithm is tested in this section with real traffic data for a road network located in London.

Data Description
The

Fuel Curve for HGVs
gradients require less fuel. Also no allowance has been made for the changing weight of the load carried. The fuel used in litres is converted to the equivalent greenhouse gas effect of CO 2 in kg by multiplying by a factor of 3.1787. This was the factor proposed by the UK Government in 2010 in its guidance for reporting emissions [6].

Results Analysis
The approach described in this paper was applied to the instances described and a summary of the results is provided in Appendix B. This shows the sequence of customers on each route, the total CO 2 emissions, the distance travelled, the time required and the total waiting time for each instance. There are two sets of results: one set for instances without customer time windows and one set for instances with time windows.
In the computational experiments, the parameter values used are and the length of the tabu list . Details concerning the parameter values can be found in Qian [25]. All the results are recorded as the better of two runs. All computations were carried out on the Lancaster high performance computer cluster with quad core Intel In the following subsections, various aspects of the problem and solution method are investigated. The first subsection investigates the impact of the speed adjustment and path selection process, waiting time and starting time based on the tests on the five customer sets without time windows. The speed adjustment and path selection process refers to the path chosen and speeds adopted on the road network in travelling between specific customers. Then, the original method is compared against a simpler method, which generates customer sequences based on static information. Finally, the effect of the time window constraint is examined.

Speed Adjustment & Path Selection
NHA, as described in [26], finds the best path between two customers and the best speeds  The path selection and speed adjustment processes reveal some contribution to the CO2e emissions reduction. On the average, about 10 kg CO2e emissions can be saved by applying both the path selection and speed adjustment processes, equivalent to a reduction of about 2-3%. In particular, the reduction for B_0 is up to 15 kg, which is 4% Scenario Instance less than 'Without PS & SA'. Compared to the speed adjustment process, the path selection process plays a more significant role in reducing fuel emissions. Around 4 kg extra CO2e is emitted on average by choosing a promising path and travelling with the fastest speeds allowed by the traffic conditions up to an optimal speed; while in the worst case B_0, another 12 kg emissions is produced by following the fastest paths at the fastest allowed speeds.
One reason that the performance of the speed adjustment process does not have a greater effect could be the shape of the fuel curves used in the experiments. As illustrated in Figure 1, the maximum valid speed is 90km/h, whose fuel efficiency is close to that of the optimal speed 65 km/h; and the fuel efficiency reduces quickly with the decrease of the speed. Therefore, it would not be wise to slow the speed any more than necessary. If the fuel emissions curve is very sensitive to the speeds, or increasing the speed makes the fuel curve steep, then the speed adjustment would be worth applying; otherwise, just carefully selecting a path, then travelling as fast as traffic conditions allow up to an optimal speed still provides a good solution.

Starting Time & Waiting Time
All experiments so far are based on the condition that all the deliveries should start at 7 am, and waiting time at each customer node as well as the depot cannot exceed 5 minutes.
Note that 7 am is the start of a peak period, when the road traffic is getting busy. Vehicles may be caught in traffic congestion, which results in extra fuel emissions. The fuel emissions may be improved by allowing longer waiting, or starting after the rush hour.
Two scenarios are tested: one is starting at 7 am, but allowing waiting up to 4 hours; and the other is starting at 10 am, and with maximum waiting time being 5 minutes.  Instance and the traffic flows in the city centre keep busy during the day time. It may also be because more travel is moved towards the evening peak when traffic is more congested.
In practice, the starting time is restricted by many factors, such as local authority regulations and the driver shift times. Attempts may be made to provide night time deliveries to avoid daytime congestion, but other problems like noise, safety and out-ofhours deliveries may mean that night time deliveries are not feasible.

Distance Based Approach
In the original method, the customer sequence is determined by approximate fuel  Table 3.   As shown in the table, DBA is able to obtain a solution close to the original method for instances A_0, C_0, D_0 and E_0. The performance of DBA for instance B_0 is not as good as for other instances, but it is only 1.6 % worse than the solution obtained by the original method. In general, minimizing distance seems to be a good criterion for obtaining a good sequence of customers, before using NHA to find the detailed paths and speeds to be used when travelling between customers.

Instances
Meanwhile, the running time for DBA is much less than the original method. The running time of DBA is about 5-6 minutes, while that of the original method is more than 20 hours. This is because NHA is only used at the end of the process, after the customers have been allocated to routes and their sequences have been determined. DBA is able to get solutions that are very close to those obtained by the original method, but it needs less than 0.5% of the running time of the original method.
For the instances without time window constraints, DBA provides good solutions in terms of running time as well as solution quality. However, DBA has the potential risk of missing time windows by using static travel time data in scheduling the customer sequence. The impact of the time window constraint will be discussed in the next section.

Time Window Constraints
All the above discussion is based on the instances without time window constraints, but only with 7 am as the earliest starting time and 5pm as the latest finish time for the drivers' working shift. In this section, the instances with time windows are tested. In order to guarantee a feasible solution for each instance, the maximum waiting time is set at 4 hours.  higher-emission strategy in order to satisfy the time window constraints or even cannot meet the time window constraints, when the time-varying speeds are taken into account.

Conclusions
In this paper, a column generation based tabu search algorithm is proposed to solve a vehicle routing problem where there are time-varying speeds and the objective is to minimize fuel emissions. The main findings can be summarized as follows:  The proposed algorithm can produce sets of routes for an urban distribution problem with a reduction in GHG emissions of about 3% compared with an approach where the objective is to minimize total time.
 Within NHA, the path selection process plays a more significant role in reducing fuel emissions, when compared with the speed adjustment process. When the speeds are allowed to be as fast as the traffic allows up to a maximum speed, the best fuel emission solution is only slightly worse than the one obtained using NHA, where the speeds are decision variables.
 By allowing a long waiting time at customer nodes, vehicles can avoid being caught in congestion, and the fuel emissions can be reduced. However, the costs of such waiting, such as drivers' salaries, have to be considered. An alternative way to attempt to avoid congestion is to depart later to avoid the morning peak period, but this may not significantly reduce the fuel emissions.
 Minimizing distance can be an effective way to obtain a good sequence of customers.
DBA, which determines the customer sequences by minimizing distance, has a good performance on instances without time windows, but there is a risk that it could fail to find a feasible solution, when it is applied with time-varying speeds. In this case, some simple repair strategy could be used to obtain a good feasible solution. 13

Results of Experiments
The results of the instances A_0, B_0, C_0, D_0 and E_0 are shown as follows. The starting time is 7 am, and the maximum waiting time is 5 minutes.