# Spatiotemporal-Dependent Vehicle Routing Problem Considering Carbon Emissions

Authors: Ziqi Liu; Yeping Chen; Jian Li; Dongqing Zhang
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

Aiming to improve the timeliness of logistics distribution and render the optimized route scheme effective under the real traffic network, we study the green vehicle routing problem with dynamic travel speed from both dimensions of time and space. A discrete formulation is proposed to calculate the travel time based on periods and arcs, which allows a vehicle to travel across an arc in multiple periods. Then, we establish a mixed-integer nonlinear programming model with minimum distribution costs including transportation costs, carbon emissions costs, and penalty costs on earliness and tardiness. A hybrid adaptive genetic algorithm with elite neighborhood search is developed to solve the problem. In the algorithm, a neighborhood search operator is employed to optimize elite individuals so that the algorithm can stimulate the intensification and avoid falling into a local optimum. Experimental instances are constructed based on benchmark instances of vehicle routing problem. The numerical results indicate that the proposed algorithm is rather effective in global convergence. Compared with the routing schemes in which travel speed merely varies with time periods or locations, the vehicle route optimized on spatiotemporal-varying speed outperforms them in terms of carbon emissions and timeliness. The research can provide a scientific and reasonable method for logistics enterprises to plan the vehicle schedule focusing on spatiotemporal-dependent speed of the road network.


---

## Body

## 1. Introduction

Most freight companies ignore the varying speed of the vehicle when planning the distribution routes (i.e., Çam and Sezen [1]), so that the preoptimized route scheme is hard for vehicles to service customers within the time windows. In order to ensure the timeliness of distribution, the vehicle route planning should apply the traffic information about the road network to calculate the travel time more precisely. At the same time, with the increasing carbon emissions, reducing vehicle fuel consumption has become the main trend of logistics distribution. Thus, the objective of distribution should additionally take carbon emissions influenced by travel speed on arcs into account. Therefore, we focus on the green vehicle routing problem with the spatiotemporal-varying vehicle speed and soft time windows, which is a variant of vehicle routing problem (VRP).In classical VRP, a fleet of vehicles at one depot must deliver known demands to a set of customers [2]. The objective is to minimize the travel cost with a set of specified routes for the vehicles under various constraints. VRP with time windows (VRPTW) can be generalized by imposing additional constraints that customers are only available during a specified time interval, called time window [3]. In most VRPTW, it is assumed that the distribution cost is proportional to the travel distance; however, the cost of carbon emission is not proportional to the travel distance. Therefore, considering the minimum distance or the minimum distribution cost cannot reduce carbon emissions. The cost of carbon emission is directly related to vehicle fuel consumption, which is associated with other factors, such as vehicle speed, current vehicle load, and road inclination [4]. With the heavy traffic congestion, vehicle speed has become an important factor that affects carbon emissions. Travel speed varies with different moments and locations, and such dynamic travel speed affects carbon emission and the arrival time at the customer [5]. Therefore, this study investigates the spatiotemporal-dependent green vehicle routing problem with time windows. The fixed costs of distribution, the variable cost of transportation, penalty costs on time windows, and carbon emission cost are considered in the objective function.We focus on the spatiotemporal-dependent green vehicle routing problem with time windows which is related to the time-dependent vehicle routing problem with time windows (TDVRPTW) and green vehicle routing problem with time windows (GVRPTW). VRP has been proven to be an NP-hard problem, so it is difficult to obtain the optimal solution in the finite amount of computation time. Liu et al. [6] applied ant colony algorithm, a metaheuristic algorithm, to solve time-dependent vehicle routing problem with time windows (TDVRPTW). Time-dependent green vehicle routing problem with time windows (TDGVRPTW) studied in our work is more complicated than TDVRPTW which needs to calculate the carbon emissions of each period and road section. Thus, we apply the genetic algorithm, a classical metaheuristic algorithm to find the near-optimal solutions of the problem. The contributions of this work are summarized as follows:(1)
Due to the real traffic network speed distribution, the calculation of travel time is a discrete function of periods and arcs, which meets the first-in and first-out rules. We adopt the fuel consumption introduced by Esteves-Booth et al. [7] with the spatiotemporal speed considered. This work establishes a mixed-integer nonlinear planning model (MINP) with the total distribution costs.(2)
A hybrid adaptive genetic algorithm with elite neighborhood search (HAGA_ELS) is proposed. In the algorithm, the elite neighborhood search is employed to avoid premature convergence and balance the intensification and diversification of the solution.The rest of the article is as follows. Section2 introduces relevant research. Section 3 describes the research problem and introduces the calculation method of travel time with spatiotemporal dependence, upon which the mixed-integer nonlinear planning model is formulated. Section 4 presents HAGA_ELS. Section 5 assesses the algorithm through modified benchmark instances and analyzes the sensitivity of the model to different parameter settings. Section 6 presents the conclusions.


## 2. Literature Review

VRPTW was developed based on classic VRP. Russell [8] was the first to propose an effective heuristic algorithm to solve the problem of travelers with time windows. Solomon [9] established a mixed-integer nonlinear planning problem for VRPTW. Most studies about VRPTW used a fixed speed, while the dynamics of speed in real traffic networks was ignored. At the same time, to reduce the increasing carbon emissions, the objective of distribution should additionally take carbon emissions influenced by travel speed into account. Therefore, on the basis of VRPTW, this work mainly reviews the time-dependent vehicle routing problem with time windows (TDVRPTW) and time-dependent green vehicle routing problem with time windows (TDGVRPTW).


### 2.1. Research on TDVRPTW

The traditional VRP assumes that the speed is constant, but this assumption is inconsistent with the real situation because the vehicle is running in the real traffic network [3]. Thus, time-dependent VRP (TDVRP) has arisen to describe more real network optimization problems. TDVRP accounts for the real traffic conditions by setting travel time based on the path between two customers and the time of day. Generally, time-dependent issues can be addressed in two ways: one is that travel time changes over time, and the other is that travel speed varies with moment.The first type of research on time-dependent travel time directly assumes that travel time is a function related to the time in a day. The travel time function can be piecewise linear form (Sun et al. [10]) or distribute function form (Hu et al. [11] and Babaei and Rajabi-Bahaabadi [12]). Sun et al. [10] applied a travel time function that is piecewise linear and continuous based on the change of speed. As for the distribution function of travel time, Hu et al. [11] built a travel time distribution function and considered the uncertainty of travel time by adding uncertainty budget into the function. Babaei and Rajabi-Bahaabadi [12] assumed that travel time for each time of day satisfies a distribution function and the travel time between specific customers can be calculated by the departure time. The travel time can also be acquired from intelligent transportation system. Saint-Guillain et al. [13] retrieved the travel time between customers by Google Map’s API before scheduling.The second research on time-dependent travel speed studies the relationship between travel speed and time to simulate the real traffic network. Travel speed is mostly a stepwise function (Xiao and Konak [14], Wu and Ma [15], Zhu and Hu [16], Duan et al. [17], Gruler et al. [18], and Gmira et al. [19]), or a continuous function (Xu et al. [20] and Fan et al. [21]), or retrieved from historical data (Franceschetti et al. [22]). Gmira et al. [19] assumed travel speed of each time and arc as a stepwise function, from which a piecewise travel time function can be derived. As for the continuous speed function, Xu et al. [20] used a sine function to simulate time-dependent speed and a nondominated sorting genetic algorithm to solve the problem. Fan et al. [21] also used the sine function in [20] to describe the travel speed but segmented it based on the road grade. Some research studies directly retrieved travel speed from historical data rather than build a function. Franceschetti et al. [22] assumed that the travel speed and time period are constants that can be extracted from historical data, and then travel time can be calculated from the speed and distance.Most of the above literature studies only consider the speed varying in time, ignoring the speed varying in space. Table1 lists the methods of time-dependent issue processing of the above literature studies.Table 1
Literature studies in TDVRPTW.
InstancesYearTravel timeTravel speedPiecewise linear functionDistribution functionHistorical dataStepwise functionContinuous functionHistorical dataSun et al. [10]2018√Hu et al. [11]2018√Babaei and Rajabi-Bahaabadi [12]2019√Saint-Guillain et al. [13]2021√Xiao and Konak [14]2016√Wu and Ma [15]2017√Zhu and Hu [16]2019√Duan et al. [17]2019√Gruler et al. [18]2020√Gmira et al. [19]2021√Xu et al. [20]2019√Fan et al. [21]2021√Franceschetti et al. [22]2013√Our workTemporal axis: travel speed varies with a stepwise function converted from continuous functionSpatial axis: travel speed obeys a distribution in accordance with historical data


### 2.2. Research on TDGVRPTW

In recent years, reducing carbon emissions has become an important issue in vehicle routing optimization. Kwon et al. [23] aimed at the green vehicle routing problem to set carbon emission as a ratio function for driving distance, thus equating the problem of reducing carbon emission to minimizing driving paths, but which ignored the impact of speed on fuel consumption.In the subsequent research studies, time-varying speed was taken into account in TDGVRPTW. Some of them only targeted carbon emissions. For example, Qian and Eglese [24] established a model that aims to minimize fuel consumption under time-dependent travel speed and used the tabu search algorithm to solve the problem. And Xiao and Konak [25] formulated a model for a fleet of heterogeneous vehicles operating, in which carbon emissions are calculated by dynamic payload weights and time-varying speed. While some research studies considered the biobjective of minimizing carbon emissions and costs (see Wang et al. [26] and Poonthalir and Nadarajan [27]), in which Poonthalir and Nadarajan [27] adopted triangular distribution to describe varying speed.Further studies focused on more factors affecting carbon emissions. Xu et al. [20] calculated the fuel consumption by the load of a vehicle and its capacity, time-varying speed, driving distance, and traffic congestion. Fan et al. [21] established a travel speed function that considers three road conditions and analyzed the impact of vehicle type, speed, load, and road gradient on fuel consumption, proposing a hybrid variable neighborhood genetic algorithm for solving the problem.Table2 summarizes and compares the formulation in the stream of TDGVRPTW.Table 2
Literature studies in TDGVRPTW.
ReferencesYearProblem featuresObjectivesSolutionTemporal speedSpatial speedFactors in fuel consumption/carbon emissionKwon et al. [23]2013DistanceOperation costs and carbon emission trading costsSimulated annealing algorithmQian and Eglese [24]2016∗SpeedCarbon emissionsTabu search algorithmXiao and Konak [25]2015∗Speed, load, and vehicle typeFuel consumptionPartial mixed-integer planning and iterated neighborhood searchWang et al. [26]2019∗Speed and distanceCarbon emission and operating costsA multiphase hybrid heuristic algorithmPoonthalir and Nadarajan [27]2018∗Speed and distanceCarbon emission and operating costsParticle swarm optimizationXu et al. [20]2019∗Speed, load, vehicle capacity, and distanceFuel consumption and customer satisfactionNondominated sorting genetic algorithm IIFan et al. [21]2021∗∗ (Divided by three road grades)Speed, load, vehicle type, and road gradientTotal costs: fixed costs, penalty costs, and fuel costsHybrid genetic algorithmOur work∗∗ (Stochastic speed distribution)Speed, load, vehicle parameters, and distanceTotal costs: transportation costs, carbon emission costs, and penalty costsHAGA_ELSFrom the literature presented above, we can conclude that most of the studies on TDGVRPTW did not consider the difference in speed or travel time on different arcs. Among them, Fan et al. [21] considered the time dependence of travel speed under different road grades and established travel speed function of three road grades. However, the stochastic speed distribution was not fully considered, and the spatial difference in speed was only reflected in the road classification. Given that a spatial speed difference also exists in the same road grade section, the treatment cannot fully reflect the spatial difference in real speed in the wide distribution problem.In the current study, a spatiotemporal-dependent speed division pattern (T&S) is used. In this pattern, each arc has different speeds in different periods, and the speed of arcs is determined by periods and arcs. The travel speed in different dimensions obeys a certain distribution. This study analyzes the limitations of the time-dependent speed division pattern (T) and the space-dependent speed division pattern (S).


## 2.1. Research on TDVRPTW

