A method for planning multiple equivalent optimal paths considering actual multiple constraints

Abstract

The present invention discloses a method for planning multiple equivalent optimal paths considering actual multiple constraints, belonging to the field of transportation science and technology. Based on the actual situation in urban logistics distribution, this method considers complex road network constraints, load constraints, time window constraints, and demand split constraints that are more in line with real life, and constructs corresponding objective functions; when solving multiple equivalent optimal paths, an improved differential evolution algorithm is designed. First, in the initialization process, a generalized reverse learning strategy is used to generate a reverse customer population, and solutions in the current space and the transformed space are searched simultaneously; and three probability methods are used to select individuals for generating differential vectors in the population, adaptively balancing the exploration ability and exploitation ability of the population, and improving the diversity of the decision space and the objective space; moreover, the individuals that generate mutations are legalized instead of being directly discarded, so that the individuals beyond the solution space range can be corrected, further ensuring that multiple equivalent optimal paths can be obtained.

Description

Technical Field

[0001] The present invention relates to a method for planning multiple equivalent optimal paths considering actual multiple constraints, and belongs to the field of transportation science and technology. Background Art

[0002] In the past few decades, the Vehicle Routing Problem (VRP) and its variants have gained wide popularity due to their ability to simulate practical applications in various fields. The Rich VRP (RVRP) considering actual multiple constraints is more relevant in practical applications than the VRP. The goal of the RVRP is to find a reasonable vehicle allocation and route planning under given constraints to minimize the cost. In the RVRP, the cost can be the number of vehicles / routes, the total distance traveled by all vehicles, the total travel time, etc.

[0003] The RVRP emerges in response to actual needs. As an NP-hard problem and a multi-objective optimization problem, the importance of its objective function varies from field to field. For example, for food delivery and the medical industry, the delay time is crucial; the freight transportation industry can consider the total itinerary as the key objective to be minimized compared to other objectives, because the fuel consumption is proportional to the travel distance, so from an economic perspective, it is important to minimize the total distance traveled by all vehicles; for small industries, minimizing the number of vehicles may be the top priority compared to other objectives. When planning vehicle routes, decision-makers hope to obtain multiple paths that meet the target requirements simultaneously to ensure the stability of decisions.

[0004] VRP and its basic variants have been widely explored, and the RVRP for practical problems has become a new research trend in recent years. Although certain research results have been achieved in the study of RVRP, there are still some problems: (1) It is overlooked that the logistics distribution process is based on a complex urban road network, and the current research on RVRP has not considered the impact of the complex urban road network on vehicle routes. (2) The solution objectives are single. Currently, the research on RVRP mainly solves the objectives from a certain perspective, such as the minimum total cost, the minimum travel time, the minimum average waiting time of customers, the minimum traveltime, etc. However, in the actual logistics distribution process, multi-objective optimization needs to be considered. (3) Most of the existing datasets of RVRP either do not consider real-world characteristics or come from classical single-objective datasets, and there are no problem instances that simulate real-world characteristics. Therefore, the relevant methods for solving RVRP are also different. (4) The resource utilization rate is not high. For most logistics enterprises or distribution centers, the number of vehicles and resources such as distribution personnel available for each transportation task are limited. Therefore, how to reduce the logistics transportation cost, maximize the use of limited resources, and improve the loading rate of transportation vehicles is crucial for relevant research on route optimization.

[0005] In the prior art, Gaurav Srivastava et al. proposed a non-dominated sorting genetic algorithm with objective-specific mutation operators (NSGA-II). By establishing a multi-objective minimization model for logistics distribution costs, INSGA-II was used to optimize and solve the model (Srivastava G, Singh A, Mallipeddi R (2021). NSGA-II with objective-specific variation operators for multi-objective vehicle routing problem with time windows. Expert Systems with Applications, 176: 114779.). Caitong Yue et al. proposed a multi-modal multi-objective differential evolution algorithm using improved crowding distance (MMODE_ICD), which simultaneously ensured the diversity and convergence of the decision space and the objective space (Yue C, Suganthan P N, Liang J, et al. (2021) Differential Evolution Using Improved Crowding Distance for Multi-modal Multi-objective Optimization. Swarm and Evolutionary Computation, 62(9): 100849.). These algorithms have solved the problem of urban logistics distribution to a certain extent, but when conducting logistics distribution for urban logistics considering actual multiple constraints, they still cannot simultaneously meet the planning requirements of different decision-makers for vehicle routes. Summary of the Invention

[0006] In order to more realistically reflect the vehicle routing problem, while meeting the planning requirements of different decision-makers for vehicle routes and seeking multiple equivalent optimal routes, the present invention provides a method for planning multiple equivalent optimal routes considering actual multiple constraints. By considering the Rich VRP (RVRP) for actual multiple constraints under four constraints, namely, more realistic complex road network constraints, load constraints, time window constraints, and demand split constraints, and regarding it as a multi-modal multi-objective optimization problem, an improved differential evolution algorithm is used to achieve the planning of multiple equivalent optimal routes for urban logistics distribution.

[0007] A multi-constraint multi-equivalent optimal path planning method considering actual situations, the method comprising:

[0008] Considering actual multi-constraint conditions and constructing a corresponding objective function; the multi-constraint conditions include complex road network constraints, load constraints, time window constraints, and demand split constraints;

