Vehicle Routing Problems with Fuel Consumption and Stochastic Travel Speeds

Conventional vehicle routing problems (VRP) always assume that the vehicle travel speed is fixed or time-dependent on arcs. However, due to the uncertainty of weather, traffic conditions, and other random factors, it is not appropriate to set travel speeds to fixed constants in advance. Consequently, we propose a mathematic model for calculating expected fuel consumption and fixed vehicle cost where average speed is assumed to obey normal distribution on each arc which is more realistic than the existing model. For small-scaled problems, we make a linear transformation and solve them by existing solver CPLEX, while, for largescaled problems, an improved simulated annealing (ISA) algorithm is constructed. Finally, instances from real road networks of England are performed with the ISA algorithm. Computational results show that our ISA algorithm performs well in a reasonable amount of time.We also find that when taking stochastic speeds into consideration, the fuel consumption is always larger than that with fixed speed model.


Introduction
The vehicle routing problem (VRP) is one of the most important and studied combinatorial optimization problems [1].It has gained great attentions from many researchers, especially in distribution and logistics management fields.More than fifty years have elapsed since Dantzig and Ramser [2] firstly introduced the problem in in 1959.They described a real-world application concerning the delivery of gasoline to service stations and proposed the first mathematical programming formulation and algorithmic approach.Following this seminal paper, hundreds of models and algorithms were proposed for the optimal and approximate solution of the different versions of VRP.Then in order to be better aligned with the real-world applications, many different VRP versions have been studied.The most common version is the CVRP [3], where each vehicle has limited load capacity.The VRP with Time Window (VRPTW) [4][5][6] aims to find optimal route sets with minimum total travel cost while serving each customer within specified time window.Other extended versions include VRP with backhauls [7], VRP with pickups and deliveries [8], and the multidepot VRP [9].
The deterministic VRP cannot cover all the situations in reality while considering stochastic VRP components.Consequently, Stochastic VRP are developed.For example, Ritzinger et al. [10] made a detailed review of the stochastic VRP.Mehrjerdi [11] combined chance constrained programming and multiple objective programming to obtain satisfactory solutions.Tas ¸et al. [12,13], Ehmke et al. [14], and Laporte et al. [15] put forward different heuristic methods for VPR with stochastic travel time and time windows.Marinakis et al. [16] developed a particle swarm algorithm for VRP with stochastic demands.For more stochastic demand results, one may refer to [17][18][19].
Based on the NP hardness of VRP, multiple heuristic algorithms have been put forward to solve this problem.Xiao and Konak [20] present simulating annealing algorithm to solve the green vehicle routing and scheduling problem with hierarchical objectives and weighted tardiness.Kondekar et al. [21] provide a mapreduce based hybrid genetic solution for solving large-scale vehicle routing problems in dynamic network with fluctuant link travel time.Neural network is also applied to solve stochastic multiconstraint problems with different time-scales; see Zhang et al. [22] and Meyer-Bäse et al. [23].Demir et al. [24] proposed an adaptive large neighborhood search algorithm (ALNS) to minimize the fuel consumption and the driving time with Pareto optimality.Garaix et al. [25] proposed a column generation algorithm for a dial-a-ride problem.
Due to the uncertain traffic factors, such as unexpected workload or bad weather, stochastic VRP has attracted more and more attentions.Cao et al. [26] proposed a partial Lagrange multiplier method.Ishigaki [27] considered a dynamic collection plan with stochastic demand and apply their search algorithm to an actual trash collection problem.Novel applications of the well-established problem can also be found in optical flow and vehicle systems [28][29][30].
In most literatures, travel speeds are assumed to be fixed or time-dependent (see, e.g., [31,32]).However, in practice, due to the uncertainty of weather, traffic conditions, and other random factors, it is not appropriate to set travel speeds to be fixed constants in advance.The interest in stochastic VRP in this paper is motivated by both its practical relevance and its considerable difficulty: large VRP instances may be solved to optimality only in particular cases.Therefore, we study vehicle routing problems with stochastic travel speeds.Moreover, as each arc has limited speed and other random factors, the average speed of the same type of vehicles on the same arc approximately obeys normal distribution.
Figure 1 is an illustration of this VPR with stochastic average travel speed.Assume that a logistics company owns a fleet of trucks and one truck delivers goods for customers 1, 2, and 3. On the first day, the truck travels with an average speed of 25 m/s on arc(0, 1).On the second day, it rains heavily when it delivers goods for customer 1.The truck has to get over the poor road conditions caused by the heavy rain; as a result, the average speed it travels on arc(0, 1) becomes 15 m/s.Several days later, a traffic accident occurred on arc(0, 1), so it may travel with an average speed of 10 m/s.Practically, during a relatively long period of time, the average speed on arc(0, 1) is not fixed but with slight fluctuation.From this point of view, the average speed is a stochastic variable.Here we assume that the average speed on each arc follows normal distribution.
On the other hand, environmental issues have become worldwide problems.In fact, fuel consumption of traffic vehicles is a great contributor to CO 2 emissions which is the main culprit of global warming [33,34].As a result, we take the sum of the expected fuel consumption and total vehicle cost as the objective to evaluate different routes.It has been studied that the fuel consumed of a vehicle traveling along a route depends on many factors, which include distance, load, speed, road conditions, and vehicle types.Models can be found in Bektas ¸and Laporte [34] and Xiao et al. [35].
Table 1 summarizes and compares the formulation of VRP in the stream of uncertain VRP and green VRP.
Two main contributions of this paper are described as follows: (1) Based on the fuel consumption model introduced by Bektas ¸and Laporte [34], the fuel consumption is a nonlinear function of travel speed, distance, and vehicle load.Due to the stochasticity of the travel speeds, we extend the model with stochastic speeds considered.
(2) An ISA (simulated annealing algorithm) is presented to solve large-scaled VRP problems.
The remainder of this paper is as follows.In Section 2, the description and formulation of this model is provided.In Section 3, we made a linearization of this problem.Then Section 4 introduces the improved simulated annealing algorithm with memory for optimizing the routing plan.In Section 5, the computational experiments are performed.Finally, Section 6 presents the conclusions and managerial insights.