The traditional VRP assumes that the speed is constant, but this assumption is inconsistent with the real situation because the vehicle is running in the real traffic network [3]. Thus, time-dependent VRP (TDVRP) has arisen to describe more real network optimization problems. TDVRP accounts for the real traffic conditions by setting travel time based on the path between two customers and the time of day. Generally, time-dependent issues can be addressed in two ways: one is that travel time changes over time, and the other is that travel speed varies with moment.The first type of research on time-dependent travel time directly assumes that travel time is a function related to the time in a day. The travel time function can be piecewise linear form (Sun et al. [10]) or distribute function form (Hu et al. [11] and Babaei and Rajabi-Bahaabadi [12]). Sun et al. [10] applied a travel time function that is piecewise linear and continuous based on the change of speed. As for the distribution function of travel time, Hu et al. [11] built a travel time distribution function and considered the uncertainty of travel time by adding uncertainty budget into the function. Babaei and Rajabi-Bahaabadi [12] assumed that travel time for each time of day satisfies a distribution function and the travel time between specific customers can be calculated by the departure time. The travel time can also be acquired from intelligent transportation system. Saint-Guillain et al. [13] retrieved the travel time between customers by Google Map’s API before scheduling.The second research on time-dependent travel speed studies the relationship between travel speed and time to simulate the real traffic network. Travel speed is mostly a stepwise function (Xiao and Konak [14], Wu and Ma [15], Zhu and Hu [16], Duan et al. [17], Gruler et al. [18], and Gmira et al. [19]), or a continuous function (Xu et al. [20] and Fan et al. [21]), or retrieved from historical data (Franceschetti et al. [22]). Gmira et al. [19] assumed travel speed of each time and arc as a stepwise function, from which a piecewise travel time function can be derived. As for the continuous speed function, Xu et al. [20] used a sine function to simulate time-dependent speed and a nondominated sorting genetic algorithm to solve the problem. Fan et al. [21] also used the sine function in [20] to describe the travel speed but segmented it based on the road grade. Some research studies directly retrieved travel speed from historical data rather than build a function. Franceschetti et al. [22] assumed that the travel speed and time period are constants that can be extracted from historical data, and then travel time can be calculated from the speed and distance.Most of the above literature studies only consider the speed varying in time, ignoring the speed varying in space. Table1 lists the methods of time-dependent issue processing of the above literature studies.Table 1
Literature studies in TDVRPTW.
InstancesYearTravel timeTravel speedPiecewise linear functionDistribution functionHistorical dataStepwise functionContinuous functionHistorical dataSun et al. [10]2018√Hu et al. [11]2018√Babaei and Rajabi-Bahaabadi [12]2019√Saint-Guillain et al. [13]2021√Xiao and Konak [14]2016√Wu and Ma [15]2017√Zhu and Hu [16]2019√Duan et al. [17]2019√Gruler et al. [18]2020√Gmira et al. [19]2021√Xu et al. [20]2019√Fan et al. [21]2021√Franceschetti et al. [22]2013√Our workTemporal axis: travel speed varies with a stepwise function converted from continuous functionSpatial axis: travel speed obeys a distribution in accordance with historical data


## 2.2. Research on TDGVRPTW

In recent years, reducing carbon emissions has become an important issue in vehicle routing optimization. Kwon et al. [23] aimed at the green vehicle routing problem to set carbon emission as a ratio function for driving distance, thus equating the problem of reducing carbon emission to minimizing driving paths, but which ignored the impact of speed on fuel consumption.In the subsequent research studies, time-varying speed was taken into account in TDGVRPTW. Some of them only targeted carbon emissions. For example, Qian and Eglese [24] established a model that aims to minimize fuel consumption under time-dependent travel speed and used the tabu search algorithm to solve the problem. And Xiao and Konak [25] formulated a model for a fleet of heterogeneous vehicles operating, in which carbon emissions are calculated by dynamic payload weights and time-varying speed. While some research studies considered the biobjective of minimizing carbon emissions and costs (see Wang et al. [26] and Poonthalir and Nadarajan [27]), in which Poonthalir and Nadarajan [27] adopted triangular distribution to describe varying speed.Further studies focused on more factors affecting carbon emissions. Xu et al. [20] calculated the fuel consumption by the load of a vehicle and its capacity, time-varying speed, driving distance, and traffic congestion. Fan et al. [21] established a travel speed function that considers three road conditions and analyzed the impact of vehicle type, speed, load, and road gradient on fuel consumption, proposing a hybrid variable neighborhood genetic algorithm for solving the problem.Table2 summarizes and compares the formulation in the stream of TDGVRPTW.Table 2
Literature studies in TDGVRPTW.
ReferencesYearProblem featuresObjectivesSolutionTemporal speedSpatial speedFactors in fuel consumption/carbon emissionKwon et al. [23]2013DistanceOperation costs and carbon emission trading costsSimulated annealing algorithmQian and Eglese [24]2016∗SpeedCarbon emissionsTabu search algorithmXiao and Konak [25]2015∗Speed, load, and vehicle typeFuel consumptionPartial mixed-integer planning and iterated neighborhood searchWang et al. [26]2019∗Speed and distanceCarbon emission and operating costsA multiphase hybrid heuristic algorithmPoonthalir and Nadarajan [27]2018∗Speed and distanceCarbon emission and operating costsParticle swarm optimizationXu et al. [20]2019∗Speed, load, vehicle capacity, and distanceFuel consumption and customer satisfactionNondominated sorting genetic algorithm IIFan et al. [21]2021∗∗ (Divided by three road grades)Speed, load, vehicle type, and road gradientTotal costs: fixed costs, penalty costs, and fuel costsHybrid genetic algorithmOur work∗∗ (Stochastic speed distribution)Speed, load, vehicle parameters, and distanceTotal costs: transportation costs, carbon emission costs, and penalty costsHAGA_ELSFrom the literature presented above, we can conclude that most of the studies on TDGVRPTW did not consider the difference in speed or travel time on different arcs. Among them, Fan et al. [21] considered the time dependence of travel speed under different road grades and established travel speed function of three road grades. However, the stochastic speed distribution was not fully considered, and the spatial difference in speed was only reflected in the road classification. Given that a spatial speed difference also exists in the same road grade section, the treatment cannot fully reflect the spatial difference in real speed in the wide distribution problem.In the current study, a spatiotemporal-dependent speed division pattern (T&S) is used. In this pattern, each arc has different speeds in different periods, and the speed of arcs is determined by periods and arcs. The travel speed in different dimensions obeys a certain distribution. This study analyzes the limitations of the time-dependent speed division pattern (T) and the space-dependent speed division pattern (S).


## 3. Problem and Mathematical Model

TDGVRPTW studied in this paper can be described as follows. In the network graphG=V,E, the set of nodes is V=0∪​V0, where 0 is the distribution center and V0=1,...,n is the set of known customers. The set of edges is E=i,j|∨i≠j∈V. The location xi,yi, nonnegative demand wi, service time STi, and delivery time of the arrival window ETi,LTi of customer i are known. The shortest path between two nodes, namely, edge i,j∈E, has a fixed distance of dij. The distribution road network has a time-varying traffic state, and the distribution time range is divided into a fixed number of time intervals L=1,2,…,h. The different arcs i,j∈E in each period l∈L have various travel speed vijl, and travel speed vijl is known in advance in accordance with the travel speed function. Distribution center has a certain number of distribution vehicles K=1,2,…,m, and all customers V0 go back to the distribution center after the completion of the distribution. Distribution vehicle k∈K has the same maximum load Q. The fuel consumption formula FUELv,q,d is given based on travel speed v, load weight q, and travel distance d, and fuel-to-CO2 coefficient ρ is known. On the premise of the constraints of vehicle loading and customer time window, the vehicle distribution routing is optimized considering the carbon emission of the vehicle. Thus, the total costs including the fixed costs of distribution, transportation costs, carbon emission costs, and penalty costs of the time window are the lowest. In the TDGVRPTW, the following constraints need to be met:(1)
The distribution center is unique, and all vehicles from the distribution center go back to the distribution center after completing the distribution task.(2)
Each customer can only be served once and by one vehicle only.(3)
Only one type of vehicle is available in the distribution center, and the number is unlimited.(4)
Customer demands for goods, location, and service time window are known in advance.(5)
Vehicles in service time do not produce fuel consumption.


### 3.1. Spatiotemporal-Dependent Travel Time Function

In a time-varying road network, the vehicle speed changes continuously with time. According to the historical traffic data, the morning rush hour is usually from 7: 00 am to 9: 00 am, and the evening rush hour is from 17: 00 pm to 19: 00 pm. Many vehicles emerge during the rush hour, and the average speed on the road decreases. For example, the relationship between vehicle speedv and time t in an entire day proposed by Xu et al. [20] can be approximated by the trigonometric function in equation (1), where φp, γp, and δp are constants related to the road conditions of the day:(1)vt=φ1sinγ1t+δ1,t∈0,T1⋮⋮φpsinγpt+δp,t∈Tp−1,Tp.TDVRPTW is an NP-hard problem, and it greatly increases the difficulty of solving the problem by using the function of continuous varying speed. To reduce the computational complexity caused by nonlinear speed change, this study transforms the trigonometric function of travel speed into a piecewise function in accordance with the period and takes the average value of velocity in each period as the constant velocity in that period (Figure1). The finer the time division is, the closer the speed segmentation function is to the continuous change in speed. The time interval is divided into periods L=1,2,…,h, the time of each period is recorded as Γll∈L, and the time interval is ΔI. The travel speed of the vehicle in any period l is(2)vl=∫ll+ΔIvtdtΔI,l∈L.Figure 1
Travel speed variation function.In some studies (e.g., Kwon et al. [23]), the travel time calculation of arc i,jdepends on the departure time of vehicles from node i and it assumes that the characteristics of arcs are the same. However, in a real road network, the traffic congestion of different arcs varies due to the influence of geographic location, traffic flow, and other factors. The all-day travel speed of each arc is not completely equal, instead the travel speed is spatiotemporal-dependent; that is, at the same time l∈L, the traveling speed vijl of different arc i,j is various. Therefore, the travel speed and travel time in the process of vehicle distribution not only depend on the time when the vehicle leaves node i but also relate to the arc from node i to the next node j. According to the historical traffic data on all-day arcs in a province of China in the study of Yang [28], the travel speed of different arcs in the same period obeys a normal distribution; that is, the travel speed vijl of any arc i,j∈E in period l obeys normal distribution Nvl,σ2, and the variance σ2 reflects the difference degree of travel speeds of different arcs. On this basis, the nonnegative random velocity matrix V˜ of all arcs in the road network in different periods can be obtained. A small network with 6 nodes, 21 arcs, and 9 periods is taken as an example, as shown in Figure 2. Figure 3 presents the mean value of the speed distribution of the small network, including 189 (21 × 9) random velocity variables.Figure 2
Small-size network.Figure 3
Velocity matrix.In the distribution process, vehicles are likely to cross from a period to the next period or even cross multiple periods when passing through the arci,j. Hill and Benton [29] and others did not consider the FIFO criterion when dealing with time-varying speeds across periods; that is, the vehicles that start first should arrive first starting from one point to the same next point and driving on the same arc [30]. Several studies (e.g., Ioannou et al. [31]) ensured the continuity of travel time by adjusting the speed across periods, but the calculation was complicated for solving problems. Therefore, this study refers to the work of Li et al. [32] to improve the calculation of travel time across periods, which satisfies the FIFO criterion. The vehicle is assumed to travel from node i to node j at time ai across S periods, so road section length dij is divided into d1,d2,…,dS+1. The speed in the segment ss=1,2,…,S+1 is vs∈V˜, and the time of division is Γs, as shown in Figure 4.Figure 4
Example of a time path without the nonstopping assumption.The travel time of road sectiond1 is Γ1−ai, and the travel time of subsequent road section ds is Γs−Γs−1s=2,3,…,S. Thus, an equation between section length and time is available:(3)d1=v1Γ1−ai,ds=vsΓs−Γs−1s=2,3,…,S.The driving distance of theS + 1 section is(4)dS+1=dij−∑s=1Sds=dij−v1Γ1−ai−∑s=2SvsΓs−Γs−1.Therefore, the arrival timetj of node j is(5)tj=ΓS+dS+1vS+1=v1vS+1ai+ΓS+dij−v1Γ1−∑s=2SvsΓs−Γs−1vS+1.According to equation (5), the time tj when the vehicle arrives at node j is an increasing function of time ai when the vehicle departs from node i, which meets the FIFO criterion.When the arci,j spans S periods, the following conditions should be met when the vehicle starts from node i: the vehicle starts from the first period and has not finished driving in arc i,j from the first period to the Sth period. The corresponding expression is as follows:(a)
ai<Γ1,(b)
v1Γ1−ai+∑s=2SvsΓs−Γs−1≤dij.The abovementioned arrangement yields(6)Γ1−dijv1+∑s=2Svsv1Γs−Γs−1≤ai<Γ1.The calculation formula of the arrival timetj of the next node j where any arc i,j crosses SS<h periods is as follows:(7)tj=ti+dijvijl,ai<Γl+1−dijvijl,l=1,…,L,vijlvijl+Sai+Γl+S+dij−vijlΓl+1−∑s=l+1l+S−1vijsΓs−Γs−1vijl+S,Γl+1−dijvijl+∑s=l+1l+S−1vijl+SvijlΓs−Γs−1≤ai<Γl+1,l=1,…,L−S.


### 3.2. Calculation of the Amount of Fuel Consumption

