The Pollution-Routing Problem with Speed Optimization and Uneven Topography

This paper considers a joint pollution-routing and speed optimization problem (PRP-SO) where fuel costs and $\textit{CO}_2e$ emissions depend on the vehicle speed, arc payloads, and road grades. We present two methods, one approximate and one exact, for solving the PRP-SO. The approximate strategy solves large-scale instances of the problem with a tabu search-based metaheuristic coupled with an efficient fixed-sequence speed optimization algorithm. The second strategy consists of a tailored branch-and-price (BP) algorithm in which speed optimization is managed within the pricing problem. We test both methods on modified Solomon benchmarks and newly constructed real-life instance sets. Our BP algorithm solves most instances with up to 50 customers and many instances with 75 and 100 customers. The heuristic is able to find near-optimal solutions to all instances and requires less than one minute of computational time per instance. Results on real-world instances suggest several managerial insights. First, fuel savings of up to 53\% are realized when explicitly taking into account arc payloads and road grades. Second, fuel savings and emissions reduction are also achieved by scheduling uphill customers later along the routes. Lastly, we show that ignoring elevation information when planning routes leads to highly inaccurate fuel consumption estimates.


Introduction
Road freight transportation is vital to the functioning of the economy and the supply chain.However, significant negative impacts on people and the environment due to excessive energy usage and considerable greenhouse emissions need to be considered.According to the International Energy Agency (IEA, 2020), global transportation is still responsible for 24% of direct CO 2 -equivalent (CO 2 e) emissions from fuel combustion.This is especially true in a city logistics context (Savelsbergh and Van Woensel, 2016).
Over the past years, these observations gave rise to the introduction of pollution-and sustainabilityrelated aspects into traditional Vehicle Routing Problems (VRPs), popularly coined in the literature as the Pollution-Routing Problem (Bektaş and Laporte, 2011).In the literature, similar problems and definitions can be found as the Emissions Minimizing VRPs (EMVRPs) (Raeesi and Zografos, 2019) or Green VRPs (Erdogan and Miller-Hooks, 2012;Moghdani et al., 2020).These models make use of the fact that the amount of transport-related greenhouse gases (GHGs) emissions is directly proportional to the fuel consumption (Kirby et al., 2000).Multiple factors are considered, including the slope (Suzuki, 2011), vehicle speed (Demir et al., 2012), the payload (Bektaş and Laporte, 2011), traffic congestion (Franceschetti et al., 2013), driver's operating habit (Bandeira et al., 2013), and the fleet size and mix (Koç et al., 2014).
Specifically, the PRP aims to build routes that minimize an objective function integrating the vehicle's routing cost (e.g., fuel consumption, the pollution aspect) and the driving costs (e.g., vehicle usage, drivers' wage, and other direct costs aspect).This paper builds upon this literature and considers modelling the driving costs as fixed costs of vehicles, which is commonly used in heterogeneous vehicle routing problems (see, e.g., Koç et al., 2016), but often disregarded and handled as duration-dependent costs in PRP.Moreover, considering this rich version of the PRP leads to interesting new methodological challenges for speed optimization.
The majority of models describe the road angle, and thus the network topography, as one of the parameters used to formulate the instantaneous engine-out emission rate (Barth and Boriboonsomsin, 2009), but do not consider this further in the model, in the solution methodology or in the results and insights.More specifically, efficient solution methods and extensive computational experiments analyzing the effect of road gradient on fuel consumption and CO 2 e emissions are missing in the literature.One notable exception in the same application domain is the paper by Brunner et al. (2019).These authors however assume that speed is constant, which is an unrealistic assumption within a city logistics context (Franceschetti et al., 2013).Concerning vehicle speed, several studies are proposing various optimization procedures.However, there is still a need for efficient and fast speed optimization algorithms to use in exact algorithms.
The contributions of this paper are threefold.
• We introduce the road gradient to the computation of fuel consumption utilizing terrain elevation information.Despite the inclusion of road gradient into the original formulation of fuel consumption, most papers assume a constant road angle for all vehicle trips.
• A number of novel solution approaches are presented, including a branch-and-price algorithm leading to optimal solutions and metaheuristic for larger instances.We develop a novel, fast and efficient algorithm for the vehicle speed optimization.
May 20, 2021 • We generate important insights based on a broad set of instances to investigate the importance of road gradient on emissions.Additionally, we propose a number of real-life instances.
The remainder of this paper is organized as follows.Section 2 provides a brief review of related scientific literature.In Section 3, we give a detailed description of our PRP model.Section 4 shows the proposed Tabu Search metaheuristic together with a detailed description of the fixed-sequence speed optimization algorithm.Later in Section 5, we present the major components of our branch-and-price algorithm.Section 6 offers extensive computational experiments conducted, and finally, we conclude this research in Section 7.