[0009] Using an improved differential evolution algorithm EIDSDE to solve for multiple equivalent optimal path plans for logistics distribution. In the population initialization stage of the improved differential evolution algorithm EIDSDE, a generalized reverse learning strategy is introduced to constrain the search space of the population; in the individual selection stage, individuals for generating differential vectors are selected from the population in three probabilistic ways, and the crowding distance and special crowding distance are calculated for all selected individuals. At the same time, the calculation method of the crowding distance is converted into the weighted sum of adjacent Euclidean distances; in the mutation stage, if the individual generated by mutation does not meet the boundary conditions, secondary mutation is performed. If the individual after secondary mutation still does not meet the boundary conditions, it is repaired according to a preset repair strategy to make it meet the boundary conditions; in the environmental selection stage, a predetermined percentage of the previous individuals are selected; each individual in the population represents a possible path for logistics distribution.

[0010] Optionally, the complex road network constraint refers to the constraint generated by traffic elements and the urban road network topology; the load constraint means that the vehicle shall not exceed the maximum load when loading; the time window constraint refers to the upper and lower limits of the logistics distribution time; the demand split constraint refers to the constraint on whether the goods to be distributed can be split.

[0011] Optionally, the objective functions constructed corresponding to the 4 constraints considered by the method are respectively:

[0012] f1 Number of vehicles:

[0013] f2 Total distance:

[0014] f3 Total distribution time:

[0015] f4 Total distribution cost:

[0016] Wherein, R is the number of vehicles required to complete the distribution task, q i is the demand of customer i, w is the maximum carrying capacity of the vehicle, n represents the number of customers, represents rounding up;

[0017] d ij is the distribution distance from customer point i to customer point j, indicates whether the vehicle on the r-th route passes through the arc (i, j), is a decision variable. When and only when the vehicle on the r-th route passes through the arc (i, j), otherwise It means that each customer is visited at least once. C = {0, 1, 2,..., n} represents the set of the distribution center and customers, and C' = C / {C0} represents the set of n customers;

[0018] Td j is the waiting time of the vehicle at customer j. v represents the speed of the distribution vehicle, and β is the time cost coefficient generated by violating the specified delivery time of the customer. It means whether the delivery task has a time requirement, b ir represents the actual time when vehicle r arrives at customer i, b or = 0 means the vehicle departs at time 0, LT i represents the latest time when the distribution vehicle is allowed to arrive at customer i;

[0019] FY is the cost matrix. The cost fy corresponding to each path a ij ∈A ij ∈FY, G is the fixed cost of the vehicle, l is the time delay cost, α is the time delay cost corresponding to the unit distance, S r represents the set of customers served in the r-th route, that is, the set of customers served by the r-th vehicle. |S r | represents the number of elements in the set, that is, the number of customers. (x i , y i ) represents the coordinates of customer i.

[0020] Optionally, in the method, each individual represents a possible path of logistics distribution. Corresponding to the population, each individual corresponds to a point, and the dimension of each point represents the number of customers passed by the path. The introduction of the generalized opposition-based learning strategy to constrain the search space of the population in the population initialization stage includes:

[0021] Suppose P is a candidate solution, representing a possible path. Suppose P = (z1, z2,..., z D ) is a point in a D-dimensional space, where z1, z2,..., z D ∈R and z m ∈[L m , U m , f(·) is the objective function fitness value of the candidate solution. Then the opposite point of P is where k = random(0, 1). If then If then it means has a better fitness value than P, and at this time, select to replace P, otherwise remain unchanged.

[0022] Optionally, in the individual selection stage of the method, individuals for generating differential vectors are selected from the population in three probability ways, including:

[0023] The first probability way p1: p1 = 1 - (G c - 1) / Max_gen, randomly select five neighborhoods in the entire population with a probability exceeding 0.5 as

[0024] The second probability way p2: p2 = (1 - p1) / 2. When selecting neighborhoods in the decision space, according to the average Euclidean distance between the current individual and the remaining solutions in the population, select a certain number of neighborhoods, and then randomly select 5 neighbors from these selected neighborhoods, and select the one with the largest crowding distance as The remaining four are

[0025] The third probability way p3: p3 = 1 - p1 - p2. When selecting neighborhoods in the objective space, according to the average Euclidean distance between the current individual and the remaining solutions in the population, select a certain number of neighborhoods, and then randomly select 5 neighbors from these selected neighborhoods, and then select the one with the largest crowding distance as The remaining four are

[0026] Among them, G c represents the current iteration generation number, and Max_gen represents the maximum iteration number;

[0027] Embed a niche method based on Euclidean distance in the differential evolution algorithm, and select r1, r2, r3, r4, r5 in the neighborhoods of the corresponding individuals according to the crowding distances of the corresponding individuals in the decision space or the objective space.

[0028] Optionally, in the mutation stage of the method, legalize the mutated individuals, including:

[0029] Use DE / rand / 2 for mutation, and the differential vector is generated as: where, v i represents the differential vector, r1, r2, r3, r4, r5 are unequal integers; F is the scaling factor for scaling the differential vector;

[0030] If the individual generated by mutation does not satisfy the boundary conditions, then:

[0031] Calculate the average Euclidean distance d between the first out-of-bounds individual and other individuals with mutation completed avg ;

[0032] In the solution space, select an individual x' whose Euclidean distance from the first out-of-bounds individual is d avg ;

[0033] Mutate the individual x' according to V i,m =x r1,m -F[(x r2,m -x r3,m )+(x r4,m -x r5,m )]; if the individual is still out-of-bounds after the second mutation, repair it according to ; where, v i,m represents the value of the i th -th individual in the m th -th dimension, and U m and L m represent the upper and lower bounds of the decision space

[0034] Optionally, in the environmental selection stage, the method selects a predetermined percentage of the previous individuals, including:

[0035] Iterate according to , where, G c is the current genetic generation, and Max_gen is the maximum genetic generation

[0036] When the genetic generation is 1, the front-row selection ratio is R (0 < R < 1), and when the genetic generation is G * Max_gen (1 < G * Max_gen < Max_gen), the ratio is 1

[0037] Optionally, the scaling factor F of the scaled differential vector gradually decreases as the number of iterations increases

[0038] This application also provides a vehicle scheduling method, which determines the vehicle scheduling plan according to the above-mentioned multiple equivalent optimal path planning method considering actual multiple constraints

[0039] This application also provides the application of the above-mentioned multiple equivalent optimal path planning method considering actual multiple constraints in the field of logistics distribution

[0040] The beneficial effects of the present invention are as follows

[0041] (1) By considering the logistics distribution under four constraints of vehicle load constraint, complex road network constraint, time window constraint, and demand split constraint, the present invention seeks multiple equivalent optimal paths and regards the RVRP as a multi-modal multi-objective optimization problem

[0042] (2) An improved differential evolution algorithm is proposed. First, in the initialization process, the generalized opposite learning (GOBL) strategy is used to generate the opposite customer population. The transformation search space technology is adopted to transform the solutions in the current space to a new space, that is, in the solution process, not only the candidate solutions in the current space but also the candidate solutions in the transformed space need to be considered. Since the solutions in both the current space and the transformed space are searched simultaneously, GOBL can quickly find the optimal solution. In this application, differential vector customer individuals are generated in the current space and the transformed space, the crowding distance and the special crowding distance are calculated for all the selected individuals, and at the same time, the calculation method of the crowding distance is converted to the weighted sum of adjacent Euclidean distances. Three probability methods are used to select the individuals for generating the differential vectors in the population, which can adaptively balance the exploration ability and the exploitation ability of the individuals and improve the diversity of the decision space and the target space. Then, the individuals that generate mutations are legalized, and all the processed individuals are sorted instead of being directly discarded, so that the individuals beyond the solution space range can be corrected, ensuring the diversity of the individuals and further ensuring that multiple equivalent optimal paths can be obtained. Then, the sorted individuals are adaptively selected at a certain ratio. Finally, multiple equivalent optimal paths are obtained. Brief Description of the Drawings

[0043] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the following will briefly introduce the drawings required for the description of the embodiments. Obviously, the drawings in the following description are only some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0044] Figure 1 It is a flowchart of the improved differential evolution algorithm provided in an embodiment of the present invention. Detailed Embodiment

[0045] To make the objectives, technical solutions, and advantages of the present invention clearer, the following will further describe the embodiments of the present invention in detail with reference to the drawings.

[0046] First, the basic theory involved in this application is introduced as follows:

[0047] 1. The running steps of the differential evolution algorithm include:

[0048] (1) Initialization

[0049] Let D be the dimension of the individual, NP be the population size, t be the evolution generation, and X(t) be the population of the t-th generation;

[0050] First, the 0-th generation population (i.e., the initial population) is randomly generated within the decision space of the problem: X(0) = {x1(0), x2(0),..., xNP (0)}, where x i (t) = (x i,1 (t), x i,2 (t),..., x i,D (t)) is used to represent the i-th individual in the t-th generation population.

