Spatiotemporal-Dependent Vehicle Routing Problem Considering Carbon Emissions

Aiming to improve the timeliness of logistics distribution and render the optimized route scheme eﬀective under the real traﬃc 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 intensiﬁcation 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 eﬀective 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 out-performs them in terms of carbon emissions and timeliness. The research can provide a scientiﬁc and reasonable method for logistics enterprises to plan the vehicle schedule focusing on spatiotemporal-dependent speed of the road network.


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.
us, the objective of distribution should additionally take carbon emissions influenced by travel speed on arcs into account. erefore, 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]. e 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. erefore, considering the minimum distance or the minimum distribution cost cannot reduce carbon emissions. e 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]. erefore, this study investigates the spatiotemporal-dependent green vehicle routing problem with time windows. e 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. us, we apply the genetic algorithm, a classical metaheuristic algorithm to find the near-optimal solutions of the problem. e 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 firstout rules. We adopt the fuel consumption introduced by Esteves-Booth et al. [7] with the spatiotemporal speed considered. is work establishes a mixedinteger 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.
e rest of the article is as follows. Section 2 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.

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. erefore, 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).

Research on TDVRPTW.
e 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]. us, 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. e first type of research on time-dependent travel time directly assumes that travel time is a function related to the time in a day. e 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. e 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. e 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. Table 1 lists the methods of time-dependent issue processing of the above literature studies.

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. Table 2 summarizes and compares the formulation in the stream of TDGVRPTW.
From 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. e travel speed in different dimensions obeys a certain distribution. is study analyzes the limitations of the time-dependent speed division pattern (T) and the space-dependent speed division pattern (S).

Problem and Mathematical Model
TDGVRPTW studied in this paper can be described as follows. In the network graph G � (V, E), the set of nodes is where 0 is the distribution center and V 0 � 1, ..., n { } is the set of known customers. e set of edges is e location (x i , y i ), nonnegative demand w i , service time ST i , and delivery time of the arrival window [ET i , LT i ] of customer i are known. e shortest path between two nodes, namely, edge (i, j) ∈ E, has a fixed distance of d ij . e distribution road network has a timevarying traffic state, and the distribution time range is divided into a fixed number of time intervals L � 1, 2, . . . , h { }. e different arcs (i, j) ∈ E in each period l ∈ L have various travel speed v l ij , and travel speed v l ij 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 V 0 go back to the distribution center after the completion of the distribution. Distribution vehicle k ∈ K has the same maximum load Q. e fuel consumption formula FUEL(v, q, d) is given based on travel speed v, load weight q, and travel distance d, and fuel-to-CO 2 coefficient ρ is known. On the premise of the constraints of vehicle loading and customer time window, Stepwise function

