A distributed collaborative parallel evolutionary optimization method for green low-carbon logistics transportation

Abstract

This invention proposes a distributed, collaborative, parallel evolutionary optimization method for green and low-carbon logistics transportation. It treats the Vehicle Routing Problem (VRPTW) with time windows as a satisfiability problem and divides the solution process into a vehicle optimization stage prioritizing the number of vehicles and a distance optimization stage under a fixed number of vehicles. This two-stage solution balances the optimization objectives of the number of vehicles used and the travel distance. A fitness penalty mechanism is introduced during the search process, penalizing infeasible solutions under both capacity and time window constraints, thereby allowing infeasible solutions to enter the population and participate in evolution, thus improving population diversity. The method uses the number of computing devices, K, as the parallel scale parameter. Based on the value of K, a corresponding number of parallel subpopulations are constructed. Each subpopulation is configured with a differentiated search mode that balances neighborhood size, waiting time, and time distortion as key control factors, allowing different subpopulations to have different search focuses, thus fully utilizing computing resources.

Description

Technical Field

[0001] This invention belongs to the field of evolutionary algorithm technology and relates to a distributed, cooperative, parallel evolutionary optimization method for green and low-carbon logistics transportation. This algorithm can be widely used in large-scale real-time vehicle routing problems under dynamic and complex environments such as logistics transportation and cargo delivery. Background Technology

[0002] The Vehicle Routing Problem (VRP), first proposed by Dantzig and Ramser (1959) in their seminal paper "The Truck Dispatching Problem," describes a traveling merchant problem in which a fleet of trucks departs from a distribution center, serves a group of customers, and returns, with the goal of minimizing the total cost of travel. It was originally developed to address the truck dispatching problem for gasoline delivery.

[0003] However, the logistics and transportation industry, as a major contributor to carbon emissions, faces unprecedented pressure to reduce emissions. Traditional Vehicle Routing Problems (VRP) focus solely on optimizing distance or cost. In modern green logistics systems, however, achieving the dual minimization of transportation costs and carbon emissions while meeting stringent customer service time windows has become a core challenge for both academia and industry. Green VRPTW (G-VRPTW) emerged to address this need. Its core lies in incorporating environmental costs (such as fuel consumption, carbon emissions, and battery degradation) into the optimization objective function and adding consideration of the physical characteristics of new energy vehicles to the constraints. Solomon (1987), in his pioneering work, not only formalized the definition of VRPTW but also established a set of benchmarks (Solomon Benchmarks) that are still widely used today, laying the foundation for evaluation standards in this field. VRPTW, in addition to meeting the constraints of total vehicle carrying capacity and minimizing the number of vehicles and total distance, introduces an additional time window constraint: each customer has a demand and a time window within which the vehicle must reach the customer. Arriving early incurs an opportunity cost, while arriving late renders the solution invalid. Violating the time window is permissible, but at a high penalty cost. The key computational challenge of VRPTW lies in optimizing vehicle routes under multiple constraints, as failing to meet the time window or total capacity limit renders the solution infeasible, making both VRP and VRPTW typical NP-hard combinatorial optimization problems. As the number of nodes increases, the number of feasible solutions explodes exponentially, making it difficult for exact algorithms to solve large-scale problems in polynomial time.

[0004] Among them, search operators are an important component of heuristic algorithms. Well-known search operators include 2-opt, 3-opt, and swap. 2-opt was proposed by Lin S (1965). This method reduces the total path distance by changing the order of customers in a path (eliminating intersections). If the cost of the path is reduced, the improved path is retained; otherwise, the path reverts to its original state. The 3-opt algorithm is an extension of the 2-opt algorithm. It not only reverses the path between two points but also allows for the reconnection of three non-adjacent vertices in the path, trying seven different connection methods and calculating the path length corresponding to each method. The connection with the shortest path length is selected as the new connection method, and this process continues until no further improvement is possible. This breaks the complex intersection cases that 2-opt cannot solve, helping the algorithm escape local optima and explore a wider solution space. Swap, proposed by Chen and Wu, is a path improvement method. When a vehicle on a path is delayed due to a customer, the "troubled" customer can be swapped with a vehicle on another path that has time leeway. By exchanging demand and different customers, the load factor of each vehicle can be dynamically balanced. These operators mainly explore the neighborhood of the current solution and make small-scale adjustments. However, due to the two constraints of vehicle capacity and customer time window in the VRPTW problem, the neighborhood solutions generated by the local search of these operators are likely to violate these constraints, which is not conducive to getting rid of local optima.

[0005] To address the difficulty of escaping local optima in local search, Shaw (1998) proposed Large Neighborhood Search (LNS). The basic idea is to randomly generate an initial solution; first, randomly or according to a certain strategy, remove a portion of elements from the current solution to introduce variation; then, based on a series of strategies, reconstruct the solution, inserting the removed elements back into the solution in an optimized manner; and finally, determine whether to choose the new solution or maintain the current solution based on the fitness of the new solution. This allows for the acceptance of suboptimal solutions to avoid getting trapped in local optima. The above steps are repeated until predefined conditions are met, such as reaching the maximum number of iterations or finding the ideal solution. While this method effectively explores the solution space and helps the algorithm escape local optima, how to handle the combination and effectiveness of operators during the search process remains a key challenge for the algorithm. Summary of the Invention

[0006] To address the above technical problems, this invention proposes a distributed, collaborative, parallel evolutionary optimization method for green and low-carbon logistics transportation, comprising the following steps: S1. Obtain a vehicle routing model with time windows. The model includes the upper limit of the number of vehicles allowed, the maximum load capacity Q of each vehicle, the coordinates of the warehouse and customer points, the earliest start time window and the latest end time window of each customer point, the service time, and the algorithm parameters. With the upper limit of the number of vehicles allowed and the total driving distance as the optimization objectives, under the conditions of satisfying the vehicle capacity constraints and customer time window constraints, search for the vehicle routing scheme, i.e. the global optimal solution. The parameters include the parallel scale parameter K, the population size N, the maximum number of iterations MaxGen, and the two-stage iteration budget allocation coefficient ŋ. S2, calculate the theoretical minimum number of vehicles Lb based on the total demand D and the single vehicle capacity Q; S3 sets the maximum number of vehicles (MaxNV) during the vehicle optimization phase. S4. Construct K parallel subpopulations based on the parallel scale parameter K, and configure differentiated search mode parameters for each subpopulation, including neighborhood size M, wait time weight waitTW, and time warp weight warpTW. S5, under the MaxNV constraint, the parallel multi-subpopulation cooperative evolution of the vehicle optimization stage is initiated. Each subpopulation independently performs evolutionary operations according to its differentiated search mode, and cross-subpopulation cooperation is achieved through a global knowledge sharing mechanism based on elite solutions. S6, determine if there is a feasible solution that simultaneously satisfies the capacity constraint and the time window constraint; if there is, proceed to S7; otherwise, jump to S10. S7. Select the solution with the smallest number of vehicles NV from all feasible solutions, and denote its number of vehicles as bestNV. S8, determine whether the vehicle optimization phase should be terminated: if bestNV equals Lb, or the number of iterations in the vehicle optimization phase reaches the allocated budget, then terminate the vehicle optimization phase and jump to S10; otherwise, execute S9. S9, update MaxNV to bestNV-1, and return to S5 to continue the search; S10, summarize the best results of each subpopulation in the vehicle optimization stage, select the solution with the smallest number of vehicles NV and the smallest total driving distance TD as the global optimal solution in the first stage, and record its number of vehicles as bestNV; S11, using bestNV as a fixed vehicle number constraint, initiates parallel multi-subpopulation cooperative evolution in the distance optimization phase; Under the constraint of a fixed number of vehicles, each subpopulation continues to evolve to optimize the total driving distance TD; S12, output the global optimal solution of the distance optimization stage as the final result; The time window constraint is handled using a time warp correction mechanism: when a vehicle arrives at a customer point later than its latest end time window, the service start time of that customer point is corrected to its latest end time window to prevent the cumulative impact of lateness on the time feasibility judgment of subsequent customer points.