Carbon emissions are directly affected by vehicle fuel consumption. Hence, this study applies vehicle fuel consumption to calculate carbon emissions. Considering that most distribution vehicles are heavy types, a carbon emission model is established in combination with the comprehensive modal emission model (CMEM) [7]. The objective function of fuel consumption when driving in different sections is as follows:(8)fuel=λα1dv2+α2dv+α3Qd.In equation (8), λ=ζ/κφ, ζ is the fuel-air mass ratio, κ is a typical diesel calorific value, φ is the fuel conversion coefficient, engine module coefficient α1=0.5CdAfρ/1000εω, speed module coefficient α2=kNV, and vehicle load module coefficient α3=ρ0+ρmax−ρ0/Q0. The meanings and values of specific parameters are presented in Appendix A.In summary, throughout the distribution, the total fuel consumption after passing through all customer nodes is(9)FUEL=∑i,j∈E∑k∈Kxijkfuelij=∑i,j∈E∑k∈Kxijk∑l∈Lλα1dijlvijl2+α2dijlvijl+α3Qijdijl.In equation (9), xijk is a binary variable indicating whether arc i,j is traveled by vehicle k (xijk=1) or not (xijk=0). fuelij is the fuel consumption of the vehicle from node i to node j. dijl is the length of the arc i,j during time l. vijl represents the travel speed of the arc i,j during time l. And Qij represents the loading weight of the vehicle from node i to node j.Figure5 shows the relationship between fuel consumption and travel speed for a fixed vehicle type. With the increase of travel speed, the fuel consumption decreases at first. However, when the travel speed reaches a certain value, here 48.14 km/h, the fuel consumption begins to increase.Figure 5
Fuel consumption function of travel speed (with 2000 kg of load).The specific process of calculating carbon emission in consideration of spatiotemporal-dependent travel speed is shown in Algorithm1.Algorithm1: Calculation of travel time and fuel consumption.
//Depart fromi to j, assume that the period is divided into L segments(1)
Inputdij,ai,Qij, Velocity matrix V˜, Time slot division Γ, parametersλ,α1,α2, α3;(2)
Γ0←ai,v←V˜arci,j,whereΓ>ai,fuel←0,d′←dij;//Judging whether to span S periods(3)
If the criteria of crossing S periods is satisfied. //i.e. The path is divided into S + 1 segments(4)
Fork = 1 to S do begin //Recursive method(5)
dk←vkΓk−Γk−1; //The distance traveled in each period(6)
fuel←fuel + λα1dkvk2+α2Γk−Γk−1+α3Qijdk;(7)
d′←d′−dk;(8)
End for(9)
dS+1←d′; //the length of the S + 1 section//Calculate the arrival time at nodej and fuel consumption of traveling arci,j(10)
tj←ΓS+dS+1/vS+1;(11)
fuel←fuel+λα1dS+1vS+12+α2tj−ΓS+α3QijdS+1;(12)
Else(13)
tj←ai+dS+1/vS+1;(14)
fuel←fuel+λα1dijv12+α2tj−ai+α3Qijdij;(15)
End if


### 3.3. Mathematical Model

On the basis of the analysis above, the TDGVRPTW model is established, and the model parameters and decision variables are shown in Table3.Table 3
Notations of sets, parameters, and variables.
Notations for customer locations and requirementsiIndex of nodes,i=0,1,2,…,n (the depot is represented by 0)nTotal number of nodes (including the depot)VSet of nodes including the depotV0Set of nodes excluding the depot,V0=V/0ESet of arcs formed by all pairs of nodes,E=i,j:∀i∈V,j∈V,i≠jdijDistance of arci,jwiThe demand of customeriSTiService time of customeriETi,LTiTime window for starting to serve customeriNotations for vehicle parameterskIndex of vehicles,k=0,1,2,…,mKSet of vehiclesQCarrying capacity of a vehicleNotations for time-varying traffic conditionslIndex of period,l=1,2,…,hhNumber of periods,h = cardLLSet of periods,l∈LΓl,Γl+1Beginning and ending time of periodlvijlTravel speed on arci,j in period ldijlDistance of arci,j in period ltijTravel time spent on arci,jNotations for fuel consumption and CO2 emissionα1Engine power module coefficientα2Engine power speed module coefficientα3Vehicle load module coefficientλCoefficient of fuel consumptionρFuel-to-CO2 coefficientNotations for total costsf1Fixed cost of vehiclef2Drive cost rateβ1Penalty coefficient of arrival before the time windowβ2Penalty coefficient of arrival after the time windowtaxCarbon taxThe binary and continuous decision variables of the model are defined as follows:xijk0,1 variable indicating whether arc i,j is traveled 1or not 0 by vehicle kyki0,1 variable Indicating whether customer i is served 1 or not 0 by vehicle kQijContinuous variable indicating the load of vehicle on arci,jtiContinuous variable indicating the arrival time (or service starting time) at nodeiIn the above decision variables,xijk is the variable related to the arc of vehicle routing. yki represents the connection between vehicles and customers. Qij and ti are variables in vehicle driving process. Qij indicates the connection between vehicle load and arcs. ti is calculated by spatiotemporal-dependent speed, which is the crucial variable considered for TDGVRPTW. The mathematical optimization model of this paper is formulated to minimize the total distribution costs, as follows:(10)F=∑i∈V0maxETi−ti,0β1+∑i∈V0max0,ti−LTiβ2+∑k∈K∑i∈V0f1yki+∑i∈V∑j∈Vxijkf2dij+∑i,j∈E∑k∈Kxijk∑l∈Lλα1dijlvijl2+α2dijlvijl+α3Qijdijlρ×tax,(11)∑i∈V0ykiwi≤Q∀k∈K,(12)∑k∈Kyki=1∀i∈V0,(13)∑i∈V0xijk−∑i∈V0xjik=0∀j∈V0,∀k∈K,(14)∑i∈V0xijk=ykj∀j∈V0,∀k∈K,(15)∑j∈V0xijk=yki∀i∈V0,k∈K,(16)∑j∈V0x0jk=∑i∈V0xi0k≤1,(17)tj≥ti+STi+tij−M1−xijk∀i∈V,j∈V0,k∈K,(18)tj≤ti+STi+tij+M1−xijk∀i∈V,j∈V0,k∈K,(19)∑i∈VQij−∑p∈VQjp=wj∀j∈V0.Equation (10) is the objective function that represents the sum of penalty costs on earliness and tardiness, vehicle fixed costs, transportation costs, and carbon emission tax. Formulation (11) represents the vehicle maximum load constraint. Equation (12) indicates that each customer must be delivered by one vehicle. Equation (13) expresses the vehicle flow balance constraint in the node. Equations (14) and (15) represent the uniqueness constraint of node distribution vehicles. Formulation (16) represents the constraint that vehicle departures and returns to the distribution center. Formulations (17) and (18) express the relationship between the arrival time of a node and that of the previous node, where M is a positive big number and tij can be deduced from equation (7). Equation (19) indicates that the node delivery quantity is equal to the difference between the loading quantity when entering and the loading quantity when leaving.


## 3.1. Spatiotemporal-Dependent Travel Time Function

In a time-varying road network, the vehicle speed changes continuously with time. According to the historical traffic data, the morning rush hour is usually from 7: 00 am to 9: 00 am, and the evening rush hour is from 17: 00 pm to 19: 00 pm. Many vehicles emerge during the rush hour, and the average speed on the road decreases. For example, the relationship between vehicle speedv and time t in an entire day proposed by Xu et al. [20] can be approximated by the trigonometric function in equation (1), where φp, γp, and δp are constants related to the road conditions of the day:(1)vt=φ1sinγ1t+δ1,t∈0,T1⋮⋮φpsinγpt+δp,t∈Tp−1,Tp.TDVRPTW is an NP-hard problem, and it greatly increases the difficulty of solving the problem by using the function of continuous varying speed. To reduce the computational complexity caused by nonlinear speed change, this study transforms the trigonometric function of travel speed into a piecewise function in accordance with the period and takes the average value of velocity in each period as the constant velocity in that period (Figure1). The finer the time division is, the closer the speed segmentation function is to the continuous change in speed. The time interval is divided into periods L=1,2,…,h, the time of each period is recorded as Γll∈L, and the time interval is ΔI. The travel speed of the vehicle in any period l is(2)vl=∫ll+ΔIvtdtΔI,l∈L.Figure 1
Travel speed variation function.In some studies (e.g., Kwon et al. [23]), the travel time calculation of arc i,jdepends on the departure time of vehicles from node i and it assumes that the characteristics of arcs are the same. However, in a real road network, the traffic congestion of different arcs varies due to the influence of geographic location, traffic flow, and other factors. The all-day travel speed of each arc is not completely equal, instead the travel speed is spatiotemporal-dependent; that is, at the same time l∈L, the traveling speed vijl of different arc i,j is various. Therefore, the travel speed and travel time in the process of vehicle distribution not only depend on the time when the vehicle leaves node i but also relate to the arc from node i to the next node j. According to the historical traffic data on all-day arcs in a province of China in the study of Yang [28], the travel speed of different arcs in the same period obeys a normal distribution; that is, the travel speed vijl of any arc i,j∈E in period l obeys normal distribution Nvl,σ2, and the variance σ2 reflects the difference degree of travel speeds of different arcs. On this basis, the nonnegative random velocity matrix V˜ of all arcs in the road network in different periods can be obtained. A small network with 6 nodes, 21 arcs, and 9 periods is taken as an example, as shown in Figure 2. Figure 3 presents the mean value of the speed distribution of the small network, including 189 (21 × 9) random velocity variables.Figure 2
Small-size network.Figure 3
Velocity matrix.In the distribution process, vehicles are likely to cross from a period to the next period or even cross multiple periods when passing through the arci,j. Hill and Benton [29] and others did not consider the FIFO criterion when dealing with time-varying speeds across periods; that is, the vehicles that start first should arrive first starting from one point to the same next point and driving on the same arc [30]. Several studies (e.g., Ioannou et al. [31]) ensured the continuity of travel time by adjusting the speed across periods, but the calculation was complicated for solving problems. Therefore, this study refers to the work of Li et al. [32] to improve the calculation of travel time across periods, which satisfies the FIFO criterion. The vehicle is assumed to travel from node i to node j at time ai across S periods, so road section length dij is divided into d1,d2,…,dS+1. The speed in the segment ss=1,2,…,S+1 is vs∈V˜, and the time of division is Γs, as shown in Figure 4.Figure 4
Example of a time path without the nonstopping assumption.The travel time of road sectiond1 is Γ1−ai, and the travel time of subsequent road section ds is Γs−Γs−1s=2,3,…,S. Thus, an equation between section length and time is available:(3)d1=v1Γ1−ai,ds=vsΓs−Γs−1s=2,3,…,S.The driving distance of theS + 1 section is(4)dS+1=dij−∑s=1Sds=dij−v1Γ1−ai−∑s=2SvsΓs−Γs−1.Therefore, the arrival timetj of node j is(5)tj=ΓS+dS+1vS+1=v1vS+1ai+ΓS+dij−v1Γ1−∑s=2SvsΓs−Γs−1vS+1.According to equation (5), the time tj when the vehicle arrives at node j is an increasing function of time ai when the vehicle departs from node i, which meets the FIFO criterion.When the arci,j spans S periods, the following conditions should be met when the vehicle starts from node i: the vehicle starts from the first period and has not finished driving in arc i,j from the first period to the Sth period. The corresponding expression is as follows:(a)
ai<Γ1,(b)
v1Γ1−ai+∑s=2SvsΓs−Γs−1≤dij.The abovementioned arrangement yields(6)Γ1−dijv1+∑s=2Svsv1Γs−Γs−1≤ai<Γ1.The calculation formula of the arrival timetj of the next node j where any arc i,j crosses SS<h periods is as follows:(7)tj=ti+dijvijl,ai<Γl+1−dijvijl,l=1,…,L,vijlvijl+Sai+Γl+S+dij−vijlΓl+1−∑s=l+1l+S−1vijsΓs−Γs−1vijl+S,Γl+1−dijvijl+∑s=l+1l+S−1vijl+SvijlΓs−Γs−1≤ai<Γl+1,l=1,…,L−S.


## 3.2. Calculation of the Amount of Fuel Consumption

Carbon emissions are directly affected by vehicle fuel consumption. Hence, this study applies vehicle fuel consumption to calculate carbon emissions. Considering that most distribution vehicles are heavy types, a carbon emission model is established in combination with the comprehensive modal emission model (CMEM) [7]. The objective function of fuel consumption when driving in different sections is as follows:(8)fuel=λα1dv2+α2dv+α3Qd.In equation (8), λ=ζ/κφ, ζ is the fuel-air mass ratio, κ is a typical diesel calorific value, φ is the fuel conversion coefficient, engine module coefficient α1=0.5CdAfρ/1000εω, speed module coefficient α2=kNV, and vehicle load module coefficient α3=ρ0+ρmax−ρ0/Q0. The meanings and values of specific parameters are presented in Appendix A.In summary, throughout the distribution, the total fuel consumption after passing through all customer nodes is(9)FUEL=∑i,j∈E∑k∈Kxijkfuelij=∑i,j∈E∑k∈Kxijk∑l∈Lλα1dijlvijl2+α2dijlvijl+α3Qijdijl.In equation (9), xijk is a binary variable indicating whether arc i,j is traveled by vehicle k (xijk=1) or not (xijk=0). fuelij is the fuel consumption of the vehicle from node i to node j. dijl is the length of the arc i,j during time l. vijl represents the travel speed of the arc i,j during time l. And Qij represents the loading weight of the vehicle from node i to node j.Figure5 shows the relationship between fuel consumption and travel speed for a fixed vehicle type. With the increase of travel speed, the fuel consumption decreases at first. However, when the travel speed reaches a certain value, here 48.14 km/h, the fuel consumption begins to increase.Figure 5
Fuel consumption function of travel speed (with 2000 kg of load).The specific process of calculating carbon emission in consideration of spatiotemporal-dependent travel speed is shown in Algorithm1.Algorithm1: Calculation of travel time and fuel consumption.
//Depart fromi to j, assume that the period is divided into L segments(1)
Inputdij,ai,Qij, Velocity matrix V˜, Time slot division Γ, parametersλ,α1,α2, α3;(2)
Γ0←ai,v←V˜arci,j,whereΓ>ai,fuel←0,d′←dij;//Judging whether to span S periods(3)
If the criteria of crossing S periods is satisfied. //i.e. The path is divided into S + 1 segments(4)
Fork = 1 to S do begin //Recursive method(5)
dk←vkΓk−Γk−1; //The distance traveled in each period(6)
fuel←fuel + λα1dkvk2+α2Γk−Γk−1+α3Qijdk;(7)
d′←d′−dk;(8)
End for(9)
dS+1←d′; //the length of the S + 1 section//Calculate the arrival time at nodej and fuel consumption of traveling arci,j(10)
tj←ΓS+dS+1/vS+1;(11)
fuel←fuel+λα1dS+1vS+12+α2tj−ΓS+α3QijdS+1;(12)
Else(13)
tj←ai+dS+1/vS+1;(14)
fuel←fuel+λα1dijv12+α2tj−ai+α3Qijdij;(15)
End if