Literature Review
The significant increasing amount of CO 2 e emissions derived from road freight transportation has certainly ignited worldwide concerns.Over the last ten years, this resulted in a large body of literature on emissions-aware transportation problems (see, e.g., the survey by Moghdani et al., 2020).The Pollution-Routing Problem (PRP) is an efficient and comprehensive formulation to address the minimization of carbon emissions.The resulting greenhouse emissions of vehicle fuel consumption are the consequence of some influential factors beyond the travel distance (Ericsson, 2001;Brundell-Freij and Ericsson, 2005).
According to Demir et al. (2014b), vehicle fuel consumption is affected by multiple factors such as speed, road gradient, road congestion, driver's operating habit, size and composition (mix) of the vehicle fleet, and payload.
As summarized in Table 1, we observe that some important factors are less studied in previous vehicle routing research, specifically the interaction of load, road gradient, and vehicle speed is missing.Road gradient.The road grade (slope) has a significant influence on both the conventional-vehicle fuel economy (Suzuki, 2011) and the electric-vehicle energy consumption (Goeke and Schneider, 2015).
The transit of a typical light-duty vehicle over a sloping road surface, with a +6 percent grade, could increase the fuel consumption by 15-20 percent (Boriboonsomsin and Barth, 2009).The influence of hilly roads is undoubtedly more significant on the fuel consumption of heavy-duty trucks.According to Davis et al. (2009), just a minor increase or decrease in road grade (1-4%) can reduce or increase fuel economy by more than 50%.In the case of electric vehicles, the energy consumption is less affected by the road gradient as opposed to conventional vehicles.However, experimental results show the impact of hilly terrain on the electric vehicle miles traveled, confirming that the electric vehicle range in mountainous landscapes is lower (Travesset-Baro et al., 2015).
The road gradient has been mainly studied in PRPs and the so-called Electric Vehicle Routing Problems (E-VRPs).Since the first mathematical formulation of PRPs (Bektaş and Laporte, 2011), the road angle was one of the parameters used to define the instantaneous engine-out emission rate.It is worth mentioning that the present work is the first study in the area to create and optimally solve a set of realistic problem instances, which include elevation information for computing the road slopes along the paths.
As mentioned earlier, we found the paper by Brunner et al. (2019) as the only contribution that studied the influence of road grade on fuel consumption.Although these authors considered the slope in their VRP formulation, they assumed a constant road grade for all arcs that form the directed graph.The above does not reflect accurately the hilly topography profile of a road, which typically connects two nodes (arc) with a sequence of multiple uphill/downhill segments.Another assumption in this contribution is related to the vehicle speed.The authors addressed a VRP without time windows, assuming that the vehicles travel through each arc with a given constant speed (input parameter).Last but not least important with respect to the paper assumptions, in Brunner et al. (2019), the payload carried by the vehicle is chosen from a prescribed set of values which means that the authors did not consider the payload as a continuous decision variable.
There are very few papers, outside the application context of PRPs, that take into account the effect of road gradients on fuel consumption.For instance, Tavares et al. (2009) utilizing an exponential regression model (COPERT-III method, introduced by Ntziachristos et al. (2000)) to estimate the minimum fuel consumption during waste collection process.They studied two realistic routing problems showing that the optimal route does not necessarily correspond to the shortest traveled distance.Their computational results demonstrated that significant fuel consumption savings are possible for longer routes with moderate road inclination.Moreover, the contribution of Suzuki (2011) proposed a linear regression model to compute the fuel consumption rate for a heavy-duty truck, based on the fuel-efficiency study of Davis et al. (2009).The author included the road-gradient factor as one of the objective function components (distance and fuel consumption), designed to formulate a traveling salesman problem with May 20, 2021 time windows.
Concerning the E-VRPs, Yang et al. (2014) performed a numerical simulation.In their paper, the authors used the electric vehicle's battery (physical model) theory to study the effects of the road's slope on electricity consumption for both uphill and downhill paths.They concluded that with the increase of the uphill's tilt angle, each electric vehicle's electricity consumption increases significantly.Goeke and Schneider (2015) also assumed not-flat terrain with grades in their energy consumption model proposed for electric vehicles.They intensely focused on the effect of load distribution on the performance of commercial electric vehicles.Later, Liu et al. (2017) also investigate the impact of road gradient on the electricity consumption of electric cars.Using 12 gradient ranges and GPS tracking data with a digital elevation map, the authors showed that, in uphill trips, the energy consumption increases almost linearly with the absolute gradient.However, the numerical results also exhibited the positive effect of regenerative braking power acquired during the downhill trips.Finally, Macrina et al. (2019) modeled a comprehensive energy consumption function considering the road gradients.Nevertheless, in their computational study, the authors set the road angle equal to zero for conventional and electric vehicles.
Vehicle speed.Bektaş and Laporte (2011) first addressed the speed of the vehicle in the area of PRPs.In the original formulation of this problem, they explicitly assumed that speed over each arc is chosen from a predefined list of possible values.Koç et al. (2014) also utilized the same discrete speed function but analyzed, for the first time, the effect of fleet size and mix in PRPs.The discretization of vehicle speed was also adopted by Eshtehadi et al. (2017).Here, the authors considered the demand and travel time uncertainty, both aspects addressed by several robust optimization techniques.Moreover, Demir et al. (2012) proposed a specialized speed optimization algorithm (SOA), which computes optimal speeds on a given path so as to minimize fuel consumption, emissions and driver costs.The authors modified the original SOA, which was initially designed for solving the tramp speed optimization problem (Norstad et al., 2011).The speed and traffic congestion were also studied by Franceschetti et al. (2013), originating the time-dependent PRP.This research considers two phases within the planning horizon, the free-transit phase, and the congested phase.One interesting fact of their computational results is that they reduced the emissions cost by waiting at specific locations (stopped vehicles) and avoiding traffic congestion.Kramer et al. (2015) also modified the original SOA providing a new speed and departure time optimization algorithm.
Finally, a huge variety of solution methods has been proposed to solve PRPs.One can find the frequent application of metaheuristics such as ALNS (Demir et al., 2012(Demir et al., , 2014a)), evolutionary and genetic algorithms (Koç et al., 2014;Rauniyar et al., 2019), and hybrid approaches (Tirkolaee et al., 2020).Table 1 shows the methods used for solving different variants of PRPs.As can be seen, the utilization of exact algorithms in PRPs is very limited.Existing exact methods often approximated or discretized the vehicle speeds in order to reduce the complexity.Fukasawa et al. (2016)  method for linearizing the original fuel consumption equation (Bektaş and Laporte, 2011), and solved the PRP instances with up to 25 customers.
To the best of our knowledge, Dabia et al. (2017) is the only previous contribution that addressed a complex variant of PRP, and developed an exact branch-and-price algorithm.To address speed optimization within the pricing subproblem, the authors introduced a ready-time function for updating the speed-dependent routing costs within a bidirectional labeling algorithm.We propose an improved branch-and-price method that uses a novel speed optimization algorithm.
In Dabia et al. (2017), a complex ready time function has to be solved recursively for determining the optimal speed that minimizes the routing cost of a partial path, which can be computationally expensive.In our proposed labeling algorithm, the optimal speed can be determined efficiently and without full "backtracking" of the partial path.In addition, we have also developed completion bounds that further improved the efficiency of solving the pricing subproblem.