[0007] Furthermore, the configuration method of the differentiated search mode parameters in step S4 is as follows: For the k-th subpopulation (k=1,2,…,K), its neighborhood size M, waiting time weight waitTW, and time warp weight warpTW are calculated as follows: In the formula, V represents the number of customer points; N represents the population size; The neighborhood reference coefficient, ; This represents the subpopulation difference offset. .

[0008] Furthermore, the evolutionary process of each subpopulation includes: (a) The subpopulation is initialized using a variety of construction strategies, including: the shortest distance strategy, the earliest latest end time window strategy, the maximum demand strategy, the distance urgency strategy, the random order strategy, the polar angle scan order strategy, and the spatial clustering grouping order strategy. (b) Calculate the fitness of an individual, wherein the fitness function is: In the formula, TD represents the total travel distance, EX represents the total capacity violation, and TW represents the total late violation calculated using the time warp correction method. and The penalty coefficient is for adaptive adjustment; (c) Detect the number of iterations in which the optimal solution of the current subpopulation has not been improved since the last update. If the number exceeds the preset stagnation threshold, reinitialize the subpopulation and return to step a; otherwise, continue to step d. (d) Use the tournament selection method to select two parent individuals from the feasible solutions first; (e) Perform crossover on the parent generation to generate offspring using the selective path exchange operator SREX; (f) Perform a local search on the offspring, the local search being based on the differentiated neighborhood set configured in step S4. Constructing a neighborhood set based on rules Iterative improvements are then made by applying a pre-defined combination of local search operators.

[0009] Furthermore, the penalty coefficient and The adaptive adjustment method is as follows: Let the proportion of individuals satisfying the capacity constraint within the statistical period be . The proportion of individuals that satisfy the time window constraint is The target ratio is ,but: , in, To punish the upward adjustment factor, The penalty is to lower the factor.

[0010] Furthermore, the combination of local search operators includes multiple types of local search operators; among them, the ExchangeUY operator is used to exchange customer points between two paths. and Starting with U and Y consecutive client point segments, different scale structural adjustments from single-point swapping to multi-client segment swapping can be achieved by setting different U and Y values; it also includes SwapTails operator and SWAP* operator. The path tail swapping operator (SwapTails) is used to swap the tail sub-paths after a specified client point in two paths, and the re-insertable swapping operator (SWAP*) is used to swap two client points and re-insert them into the optimal position of the other path.

[0011] Furthermore, in the parallel multi-subpopulation co-evolution process of steps S5 and / or S11, cross-subpopulation cooperation is achieved through the following global knowledge-sharing mechanism based on elite solutions: Each subpopulation periodically compares its internal best feasible solution with the global elite solution: if the number of vehicles is less, or the number of vehicles is the same and the total travel distance improves by ≥5% in the early stage of evolution and ≥1% in the later stage, then the global elite solution is updated; Every preset number of generations, each subpopulation receives the current global elite solution. If the solution is better than its own optimal solution, it is added to the subpopulation.

[0012] Furthermore, it also includes a population restart mechanism: when the number of consecutive times the optimal solution within a subpopulation remains unimproved exceeds the stagnation threshold ST, the subpopulation is reinitialized.

[0013] Furthermore, it also includes a population size control mechanism: when the subpopulation size exceeds the threshold PT, it sorts individuals by fitness and removes duplicate individuals and individuals with poor fitness, so that the population size is restored to N.

[0014] Furthermore, the two-stage iterative budget allocation satisfies: the maximum number of iterations in the vehicle optimization stage. Minimum number of iterations from the optimization phase If the vehicle optimization phase terminates prematurely, its remaining iterations are automatically incorporated into the available iteration budget of the distance optimization phase.

[0015] Compared with the prior art, the beneficial effects of the present invention are: 1. By employing a two-stage hierarchical optimization mechanism, VRPTW is modeled as a satisfiability problem, divided into a vehicle number optimization stage and a distance optimization stage, achieving dual-objective collaborative optimization of both vehicle number and travel distance. This method, while minimizing the number of vehicles, further optimizes the travel distance, overcoming the problem in traditional methods where both are difficult to simultaneously achieve.

[0016] 2. By introducing a fitness evaluation mechanism with penalties, infeasible solutions are allowed to participate in evolution, which effectively improves population diversity, avoids the algorithm getting stuck in local optima in the early stages, and enhances global search capabilities.

[0017] 3. A multi-subpopulation parallel evolution strategy based on the number of computing devices is adopted. Each subpopulation is equipped with differentiated neighborhood construction and time window processing strategies to form complementary search behavior, make full use of computing resources, and significantly improve the solution efficiency of large-scale problems.

[0018] 4. By adjusting the adaptive penalty coefficient, cross-population migration based on elite solutions, and global knowledge sharing, we can accelerate convergence while maintaining population diversity, avoid premature convergence, and improve the robustness and optimization quality of the algorithm.

[0019] 5. Design a variety of local search operators and a stagnation restart mechanism to enhance the local optimization ability and the ability to escape local optima, ensuring that the algorithm can continue to improve under complex constraints. Attached Figure Description

