Vehicle path optimization method of multi-objective ant colony robust optimization algorithm driven by time-space attenuation factors

Abstract

Aiming at the problem of time-dependent multi-target green vehicle path optimization, the invention designs a space-time attenuation factor-driven urban logistics low-carbon robust optimization method, takes the total vehicle distribution time and carbon emission as optimization targets, and adopts double-layer ant colony pheromones to guide ant colony search, so that the optimal path optimization is realized. The adaptive capacity of ants in a high-disturbance uncertain environment is improved, the convergence of solutions is enhanced, more solution sets balancing robustness and optimality are excavated through solution robustness evaluation and a feedback mechanism thereof, the diversity of optimal solutions is improved, an optimal vehicle path robust optimization scheme is obtained, the carbon emission of urban road network distribution is reduced, and the urban road network distribution efficiency is improved. And the distribution efficiency of logistics is improved.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of robust optimization and intelligent path planning, and relates to a spatiotemporal attenuation factor driven multi-objective ant colony robust optimization method, aiming to balance the robustness and efficiency of the distribution process under an uncertain time-dependent environment. BACKGROUND

[0002] In the process of green and sustainable development of urban road logistics, the time-varying green vehicle routing problem with time windows is one of the key challenges. Due to the complex environmental disturbances in the distribution process, such as traffic congestion, weather changes, etc., it becomes extremely difficult to control carbon emissions while ensuring optimal distribution time. This leads to low optimization efficiency and insufficient robustness, which cannot meet the dual requirements of efficiency and environmental protection in actual logistics distribution.

[0003] Levy flight, as a random walk process, can effectively cover a larger search space and perform well in exploring uncertain environments or solving complex optimization problems. The multi-objective ant colony optimization (MOACO) algorithm can effectively balance the trade-off relationship between different objectives. By combining Levy flight and the improved MOACO algorithm, the simultaneous processing of multi-objective optimization in the distribution process can be realized, and the robustness of the algorithm in complex uncertain scenarios can be improved.

[0004] The application provides a spatiotemporal attenuation factor driven robust optimization method. First, dynamic waiting time window constraints are established to enhance the adaptability of the ant colony to time window constraints in a high disturbance environment. By constructing a double-layer pheromone to enhance the ability to capture uncertain factors, then in the closed-loop optimization of robust solution uniformity evaluation, according to the evaluation results, dynamic Levy flight is performed to adjust the pheromone distribution, and finally through an adaptive feedback mechanism, the relative important parameters are updated to realize a balanced distribution process of robustness and efficiency in a high disturbance environment. SUMMARY

[0005] The application provides a spatiotemporal attenuation factor driven robust optimization method, aiming to balance the robustness and efficiency optimization of the distribution process under an uncertain environment. The robust optimization is realized by the following steps:

[0006] The method steps of the application are as follows:

[0007] Step 1: Set the objective function of the multi-objective time-dependent vehicle routing problem, including the minimization of carbon emissions and total travel time.

[0008] Step 2: Dynamic waiting time window design.

[0009] In the initial stages of delivery, vehicles incur excessive waiting time while searching for the optimal departure time to minimize carbon emissions. This leads to reduced delivery efficiency for subsequent customers due to violations of time window constraints. To balance this conflicting goal, this invention dynamically optimizes the vehicle's waiting time at each customer point, thereby satisfying load limits and time window constraints as much as possible.

[0010]

[0011] T v =V ct1 T n +V ct2 (2)

[0012] Where the waiting time w ijkr It is based on the number K of customers whose deliveries have been completed on the current route. T and the current time T n Determined, serv i The time required to serve customer point i. Here, K a V ct1 V ct2 and T v These are adjustable parameters related to the expected number of customers to be delivered and the delivery time. The constraint has a strong effect on reducing waiting time and is suitable for the initial stage of delivery when time constraints are tight and there are many target customers. Similarly, the constraint has a weak effect on reducing waiting time and is suitable for the later stage of delivery when delivery pressure is lower.

[0013] Step 3: Initial solution generation.

[0014] The multi-objective ant colony optimization algorithm is used to initially solve the problem and obtain an initial solution set. The decision function P of the ants is... ij As shown below:

[0015]

[0016] Among them, P ij This refers to the transition probability from node i to node j. τ B ijt τ represents the average basic pheromone content along the path from i to j. R ijt This represents the average robust pheromone content along the same edge. I∈N represents the reachable client nodes under the constraints. ijkr It is the reciprocal of the waiting time at node k. η is the reciprocal of the distance between clients. α c β c σ and γ c These are all weighting factors used to adjust the influence of each parameter.

[0017] Step 3.1: Constructing the Elite Solution Database.

[0018] The best solution with the front of 1 after each iteration is stored in the constructed elite solution database by fast Pareto sorting, the database stores 30 best solutions, the database is updated by fast Pareto sorting, the overall optimal state is maintained, and the subsequent information update is facilitated.

[0019] Step 4: Ant colony pheromone update.

[0020] The update of the ant colony pheromone can guide the subsequent ants and enhance the identification ability of the optimal path. At the same time, the update of the pheromone considers the capture of uncertain speed characteristics and the change of the objective function value of each generation optimal solution, so that the ant colony can cope with uncertain speed and maintain good solving performance.

[0021] Step 4.1: Ant colony pheromone evaporation design.

[0022] In order to maintain the optimal optimization ability of vehicle distribution under the influence of vehicle speed fluctuation under different external environmental interference intensity, the evaporation intensity of ant colony pheromone is adjusted by capturing the uncertain distribution speed characteristics of the vehicle, and the formula is as follows:

[0023]

[0024] Where, σ v and μ v are the variance and mean value of the uncertain distribution speed, respectively, and ρ v0 is an adjustable control parameter.

[0025] Step 4.2: Double-layer ant colony pheromone driven by space-time decay factor.

[0026] In order to realize intelligent ant colony pheromone update, so as to effectively guide the search direction of ant colony, balance the robustness and optimality of the solution. Therefore, based on the total distribution time and total carbon emission, the space-time decay factor is designed, and the space decay factor driven basic ant colony pheromone update formula is defined as follows:

[0027]