## 3.3. Mathematical Model

On the basis of the analysis above, the TDGVRPTW model is established, and the model parameters and decision variables are shown in Table3.Table 3
Notations of sets, parameters, and variables.
Notations for customer locations and requirementsiIndex of nodes,i=0,1,2,…,n (the depot is represented by 0)nTotal number of nodes (including the depot)VSet of nodes including the depotV0Set of nodes excluding the depot,V0=V/0ESet of arcs formed by all pairs of nodes,E=i,j:∀i∈V,j∈V,i≠jdijDistance of arci,jwiThe demand of customeriSTiService time of customeriETi,LTiTime window for starting to serve customeriNotations for vehicle parameterskIndex of vehicles,k=0,1,2,…,mKSet of vehiclesQCarrying capacity of a vehicleNotations for time-varying traffic conditionslIndex of period,l=1,2,…,hhNumber of periods,h = cardLLSet of periods,l∈LΓl,Γl+1Beginning and ending time of periodlvijlTravel speed on arci,j in period ldijlDistance of arci,j in period ltijTravel time spent on arci,jNotations for fuel consumption and CO2 emissionα1Engine power module coefficientα2Engine power speed module coefficientα3Vehicle load module coefficientλCoefficient of fuel consumptionρFuel-to-CO2 coefficientNotations for total costsf1Fixed cost of vehiclef2Drive cost rateβ1Penalty coefficient of arrival before the time windowβ2Penalty coefficient of arrival after the time windowtaxCarbon taxThe binary and continuous decision variables of the model are defined as follows:xijk0,1 variable indicating whether arc i,j is traveled 1or not 0 by vehicle kyki0,1 variable Indicating whether customer i is served 1 or not 0 by vehicle kQijContinuous variable indicating the load of vehicle on arci,jtiContinuous variable indicating the arrival time (or service starting time) at nodeiIn the above decision variables,xijk is the variable related to the arc of vehicle routing. yki represents the connection between vehicles and customers. Qij and ti are variables in vehicle driving process. Qij indicates the connection between vehicle load and arcs. ti is calculated by spatiotemporal-dependent speed, which is the crucial variable considered for TDGVRPTW. The mathematical optimization model of this paper is formulated to minimize the total distribution costs, as follows:(10)F=∑i∈V0maxETi−ti,0β1+∑i∈V0max0,ti−LTiβ2+∑k∈K∑i∈V0f1yki+∑i∈V∑j∈Vxijkf2dij+∑i,j∈E∑k∈Kxijk∑l∈Lλα1dijlvijl2+α2dijlvijl+α3Qijdijlρ×tax,(11)∑i∈V0ykiwi≤Q∀k∈K,(12)∑k∈Kyki=1∀i∈V0,(13)∑i∈V0xijk−∑i∈V0xjik=0∀j∈V0,∀k∈K,(14)∑i∈V0xijk=ykj∀j∈V0,∀k∈K,(15)∑j∈V0xijk=yki∀i∈V0,k∈K,(16)∑j∈V0x0jk=∑i∈V0xi0k≤1,(17)tj≥ti+STi+tij−M1−xijk∀i∈V,j∈V0,k∈K,(18)tj≤ti+STi+tij+M1−xijk∀i∈V,j∈V0,k∈K,(19)∑i∈VQij−∑p∈VQjp=wj∀j∈V0.Equation (10) is the objective function that represents the sum of penalty costs on earliness and tardiness, vehicle fixed costs, transportation costs, and carbon emission tax. Formulation (11) represents the vehicle maximum load constraint. Equation (12) indicates that each customer must be delivered by one vehicle. Equation (13) expresses the vehicle flow balance constraint in the node. Equations (14) and (15) represent the uniqueness constraint of node distribution vehicles. Formulation (16) represents the constraint that vehicle departures and returns to the distribution center. Formulations (17) and (18) express the relationship between the arrival time of a node and that of the previous node, where M is a positive big number and tij can be deduced from equation (7). Equation (19) indicates that the node delivery quantity is equal to the difference between the loading quantity when entering and the loading quantity when leaving.


## 4. Problem-Solving Method

TDGVRPTW is an NP-hard problem. In this study, the spatiotemporal division of travel speed is added to the classical VRP, thus making the solution increasingly difficult. A metaheuristic algorithm is usually used to solve this kind of problem. The genetic algorithm has strong global search ability and robustness, and it is mature enough for application in vehicle scheduling problems. However, it easily falls into local convergence. Therefore, this study designs a hybrid adaptive genetic algorithm with elite local search (HAGA_ELS). The basic idea is to obtain the initial population by using the Clarke and Wright Saving (CWS) heuristic in the study of Marković et al. [33] and then perform crossover and mutation on individual populations in the evolution process. The elite individuals in the original population are optimized by neighborhood search and reinserted into the new population to keep the population size unchanged and increase the population diversity. Meanwhile, the crossover and mutation rates are adaptively adjusted according to the iteration times and individual fitness changes. The operation flow of HAGA_ELS is shown in Figure 6.Figure 6
HAGA_ELS algorithmic flow.


### 4.1. Initial Solution

In this study, population size is set asN, and each chromosome is composed of n (n is the number of customers) integers. The sequence of gene bits indicates the order in which distribution vehicles visit customers. To improve the quality of the initial population, this study uses the CWS heuristic and random permutation operator to generate the initial population. The steps of the random permutation operator are as follows: an integer sequence containing customers is randomly generated, and vehicles are assigned to the route in accordance with the vehicle load constraint. When the vehicle cannot meet the next customer demand, another vehicle is called, and the current gene of the chromosome is set to 0. The driving sequence of a vehicle appears as subrouting. If the number of required vehicles exceeds the number of existing vehicles, then the route is an infeasible solution. The steps above are repeated until the chromosome number reaches population size N.


### 4.2. Fitness Function

The fitness of each chromosome in the population is constructed by the model objective function equation (10). In this study, the pressure principle is used to deepen the gap of fitness between individuals, thus increasing the probability of excellent individuals entering the next generation [34]. The fitness function is shown in equation (20), where a is a constant greater than 1, fi is the fitness of chromosome i, and zi is the objective function value corresponding to chromosome i:(20)fi=amaxzin−zi.


### 4.3. Selection Operator

In this study, elite and roulette strategies are combined to implement the selection operation. The specific steps are as follows: with the descending order of chromosome fitness value,N1−Φ of the top individuals is selected for elite evolution operation (Φ is the generation gap of population evolution) to ensure that the outstanding individuals in the parent will not be lost due to crossover and mutation. Then, the roulette strategy is adopted to select NΦ individuals from their parents. The selection probability of each chromosome is pi=fi/∑i∈Nfi, and the greater the fitness is, the greater the probability of being selected is. The selected individuals are subjected to crossover mutation operation, and the elite evolutionary individuals are reinserted into the population after crossover and mutation so that the population size is always N.


### 4.4. Improved Crossover and Mutation Operators

For the traditional adaptive genetic algorithm, when the larger fitness value of the crossover operation individual is equal to the maximum fitness value of the population, the crossover and mutation rates are zero; additionally, the influence of evolution time on crossover and mutation rate is not considered. Therefore, the improved method of cloud adaptive adjustment proposed by Dai et al. [35] is adopted. An adaptive adjustment algorithm of crossover and mutation rate is designed based on the robustness and efficiency of cloud theory. In the early stage of evolution, large crossover and mutation rates make the individuals with low fitness participate in the crossover and mutation operation with high probability to produce excellent individuals. In the late stage of evolution, small crossover and mutation rates make individuals with high fitness participate in the crossover and mutation operation with low probability to protect outstanding individuals from being destroyed. The algorithm for generating crossover rate Pc and mutation rate Pm is as follows:(21)Pc=k1e−f−favg2/2Normfmax−favg/C1,fmax−favg/C1C22,f≥favg,k3,f≤favg,Pm=k2e−f′−favg2/2Normfmax−favg/C3,fmax−favg/C3C42,f′≥favg,k4,f′≤favg,where f is the larger fitness value among the crossover individuals, f′ is the fitness value of the mutant individual, fmax is the maximum fitness value in the population, and favg is the average fitness value. Norm is a normal random number generator. C1∼C4 are the control coefficients (C1=C3=T+1/6N, C2=C4=15−T−N/22), and k1∼k4 are constants in the range of 0,1. The control coefficient in the algorithm is appropriately adjusted in accordance with evolutional generationT, search accuracy, and other factors.


#### 4.4.1. Hybrid Crossover Operator

To avoid the ineffectiveness of crossover caused by the same crossover genomes, two hybrid crossover operators are used for crossover operations to increase individual diversity. When the crossover genomes are different, the Partial-Mapped Crossover (PMX) crossover operator is employed to exchange the crossover genomes of parents A and B, and a mapping relationship is established. The conflicting genes in offspring A and B are replaced in accordance with the mapping relationship to ensure that the coding of offspring obtained by crossover is nonrepetitive. When the cross genomes are the same, we adopt the improved Order Crossover (OX) crossover operator in the work of Zhang and Li [34]. The crossover regions of parents A and B are placed at the front and back of the other parent, respectively, and the genes in parents A and B that are duplicated in the inserted cross genomes are deleted to form two new children. The crossover process is shown in Figure 7.Figure 7
Crossover operator. (a) PMX crossover operator. (b) Improved OX crossover operator.
(a)(b)


#### 4.4.2. Reverse Mutation Operator

In this study, the mutation operation is improved by the reversal operator adopted by Barma et al. [36]. Two nodes in the mutated chromosome are randomly selected and then inserted reversely after the selected two nodes break. If the fitness value of the new chromosome after reversal mutation is better than the fitness value of the original chromosome, the reversal mutation operation is effective. The new chromosome after mutation is retained, and the original chromosome is deleted so that the individual is always mutated in an improved direction. The specific mutation operation is shown in Figure 8.Figure 8
Reverse mutation operator.


### 4.5. Elite Local Search

To avoid the premature convergence by the traditional elite retention strategy, the neighborhood search operator is designed to further optimize elite individuals and improve the local optimum caused by direct retention of elite individuals. To search the new solution space effectively, the four neighborhood operators in the study of Xu et al. [37] are employed to search the elite individuals. These four operators are cyclic operator, 1-0 exchange, 2-opt∗ exchange, and 3-opt exchange, respectively. If the new individual after the neighborhood search is better than the original elite individual, then the elite optimization operation is effective. Afterward, the improved elite individual is recombined with the individuals after crossover and mutation operations to ensure that the optimal individual of the new population is always better than or equal to the optimal individual of the previous generation (see Algorithm 2).(1)
Cyclic operator: two node positions in a subrouting are randomly selected, and the genes between the two node positions are circularly moved. As shown in Figure9(a), the positions of nodes 6 and 5 are selected, and the intermediate part6¯2¯5¯ is sequentially moved forward, i.e., 2¯5¯6¯.(2)
1-0 exchange: it is also known as swap move, which means inserting a node in a subrouting into another subrouting. The exchange operation is illustrated in Figure9(b).(3)
2-opt∗ exchange: a node position in the two subroutings is selected, and the tails of the two subroutings are exchanged in accordance with the node position, as shown in Figure9(c).(4)
3-opt exchange: four continuous pointsi,i+1,j,j+1 in a certain subrouting are transformed into i+1,i,j+1,j, as shown in Figure 9(d).Algorithm2: Procedure of ELS.
(1)InputXGood,fitness(2)Fori = 1 to size (XGood,1) do begin(3)XGood’←XGoodi,: //Update XGood’after each step of operation//Cyclic operator. Choosingr1 and r2 at random, which are in the same subrouting(4)New_XGoodr2←XGood’r1(5)Fork = r1 to r2−1 do begin(6)New_XGoodk←XGood’k+1(7)End for//1-0 swap. Choosingr1 andr2 in different subroutings(8) RemoveNew_XGoodi,r1(9)New_XGoodi,r2←XGood’i,r1;New_XGoodi,r2+1:end←XGood’i,r2:end//2-opt swap. Choosingr1 and r2 in different subroutings(10)New_XGoodi,r1:where node is0←XGood’i,r2where node is0)(11) Similarly for the other subrouting wherer2 is located//3-opt swap. Choosingr1 andr2 in the same subrouting(12)New_XGoodi,r1←XGood’i,r1+1;New_XGoodi,r1+1←XGood’i,r1(13)New_XGoodi,r2←XGood’i,r2+1;New_XGoodi,r2+1←XGood’i,r2//Judging the effectiveness of elite evolution(14)new_fitness←FNew_XGood(15)Ifnew_fitness>fitnessi(16)XGoodi,:←New_XGood(17)End if(18)End forFigure 9
Elite neighborhood search operation. (a) Cyclic operator. (b) 1-0 swap. (c) 2-opt∗ swap. (d) 3-opt swap.
(a)(b)(c)(d)