[0020] Figure 1 This is a flowchart of the main overall method of the present invention; Figure 2 This is a flowchart illustrating the internal evolution of the subpopulation in this invention. Figure 3 This is a schematic diagram illustrating a problem of the present invention; Figure 4 for Figure 3 The following image; Figure 5 for Figure 4 The following image; Figure 6 This is a comparison chart of the experimental results. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of the embodiments of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0022] Example A distributed, collaborative, parallel evolutionary optimization method for green and low-carbon logistics transportation is applied to a delivery service area consisting of one warehouse, V customers, and multiple vehicles. Taking the warehouse as the starting and ending point, and the number of vehicles (NV) and total travel distance (TD) as the main optimization objectives, the method solves for a delivery plan that simultaneously satisfies the time window requirements and capacity demands of all customers, while satisfying vehicle capacity constraints and customer time window constraints. The method includes the following steps: S1, First, obtain the path planning model with time windows. Path planning models with time windows Figures 3-5 Taking C121 as an example, the figure shows the basic data of the model. The maximum number of vehicles allowed is 50, the maximum load capacity Q of each vehicle is 200, the warehouse and V customer numbers are in the first column of the table, where customer number 0 represents the warehouse, and the rest are customer numbers; the x and y coordinates of the i-th customer are in the second and third columns respectively, the demand of the i-th customer is in the fourth column, the earliest start time window and the latest end time for the i-th customer are in the fifth and sixth columns, the service time of the i-th customer is in the seventh column, and all vehicles start and end at the warehouse.

[0023] After completing the instance data reading of the path planning model with time windows, the parallel scale parameter K, subpopulation size N, stopping condition MaxGen, and allocation coefficient ŋ for two-stage optimization are further determined. K represents the number of available computing devices. The method determines the number of parallel subpopulations based on the value of K. N specifies the population size of each subpopulation. MaxGen limits the maximum number of iterations. ŋ allocates the iteration budget for the vehicle optimization stage and the distance optimization stage. The maximum number of iterations in the vehicle optimization stage is GenV, and the minimum number of iterations in the distance optimization stage is GenD. The value of K is determined by the computing devices. N is set to 25, MaxGen to 100,000, and ŋ to 20%. In the vehicle optimization stage, GenV can use up to 80% of MaxGen. The first attempt uses 10% of GenV, and every two subsequent attempts increase by 5000 iterations. GenD uses at least 20% of MaxGen.

[0024] S2, calculate the lower bound Lb of the theoretical number of vehicles based on capacity constraints. First, sum up the demand of all customers to obtain the total demand D, and denote the maximum load capacity of a single vehicle as Q. Since in any feasible delivery plan, the load of each vehicle cannot exceed Q, and all customer demands must be fully delivered, when only considering the necessary condition of capacity constraints, the minimum number of vehicles required to cover the total demand is Lb, which is equal to the ratio of the total demand rounded up to the single vehicle capacity. The formula for calculating Lb is as follows (1): (1) In the formula, Lb represents the theoretical lower limit of the number of vehicles determined solely by capacity constraints, serving as the basis for determining whether to terminate the subsequent vehicle optimization phase. Next, we proceed to the first phase: the vehicle optimization phase.

[0025] S3 determines the maximum number of vehicles (MaxNV) available for vehicle optimization. In this stage, the maximum number of available vehicles given by the specific problem instance in S1 is initially set as MaxNV. Subsequently, during iterations, MaxNV is dynamically adjusted based on the update results obtained in S9, serving as a constraint on the upper limit of the number of vehicles used in subsequent vehicle optimization stages.

[0026] S4, based on the parallel scale parameter K, subpopulation size N, and number of customer points V obtained in S1, construct K parallel subpopulations. Each subpopulation Configure a set of different search parameters, and configure different neighborhood sets according to the corresponding search parameters. Develop rules to guide the creation of dedicated neighborhood sets. The construction of the neighborhood set used by each subpopulation when executing the local search operator. The difference makes the neighborhood set of the subpopulation with smaller subpopulation number k different. More inclined to exclude customers who are more likely to cause lateness or breach of contract. This makes the local search more biased towards paths that satisfy the latest end time window and improves the efficiency of advancing the time window feasibility. The neighborhood set of the subpopulation with the larger subpopulation number k is also considered. More inclined to retain customers with smaller waiting times This makes local searches more likely to form paths that better match the arrival time with the customer's available service start time, thereby reducing waiting accumulation and improving the compactness of the schedule, ultimately forming complementary search strategies in the process of parallel evolution.

[0027] The search parameters include: neighborhood size M, wait time weight waitTW, and time warp weight warpTW; Neighborhood set It is the range of candidate points for the local search operator operation, the neighborhood set. The construction rules include: each customer point neighborhood set Selected from customer points Other customer points The customer points corresponding to the top M nearest neighbors after sorting in ascending order. The nearest neighbor value is determined by the distance between two points and the waiting time waitT. With Time WarpT They are jointly constituted and weighted by waitTW and warpTW respectively, as shown in Equation 2. The smaller the proximity value, the more matched the two points are in terms of time arrangement and the less likely it is to introduce waiting and lateness violations.

[0028] (2) in, Indicates customer point With customer points The distance between them.

[0029] The neighborhood size M is mainly used to constrain the neighborhood set of local search operators. The number of candidate points is determined to control the computational scale while ensuring search effectiveness. The neighborhood size M is calculated using Equation 3, where... To serve as a baseline order of magnitude for determining the neighborhood size, M is made to adapt adaptively to the number of customer points V and the subpopulation size N. Both methods are used to generate stable difference amplitudes for different subpopulations. Both aim to avoid situations where the neighborhood of some subpopulations is too large, resulting in excessive computational overhead, or the neighborhood is too small, resulting in insufficient search coverage, under different scale problems. The calculation formula is shown in Equation 4. The definition is shown in Equation 5.

[0030] (3) (4) (5) The wait time weight (waitTW) is used to control the neighborhood set. The strategy is constructed to assess its sensitivity to early arrival waiting, where the waiting time is waitT. Indicates that the vehicle departed from the customer's location. Head straight to the customer service point The minimum possible waiting time, reflecting the arrival time of the vehicle. The strategy treats the waiting caused by early arrival as a constraint violation of underutilization of the time window, and weights it using waitTW. The larger the waitTW value, the stronger the amplification effect of the nearest neighbor value on the waiting item, leading to... More inclined to retain customers with smaller waiting times This allows local searches to more readily generate path structures that better match the arrival time with the customer's available service start time, thereby reducing wait accumulation and improving the compactness of scheduling. Waiting time (waitT) The calculation method is shown in Equation 6, and the calculation method of the waiting time weight waitTW is shown in Equation 7.

[0031] (6) In the formula, Indicates customer point The earliest start time window, i.e. the earliest time when service is allowed to begin; Indicates customer point The latest end time window, i.e. the latest time when service is allowed to begin; Indicates customer point With customer points The travel time between them is equivalent to the distance between the two points. ; Indicates at the customer point Required service time.