Problem Statement and Model Formulation
Generally speaking, logistics companies tend to make budget decisions for a planning period; for example, they figure out the amount of the fuel consumption, the cost of keeping, and maintaining a fleet of vehicles in advance.Thus, we try to develop a fuel consumption model and algorithm to obtain more precise budgets.We use a digraph to describe the vehicle routing problems with fuel consumption and stochastic speeds (VRPFSV).Let a complete connected digraph  = (, ) be the logistics network with a node set  = {0, 1, . . ., } and an arc set  = {(, ) | ,  ∈ ,  ̸ = }.Arc(, ) represents the path from node  to node .The depot is denoted by node 0 and the other nodes in /{0} represent  customers with nonnegative demand   .The depot owns enough homogeneous vehicles with limited capacity , so the total demands of customers assigned to the same vehicle must be less than or equal to .Each route is finished by only one vehicle and each customer is served only once.The average travel speed V  on arc(, ) is stochastic.After serving all the assigned customers, each vehicle has to return to the depot.

Assumptions and Notations.
In this section, the assumptions and notations used in the model formulation are listed.First, we assume the following: (1) Each demand must be satisfied and each customer is served only once.
(3) Each vehicle must departure from the depot and after having served its customers it must return to the depot.
(4) Time window constraints are not considered here.construct the VRPFSV.The VRPFSV is an extension of the classical VRP and assumes that average travel speed V  on arc(, ) is a nonnegative random variable.Furthermore, the average travel speed V  on arc(, ) is assumed to obey normal distribution with mean   and variance  2  .The fuel consumption model applied in this paper is based on the comprehensive model formed by Barth and Boriboonsomsin [33].It is given by where  = / and  = 1/1000   are constants decided by different fuel properties,  = 0.5  , and  are vehicle characteristics related constants.  is a constant associated with the road characteristics and acceleration,   =   +  sin   +   cos   , where   is acceleration and   is the road angle inclination.The reference values of all parameters are given in Table 2. Assuming that acceleration and road inclination are zeros, then the fuel consumption on arc(, ) can be rewritten as where  1 = ,  2 =   ,  3 =   ,  4 = , and the price of the fuel  = 1.4m/.Equation ( 2) establishes a good correlation between fuel consumption and travel speed, weight, and distance.
In our research, the service cost includes two components: the expected fuel consumption and the fixed vehicle cost.
Consider arc(, ); if it is traveled by a vehicle with average speed of V  , then it has an expected fuel consumption; according to (2), the expected fuel consumption is given in Consequently, the mathematical formulation of VRPFSV can be expressed as follows: ≤   , ∀(, ) ∈  (8) ≥ 0, ∀ (, ) ∈ .
In this optimization model, the two variable sets are   and   .Equation (4) is the objective function which consists of the expected fuel consumption and the total vehicle cost. 0 = 1 indicates a departure of a newly used vehicle from the depot; as any used vehicle will trigger a fixed cost, the sum of  0 represents the total vehicle cost.Constraint (5) represents that the vehicle visits each node (except the depot) only once.Equation ( 6) represents the conservation of flow constraints.Equation (7) ensures that the demand of each customer must be met and it also eliminates subtours.Equation (8) guarantees that the vehicle load cannot exceed the vehicle capacity .Equation (9) indicates that   is a binary variable.Equation (10) is the nonnegative constraint of vehicle load.

