# Vehicle Routing Problems with Fuel Consumption and Stochastic Travel Speeds

Authors: Yanling Feng; Ren-Qian Zhang; Guozhu Jia
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

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 large-scaled 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.


---

## Body

## 1. 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–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. Taş 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–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–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.Figure1 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.Figure 1
Illustration of a delivery.On the other hand, environmental issues have become worldwide problems. In fact, fuel consumption of traffic vehicles is a great contributor to CO2emissions 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 Bektaş and Laporte [34] and Xiao et al. [35].Table1 summarizes and compares the formulation of VRP in the stream of uncertain VRP and green VRP.Table 1
A comparison of existing vehicle routing models.
Reference Stochastic travel speed Stochastic service time Stochastic demand Objective function Weight Time window Independent decision variables Solution approach Uncertain VRP Figliozzi [36] No Yes No Vehicle cost; travel time No Yes Routes; service time New metaheuristic Errico et al. [37] No Yes No Cost No Yes Routes; service time Column generation; branch-price-and-cut algorithms Marinaki and Marinakis [38] No No Yes Expected length Yes No Route Glowworm swarm optimization Sarasola et al. [39] No No Yes Total distance; lateness penalty No No Route Variable neighborhood search Zhang et al. [19] No No Yes Travel time; earliness penalty; lateness penalty No Yes Route Dynamic programming  Green VRP Bektaş and Laporte [34] No No No Driver cost; operational cost; emissions cost Yes Yes Route; start time CPLEX solver Jabali et al. [40]; No Yes No Driver cost; fuel cost; emissions cost No Yes Routes; speed; route start time Tabu search Xiao et al. [35]; Kuo [31] No No No Fuel consumption Yes No Routes Simulated annealing Franceschetti et al. [41] No Only two fixed periods No Emissions cost; driver wage Yes Yes Routes; speed; departure time MIP & DSOP algorithm Xiao and Abdullah [42] No (time-dependent speed) No No CO2 emissions; travel time; travel distance Yes Yes Routes; travel distance in each period; time Dynamic programming This paper Yes No No Fuel consumption; vehicle cost Yes No Routes; vehicle load on each arc MIP & improved simulated annealingTwo main contributions of this paper are described as follows:(1)
Based on the fuel consumption model introduced by Bektaş 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 Section2, 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.


## 2. 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 digraphG=(N,A) be the logistics network with a node set N=0,1,…,n and an arc set A=(i,j)∣i,j∈N,i≠j. Arc(i,j) represents the path from node i to node j. The depot is denoted by node 0 and the other nodes in N/{0} represent n customers with nonnegative demand qi. The depot owns enough homogeneous vehicles with limited capacity Q, so the total demands of customers assigned to the same vehicle must be less than or equal to Q. Each route is finished by only one vehicle and each customer is served only once. The average travel speed vij on arc(i,j) is stochastic. After serving all the assigned customers, each vehicle has to return to the depot.