[0028] Where, ρ is the pheromone volatilization factor, x i -x j is the distance between i and j points in the current optimal distribution path, D max is the farthest distance among all customer points, Q is a finer parameter of pheromone, f1 best is the current optimal carbon emission value.

[0029] Then, the robust ant colony pheromone update formula driven by time decay factor is defined as follows:

[0030]

[0031] where w z ijk -w zb ijk is the difference between the departure time of the current delivery vehicle to complete the delivery from customer point i to customer point j in z period and the historical optimal time, U-L is the earliest delivery start time and the latest end time of the distribution center, f2 best is the current optimal carbon emission value, p T is an adjustable parameter, k v and ξ v are the kurtosis and skewness of the uncertain delivery speed of the vehicle, respectively.

[0032] Step 5: Update based on the feedback of solution robustness.

[0033] By evaluating the robustness of the solution, the feedback parameter is transmitted to the pheromone adjustment mechanism based on Levy flight to expand the ant colony exploration space, which helps to jump out of the local optimum and guide the ant colony to find solutions with more robustness.

[0034] Step 5.1: Solution robustness evaluation

[0035] By evaluating the uniformity of the solution distribution in the generation, it is determined whether the result of the ant colony search has reached the best balance between robustness and optimality. The solution robustness evaluation formula is as follows:

[0036]

[0037] The uniformity of the solution is evaluated by the normalized average nearest neighbor distance between the F solutions in the cth generation. Overall, the uniformity of the solution distribution in the iteration process is evaluated by comparing the sliding window to determine the direction of uniformity change. In formula (10), D c represents the uniformity index of the cth generation, 2 represents the square of the difference between the average distance of the ith solution and the F solutions in the cth generation. In formula (11), D f is the average index value of the previous K c generations, and D t represents the average distance of the previous t generations of solutions. By tracking the trend of uniformity change, the improvement degree and convergence situation are quantified.

[0038] Step 5.2: Pheromone adjustment based on Levy flight

[0039] According to the evaluation result of the uniformity of the previous solution distribution, i.e. D fAfterwards, by combining the change direction of the distribution uniformity of the historical solution, the adjustment strength of the Lévy flight is controlled adaptively, and by adjusting the distribution of the global pheromone, the ant colony is guided to jump out of the local optimum. The formula of the pheromone adjustment is as follows:

[0040]

[0041] Wherein, ω is the disturbance strength control parameter of the Lévy flight, and ω0 is the basic disturbance strength control parameter of the Lévy flight, θ f Generally, β is set to 1, and β is the basic step length of the Lévy flight, which is dynamically controlled by c max The maximum iteration number and c current iteration number, so the global pheromone adjustment trend is: in the early iteration, the basic step length is large, which can effectively assist the ant colony search process, expand the search space, and mine the potential robustness and optimal solution. In the later iteration, the basic step length is small, which can assist in improving the convergence ability of the ant colony and focus on searching in the optimal direction. In the whole search process, the total strength of the global pheromone adjustment is controlled by the change direction of the solution space distribution, so that the pheromone adjustment has strong pertinence.

[0042] Step 5.3 Feedback update of relative important factors of double-layer pheromone heuristic weight The double-layer pheromone heuristic weight is updated or reset according to the improvement of the solution after each iteration. When the solution is optimized, the update is performed, otherwise the weight is restored to the default value. The update mechanism is as follows:

[0043]

[0044] γ c+1 = γ0- α c+1 . (16)

[0045] In the application, the solution improvement judgment mode is: the objective function of the best elite solution obtained by the current c generation is quickly Pareto sorted with the solution of the c-1 generation, if the c generation is better than the c-1 generation, it means that the performance of the solution is improved, and α c0 The initial value is set to 2.5, and in addition, the important factors of the basic pheromone and the robust pheromone are updated synchronously, that is, β c+1 = α c+1 , the importance factor of the pheromone and the heuristic information γ c+1 Always keep complementary balance relationship with γ0 value.

[0046] As described above, the advantages of the application are:

[0047] By constructing double-layer ant colony pheromone, a time-space attenuation factor driven ant colony pheromone mechanism is proposed, which can clearly determine the change direction of the objective function of the ant colony search process, capture the distribution characteristics of uncertain speed, and quickly respond to uncertain distribution scenarios. At the same time, by establishing a feedback mechanism based on solution robustness, real-time monitoring and analysis of the distribution of the solution space is realized, which improves the global search ability of the ant colony and effectively balances the robustness and optimality of the distribution process. In the urban road network, the distribution speed of the vehicle is affected by the uncertain fluctuations of the complex environmental factors, and the above method can adapt to and significantly improve the vehicle scheduling ability, effectively reduce the urban carbon emissions, and promote the sustainable development of green city logistics, which has wide application prospects. BRIEF DESCRIPTION OF DRAWINGS

[0048] In order to clearly describe the technical solutions of the embodiments of the present application and the prior art, the drawings required to be cited in the description of the embodiments and the prior art will be briefly described below. Obviously, the following drawings are only part of the embodiments of the present application, and those skilled in the art can further obtain other related drawings based on these drawings without creative labor

[0049] Figure 1 The time-dependent uncertain speed disturbance data provided by the embodiment of the present application.

[0050] Figure 2 、 3 The distribution of the solution and its convergence performance in the high disturbance distribution scenario provided by the embodiment of the present application.

[0051] Figure 4 The optimization result graph of the time-space attenuation factor driven ant colony optimization algorithm provided by the embodiment of the present application. DETAILED DESCRIPTION

[0052] In order to make the purpose, technical scheme and advantages of the embodiments of the present application clearer, the technical scheme in the embodiments of the present application will be described clearly and completely below in combination with the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present application.

[0053] The time-dependent vehicle routing problem model of the present application and its ant colony optimization algorithm application process are described below, which specifically includes the following steps:

[0054] Step 1: Determine the optimization goal of the vehicle routing problem:

[0055] To better define the multi-objective optimization delivery model, it is necessary to clarify the characteristics of the time-varying green vehicle routing problem with time windows. Since vehicle speed changes dynamically with time, the travel time between customer points i and j is affected by the departure time. To characterize this property, we first need to define a directed graph G = (A, C), where C is {(i, j) | i, j ∈ A and i ≠ j}, representing all feasible paths from node i to node j. The node set A = {0, 1, ..., n+1} satisfies the following: node 0 represents the delivery center, serving as the vehicle's starting point; the set N = {1, 2, ..., n} represents the customer nodes to be served; and node n+1 is the final return node, corresponding to the delivery center. All paths must satisfy the following requirement q for each customer node i ∈ A. i and the stipulated delivery time window [e i ,l i ]. Where e i and l i s represents the earliest and latest times when node i can start providing services, respectively. i For service duration, the coefficient d ij This represents the travel distance along path (i,j). Since vehicle speeds on arc (i,j)∈C are time-varying, the travel time and carbon emissions between two client nodes depend on the departure time g of node i∈{0}∪N. i The carbon emission calculation model used in the time-dependent green vehicle route model is a comprehensive modal emission model, an analytical framework for quantifying pollutant emissions from transportation vehicles under different operating conditions. It is commonly used in transportation planning, green logistics optimization, and environmental impact assessment. Its calculation formula is shown below:

[0056] F(t,v,f)=ψ e t+ψ s v 3 t+ψ ω d(μ k +f). (17)

[0057] Where, μ k Let represent the unloaded weight of vehicle k, which is set to 1592 kg in this problem, and ψ e ψ s and ψ ω The coefficients for the engine module, speed module, and load module are respectively. The vehicle is a default logistics delivery van. Under urban conditions, the engine module coefficient is 0.3, the speed module coefficient is 0.06, and the load module coefficient is 0.008 when the vehicle speed is approximately 60 km / h. t is the total delivery time, v is the real-time vehicle speed, f is the cargo being transported, and d is the distance traveled for a cargo weight of f kg. Therefore, the carbon emission model design of this invention is as shown in formula (17).

[0058] The model of the application takes carbon emissions and total travel time of vehicles as optimization objectives:

[0059]