Linearization and Solution by CEPLEX
Since the VRP is NP-hard and this VRPFSV model includes stochastic travel speed, it is at least NP-hard.In this section, we transform the original model into a linear one, and then small instances can be solved easily by CEPLEX 12.6.2.
This VRPFSV model is nonlinear due to the existence item     in the objective.In fact, as the Hessian matrix of the objective is nonpositive, the VRPFSV model is even not convex.Fortunately, we manage to make the linearization and obtain the linear model which is equivalent to the original one.
For any given arc(, ), consider its expected fuel consumption: In fact, the total fuel consumption can be reduced as follows: Property 1.In the VRPFSV, ( 11) is equivalent to (12), subjecting to constraint (8).
Proof of Property 1.In order to prove this property, we only need to prove that the product of   and   is equivalent to   .In fact,   is a binary variable; if   = 1, then     =   , else if   = 0, due to constraints ( 8) and ( 10), we get   = 0. Consequently,     = 0 =   .Thus,     is equivalent to   which indicates that ( 11) is equivalent to (12).
According to (12), objective ( 4) is converted to the equivalent linear form: Let As a result, objective (4) can be transformed into , which is a linear objective function about   and   now.Moreover, constraints (5)- (10) are linear about the decision variables, so we established a linear mixed integer programming model for VRPFSV.
Then the following intractable problem is how to calculate  1 and  2 , as it is unable to get the closed form of the integral in (14).In the following section, we use numerical integration to obtain and store the value of  1 and  2 on every arc in advance.Now all the expressions in the model are linear, so we got the coefficient value in the objective.It is possible to optimally solve problem VRPFSV by the existing solver CPLEX (version 12.6.2) for small sized problem.The linearization plays a significant role in finding optimal solution and reducing computation time.