### 2.1. 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.(2)
The depot owns enough homogenous vehicles.(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.Then the notations used in the VRPFSV formulation are listed as follows:n: total number of customers.i: index of nodes (the depot is represented by 0).f: fixed cost of using a vehicle.Q: the capacity of a vehicle.N: set of nodes including the depot, i∈N.A: set of arcs formed by all pairs of nodes, (i,j)∈A,∀i∈N,j∈N,i≠j.lij: distance between node i and node j.vij: stochastic average travel speed on arc(i,j).Two decision variables are as follows:xij: binary variable indicating whether arc(i,j) is traveled. If arc(i,j) is traveled by a vehicle, then xij=1; otherwise, xij=0.Qij: load carried from node i to node j.


### 2.2. Formulation of the Model of VRPFSV

In this section, a mixed integer linear programming model is developed to construct the VRPFSV. The VRPFSV is an extension of the classical VRP and assumes that average travel speedvij on arc(i,j) is a nonnegative random variable. Furthermore, the average travel speed vij on arc(i,j) is assumed to obey normal distribution with mean uij and variance σij2.The fuel consumption model applied in this paper is based on the comprehensive model formed by Barth and Boriboonsomsin [33]. It is given by (1)Zijvij=λkNV+wγαijvij+γαijfijvij+βγvij3dijvij,where λ=ξ/κψ and γ=1/1000ηtfη are constants decided by different fuel properties, β=0.5CdρA, and w are vehicle characteristics related constants. αij is a constant associated with the road characteristics and acceleration, αij=τij+gsin⁡θij+gCrcos⁡θij, where τij is acceleration and θij 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(i,j) can be rewritten as (2)Zijvij=ω1vij+ω2+ω3fij+ω4vij2dij,where ω1=λkNV,ω2=λwγgCr,ω3=λγgCr,ω4=λβγ, and the price of the fuel c=1.4£/L. Equation (2) establishes a good correlation between fuel consumption and travel speed, weight, and distance.Table 2
Values and notations of parameters.
Notation Description Typical values W Curb-weight (kilogram) 6350 Ξ Fuel-to-air mass ratio 1 K Engineer friction factor (kilojoule/rev/liter) 0.2 N Engineer speed (meter/second2) 33 V Engineer displacement (liters) 5 G Gravitational constant (meter/second2) 9.81 C d Coefficient of aerodynamic drag 0.7 P Air density (kilogram/meter3) 1.204 A Frontal surface area (meter2) 3.912 C r Coefficient of rolling resistance 0.01 n t f Vehicle drive train efficiency 0.4 H Efficiency parameter for diesel engines 0.9 c Fuel and CO2 emissions cost per liter (£) 1.4 K Heating value of a typical diesel fuel (kilojoule/gram) 44 Ψ Conversion factor (gram/second to liter/second) 737 v m i n Lower speed limit (meter/second) 5 v m a x Upper speed limit (meter/second) 25 Source: Demir et al., 2012 [24].In our research, the service cost includes two components: the expected fuel consumption and the fixed vehicle cost.Consider arc(i,j); if it is traveled by a vehicle with average speed of vij, then it has an expected fuel consumption; according to (2), the expected fuel consumption is given in (3)EFijvij=∫vijminvijmaxlijw1vij+w2+w3Qij+w4vij212πσije-vij-uij2/2σij2dvij.Consequently, the mathematical formulation of VRPFSV can be expressed as follows:(4)min∑i=1n∑j=1ncEFijvijxij+f∑j=1nx0j(5)s.t.∑j=1nxij=1,∀i∈N/0(6)∑j=1nxij-∑j=1nxji=0,∀i∈N(7)∑j=1,j≠inQji-∑j=1,j≠inQij=Di,∀i∈N/0(8)Qij≤Qxij,∀i,j∈A(9)xij∈0,1,∀i,j∈A(10)Qij≥0,∀i,j∈A.In this optimization model, the two variable sets arexij and Qij. Equation (4) is the objective function which consists of the expected fuel consumption and the total vehicle cost. x0j=1 indicates a departure of a newly used vehicle from the depot; as any used vehicle will trigger a fixed cost, the sum of x0j 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 Q. Equation (9) indicates that xij is a binary variable. Equation (10) is the nonnegative constraint of vehicle load.


## 2.1. 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.(2)
The depot owns enough homogenous vehicles.(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.Then the notations used in the VRPFSV formulation are listed as follows:n: total number of customers.i: index of nodes (the depot is represented by 0).f: fixed cost of using a vehicle.Q: the capacity of a vehicle.N: set of nodes including the depot, i∈N.A: set of arcs formed by all pairs of nodes, (i,j)∈A,∀i∈N,j∈N,i≠j.lij: distance between node i and node j.vij: stochastic average travel speed on arc(i,j).Two decision variables are as follows:xij: binary variable indicating whether arc(i,j) is traveled. If arc(i,j) is traveled by a vehicle, then xij=1; otherwise, xij=0.Qij: load carried from node i to node j.


## 2.2. Formulation of the Model of VRPFSV

In this section, a mixed integer linear programming model is developed to construct the VRPFSV. The VRPFSV is an extension of the classical VRP and assumes that average travel speedvij on arc(i,j) is a nonnegative random variable. Furthermore, the average travel speed vij on arc(i,j) is assumed to obey normal distribution with mean uij and variance σij2.The fuel consumption model applied in this paper is based on the comprehensive model formed by Barth and Boriboonsomsin [33]. It is given by (1)Zijvij=λkNV+wγαijvij+γαijfijvij+βγvij3dijvij,where λ=ξ/κψ and γ=1/1000ηtfη are constants decided by different fuel properties, β=0.5CdρA, and w are vehicle characteristics related constants. αij is a constant associated with the road characteristics and acceleration, αij=τij+gsin⁡θij+gCrcos⁡θij, where τij is acceleration and θij 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(i,j) can be rewritten as (2)Zijvij=ω1vij+ω2+ω3fij+ω4vij2dij,where ω1=λkNV,ω2=λwγgCr,ω3=λγgCr,ω4=λβγ, and the price of the fuel c=1.4£/L. Equation (2) establishes a good correlation between fuel consumption and travel speed, weight, and distance.Table 2
Values and notations of parameters.
Notation Description Typical values W Curb-weight (kilogram) 6350 Ξ Fuel-to-air mass ratio 1 K Engineer friction factor (kilojoule/rev/liter) 0.2 N Engineer speed (meter/second2) 33 V Engineer displacement (liters) 5 G Gravitational constant (meter/second2) 9.81 C d Coefficient of aerodynamic drag 0.7 P Air density (kilogram/meter3) 1.204 A Frontal surface area (meter2) 3.912 C r Coefficient of rolling resistance 0.01 n t f Vehicle drive train efficiency 0.4 H Efficiency parameter for diesel engines 0.9 c Fuel and CO2 emissions cost per liter (£) 1.4 K Heating value of a typical diesel fuel (kilojoule/gram) 44 Ψ Conversion factor (gram/second to liter/second) 737 v m i n Lower speed limit (meter/second) 5 v m a x Upper speed limit (meter/second) 25 Source: Demir et al., 2012 [24].In our research, the service cost includes two components: the expected fuel consumption and the fixed vehicle cost.Consider arc(i,j); if it is traveled by a vehicle with average speed of vij, then it has an expected fuel consumption; according to (2), the expected fuel consumption is given in (3)EFijvij=∫vijminvijmaxlijw1vij+w2+w3Qij+w4vij212πσije-vij-uij2/2σij2dvij.Consequently, the mathematical formulation of VRPFSV can be expressed as follows:(4)min∑i=1n∑j=1ncEFijvijxij+f∑j=1nx0j(5)s.t.∑j=1nxij=1,∀i∈N/0(6)∑j=1nxij-∑j=1nxji=0,∀i∈N(7)∑j=1,j≠inQji-∑j=1,j≠inQij=Di,∀i∈N/0(8)Qij≤Qxij,∀i,j∈A(9)xij∈0,1,∀i,j∈A(10)Qij≥0,∀i,j∈A.In this optimization model, the two variable sets arexij and Qij. Equation (4) is the objective function which consists of the expected fuel consumption and the total vehicle cost. x0j=1 indicates a departure of a newly used vehicle from the depot; as any used vehicle will trigger a fixed cost, the sum of x0j 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 Q. Equation (9) indicates that xij is a binary variable. Equation (10) is the nonnegative constraint of vehicle load.


## 3. 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 itemxijQij 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(i,j), consider its expected fuel consumption:(11)EFijvij=xij∫vijminvijmaxlijω1vij+ω2+ω3Qij+ω4vij212πσije-vij-uij2/2σij2dvij=xij∫vijminvijmaxlijω1vij+ω2+ω4vij212πσije-vij-uij2/2σij2dvij+xij∫vijminvijmaxlijω3Qij12πσije-vij-uij2/2σij2dvij=xij∫vijminvijmaxlijω1vij+ω2+ω4vij212πσije-vij-uij2/2σij2dvij+Qijxij∫vijminvijmaxlijω312πσije-vij-uij2/2σij2dvij.In fact, the total fuel consumption can be reduced as follows:(12)EFijvij=xij∫vijminvijmaxlijw1vij+w2+w4vij212πσije-vij-uij2/2σij2dvij+Qij∫vijminvijmaxlijw312πσije-vij-uij2/2σij2dvij.Property 1.
In the VRPFSV, (11) is equivalent to (12), subjecting to constraint (8).Proof of Property1.
In order to prove this property, we only need to prove that the product ofQij and xij is equivalent to Qij. In fact, xij is a binary variable; if xij=1, then xijQij=Qij, else if xij=0, due to constraints (8) and (10), we get Qij=0. Consequently, xijQij=0=Qij. Thus, xijQij is equivalent to Qij which indicates that (11) is equivalent to (12).
According to (12), objective (4) is converted to the equivalent linear form:(13)∑i=1n∑j=1ncEFijvijxij+f∑j=1nx0j=∑i=0n∑j=0ncxij∫vijminvijmaxlijω1vij+ω2+ω4vij212πσije-vij-uij2/2σij2dvij+∑i=1n∑j=1ncQij∫vijminvijmaxlijω312πσije-vij-uij2/2σij2dvij+f∑j=1nx0j.
Let(14)C1=clij∫vijminvijmaxω1vij+ω2+ω4vij212πσije-vij-uij2/2σij2dvij(15)C2=clij∫vijminvijmaxω312πσije-vij-uij2/2σij2dvij.
As a result, objective (4) can be transformed into ∑i=1n∑j=1ncE[Fij(vij)]xij+f∑j=1nx0j=C1∑i=0n∑j=0nxij+C2∑i=1n∑j=1nQij+f∑j=1nx0j, which is a linear objective function about xij and Qij 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 calculateC1 and C2, 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  C1 and C2 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.


## 4. 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].


### 4.1. Construction of Initial Solution

In the string model, the unique depot is encoded as 0 and the customers are encoded as a series of positive numbers. The string vector which represents the initial route is formatted asS=0index⋯index0 and the length of S is n+2.Then makeS feasible for VRPFSV; that is, assign n 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 S is transformed to feasible solution with the following procedures. First, start from the second character of S, 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 S before the current character and then reset the cumulative demands to 0. Repeat the above steps until the end of S. Take n=10 as an example, assuming that S=0341759810620 is the initial string. As we do not know how many zeros are needed to be inserted into S 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 S; then we can obtain VRPS={0341075098100620000000}, which represents ten customers being visited by four vehicles.


### 4.2. Generate Neighborhood Solutions

Exchange rule: three commonly used exchange rules swap, relocation, and 2-opt (see Figure2) 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.Figure 2
Three exchange rules.
(a)
Swap (b)
Relocation (c)
2-opt


### 4.3. 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 arrayVRPSbest in the algorithm to record the best solution. Only when a better solution appears, will VRPSbest 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 AppendixB.This algorithm adopts adaptive cooling processing, and the temperature cooling coefficientτ can be set to A/A+loop, where A is the accepted solution numbers of the current temperature andloop 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 demandD and distance matrix L, vehicle capacity Q, fixed single vehicle cost f, speed distribution parameters of each arc u and σ2, the end temperature Tend, and the number of inner loopsloop.Step 2.
Generate initial stringS according to the method introduced in Section 4.1.Step 3.
Generate initial temperature. Randomly changeS 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 Step8. 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 Step9.Step 8.
Whether the end temperature is reached, if so, go to Step9. 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 AppendixA.In Step3, 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 Step5, let Ci and Cnow be the objective value of the ith iteration and the best objective of the accepted routes, and ΔC=Ci-Cnow. VRPSnow and VRPSnew denote, respectively, the minimum cost route among the accepted ones and the ith route. If ΔC<0, then replace VRPSnow with VRPSnew. Else, replace VRPSnow with VRPSnew if random(0,1)<Tnow-Tend/T0-Tend.The algorithm terminates on two conditions: the current temperature is below or equal to the end temperatureTend, 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 isOn.Proof.
In the outer temperature loop, the number of temperatures in the cooling procedure islogτT0/Tend; in the internal loop, the number of loops is loop, and each loop has a calculation of cost function C with time complexity O(n). Thus, the total time complexity becomes Onlooplogτ⁡T0/Tend. As loop, τ, T0, and Tend are all given constants, Theory 2 is proved.


## 4.1. Construction of Initial Solution

In the string model, the unique depot is encoded as 0 and the customers are encoded as a series of positive numbers. The string vector which represents the initial route is formatted asS=0index⋯index0 and the length of S is n+2.Then makeS feasible for VRPFSV; that is, assign n 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 S is transformed to feasible solution with the following procedures. First, start from the second character of S, 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 S before the current character and then reset the cumulative demands to 0. Repeat the above steps until the end of S. Take n=10 as an example, assuming that S=0341759810620 is the initial string. As we do not know how many zeros are needed to be inserted into S 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 S; then we can obtain VRPS={0341075098100620000000}, which represents ten customers being visited by four vehicles.


## 4.2. Generate Neighborhood Solutions

Exchange rule: three commonly used exchange rules swap, relocation, and 2-opt (see Figure2) 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.Figure 2
Three exchange rules.
(a)
Swap (b)
Relocation (c)
2-opt


## 4.3. 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 arrayVRPSbest in the algorithm to record the best solution. Only when a better solution appears, will VRPSbest 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 AppendixB.This algorithm adopts adaptive cooling processing, and the temperature cooling coefficientτ can be set to A/A+loop, where A is the accepted solution numbers of the current temperature andloop 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 demandD and distance matrix L, vehicle capacity Q, fixed single vehicle cost f, speed distribution parameters of each arc u and σ2, the end temperature Tend, and the number of inner loopsloop.Step 2.
Generate initial stringS according to the method introduced in Section 4.1.Step 3.
Generate initial temperature. Randomly changeS 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 Step8. 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 Step9.Step 8.
Whether the end temperature is reached, if so, go to Step9. 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 AppendixA.In Step3, 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 Step5, let Ci and Cnow be the objective value of the ith iteration and the best objective of the accepted routes, and ΔC=Ci-Cnow. VRPSnow and VRPSnew denote, respectively, the minimum cost route among the accepted ones and the ith route. If ΔC<0, then replace VRPSnow with VRPSnew. Else, replace VRPSnow with VRPSnew if random(0,1)<Tnow-Tend/T0-Tend.The algorithm terminates on two conditions: the current temperature is below or equal to the end temperatureTend, 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 isOn.Proof.
In the outer temperature loop, the number of temperatures in the cooling procedure islogτT0/Tend; in the internal loop, the number of loops is loop, and each loop has a calculation of cost function C with time complexity O(n). Thus, the total time complexity becomes Onlooplogτ⁡T0/Tend. As loop, τ, T0, and Tend are all given constants, Theory 2 is proved.


## 5. 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 different-scale 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 n 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 thatvij is the average travel speed on arc(i,j) and it obeys normal distribution with mean parameter uij and variance parameter σij.In these examples, speed limit of all arcs is 5 m/s~25 m/s. For each arc, the integer mean parametersuij are generated randomly from interval [5, 25] and set σij=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

Tables2 and 3 compared results of the optimal solutions found by CPLEX, SA, and ISA algorithm when n 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=ObjISAbest-Obji/Obji×100%,i=CEPLXorSAbest. Avg-CPU time is the average run time (CPU) of the SA or ISA algorithm.Table 3
Comparison results whenn=10.
Instance ID CPLEX Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from CPLEX Dev. from SA Obj (£) CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK10_01 99.63 1.39 99.63 99.63 33.01 99.96 99.63 13.61 0.00% 0.00% UK10_02 124.46 0.33 124.46 124.46 22.46 124.46 124.46 17.17 0.00% 0.00% UK10_03 100.81 0.47 100.97 100.97 32.81 101.35 100.97 5.36 0.15% 0.00% UK10_04 100.18 0.81 106.07 106.07 33.37 106.07 106.07 4.57 5.88% 0.00% UK10_05 107.19 1.22 109.69 109.69 42.67 109.69 109.69 22.07 2.33% 0.00%


#### 5.1.1. 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 ofloop and Tend. 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 loop or a smaller Tend may improve the solution quality a little but increase the computation time significantly.Figure 3
Sensitive analysis ofloop and Tend.By trading off the effectiveness and the efficiency of the algorithm, we found the values ofloop=100 and Tend=0.1 are appropriate parameters combination.As analyzed above, we setloop=100, Tend=0.1, c=1.4, and f=1.ω1, ω2, ω3, ω4 are set to the same value as mentioned in Section 2.2; that is, ω1=1.0176×10-3, ω2=5.3360×10-5, ω3=8.4032×10-9, and ω4=1.4122×10-7. In order to compare the results obtained by SA and ISA, we set the same parameter with the same value both in SA and in ISA. Thus, the parameters in SA is as follows: T0=1000,loop=100, and Tend=0.1.In Table3, 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.Table 4
Comparison results whenn=15.
Instance ID CPLEX Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from CPLEX Dev. from SA Obj (£) CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK15_01 161.86 18.42 159.59 159.59 33.01 160.16 159.59 13.51 1.40% 0.00% UK15_02 109.19 0.8 109.87 109.87 22.46 109.89 109.87 25.05 0.62% 0.00% UK15_03 167.96 46.66 167.96 167.96 32.81 170.22 167.96 6.43 0.00% 0.00% UK15_04 172.26 83.89 172.55 172.26 33.37 174.52 172.26 7.91 0.00% 0.00% UK15_05 166.73 0.88 166.73 166.73 42.67 166.95 166.73 7.25 0.00% 0.00%It is shown in Table5 that when n=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.Table 5
Comparison results whenn=25.
Instance ID CPLEX Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from CPLEX Dev. from SA Obj (£) CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK25_01 153.24 4.16 155.25 155.11 61.85 156 153.24 22.63 0.00% - 1.21 % UK25_02 185.74 41.39 189.7 187.6 52.25 189.6 185.74 25.05 0.00% - 1.01 % UK25_03 105.6 16.62 107.46 106.44 51.69 108.2 105.59 20.25 0.01% - 0.80 % UK25_04 134.63 9.73 138.38 136.26 51.90 137.6 136.59 31.36 1.46% 0.24 % UK25_05 187.95 24.97 189.68 189.15 53.63 192.2 187.95 76.98 0.00% - 0.63 %In Table6, 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.Table 6
Comparison results whenn=50.
Instance ID CPLEX Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from CPLEX Dev. from SA Obj (£) CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK50_01 314.75 55276.10 347.99 340.31 87.54 353.42 329.05 55.41 4.54% - 3.31 % UK50_02 317.82 30042.20 363.70 356.72 86.72 356.37 334.5 39.71 5.25% - 6.23 % UK50_03 325.77 35527.60 354.47 349.60 86.05 359.56 338.58 41.7 3.93% - 3.15 % UK50_04 418.24 55624.80 444.30 437.77 87.65 436.86 425.8 28.03 1.81% - 2.73 % UK50_05 352.95 46121.40 371.80 366.43 84.67 386.03 359.88 31.88 1.96% - 1.79 %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, Tables7 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.Table 7
Results ofn=100.
Instance ID Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from SA Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK100_01 786.44 774.13 155.56 765.4 718.92 99.39 - 7.13 % UK100_02 774.64 744.66 153.76 727.98 709.02 57.02 - 4.79 % UK100_03 690.51 685.67 154.38 682.02 647.96 93.33 - 5.50 % UK100_04 663.71 659.64 157.92 683.57 642.06 27.75 - 2.67 % UK100_05 629.12 617.56 156.52 612.15 571.8 92.49 - 7.41 %Table 8
Results ofn=200.
Instance ID Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from SA Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK200_01 1413.97 1387.87 327.66 1388.82 1307.69 40.31 - 6.13 % UK200_02 1372.61 1326.84 284.24 1303.79 1249.02 57.18 - 6.23 % UK200_03 1387.10 1363.50 371.36 1336.83 1264.7 48.27 - 7.81 % UK200_04 1316.76 1244.31 292.16 1295.71 1242.86 36.85 - 0.12 % UK200_05 1524.22 1503.30 293.40 1497.39 1381.58 42.39 - 8.81 %


#### 5.1.2. 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.Taken=200 as an example. First let c=1, and variable c values are tested over the ranges f∈0.01,0.1,10,100,1000. The total cost and number of vehicles are listed in Table 9.Table 9
Test results of different fixed vehicle cost.
Values off 1000 100 10 0.1 0.01 N 28 28 29 29 29 Total cost (£) 28975.39 4459.221 1474.35 1154.432 968.6605From Table9, 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 f and c 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.


### 5.2. 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 whenn=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.Table 10
Compared results ofn=10.
Given speed value (m/s) Objective value (£) Stochastic speed model Fixed speed model 5 92.42 201.39 7 158.4 160.77 9 141.34 140.16 11 129.78 129.1 13 124.03 123.57 15 124.46 121.7 17 122.74 122.48 19 125.56 125.35 21 130.13 129.95 23 132.95 136.04 25 71.21 143.47As can be seen from Table10, 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 (Figure4) 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.Figure 4
Relationship between fuel consumption and travel speed.


## 5.1. Quality and Efficiency of the ISA Algorithm

Tables2 and 3 compared results of the optimal solutions found by CPLEX, SA, and ISA algorithm when n 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=ObjISAbest-Obji/Obji×100%,i=CEPLXorSAbest. Avg-CPU time is the average run time (CPU) of the SA or ISA algorithm.Table 3
Comparison results whenn=10.
Instance ID CPLEX Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from CPLEX Dev. from SA Obj (£) CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK10_01 99.63 1.39 99.63 99.63 33.01 99.96 99.63 13.61 0.00% 0.00% UK10_02 124.46 0.33 124.46 124.46 22.46 124.46 124.46 17.17 0.00% 0.00% UK10_03 100.81 0.47 100.97 100.97 32.81 101.35 100.97 5.36 0.15% 0.00% UK10_04 100.18 0.81 106.07 106.07 33.37 106.07 106.07 4.57 5.88% 0.00% UK10_05 107.19 1.22 109.69 109.69 42.67 109.69 109.69 22.07 2.33% 0.00%


### 5.1.1. 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 ofloop and Tend. 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 loop or a smaller Tend may improve the solution quality a little but increase the computation time significantly.Figure 3
Sensitive analysis ofloop and Tend.By trading off the effectiveness and the efficiency of the algorithm, we found the values ofloop=100 and Tend=0.1 are appropriate parameters combination.As analyzed above, we setloop=100, Tend=0.1, c=1.4, and f=1.ω1, ω2, ω3, ω4 are set to the same value as mentioned in Section 2.2; that is, ω1=1.0176×10-3, ω2=5.3360×10-5, ω3=8.4032×10-9, and ω4=1.4122×10-7. In order to compare the results obtained by SA and ISA, we set the same parameter with the same value both in SA and in ISA. Thus, the parameters in SA is as follows: T0=1000,loop=100, and Tend=0.1.In Table3, 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.Table 4
Comparison results whenn=15.
Instance ID CPLEX Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from CPLEX Dev. from SA Obj (£) CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK15_01 161.86 18.42 159.59 159.59 33.01 160.16 159.59 13.51 1.40% 0.00% UK15_02 109.19 0.8 109.87 109.87 22.46 109.89 109.87 25.05 0.62% 0.00% UK15_03 167.96 46.66 167.96 167.96 32.81 170.22 167.96 6.43 0.00% 0.00% UK15_04 172.26 83.89 172.55 172.26 33.37 174.52 172.26 7.91 0.00% 0.00% UK15_05 166.73 0.88 166.73 166.73 42.67 166.95 166.73 7.25 0.00% 0.00%It is shown in Table5 that when n=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.Table 5
Comparison results whenn=25.
Instance ID CPLEX Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from CPLEX Dev. from SA Obj (£) CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK25_01 153.24 4.16 155.25 155.11 61.85 156 153.24 22.63 0.00% - 1.21 % UK25_02 185.74 41.39 189.7 187.6 52.25 189.6 185.74 25.05 0.00% - 1.01 % UK25_03 105.6 16.62 107.46 106.44 51.69 108.2 105.59 20.25 0.01% - 0.80 % UK25_04 134.63 9.73 138.38 136.26 51.90 137.6 136.59 31.36 1.46% 0.24 % UK25_05 187.95 24.97 189.68 189.15 53.63 192.2 187.95 76.98 0.00% - 0.63 %In Table6, 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.Table 6
Comparison results whenn=50.
Instance ID CPLEX Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from CPLEX Dev. from SA Obj (£) CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK50_01 314.75 55276.10 347.99 340.31 87.54 353.42 329.05 55.41 4.54% - 3.31 % UK50_02 317.82 30042.20 363.70 356.72 86.72 356.37 334.5 39.71 5.25% - 6.23 % UK50_03 325.77 35527.60 354.47 349.60 86.05 359.56 338.58 41.7 3.93% - 3.15 % UK50_04 418.24 55624.80 444.30 437.77 87.65 436.86 425.8 28.03 1.81% - 2.73 % UK50_05 352.95 46121.40 371.80 366.43 84.67 386.03 359.88 31.88 1.96% - 1.79 %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, Tables7 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.Table 7
Results ofn=100.
Instance ID Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from SA Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK100_01 786.44 774.13 155.56 765.4 718.92 99.39 - 7.13 % UK100_02 774.64 744.66 153.76 727.98 709.02 57.02 - 4.79 % UK100_03 690.51 685.67 154.38 682.02 647.96 93.33 - 5.50 % UK100_04 663.71 659.64 157.92 683.57 642.06 27.75 - 2.67 % UK100_05 629.12 617.56 156.52 612.15 571.8 92.49 - 7.41 %Table 8
Results ofn=200.
Instance ID Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from SA Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK200_01 1413.97 1387.87 327.66 1388.82 1307.69 40.31 - 6.13 % UK200_02 1372.61 1326.84 284.24 1303.79 1249.02 57.18 - 6.23 % UK200_03 1387.10 1363.50 371.36 1336.83 1264.7 48.27 - 7.81 % UK200_04 1316.76 1244.31 292.16 1295.71 1242.86 36.85 - 0.12 % UK200_05 1524.22 1503.30 293.40 1497.39 1381.58 42.39 - 8.81 %


### 5.1.2. 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.Taken=200 as an example. First let c=1, and variable c values are tested over the ranges f∈0.01,0.1,10,100,1000. The total cost and number of vehicles are listed in Table 9.Table 9
Test results of different fixed vehicle cost.
Values off 1000 100 10 0.1 0.01 N 28 28 29 29 29 Total cost (£) 28975.39 4459.221 1474.35 1154.432 968.6605From Table9, 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 f and c 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.


## 5.1.1. 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 ofloop and Tend. 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 loop or a smaller Tend may improve the solution quality a little but increase the computation time significantly.Figure 3
Sensitive analysis ofloop and Tend.By trading off the effectiveness and the efficiency of the algorithm, we found the values ofloop=100 and Tend=0.1 are appropriate parameters combination.As analyzed above, we setloop=100, Tend=0.1, c=1.4, and f=1.ω1, ω2, ω3, ω4 are set to the same value as mentioned in Section 2.2; that is, ω1=1.0176×10-3, ω2=5.3360×10-5, ω3=8.4032×10-9, and ω4=1.4122×10-7. In order to compare the results obtained by SA and ISA, we set the same parameter with the same value both in SA and in ISA. Thus, the parameters in SA is as follows: T0=1000,loop=100, and Tend=0.1.In Table3, 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.Table 4
Comparison results whenn=15.
Instance ID CPLEX Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from CPLEX Dev. from SA Obj (£) CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK15_01 161.86 18.42 159.59 159.59 33.01 160.16 159.59 13.51 1.40% 0.00% UK15_02 109.19 0.8 109.87 109.87 22.46 109.89 109.87 25.05 0.62% 0.00% UK15_03 167.96 46.66 167.96 167.96 32.81 170.22 167.96 6.43 0.00% 0.00% UK15_04 172.26 83.89 172.55 172.26 33.37 174.52 172.26 7.91 0.00% 0.00% UK15_05 166.73 0.88 166.73 166.73 42.67 166.95 166.73 7.25 0.00% 0.00%It is shown in Table5 that when n=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.Table 5
Comparison results whenn=25.
Instance ID CPLEX Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from CPLEX Dev. from SA Obj (£) CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK25_01 153.24 4.16 155.25 155.11 61.85 156 153.24 22.63 0.00% - 1.21 % UK25_02 185.74 41.39 189.7 187.6 52.25 189.6 185.74 25.05 0.00% - 1.01 % UK25_03 105.6 16.62 107.46 106.44 51.69 108.2 105.59 20.25 0.01% - 0.80 % UK25_04 134.63 9.73 138.38 136.26 51.90 137.6 136.59 31.36 1.46% 0.24 % UK25_05 187.95 24.97 189.68 189.15 53.63 192.2 187.95 76.98 0.00% - 0.63 %In Table6, 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.Table 6
Comparison results whenn=50.
Instance ID CPLEX Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from CPLEX Dev. from SA Obj (£) CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK50_01 314.75 55276.10 347.99 340.31 87.54 353.42 329.05 55.41 4.54% - 3.31 % UK50_02 317.82 30042.20 363.70 356.72 86.72 356.37 334.5 39.71 5.25% - 6.23 % UK50_03 325.77 35527.60 354.47 349.60 86.05 359.56 338.58 41.7 3.93% - 3.15 % UK50_04 418.24 55624.80 444.30 437.77 87.65 436.86 425.8 28.03 1.81% - 2.73 % UK50_05 352.95 46121.40 371.80 366.43 84.67 386.03 359.88 31.88 1.96% - 1.79 %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, Tables7 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.Table 7
Results ofn=100.
Instance ID Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from SA Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK100_01 786.44 774.13 155.56 765.4 718.92 99.39 - 7.13 % UK100_02 774.64 744.66 153.76 727.98 709.02 57.02 - 4.79 % UK100_03 690.51 685.67 154.38 682.02 647.96 93.33 - 5.50 % UK100_04 663.71 659.64 157.92 683.57 642.06 27.75 - 2.67 % UK100_05 629.12 617.56 156.52 612.15 571.8 92.49 - 7.41 %Table 8
Results ofn=200.
Instance ID Solutions found by SA in 10 runs Solutions found by ISA in 10 runs Dev. from SA Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) Avg-Obj (£) Best-Obj (£) Avg-CPU time (s) UK200_01 1413.97 1387.87 327.66 1388.82 1307.69 40.31 - 6.13 % UK200_02 1372.61 1326.84 284.24 1303.79 1249.02 57.18 - 6.23 % UK200_03 1387.10 1363.50 371.36 1336.83 1264.7 48.27 - 7.81 % UK200_04 1316.76 1244.31 292.16 1295.71 1242.86 36.85 - 0.12 % UK200_05 1524.22 1503.30 293.40 1497.39 1381.58 42.39 - 8.81 %


## 5.1.2. 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.Taken=200 as an example. First let c=1, and variable c values are tested over the ranges f∈0.01,0.1,10,100,1000. The total cost and number of vehicles are listed in Table 9.Table 9
Test results of different fixed vehicle cost.
Values off 1000 100 10 0.1 0.01 N 28 28 29 29 29 Total cost (£) 28975.39 4459.221 1474.35 1154.432 968.6605From Table9, 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 f and c 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.


## 5.2. 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 whenn=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.Table 10
Compared results ofn=10.
Given speed value (m/s) Objective value (£) Stochastic speed model Fixed speed model 5 92.42 201.39 7 158.4 160.77 9 141.34 140.16 11 129.78 129.1 13 124.03 123.57 15 124.46 121.7 17 122.74 122.48 19 125.56 125.35 21 130.13 129.95 23 132.95 136.04 25 71.21 143.47As can be seen from Table10, 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 (Figure4) 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.Figure 4
Relationship between fuel consumption and travel speed.


## 6. 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 large-scale 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.



---

*Source: 6329203-2017-01-29.xml*