Problem Formulation
In Section 3.1, we formulate the carbon dioxide equivalent (CO 2 e) emission using the comprehensive modal emissions model (CMEM).In Section 3.2, we introduce a mixed-integer linear programming formulation for the joint Pollution Routing and Speed Optimization problem (PRP-SO).

Modelling CO 2 e emissions
We model CO 2 e emissions using the CMEM (refers to e.g.Barth and Boriboonsomsin (2009); Boriboonsomsin and Barth (2009); Demir et al. (2012)).The parameters and the values for light, medium, and heavy-duty vehicles (denoted respectively as LDV, MDV, HDV) we used in our experiments are shown in Table 2.According to the CMEM model, the instantaneous fuel use rate of a vehicle k when traveling at a May 20, 2021 constant speed ν k with payload x on a path with the road angle φ is given by When traversing a distance of d meters, the amount of fuel consumption is therefore given by Now that, there are a sequence of road segments associated on each arc e which are denoted by S e .
Let d es denote the travel distance of segment s ∈ S e , and x ke denote the payload of vehicle k when traversing arc e.For traversing the sequence of road segments associated on arc e, the amount of fuel consumption of vehicle k can then be determined by where To speed up the CO 2 e emission calculations, we precompute the parameters α ke , β ke , and γ ke for all the vehicles k and arcs e.These parameters are used in formulating the mathematical model in Section 3.2 and developing the solution approaches in Section 4 and Section 5.

Mathematical model
Let K be the set of vehicles.Let G(N , A) be the underlying directed graph.A set of nodes N = {0, 1, ..., n} contains n customers and a depot (represented by node 0).The set of arcs A, defined as  to the following constraints: i) the total demand in a vehicle does not exceed the vehicle capacity; ii) every route starts and ends at the vehicle's home depot; iii) every customer is visited exactly once by exactly one vehicle; iv) all vehicles should return to its home depot within a time limit; v) every vehicle travels at a speed within the speed limits.
The joint Pollution Routing and Speed Optimization problem (PRP-SO) can be formulated as the following nonlinear integer programming model.
Constraints ( 8) are the flow conservation constraints.Constraints (9) ensure that the payload does not exceed vehicle capacity.Constraints ( 10) -( 12) ensure that customers are visited within the given time windows.Constraints (13) ensure that vehicles travel at a speed within the limits.

Metaheuristic
Metaheuristic approaches require frequently evaluating solutions with fixed vehicle routes.We present an efficient algorithm for finding the optimal vehicle speed of a given route so that we can compute the costs of CO 2 e emissions and fuel consumption efficiently.In Section 4.1, we present a novel polynomialtime algorithm for solving the fixed-sequence speed optimization subproblem -determine the optimal speed of a vehicle when the customer sequence is fixed.To demonstrate the effectiveness of the speed optimization algorithm, it is embedded into a Tabu Search (TS) metaheuristic for our experiments.
TS is originally proposed by Glover (1986) as a synthesis of the perspectives of operations research and artificial intelligence.TS has been widely used and has been shown to be effective for finding near-optimal solutions to vehicle routing problems, see, e.g.Gendreau et al. (1994), Toth and Vigo (2003), Lai et al. (2016).Review on the recent development of TS can refer to Gendreau and Potvin (2019).Major components of the proposed TS include the fixed-sequence speed optimization subproblem in Section 4.1, the penalized objective function in Section 4.2, the initial solutions in Section 4.3, the neighborhood structure in Section 4.4, the intra-route improvement procedure in Section 4.5, and the search procedure in Section 4.6.