### 4.6. Termination Condition

As the process of the metaheuristic algorithm is uncertain, it is difficult to determine a number of terminal iterations to stop the calculation. Considering the change of the fitness value would get smaller, it is reasonable to stop the process when the best solution is not improved in successive generationsΩ.


## 4.1. Initial Solution

In this study, population size is set asN, and each chromosome is composed of n (n is the number of customers) integers. The sequence of gene bits indicates the order in which distribution vehicles visit customers. To improve the quality of the initial population, this study uses the CWS heuristic and random permutation operator to generate the initial population. The steps of the random permutation operator are as follows: an integer sequence containing customers is randomly generated, and vehicles are assigned to the route in accordance with the vehicle load constraint. When the vehicle cannot meet the next customer demand, another vehicle is called, and the current gene of the chromosome is set to 0. The driving sequence of a vehicle appears as subrouting. If the number of required vehicles exceeds the number of existing vehicles, then the route is an infeasible solution. The steps above are repeated until the chromosome number reaches population size N.


## 4.2. Fitness Function

The fitness of each chromosome in the population is constructed by the model objective function equation (10). In this study, the pressure principle is used to deepen the gap of fitness between individuals, thus increasing the probability of excellent individuals entering the next generation [34]. The fitness function is shown in equation (20), where a is a constant greater than 1, fi is the fitness of chromosome i, and zi is the objective function value corresponding to chromosome i:(20)fi=amaxzin−zi.


## 4.3. Selection Operator

In this study, elite and roulette strategies are combined to implement the selection operation. The specific steps are as follows: with the descending order of chromosome fitness value,N1−Φ of the top individuals is selected for elite evolution operation (Φ is the generation gap of population evolution) to ensure that the outstanding individuals in the parent will not be lost due to crossover and mutation. Then, the roulette strategy is adopted to select NΦ individuals from their parents. The selection probability of each chromosome is pi=fi/∑i∈Nfi, and the greater the fitness is, the greater the probability of being selected is. The selected individuals are subjected to crossover mutation operation, and the elite evolutionary individuals are reinserted into the population after crossover and mutation so that the population size is always N.


## 4.4. Improved Crossover and Mutation Operators

For the traditional adaptive genetic algorithm, when the larger fitness value of the crossover operation individual is equal to the maximum fitness value of the population, the crossover and mutation rates are zero; additionally, the influence of evolution time on crossover and mutation rate is not considered. Therefore, the improved method of cloud adaptive adjustment proposed by Dai et al. [35] is adopted. An adaptive adjustment algorithm of crossover and mutation rate is designed based on the robustness and efficiency of cloud theory. In the early stage of evolution, large crossover and mutation rates make the individuals with low fitness participate in the crossover and mutation operation with high probability to produce excellent individuals. In the late stage of evolution, small crossover and mutation rates make individuals with high fitness participate in the crossover and mutation operation with low probability to protect outstanding individuals from being destroyed. The algorithm for generating crossover rate Pc and mutation rate Pm is as follows:(21)Pc=k1e−f−favg2/2Normfmax−favg/C1,fmax−favg/C1C22,f≥favg,k3,f≤favg,Pm=k2e−f′−favg2/2Normfmax−favg/C3,fmax−favg/C3C42,f′≥favg,k4,f′≤favg,where f is the larger fitness value among the crossover individuals, f′ is the fitness value of the mutant individual, fmax is the maximum fitness value in the population, and favg is the average fitness value. Norm is a normal random number generator. C1∼C4 are the control coefficients (C1=C3=T+1/6N, C2=C4=15−T−N/22), and k1∼k4 are constants in the range of 0,1. The control coefficient in the algorithm is appropriately adjusted in accordance with evolutional generationT, search accuracy, and other factors.


### 4.4.1. Hybrid Crossover Operator

To avoid the ineffectiveness of crossover caused by the same crossover genomes, two hybrid crossover operators are used for crossover operations to increase individual diversity. When the crossover genomes are different, the Partial-Mapped Crossover (PMX) crossover operator is employed to exchange the crossover genomes of parents A and B, and a mapping relationship is established. The conflicting genes in offspring A and B are replaced in accordance with the mapping relationship to ensure that the coding of offspring obtained by crossover is nonrepetitive. When the cross genomes are the same, we adopt the improved Order Crossover (OX) crossover operator in the work of Zhang and Li [34]. The crossover regions of parents A and B are placed at the front and back of the other parent, respectively, and the genes in parents A and B that are duplicated in the inserted cross genomes are deleted to form two new children. The crossover process is shown in Figure 7.Figure 7
Crossover operator. (a) PMX crossover operator. (b) Improved OX crossover operator.
(a)(b)


### 4.4.2. Reverse Mutation Operator

In this study, the mutation operation is improved by the reversal operator adopted by Barma et al. [36]. Two nodes in the mutated chromosome are randomly selected and then inserted reversely after the selected two nodes break. If the fitness value of the new chromosome after reversal mutation is better than the fitness value of the original chromosome, the reversal mutation operation is effective. The new chromosome after mutation is retained, and the original chromosome is deleted so that the individual is always mutated in an improved direction. The specific mutation operation is shown in Figure 8.Figure 8
Reverse mutation operator.


## 4.4.1. Hybrid Crossover Operator

To avoid the ineffectiveness of crossover caused by the same crossover genomes, two hybrid crossover operators are used for crossover operations to increase individual diversity. When the crossover genomes are different, the Partial-Mapped Crossover (PMX) crossover operator is employed to exchange the crossover genomes of parents A and B, and a mapping relationship is established. The conflicting genes in offspring A and B are replaced in accordance with the mapping relationship to ensure that the coding of offspring obtained by crossover is nonrepetitive. When the cross genomes are the same, we adopt the improved Order Crossover (OX) crossover operator in the work of Zhang and Li [34]. The crossover regions of parents A and B are placed at the front and back of the other parent, respectively, and the genes in parents A and B that are duplicated in the inserted cross genomes are deleted to form two new children. The crossover process is shown in Figure 7.Figure 7
Crossover operator. (a) PMX crossover operator. (b) Improved OX crossover operator.
(a)(b)


## 4.4.2. Reverse Mutation Operator

In this study, the mutation operation is improved by the reversal operator adopted by Barma et al. [36]. Two nodes in the mutated chromosome are randomly selected and then inserted reversely after the selected two nodes break. If the fitness value of the new chromosome after reversal mutation is better than the fitness value of the original chromosome, the reversal mutation operation is effective. The new chromosome after mutation is retained, and the original chromosome is deleted so that the individual is always mutated in an improved direction. The specific mutation operation is shown in Figure 8.Figure 8
Reverse mutation operator.


## 4.5. Elite Local Search

To avoid the premature convergence by the traditional elite retention strategy, the neighborhood search operator is designed to further optimize elite individuals and improve the local optimum caused by direct retention of elite individuals. To search the new solution space effectively, the four neighborhood operators in the study of Xu et al. [37] are employed to search the elite individuals. These four operators are cyclic operator, 1-0 exchange, 2-opt∗ exchange, and 3-opt exchange, respectively. If the new individual after the neighborhood search is better than the original elite individual, then the elite optimization operation is effective. Afterward, the improved elite individual is recombined with the individuals after crossover and mutation operations to ensure that the optimal individual of the new population is always better than or equal to the optimal individual of the previous generation (see Algorithm 2).(1)
Cyclic operator: two node positions in a subrouting are randomly selected, and the genes between the two node positions are circularly moved. As shown in Figure9(a), the positions of nodes 6 and 5 are selected, and the intermediate part6¯2¯5¯ is sequentially moved forward, i.e., 2¯5¯6¯.(2)
1-0 exchange: it is also known as swap move, which means inserting a node in a subrouting into another subrouting. The exchange operation is illustrated in Figure9(b).(3)
2-opt∗ exchange: a node position in the two subroutings is selected, and the tails of the two subroutings are exchanged in accordance with the node position, as shown in Figure9(c).(4)
3-opt exchange: four continuous pointsi,i+1,j,j+1 in a certain subrouting are transformed into i+1,i,j+1,j, as shown in Figure 9(d).Algorithm2: Procedure of ELS.
(1)InputXGood,fitness(2)Fori = 1 to size (XGood,1) do begin(3)XGood’←XGoodi,: //Update XGood’after each step of operation//Cyclic operator. Choosingr1 and r2 at random, which are in the same subrouting(4)New_XGoodr2←XGood’r1(5)Fork = r1 to r2−1 do begin(6)New_XGoodk←XGood’k+1(7)End for//1-0 swap. Choosingr1 andr2 in different subroutings(8) RemoveNew_XGoodi,r1(9)New_XGoodi,r2←XGood’i,r1;New_XGoodi,r2+1:end←XGood’i,r2:end//2-opt swap. Choosingr1 and r2 in different subroutings(10)New_XGoodi,r1:where node is0←XGood’i,r2where node is0)(11) Similarly for the other subrouting wherer2 is located//3-opt swap. Choosingr1 andr2 in the same subrouting(12)New_XGoodi,r1←XGood’i,r1+1;New_XGoodi,r1+1←XGood’i,r1(13)New_XGoodi,r2←XGood’i,r2+1;New_XGoodi,r2+1←XGood’i,r2//Judging the effectiveness of elite evolution(14)new_fitness←FNew_XGood(15)Ifnew_fitness>fitnessi(16)XGoodi,:←New_XGood(17)End if(18)End forFigure 9
Elite neighborhood search operation. (a) Cyclic operator. (b) 1-0 swap. (c) 2-opt∗ swap. (d) 3-opt swap.
(a)(b)(c)(d)


## 4.6. Termination Condition

As the process of the metaheuristic algorithm is uncertain, it is difficult to determine a number of terminal iterations to stop the calculation. Considering the change of the fitness value would get smaller, it is reasonable to stop the process when the best solution is not improved in successive generationsΩ.


## 5. Computational Experiment and Analyses

In this section, the experimental design is described in Section5.1. Parameter setting is provided in Section 5.2. The performance of the algorithm is evaluated in Section 5.3. In Sections 5.4 and 5.5, three speed patterns are compared and the sensitivity of the model is assessed. All experiments are implemented in Matlab R2018b in a PC with Intel Core i5 2.4 GHz processor, 8 GB RAM, and Windows 10 professional operating system.


### 5.1. Experimental Design

Nine benchmark instances are used to assess the effectiveness of HAGA_ELS (see, e.g., Augerat et al. [38], Breedam [39], and Christofides and Eilon [40]). In these instances, the geographical distribution of customers in E instances is randomly and evenly distributed, the geographical location of customers in A instances has semiclustering characteristics, and the geographical location of customers in B instances is in accordance with the clustering characteristics. The time windows are set to [60,600], which allows early arrivals but does not allow delayed arrivals. n represents the number of nodes in the instances, NV represents the number of vehicles used, Capacity represents the maximum load of vehicles, and r0 represents the ratio of the total demand to the total capacity. Table 4 lists the characteristics and basic information of the instances.Table 4
Characteristics of benchmark instances.
No.InstancesnNVCapacityr01E_n33_k433480000.9182E_n45_k445420100.8983E_n51_k55151600.9714A_n37_k63761000.9505A_n55_k95591000.9326A_n65_k96591000.9747B_n35_k53551000.8748B_n56_k75671000.8809B_n78_k1078101000.937The instance for sensitivity analysis is designed following specific distributions.nn∈10,20,40 customers are generated randomly, where the location x,y obeys the uniform distribution of U−50,50, customer demand Di obeys the uniform distribution of U0,10, and the customer time window ETi,LTi is set to (60,600). On the basis of travel speed function vt of equation (1), travel speed vijt of different arcs in different periods is randomly generated so that vijt obeys Normalvt,σ2.


### 5.2. Parameter Setting