ISA Algorithm for Large-Scaled VRPFSV
For large-scaled instances, the existing solver will go expired, as the computation process will cost a very large amount of time.We develop an algorithm based on SA algorithm.This ISA algorithm includes four parts, construction of the initial solution, generating neighborhood solutions, local search, and replacement of the current best solution.And theoretically, the SA algorithm is guaranteed to converge to the global optimization solution with probability one [43].Then make  feasible for VRPFSV; that is, assign  customers to appropriate number of vehicles.Usually, vehicle load contributes a lot to the fuel consumption.For example, if a vehicle carries a heavy load starting from the depot, it manages to serve more customers.On the contrary, if the initial load is small, more vehicles will be needed because of the capacity limit.According to vehicle capacity, string  is transformed to feasible solution with the following procedures.First, start from the second character of , and then accumulate the demands of nodes successively.If the cumulative demands exceed the vehicle capacity, start a new vehicle; that is, insert a 0 into string  before the current character and then reset the cumulative demands to 0. Repeat the above steps until the end of .Take  = 10 as an example, assuming that  = {0 3 4 1 7 5 9 8 10 6 2 0} is the initial string.As we do not know how many zeros are needed to be inserted into  initially, we set the length of VRPX to the largest likely length 21.Because the cumulative demands of 3, 4, and 1 are less than the capacity and the cumulative demands of 3, 4, 1, and 7 exceed the capacity, 0 is inserted between 1 and 7. Continue checking until the end of ; then we can obtain VRPS = {0 3 4 1 0 7 5 0 9 8 10 0 6 2 0 0 0 0 0 0 0}, which represents ten customers being visited by four vehicles.

Generate Neighborhood Solutions.
Exchange rule: three commonly used exchange rules swap, relocation, and 2opt (see Figure 2) are used to configure the neighborhood solutions.After exchanging, we get a neighborhood solution, and then make it feasible following steps of Section 4.1.Finally, the neighborhood solutions are obtained.

Local Search and Update the Current Best Solution.
The objective function value is the criterion to evaluate which solution is better.If the total cost becomes smaller, the newly generated solution is accepted.Otherwise, the solution is accepted with certain probability.As poor solutions maybe accepted at certain probability, so the best of the output cannot be guaranteed.Consequently, we use a memory array VRPS best in the algorithm to record the best solution.Only when a better solution appears, will VRPS best be updated.After continuous improvement, we can get the best solution in all searched neighborhoods.Furthermore, it is not so easy to calculate the vehicle load on each arc, so once a route is generated, the load will be calculated in a reverse order.More details are listed in in Appendix B.
This algorithm adopts adaptive cooling processing, and the temperature cooling coefficient  can be set to /( + ), where  is the accepted solution numbers of the current temperature and loop is the total loop numbers.In this way, if the number of accepted solutions is small, that is, the current solution is close to optimal solution, the coefficient may lead to smaller search scope which may cost less time.On the other hand, if the number of accepted solutions is big, a relatively big cooling coefficient will help expand the search scope.
The main steps of this improved simulated annealing algorithm are as follows.
Step 1. Initialize parameters.Input customers with demand  and distance matrix , vehicle capacity , fixed single vehicle cost , speed distribution parameters of each arc  and  2 , the end temperature  end , and the number of inner loops loop.
Step 2. Generate initial string  according to the method introduced in Section 4.1.
Step 3. Generate initial temperature.Randomly change  to its neighbor for 1000 times; the maximum objective deviation is selected as the initial temperature.
Step 4. Whether the inner loop number is reached, if it is, go to Step 8. Otherwise, generate new solutions and transform it into feasible solution.
Step 5. Calculate the deviation between the current cost and the former one and accept or reject it according to Metropolis rule.
Step 6. Update the best so far solution.
Step 7. Whether the improvement is less than 0.01, if it is, go to Step 9.
Step 8. Whether the end temperature is reached, if so, go to Step 9. Otherwise, drop the temperature and go to Step 4.
Step 9. Stop and output the best solution.The pseudocode of the improved annealing algorithm is listed in Appendix A.
In Step 3, we use self-adaptive method to obtain the best initial temperature.As too high initial temperature will increase the computational time and too low initial temperature will trap in local optimal solution.Here, we generate the maximum temperature adaptively according to the problem; that is, randomly change the incumbent solution for 1000 times.The maximum cost deviation between two neighboring solution is selected as the initial temperature.
Moreover, the temperature drops adaptively, under certain temperature; if the number of solutions being accepted is large, then the temperature drops very quickly.Else, if the number of accepted solutions is small, the temperature drops relatively slowly.
In Step 5, let   and  now be the objective value of the th iteration and the best objective of the accepted routes, and Δ =   −  now .VRPS now and VRPS new denote, respectively, the minimum cost route among the accepted ones and the th route.If Δ < 0, then replace VRPS now with VRPS new .Else, replace VRPS now with VRPS new if random(0, 1) < ( now −  end )/( 0 −  end ).
The algorithm terminates on two conditions: the current temperature is below or equal to the end temperature  end , or the improvement of the best solution is less than 0.01.Theory 2. In the VRPFSV, the time complexity of the proposed SA algorithm is ().
Proof.In the outer temperature loop, the number of temperatures in the cooling procedure is log  ( 0 / end ); in the internal loop, the number of loops is loop, and each loop has a calculation of cost function  with time complexity ().Thus, the total time complexity becomes ( loop log  ( 0 / end )).As loop, ,  0 , and  end are all given constants, Theory 2 is proved.