Fixed-sequence speed optimization subproblem
The fixed-sequence speed optimization subproblem determines the optimal speed of a vehicle for a given customer sequence.Let R = (v 1 , v 2 , ..., v n ) denote a fixed sequence of nodes in a vehicle route where both v 1 and v n represent the home depot.
Without time-window constraints, vehicles should always travel at a speed that is most cost-efficient and within the speed limits of the vehicle.The most cost-efficient speed of vehicle k when there are no May 20, 2021 time window constraints is given by Since vk is independent of the vehicle route when there are no time window constraints, the optimal speed of vehicle k within speed limits [a k , b k ] can be computed by In the presence of time window constraints, a vehicle should travel at a speed that satisfies all time window constraints and incurs a minimal cost.Let σ(R) denote the lowest vehicle speed without violating any time windows associated on the nodes of route R. When a vehicle travels at the speed of σ(R) in route R, all time window constraints will be satisfied.We will show in Proposition 1 that σ(R) can be computed efficiently by where denotes the total distance from the depot to node i along vehicle route R, S ij = j−1 k=i s v k denotes the total service time from node v i to node v j−1 along vehicle route R, and that [e v , l v ] is the time window associated on node v.If σ(R) ≤ 0, no feasible speed exists due to conflicting time window constraints.
Since a vehicle can wait at the customer node if the vehicle arrives earlier than the lower bound of a time window, any speed that is greater than σ(R) also satisfies the time-window constraints in route R.
Thus, the optimal speed of vehicle k when traversing on route R is given by max The total cost (as defined in the objective function ( 5)) of a given solution can be evaluated straightforward when the optimal speeds of the vehicle routes have been found by using ( 19) and ( 20).
Proposition 1.The minimum speed without violating any time window constraints of a vehicle route R is given by σ(R) when σ(R) > 0, and can be obtained in a O(n 2 ) complexity where n is the number of nodes in route R.
Proof.Let R = (v 1 , v 2 , ..., v n ) denote the sequence of nodes in a vehicle route R with both v 1 and v n representing the home depot, and there is at least one customer node in R i.e. n ≥ 3. Furthermore, we assume that the nodes in R are not all located at the same location.For i = 2, ..., n, let ∆(i) = i−1 l=1 d v l ,v l+1 denote the total distance from node v 1 (the depot) to node v i along vehicle route R. Set ∆(1) = 0.For all i = 1, 2, ..., n and j = 2, ..., n, with i < j, let S ij = j−1 k=i s v k denote the total service time from node v i to node v j−1 .For every distinct pair of nodes i, j ∈ {1, 2, ..., n}, with i < j, the May 20, 2021 following lower bound of σ(R) can be derived: Any vehicle speed that is higher than the above lower bound induced from nodes i and j would violate one of the time window constraints associated on nodes v i and v j .
The minimum speed without violating any time window constraints in R is given by the maximum lower bounds for all distinct pairs of nodes, and thus we have Since there are n(n−1) 2 distinct pairs of nodes in a vehicle route of length n, σ(R) can be computed in To complete the proof, we will show that a pair of conflicting time window constraints exist in route R if and only if σ(R) ≤ 0. Consider two cases: i) σ(R) < 0; ii) σ(R) = 0. Case i: By definition we have ∆(j) − ∆(i) ≥ 0 for all i, j = 1, ..., n with i < j.Therefore, σ(R) < 0 iff there exist v i and v j with l vj − e vi < S ij which implies a violation on one of the time window constraints associated on nodes v i and v j .Case ii: Suppose the contrary, we consider a vehicle with a speed equals to zero and that all the time window constraints in R are satisfied.Since the due date l v = ∞ for all v, we have σ(R) = 0 iff ∆(j) − ∆(i) = 0 for all i, j = 1, ..., n with i < j.Since the nodes in R are not all located at the same location, there exist two distinct nodes v i and v j in R with d vi,vj > 0 and i < j.This implies that, there exist v i and v j in R with i < j and ∆(j) − ∆(i) > 0 which leads to a contradiction.

Penalized objective function
We allow infeasible solutions in the search space.It is implemented as a penalized objective function that is obtained by relaxing some of the constraints and incorporating them into the objective function with the use of self-adjusting penalty parameters.
If a vehicle k ∈ K is assigned to a non-empty route R, and is travelling at a speed of ν k , the travel cost c(R) can be written as where xke is the payload of vehicle k on arc e and A(R) denotes the arcs on route R. The overload P(R) representing the violation of the vehicle capacity constraints is defined as P denotes the customer nodes on route R.After incorporating the penalties for possible violations, the penalized objective function value is computed by z(R) = c(R) + ρP(R) where ρ ∈ R + is the penalty weight that is self-adjusting in the search.The penalty weight ρ is initialized as 1, and updated in every δ iterations as follows.If P(R) = 0 for all vehicle routes R, then update ρ to 0.5ρ; otherwise, update ρ to 2ρ.

Initial solution
The initial solution is created by first applying a stochastic insertion-based heuristic and then improved by using the TS procedure described in Algorithm 1 with a limited number of iterations.The number of iterations is set to I 1 n where I 1 is a user-controlled parameter and n is the number of customers.After May 20, 2021 creating ten such initial solutions, the best one is returned as the initial solution for the main search procedure described in Section 4.6.
The following stochastic insertion-based heuristic is applied.It begins by assigning exactly one arbitrarily selected customer to each of the randomly selected l vehicles, where l is a randomly generated integer between 1 and the total number of vehicles available.Then, the remaining customers are considered one by one, following a randomized order.For each customer, all possible locations in all routes are evaluated for insertion, and the customer is subsequently inserted into a vehicle route at the position that minimizes the insertion cost (the incremental change in the penalized objective function value).
We determine the insertion costs by using the algorithm described in Section 4.1 that performs speed optimization for a fixed sequence of nodes on a route.

Neighborhood structure
Let X denote the set of all feasible solutions for the instance.We define the solutions in the neighborhood of a given solution x ∈ X as N(x).At each iteration, all possible move operations for all customers are evaluated and the best one is subsequently performed.The move operation involves relocating a customer from its current route to another route at the location that minimizes the insertion cost.The insertion cost is determined by applying the speed optimization algorithm described in Section 4.1.The best move operation is the one that leads to the lowest total penalized objective function value and the following diversification penalty.For a given solution x ∈ X , we define the diversification penalty as φ(x) = λc(x) √ nϑ ir where n is the number of customers, ϑ ir counts the number of times customer i has been moved to route r so far in the search, and λ is a positive parameter that controls the intensity of diversification.Readers can refer to Soriano and Gendreau (1996) for extensive discussion on diversification schemes.
To prevent cycling, if a customer has been moved from route r to route s in a given iteration, then moving the same customer back to route r is declared tabu which implies that this reverse operation is forbidden for the next h log 10 (n) iterations where h is a user-controlled parameter and n is the number of customers.To prevent the search from stagnating, the following aspiration criterion is applied: a Tabu move is allowed only when the resulting solution is feasible and has an objective function value that is better than that of the current best feasible solution found by the search.