[0032] (7) The time warp weight (warpTW) is used to control the sensitivity of neighborhood construction to late violations, where the time warp weight (warpTW) is... Indicates that the vehicle departed from the customer's location. Head straight to the customer service point The minimum amount of delay violation that may occur when the vehicle arrives reflects the minimum amount of delay violation that may occur when the vehicle arrives. The degree to which it is later than its latest closing time window. Time warping uses the idea of ​​arrival time correction, that is, if the vehicle arrives at the customer's point... The customer is already late, so the order is placed... The service start time is no longer calculated based on the current arrival time, but rather on... The latest end time window is used as the service start time. This way, when a lateness occurs at a customer point, the lateness is not propagated to subsequent customers in the original way. Instead, a correction mechanism reduces the collateral impact of the lateness on the feasibility assessment of subsequent customers, thus avoiding misjudgments of subsequent sub-paths due to accumulated lateness and protecting potentially well-structured subsequent sub-paths. The larger the warpTW value, the stronger the amplification effect of the nearest neighbor value on the time warp term, leading to... More inclined to exclude customers who are more likely to cause lateness violations. This makes local searches more biased towards satisfying the latest end time window and improves the efficiency of advancing time window feasibility. (Time WarpT) The calculation method is shown in Equation 8, and the calculation method of the time warp weight warpTW is shown in Equation 9.

[0033] (8) In the formula, Indicates customer point The earliest start time window, i.e. the earliest time when service is allowed to begin; Indicates customer point The latest end time window, i.e. the latest time when service is allowed to begin; Indicates customer point With customer points The travel time between them is equivalent to the distance between the two points. ; Indicates at the customer point Required service time.

[0034] (9) In S5, under the constraint of MaxNV, the maximum number of vehicles that can be used for vehicle optimization determined in S3, the parallel co-evolution of K subpopulations in the vehicle optimization phase is initiated. Each subpopulation is based on its dedicated neighborhood set determined in S4. Neighborhood set constructed by rules Each subpopulation independently performs evolutionary operations based on its differentiated search pattern, and achieves cross-subpopulation collaboration through global knowledge sharing based on elite solutions during the iteration process, continuously updating the current best result of the vehicle optimization stage; the number of iterations does not exceed the maximum number of iterations GenV allocated in the vehicle optimization stage in S1, and the evolutionary process within each subpopulation is described in detail in S501~S513.

[0035] S6: Determine if there exists a feasible solution that simultaneously satisfies both the capacity constraint and the time window constraint; if so, proceed to S7; otherwise, jump to S10.

[0036] After completing one parallel evolution in the vehicle optimization phase, the optimal solutions of each subpopulation are summarized and their feasibility is determined. If a feasible solution that simultaneously satisfies both the capacity constraint and the time window constraint is found, it indicates that there is still room for further reduction in the number of vehicles under the current maximum vehicle number constraint (MaxNV), and the process proceeds to S7. If no feasible solution that simultaneously satisfies both the capacity constraint and the time window constraint is found, it indicates that further optimization of the number of vehicles is not possible under the current MaxNV constraint, and the vehicle number reduction process terminates here. At this point, the method takes the minimum number of vehicles corresponding to the most recently successfully found feasible solution as the best result of the vehicle optimization phase and directly jumps to S10 for execution.

[0037] S7: Select the solution with the minimum number of vehicles (NV) from all feasible solutions and denote its vehicle count as bestNV. If and only if S6 determines this, the method will aggregate all feasible solutions and select the solution with the minimum NV based on the vehicle count (NV) as the priority criterion, denoteing it as bestNV. This bestNV will be used to update the subsequent vehicle count upper limit MaxNV and will serve as the termination criterion for S8.

[0038] S8, determine whether the vehicle optimization phase should be terminated: if bestNV equals Lb, or the number of iterations in the vehicle optimization phase reaches the allocated budget, then terminate the vehicle optimization phase and jump to S10; otherwise, execute S9. Lb is the lower bound of the theoretical number of vehicles obtained in S2. If bestNV equals Lb, it indicates that the number of vehicles has reached the theoretical lower bound under the capacity constraint, and there is no longer room for further reduction in the number of vehicles in the vehicle optimization stage. The process then proceeds to S10 for the next step. If bestNV is greater than Lb, it indicates that there is still potential for further optimization of the number of vehicles. If the current iteration count reaches the maximum iteration count GenV allocated to the vehicle optimization stage in S1, it indicates that the iteration budget for the vehicle optimization stage has been exhausted, and vehicle optimization will no longer continue, proceeding to S10. If the current iteration count has not reached the maximum iteration count GenV allocated to the vehicle optimization stage in S1, it indicates that the iteration budget for the vehicle optimization stage has not yet been exhausted. The two decision conditions are connected by an OR logic. The vehicle optimization stage terminates and proceeds to S10 if either condition is true. Only when both conditions are false does the process proceed to S9 to continue shrinking the upper limit of the number of vehicles.

[0039] S9 updates the maximum number of vehicles usable during the vehicle optimization phase, MaxNV. If S8 determines otherwise, using the current optimal number of vehicles (bestNV) obtained in S7 as a baseline, MaxNV is updated to bestNV minus one, i.e.: This allows subsequent vehicle optimization phases to continue searching for feasible solutions under a stricter upper limit on the number of vehicles. The aim is to verify whether there are feasible path schemes using fewer vehicles, thereby achieving a progressively decreasing search for the number of vehicles and providing a new upper limit on the number of vehicles for the next round of parallel evolution.

[0040] S10, summarize the best results of each subpopulation in the vehicle optimization stage, select the solution with the smallest number of vehicles NV and the smallest total driving distance TD as the global optimal solution in the first stage, and record its number of vehicles as bestNV; When S6 determines that no feasible solution has been found or S8 determines that the vehicle optimization phase has terminated, the optimal results of each parallel subpopulation in the vehicle optimization phase are summarized and compared. The global optimal solution is selected according to the following two-level ranking criteria: the global optimal solution of the first phase is selected based on minimizing the number of vehicles used (NV); if NV is the same, the global optimal solution of the first phase is selected based on minimizing the total travel distance (TD). The number of vehicles used in this solution is recorded as bestNV. Subsequently, this global optimal solution is passed into the distance optimization phase as an elite solution, and bestNV is used as a fixed vehicle number constraint for the distance optimization phase to initialize and guide the distance optimization process of the second phase.

[0041] In S11, using bestNV as a constraint on the fixed number of vehicles, the parallel multi-subpopulation co-evolutionary process for distance optimization is initiated. In this stage, under the constraint of bestNV determined in S10, K subpopulations are run simultaneously, each subpopulation based on its dedicated neighborhood set determined in S4. Neighborhood set constructed by rules Each subpopulation independently performs evolutionary operations based on its differentiated search pattern, and achieves cross-subpopulation collaboration through global knowledge sharing based on elite solutions during the iteration process, continuously updating the current best result of the vehicle optimization stage; the number of iterations does not exceed the maximum number of iterations GenV allocated in the vehicle optimization stage in S1, and the evolutionary process within each subpopulation is described in detail in S501~S513.

[0042] Compared to S5, the maximum number of vehicles MaxNV in this optimization phase has been updated to bestNV, and the optimization objective for this phase has changed from the number of vehicles to the total driving distance.

[0043] S12, output the global optimal solution for the second stage distance optimization and end.