Continuous function
Historical data Sun 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 work Temporal axis: travel speed varies with a stepwise function converted from continuous function Spatial axis: travel speed obeys a distribution in accordance with historical data the vehicle distribution routing is optimized considering the carbon emission of the vehicle. us, 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) e 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.

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 speed v 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 , c p , and δ p are constants related to the road conditions of the day: 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 ( Figure 1). e finer the time division is, the closer the speed segmentation function is to the continuous change in speed. e time interval is divided into periods L � 1, 2, . . . , h { }, the time of each period is recorded as Γ l (l ∈ L), and the time interval is ΔI. e travel speed of the vehicle in any period l is In some studies (e.g., Kwon et al. [23]), the travel time calculation of arc (i, j) depends 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. e allday 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 v l ij of different arc (i, j) is various. erefore, the travel speed and travel time in the Discrete Dynamics in Nature and Society 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 v l ij of any arc (i, j) ∈ E in period l obeys normal distribution N(v l , σ 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.
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 arc (i, 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. erefore, 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. e vehicle is assumed to travel from node i to node j at time a i across S periods, so e driving distance of the S + 1 section is According to equation (5), the time t j when the vehicle arrives at node j is an increasing function of time a i when the vehicle departs from node i, which meets the FIFO criterion.
When the arc (i, 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. e corresponding expression is as follows: e calculation formula of the arrival time t j of the next node j where any arc (i, j) crosses S(S < h) periods is as follows:  Discrete Dynamics in Nature and Society

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]. e objective function of fuel consumption when driving in different sections is as follows: In equation (8), λ � ζ/κφ, ζ is the fuel-air mass ratio, κ is a typical diesel calorific value, φ is the fuel conversion coefficient, engine module coefficient e 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 In equation (9), x ijk is a binary variable indicating whether arc (i, j) is traveled by vehicle k (x ijk � 1) or not (x ijk � 0). fuel ij is the fuel consumption of the vehicle from node i to node j. d l ij is the length of the arc (i, j) during time l. v l ij represents the travel speed of the arc (i, j) during time l. And Q ij represents the loading weight of the vehicle from node i to node j. Figure 5 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. e specific process of calculating carbon emission in consideration of spatiotemporal-dependent travel speed is shown in Algorithm 1.

Mathematical Model.
On the basis of the analysis above, the TDGVRPTW model is established, and the model parameters and decision variables are shown in Table 3.
In the above decision variables, x ijk is the variable related to the arc of vehicle routing. y ki represents the connection between vehicles and customers. Q ij and t i are variables in vehicle driving process. Q ij indicates the connection between vehicle load and arcs. t i is calculated by  Discrete Dynamics in Nature and Society spatiotemporal-dependent speed, which is the crucial variable considered for TDGVRPTW. e mathematical optimization model of this paper is formulated to minimize the total distribution costs, as follows: i∈V 0  //Depart from i to j, assume that the period is divided into L segments //Judging whether to span S periods (3) If the criteria of crossing S periods is satisfied. //i.e. e path is divided into S + 1 segments d′←d′ − d k ; (8) End for (9) d S+1 ←d′; //the length of the S + 1 section //Calculate the arrival time at node j and fuel consumption of traveling arc(i, j) i∈V 0 i∈V 0 j∈V 0 j∈V 0 Set of nodes excluding the depot, Travel time spent on arc (i, j) Notations for fuel consumption and CO 2 emission α 1 Engine power module coefficient α 2 Engine power speed module coefficient α 3 Vehicle load module coefficient λ Coefficient of fuel consumption ρ Fuel-to-CO 2 coefficient Notations for total costs f 1 Fixed cost of vehicle f 2 Drive cost rate β 1 Penalty coefficient of arrival before the time window β 2 Penalty coefficient of arrival after the time window tax Carbon tax e binary and continuous decision variables of the model are defined as follows: 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 t ij 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.

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. e 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. erefore, this study designs a hybrid adaptive genetic algorithm with elite local search (HAGA_ELS). e 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. e 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. e operation flow of HAGA_ELS is shown in Figure 6.

Initial Solution.
In this study, population size is set as N, and each chromosome is composed of n (n is the number of customers) integers. e 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. e 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. e 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. e steps above are repeated until the chromosome number reaches population size N.

Fitness Function.
e 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].
e fitness function is shown in equation (20), where a is a constant greater than 1, f i is the fitness of chromosome i, and z i is the objective function value corresponding to chromosome i:

Selection Operator.
In this study, elite and roulette strategies are combined to implement the selection operation. e specific steps are as follows: with the descending order of chromosome fitness value, N(1 − Φ) 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. en, the roulette strategy is adopted to select NΦ individuals from their parents. e selection probability of each chromosome is p i � f i /( i∈N f i ), and the greater the fitness is, the greater the probability of being selected is. e 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.

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. erefore, 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. e algorithm for generating crossover rate P c and mutation rate P m is as follows: Discrete Dynamics in Nature and Society where f is the larger fitness value among the crossover individuals, f ′ is the fitness value of the mutant individual, f max is the maximum fitness value in the population, and f avg is the average fitness value. Norm is a normal random number generator. C 1 ∼ C 4 are the control coefficients

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. e 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]. e 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. e crossover process is shown in Figure 7.

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.
e new chromosome after mutation is retained, and the original chromosome is deleted so that the individual is always mutated in an improved direction. e specific mutation operation is shown in Figure 8.

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.
ese 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 Figure 9(a), the positions of nodes 6 and 5 are selected, and the intermediate part 6 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. e exchange operation is illustrated in Figure 9(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 Figure 9(c). (4) 3-opt exchange: four continuous points (i, i + 1, j, j + 1) in a certain subrouting are transformed into (i + 1, i, j + 1, j), as shown in Figure 9(d).

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 Ω.

Computational Experiment and Analyses
In this section, the experimental design is described in Section 5.1. Parameter setting is provided in Section 5.2. e 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.

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. e 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 r 0 represents the ratio of the total demand to the total capacity. Table 4 lists the characteristics and basic information of the instances. e instance for sensitivity analysis is designed following specific distributions. n(n ∈ 10, 20, 40 { }) customers are generated randomly, where the location (x, y) obeys the uniform distribution of U(− 50, 50), customer demand D i obeys the uniform distribution of U(0, 10), and the customer time window (ET i , LT i ) is set to (60,600). On the basis of travel speed function v(t) of equation (1), travel speed v t ij of different arcs in different periods is randomly generated so that v t ij obeys Normal(v(t), σ 2 ).

Parameter
Setting. e parameters, for instance, are set as follows: fixed vehicle cost f 1 � 110 yuan/vehicle; time penalty cost β 1 � 50 yuan/hour and β 2 � 100 yuan/hour; fuel Discrete Dynamics in Nature and Society cost rate f 2 � 0.8 yuan/L; fuel consumption parameters are set as a 1 � 0.111, a 2 � 36.0302, and a 3 � 0.000059 [41]; fuelto-CO 2 coefficient ρ � 2.61 kg/L; and carbon tax � 1.5 yuan/ t. e parameters for the CMEM model are given in Appendix A. e parameters to be determined for the algorithm include the cloud model k 1 ∼ k 4 , the population size N, and successive generations Ω. For the parameters k 1 ∼ k 4 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 k 1 , k 3 ranging from 0 to 1 with an increment of 0.5 and mutation related probabilities k 2 , k 4 ranging from 0 to 1 with an increment of 0.5. us, 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. e algorithm has been run 10 times for each combination within fixed generations of 500. Detailed results are coming in Appendix B. e 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 (k 1 � k 3 � 1), mutation related probabilities of 0.5 (k 2 � k 4 � 0.5), and successive generations of 50 (for HAGA_ELS).

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 Section 4.5. (2) EGA adopts elite local search in HAGA_ELS, but the crossover and mutation are the same as the standard genetic algorithm (SGA). Table 5 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 � ST D/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% � (SAGA Avg − HAGA ELS Avg )/ HAGA ELS Avg * 100.
As shown in Table 5, the average optimization percentage compared with SAGA and EGA is 2.59% and 1.67%, respectively. e 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. e result shows that HAGA_ELS is rather effective and stable. Table 6 lists the results of the three metaheuristic algorithms and the known optimal solution (OPT). e 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 ELS Best − 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. erefore, HAGA_ELS has firm solving capability for instances with different characteristics and sizes. Figure 10 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 (1)Input XGood, fitness (2)For i � 1 to size (XGood,1) do begin (3) XGood' ←XGood(i, : ) //Update XGood' after each step of operation //Cyclic operator. Choosing r 1 and r 2 at random, which are in the same subrouting (4) New XGoo d(r 2 )←XGoo d'(r 1 ) (5) For k � r 1 to r 2 − 1 do begin (6) New XGoo d(k)←XGoo d'(k + 1) (7) End for //1-0 swap. Choosing r 1 and r 2 in different subroutings (8) Remove New XGoo d(i, r 1 ) (9) New XGoo d(i, r 2 )←XGoo d'(i, r 1 ); New XGoo d(i, r 2 + 1: end)←XGoo d'(i, r 2 : end) //2-opt swap. Choosing r 1 and r 2 in different subroutings (10) New XGoo d(i, r 1 : where node is 0)←XGoo d'(i, r 2 where node is 0)) (11) Similarly for the other subrouting where r 2 is located //3-opt swap. Choosing r 1 and r 2 in the same subrouting Discrete Dynamics in Nature and Society    instances, respectively. According to the graph, HAGA_ELS can obtain more qualified solution within more iterations. e computational results of HAGA_ELS and classical constructive heuristics are shown in Table 7. Gap h% is much smaller than Gap si% and Gap cws%, which means HAGA_ELS performs better than SI and CWS in any instances. e average gap of HAGA_ELS is 3.87%, while the average gaps of SI and CWS are all over 15%. us, compared with constructive heuristics, HAGA_ELS has ability to obtain the better solution.

Comparison of ree Patterns.
e 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. Table 8 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. e 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 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.
is finding shows that applying the spatiotemporal-varying speed facilitates to reduce the default of the customer required. Figure 11 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. e 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. e average travel speed of arcs can be seen in Figure 11. e 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. e T&S pattern improves travel speed by avoiding congestion, thus reducing the carbon emission cost.
As shown in Figure 12, the T&S pattern mainly selects the route to avoid congestion through the node sequence.
us, 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. is 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.

Comparison of Different Factor Impacts on ree Patterns.
To assess the effectiveness of the T&S pattern further, we simulate and analyze different customer sizes in Section 5.5.1 and different speed distributions in Section 5.5.2.

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 Table 9, 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). e Gap T% in cost is over 10%, and that of Gap S% is over 15%. e routing optimization under the T&S pattern is superior to those under T and S patterns in terms of cost. Figure 13(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%. erefore, 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.

Analysis of Different Speed Deviation.
ree different standard deviations of speed distribution (σ � 10, σ � 20, and σ � 40) are selected to randomly generate a spatiotemporal speed matrix. e travel speed v l ij of different arcs in different periods obeys the normal distribution of Normal(v(t), σ 2 ). (v(t) 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. e 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. e numerical results of each group of experiments are shown in Table 10, which illustrates that Gap T% and Gap S% are greater than or equal to 0 under the velocity matrix generated by different variances σ 2 .
us, 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. Figure 14(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. e total cost deviations keep growing, Gap T% increases from 0.00% to 25.31%, and Gap S% increases from 3.14% to 25.31%. e 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. e 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.

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. e 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. e 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. e rationality specifically reflects in less time penalty cost and carbon emission cost. (3) e 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.
e 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.