Intra-route improvement procedure
The following procedure is applied for improving the solution by modifying customers' position within the same route: a customer is randomly picked and then reinserted into the best location of the same route, until no further improvement is possible.The intra-route improvement procedure is invoked after the selected move operation is performed, the parameter for overload penalty is updated, or when the best feasible solution is improved.

Search procedure
The TS routine is presented in Algorithm 1.The search procedure consists of two phases.The first phase of the procedure starts by constructing an initial solution as described in Section 4.3.Next, the May 20, 2021 best solution found in the first phase is improved in the second phase by executing I 2 iterations of the TS routine.
3: determine z(x) and φ(x) for all x ∈ N (x) 4: while stopping condition is not satisfied do 5: -minimizes z(x) + φ(x) 7: -x is non-tabu or it satisfies the aspiration criteria 8: set the reverse move tabu for θ iterations 9: perform the intra-route improvement procedure on x. 10: if x is feasible and z(x) < z * then set x * = x and z * = z(x). 11: set x = x, and update the penalty weight of overload for every δ iterations 12: update z(x) and φ(x) for all x ∈ N(x) 13: return x * Algorithm 1 shows the TS procedure.In line 2, the algorithm starts with an initial solution x 0 which is obtained by running the TS routine with I 1 iterations as described in Section 4.3.In line 3, the penalized objective function and diversification penalty associated on all the move operations are updated.The loop in lines 4-12 is then invoked and stops until after I 2 iterations.In line 5-7, the best neighborhood solution is picked which takes into account the diversification penalty, overload penalty, solutions in the Tabu list, and the aspiration criteria.As described in Section 4.4, the diversification penalty φ(x) depends on the intensification parameter λ.In line 8, the reverse move is forbidden for the next h log 10 (n) iterations.As shown in lines 9-10, the intra-route improvement procedure described in Section 4.5 is performed on the two routes that are modified in the best neighborhood solution.The incumbent is updated if a new best feasible solution is identified.In lines 11-12, the search moves to the selected neighboring solution.The penalty weight of overload is updated in every δ iterations.Whenever the search moves to a neighboring solution, the penalized objective function and diversification penalty associated on each of the move operations are again updated by using the speed optimization algorithm described in Section 4.1 and the penalized objective function defined in Section 4.2.The computational time is manageable since we only need to update the costs of the move operations that are involved with the modified routes.The values for the algorithmic parameters in the search procedure include I 1 , I 2 , λ, h and δ which are set according to the parameter tuning experiment described in Section 6.3.

Branch-and-price Algorithm
Branch-and-price is a successful exact method for solving several VRP variants.It relies on reformulating the problem as a set-partitioning (SP) problem and applying branch-and-bound to solve the SP formulation, which is known to provide tight linear bounds.As the number of variables in the SP model grows exponentially with the instance size, column generation is applied to identify profitable variables and add them to the model dynamically.For a review of branch-and-price applied to vehicle routing we refer to Costa et al. (2019).In this section, we present the SP-based PRP-SO formulation and describe our branch-and-price solution method, with a focus on the specialized algorithm for solving the non-linear PRP with Speed Optimization and Uneven Topography May 20, 2021 column generation subproblem.

Set-partitioning formulation
A route p is feasible for vehicle k when i) there exists a value ν ∈ [a k , b k ] such that no time-window constraint is violated when vehicle k executes route p at speed ν; and ii) the total demand of the customers served along p does not exceed the vehicle capacity Q k .Let Ω k be the set of feasible routes for vehicle k.Furthermore, for each p ∈ Ω k , let c k p be the cost incurred when vehicle k executes route p at the optimal speed, as defined in Section 4.1.Finally, for each route p ∈ Ω k and customer i ∈ N \ {0}, let We define binary decision variables λ k p for all k ∈ K and p ∈ Ω k , such that λ k p = 1 if and only if vehicle k executes route p in the solution.Then, the PRP-SO is formulated as the following SP problem: (SP) : min The objective function (21) minimizes the total costs.Constraints (22) ensure that each customer is visited exactly once, and constraints (23) enforce a maximum of one route per vehicle.

Column generation subproblem
The restricted master problem (RMP) refers to the linear relaxation of ( 21)-( 24) with a restricted set of vehicle routes.Let π i , i ∈ N \ {0}, be the dual prices associated with constraints (22) after solving the RMP.The pricing problem for vehicle k is defined as follows: (PP) : min s.t.
e∈δ + (i) x e − e∈δ + (i) x e = q i z i , ∀i ∈ N \ {0}, (28) May 20, 2021 The objective function ( 25) minimizes the route-dependent costs (including vehicle fixed cost and CO 2 e emissions cost) and the dual prices of the RMP.Constraints ( 26)-( 27) are the flow conservation constraints.Constraints ( 28)-( 29) keep track of the payload in the vehicle along each arc traversed.
Finally, constraint (34) ensures that the vehicle travels at a speed within the prescribed limits.