[0044] When the number of iterations in the distance optimization stage reaches the stopping condition MaxGen set by S1, the optimal results obtained by each parallel subpopulation under the constraint of a fixed number of vehicles are summarized and compared. The solution with the smallest total travel distance is selected as the global optimal solution in the second stage, and it is used as the final output result of the method of this invention.

[0045] The number of vehicles used is the theoretically achievable minimum, at which the total operating cost is lowest, satisfying the time window and vehicle capacity constraints for all customers.

[0046] Within each subpopulation, the evolutionary process of the vehicle optimization phase and the distance optimization phase is consistent. Therefore, the following discussion focuses on the subpopulations. Taking this as an example, the evolutionary process of the vehicle optimization stage and the distance optimization stage, S501~S513, is introduced: S501, Initialize the population .

[0047] Input the basic data of the path planning model with time window, the maximum number of vehicles allowed in this optimization (MaxNV), and the iteration budget. The basic data includes the number of customer points (V), the maximum number of iterations in the vehicle optimization phase (GenV), and the minimum number of iterations in the distance optimization phase (GenD). If the vehicle optimization phase terminates prematurely before exhausting GenV, the remaining iterations will be incorporated into GenD for distance optimization.

[0048] subpopulation Seven initialization strategies were employed to obtain the initial solution set, including four heuristic greedy construction strategies and three random or geometric order construction strategies. Greedy construction was guided by the shortest distance, earliest latest end time window, maximum demand, and distance urgency, respectively. Random or geometric order construction first generated the customer visit sequence and then constructed paths under capacity and time window constraints. To ensure structural differences within the initial population, subpopulations... The initialization strategy index used by the nth individual is the remainder of n modulo 7. Each initialization strategy index corresponds to an initialization strategy, and an initial solution is generated according to it. The above process is repeated until the subpopulation size reaches N.

[0049] When constructing individual p, the shortest distance strategy, excluding the warehouse... In addition, all customer points All were added to the unassigned customer set UV, with vehicle route number as follows: , During initialization, each vehicle path of individual p... Initially, it only contains warehouses. Then, without violating time window and capacity constraints, the customer point closest to the end node of the path is selected from the UVs. Add to The tail end and remove it from the UV; when After the tail node is updated, the selection process is repeated based on the new tail customer point until... If it is impossible to insert more customer points without violating the constraints, then confirm. The final route is calculated and the number of vehicle paths NV for individual p is incremented. The next vehicle path is then constructed until UV is empty or the number of vehicle paths reaches MaxNV. After construction according to the above process, if UV is still not empty, the remaining customer points in UV are sequentially added to the end of the vehicle path that causes the smallest distance increment, where the distance increment refers to the increase in total travel distance caused by customer insertion. Let the customer before the insertion position be denoted as p. The latter customer is referred to as Then the customer to be inserted The formula for calculating the distance increment is shown in Equation 10: (10) Here, the customer point is inserted at the end of the vehicle path. Can be considered empty; symbol This represents a 2-norm, which in this embodiment represents the Euclidean distance between vectors.

[0050] The "latest end time window earliest" strategy follows the basic process of the "shortest distance" strategy when constructing individual p. The only difference lies in the selection criterion for determining the next service customer point from the unassigned customer set (UV) each time. The "latest end time window earliest" strategy selects the customer point with the earliest latest end time window from the UV each time. Add to Remove the tail end from the UV, and the rest of the operations are the same as the shortest distance strategy.

[0051] The demand-maximization strategy follows the basic process of the shortest distance strategy when constructing individual p. The only difference lies in the selection criterion for determining the next service customer point from the unallocated customer set UV each time. The demand-maximization strategy selects the customer point with the highest demand from the UV each time. Add to Remove the tail end from the UV, and the rest of the operations are the same as the shortest distance strategy.

[0052] The distance urgency strategy, when constructing an individual p, follows the basic process of the shortest distance strategy. The only difference lies in the selection criterion for determining the next service customer point from the unassigned customer set (UV) each time. The distance urgency strategy selects the customer point with the smallest urgency value from the UV each time. Add to Remove the tail end from the UV; all other operations are the same as the shortest distance strategy. Urgency value This is used to simultaneously reflect spatial proximity and time window urgency, ensuring that the construction process maintains a compact path space while prioritizing customer points with tighter latest completion windows to avoid delays. Customer Point With customer points The urgency value between The calculation method is as shown in Equation 11: (11) In the formula, Indicates customer point With customer points The distance between them; Indicates customer point The latest end time window, i.e. the latest time when service is allowed to begin; Indicates vehicle arrival The arrival time.

[0053] When constructing an individual p, the random order strategy first processes the individual p excluding the warehouse. All customer points outside Generate a random access sequence and use this sequence as the traversal order for UVs; then attempt to access customer points sequentially according to this traversal order. Add to the current vehicle path without violating time window constraints and capacity constraints. The tail end is removed from the UV, if If it is no longer possible to accommodate more customer points without violating the constraints, then it is determined that... The final route is determined and the number of vehicle paths NV for individual p is incremented. The next vehicle path is then constructed until UV is empty or the number of vehicle paths reaches MaxNV. After the construction is completed according to the above process, if UV is still not empty, the remaining customer points in UV are added sequentially to the end of the vehicle path that causes the smallest distance increment.

[0054] The polar angle scanning sequence strategy, when constructing individual p, first uses the warehouse... To calculate the polar angle of each customer point relative to the positive x-axis, the positive x-axis direction is defined as 0°. The points are sorted in ascending order of polar angle. A starting position is randomly selected on the sorted sequence to form a cyclic access sequence, and this sequence is used as the traversal order of UVs. All other operations are the same as the random order strategy.

[0055] The spatial clustering grouping strategy first performs spatial clustering based on the coordinates of customer points when constructing individual p, with the number of clusters set to MaxNV. Then, MaxNV customer points are randomly selected from the customer set as initial cluster centers, and the remaining customer points are then clustered. Calculate its Euclidean distance to each cluster center, and then... Customers are assigned to the clusters corresponding to the nearest centers, thus forming various customer clusters. After clustering, the customer points within each cluster are first randomly shuffled to form the intra-cluster access order. Then, according to the selection order of the cluster centers, each cluster is taken out in sequence and concatenated to form the overall customer access sequence. This sequence is used as the traversal order of UVs. All other operations are the same as the random order strategy.

[0056] S503, calculate the individual fitness value.

[0057] The fitness function is as shown in Equation 14: (14) In the formula, TD represents the total travel distance, EX represents the total capacity violation, and TW represents the total late violation calculated using the time warp correction method. and The penalty coefficient is for adaptive adjustment; This invention models the VRPTW problem as a satisfiability problem, treating capacity constraints and time window constraints as mandatory conditions, with the upper limit of the number of vehicles (MaxNV) as an externally given constraint. When a solution can be found that simultaneously satisfies both capacity and time window constraints under the MaxNV constraint, the problem is considered satisfiable under the upper limit of the number of vehicles; when no such solution can be found, the problem is considered unsatisfiable under the upper limit of the number of vehicles. Based on this modeling approach, the vehicle optimization phase approximates the minimum number of vehicles by continuously tightening MaxNV and repeatedly determining satisfiability, while the distance optimization phase continues to optimize the travel distance under a fixed number of vehicles.