[0051] The value of each dimension of the individual can be generated according to x i,m (0) = L m + rand i,m [0, 1](U m - L m ), where: 1 ≤ i ≤ NP, 1 ≤ m ≤ D, [L m , U m is the value range of the m-th dimension, and rand i,m [0, 1] is a uniformly distributed random number between 0 and 1.

[0052] (2) Mutation

[0053] From a biological perspective, mutation means that the gene order in the chromosome changes. In the field of evolutionary computation, mutation is regarded as a change in a certain element. In the differential evolution algorithm, taking the simplest mutation operation (DE / rand / 1) as an example, for the i th benchmark vectors , the mutation operation is: randomly select three vectors x r1 (t), x r2 (t), x r3 (t) from the current population. It is required that r1, r2, r3 are randomly selected and distinct integers from the set {1, 2,..., NP}\{i}. The difference between two vectors is scaled and added to the third vector to obtain a mutant vector V i (t): V i (t) = x r1 (t) + F(x r2 (t) - x r3 (t)), where F is the scaling factor of DE, and its value range is [0, 1].

[0054] Obviously, the smaller the difference vector between x r2 (t) and x r3 (t), the smaller the perturbation. This means that in the initial stage of the algorithm, due to the large difference between individuals, the perturbation is large, and the algorithm searches in a relatively large range. In the later stage of the algorithm, because the population is close to the optimal individual, the perturbation value is small, and the algorithm searches in a small range.

[0055] (3) Crossover

[0056] To improve the diversity of the population, the differential evolution algorithm introduces a discrete crossover operator. Different from the crossover operators in other evolutionary algorithms that exchange genes based on multiple reference vectors from the parent generation, the crossover operator in the differential evolution algorithm operates on the reference vector and the mutant vector.

[0057] (4) Selection

[0058] The selection operation of the differential evolution algorithm is a greedy selection mechanism that retains the best between the target vector and its corresponding trial vector, making the fitness value of the offspring individuals always better than that of the parent individuals. This causes the population to always evolve towards the optimal solution position and gradually focus on the optimal solution position or the satisfactory solution position.

[0059] Example 1:

[0060] This embodiment provides a method for planning multiple equivalent optimal paths considering actual multiple constraints. Refer to Figure 1 , the method includes:

[0061] (1) Consider multiple practical constraints

[0062] For urban logistics distribution, consider more practical complex road network constraints, load constraints, time window constraints, and demand split constraints, and consider multi-objective optimization problems;

[0063] Load constraint: When loading goods, the vehicle shall not exceed its maximum load.

[0064] Complex road network constraint: The complexity of the road network is mainly manifested in two aspects: detailed and complex traffic elements and complex topological structures. The complex traffic elements mainly include traffic lights, accidents, peak hours, etc., as well as the restrictions of urban road intersections and road sections; the complex urban road network topological structure is mainly reflected in that there are multiple choices for the connection paths between two points.

[0065] Time window constraint: The limitation of the time window indicates the upper and lower bounds within which customers can accept the delivery service. Delivering earlier or later than this bound requires a certain penalty cost, which is more in line with the actual delivery business scenario.

[0066] Demand split constraint: When the distribution center delivers orders to customer points, if the customer's demand allows for split delivery, it is better in terms of both the driving distance of the vehicle and the number of vehicles used than the situation where the customer's demand can only be served by one vehicle once.

[0067] (2) Construct the objective function, that is, consider the multi-objective function under the four constraints;

[0068] f1 number of vehicles:

[0069] Total distance of f2:

[0070] Total delivery time of f3:

[0071] Total delivery cost of f4:

[0072] Wherein, G is the fixed cost of the vehicle, l is the time delay cost,

[0073] n is the number of customers, (x i , y i ) represents the coordinates of customer i, R is the number of vehicles required to complete the delivery task. Since in actual delivery, one delivery route requires one vehicle, so r ∈ R, represents rounding up, w is the maximum carrying capacity of the vehicle, q i is the demand of customer i.

[0074] d ij is the delivery distance from customer point i to customer point j, means that each customer is visited at least once. C = {0, 1, 2,..., n} represents the set of the distribution center and customers, and C' = C / {C0} represents the set of n customers.

[0075] α is the time delay cost corresponding to the unit distance, and β is the time cost coefficient generated by violating the delivery time specified by the customer.

[0076] Td j is the waiting time of the vehicle at customer j, b ir represents the actual time when vehicle r arrives at customer i, LT i represents the latest time when the delivery vehicle is allowed to arrive at customer i, b or = 0 means the vehicle departs at time 0, and v represents the speed of the delivery vehicle.

[0077] S r represents the set of customers served in the r-th route, that is, the set of customers served by the r-th vehicle, |S r | represents the number of elements in the set, that is, the number of customers.

[0078] indicates whether the delivery task has time requirements.

[0079] FY is the cost matrix, and the cost fy ij ∈ A corresponding to each path a ij ∈ FY,

[0080] is a decision variable. When and only when the vehicle in the r-th route passes through the arc (i, j), Otherwise

[0081] (3) Population generation based on generalized opposition-based learning

[0082] Introduce the generalized opposition-based learning (GOBL) strategy in the initialization process of the evolutionary algorithm to constrain the search space of the population and accelerate the convergence speed and efficiency of the algorithm:

[0083] Suppose \(P=(z_1,z_2,\cdots,z\) D ) is a point in a \(D -\)dimensional space (assuming \(P\) is a candidate solution), where \(z_1,z_2,\cdots,z\) D \(\in R\) and \(z\) m \(\in[L\) m ,U m , Let \(f(\cdot)\) be the objective function fitness value of the candidate solution, then the opposite point of \(P\) is where \(k = random(0,1)\), if then if then it means has a better fitness value than \(P\). At this time, select to replace \(P\), otherwise keep it unchanged.

[0084] (4) Initial individual selection

[0085] In order to adaptively balance the exploration ability and exploitation ability of individuals, improve the diversity of the decision space and the target space, use three probability methods to select individuals for generating differential vectors in the population, and calculate the crowding distance between individuals;

[0086] Use three probability methods to select individuals for generating differential vectors in the population:

[0087] ① \(p_1\): \(p_1 = 1-(G\) c \(-1) / Max\_gen\), randomly select five neighborhoods in the entire population with a relatively high probability as In practical applications, random selection can be performed with a probability exceeding \(0.5\).

[0088] ② \(p_2\): \(p_2=(1 - p_1) / 2\). When selecting neighborhoods in the decision space, according to the average Euclidean distance between the current individual and the remaining solutions in the population, select a certain number of neighborhoods (the size of this neighborhood is determined according to the population size. In this case, the neighborhood size is set to 12), then randomly select 5 neighbors from the determined number of neighborhoods, and then select the one with the largest crowding distance as The remaining four are

[0089] ③p3: p3 = 1-p1-p2. When selecting neighbors in the target space, a certain number of neighbors are selected according to the average Euclidean distance between the current individual and the remaining solutions in the population (the meaning of "a certain number of neighbors" here is also determined according to the population size. In this case, the neighborhood size is set to 12). Then, 5 neighbors are randomly selected from the determined number of neighbors, and the one with the largest crowding distance is selected as The remaining four are

[0090] Among them, G c Represents the number of the current iteration, and Max_gen represents the maximum number of iterations.

[0091] A niche method based on Euclidean distance is embedded in the differential evolution algorithm to select r1, r2, r3, r4, and r5 in the neighborhood of the corresponding individuals according to the crowding distance of the corresponding individuals in the decision space or the target space.

[0092] (5) Variation and legalization

[0093] Using DE / rand / 2 for mutation, the difference vector is generated as: Among them, v i represents the difference vector, r1, r2, r3, r4, r5 are unequal integers. F is the scaling factor used to scale the difference vector. The mutated vector individuals may be infeasible, that is, the mutated individuals do not meet the boundary conditions and directly exceed the solution space range, that is, fall outside the search space, so it is necessary to legalize the infeasible solution.

[0094] (5.1) Calculate the average Euclidean distance d between the first crossing individual and other mutation-completed individuals avg ;

[0095] (5.2) In the solution space, the Euclidean distance from the first individual crossing the boundary is d avg Individual x';

[0096] (5.3) Individual x' is transformed into i,m =x r1,m -F[(x r2,m -x r3,m )+(x r4,m -x r5,m )] to perform secondary mutation;

[0097] If the individual still crosses the boundary after the second mutation, Make repairs;

[0098] Among them, vi,m Denote the value of the \(i\)-th individual in the \(m\)-dimensional space as \(U\) and \(L\) represents the upper and lower bounds of the decision space. th -th individual in the \(m\)-dimensional space, where \(U\) and \(L\) represent the upper and lower bounds of the decision space. th -dimensional value of the \(i\)-th individual, and \(U\) and \(L\) represent the upper and lower bounds of the decision space. m and \(L\) m represent the upper and lower bounds of the decision space.

[0099] (6) Individual crossover

[0100] (7) Environmental selection operation

[0101] In the new environmental selection method, only a certain percentage of the prefix elements are selected.

[0102] The ratio of the selected individuals is iterated according to , where \(G\) is the current generation number of the genetic algorithm and \(Max\_gen\) is the maximum generation number of the genetic algorithm. When the generation number \(G = 1\), the front-row selection ratio is \(R(0\lt R\lt1)\); when the generation number \(G = G\times Max\_gen(1\lt G\times Max\_gen\lt Max\_gen)\), the ratio is \(1\). c where \(G\) is the current generation number of the genetic algorithm and \(Max\_gen\) is the maximum generation number of the genetic algorithm. When \(G = 1\), the front-row selection ratio is \(R(0\lt R\lt1)\); when \(G = G\times Max\_gen(1\lt G\times Max\_gen\lt Max\_gen)\), the ratio is \(1\).

[0103] After initialization, mutation, crossover, and environmental selection, multiple equivalent optimal paths for urban logistics distribution are obtained.

[0104] The effects of the present invention can be further illustrated by the following simulation experiments.

[0105] 1. Simulation conditions

[0106] This embodiment conducts simulation experiments using MATLAB 2020a, and the experimental device is a computer with an i7-7500 2.90GHz CPU. Since the Solomon dataset (the Solomon dataset is a publicly available dataset) cannot represent the ideal benchmark scenario of the RVRP proposed in the present invention, and the time window constraint has a certain dominant position in the RVRP proposed in the present invention, the present invention uses the same test dataset as Zhou and Wang (Zhou Y, Wang J. A Local Search-Based Multi-objective Optimization Algorithm for Multi-objective Vehicle Routing Problem With Time Windows(2017). IEEE Systems Journal, 9(3):1100-1113.) and Castro Gutierrez et al. (JPC Gutiérrez, Landa-Silva D, JA Moreno-Pérez. Nature of real-world multi-objective vehicle routing with evolutionary algorithms(2011). IEEE International Conference on Systems. IEEE, 257-264.). This dataset consists of 45 real VRP data instances, which are composed of customers of 3 different scales, 5 time window configurations, and 3 vehicle models with different capacities, that is, 3 groups * 5 time windows * 3 vehicle capacities = 45 VRP instances. This dataset and configuration file can be downloaded from https: / / github.com / psxjpc / .

[0107] This embodiment compares the method EIDSDE of the present application with the existing INSGA_II and MMODE_ICD. These two methods can be referred to:

[0108] INSGA_II can be referred to "Srivastava G, Singh A, Mallipeddi R. NSGA-II with objective-specific variation operators for multi-objective vehicle routing problem with time windows[J]. Expert Systems with Applications, 2021, 176:114779."

[0109] MMODE_ICD can refer to "Yue C, Suganthan P N, Liang J, et al. Differential Evolution Using Improved Crowding Distance for Multi-modal Multi-objective Optimization[J]. Swarm and Evolutionary Computation, 2021, 62(9): 100849."

[0110] 2. Simulation Results

[0111] A single performance metric cannot comprehensively measure the performance of a multi-objective optimization algorithm. Therefore, four metrics are used in this invention:

[0112] (1) Inverted generational distance (IGD): The smaller the IGD value, the closer the non-dominated solution set is to the true Pareto front and the more evenly distributed it is, indicating better convergence and diversity of the solution set.

[0113] (2) Coverage metric (C-metric): Its value is calculated by C(X, Y), where C(X, Y) represents the proportion of solutions in solution set Y that are dominated by a certain solution in solution set X.

[0114] (3) 1 / HV: HV (hyper volume) is the hyper volume. 1 / HV can measure performance in the decision space, and the smaller the value, the better the performance.

[0115] (4) 1 / PSP: PSP is the Pareto Sets Proximity (PSP). PSP reflects the overlap rate and distance between the true PS and the obtained PS. 1 / PSP can measure performance in the objective space, and the smaller the value, the better the performance.

[0116] "Number1-Number2-Number3" represents the name of the instance, where Number1 represents the number of customers, Number2 represents different types of indices regarding vehicle capacity, and Number3 represents the index of time window configuration. 'B / S / W' means that the proposed algorithm is significantly better than / essentially similar to / significantly worse than the current algorithm.

[0117] In the following Table 1 and Table 2, the abbreviated meanings are as follows:

[0118] C(IN, EI): C(INSGA-II, EIDSDE)

[0119] C(EI,IN): C(EIDSDE, INSGA-II)

[0120] C(MM, EI): C(MMODE_ICD, EIDSDE)

[0121] C(EI, MM): C(EIDSDE, MMODE_ICD)

[0122] Table 1 Comparison of the average values of four indicators between EIDSDE and INSGA_II

[0123]

[0124] Table 2 Comparison of the average values of four indicators between EIDSDE and MMODE_ICD

[0125]

[0126] It can be seen from Table 1 and Table 2 that:

[0127] In terms of IGD, the values of the method EIDSDE of the present application in all 45 instances are significantly smaller than those of INSGA-II, and the values in 32 instances are smaller than those of MMODE_ICD. And the smaller the IGD value, the closer the non-dominated solution set is to the true Pareto front and the more uniform the distribution, and the better the convergence and diversity of the solution set.

[0128] In terms of C-metric, the values of the method EIDSDE of the present application in 30 instances are significantly better than those of INSGA-II, and the values in 28 instances are significantly better than those of MMODE_ICD.

[0129] In terms of 1 / HV and 1 / PSP, the method EIDSDE of the present application also shows good effects.

[0130] It can be obviously seen from Table 1 and Table 2 that the difficulty of the problem increases with the increase in the number of customers and the decrease in vehicle capacity. The reason is that the problem with more customers and smaller capacity vehicles will have more route planning solutions. Therefore, it will be more difficult to converge on all objectives. Moreover, any solution with the minimum value in any one objective is a non-dominated solution, regardless of the values of other objectives.

[0131] In a real - world scenario, the preference for one target may be higher compared to other targets. In the RVRP with 4 targets proposed in the present invention, from a certain perspective, the target f2 (total travel distance) may be more important than other targets because the travel distance is proportional to the fuel consumption, and thus, it has a direct impact on environmental pollution. In addition, for a logistics company, the targets f4 (total cost) and f3 (total travel time) are equally important because it is necessary to ensure service to customers within the specified time to improve customer satisfaction and also to ensure reduction of the total distribution cost during the entire logistics distribution process. Therefore, it is crucial for the RVRP to find the optimal value among all targets through a method because in most cases, the decision - maker's preference is a priori unknown. Table 3 analyzes this.

[0132] Table 3 Comparison of 30 experiments between EIDSDE and MMODE_ICD / INSGA_II on four targets

[0133]

[0134] Example Two

[0135] In this example, logistics distribution is carried out in a certain area on the dataset selected in the simulation experiment of the above - mentioned Example One. Among them, the number of customers is 50, the vehicle load is 200 tons, the vehicle traveling speed is 60 km / h, the vehicle fixed cost is 500 yuan, the unit transportation cost is 40 yuan / km, the earliest service time that the customer can accept can be 5 - 15 minutes in advance, and the latest service time that the customer can accept can be postponed by 5 - 15 minutes.

[0136] Multiple equivalent optimal paths obtained by using the method of the present invention are as follows:

[0137] 0→1→7→9→15→20→25→33→36→39→45→0

[0138] 0→2→10→14→26→30→34→37→41→49→50→0

[0139] 0→3→6→11→17→23→27→31→38→42→48→0

[0140] 0→4→8→12→18→21→24→28→35→44→47→0

[0141] 0→5→13→16→19→22→29→32→40→43→46→0

[0142] Or

[0143] 0→1→6→9→17→20→26→30→34→37→41→45→0

[0144] 0 → 2 → 10 → 15 → 21 → 25 → 31 → 36 → 39 → 49 → 50 → 0

[0145] 0 → 3 → 7 → 14 → 19 → 23 → 27 → 33 → 38 → 42 → 48 → 0

[0146] 0 → 4 → 11 → 12 → 18 → 24 → 28 → 35 → 44 → 47 → 0

[0147] 0 → 5 → 8 → 13 → 16 → 22 → 29 → 32 → 40 → 43 → 46 → 0

[0148] or

[0149] 0 → 1 → 9 → 14 → 17 → 23 → 28 → 35 → 39 → 44 → 47 → 0

[0150] 0 → 2 → 6 → 12 → 18 → 24 → 27 → 33 → 38 → 42 → 48 → 0

[0151] 0 → 3 → 7 → 11 → 20 → 25 → 29 → 31 → 36 → 45 → 50 → 0

[0152] 0 → 4 → 8 → 13 → 16 → 22 → 30 → 34 → 37 → 41 → 49 → 0

[0153] 0 → 5 → 10 → 15 → 19 → 21 → 26 → 32 → 40 → 43 → 46 → 0

[0154] The minimum number of vehicles required under the equivalent optimal path is 5, the shortest path length is 264.36 kilometers, the minimum cost is 9,866.65 yuan, and the shortest total travel time is 4.41 hours (excluding the service time of the vehicles). The above are the paths for the required 5 vehicles, where the numbers represent customer nodes.

[0155] The unique path obtained using MMODE_ICD is:

[0156] 0 → 1 → 7 → 10 → 15 → 25 → 33 → 38 → 45 → 49 → 0

[0157] 0 → 2 → 11 → 14 → 20 → 26 → 30 → 34 → 41 → 0

[0158] 0 → 3 → 6 → 12 → 18 → 23 → 31 → 37 → 42 → 48 → 0

[0159] 0 → 4 → 8 → 13 → 21 → 27 → 35 → 40 → 44 → 0

[0160] 0 → 5 → 16 → 22 → 29 → 32 → 46 → 39 → 47 → 0

[0161] 0 → 9 → 17 → 19 → 24 → 28 → 36 → 43 → 50 → 0

[0162] The minimum number of vehicles required for the optimal path is 6, the shortest path length is 375.82 kilometers, the minimum cost is 122,045.64 yuan, and the shortest total travel time is 6.26 hours (excluding the service time of the vehicles).

[0163] The only path obtained by using INSGA_II is as follows:

[0164] 0 → 1 → 7 → 25 → 33 → 45 → 0

[0165] 0 → 2 → 14 → 26 → 34 → 41 → 0

[0166] 0 → 3 → 6 → 23 → 31 → 42 → 0

[0167] 0 → 4 → 8 → 21 → 35 → 44 → 0

[0168] 0 → 5 → 16 → 22 → 32 → 46 → 0

[0169] 0 → 9 → 19 → 24 → 36 → 50 → 0

[0170] 0 → 10 → 18 → 28 → 40 → 43 → 0

[0171] 0 → 11 → 15 → 30 → 37 → 48 → 0

[0172] 0 → 12 → 17 → 27 → 38 → 49 → 0

[0173] 0 → 13 → 20 → 29 → 39 → 47 → 0

[0174] The minimum number of vehicles required for the optimal path is 10, the shortest path length is 653.47 kilometers, the minimum cost is 23,568.76 yuan, and the shortest total travel time is 10.89 hours (excluding the service time of the vehicles).

[0175] The experimental results show that the method proposed in the present invention can meet the planning requirements of different decision-makers for vehicle paths, and at the same time satisfy the simultaneous optimization of multiple objectives, obtaining multiple equivalent optimal paths, while the existing MMODE_ICD and INSGA_II can only obtain a unique optimal path and cannot meet the planning requirements of different decision-makers for vehicle paths at the same time.

[0176] Some steps in the embodiments of the present invention can be implemented by software, and the corresponding software program can be stored in a readable storage medium, such as an optical disc or a hard disk, etc.

[0177] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A method for planning multiple equivalent optimal paths considering actual multiple constraints, characterized in that, The method includes: Considering practical multi-constraints and constructing a corresponding objective function; the multi-constraints include complex road network constraints, load constraints, time window constraints, and demand split constraints; Using the improved differential evolution algorithm EIDSDE to solve for multiple equivalent optimal path plans for logistics distribution. In the population initialization stage of the improved differential evolution algorithm EIDSDE, a generalized reverse learning strategy is introduced to constrain the search space of the population. In the individual selection stage, individuals for generating differential vectors are selected from the population in three probabilistic ways, and the crowding distance and special crowding distance are calculated for all selected individuals. At the same time, the calculation method of the crowding distance is converted into the weighted sum of adjacent Euclidean distances. In the mutation stage, if the individual generated by mutation does not meet the boundary conditions, secondary mutation is performed. If the individual after secondary mutation still does not meet the boundary conditions, it is repaired according to a preset repair strategy to meet the boundary conditions. In the environmental selection stage, a predetermined percentage of the previous individuals are selected; each individual in the population represents a possible path for logistics distribution; The corresponding objective functions constructed under the four constraints considered in the method are: f1 Number of vehicles: Total distance of f2: f3 Total delivery time: Total distribution cost of f4: where \(R\) is the number of vehicles required to complete the delivery task, \(q\) i is the demand of customer \(i\), \(w\) is the maximum carrying capacity of the vehicle, \(n\) represents the number of customers, denotes rounding up; d ij The delivery distance from customer point i to customer point j Indicates whether the vehicle on the r-th route passes through the arc (i, j). r ∈ R, i, j ∈ C' are decision variables. If and only if the vehicle on the r-th route passes through the arc (i, j), otherwise Indicates that each customer is visited at least once. C = {0, 1, 2,..., n} represents the set of the distribution center and customers, and C' = C / {C0} represents the set of n customers; Td j is the waiting time of the vehicle at customer j, v represents the speed of the delivery vehicle, and β is the time cost coefficient incurred for violating the delivery time stipulated by the customer. indicates whether the delivery task has a time requirement, b ir represents the actual time when vehicle r arrives at customer i, b or = 0 indicates that the vehicle departs at time 0, LT i represents the latest time when the delivery vehicle is allowed to arrive at customer i; FY is the cost matrix, and for each path a ij ∈A, the corresponding cost fy ij ∈FY. G is the fixed cost of the vehicle, l is the time delay cost, α is the time delay cost per unit distance, and S r represents the set of customers served in the r-th route, that is, the set of customers served by the r-th vehicle. |S r | represents the number of elements in the set, that is, the number of customers. (x i ,y i ) represents the coordinates of customer i; In the method, each individual represents a possible path for logistics distribution. Corresponding to the population, each individual corresponds to a point, and the dimension of each point represents the number of customers passed by the path. The introduction of a generalized reverse learning strategy to constrain the search space of the population in the population initialization stage includes: Suppose \(P\) is a candidate solution, representing a possible path. Suppose \(P=(z_1,z_2,\cdots,z\) D ) is a point in a \(D\)-dimensional space, where \(z_1,z_2,\cdots,z\) D \(\in R\) and \(z\) m \(\in [L\) m ,U m , If \(f(\cdot)\) is the objective function fitness value of the candidate solution, then the reverse point of \(P\) is where \(k = random(0,1)\). If then If then it means has a better fitness value than \(P\). At this time, select to replace \(P\), otherwise remain unchanged; In the individual selection stage of the method, individuals for generating differential vectors are selected from the population in three probabilistic ways, including: The first probability method p1 is p1 = 1 - (G c - 1) / Max_gen, and five neighborhoods in the entire population are randomly selected with a probability exceeding 0.5 as The second probability method p2 is p2 = (1 - p1) / 2. When selecting neighborhoods in the decision space, a certain number of neighborhoods are selected according to the average Euclidean distance between the current individual and the remaining solutions in the population. Then, 5 neighbors are randomly selected from the selected number of neighborhoods, and the one with the largest crowding distance is selected as The remaining four are The third probability method p3 is p3 = 1 - p1 - p2. When selecting neighborhoods in the target space, a certain number of neighborhoods are selected according to the average Euclidean distance between the current individual and the remaining solutions in the population. Then, 5 neighbors are randomly selected from the determined number of neighborhoods, and then the one with the largest crowding distance is selected as The remaining four are Among them, G c represents the generation number of the current iteration, and Max_gen represents the maximum number of iterations; Embedding a niche method based on Euclidean distance in the differential evolution algorithm, and selecting r1, r2, r3, r4, r5 in the neighborhood of the corresponding individual according to the crowding distance of the corresponding individual in the decision space or the objective space; In the mutation stage of the method, the mutated individuals are legalized, including: Mutation is performed using DE / rand / 2, and the difference vector is generated as follows: where v i represents the difference vector, r1, r2, r3, r4, r5 are unequal integers; F is the scaling factor used to scale the difference vector; If the individual generated by mutation does not meet the boundary conditions, then: Calculate the average Euclidean distance d between the first out-of-bounds individual and other individuals with completed mutations avg ; In the solution space, an individual x' with an Euclidean distance of d from the first out-of-bounds individual is taken; avg ​ Mutate the individual x' according to V i,m = x r1,m - F[(x r2,m - x r3,m ) + (x r4,m - x r5,m )]; If the individual is still out of bounds after the second mutation, repair it according to ; Among them, v i,m represents the value of the i th -th individual in the m th -th dimension, and U m and L m represent the upper and lower bounds of the decision space.

2. The method according to claim 1, wherein The complex road network constraint refers to the constraint generated by traffic elements and the urban road network topology; the load constraint means that the vehicle shall not exceed the maximum load when loading; the time window constraint refers to the upper and lower limits of the logistics distribution time; the demand split constraint refers to the constraint on whether the goods to be distributed can be split.

3. The method according to claim 2, wherein In the environmental selection stage of the method, a predetermined percentage of the previous individuals are selected, including: When the generation number of the genetic algorithm is 1, the front row selection ratio is R, and 0 < R < 1. When the generation number of the genetic algorithm is G*Max_gen, and 1 < G*Max_gen < Max_gen, the ratio is 1.

4. The method according to claim 3, wherein The scaling factor F of the scaled differential vector gradually decreases with the increase of the number of iterations.

5. A vehicle scheduling method, characterized in that, The method determines the vehicle scheduling plan according to the multiple equivalent optimal path planning method considering actual multi-constraints described in any one of claims 1-4.

6. Application of the multiple equivalent optimal path planning method considering actual multi-constraints described in any one of claims 1-4 in the field of logistics distribution.

Images