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.
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. 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.
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. Table 1 lists the methods of time-dependent issue processing of the above literature studies.
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.
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. 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 graph , the set of nodes is , where 0 is the distribution center and is the set of known customers. The set of edges is . The location , nonnegative demand , service time , and delivery time of the arrival window of customer are known. The shortest path between two nodes, namely, edge , has a fixed distance of . The distribution road network has a time-varying traffic state, and the distribution time range is divided into a fixed number of time intervals . The different arcs in each period have various travel speed , and travel speed is known in advance in accordance with the travel speed function. Distribution center has a certain number of distribution vehicles , and all customers go back to the distribution center after the completion of the distribution. Distribution vehicle has the same maximum load . The fuel consumption formula is given based on travel speed , load weight , and travel distance , 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 speed and time in an entire day proposed by Xu et al. [20] can be approximated by the trigonometric function in equation (1), where , , and 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). 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 , the time of each period is recorded as , and the time interval is . The travel speed of the vehicle in any period is
In some studies (e.g., Kwon et al. [23]), the travel time calculation of arc depends on the departure time of vehicles from node 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 , the traveling speed of different arc 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 but also relate to the arc from node to the next node . 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 of any arc in period obeys normal distribution , and the variance reflects the difference degree of travel speeds of different arcs. On this basis, the nonnegative random velocity matrix 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 . 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 to node at time across periods, so road section length is divided into . The speed in the segment is , and the time of division is , as shown in Figure 4.
The travel time of road section is , and the travel time of subsequent road section is . Thus, an equation between section length and time is available:
The driving distance of the S + 1 section is
Therefore, the arrival time of node is
According to equation (5), the time when the vehicle arrives at node is an increasing function of time when the vehicle departs from node , which meets the FIFO criterion.
When the arc spans periods, the following conditions should be met when the vehicle starts from node : the vehicle starts from the first period and has not finished driving in arc from the first period to the th period. The corresponding expression is as follows:(a),(b)
The abovementioned arrangement yields
The calculation formula of the arrival time of the next node where any arc crosses periods is as follows:
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:
In equation (8), , is the fuel-air mass ratio, is a typical diesel calorific value, is the fuel conversion coefficient, engine module coefficient , speed module coefficient , and vehicle load module coefficient . 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
In equation (9), is a binary variable indicating whether arc is traveled by vehicle () or not (). is the fuel consumption of the vehicle from node to node . is the length of the arc during time . represents the travel speed of the arc during time . And represents the loading weight of the vehicle from node to node .
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.
The specific process of calculating carbon emission in consideration of spatiotemporal-dependent travel speed is shown in Algorithm 1.
|
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 Table 3.
In the above decision variables, is the variable related to the arc of vehicle routing. represents the connection between vehicles and customers. and are variables in vehicle driving process. indicates the connection between vehicle load and arcs. 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:
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 is a positive big number and 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.
4.1. Initial Solution
In this study, population size is set as , and each chromosome is composed of ( 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 .
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 is a constant greater than 1, is the fitness of chromosome , and is the objective function value corresponding to chromosome :
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, 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 individuals from their parents. The selection probability of each chromosome is , 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 .
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 and mutation rate is as follows:where is the larger fitness value among the crossover individuals, is the fitness value of the mutant individual, is the maximum fitness value in the population, and is the average fitness value. is a normal random number generator. are the control coefficients (, ), and are constants in the range of . The control coefficient in the algorithm is appropriately adjusted in accordance with evolutional generation, 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.
(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.
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 Figure 9(a), the positions of nodes 6 and 5 are selected, and the intermediate part is sequentially moved forward, i.e., .(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 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 in a certain subrouting are transformed into , as shown in Figure 9(d).
| ||||||||||||||||||||||||||||||||||||||||||||||
(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 Section 5.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. represents the number of nodes in the instances, represents the number of vehicles used, represents the maximum load of vehicles, and represents the ratio of the total demand to the total capacity. Table 4 lists the characteristics and basic information of the instances.
The instance for sensitivity analysis is designed following specific distributions. customers are generated randomly, where the location obeys the uniform distribution of , customer demand obeys the uniform distribution of , and the customer time window is set to (60,600). On the basis of travel speed function of equation (1), travel speed of different arcs in different periods is randomly generated so that obeys .
5.2. Parameter Setting
The parameters, for instance, are set as follows: fixed vehicle cost yuan/vehicle; time penalty cost yuan/hour and yuan/hour; fuel cost rate yuan/L; fuel consumption parameters are set as [41]; fuel-to-CO2 coefficient ; and carbon yuan/t. The parameters for the CMEM model are given in Appendix A.
The parameters to be determined for the algorithm include the cloud model , the population size , and successive generations . For the parameters and , 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 ranging from 0 to 1 with an increment of 0.5 and mutation related probabilities 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 (), mutation related probabilities of 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 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 () under each algorithm. and , respectively, represent the optimized percentage of the average values of SAGA and EGA, for example, .
As shown in Table 5, 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.
Table 6 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 , , and , respectively, represent the gap between the solution of three algorithms and the known optimal solution. For example, . 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.
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 instances, respectively. According to the graph, HAGA_ELS can obtain more qualified solution within more iterations.
(a)
(b)
(c)
The computational results of HAGA_ELS and classical constructive heuristics are shown in Table 7. is much smaller than and , 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.
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. 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 and , 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 . 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. This 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. 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 Figure 12, 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 Section 5.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 (, , and ) 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 , 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 , . and are greater than zero under different customer sizes ; that is, the best result is obtained by using the spatiotemporal travel speed pattern (T&S). The in cost is over 10%, and that of is over 15%. The 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 appears an upward trend with customers increasing. increases from 7.33% to 12.92%, and 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.
(a)
(b)
5.5.2. Analysis of Different Speed Deviation
Three different standard deviations of speed distribution (, , and ) are selected to randomly generate a spatiotemporal speed matrix. The travel speed of different arcs in different periods obeys the normal distribution of . ( 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 (), 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 Table 10, which illustrates that and are greater than or equal to 0 under the velocity matrix generated by different variances . 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, , 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 , increases continuously and reaches at 17.53% under . Meanwhile, 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. The total cost deviations keep growing, increases from 0.00% to 25.31%, and 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 and 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.
(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.
Appendix
(A). CMEM Model Parameters in Table 11
(B). Algorithm Performance for Various Population Sizes
In Table 12, the obtained average optimal solution values and the computational time for considered parameter combinations are presented. Parameter combinations consists of the population size and the cloud model settings or .
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Authors’ Contributions
Ziqi Liu and Yeping Chen contributed to the work equally.
Acknowledgments
This work was supported by the National Key Research and Development Program of China (2018YFB1403100) and the Humanities and Social Sciences Youth Foundation of Ministry of Education of China (19YJC630079).