[0058] To enable infeasible solutions to participate in evolution and improve population diversity, this invention introduces a penalty mechanism in fitness evaluation. First, the relevant parameters of individual p are defined. Individual p consists of several vehicle paths, the number of vehicles used is denoted as NV, and the total travel distance is denoted as TD. To determine the degree of capacity constraint violation, the cargo capacity of the nth vehicle path is denoted as the sum of the customer demand served by that path, and the maximum load capacity of a single vehicle is denoted as Q. The capacity violation is then defined as the sum of the overload amounts of all paths, denoted as EX, and its specific calculation formula is shown in Equation 12. To determine the degree of time window constraint violation, the customer point... The earliest start time window is The latest end time window is Service hours are The travel time between the two points is The demand is ,arrive The time is For any path, the arrival time and service start time of each customer point on the path are recursively calculated according to the access order, and the late violation amount, denoted as TW, is calculated using a time warp correction method. The specific calculation formula is shown in Equation 13. An individual p is considered feasible if and only if EX is 0 and TW is 0 (S502).

[0059] (12) (13) The fitness of individual p is composed of the total travel distance TD and the penalties EX and TW for violating the two types of constraints, forming a fitness with penalty values ​​used for ranking and selecting infeasible solutions. The capacity penalty coefficient is denoted as... The time window penalty coefficient is denoted as This penalty term ensures that feasible solutions are judged solely by distance when the violation is zero, while infeasible solutions are penalized based on the degree of violation. This avoids insufficient exploration of the solution space in the early stages of the search due to the strict rejection of infeasible solutions.

[0060] This step employs a penalty management mechanism. and Adaptive adjustment is performed based on the population's satisfaction of capacity and time window constraints during the search process. When a constraint becomes more difficult to satisfy, the corresponding penalty coefficient is increased; when a constraint becomes easier to satisfy, the corresponding penalty coefficient is decreased. Since a larger penalty coefficient results in a stronger penalty for individuals violating the corresponding constraint in terms of fitness, making them less likely to be retained in the selection process, increasing the penalty when a constraint is more difficult to satisfy enhances the suppression of such violations and promotes the population's convergence towards the feasible region. Conversely, decreasing the penalty when a constraint is easier to satisfy avoids excessive penalty leading to premature contraction of the search scope. Let's assume that, within the current statistical period, the proportion of individuals satisfying the capacity constraint out of the total number of individuals being evaluated is... The proportion of individuals that meet the time window constraint out of the total number of individuals evaluated is And set the target ratio as The penalty adjustment factor is The penalty reduction factor is Then, the two types of penalties are updated independently according to Equation 15, thereby allowing a certain proportion of infeasible solutions to participate in evolution to maintain diversity, while continuously controlling the degree of violation and guiding the population to gradually enter and stably remain near the feasible region. The statistical period is set to 500 iterations, and the target proportion... The penalty upregulation factor is 0.9. The penalty downregulation factor is 1.25. It is 0.85.

[0061] (15) In the formula, Corresponding to or , Corresponding to or This enables independent adaptive adjustment of capacity penalty and time window penalty.

[0062] S504, at this stage, in the population In each iteration, a stall state detection is performed; Let pc be the number of consecutive unimproved iterations of the current population's optimal solution since the last update, and ST be the preset stagnation threshold, with the initial value of ST set to 2,0000; where the pc value is updated in S508; If the pc exceeds ST, it indicates that the population has entered a long-term stagnation and will proceed to S505 for restart; otherwise, proceed to S506.

[0063] S505, if the determination of S504 is yes, then for the population Perform a restart, i.e., reinitialize the population, enter S501, and after the restart is complete, reset the stagnation count (PC) to zero and use the updated... The subpopulation continues to evolve as input for subsequent iterations.

[0064] S506, tournament method to select two parent individuals.

[0065] In this stage, parent generation is preferentially extracted from the feasible solution set. Each time, two feasible individuals are randomly selected for fitness comparison, and the one with the better fitness is selected as the parent. When the number of feasible solutions is insufficient to complete the extraction, additional individuals are randomly selected from the infeasible solution set, and the selection is performed in the same manner. If only one feasible solution exists in the current population, that feasible solution is directly designated as the first parent. The second parent is selected by randomly selecting two individuals from the infeasible solution set for comparison, and the one with the better fitness is chosen, ensuring that parent selection prioritizes feasible structures. This stage obtains two parents sequentially according to the above process, avoiding the two parents being the same solution during the selection process, thus forming the parent input set for subsequent crossovers.

[0066] S507 uses the selective path exchange operator SREX to perform crossover on two parent generations to generate child solutions, and evaluates the child solutions according to the fitness calculation method defined in S503.

[0067] S508 performs a local search on the child solution.

[0068] This stage takes the offspring solutions generated in S507 as input and configures a differentiated neighborhood set based on the specific search parameters of the subpopulation numbers in S4. Establish rules, and then construct a dedicated neighborhood set. Restricting any customer point It can only be used with sets located in its neighborhood. Other customer points within Various local search operations are performed to ensure effective search coverage while controlling the computational scale. To achieve differentiated enhancement and complementary search, the method of this invention pre-defines ten combinations of local search operators and assigns specific operator combinations to each subpopulation based on the modulo operation result of the subpopulation number k. These local operator combinations include node operators and path operators. After determining the operator combinations, the local search is executed in the order of node operators first, then path operators: First, node operators are repeatedly called to perform small-scale structural adjustments using the internal customer points as the operation unit. As long as the fitness of the individual after the move is better than the current solution, the move is accepted and the solution is updated, until node operators can no longer produce improvement. Then, path operators are called to implement stronger structural adjustments using the internal vehicle paths as the operation unit, again adopting an iterative strategy of accepting improvements, until neither type of operator can further improve the current solution. Finally, the sub-solution enhanced by the local search is output and compared with... The internal optimal solution is compared to determine whether it has been updated. If the optimal solution has been updated, the population stagnation counter pc is reset to 0; otherwise, pc is incremented.

[0069] The ten local search operator combinations are: {Exchange10, Exchange11, SwapTails, SWAP*}; {Exchange10, Exchange20, SwapTails}; {Exchange10, Exchange11, Exchange20 , SwapTails, SWAP*}; {Exchange10, Exchange21, Exchange11, SwapTails, SWAP*}; {Exchange10, Exchange30, Exchange20, SwapTails, SWAP*}; { Exchange10, Exchange22, Exchange11, SwapTails, SWAP*}; {Exchange10, Exchange31, Exchange21, SwapTails, SWAP*}; {Exchange10, Exchange 32, Exchange33, SwapTails, SWAP*}; {Exchange10, SwapTails, SWAP*}; {Exchange10, Exchange11, Exchange20, Exchange21, SwapTails, SWAP*}.