The parameters, for instance, are set as follows: fixed vehicle costf1=110 yuan/vehicle; time penalty cost β1=50 yuan/hour and β2=100 yuan/hour; fuel cost rate f2=0.8 yuan/L; fuel consumption parameters are set as a1=0.111,a2=36.0302,anda3=0.000059 [41]; fuel-to-CO2 coefficient ρ=2.61kg/L; and carbon tax=1.5 yuan/t. The parameters for the CMEM model are given in Appendix A.The parameters to be determined for the algorithm include the cloud modelk1∼k4, the population size N, and successive generations Ω. For the parameters k1∼k4 and N, the trial and error are executed to achieve the best configuration of HAGA_ELS in a set of experiments. We test the population size of 10, 50, 100, and 200 and use crossover related probabilities k1,k3 ranging from 0 to 1 with an increment of 0.5 and mutation related probabilities k2,k4 ranging from 0 to 1 with an increment of 0.5. Thus, we have a 3D matrix holding all the possible combinations. Because finding the appropriate configuration over all the problem instances is time consuming, we opt to perform on a particular instance B_n78_k10 with maximum nodes. The algorithm has been run 10 times for each combination within fixed generations of 500. Detailed results are coming in Appendix B. The results show that when the population size increases to 200, the solution does not improve significantly in general, but the CPU time increases by 2-3 times, so the population size could be defined less than 200. We also find that in most experiments the difference between solutions in last consecutive iterations of 50 final is insignificant. From all the test runs, the best configuration has a population size of 100, crossover related probabilities of 1 (k1=k3=1), mutation related probabilities of 0.5 (k2=k4=0.5), and successive generations of 50 (for HAGA_ELS).


### 5.3. Evaluating the Performance of HAGA_ELS on TDVRPTW

To evaluate the proposed algorithm HAGA_ELS, two GA-based algorithms and two classical constructive heuristics are used for comparison. Classical constructive heuristics, sequential insertion algorithm (SI), and Clarke and Wright Saving (CWS) heuristic are, respectively, discussed in Avdoshin and Beresneva [42] and Marković et al. [33]. GA-based algorithms, standard adaptive genetic algorithm (SAGA), and elite strategy genetic algorithm (EGA) were coded as follows:(1)
SAGA is the same as HAGA_ELS, except for the elite local search described in Section4.5.(2)
EGA adopts elite local search in HAGA_ELS, but the crossover and mutation are the same as the standard genetic algorithm (SGA).Table5 lists the results obtained by the three metaheuristic algorithms. Each instance is tested on the three metaheuristic algorithms 10 times. Avg represents the average value of 10 calculation results, CV represents the deviation coefficient (CV=STD/Mean) under each algorithm. Gap_s% and Gap_e%, respectively, represent the optimized percentage of the average values of SAGA and EGA, for example, Gap_s%=SAGAAvg−HAGA_ELSAvg/HAGA_ELSAvg∗100.Table 5
Comparison of three GA-based algorithms.
InstancesHAGA_ELSSAGAEGAGap_s%Gap_e%AvgCVAvgCVAvgCVE_n33_k4863.564.22876.983.60879.912.441.551.89E_n45_k4740.763.02755.057.24753.276.231.931.69E_n51_k5579.686.32610.588.59597.287.375.333.04A_n37_k6986.344.031009.915.401004.585.832.391.85A_n55_k91125.338.141154.335.971147.965.172.582.01A_n65_k91276.165.361330.866.821304.318.674.292.21B_n35_k5962.902.79966.106.02965.393.900.330.26B_n56_k7724.173.96736.267.02730.386.571.670.86B_n78_k101344.556.321388.136.261360.617.623.241.19Average955.944.91980.916.32971.525.982.591.67As shown in Table5, the average optimization percentage compared with SAGA and EGA is 2.59% and 1.67%, respectively. The average results of the nine instances obtained by HAGA_ELS are better than those obtained by the two other algorithms. HAGA_ELS has an average deviation coefficient of 4.91% for nine instances, which is much smaller than other two algorithms compared to 6.32% of SAGA and 5.98% of EGA. The result shows that HAGA_ELS is rather effective and stable.Table6 lists the results of the three metaheuristic algorithms and the known optimal solution (OPT). The iterations represent the number of average iterations to obtain the optimal solution in 10 experiments, and Gap_h%, Gap_s%, and Gap_e%, respectively, represent the gap between the solution of three algorithms and the known optimal solution. For example, Gap_h%=HAGA_ELSBest−opt/opt∗100. As shown in Table 6, the average gap between HAGA_ELS and OPT is 3.87%, while the average gap of SAGA and EGA is 6.61% and 5.54%, respectively. Given that a time window is added to the calculation instance in this study, the deviation is still acceptable. Therefore, HAGA_ELS has firm solving capability for instances with different characteristics and sizes.Table 6
Comparison of optimal solutions of GA-based algorithms.
InstancesOPTHAGA_ELSSAGAEGAGap_h%Gap_s%Gap_e%BestIterationsBestIterationsBestIterationsE_n33_k4835853.44208871.79201874.731732.214.414.76E_n45_k4728737.18453747.40420735.444481.262.661.02E_n51_k5521575.58375601.58342590.7431110.4815.4713.39A_n37_k6949979.712021001.311681001.311343.245.515.51A_n55_k910731114.102481145.742071143.641893.836.786.58A_n65_k911741262.024171322.443291286.072707.5012.649.55B_n35_k5955957.11204959.08197961.721760.220.430.70B_n56_k7707715.06295728.36289721.132801.143.022.00B_n78_k1012661328.635251373.953811346.383644.958.536.35Average912946.98325972.41282962.352603.876.615.54Figure10 shows the convergence graph of the three metaheuristic algorithms for large instances. Figures 10(a)–10(c) response to E_n51_k5, A_n65_k9, and B_n78_k10 instances, respectively. According to the graph, HAGA_ELS can obtain more qualified solution within more iterations.Figure 10
Convergence comparison.
(a)(b)(c)The computational results of HAGA_ELS and classical constructive heuristics are shown in Table7. Gap_h% is much smaller than Gap_si% and Gap_cws%, which means HAGA_ELS performs better than SI and CWS in any instances. The average gap of HAGA_ELS is 3.87%, while the average gaps of SI and CWS are all over 15%. Thus, compared with constructive heuristics, HAGA_ELS has ability to obtain the better solution.Table 7
Comparison of optimal solutions of HAGA_ELS and constructive heuristics.
InstancesOPTHAGA-ELSSICWSGap_h%Gap_si%Gap_cws%E_n33_k4835853.44888.82937.772.216.4512.31E_n45_k4728737.18881.00964.601.2621.0232.50E_n51_k5521575.58638.28771.8510.4822.5148.15A_n37_k6949979.711175.901065.503.2423.9112.28A_n55_k910731114.101394.101257.703.8329.9317.21A_n65_k911741262.021571.701639.607.5033.8839.66B_n35_k5955957.111007.501245.100.225.5030.38B_n56_k7707715.06869.39861.501.1422.9721.85B_n78_k1012661328.631424.401704.204.9512.5134.61Average912946.981094.571160.873.8719.8527.66


### 5.4. Comparison of Three Patterns

The second experiment is a small-size experiment with 10 customers to verify the validity of the model, and it is carried out under three patterns (T&S, T, and S). Each experiment is run 10 times. Table8 shows a cost comparison for the three patterns. Given that the fixed costs are the same, they are omitted in the Table 8, in which Gap_T% and Gap_S%, respectively, represent the difference between T&S pattern and T and or S patterns. The optimization percentage of each kind of costs between S pattern and T&S pattern is calculated as Gap%=cost′−cost/cost. Figure 11 shows the distribution scheme under the three patterns with the corresponding carbon emission costs and average travel speed. Figure 12 shows the travel speed of each road section under the three patterns.Table 8
Cost item deviations of the different patterns at a small-size instance.
Cost itemT&STSGap_T%Gap_S%Total cost670.74741.21851.5410.5126.96Transportation cost203.66197.52191.27−3.01−6.08Carbon cost170.58182.19199.936.8117.21Time penalty cost76.49141.5240.3484.99214.21Figure 11
Comparison of paths under three patterns.Figure 12
Speeds of small-size instance under three patterns.Table8 shows the comparison of costs for the three patterns. Although the transportation cost under the T&S pattern is slightly higher, the total cost of T and S patterns is, respectively, 10.51% and 26.96% higher than that of the T&S pattern, mainly because the time penalty cost is greatly reduced under the T&S pattern. This finding shows that applying the spatiotemporal-varying speed facilitates to reduce the default of the customer required.Figure11 shows the distribution scheme under the three patterns, and Figure 12 shows the average travel speed of each arc under the three patterns. As shown in Figure 11, the scheme obtained by the T&S pattern has the lowest carbon emission cost and the fastest average travel speed. The average travel speeds under T&S, T, and S patterns are 26.77, 24.38, and 23.17 km/h, respectively. Figure 12 indicates that under the T&S pattern, the travel routing can avoid traffic congestion, so that the distribution has a high travel speed. The average travel speed of arcs can be seen in Figure 11. The maximum travel speed under the three patterns does not reach the inflection point of the carbon emission cost function (Equation (9)). Hence, the carbon emission cost drops with the increase in speed. The T&S pattern improves travel speed by avoiding congestion, thus reducing the carbon emission cost.As shown in Figure12, the T&S pattern mainly selects the route to avoid congestion through the node sequence. Thus, its travel speed is generally at a high level. Detailed analysis of Figure 11 indicates that the T&S pattern avoids the congestion in the fourth arc (a⟶c) and the fifth arc (c⟶b). In these arcs, the order of node arrival is c⟶a⟶b and b⟶a⟶c in T and S patterns, respectively, and the order of node arrival is changed to a⟶c⟶b in the T&S pattern. This finding shows that the T&S pattern can help vehicles plan the node arrival sequence and avoid the traffic congestion, thus reducing carbon emission cost and improving distribution speed.


### 5.5. Comparison of Different Factor Impacts on Three Patterns

To assess the effectiveness of the T&S pattern further, we simulate and analyze different customer sizes in Section5.5.1 and different speed distributions in Section 5.5.2.


#### 5.5.1. Analysis of Different Customer Size

We select three customer groups with different sizes (n=10, n=20, and n=40) for testing and adopt three dynamic speed patterns, namely, T&S, T, and S. Nine groups of experiments are performed, and each group of experiments is run 10 times. We obtain the best results of the three dynamic speed patterns with different sizes and compare them in terms of cost, travel time, and speed. Under σ=5, we generate speed data that obey a normal distribution under different arcs in the same period to ensure a random difference of speed in road network. At the same time, the travel speed on each arc obeys the function law of the period, which reflects the real traffic network to a certain extent.We report the computational results of the nine groups of experiments in Table9, in which Gap_T%=Item_T−Item_T&S/Item_T&S×100, Gap_S%=Item_S−Item_T&S/Item_T&S×100. Gap_T% and Gap_S% are greater than zero under different customer sizes n; that is, the best result is obtained by using the spatiotemporal travel speed pattern (T&S). The Gap_T% in cost is over 10%, and that of Gap_S% is over 15%. The routing optimization under the T&S pattern is superior to those under T and S patterns in terms of cost.Table 9
Numerical results of different patterns under different sizes.
Customer sizeItemT&STSGap_T%Gap_S%n=10Cost670.74741.21851.5410.5126.96Travel time1168.381198.301250.512.567.03Average speed26.7724.3823.178.9313.45Average———7.3315.81n=20Cost1026.021226.81222.5219.5719.15Travel time2245.092430.232456.088.259.40Average speed28.3925.8025.309.1210.88Average———12.3113.14n=40Cost2730.463254.573650.3419.1933.69Travel time5332.415813.85855.059.039.80Average speed31.0327.7627.4110.5411.67Average———12.9218.39Figure13(a) shows the total cost deviation between T&S pattern and T (or S) pattern under different customer sizes, and Figure 13(b) presents the total travel time deviation between T&S pattern and T (or S) pattern. With the increase in customer size, the total cost deviation and total travel time deviation also increase. Table 9 also indicates that the average optimization percentage Gap% appears an upward trend with customers increasing. Gap_T% increases from 7.33% to 12.92%, and Gap_S% increases from 15.81% to 18.39%. Therefore, with the increase in customer size, the T&S pattern performs much better than T and S patterns by reducing the distribution cost and improving the distribution efficiency.Figure 13
Total cost and total time deviation of different patterns under different sizes.
(a)(b)


#### 5.5.2. Analysis of Different Speed Deviation