[0060] wherein f1 is the carbon emissions generated by all vehicles constructing paths in the scheduling process, and f2 is the total travel time of vehicles; the number of vehicles k satisfies 0 < k ≤ K, K is the maximum vehicle quantity limit, and K is usually set to 50 when the number of customers is 50, that is, each vehicle that departs must be able to serve at least one customer. C represents all feasible paths from node i to node j, and in the distribution process, the speed of the vehicle changes with time, and the change of the speed is shown in the speed-time graph as follows: the speed change process between adjacent time instants is a piecewise function with a certain slope, and the vehicle distribution time is divided into travel time piecewise functions according to time change instants [M 1 ij …M z ij ,M z+1 ij ] wherein z ∈ Z represents the zth time instant of the piecewise function, and ij represents the starting point and the end point of the distribution in the time interval, and in this problem, the starting point of the distribution time axis of the application is 8:00, and the end point is 20:00, and the speed fluctuates with time every 15 minutes, so the speed changes with time 48 times in total, so the value range of Z is [1, 48], and then r ∈ R, R represents 49 time intervals. In addition, most distribution processes are not completed within one time interval, and the distribution is usually completed across time intervals.v ijk r represents the distribution process of vehicle k from point i to point j in time interval r, and T ijk r (w) represents the vehicle distribution time from point i to point j at speed v ijk r the piecewise linear travel time function of vehicle k from point i to point j at time w, σ zr ijk the slope of the piecewise linear travel time function of vehicle k from point i to point j at time w, θ z ij the intercept of the piecewise linear travel time function of vehicle k from point i to point j at time w, and z+1 ij the piecewise linear function T ijk r (w) of vehicle k from point i to point j in interval [M zr ijk (w) of vehicle k from point i to point j in interval [M z ij (w) of vehicle k from point i to point j in interval [M z+1 ij (w) of vehicle k from point i to point j in interval [M ijk rSlope of (w) ijk w represents the load of vehicle k between point i and point j z ijk w represents the load of vehicle k between point i and point j z ij x represents the departure time of vehicle k from point i to point j z+1 ij x represents the departure time of vehicle k from point i to point j z ijk d represents whether vehicle k is on the optimal path from point i to point j z ij d represents whether vehicle k is on the optimal path from point i to point j z+1 ij d represents whether vehicle k is on the optimal path from point i to point j ij x represents the delivery distance from point i to point j ijk d represents whether vehicle k is on the optimal path from point i to point j

[0061] Each vehicle needs to satisfy the following constraints during delivery, each customer point is visited only once and is on the optimal path:

[0062]

[0063] In equations (21) and (22), the vehicle flow and demand quantity are conserved, where in equation (22), f ij and f jh represent the real-time load of vehicle k from point i to point j and from point j to point h, respectively, and the difference between the two must be equal to the customer demand q j at point j. A and N represent the node set and the delivery node set that satisfies the constraints, respectively.

[0064]

[0065] Equation (23) ensures that the total delivery quantity of each vehicle does not exceed its maximum load Q k , which represents the maximum load of the vehicle. Since the customer demand q i in the example is generally 10, 20, 30, the maximum load of the vehicle Q k is 200.

[0066]

[0067] Equation (24) limits the vehicle between two customer nodes to have multiple time nodes z, requiring only one departure time to be selected x z ijk d represents whether vehicle k is on the optimal path from point i to point j z ij d represents whether vehicle k is on the optimal path from point i to point j z+1 ij d represents whether vehicle k is on the optimal path from point i to point jijk denotes whether vehicle k is on the optimal path i to j, C denotes the set of feasible paths from i to j, and K is the set of vehicles:

[0068]

[0069] Equation (25) requires each vehicle to return to the distribution center, K is the set of distribution vehicles, N is the set of customer points that satisfy the constraints, where x n+1,0，k is a 0-1 variable that denotes whether vehicle k returns from node n+1 to node 0 (distribution center):

[0070]

[0071] Equation (26) is a hard time window constraint that limits the arrival time of a vehicle to not exceed the right boundary of the time window of the next customer node, where b j is the right time window of customer point j, Equation (27) requires that a vehicle must serve each customer point within the specified time window, where a i and b i are the left and right time windows of customer point i, respectively, s ik is the time at which vehicle k starts service at point i, and T ijk r (w) denotes a piecewise linear travel time function with speed v ijk r from point i to point j starting at time w, K and N denote the set of vehicles and the set of distribution nodes that satisfy the constraints, respectively.

[0072]

[0073]

[0074] Equation (28) is a service time constraint for a vehicle, s ik is the time at which vehicle k starts service at point i, K and N denote the set of vehicles and the set of distribution nodes that satisfy the constraints, respectively.

[0075]

[0076] The vehicle generates excessive waiting time at the beginning of the distribution due to searching for the optimal carbon emission departure time point, which leads to a decrease in distribution efficiency due to violation of the time window constraint by subsequent distribution customers. In order to balance the conflict between the two goals, the present application dynamically optimizes the waiting time of the vehicle at each customer point to as far as possible under the conditions of satisfying the load limit and the time window constraint:

[0077]

[0078] T v = Vct1 T n +V ct2 . (30)

[0079] where w ijkr is determined according to the number of customers K T that have been served on the current route and the current time T n , s ik is the time at which vehicle k starts serving point i, and serv i is the time required to serve customer point i. Here, K a , V ct1 , V ct2 and T v are adjustable parameters related to the expected number of delivery customers and time, and in the results of multiple experiments, the average number of delivery customers K T per vehicle is about 5, so K a is set to 8, and when V ct1 , V ct2 are -0.001 and 1.27 respectively, the reduction rate interval of the waiting time window of the vehicle according to the number of delivery customers and the current time is 0.3 to 0.95. Experiments have shown that 0.3 represents a strong waiting time reduction effect, suitable for the early stage of delivery when time constraints are tight and the target customers are many, and 0.95 represents a weak waiting time reduction effect, suitable for the later stage of delivery when the delivery pressure is smaller.

[0080] Step 2 Ant colony algorithm initialization and initial solution generation:

[0081] Step 2.1: Initialization:

[0082] 2.1.1 Read instance data and convert

[0083] Store the x and y coordinates of the customer points in the dataset as vectors, and store the distances between the customer points in the array d according to the Euclidean distance calculation formula, with a unit of kilometers. In addition, there are customer demand, customer left and right time window, and the delivery time window for one day is defined as 8:00 to 20:00.

[0084] 2.1.2 Uncertain time-varying delivery speed generation

[0085] Environmental disturbances can cause uncertainty in vehicle speed, which in turn affects urban carbon emissions and delivery timeliness. The same data characteristics are used for delivery speeds with uncertainty, and the vehicle speed uncertainty set v t is constructed as:

[0086] v z = { | μ z - τσ z ≤ v z ≤ μz + τs z , | | v z - v z-1 | | ≤ Δ max}. (31)

[0087] where μ z = [μ1, μ2, …, μ n ] and μ z ≥ 0 represents the velocity vector within the time period z, σ z = [σ1, σ2, …, σ n ] (σ z ≥ 0) and τ ≥ 1 represent the perturbation reference value of the velocity and the scaling factor of the uncertainty set respectively, Δ m axcorresponds to the physical constraint condition of the velocity, in order to simulate the real urban logistics distribution environment, the value range of the expected velocity μ z is 50-70 km / h, secondly, due to the uncertain interference of multiple factors in the distribution process, the perturbation reference value of the velocity is set to 10 km / h, in addition, the design of the scaling factor of the uncertainty set in the patent is as follows.

[0088]

[0089] wherein the scaling factor τ of the uncertainty set is dynamically adjusted based on the difference between the current time T and the peak time T peak , wherein τ0 is the benchmark fluctuation parameter, the value range of the current time T is 480-1200, T max is the maximum distribution time, which is 1200, in minutes, the peak time T peak is taken from the early and late two peak points, which are 8:00 (480) in the morning and 5:30 (1080) in the evening, in the calculation, the nearest peak period is adaptively selected for calculation, τ0 is taken as -0.15, in order to more truly reflect the speed fluctuation of urban logistics vehicle distribution, the closer the current time T to the peak time T peak , the greater the value of the scaling factor τ, indicating that the speed fluctuation is more intense, on the contrary, when it is away from the peak period, the value is relatively small, the value range of τ is 1.35-1.0, when T = T peak (peak time), (T-T peak ) 2 = 0, the exp term is 1, at this time τ = -0.15 + 1.5 = 1.35 (maximum value), corresponding to the most intense speed fluctuation in the peak period, which matches the setting of "peak period speed perturbation reference value 10 km / h", when T is away from the peak (such as 13:00, i.e. 780 minutes), (T-T peak ) / T maxThe value of exp term tends to 0, and τ ≈τ0+0=-0.15. However, considering the fact that there are still basic fluctuations in off-peak hours, the exponential function eventually makes τ stable at 1.0 (the minimum value), which corresponds to the relatively flat speed fluctuations in off-peak hours, consistent with the reality that off-peak traffic is more stable.

[0090] Step 2.2: Initial solution generation:

[0091] The multi-objective ant colony optimization algorithm is used to solve the problem and obtain the initial solution set. The decision function P ij is as follows:

[0092]

[0093] where P ij is the transition probability from node i to node j. τ B ijt represents the average amount of basic pheromone on the i-j path, τ R ijt represents the average amount of robust pheromone on the same edge. I∈N is the customer node that can be reached under the constraint. w ijkr is the inverse of the waiting time at node k. η is the inverse of the distance between customers. α c , β c , σ and γ c are weight factors that adjust the influence of each parameter, set to 1.25, 2, 0.2, and 1.25 respectively in the actual problem. Among them, α c and γ c are dynamically adjusted in the subsequent algorithm iteration process, and 1.25 is the initial value of the double-layer pheromone. Due to the previous dynamic time window constraint, in order to balance the pressure of the time window constraint under the strong interference of speed, the weight factor value of w ijkr is low. In addition, if there is no waiting time, w ijkr takes a value of 0.1.

[0094] Step 2.3: Elite solution database construction:

[0095] 2.3.1 Construction of elite solution library

[0096] The best solution with a front of 1 obtained after each iteration is stored in the constructed elite solution database, which can store a maximum of 30 best solutions.

[0097] 2.3.2 Update of elite solution library

[0098] When the database is full of 30 best solutions, the new solution and the original solution of the database are updated by fast pareto sorting, the best solution with the front of 1 is retained, and 30 elite solutions are dynamically maintained to facilitate subsequent information updating.

[0099] Step 3 Ant colony pheromone update

[0100] Step 3.1 Ant colony pheromone evaporation mechanism

[0101] In order to maintain the optimal optimization ability of vehicle distribution under the influence of vehicle speed fluctuation under different external environmental interference intensity, the evaporation intensity of ant colony pheromone is adjusted by capturing the uncertain distribution speed characteristics of the vehicle, and the pheromone volatilization factor formula is as follows:

[0102]

[0103] Where, σ v and μ v are the variance and mean value of the uncertain distribution speed, and ρ v0 is an adjustable control parameter, which is set to 0.3, so that the evaporation intensity of ant colony pheromone changes in the range of 0.98-0.92 under the condition that the speed presents different fluctuation amplitudes at different times. When the speed fluctuation is strong, the ant colony search process is more affected by the time window constraint and more unstable, and the search direction is affected, so the evaporation intensity value is relatively low.

[0104] 3.2 Double-layer ant colony pheromone design To realize intelligent ant colony pheromone update, so as to effectively guide the search direction of ant colony and balance the robustness and optimality of the solution. Therefore, the time and space decay factors are designed based on the total distribution time and total carbon emission.

[0105] 3.2.1 Space decay factor driven basic ant colony pheromone The distance between the current distribution vehicle and the data set is calculated by the distance between the two points, and the core idea mainly lies in that the farther the distribution distance of the vehicle between the two points, the thinner the amount of pheromone update, and the distribution distance is too far, which is more likely to be affected by the speed, and the speed fluctuation amplitude is large, which is more likely to violate the time constraint and lead to the failure to meet the requirements of the model. The space decay factor driven basic ant colony pheromone update formula is defined as follows:

[0106]

[0107] Where, ρ is the pheromone volatilization factor, ||x i -x j || is the Euclidean distance between i and j points in the current optimal distribution path, and D maxThe farthest distance between all customer points, Δτ B ijt The increment of pheromone for base information, which is calculated by formula (36), Q is the pheromone update parameter, which is set to 120 in this experiment, f1 best The current optimal carbon emission value, in order to balance the value of double-layer ant colony pheromone update as much as possible, the adjustment parameter ρ c Set to 3.4*10 -4 .

[0108] 3.2.2 Time decay factor driven robust ant colony pheromone

[0109] According to the proportion of the delivery time point of the current delivery vehicle and the total delivery time, participate in the calculation of pheromone update, in the high disturbance environment, the delivery time of the vehicle delivery process is affected by the uncertain speed, the path with long delivery time and less customer points is affected by the larger, and its anti-interference ability is poor. Therefore, the time decay factor driven robust ant colony pheromone update formula is defined as follows:

[0110]

[0111] Wherein, |w z ijk -w zb ijk | is the absolute value of the difference between the departure time of the current delivery vehicle in z period to complete the delivery from customer point i to customer point j and the historical optimal time, the time decay factor T d Defined by formula (38), wherein, the maximum delivery time T max Take 1200, ρ T Adjustment parameter, f2 best The current optimal carbon emission value, the increment of robust pheromone Δτ R ijt The calculation of is defined by formula (39), and in order to balance the magnitude difference between the two objectives, in the experiment, ρ T Parameter is set to 250, and in the calculation of the increment of robust pheromone Δτ R ijt Q is also taken as 120, κ v And ξ v The kurtosis and skewness of the uncertain delivery speed of the vehicle.

[0112] Step 4: Feedback update based on solution robustness.

[0113] Step 4.1 Evaluation of solution robustness

[0114] The uniformity of the solution distribution is evaluated to determine whether the robustness and optimality of the ant colony search have reached the best balance. The robustness of the solution is evaluated according to the following formula (40) and (41):

[0115]

[0116] The uniformity of the solution is evaluated by the normalized average nearest neighbor distance between the F solutions in the cth generation. Overall, the uniformity of the solution distribution is evaluated by comparing the sliding window to determine the direction of change in uniformity. In formula (40), D c represents the uniformity index of the cth generation, 2 represents the square of the difference between the average distance of the ith solution and the F solutions in the cth generation. In formula (41), D f is the average sliding window index of the previous K c generations, and D t represents the average distance of the previous t generations of solutions. By tracking the trend of uniformity, the degree of improvement and convergence can be quantified.

[0117] Step 4.2 Pheromone adjustment based on Levy flight

[0118] According to the evaluation results of the uniformity of the solution distribution, i.e., D f , the adjustment strength of the Levy flight is adapted by combining the direction of change in the uniformity of the distribution of historical solutions, and the distribution of global pheromones is adjusted to guide the ants to escape from local optimality. The formula for pheromone adjustment is as follows:

[0119]

[0120] where ω is the disturbance strength control parameter of the Levy flight disturbance control L(β), which adaptively adjusts the global pheromone in the global range. In the Levy flight disturbance control L(β), the disturbance is deployed by the random walking strategy of the heavy-tailed distribution of the global pheromone by the probability distribution of the Levy flight basic step size β, ω0 is the basic disturbance strength control parameter of the Levy flight, and is set to 1.47. The strength of ω is controlled in the range of 0.5-1.5, and θ f is usually set to 1. β is the basic step size of the Levy flight, which is determined by c maxThe maximum iteration number and the current iteration number c are dynamically controlled, in order to reasonably control the adjustment process of the global pheromone by the Levy flight, if the adjustment intensity is too large, the updating effect of the ant colony pheromone is reduced, and if the adjustment intensity is too small, the search space of the ant colony is limited, therefore, the step constant β0 is set to 0.03, the variation range of the perturbation step is 0.03-0.0003, when the perturbation step is less than 0.005, the current iteration number c is greater than 167 generations, which belongs to the later iteration, the perturbation of the Levy flight is gradually reduced, the ant colony is guided to quickly converge to the local optimum, and therefore better optimization results are obtained.

[0121] In step 4.3, the feedback updating double-layer pheromone heuristic weight of the relative important factor is updated or reset according to the improvement of the solution after each iteration. When the solution is optimized, the updating is performed, otherwise the weight is restored to the default value. The updating mechanism is as follows:

[0122]

[0123] γ c+1 = γ0- α c+1 . (46)

[0124] In the application, the improvement judgment mode of the solution is that the best elite solution target function obtained by the current c generation is quickly Pareto sorted with the solution of the c-1 generation, if the c generation is better than the c-1 generation, it indicates that the performance of the solution is improved, α c0 The initial value is set to 2.5, in addition, the important factors of the basic pheromone and the robust pheromone are updated synchronously, that is, β c+1 = α c+1 , the importance factor of the pheromone and the heuristic information γ c+1 is always kept in a complementary balance relationship with the value of γ0, and the specific value is 5, α0 is set to 1.25, and the relative important factors of the double-layer pheromone are kept consistent, and the purpose is to balance the two directions of optimization and robustness in optimization.

[0125] When the iteration number reaches 200 times, the calculation is stopped and the final distribution result is output, including carbon emission (kg) and total distribution time (min), table 1 and table 2 show the optimization data of two examples, wherein C201 represents an example with more stringent time window constraints, from the normal working time 8:00 to 20:00, the specific setting is that the start time CT is 480, the end time is 1200, the vehicle capacity is 200, the left and right time windows and the service time of the customer are converted into corresponding distribution working time data, and table 3 and table 4 give the route table of the final vehicle path planning.

[0126] Table 1. Vehicle path optimization data of C101

[0127]

[0128]

[0129] Table 2. C201 vehicle routing optimization data

[0130]

[0131]

[0132] Table 3. C101 optimal vehicle routing solution

[0133]

[0134] Table 4. C201 optimal vehicle routing solution

[0135]

[0136]

Claims

1. A robust optimization method for multi-objective ant colonies driven by a spatiotemporal decay factor, characterized in that, Includes the following steps: Step 1: Define the optimization objective for the vehicle routing problem: Because vehicle speed changes dynamically over time, the travel time between customer points i and j is affected by the departure time. To characterize this property, we first need to define a directed graph G = (A, C), where C is {(i, j) | i, j ∈ A and i ≠ j}, representing all feasible paths from node i to node j. The node set A = {0, 1, ..., n+1} satisfies the following: node 0 represents the distribution center, serving as the vehicle's starting point; the set N = {1, 2, ..., n} represents the customer nodes to be served; and node n+1 is the final return node, corresponding to the distribution center. All paths must satisfy the following requirement q for each customer node i ∈ A. i and the stipulated delivery time window [e i ,l i ]; where e i and l i s represents the earliest and latest times when node i can start providing services, respectively. i For service duration, the coefficient d ij Let g represent the travel distance along path (i,j). Since the vehicle speed on arc (i,j)∈C is time-varying, the travel time and carbon emissions between two client nodes depend on the departure time g of node i∈{0}∪N. i The carbon emission calculation formula used in the time-dependent green vehicle path model is shown below: F(t,v,f)=ψ e t+ψ s v 3 t+ψ ω d(μ k +f). (1) Where, μ k Let represent the unloaded weight of vehicle k, which is set to 1592 kg in this problem, and ψ e ψ s and ψ ω The coefficients are for the engine module, speed module, and load module, respectively. The vehicle is a small truck used for logistics delivery by default. Under urban conditions, the engine module coefficient is 0.3, the speed module coefficient is 0.06, and the load module coefficient is 0.008 when the vehicle speed is about 60km / h. t is the total delivery time, v is the real-time speed of the vehicle, f is the cargo being transported, and d is the distance traveled by the cargo with a weight of fkg. Therefore, the carbon emission model design is as shown in formula (2). The model of this invention aims to minimize carbon emissions and total vehicle driving time. Among them, f1 is the carbon emission generated by all vehicles constructing paths during the scheduling process, and f2 is the total vehicle travel time; the number of vehicles k satisfies 0 < k ≤ K, where K is the maximum vehicle number limit, that is, each dispatched vehicle must be able to serve at least one customer; C represents all feasible paths from node i to node j. During the distribution process, the vehicle speed changes with time, and the change in speed is manifested in the speed-time graph as follows: the speed change process between adjacent moments is a piecewise function with a certain slope. The vehicle distribution time is divided into [M 1 ij …M z ij ,M z+1 ij of the travel time piecewise function, where z ∈ Z represents the z-th moment of the piecewise function, Z is the moment of the piecewise function, and ij represents the starting and ending points of the distribution within this time interval. The starting point of the distribution time axis is 8:00, and the ending point is 20:

00. The speed fluctuates with time every 15 minutes. Therefore, the speed changes 48 times in total. So the value range of Z is [1, 48]. Then, in r ∈ R, R represents 49 time periods; Most delivery processes are not completed within a single time frame; they typically span multiple time periods. ijk r This represents the delivery process of vehicle k from point i to point j within time period r, therefore T ijk r (w) indicates that the vehicle travels at a speed of v ijk r A piecewise linear travel time function with departure time w from point i to point j, σ zr ijk In the interval [M] z ij M z+1 ij The piecewise linear function T for delivery from point i to point j within the range. ijk r The intercept of (w), θ zr ijk Indicates the interval [M] z ij M z+1 ij The piecewise linear function T for delivery from point i to point j within the range. ijk r The slope of (w), f ijk w represents the load of vehicle k between point i and point j. z ijk Indicates vehicle k in M z ij To M z +1 ij The departure time from point i to point j between moments, x z ijk Indicates whether vehicle k is in M z ij To M z+1 ij The value of d is 1 if the journey from point i to point j occurs within a given time interval, and 0 otherwise. ij x represents the delivery distance from point i to point j; ijk This indicates whether vehicle k is on the optimal path from i to j; ψ is 1 if it is, and 0 otherwise. e ψ s and ψ ω The coefficients for the engine module, speed module, and load module are set to 0.3, 0.06, and 0.008, respectively. Each vehicle must meet the following constraints during delivery: each customer point is visited only once and always on the optimal route: In formulas (5) and (6), traffic flow and demand for goods are conserved. In formula (6), f ij with f jh This represents the real-time load of vehicle k from point i to point j and from point j to point h. The difference between the two must equal the customer demand q at point j. j A and N represent the set of nodes and the set of delivery nodes that satisfy the constraints, respectively; Formula (7) ensures that the total delivery volume of each vehicle does not exceed its maximum load capacity, Q. k This indicates the vehicle's maximum load capacity, due to customer requirements q. i The maximum load capacity Q of the vehicle is 10, 20, or 30. k It is 200; Formula (8) restricts the vehicle to multiple time points z between two customer nodes, requiring only one departure time x to be selected. z ijk Indicates whether vehicle k is in M z ij To M z+1 ij Between points i and j, x ijk This indicates whether vehicle k is on the optimal path from i to j, where C represents the set of feasible paths from i to j, and K is the set of vehicles. Formula (9) requires that each vehicle must return to the distribution center, where K is the set of delivery vehicles and N is the set of customer points that satisfy the constraints, where x n+1,0，k Let z be a 0-1 variable, representing whether vehicle k returns from node n+1 to node 0 (the distribution center), and z∈Z represent the z-th time step of the piecewise function. Formula (10) is a hard time window constraint, which restricts the arrival time of the vehicle from exceeding the right boundary of the time window of the next customer node, where b j For the right time window of customer point j, formula (11) requires that the vehicle must provide service within the specified time window of each customer point, where a i With b i These represent the left and right time windows for customer point i, respectively. ik T is the time when vehicle k starts service at point i. ijk r (w) indicates that the vehicle travels at a speed of v ijk r The piecewise linear travel time function from point i to point j with departure time w in time period r, where K and N represent the set of vehicles and the set of delivery nodes that satisfy the constraints, respectively. Formula (12) represents the service time constraint for the vehicle, s ik Let K be the time when vehicle k starts serving at point i, and let K and N represent the set of vehicles and the set of delivery nodes satisfying the constraints, respectively. In the initial stages of delivery, vehicles incur excessive waiting time while searching for the optimal departure time for carbon emissions. This leads to reduced delivery efficiency for subsequent customers due to violations of time window constraints. To balance the conflict between these two objectives, the waiting time of vehicles at each customer point is dynamically optimized to satisfy load limits and time window constraints as much as possible. T v =V ct1 T n +V ct2 . (14) Where the waiting time w ijkr It is based on the number K of customers whose deliveries have been completed on the current route. T and the current time T n Determined, s ik For the time when vehicle k starts service at point i, serv i The time required to serve customer point i; here, K a V ct1 V ct2 and T v These are adjustable parameters related to the expected number of customers to be delivered and the time. In the results of multiple experiments, the average number of customers delivered per vehicle, K... T If it is approximately 5, then K a Set to 8, when V ct1 V ct2 When the values ​​are -0.001 and 1.27 respectively, the vehicle waiting time window has a reduction factor ranging from 0.3 to 0.95, depending on the number of vehicles delivering and the current time. Step 2: Ant Colony Algorithm Initialization and Initial Solution Generation Step 2.1: Initialization: 2.1.1 Read and transform instance data The x and y coordinates of customer points in the dataset are stored as vectors. Based on the Euclidean distance calculation formula, the distance between customer points is stored in array d, with the unit being kilometers. In addition, there are customer needs and customer left and right time windows. The daily delivery time window is defined as 8:00 to 20:

00. 2.1.2 Generation of Uncertain Time-Varying Delivery Speed Environmental disturbances can introduce uncertainty in vehicle speed, thereby impacting urban carbon emissions and delivery timeliness. The same data characteristics are applied to delivery speeds with inherent uncertainties, and the vehicle speed uncertainty set v is defined as follows: t Build as: v z ={|μ z -ts z ≤v z ≤μ z +ts z ,||v z -v z-1 ||≤△ max }. (15) Where, μ z =[μ1,μ2,…,μ n And μ z ≥0 represents the velocity vector within time interval z, used as a reference value, σ z =[σ1,σ2,…,σ n ](σ z ≥0) and τ≥1 represent the velocity disturbance reference value and the scaling factor of the uncertainty set, respectively, Δ max Given the physical constraints on velocity, and in order to simulate a real-world urban logistics and distribution environment, the desired velocity μ is... z The value range is 50-70km / h. Secondly, due to the uncertain interference of multiple factors during the delivery process, the speed disturbance reference value is set to 10km / h. In addition, the design of the scaling factor of the uncertainty set of this patent is as follows. The scaling factor τ of the uncertainty set is based on the current time T and the peak time T. peak The difference is dynamically adjusted, where τ0 is the baseline fluctuation parameter, and the current time T ranges from 480 to 1200. max The maximum delivery time is set at 1200, in minutes, with peak time T. peak The values ​​are taken from two peak times, 8:00 AM (480) and 5:30 PM (1080). The calculation will adaptively select the nearest peak time, with τ0 set to -0.

15. The closer the current time T is to the peak time T... peak When T = T, a larger scaling factor τ indicates more severe speed fluctuations, and conversely, a smaller value indicates a lower speed fluctuation when far from the peak. The value of τ ranges from 1.35 to 1.

0. peak (Peak hours), (TT) peak ) 2 =0, exp term is 1, at this time τ = -0.15 + 1.5 = 1.35, corresponding to the most violent speed fluctuation during the peak period, which matches the setting of "peak period speed disturbance reference value 10km / h". When T is far away from the peak, (TT) peak ) / T max The value of τ increases significantly, and the exp term approaches 0. At this time, τ≈τ0+0=-0.

15. However, considering the characteristic of "basic fluctuations still exist during off-peak hours" in actual traffic, the decay characteristics of the exponential function eventually stabilize τ at 1.0, which corresponds to a relatively gentle speed fluctuation during off-peak hours, consistent with the reality that "off-peak traffic conditions are more stable." Step 2.2: Initial solution generation: The multi-objective ant colony optimization algorithm is used to initially solve the problem and obtain an initial solution set. The decision function P of the ants is... ij As shown below: Among them, P ij τ refers to the transition probability from node i to node j; B ijt τ represents the average basic pheromone content along the path from i to j. R ijt Then, it represents the average robust pheromone content on the same edge; I∈N represents the client nodes that can be reached under the constraints; w ijkr η is the reciprocal of the waiting time at node k; η is the reciprocal of the distance between customers; α c β c σ and γ c These are all weighting factors used to adjust the influence of each parameter, and in practical problems, they are set to 1.25, 2, 0.2, and 1.25 respectively; where α c With γ c The values ​​are dynamically adjusted during subsequent algorithm iterations. 1.25 is the initial value for each of the two pheromones. If no waiting time is generated, then w... ijkr The value is 0.1; Step 2.3: Elite Solution Database Construction: 2.3.1 Elite Library Decomposition and Construction The best solutions with a frontier of 1 are obtained by sorting them by Pareto after each iteration and stored in the constructed elite solution database, which can store a maximum of 30 best solutions. 2.3.2 Elite Database Decomposition Update Once the database contains 30 optimal solutions, the new solutions and the original solutions in the database will be updated using fast Pareto sorting, retaining the optimal solutions with a leading edge of 1 and dynamically maintaining 30 elite solutions. Step 3: Ant colony pheromone update Step 3.1 Ant colony pheromone evaporation mechanism By capturing the uncertain delivery speed characteristics of vehicles, the intensity of ant colony pheromone evaporation can be adjusted in a targeted manner. The formula for the pheromone evaporation factor is shown below: Where, σ v and μ v ρ represents the variance and mean of the uncertain delivery speed, respectively. v0 As the adjustable control parameter is set to 0.3, the variation range of ant colony pheromone evaporation intensity is 0.98-0.92 when the speed exhibits different fluctuations at different times. 3.2 Two-layer ant colony pheromone design In order to achieve intelligent ant colony pheromone updates and thus effectively guide the ant colony search direction, balancing the robustness and optimality of the solution, a temporal and spatial decay factor was designed based on the total delivery time and total carbon emissions. 3.2.1 Basic Ant Colony Pheromone Driven by Spatial Decay Factor The basic ant colony pheromone update formula driven by the spatial decay factor is defined as follows: The formula is calculated by proportionally comparing the distance traveled by the delivery vehicle from point i to point j with the maximum distance between two points in the dataset. Where ρ is the pheromone evaporation factor, ||x i -x j || represents the Euclidean distance between points i and j in the current optimal delivery route, D. max Δτ represents the furthest distance between all customer locations. B ijt The increment of the basic pheromone is calculated using formula (20), where Q is the pheromone update parameter, which is set to 120 in this experiment. best To achieve the current optimal carbon emission value, and to balance the pheromone update values ​​of the two-layer ant colony as much as possible, the parameter ρ is adjusted. c Set to 3.4*10 -4 ; 3.2.2 Time-decay factor-driven robust ant colony pheromones The pheromone update calculation is based on the ratio of the current delivery time to the total delivery time. In a highly disturbed environment, the delivery time is affected by uncertain speed, and routes with long delivery times and few customer points are more significantly affected, exhibiting poor anti-interference capabilities. Therefore, the robust ant colony pheromone update formula driven by the time decay factor is defined as follows: In formula (21), |w z ijk -w zb ijk | represents the absolute value of the difference between the departure time of the current delivery vehicle completing the delivery from customer point i to customer point j in time period z and the historical best time, with a time decay factor T. d Defined by formula (22), where the maximum delivery time T max Take 1200, ρ T To adjust the number of deaths, f2 best For the current optimal carbon emission value, the robust pheromone increment Δτ R ijt The calculation is defined by formula (23). Similarly, in order to balance the difference in magnitude between the two targets, ρ T The parameter is set to 250, and similarly, the robust pheromone increment Δτ R ijt In the calculation, Q is also taken as 120, κ v and ξ v These represent the kurtosis and skewness of the vehicle's uncertain delivery speed, respectively.

2. The method according to claim 1, characterized in that... Also includes: Step 4: Feedback update based on solution robustness; Step 4.1 Evaluation of robustness By evaluating the uniformity of solution distribution in the solution distribution space across generations, we can determine whether the ant colony search results have achieved the best balance between robustness and optimality. The robustness evaluation calculation process is shown in formulas (24) and (25): The uniformity of the solution scheme is evaluated by the normalized average nearest neighbor distance among the F solutions in the c-th generation; while the overall uniformity of the solution distribution during the iteration process is evaluated by comparison through a sliding window to determine the direction of uniformity change; where D in formula (24) c This represents the uniformity index of generation c. D represents the square of the difference between the average distance between the i-th solution and the F solutions of the c-th generation. In formula (25), D f For the first K c The index value calculated by the average sliding window of the generation, D t This represents the average distance between the first t solutions; Step 4.2 Pheromone Adjustment Based on Levi's Flight Based on the previous evaluation of the uniformity of solution distribution, i.e., D f Then, by combining the direction of change in the uniformity of historical solution distribution, the adjustment intensity of the adaptive control Levy flight is used to guide the ant colony out of the local optimum by adjusting the global pheromone distribution; the formula for pheromone adjustment is shown below: Wherein, ω serves as the disturbance intensity control parameter of the Lévy flight disturbance control L(β), adaptively adjusting the two-layer pheromone globally. In the Lévy flight disturbance control L(β), the deployment of disturbances is a random walk strategy that uses the probability distribution of the Lévy flight base step size β to distribute the global pheromone with a heavy tail. ω0 is the basic disturbance intensity control parameter of the Lévy flight, with a value of 1.

47. The intensity range of ω is controlled between 0.5 and 1.

5. θ f Typically set to 1, β serves as the base step size for Lévy's flight, determined by c. max The maximum number of iterations and the current number of iterations c are dynamically controlled. The step size constant β0 is set to 0.03, and the perturbation step size varies from 0.03 to 0.0003. When the perturbation step size is less than 0.005, the current number of iterations c has exceeded 167 generations, which is in the late stage of iteration. Step 4.3 Feedback Update on Relatively Important Factors The two-layer pheromone heuristic weights are updated or reset based on the improvement of the solution after each iteration; the update is performed when the solution is optimized, otherwise the weights are restored to their default values; the update mechanism is as follows: c c+1 =γ0-α c+1 . (30) The method for judging the improvement of the solution is as follows: The objective function of the best elite solution obtained in the c-th generation is compared with the solution in the (c-1)-th generation using a fast Pareto sort. If the ranking of the c-th generation is better than that of the (c-1)-th generation, it indicates that the performance of the solution has improved, α. c0 The initial value was set to 2.

5. Furthermore, key factors of both the basal pheromone and robust pheromone were updated synchronously, namely β. c+1 =α c+1 The importance factor of pheromones and heuristic information γ c+1 The complementary balance is always maintained with a value of γ0, which is 5, and α0 is set to 1.25.

Images