[0070] Among them, the node operator is based on the ExchangeUY series of exchange operators and adopts a continuous segment exchange mechanism: the client point is first selected among the individuals. Then in Selected customer points And determine the positions of the two on their respective paths, from Extract from the path A continuous U customer point segment starting from [the origin] ,from Extract from the path A continuous segment of Y customer points starting from [the point name]. The two segments are then swapped to form a new path connection. Different combinations of U and Y can cover local adjustments at different scales: when both U and Y are 1, the two points are swapped; when Y is 0, it degenerates into a segment. Remove and Insert Next, the U consecutive client point segments starting from vi will be removed from the original path and inserted into the new path. Subsequently, when U or Y takes a larger value, it is possible to move or swap multiple consecutive segments of customers, thereby changing the local path structure with different perturbation amplitudes and improving the ability to escape local optima.

[0071] Furthermore, the method of this invention introduces two types of path-structure-oriented enhancement operators to improve efficiency. One is the path tail swapping operator SwapTails, given two customer points. and ,remember and They are respectively and At subsequent customer points along their respective paths, the operator will convert the original arc... and Replace with and This is equivalent to exchanging the path between two paths. and The subsequent tail sub-paths quickly change the tail structure of both paths at a relatively low cost. The second is the SWAP* operator, which operates on two customer points. and When making an exchange, it is not required that the item must land in the opponent's original position; instead, it should first be placed in the opponent's original position. and Remove from their respective paths, and then separately Reinsertion Location path, will Reinsertion The path is defined, and the possible insertion positions within the target path are enumerated. Priority is given to the position that satisfies the capacity and time window constraints and causes the minimum distance increment. If no feasible insertion position exists, the position that minimizes the penalty value increment is selected, so that the swap and re-insertion can be completed in one operation.

[0072] S509, global knowledge sharing based on elite solutions is carried out among subpopulations.

[0073] To mitigate excessive replacements caused by frequent early improvements, the method divides the evolutionary process into early and late stages based on iteration progress, with the dividing point set at 50% of the maximum number of iterations. In this stage, the subpopulation... The current best feasible solution is compared with the global elite solution: if the current best solution is better in terms of the number of vehicles, it is directly used to update the global elite solution; only when the number of vehicles is the same, the total travel distance is further compared, with an alternative condition of at least 5% improvement in the early stage and at least 1% improvement in the later stage, to avoid frequent updates due to small fluctuations when the number of vehicles is the same. The improvement threshold is also appropriately lowered considering that the population tends to stabilize and the improvement difficulty increases in the later stage. If the global elite solution is still empty, it is directly updated. The optimal feasible solution is written as the global elite solution, where the initial global elite solution in the distance optimization phase is obtained from S10. Furthermore, the method performs a reception decision every 100 generations, comparing the global elite solution with... The current optimal feasible solution is compared, and if the global elite solution is better, it is added. This enables global elite solution propagation and parallel collaboration across subpopulations; the above update and reception processes only allow feasible solutions to participate, in order to ensure the feasibility and quality consistency of the migrated solutions.

[0074] S510, testing Has the population size reached the threshold PT? To avoid subpopulation... During the iterative process, the addition of factors and migration reception cause the population size to continuously expand. The method of this invention is... Set a subpopulation size threshold PT and perform statistics during this phase. The current subpopulation size N is determined. When N reaches PT, the process proceeds to S511 to perform population selection; when N does not reach PT, the process directly proceeds to S512 to execute subsequent steps. The threshold PT is 65.

[0075] It's worth noting that the PT is set to 65 to allow space for generated offspring and received global elite solutions. Offspring generation and elite reception cause the population to expand temporarily; retaining more solutions first helps maintain diversity and avoid premature elimination. If the PT is too small, selection will be triggered frequently, so setting it to 65 also balances computational overhead and exploration capability.

[0076] S511, for Make a selection. When S510 determines... When the subpopulation size reaches the threshold PT, for All individuals are sorted according to their fitness values ​​as defined in S503. Duplicate individuals in the population are removed first, and then the top N individuals with the best fitness are retained as a new subpopulation. .

[0077] In S512, it is determined whether the current iteration number gen is less than the stage iteration limit MaxGen. MaxGen is obtained in the main process and input as a parameter into S501. If the determination is yes, the process returns to continue executing the next round of offspring generation and evolution. If the determination is no, it proves that the iteration budget has been exhausted and the process proceeds to the subsequent termination output process.

[0078] S513, Output subpopulation Internal optimal solution. When S512 determines that the iteration budget for this stage has been exhausted, the output of this stage is... The current best individual is taken as the output result and returned to the main process.

[0079] As described above, this invention treats the Vehicle Routing Problem with Time Windows (VRPTW) as a satisfiability problem and divides the solution process into a vehicle optimization stage prioritizing the number of vehicles and a distance optimization stage under a fixed number of vehicles. This two-stage solution balances the optimization objectives of the number of vehicles used and the travel distance. To enhance the solution space exploration capability, this invention introduces a fitness penalty mechanism during the search process, penalizing infeasible solutions in terms of both capacity and time window constraints. This allows infeasible solutions to enter the population and participate in evolution, thereby increasing population diversity. To achieve rapid convergence and improve the solution efficiency of large-scale problems, this invention employs a multi-population collaborative parallel computing approach. Using the number of available computing devices K as the parallel scale parameter, a corresponding number of parallel subpopulations are constructed based on the value of K. Each subpopulation is configured with a differentiated search mode that balances neighborhood size, waiting time, and time distortion as key control factors, allowing different subpopulations to have different search focuses, thus fully utilizing computing resources. Meanwhile, during the iteration process, a cross-subpopulation migration and global knowledge sharing mechanism based on elite solutions is set up to achieve parallel collaboration and diversity maintenance, which accelerates convergence while ensuring the differentiated advantages of each subpopulation's search mode.

[0080] This invention was tested on the Gehring & Homberger benchmark standard test set, using 10 instances, each with 200 clients. The client distribution was roughly divided into two categories: a C-class test set with a concentrated client distribution, and an R-class test set with a uniformly random client distribution. Since this algorithm optimizes the number of vehicles before optimizing the distance, a smaller NV value indicates a higher solution quality; for solutions with the same NV, a smaller TD value indicates better performance. Simultaneously, this invention also conducted ablation experiments on local search of offspring based on a specific search pattern that balances neighborhood size, waiting time, and time warp, and on global knowledge sharing based on elite solutions. Parallel Wait–Warp Tuned Evolution with Elite Migration and Sharing (PWT-EMS) is the English name of this invention. PWT represents removing the global knowledge sharing part based on elite solutions, and P-EMS represents removing the specific search pattern based on a trade-off between neighborhood size, waiting time, and time warp. P-NWTEMS represents removing the specific search pattern based on a trade-off between neighborhood size, waiting time, and time warp, and global knowledge sharing based on elite solutions. Experimental results are as follows Figure 4 As shown.