Pricing algorithm
The pricing algorithm identifies columns with a negative reduced cost, that is, it finds solutions to (PP) with a negative objective value.As common in branch-and-price for vehicle routing, we solve the pricing problem with a labeling algorithm.A label represents a partial route from the depot to a customer.Label extensions are created by extending labels to all feasible customers.Table 4 summarizes the attributes of a label alongside their corresponding initialization values and updating rules.A key component of our pricing algorithm is the set of label extension procedures that do not require May 20, 2021 full backtracking to determine the cost of a route under the optimal speed.The idea is formalized with the following definition and proposition.
Definition 2 (Cost of a path).Let Q be a path with arcs (e 1 , e 2 , ..., e L ) and nodes (v 0 , v 1 , v 2 , ..., v L ) where v 0 = v L = 0 is the depot.The cost of path Q when travelling at a speed of ν is defined as where Proposition 3. Let P be a complete path that ends at the depot, and let p be the route induced by P .
Then, the cost of p executed under the optimal vehicle speed, given by c p , is equal to the cost of P as determined by ( 39) where ν = ν P .
Proof.Let Q be a partial path of vehicle k with arcs (e 1 , e 2 , ..., e L ) and nodes (v 0 , v 1 , v 2 , ..., v L ) where v 0 = v L = 0 is the depot with demand q 0 set to 0. The travel cost of the route induced by Q, according to the objective function ( 25), is given by where x e l = L−1 m=l q m is the sum of the demand of the customers in the remaining part of the route, that is, the payload of arc e l .As shown below, the cost calculation by ( 40) is equivalent to (39).
The travel cost (40) can be rewritten as δ Q , and hence equivalent to the routing costs defined in the objective function ( 25).
The pricing problem can be considered as a variant of the resource-constrained shortest-path (RCSP) problem (see e.g.Feillet et al., 2004).Typically, one finds the RCSP with a labeling procedure where dominance rules are employed to discard non-promising partial paths.In our case, however, dominance rules are likely to be ineffective because of the large number of resources involved.More specifically, in addition to the usual resources to handle vehicle capacity and time windows, in our case one must also observe dominance conditions on the allowed vehicle speed range and each cost component individually.
Therefore, in our pricing algorithm, we decide to control the combinatorial growth of labels exclusively with completion bounds, which are detailed next.

Completion bounds
A completion bound is a lower bound on the reduced cost of all routes that can be generated from a label.
Completion bounds accelerate the solution of the pricing problem since partial paths with nonnegative bounds are discarded during the labeling procedure.
Consider a partial path P with arcs (e 1 , e 2 , ..., e L ) and nodes (v 0 , v 1 , v 2 , ..., v L ).The precise cost along P depends on the customers visited after v L , since the cost along each arc depends on the arc payload.
A lower bound on the cost along P , however, can be computed as follows: where Equation ( 42) considers a best-case (i.e., cost-minimizing) scenario concerning the payload along each arc.If the corresponding β is nonnegative, the vehicle is assumed to travel along arc e as light as possible.
Otherwise, the vehicle is assumed to travel as loaded as possible.
Given the lower bound (41), we propose two completion bounds for a label P .The first bound is based on a RCSP and explores the capacity resource to find a lower bound on the reduced cost of any extension of P .The second bound is based on a knapsack problem and explores not only the capacity resource but also the "timing" resources, that is, the fact that customers cannot be visited after their time windows and the vehicle must return to the depot no later than instant l 0 .
We start with the RCSP-based completion bound, which adapts the bound proposed by Florio et al. (2021) for solving the elementary RCSP.First, we associate to each arc e = (i, j) ∈ A a lower bound φe on the reduced cost change when partial path P is extended along e = (i, j): where x e = 0 if β e ≥ 0 and x e = Q k otherwise.We denote by S * i (Q) the value of the RCSP from node i ∈ N \ { } to node 0 in a graph with arc costs given by ( 43), in which the initial resource limit is Q and an amount q j of resource is consumed each time node j = 0 is visited.Then, the RCSP-based completion bound is given by: Equation ( 44) yields a valid bound because Φ P − i∈N P π i is a lower bound on the reduced cost of partial path P , and S * n P (Q k − q P ) is a lower bound on the reduced cost of any feasible extension to P .To enable evaluations of (44) in constant time for any label P , at the beginning of an iteration of the pricing problem we pre-compute S * i (Q) for all i ∈ N \ { } and Q ∈ {0, . . ., Q k }.The (non-elementary) RCSPs can be solved efficiently by dynamic programming.
While the RCSP bound is computationally efficient, it does not explore time window constraints nor the fact that we price elementary routes.In the knapsack bound, we set up a {0, 1}-knapsack problem PRP with Speed Optimization and Uneven Topography May 20, 2021 with capacity of l 0 − τ P , which corresponds to the maximum remaining routing time after customer n P is served.Then, we define a set of knapsack items Each element of I is a customer that can be visited after n P , considering time window and vehicle capacity constraints.With each item i ∈ I a value v(i) and a weight w(i) are associated: The value and weight of an item i correspond to the maximum reduced cost decrease and minimum amount of the time resource consumed, respectively, when customer i is visited in an extension of partial path P .We let K * P be the optimal solution value of the knapsack problem defined above.Then, the knapsack completion bound is given by: Each time a label P is generated, we evaluate (45) and discard P if the bound is nonnegative.This evaluation requires solving the knapsack problem to optimality, which can also be achieved efficiently by dynamic programming.

Branch-and-bound
The implemented branch-and-bound framework finds an optimal integer solution to (SP) by branching on variables y e , e ∈ A, that take fractional values.We apply a semi-strong branching rule where each potential branching variable is evaluated under the current pool of columns, and the variable on which branching leads to the highest lower bound is chosen.More precisely, at a given branch-and-bound node, we let Ω R be the set of all columns generated and Y the set of arc variables that assume fractional values in the solution to the RMP.Then, we evaluate for each y e ∈ Y, where RMP(Ω, Y , Y ∞ ) corresponds to the optimal solution value of the RMP restricted to columns Ω and enforcing constraints y e = 0 for all y e ∈ Y and y e = 1 for all y e ∈ Y ∞ in addition to the branching constraints of the parent branch-and-bound node.Finally, we branch on the variable y e ∈ Y such that (46) is maximum.Note that evaluating (46) for each potential branching variable can be implemented efficiently by loading a single linear program and (re)solving it for each variable after adjusting the cost vector accordingly, by penalizing the cost of routes that do not comply with the candidate branching decision.