Three different standard deviations of speed distribution (σ=10, σ=20, and σ=40) are selected to randomly generate a spatiotemporal speed matrix. The travel speed vijl of different arcs in different periods obeys the normal distribution of Normalvt,σ2. (vt is the time travel speed function (equation (1)). By changing the variance of the normal distribution, the difference of travel speed in arcs can be adjusted. The greater the variance is, the greater the difference in travel speed is among different arcs. In the same manner, three different spatiotemporal-varying speed matrices are tested with the three dynamic speed patterns. To avoid the compound effect caused by the change in customer group size, we set the customer group size to be constant (n=20), and each group of experiments is run 10 times to obtain the optimal results of three speed patterns under different speed matrices. A comparative analysis is performed in terms of cost, travel time, and speed.The numerical results of each group of experiments are shown in Table10, which illustrates that Gap_T% and Gap_S% are greater than or equal to 0 under the velocity matrix generated by different variances σ2. Thus, the spatiotemporal speed pattern (T&S) is the best method under the velocity matrix generated by different variances. When the variance of the distribution is small, that is, σ=2, the difference in spatial speed is small. In this case, the running results of the T&S and T patterns are consistent, so there is no obvious difference whether speed spatiality is considered or not. However, with the increase in speed distribution variance σ2, Gap_T% increases continuously and reaches at 17.53% under σ=8. Meanwhile, Gap_S% increases from 3.55% to 17.53%. When the difference of spatial speed amplifies, the T&S pattern becomes superior to T and S patterns, and the deficiency of the T pattern that does not consider speed spatiality becomes increasingly apparent.Table 10
Numerical results of different patterns under different standard deviation of speeds.
Standard deviationItemT&STSGap_T%Gap_S%σ=2Cost998.46998.461029.800.003.14Travel time2236.652236.652361.430.005.58Average speed28.5728.5728.020.001.93Average———0.003.55σ=5Cost1026.021226.81222.5219.5719.15Travel time2245.092430.232456.088.259.40Average speed28.3925.8025.309.1210.88Average———12.3113.14σ=8Cost1017.011274.391274.3925.3125.31Travel time2202.182467.532467.5312.0512.05Average speed31.1926.4426.4415.2315.23Average———17.5317.53Figure14(a) shows a comparison of the total cost of the three patterns under different speed distribution variances. In the same customer group, when the speed difference changes, the total cost of distribution of the T&S pattern changes slightly, but the total cost of T and S patterns increases drastically. The total cost deviations keep growing, Gap_T% increases from 0.00% to 25.31%, and Gap_S% increases from 3.14% to 25.31%. The change in the total cost of carbon emissions is similar to the change trend of the total cost, as indicated in Figure 14(b). With the increase in the speed difference, the carbon cost deviation Gap_T% and Gap_S% also increases. The carbon cost of the T&S pattern is always lower than those of T and S patterns, and the greater σ is, the greater the advantage of the T&S pattern is. According to Figure 14 and Table 10, compared with T and S patterns, the T&S pattern can effectively reduce the carbon emission cost and improve the timeliness of distribution.Figure 14
Cost item differences of different patterns under different standard deviation.
(a)(b)


## 5.1. Experimental Design

Nine benchmark instances are used to assess the effectiveness of HAGA_ELS (see, e.g., Augerat et al. [38], Breedam [39], and Christofides and Eilon [40]). In these instances, the geographical distribution of customers in E instances is randomly and evenly distributed, the geographical location of customers in A instances has semiclustering characteristics, and the geographical location of customers in B instances is in accordance with the clustering characteristics. The time windows are set to [60,600], which allows early arrivals but does not allow delayed arrivals. n represents the number of nodes in the instances, NV represents the number of vehicles used, Capacity represents the maximum load of vehicles, and r0 represents the ratio of the total demand to the total capacity. Table 4 lists the characteristics and basic information of the instances.Table 4
Characteristics of benchmark instances.
No.InstancesnNVCapacityr01E_n33_k433480000.9182E_n45_k445420100.8983E_n51_k55151600.9714A_n37_k63761000.9505A_n55_k95591000.9326A_n65_k96591000.9747B_n35_k53551000.8748B_n56_k75671000.8809B_n78_k1078101000.937The instance for sensitivity analysis is designed following specific distributions.nn∈10,20,40 customers are generated randomly, where the location x,y obeys the uniform distribution of U−50,50, customer demand Di obeys the uniform distribution of U0,10, and the customer time window ETi,LTi is set to (60,600). On the basis of travel speed function vt of equation (1), travel speed vijt of different arcs in different periods is randomly generated so that vijt obeys Normalvt,σ2.


## 5.2. Parameter Setting

The parameters, for instance, are set as follows: fixed vehicle costf1=110 yuan/vehicle; time penalty cost β1=50 yuan/hour and β2=100 yuan/hour; fuel cost rate f2=0.8 yuan/L; fuel consumption parameters are set as a1=0.111,a2=36.0302,anda3=0.000059 [41]; fuel-to-CO2 coefficient ρ=2.61kg/L; and carbon tax=1.5 yuan/t. The parameters for the CMEM model are given in Appendix A.The parameters to be determined for the algorithm include the cloud modelk1∼k4, the population size N, and successive generations Ω. For the parameters k1∼k4 and N, the trial and error are executed to achieve the best configuration of HAGA_ELS in a set of experiments. We test the population size of 10, 50, 100, and 200 and use crossover related probabilities k1,k3 ranging from 0 to 1 with an increment of 0.5 and mutation related probabilities k2,k4 ranging from 0 to 1 with an increment of 0.5. Thus, we have a 3D matrix holding all the possible combinations. Because finding the appropriate configuration over all the problem instances is time consuming, we opt to perform on a particular instance B_n78_k10 with maximum nodes. The algorithm has been run 10 times for each combination within fixed generations of 500. Detailed results are coming in Appendix B. The results show that when the population size increases to 200, the solution does not improve significantly in general, but the CPU time increases by 2-3 times, so the population size could be defined less than 200. We also find that in most experiments the difference between solutions in last consecutive iterations of 50 final is insignificant. From all the test runs, the best configuration has a population size of 100, crossover related probabilities of 1 (k1=k3=1), mutation related probabilities of 0.5 (k2=k4=0.5), and successive generations of 50 (for HAGA_ELS).


## 5.3. Evaluating the Performance of HAGA_ELS on TDVRPTW

To evaluate the proposed algorithm HAGA_ELS, two GA-based algorithms and two classical constructive heuristics are used for comparison. Classical constructive heuristics, sequential insertion algorithm (SI), and Clarke and Wright Saving (CWS) heuristic are, respectively, discussed in Avdoshin and Beresneva [42] and Marković et al. [33]. GA-based algorithms, standard adaptive genetic algorithm (SAGA), and elite strategy genetic algorithm (EGA) were coded as follows:(1)
SAGA is the same as HAGA_ELS, except for the elite local search described in Section4.5.(2)
EGA adopts elite local search in HAGA_ELS, but the crossover and mutation are the same as the standard genetic algorithm (SGA).Table5 lists the results obtained by the three metaheuristic algorithms. Each instance is tested on the three metaheuristic algorithms 10 times. Avg represents the average value of 10 calculation results, CV represents the deviation coefficient (CV=STD/Mean) under each algorithm. Gap_s% and Gap_e%, respectively, represent the optimized percentage of the average values of SAGA and EGA, for example, Gap_s%=SAGAAvg−HAGA_ELSAvg/HAGA_ELSAvg∗100.Table 5
Comparison of three GA-based algorithms.
InstancesHAGA_ELSSAGAEGAGap_s%Gap_e%AvgCVAvgCVAvgCVE_n33_k4863.564.22876.983.60879.912.441.551.89E_n45_k4740.763.02755.057.24753.276.231.931.69E_n51_k5579.686.32610.588.59597.287.375.333.04A_n37_k6986.344.031009.915.401004.585.832.391.85A_n55_k91125.338.141154.335.971147.965.172.582.01A_n65_k91276.165.361330.866.821304.318.674.292.21B_n35_k5962.902.79966.106.02965.393.900.330.26B_n56_k7724.173.96736.267.02730.386.571.670.86B_n78_k101344.556.321388.136.261360.617.623.241.19Average955.944.91980.916.32971.525.982.591.67As shown in Table5, the average optimization percentage compared with SAGA and EGA is 2.59% and 1.67%, respectively. The average results of the nine instances obtained by HAGA_ELS are better than those obtained by the two other algorithms. HAGA_ELS has an average deviation coefficient of 4.91% for nine instances, which is much smaller than other two algorithms compared to 6.32% of SAGA and 5.98% of EGA. The result shows that HAGA_ELS is rather effective and stable.Table6 lists the results of the three metaheuristic algorithms and the known optimal solution (OPT). The iterations represent the number of average iterations to obtain the optimal solution in 10 experiments, and Gap_h%, Gap_s%, and Gap_e%, respectively, represent the gap between the solution of three algorithms and the known optimal solution. For example, Gap_h%=HAGA_ELSBest−opt/opt∗100. As shown in Table 6, the average gap between HAGA_ELS and OPT is 3.87%, while the average gap of SAGA and EGA is 6.61% and 5.54%, respectively. Given that a time window is added to the calculation instance in this study, the deviation is still acceptable. Therefore, HAGA_ELS has firm solving capability for instances with different characteristics and sizes.Table 6
Comparison of optimal solutions of GA-based algorithms.
InstancesOPTHAGA_ELSSAGAEGAGap_h%Gap_s%Gap_e%BestIterationsBestIterationsBestIterationsE_n33_k4835853.44208871.79201874.731732.214.414.76E_n45_k4728737.18453747.40420735.444481.262.661.02E_n51_k5521575.58375601.58342590.7431110.4815.4713.39A_n37_k6949979.712021001.311681001.311343.245.515.51A_n55_k910731114.102481145.742071143.641893.836.786.58A_n65_k911741262.024171322.443291286.072707.5012.649.55B_n35_k5955957.11204959.08197961.721760.220.430.70B_n56_k7707715.06295728.36289721.132801.143.022.00B_n78_k1012661328.635251373.953811346.383644.958.536.35Average912946.98325972.41282962.352603.876.615.54Figure10 shows the convergence graph of the three metaheuristic algorithms for large instances. Figures 10(a)–10(c) response to E_n51_k5, A_n65_k9, and B_n78_k10 instances, respectively. According to the graph, HAGA_ELS can obtain more qualified solution within more iterations.Figure 10
Convergence comparison.
(a)(b)(c)The computational results of HAGA_ELS and classical constructive heuristics are shown in Table7. Gap_h% is much smaller than Gap_si% and Gap_cws%, which means HAGA_ELS performs better than SI and CWS in any instances. The average gap of HAGA_ELS is 3.87%, while the average gaps of SI and CWS are all over 15%. Thus, compared with constructive heuristics, HAGA_ELS has ability to obtain the better solution.Table 7
Comparison of optimal solutions of HAGA_ELS and constructive heuristics.
InstancesOPTHAGA-ELSSICWSGap_h%Gap_si%Gap_cws%E_n33_k4835853.44888.82937.772.216.4512.31E_n45_k4728737.18881.00964.601.2621.0232.50E_n51_k5521575.58638.28771.8510.4822.5148.15A_n37_k6949979.711175.901065.503.2423.9112.28A_n55_k910731114.101394.101257.703.8329.9317.21A_n65_k911741262.021571.701639.607.5033.8839.66B_n35_k5955957.111007.501245.100.225.5030.38B_n56_k7707715.06869.39861.501.1422.9721.85B_n78_k1012661328.631424.401704.204.9512.5134.61Average912946.981094.571160.873.8719.8527.66


## 5.4. Comparison of Three Patterns

The second experiment is a small-size experiment with 10 customers to verify the validity of the model, and it is carried out under three patterns (T&S, T, and S). Each experiment is run 10 times. Table8 shows a cost comparison for the three patterns. Given that the fixed costs are the same, they are omitted in the Table 8, in which Gap_T% and Gap_S%, respectively, represent the difference between T&S pattern and T and or S patterns. The optimization percentage of each kind of costs between S pattern and T&S pattern is calculated as Gap%=cost′−cost/cost. Figure 11 shows the distribution scheme under the three patterns with the corresponding carbon emission costs and average travel speed. Figure 12 shows the travel speed of each road section under the three patterns.Table 8
Cost item deviations of the different patterns at a small-size instance.
Cost itemT&STSGap_T%Gap_S%Total cost670.74741.21851.5410.5126.96Transportation cost203.66197.52191.27−3.01−6.08Carbon cost170.58182.19199.936.8117.21Time penalty cost76.49141.5240.3484.99214.21Figure 11
Comparison of paths under three patterns.Figure 12
Speeds of small-size instance under three patterns.Table8 shows the comparison of costs for the three patterns. Although the transportation cost under the T&S pattern is slightly higher, the total cost of T and S patterns is, respectively, 10.51% and 26.96% higher than that of the T&S pattern, mainly because the time penalty cost is greatly reduced under the T&S pattern. This finding shows that applying the spatiotemporal-varying speed facilitates to reduce the default of the customer required.Figure11 shows the distribution scheme under the three patterns, and Figure 12 shows the average travel speed of each arc under the three patterns. As shown in Figure 11, the scheme obtained by the T&S pattern has the lowest carbon emission cost and the fastest average travel speed. The average travel speeds under T&S, T, and S patterns are 26.77, 24.38, and 23.17 km/h, respectively. Figure 12 indicates that under the T&S pattern, the travel routing can avoid traffic congestion, so that the distribution has a high travel speed. The average travel speed of arcs can be seen in Figure 11. The maximum travel speed under the three patterns does not reach the inflection point of the carbon emission cost function (Equation (9)). Hence, the carbon emission cost drops with the increase in speed. The T&S pattern improves travel speed by avoiding congestion, thus reducing the carbon emission cost.As shown in Figure12, the T&S pattern mainly selects the route to avoid congestion through the node sequence. Thus, its travel speed is generally at a high level. Detailed analysis of Figure 11 indicates that the T&S pattern avoids the congestion in the fourth arc (a⟶c) and the fifth arc (c⟶b). In these arcs, the order of node arrival is c⟶a⟶b and b⟶a⟶c in T and S patterns, respectively, and the order of node arrival is changed to a⟶c⟶b in the T&S pattern. This finding shows that the T&S pattern can help vehicles plan the node arrival sequence and avoid the traffic congestion, thus reducing carbon emission cost and improving distribution speed.


## 5.5. Comparison of Different Factor Impacts on Three Patterns

To assess the effectiveness of the T&S pattern further, we simulate and analyze different customer sizes in Section5.5.1 and different speed distributions in Section 5.5.2.


### 5.5.1. Analysis of Different Customer Size

We select three customer groups with different sizes (n=10, n=20, and n=40) for testing and adopt three dynamic speed patterns, namely, T&S, T, and S. Nine groups of experiments are performed, and each group of experiments is run 10 times. We obtain the best results of the three dynamic speed patterns with different sizes and compare them in terms of cost, travel time, and speed. Under σ=5, we generate speed data that obey a normal distribution under different arcs in the same period to ensure a random difference of speed in road network. At the same time, the travel speed on each arc obeys the function law of the period, which reflects the real traffic network to a certain extent.We report the computational results of the nine groups of experiments in Table9, in which Gap_T%=Item_T−Item_T&S/Item_T&S×100, Gap_S%=Item_S−Item_T&S/Item_T&S×100. Gap_T% and Gap_S% are greater than zero under different customer sizes n; that is, the best result is obtained by using the spatiotemporal travel speed pattern (T&S). The Gap_T% in cost is over 10%, and that of Gap_S% is over 15%. The routing optimization under the T&S pattern is superior to those under T and S patterns in terms of cost.Table 9
Numerical results of different patterns under different sizes.
Customer sizeItemT&STSGap_T%Gap_S%n=10Cost670.74741.21851.5410.5126.96Travel time1168.381198.301250.512.567.03Average speed26.7724.3823.178.9313.45Average———7.3315.81n=20Cost1026.021226.81222.5219.5719.15Travel time2245.092430.232456.088.259.40Average speed28.3925.8025.309.1210.88Average———12.3113.14n=40Cost2730.463254.573650.3419.1933.69Travel time5332.415813.85855.059.039.80Average speed31.0327.7627.4110.5411.67Average———12.9218.39Figure13(a) shows the total cost deviation between T&S pattern and T (or S) pattern under different customer sizes, and Figure 13(b) presents the total travel time deviation between T&S pattern and T (or S) pattern. With the increase in customer size, the total cost deviation and total travel time deviation also increase. Table 9 also indicates that the average optimization percentage Gap% appears an upward trend with customers increasing. Gap_T% increases from 7.33% to 12.92%, and Gap_S% increases from 15.81% to 18.39%. Therefore, with the increase in customer size, the T&S pattern performs much better than T and S patterns by reducing the distribution cost and improving the distribution efficiency.Figure 13
Total cost and total time deviation of different patterns under different sizes.
(a)(b)


### 5.5.2. Analysis of Different Speed Deviation

Three different standard deviations of speed distribution (σ=10, σ=20, and σ=40) are selected to randomly generate a spatiotemporal speed matrix. The travel speed vijl of different arcs in different periods obeys the normal distribution of Normalvt,σ2. (vt is the time travel speed function (equation (1)). By changing the variance of the normal distribution, the difference of travel speed in arcs can be adjusted. The greater the variance is, the greater the difference in travel speed is among different arcs. In the same manner, three different spatiotemporal-varying speed matrices are tested with the three dynamic speed patterns. To avoid the compound effect caused by the change in customer group size, we set the customer group size to be constant (n=20), and each group of experiments is run 10 times to obtain the optimal results of three speed patterns under different speed matrices. A comparative analysis is performed in terms of cost, travel time, and speed.The numerical results of each group of experiments are shown in Table10, which illustrates that Gap_T% and Gap_S% are greater than or equal to 0 under the velocity matrix generated by different variances σ2. Thus, the spatiotemporal speed pattern (T&S) is the best method under the velocity matrix generated by different variances. When the variance of the distribution is small, that is, σ=2, the difference in spatial speed is small. In this case, the running results of the T&S and T patterns are consistent, so there is no obvious difference whether speed spatiality is considered or not. However, with the increase in speed distribution variance σ2, Gap_T% increases continuously and reaches at 17.53% under σ=8. Meanwhile, Gap_S% increases from 3.55% to 17.53%. When the difference of spatial speed amplifies, the T&S pattern becomes superior to T and S patterns, and the deficiency of the T pattern that does not consider speed spatiality becomes increasingly apparent.Table 10
Numerical results of different patterns under different standard deviation of speeds.
Standard deviationItemT&STSGap_T%Gap_S%σ=2Cost998.46998.461029.800.003.14Travel time2236.652236.652361.430.005.58Average speed28.5728.5728.020.001.93Average———0.003.55σ=5Cost1026.021226.81222.5219.5719.15Travel time2245.092430.232456.088.259.40Average speed28.3925.8025.309.1210.88Average———12.3113.14σ=8Cost1017.011274.391274.3925.3125.31Travel time2202.182467.532467.5312.0512.05Average speed31.1926.4426.4415.2315.23Average———17.5317.53Figure14(a) shows a comparison of the total cost of the three patterns under different speed distribution variances. In the same customer group, when the speed difference changes, the total cost of distribution of the T&S pattern changes slightly, but the total cost of T and S patterns increases drastically. The total cost deviations keep growing, Gap_T% increases from 0.00% to 25.31%, and Gap_S% increases from 3.14% to 25.31%. The change in the total cost of carbon emissions is similar to the change trend of the total cost, as indicated in Figure 14(b). With the increase in the speed difference, the carbon cost deviation Gap_T% and Gap_S% also increases. The carbon cost of the T&S pattern is always lower than those of T and S patterns, and the greater σ is, the greater the advantage of the T&S pattern is. According to Figure 14 and Table 10, compared with T and S patterns, the T&S pattern can effectively reduce the carbon emission cost and improve the timeliness of distribution.Figure 14
Cost item differences of different patterns under different standard deviation.
(a)(b)


## 5.5.1. Analysis of Different Customer Size

We select three customer groups with different sizes (n=10, n=20, and n=40) for testing and adopt three dynamic speed patterns, namely, T&S, T, and S. Nine groups of experiments are performed, and each group of experiments is run 10 times. We obtain the best results of the three dynamic speed patterns with different sizes and compare them in terms of cost, travel time, and speed. Under σ=5, we generate speed data that obey a normal distribution under different arcs in the same period to ensure a random difference of speed in road network. At the same time, the travel speed on each arc obeys the function law of the period, which reflects the real traffic network to a certain extent.We report the computational results of the nine groups of experiments in Table9, in which Gap_T%=Item_T−Item_T&S/Item_T&S×100, Gap_S%=Item_S−Item_T&S/Item_T&S×100. Gap_T% and Gap_S% are greater than zero under different customer sizes n; that is, the best result is obtained by using the spatiotemporal travel speed pattern (T&S). The Gap_T% in cost is over 10%, and that of Gap_S% is over 15%. The routing optimization under the T&S pattern is superior to those under T and S patterns in terms of cost.Table 9
Numerical results of different patterns under different sizes.
Customer sizeItemT&STSGap_T%Gap_S%n=10Cost670.74741.21851.5410.5126.96Travel time1168.381198.301250.512.567.03Average speed26.7724.3823.178.9313.45Average———7.3315.81n=20Cost1026.021226.81222.5219.5719.15Travel time2245.092430.232456.088.259.40Average speed28.3925.8025.309.1210.88Average———12.3113.14n=40Cost2730.463254.573650.3419.1933.69Travel time5332.415813.85855.059.039.80Average speed31.0327.7627.4110.5411.67Average———12.9218.39Figure13(a) shows the total cost deviation between T&S pattern and T (or S) pattern under different customer sizes, and Figure 13(b) presents the total travel time deviation between T&S pattern and T (or S) pattern. With the increase in customer size, the total cost deviation and total travel time deviation also increase. Table 9 also indicates that the average optimization percentage Gap% appears an upward trend with customers increasing. Gap_T% increases from 7.33% to 12.92%, and Gap_S% increases from 15.81% to 18.39%. Therefore, with the increase in customer size, the T&S pattern performs much better than T and S patterns by reducing the distribution cost and improving the distribution efficiency.Figure 13
Total cost and total time deviation of different patterns under different sizes.
(a)(b)


## 5.5.2. Analysis of Different Speed Deviation

Three different standard deviations of speed distribution (σ=10, σ=20, and σ=40) are selected to randomly generate a spatiotemporal speed matrix. The travel speed vijl of different arcs in different periods obeys the normal distribution of Normalvt,σ2. (vt is the time travel speed function (equation (1)). By changing the variance of the normal distribution, the difference of travel speed in arcs can be adjusted. The greater the variance is, the greater the difference in travel speed is among different arcs. In the same manner, three different spatiotemporal-varying speed matrices are tested with the three dynamic speed patterns. To avoid the compound effect caused by the change in customer group size, we set the customer group size to be constant (n=20), and each group of experiments is run 10 times to obtain the optimal results of three speed patterns under different speed matrices. A comparative analysis is performed in terms of cost, travel time, and speed.The numerical results of each group of experiments are shown in Table10, which illustrates that Gap_T% and Gap_S% are greater than or equal to 0 under the velocity matrix generated by different variances σ2. Thus, the spatiotemporal speed pattern (T&S) is the best method under the velocity matrix generated by different variances. When the variance of the distribution is small, that is, σ=2, the difference in spatial speed is small. In this case, the running results of the T&S and T patterns are consistent, so there is no obvious difference whether speed spatiality is considered or not. However, with the increase in speed distribution variance σ2, Gap_T% increases continuously and reaches at 17.53% under σ=8. Meanwhile, Gap_S% increases from 3.55% to 17.53%. When the difference of spatial speed amplifies, the T&S pattern becomes superior to T and S patterns, and the deficiency of the T pattern that does not consider speed spatiality becomes increasingly apparent.Table 10
Numerical results of different patterns under different standard deviation of speeds.
Standard deviationItemT&STSGap_T%Gap_S%σ=2Cost998.46998.461029.800.003.14Travel time2236.652236.652361.430.005.58Average speed28.5728.5728.020.001.93Average———0.003.55σ=5Cost1026.021226.81222.5219.5719.15Travel time2245.092430.232456.088.259.40Average speed28.3925.8025.309.1210.88Average———12.3113.14σ=8Cost1017.011274.391274.3925.3125.31Travel time2202.182467.532467.5312.0512.05Average speed31.1926.4426.4415.2315.23Average———17.5317.53Figure14(a) shows a comparison of the total cost of the three patterns under different speed distribution variances. In the same customer group, when the speed difference changes, the total cost of distribution of the T&S pattern changes slightly, but the total cost of T and S patterns increases drastically. The total cost deviations keep growing, Gap_T% increases from 0.00% to 25.31%, and Gap_S% increases from 3.14% to 25.31%. The change in the total cost of carbon emissions is similar to the change trend of the total cost, as indicated in Figure 14(b). With the increase in the speed difference, the carbon cost deviation Gap_T% and Gap_S% also increases. The carbon cost of the T&S pattern is always lower than those of T and S patterns, and the greater σ is, the greater the advantage of the T&S pattern is. According to Figure 14 and Table 10, compared with T and S patterns, the T&S pattern can effectively reduce the carbon emission cost and improve the timeliness of distribution.Figure 14
Cost item differences of different patterns under different standard deviation.
(a)(b)


## 6. Conclusion

Complex traffic networks and ever-changing traffic flow in cities increase the uncertainty in urban logistics distribution businesses. Route optimization considering traffic changes has become an important measure to reduce the cost and carbon emission of freight enterprises. In this study, we establish a mixed-integer nonlinear programming model that considers carbon emission and propose a spatiotemporal-dependent travel time calculation method. A hybrid adaptive genetic algorithm is developed, in which an elite neighborhood search operator is applied to ensure the obtainment of the best solution and to increase individual diversity, thus avoiding falling into the local optimum. The effectiveness of this algorithm is evaluated by comparing it with other algorithms on benchmark instances. A sensitivity analysis based on different sizes and speed variances assesses the effectiveness of this pattern. The following conclusions are drawn:(1)
Compared with previous algorithms, HAGA_ELS can obtain a high-quality solution. It can avoid falling into the local optimal solution and the premature problem.(2)
Compared with those patterns that only consider time or space dependence of speed, it will obtain more reasonable routings that consider spatiotemporal changes in speed. The rationality specifically reflects in less time penalty cost and carbon emission cost.(3)
The customer size and speed differences have a significant influence on the effect of spatiotemporal pattern. With the increase of customer size or speed difference, this pattern will bring lower total cost and shorter travel time compared with other patterns that only consider time or space dependence of speed.Based on the comparison of three patterns, we give the following recommendations to freight enterprises: freight enterprises should make use of the intelligent transportation system to scheme preoptimized route, which provides the possibility of obtaining traffic information about the road network. Especially in the case of scattered customers resulting in considerable speed difference, preoptimized scheme under spatiotemporal traffic information obtains performance gain.The future research will focus on speeding up the algorithm. For further research, one direction is to speed up the algorithm to solve the large size problem more efficiently.



---

*Source: 9729784-2021-09-24.xml*