[0081] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A distributed, collaborative, parallel evolutionary optimization method for green and low-carbon logistics transportation, characterized in that, Includes the following steps: S1. Obtain a vehicle routing model with time windows. The model includes the upper limit of the number of vehicles allowed, the maximum load capacity Q of each vehicle, the coordinates of the warehouse and customer points, the earliest start time window and the latest end time window of each customer point, the service time, and the algorithm parameters. With the upper limit of the number of vehicles allowed and the total driving distance as the optimization objectives, under the conditions of satisfying the vehicle capacity constraints and customer time window constraints, search for the vehicle routing scheme, i.e. the global optimal solution. The parameters include the parallel scale parameter K, the population size N, the maximum number of iterations MaxGen, and the two-stage iteration budget allocation coefficient ŋ. S2, calculate the theoretical minimum number of vehicles Lb based on the total demand D and the single vehicle capacity Q; S3 sets the maximum number of vehicles (MaxNV) during the vehicle optimization phase. S4. Construct K parallel subpopulations based on the parallel scale parameter K, and configure differentiated search mode parameters for each subpopulation, including neighborhood size M, wait time weight waitTW, and time warp weight warpTW. S5, under the MaxNV constraint, the parallel multi-subpopulation cooperative evolution of the vehicle optimization stage is initiated. Each subpopulation independently performs evolutionary operations according to its differentiated search mode, and cross-subpopulation cooperation is achieved through a global knowledge sharing mechanism based on elite solutions. S6, determine whether there is a feasible solution that simultaneously satisfies the capacity constraint and the time window constraint; If it exists, proceed to S7; otherwise, jump to S10. S7. Select the solution with the smallest number of vehicles NV from all feasible solutions, and denote its number of vehicles as bestNV. S8, determine whether the vehicle optimization phase should be terminated: if bestNV equals Lb, or the number of iterations in the vehicle optimization phase reaches the allocated budget, then terminate the vehicle optimization phase and jump to S10; Otherwise, execute S9; S9, update MaxNV to bestNV−1, and return to S5 to continue the search; S10, summarize the best results of each subpopulation in the vehicle optimization stage, select the solution with the smallest number of vehicles NV and the smallest total driving distance TD as the global optimal solution in the first stage, and record its number of vehicles as bestNV; S11, using bestNV as a fixed vehicle number constraint, initiates parallel multi-subpopulation cooperative evolution in the distance optimization phase; Under the constraint of a fixed number of vehicles, each subpopulation continues to evolve to optimize the total driving distance TD; S12, output the global optimal solution of the distance optimization stage as the final result; The time window constraint is handled using a time warp correction mechanism: when a vehicle arrives at a customer point later than its latest end time window, the service start time of that customer point is corrected to its latest end time window to prevent the cumulative impact of lateness on the time feasibility judgment of subsequent customer points.

2. The method according to claim 1, characterized in that, The configuration method for the differentiated search mode parameters in step S4 is as follows: For the k-th subpopulation (k=1,2,…,K), its neighborhood size M, waiting time weight waitTW, and time warp weight warpTW are calculated as follows: In the formula, V represents the number of customer points; N represents the population size; The neighborhood reference coefficient, ; This represents the subpopulation difference offset. .

3. The method according to claim 1, characterized in that, The evolutionary process of each subpopulation includes: (a) The subpopulation is initialized using a variety of construction strategies, including: the shortest distance strategy, the earliest latest end time window strategy, the maximum demand strategy, the distance urgency strategy, the random order strategy, the polar angle scan order strategy, and the spatial clustering grouping order strategy. (b) Calculate the fitness of an individual, wherein the fitness function is: In the formula, TD represents the total travel distance, EX represents the total capacity violation, and TW represents the total late violation calculated using the time warp correction method. and The penalty coefficient is for adaptive adjustment; (c) Detect the number of iterations in which the optimal solution of the current subpopulation has not been improved since the last update. If the number exceeds the preset stagnation threshold, reinitialize the subpopulation and return to step a; otherwise, continue to step d. (d) Use the tournament selection method to select two parent individuals from the feasible solutions first; (e) Perform crossover on the parent generation to generate offspring using the selective path exchange operator SREX; (f) Perform a local search on the offspring, the local search being based on the differentiated neighborhood set configured in step S4. Constructing a neighborhood set based on rules Iterative improvements are then made by applying a pre-defined combination of local search operators.

4. The method according to claim 3, characterized in that, The penalty coefficient and The adaptive adjustment method is as follows: Let the proportion of individuals satisfying the capacity constraint within the statistical period be . The proportion of individuals that satisfy the time window constraint is The target ratio is ,but: , in, To punish the upward adjustment factor, The penalty is to lower the factor.

5. The method according to claim 3, characterized in that, The combination of local search operators includes various types of local search operators; among them, the ExchangeUY operator is used to exchange customer points between two paths. and Starting with U and Y consecutive client point segments, different scale structural adjustments from single-point swapping to multi-client segment swapping can be achieved by setting different U and Y values; it also includes SwapTails operator and SWAP* operator. The path tail swapping operator (SwapTails) is used to swap the tail sub-paths after a specified client point in two paths, and the re-insertable swapping operator (SWAP*) is used to swap two client points and re-insert them into the optimal position of the other path.

6. The method according to any one of claims 1, 3, or 5, characterized in that, In the parallel multi-subpopulation co-evolution process of steps S5 and / or S11, cross-subpopulation cooperation is achieved through the following global knowledge-sharing mechanism based on elite solutions: Each subpopulation periodically compares its internal best feasible solution with the global elite solution: if the number of vehicles is less, or the number of vehicles is the same and the total travel distance improves by ≥5% in the early stage of evolution and ≥1% in the later stage, then the global elite solution is updated; Every preset number of generations, each subpopulation receives the current global elite solution. If the solution is better than its own optimal solution, it is added to the subpopulation.

7. The method according to claim 3, characterized in that, It also includes a population restart mechanism: when the number of consecutive times the optimal solution within a subpopulation remains unimproved exceeds the stagnation threshold ST, the subpopulation is reinitialized.

8. The method according to claim 3, characterized in that, It also includes a population size control mechanism: when the subpopulation size exceeds the threshold PT, it sorts by fitness and removes duplicate individuals and individuals with poor fitness, so that the population size is restored to N.

9. The method according to claim 1, characterized in that, The two-stage iterative budget allocation satisfies: the maximum number of iterations in the vehicle optimization stage. Minimum number of iterations from the optimization phase ; If the vehicle optimization phase terminates prematurely, its remaining iterations are automatically incorporated into the available iteration budget of the distance optimization phase.

Images