Test instances using real-world data
Another dataset is constructed based on the distribution network of a large international health and beauty retailer.There are 248 customers considered in our experiments which are representing the retailer's stores in Hong Kong.Geometric information is obtained by using Google Maps APIs, including the coordinates, elevations, suggested paths between the stores.Figure 1 illustrates the geographical locations of these stores.The landscape varies from fairly hilly to mountainous with steep slopes.The stores' locations are clustered, densely populated in the central areas.We will use this dataset for evaluating the potential benefits of using elevation information in optimizing the vehicle routes.Terminal which is the busiest port in Hong Kong.

Parameter tuning
The best parameters amongst the values specified in Table 6 are chosen for each instance class.

Efficiency of the approaches
Table 7 summarizes the computational results of BP and TS on the instances described in Section 6.1 with 25, 50, 75 and 100 customers respectively.Results on individual instances are shown in Appendix A. Column NI is the number of instances in the dataset class excluding the ones that no feasible solution can be found by using BP within 3 hours.Column NO is the number of instances that can be solved to optimality by using BP within the time limit.Column NV is average number of vehicle routes, column TD is average total travel cost, column AT is the average cpu time (in seconds), and column AG is the average optimality gap (in percentage).When reporting the average values, we excluded the instances that no feasible solution can be found by using BP within 3 hours.
As shown in Table 7, BP can solve instances up to 100 customers, with 61 instances (53%) solved to optimality.BP can find all the optimal solutions of the dataset with 25 customers within a reasonable time, and most optimal solutions of the dataset with 50 customers within 3 hours.As compared to TS, BP can determine better solutions on smaller instances, but at the expense of significantly more CPU time.TS can find near-optimal solutions for all instances within one minute, and outperforms BP on solution quality for the dataset with 100 customers.

The value of using elevation information
For evaluating the potential benefits of using elevation information in optimizing the vehicle routes, the real-world instances described in Section 6.2 are solved by using TS.Table 8 summaries the total distance (in km), average speed (in km/h), total fixed cost (in EUR), total fuel cost (in EUR), and the total elevation (in meters).Figure 2 illustrates an example vehicle route.The same dataset is solved again with which elevation information is ignored when planning the vehicle routes.The results are reported in Table 9.This is achieved in our experiment by setting the angles of all slopes to zero when optimizing the vehicle routes by using TS, and evaluating the solutions with the correct elevations and slopes after the vehicle routes have been decided.We can observe the impacts on solutions when elevation information is used for planning the vehicle routes by comparing the results in Table 8 and Table 9.As shown in the experimental results, when elevation information is considered in planning, fuel consumption decreased by 54% on average while the average vehicle speed and elevation increased slightly, and the total travel distance increased by 31%.Larger travel distances and lower fuel consumption seem to be contradictory for typical vehicle routing problems.This is because differs from typical vehicle routing problems, the fuel consumption now May 20, 2021 For typical vehicle routing problems or when vehicle routes are planned manually, travel distance is usually minimized in the objective function.Our experimental result reveals that this can lead to a suboptimal solution in reality.With the elevation information, the optimal solution can balance the tradeoff between the energy consumption due to longer distances and the higher payload when going uphill.To avoid high fuel consumption when vehicles going uphill, heavier items tend to be delivered first before going uphill.Customer time windows have to be taken into account so that vehicles do not need to speed up to meet the due times which can result in high fuel consumption.If elevation information is ignored when planning the vehicle routes, fuel consumption due to payload is underestimated when vehicles are going uphill, which leads to poor solutions.With the significant savings we observed from the experimental results, logistic service providers should consider using elevation data for planning their vehicle routes in practice.

Impact of payloads and slopes
For evaluating the impact of payloads and slopes on the optimized solutions, the real-world instances described in Section 6.2 are modified by scaling the payloads by a factor r 1 ∈ {0, 0.1, ..., 1} when calculating the costs, and scaling the slopes by a factor r 2 ∈ {0, 0.1, ..., 1}.In our experiments, we replace the payload x ke in (1) by a factor r 1 x ke when calculating the costs, and replace the slope φ es in (3) by r 2 φ es when preprocessing the data.A higher value of the payload factor (r 1 ) represents scenarios when relatively heavier items are shipped, and a higher value of the slope factor (r 2 ) represents more hilly areas where steep slopes commonly appear.There are 1089 instances tested in total, and each of them is solved by using TS with a timelimit of 300 seconds.Table 10 summarizes the average costs (and refer to appendix B for the complete results).Figure 3 shows the average cost with varying payloads (r 1 ) and slopes (r 2 ) respectively.From the experimental results, we can see that the shipping cost increases with the payload (r 1 ) and decreases with the slope (r 2 ).It is more costly to ship with a heavier payload (higher value of r 1 ) with a percentage increase in total costs up to 3.49% on average regardless of the slopes.It is less costly to ship in more hilly areas (higher value of r 2 ) with a saving up to 11.71% of the total costs on average regardless of the payloads.The impact due to slopes is significantly higher than the impact due to payloads, which reveals the importance of taking into account slopes when optimizing the vehicle routes.