Computational Experiments of the ISA Algorithm
In this section, the results of computational applications of ISA algorithm are presented.First, we introduce the Pollution-Routing Problem Instance Library (PRPLIB) (http://www.apollo.management.soton.ac.uk/prplib.htm).This library contains nine differentscale groups of instances of the Pollution-Routing Problem (PRP).Each group consists of 20 different instances.These instances are based on real distances collected from UK cities.The first number in the file name after UK shows the number of nodes contained in the instance.The second is the order number of the instance within the group.The problem consists of  nodes, and the format of each file includes data about number of customers, vehicle capacity, city name and demand, distance matrix, and minimum and maximum speed level.
As there are no suitable benchmark problems of this model, we modified the benchmark problems of the PRPLIB to test our algorithm.Suppose that V  is the average travel speed on arc(, ) and it obeys normal distribution with mean parameter   and variance parameter   .
In these examples, speed limit of all arcs is 5 m/s∼25 m/s.For each arc, the integer mean parameters   are generated randomly from interval [5,25] and set   = 1.At the same time, the other parameters stay unchanged.The ISA algorithm was implemented in MATLAB and executed on an Intel 2.0 GHz processor with 1.59 G of RAM.
5.1.Quality and Efficiency of the ISA Algorithm.Tables 2  and 3 compared results of the optimal solutions found by CPLEX, SA, and ISA algorithm when  is small.The SA and ISA algorithm are performed 10 times for each instance.The average and best objective value of the 10 runs are listed in the tables, respectively.We label the metrics associated with each objective as follows: the first column is the instance ID, and the second and third columns are the optimal objective value and computation time of solutions solved by the solver CPLEX.Avg-Obj is the average total cost found by our SA or ISA algorithm in 10 runs.Best-Obj is the best objective of the 10 runs.Dev. is the relative deviations between the best objective obtained in 10 runs and the optimal cost solved by CPLEX or SA; that is, Dev = ((Obj ISA best − Obj  )/Obj  ) × 100%,  = CEPLX or SA best.Avg-CPU time is the average run time (CPU) of the SA or ISA algorithm.

Parameter Settings.
As different parameter settings may influence the performance of ISA heuristic, we test different parameter values in advance.Furthermore, the initial temperature and temperature drop factor are generated adaptively; we only need to determine the appropriate values of  and  end .We change one parameter at a time to observe how that parameter affects the solution.It can be noticed from Figure 3 that using a larger  or a smaller  end may improve the solution quality a little but increase the computation time significantly.
By trading off the effectiveness and the efficiency of the algorithm, we found the values of  = 100 and  end = 0.1 are appropriate parameters combination.
In Table 3, the CPLEX is always superior to SA and ISA in both objective and computation time, and the largest relative deviation between ISA and CPLEX is 5.88%.Table 4 depicts in most cases that the ISA manage to obtain the optimal solution and the computation time turns out to be less than CPLEX.Bold numbers indicate that, in four instances, the ISA algorithm finds the optimal solution.In all cases, ISA consumes much less time than classic SA while obtaining the same solution.
It is shown in Table 5 that when  = 25, the largest deviation of the ISA and CPLEX is 1.46%, and the ISA obtains  the near optimal solution in less time than both CPLEX and ISA.For three instances the ISA managed to get the optimal solution.
In Table 6, results from the ISA outperform the solver CPLEX.As can be seen, the biggest deviation between ISA and CPLEX is 5.25% and the longest running time of ISA is less than 60 seconds, while solver CPLEX costs more than 8 hours.On the other hand, ISA costs less time and obtain better solutions than SA.These results highlight the high quality and time advantage of the proposed SA algorithm for large-scaled problems.
For large-scaled problems, if we set the time limit of CPLEX solver to 12 hours, all the computing process terminated without any solution.Thus, Tables 7 and 8 only perform the ISA and SA results.As can be seen, the longest running time of ISA is less than 100 seconds, which illustrates that the proposed algorithm performs well especially for large-scaled problems.

The Effect of Fixed Cost and Variable Cost.
As our interest lies in the robustness of the final approach of the problem, we consider different fixed cost and variable cost values based on average over ten runs for each instance.
Take  = 200 as an example.First let  = 1, and variable  values are tested over the ranges  ∈ {0.01, 0.1, 10, 100, 1000}.The total cost and number of vehicles are listed in Table 9.
From Table 9, we notice that the total cost increase as the fixed cost increases.The number of vehicles increases as the vehicle cost increases, but the increase is not obvious.It only changes slightly when the difference between  and  is large enough.The reason is that our heuristic algorithm always tends to accept the solution with the minimum number of vehicles with a given depots visiting sequence.

Compared Results between Models with and without
Stochastic Speed.In order to observe the effects brought by  the stochastic speeds, we focus on the compared results when  = 10.In the stochastic model, we set the means of average speed on all arcs to the given value in the first column of Table 10.In the fixed speed model, the speed of each arc is set to the given value, that is, the corresponding mean in the stochastic model.Table 10 lists the objectives obtained by the ISA algorithm.
As can be seen from Table 10, over the speed interval [5,25], cost of model with fixed speed declines first; when it reaches 15 m/s, the minimum cost is obtained; then the cost begins to increase until the end of the speed interval.
In fact, fuel consumption takes up a large proportion in the total cost, so the trend of it (Figure 4) is in line with the trend of the total cost.Generally, costs obtained by the stochastic speed model is larger than that with fixed speed over the speed interval (the bold numbers in Table 10 are relatively high cost of the stochastic speed model), but when the speed is close to the end points of the interval, model with fixed speed will obtain smaller cost.The reason of this is that the integral interval will be ineffective once speeds exceed the two interval ends.

Conclusion
A new variant of the VRP with fuel consumption is presented, where the travel speed is considered to be stochastic.To solve this problem, an improved simulated annealing algorithm is proposed and presented in this paper.The proposed ISA is effective for solving the fuel VRP especially for largescale problems.Experimental results showed that when the value of the expected speed is not close to the end points of speed limit, total travel cost of the stochastic model is always larger than that of the fixed speed model.In order to make reliable logistics routing decisions, managers should take the stochasticity of speed into consideration.Several extensions are possible for further research.One worth mentioning extension is when average travel speed follows other forms of distribution, and another one would be the time window version or traffic congestion version of this problem.

2 Mathematical
Problems in Engineering

Figure 1 :
Figure 1: Illustration of a delivery.

Figure 3 :
Figure 3: Sensitive analysis of loop and  end .

Figure 4 :
Figure 4: Relationship between fuel consumption and travel speed.

Table 1 :
A comparison of existing vehicle routing models.

Table 2 :
Values and notations of parameters.

Table 9 :
Test results of different fixed vehicle cost.