Conclusions
This paper formulates and proposes efficient solution methods for a joint Pollution Routing and Speed Optimization Problem (PRP-SO), where the total travel cost is a function of fuel consumption and CO 2 e emissions and depends, simultaneously, on road grades, arc payloads, and vehicle speed.The introduction of the vehicle speed as continuous decision variables results in more complicated optimization subproblems in the presence of time window constraints.For the fixed-sequence speed optimization problem where the vehicle route is known, the proposed approach is conceptually simple and computes the optimal vehicle speed (with and without time windows) in quadratic time.In the speed optimization with variable routes, we introduced a novel labeling algorithm, without full backtracking searching, that efficiently determines the cost of a vehicle route that travels with optimal speed.Based on the proposed speed optimization algorithms, we present two general solution approaches for solving the PRP-SO.An approximate solution strategy aims to solve large instances in a short computational time.For this purpose, we integrated the fixed-sequence speed optimization algorithm in a tabu search metaheuristic.The second approach consists of an exact branch-and-price algorithm, in which the variable-route speed optimization is managed within the pricing problem.
We carried out extensive computational experiments on modified Solomon benchmarks and newly constructed real-life instances.Numerical results show that the exact solution methodology performs very well in terms of solution quality: 61 out of 116 instances are solved to optimality.Contrary to the computational outcomes presented by Dabia et al. (2017), in which several instances with only 25 customers cannot be solved, we are able to solve all small-scale problem instances (25 customers) within a reasonable time.Our BP algorithm solved most of the problem instances with 50 customers and reached optimal solutions for some larger benchmark instances with up to 100 customers.We show that our metaheuristic works very effectively for all instances solved.The heuristic consumes less than one minute to find near-optimal solutions in all instances and improved best-known solutions where the exact algorithm did not reach optimality.
Our computational results on real-world instances provide sufficient evidence to suggest some essential managerial insights.First, significant savings (53%) in fuel consumption and CO 2 e emissions are observed, especially when shipping heavy items in hilly areas.Second, vehicle routes included a larger number of customer visits (located at flatter terrain) before going to uphill destinations, also significantly reducing fuel consumption.Third, if elevation information is ignored when planning vehicle routes, fuel consumption estimation is inaccurate.
Finally, we do believe that pending issues are requiring future research.From a modeling perspective, the route security (uphill and downhill paths) in hilly topographic cities is a challenging issue to be included.Future studies can now be focused on generalizing the proposed methodological approaches to a new setting where two metrics of performance, fuel consumption, and vehicle route security should be optimized.
May 20, 2021 represents the paths between the nodes.Each node i ∈ N is associated with a demand q i and a time window [e i , l i ].Each vehicle k ∈ K is associated with a fixed cost f k , a variable cost c k (including the costs for CO 2 e emission and fuel consumption), a vehicle capacity Q k , vehicle curb weight w k , and speed limits [a k , b k ].Each arc e ∈ A is associated with a distance d e , and the parameters α ke , β ke and γ ke described in Section 3.1 for estimating fuel consumption.For all k ∈ K and i ∈ N , let z ki be a binary decision variable, with z ki = 1 if and only if customer i ∈ N \ {0} is served by vehicle k; and with z k0 = 1 if and only if vehicle k is in use.For all k ∈ K and e ∈ A, let y ke be a binary decision variable with y ke = 1 if and only if vehicle k traverses arc e, and let x ke denote the corresponding payload when vehicle k traverses arc e.For all k ∈ K and i ∈ N , let t ki denote the time at which vehicle k starts serving customer i ∈ N \ {0}; and let t k0 denote the time vehicle k returns to the depot.For all k ∈ K, let ν k denote the speed of vehicle k.The objective is to minimize the total CO 2 e emission costs, fuel costs, and vehicle fixed costs, subject PRP with Speed Optimization and Uneven Topography May 20, 2021 = αP + α ke βQ = βP + β ke γQ = γP + γ ke δQ = δP + βQqj a Assuming label P ; b Assuming vehicle k and a label representing the partial route along arc e = (0, i); c Assuming label Q obtained by extending P along arc e = (i, j).

Figure 2 :
Figure 2: An example vehicle route

Figure 3 :
Figure 3: Impact of payloads and slopes on routing costs

Table 1
clearly shows that most PRP-previous contributions included the road angle in the mathematical formulation of fuel use rate.This table also stands out two major gaps in the area that the proposed work is trying to address: (i) despite the inclusion of road gradient into the original formulation of fuel consumption, most of the papers assume that the road angle remains constant throughout all vehicle trips, and (ii) the majority of previously generated problem instances have not yet considered the road gradient.

Table 3 :
Notation for the MINLP Model ke A constant for estimating the fuel consumption of vehicle k on arc e R + β ke A constant for estimating the fuel consumption of vehicle k on arc e R + γ ke A constant for estimating the fuel consumption of vehicle k on arc e R + z ki Decision variable, to allocate customers to vehicles B y ke Decision variable, to allocate arcs to vehicles B x ke Decision variable, payload of vehicle k when traversing on arc e R + t ki Decision variable, time at which vehicle k starts serving customer i R + ν k Decision variable, speed of vehicle k R +

Table 5 :
Vehicle Information For example, the vehicles in the Solomon's instances with a capacity of 200 units correspond to the LDV vehicles of the PRP-SO instances, and therefore a demand of 20 units in those instances represents a demand of 120 kg (given by 20 × 6) in the PRP-SO instances.All the vehicles have a fixed cost of 100 EUR, a fuel cost of 1.42 EUR per liter, and a maximum speed of 80 km per hour.Other parameters used for the CO 2 e consumption calculations are summarized in

Table 7 :
Computational Results of the PRP-SO Instances

Table 8 :
Results on Real-world Instances Instance Distance Speed Fixed cost Fuel cost Elevation

Table 9 :
Results on Real-world Instances when Slopes are Ignored Instance Distance Speed Fixed cost Fuel cost Elevation depends not only on the distance, but also on the slopes, payload, and vehicle speeds.Optimal vehicle routes should therefore save fuel costs by avoiding going uphill at a high speed with a large payload.As a result, vehicles tend to visit more customers first before going uphill.Although this will increase travel distance, fuel consumption can be saved from going uphill with less payload.

Table 10 :
Impact of Payloads and Slopes on Routing Costs

Table 12 :
Results on dataset A with 50 customers