A vehicle routing optimization method based on adaptive multi-constrained differential evolution algorithm

Abstract

The present invention aims at the multi-constraint and multi-objective vehicle path optimization problem and designs a vehicle path optimization method based on an adaptive multi-constraint differential evolution algorithm. The number of vehicles used and the total driving distance are taken as optimization objectives, a multi-objective differential evolution algorithm is adopted, multiple evolution mechanisms are adaptively adjusted, a differential mutation strategy guided by feasible solutions is utilized to improve the convergence of the solution set, an evaluation mechanism based on constraint dominance criteria is designed for selection operations, the diversity of optimal solutions is improved, the best vehicle path optimization scheme is obtained, the enterprise transportation cost is reduced, and the logistics distribution efficiency is increased.

Description

Technical Field

[0001] Aiming at the problem of vehicle path optimization from logistics center to customers in the field of supply chain logistics, the present invention designs a vehicle path optimization method based on an adaptive multi-constraint differential evolution algorithm, obtains a multi-objective vehicle path optimization scheme under multiple constraints, achieves simultaneous optimization of the number of vehicles and the total driving distance of vehicles, effectively improves the logistics distribution efficiency, and reduces the logistics transportation cost of enterprises, which belongs to both the field of intelligent information and the field of logistics. Background Art

[0002] Vehicle path optimization has important practical significance for saving the logistics and transportation costs of enterprises and promoting the development of the logistics and transportation industry. With the expansion of the scale of the logistics and distribution industry, the vehicle path problem has gradually shown the characteristics of multiple objectives and multiple constraints in the supply chain logistics scenario. First, in order to improve resource utilization, it is necessary to take into account the number of vehicles used and the driving distance, and keep the two objectives optimal at the same time. Secondly, in addition to the vehicle load limit, the vehicle path optimization in the specific logistics and transportation scenario is subject to more constraints, such as the need to deliver goods within the service time window specified by the customer. Therefore, under multiple constraints, how to determine the number of vehicles used and coordinate the distribution routes of each vehicle, and obtain a vehicle scheduling plan that meets all constraints, is a major challenge in the field of supply chain logistics. The adaptive multi-constraint differential evolution algorithm has good convergence and diversity, and can efficiently solve multi-constraint and multi-objective optimization problems, providing a new method for solving the above-mentioned vehicle path optimization problem. Therefore, the present invention has broad application prospects in vehicle path optimization problems.

[0003] The present invention proposes a vehicle path optimization method based on an adaptive multi-constraint differential evolution algorithm. The method takes the number of vehicles and the total driving distance of vehicles as optimization targets, considers multiple constraints, adaptively adjusts multiple evolution strategies, optimizes the customer service sequence, and obtains the optimal vehicle path solution, providing an effective method for improving logistics distribution efficiency and service quality. Summary of the invention

[0004] The present invention obtains a vehicle path optimization method based on an adaptive multi-constraint differential evolution algorithm. The method takes minimizing the number of vehicles and the total driving distance as the optimization goal, determines a multi-constraint vehicle path optimization model, uses feasible solutions to guide the differential mutation strategy, and uses high-quality infeasible solutions to improve the diversity of the solution set, so as to obtain the optimal path plan, improve the utilization rate of logistics transportation resources, and save the logistics costs of enterprises.

[0005] The present invention adopts the following technical solutions and implementation steps:

[0006] 1. A vehicle path optimization method based on an adaptive multi-constrained differential evolution algorithm, characterized in that the optimization target of the vehicle path optimization problem is determined, and an adaptive multi-constrained differential evolution algorithm is used to obtain a vehicle scheduling solution, solve the problem that it is difficult to obtain a feasible solution and the diversity of the solution set is poor in the vehicle path optimization process under the customer time window constraint, and achieve the simultaneous optimization of the number of vehicles and the total driving distance, including the following steps:

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

[0008] Assume that K max All transport vehicles of the same model are parked in a logistics center, ready to provide one-way delivery services to N customers. The vehicle load limit is Q, and the logistics center opening and closing times are L and E respectively; each customer corresponds to an order data, including customer number i, customer location coordinates (x i ,y i ), the cargo weight q of the i-th customer i , the earliest service start time for the i-th customer l i , the latest service start time e for the i-th customer i and the service time s of the i-th customer i , where i = 0 represents the logistics center; 0 <K max <50,0 <N<200,0<Q<300,0<L<200,0<E<200,i=0,…,N,0<x i <100,0 <y i <100,0 i <300,0 <l i <1000,0 <e i <1440,0 i <200;

[0009] The optimization goal of the vehicle routing optimization problem is to minimize the following objective function:

[0010] f 1 =K (1)

[0011]

[0012] Among them, f 1 The total number of vehicles used in the distribution process K, f 2 is the total distance traveled by the vehicle; 0 <K<36,d ij is the distance between customer i and customer j, x ijk Indicates whether the kth vehicle travels from i to j, x ijk =1 for yes, x ijk =0 means no, j=0,…,N, k=1,…,K; ​​

[0013] Each vehicle must meet the following conditions during the delivery process, and each customer will only be delivered by one vehicle:

[0014]

[0015]

[0016] Each vehicle departs from the logistics center and returns to the logistics center after completing the delivery task according to the corresponding customer service order:

[0017]

[0018]

[0019] The total customer demand for the kth car should not exceed Q:

[0020]

[0021] The time a that the vehicle arrives at the location of customer i i It should be within its service time window, and the vehicle should return to the logistics center before the logistics center closes:

[0022] l i ≤a i ≤e i (8)

[0023] L≤a 0 ≤E (9)

[0024] The time a for a vehicle to pass through i and reach j j for:

[0025]

[0026] Among them, w i ≥0 indicates the time the vehicle waits for customer i, t 0j represents the time it takes to travel from the logistics center to j, t ij represents the time it takes to travel from i to j, t i0 represents the time it takes to travel from i to the logistics center;

[0027] Step 2: Get the distance and time data for the vehicle routing optimization problem:

[0028] Calculate the distance between customer i and customer j:

[0029]

[0030] Among them, x i and i Represents the horizontal and vertical coordinates of customer i’s location, xj and j Represent the horizontal and vertical coordinates of the customer j's position respectively; according to the average vehicle speed v, calculate the corresponding to d ij The travel time t ij :

[0031]

[0032] Step 3: Design an adaptive multi-constrained differential evolution algorithm:

[0033] Step 3.1 Parameter setting and population initialization:

[0034] The maximum evolutionary number of the multi-objective differential evolution algorithm is set to T, where T is an integer in the range of [250, 2000], the mutation rate is F, 0≤F≤1, the crossover rate is Cr, 0≤Cr≤1, the population size is H, H is an integer in the range of [10, 500], the neighborhood size is c, c is an integer in the range of [1, N], each solution represents a vehicle scheduling plan, and the initial population P is randomly generated. 0 :

[0035]

[0036] Among them, p b,0 is the b-th solution of the initial population, b=1,…,H;

[0037] Step 3.2 Determination of the number of feasible solutions:

[0038] The tth generation population P t It is expressed as:

[0039]

[0040]

[0041] Among them, p b,t is the bth solution of the tth generation population, p bu,t is the u-th dimension component of the b-th solution of the t-th generation population; t=0,…,T,p bu,t =1,…,N+K max ,u=1,…,N+K max ;p bu,t =1,…,N indicates the customer, p bu,t =N+1,…,N+K max represents a logistics center; between two logistics centers, the customer order is sorted in ascending order according to the time when the customer starts service. The time when the customer starts service is read from the database, and K delivery routes r are obtained. k , k is the route number, k=1,…,K,K≤K max ;

[0042] Compute the constraint violation δ for each solution:

[0043]

[0044] Among them, q i is the cargo weight of the i-th customer, K is the total number of routes, i.e. the total number of vehicles used, Q is the vehicle load limit, and a i is the arrival time of the ith customer, which is calculated by formula (10), e i is the latest service start time for the i-th customer, which is read from the database; set the solution satisfying δ=0 as a feasible solution, and the solution satisfying δ>0 as an infeasible solution, and count P t The total number of feasible solutions in ;

[0045] Set P t The set of all feasible solutions in is F t , β represents the total number of feasible solutions, which is an integer in the range [0,H]. If β=0, execute steps 3.3-3.4 and then step 3.6. If β≥2, execute steps 3.4-3.5 and then step 3.6. If β=1, directly execute step 3.6.

[0046] Step 3.3 Differential evolution based on DE / rand / 1 mutation strategy:

[0047]

[0048]

[0049]

[0050] Among them, v b,t is the b-th mutation solution of the t-th generation population, p rand1,t is the rand of the t-th generation population 1 solution, p rand2,t is the rand of the t-th generation population 2 solution, rand 1 and rand 2 are two different integers randomly selected from [1,H] that are different from b; p bu,t+1 is the u-th dimension component of the b-th solution in the t+1-th generation population, rand bu is the random number taken by the b-th solution on the i-th dimension component, 0≤rand bu ≤1;

[0051] P b,t+1 Each dimension is rounded down:

[0052]

[0053] in, Indicates p b,t+1 The value of the u-th dimension component rounded down;

[0054] If N+1≤p bu,t+1 ≤N+K max ,N+1≤p b(u+1),t+1 ≤N+K max , and satisfy 1≤p b(u+2),t+1 ≤N, that is, p b,t+1 When there is an empty route, p b(u+1),t+1 Remove the route group, otherwise set u = u + 1 and continue to check empty routes until p b,t+1 Delete all empty routes in the .

[0055] Step 3.4 Fitness value calculation:

[0056] Calculate P t and P t+1 The dominance strength of each solution in the candidate solution set is:

[0057] If p is a feasible solution and q is an infeasible solution, then p dominates q, and the dominance strength of p increases by 1;

[0058] If p and q are both feasible solutions, f 1 (p)≤f 1 (q), f 2 (p)≤f 2 (q), then p dominates q, and the dominance strength of p increases by 1;

[0059] If both p and q are infeasible solutions, δ(p)<δ(q), then p dominates q, and the dominance strength of p increases by 1;

[0060] Among them, p and q are two different solutions in the candidate solution set, f 1 (p) and f 1 (q) are the total number of vehicles corresponding to p and the total number of vehicles corresponding to q, respectively, and f 2 (p) and f 2 (q) are the total vehicle travel distance corresponding to p and the total vehicle travel distance corresponding to q, respectively. δ(p) and δ(q) are the constraint violation degree corresponding to p and the constraint violation degree corresponding to q, respectively.

[0061] Calculate P t and P t+1 The rank of each solution in the candidate solution set is:

[0062] For the solution with δ = 0, i.e., a feasible solution, its rank is set to 1;

[0063] For solutions with δ>0, i.e. infeasible solutions, they are sorted in ascending order according to the size of δ, and the rank of each infeasible solution is set to the sum of the total number of feasible solutions and the corresponding δ ranking;

[0064] Add the dominance strength and the level value of each solution to obtain the fitness value of each solution;

[0065] Step 3.5: Differential mutation operation guided by feasible solution:

[0066] ① According to the fitness value, F t Sort in ascending order to get F t =(p 1 ,...,p β );

[0067] ② Calculate p 1 and p 2 The similarity S is:

[0068]

[0069] Among them, r p1,rand and r p2,rand Respectively represent from p 1 A delivery route randomly selected from p 2 A delivery route randomly selected from max ,|r p1,rand | and |r p2,rand |represents r p1,rand and r p2,rand The number of customers, p1,rand ∩r p2,rand | indicates that at the same time p1,rand and r p2,rand The number of customers that appear in |r p1,rand ∩r p2,rand |≥0;

[0070] ③ Adjust F t The order:

[0071] p 2,t Adjust to p β , and (p 3 ,...,p β ) is adjusted to (p 2 ,...,p β-1 );

[0072] ④If S≤0.5, execute step 3.5⑤; if S>0.5, return to step 3.5②;

[0073] ⑤Select p in the order of route number 1 A route p1,k1 , k1 Yes 1 The corresponding route number, 1≤k 1 ≤K max , set c(1,r p1,k1 ) indicates p 1 k 1 The first customer of the delivery route, in p 2 Find c(1,r p1,1 )Route p2,k2 , k 2 Yes 2 The corresponding route number, 1≤k 2 ≤K max , determine r p1,k1 and r p2,k2 The difference vector pair of :

[0074]

[0075] Among them, r p1,k1 -r p1,k1 ∩r p2,k2 Indicates that p1,k1 Rather than p2,k2 The complement of the delivery routes on p2,k2 -r p1,k1 ∩r p2,k2 Indicates that p2,k1 Rather than p1,k2 The complement of the delivery routes on the p1,k1 -r p1,k1 ∩r p2,k2 and r p2,k2 -r p1,k1 ∩r p2,k2 first customer;

[0076] ⑥Let the bth solution p of the t+1th generation population b,t+1 =p 1 , in p 1 Find the route r where n is located p1,k3 , k 3 Yes 1 Different from k 1 Route number, 1≤k 3 ≤K max , k 3 ≠k 1 , change m and n from r p1,k1 and r p1,k3 Removed:

[0077] Under the premise of satisfying formula (3)-formula (10), find r p1,k3 All positions of m can be inserted so that δ(p b,t+1) = 0 holds, calculate the route increments corresponding to all insertion positions, set the position corresponding to the minimum route increment as the insertion position of m, and insert m at this position;

[0078] On the premise of satisfying formulas (3) - (10), find all positions in r where n can be inserted such that δ(p p1,k3 ) = 0 holds, calculate the route increments corresponding to all insertion positions, set the position corresponding to the minimum route increment as the insertion position of n, and insert n at this position; b,t+1 ) = 0 holds, calculate the route increments corresponding to all insertion positions, set the position corresponding to the minimum route increment as the insertion position of n, and insert n at this position;

[0079] If there is no insertable position that makes δ(p b,t+1 ) = 0 hold, set m or n as the first customer of r p1,k3 or r p1,k1 to complete the construction of p b,t+1 ;

[0080] ⑦ Let k 1 = k 1 + 1, count the size |P t+1 | of the (t + 1)-th generation population. If k 1 = K 1 , where K 1 is the number of vehicles corresponding to p 1 , and |P t+1 | = H, where H is the population size, then go to step 3.6. If k 1 = K 1 , and |P t+1 | < H, then return to step 3.5②; If k 1 < K 1 , and |P t+1 | < H, then return to step 3.5⑤;

[0081] Step 3.6 Large-scale neighborhood search operation:

[0082] Randomly select a one-dimensional component p b,w = i, w = t or t + 1, 1 < p bu,w < N in p bu,w . According to the distance d ij between customer i and customer j, select c customers closest to customer i, where c is the neighborhood size. Find the corresponding c-dimensional component p b,w in p bv,w , 1 < p bv,w < N, v ≠ u, and remove p bv,w and the selected c-dimensional component from p b,t ;

[0083] On the premise of satisfying formulas (3) - (10), find all positions in p b,w where p bu,w can be inserted.at all positions such that δ(p b,w ) = 0 holds, calculate the route increments for all insertion positions, and set the position corresponding to the minimum route increment as p bu,w 's insertion position, and insert p bu,w into this position;

[0084] Under the premise of satisfying formulas (3) - (10), find all positions in p b,w where p bv,w can be inserted, such that δ(p b,w ) = 0 holds, calculate the route increments for all insertion positions, and set the position corresponding to the minimum route increment as p bv,w 's insertion position, and insert p bv,w into this position;

[0085] Step 3.7 Selection operation:

[0086] If P t+1 is empty, that is, there is no (t + 1)-th generation population, go to Step 3.8; otherwise, sort the candidate solution set in ascending order according to the fitness values of the solutions, and select the top H solutions to construct P t+1 ;

[0087] Step 3.8 Judge the termination condition of the algorithm:

[0088] If t < T, then let t = t + 1, and repeat Steps 3.2 - 3.7; if t = T, then stop the calculation and output the optimal vehicle scheduling plan.

[0089] The innovation of the present invention lies in:

[0090] (1) In view of the problem that it is difficult to obtain the optimal solution when considering multiple constraint conditions in the vehicle path optimization process, the present invention uses a multi-objective differential evolution algorithm to optimize both the number of vehicles and the total driving distance at the same time, adaptively adjusts multiple evolutionary strategies, defines the optimization solution with a zero constraint violation degree as a feasible solution, and uses the feasible solution to guide the differential mutation operation to improve the convergence of the optimal solution and quickly obtain the vehicle path optimization plan;

[0091] (2) In view of the problem of poor diversity of the optimal solution in the vehicle path optimization process, based on the constraint domination relationship, the present invention quantitatively describes the constraint violation degree, designs a constraint violation domination criterion, selects good-quality infeasible solutions to enter the optimal solution set, makes full use of the useful information in the infeasible solutions, improves the diversity of the optimal solution, and obtains a vehicle path optimization plan with good diversity.

[0092] It should be noted in particular that the present invention uses the number of vehicles and the total driving distance as optimization objectives, uses a constrained multi-objective differential evolution algorithm, designs a differential mutation strategy based on feasible solution guidance, and a selection strategy based on constraint dominance criteria to perform multi-constraint and multi-objective vehicle path optimization. As long as the optimization objectives and optimization algorithms of the present invention are used to conduct research on vehicle path optimization, it should fall within the scope of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0093] Figure 1 This is a schematic diagram of the differential mutation strategy guided by the feasible solution of the present invention.

[0094] Figure 2 It is the vehicle dispatch route map in the optimization result of the present invention DETAILED DESCRIPTION

[0095] The data sample used in the present invention is C101 belonging to the Solomon data set, as shown in Table 1, the first to seventh columns represent the customer number i, the horizontal coordinate x of the customer location, and the i and the ordinate y i , the customer's demand for goods q i , Earliest service start time i , the latest service start time i and service time i , where the data corresponding to i=0 is the logistics center information.

[0096] The present invention selects the number of vehicles and the total driving distance of vehicles as optimization targets; and adopts the following technical solutions and implementation steps:

[0097] 1 A vehicle path optimization method based on adaptive multi-constraint differential evolution algorithm The specific steps are as follows:

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

[0099] Assume that there are 36 transport vehicles of the same model parked in a logistics center, ready to provide one-way delivery services to 100 customers. The vehicle load limit is 200. The opening and closing times of the logistics center are 0 and 1236 respectively. Each customer corresponds to an order data, including customer number i=0,…,N, customer location (x i ,y i ), the cargo weight q of the i-th customer i , the earliest service start time for the i-th customer l i , the latest service start time e for the i-th customer i and the service time s of the i-th customer i , where i = 0 represents the logistics center;

[0100] The objective function of the vehicle routing optimization problem is:

[0101] f 1 =K (1)

[0102]

[0103] Among them, f 1 The total number of vehicles used in the distribution process K, f 2 is the total distance traveled by the vehicle; 0 <K<K max , d ij is the distance between customer i and customer j, x ijk Indicates whether the kth vehicle travels from i to j, x ijk =1 for yes, x ijk =0 means no, j=0,…,N, k=1,…,K;

[0104] Each vehicle must meet the following conditions during the delivery process, and each customer will only be delivered by one vehicle:

[0105]

[0106]

[0107] Each vehicle departs from the logistics center and returns to the logistics center after completing the delivery task according to the corresponding customer service order:

[0108]

[0109]

[0110] The total customer demand for the kth car should not exceed 200:

[0111]

[0112] The time a that the vehicle arrives at the location of customer i i It should be within its service time window, and the vehicle should return to the logistics center before the logistics center closes:

[0113] l i ≤a i ≤e i (8)

[0114] 0≤a 0 ≤1236 (9)

[0115] The time a for a vehicle to pass through i and reach j j for:

[0116]

[0117] Among them, w i ≥0 indicates the time the vehicle waits for customer i, t 0j represents the time it takes to travel from the logistics center to j, t ij represents the time it takes to travel from i to j, t i0 represents the time it takes to travel from i to the logistics center;

[0118] Step 2: Get the distance and time data for the vehicle routing optimization problem:

[0119] Calculate the distance between customer i and customer j:

[0120]

[0121] Among them, x i and i Represents the horizontal and vertical coordinates of customer i’s location, x j and j Represent the horizontal and vertical coordinates of the customer j's position respectively; according to the average vehicle speed v, calculate the corresponding to d ij The travel time t ij :

[0122]

[0123] Step 3: Design an adaptive multi-constrained differential evolution algorithm:

[0124] Step 3.1 Parameter setting and population initialization:

[0125] The maximum evolutionary number of the multi-objective differential evolution algorithm is set to 1000, the mutation rate is 0.5, the crossover rate is 1, the population size is 50, the neighborhood size is 14, and the initial population P is randomly generated. 0 :

[0126]

[0127] Among them, p b,0 is the bth solution of the initial population, b=1,…,50;

[0128] Step 3.2 Determination of the number of feasible solutions:

[0129] The tth generation population P t It is expressed as:

[0130]

[0131] p b,t =[p b1,t ,...,p bu,t ,…,p b136,t ] (15)

[0132] Among them, p b,t is the bth solution of the tth generation population, p bu,t is the u-th dimension component of the b-th solution of the t-th generation population; t=0,…,1000, p bu,t =1,…,136,u=1,…,136;p bu,t =1,…,100 means customer, p bu,t =101,…,136 represents a logistics center; between two logistics centers, the customer order is sorted in ascending order according to the time when the customer starts the service. The time when the customer starts the service is read from the database, and K delivery routes r are obtained. k , k is the route number, k=1,…,K,K≤36;

[0133] Compute the constraint violation δ for each solution:

[0134]

[0135] Among them, q i is the cargo weight of the i-th customer, K is the total number of routes, i.e. the total number of vehicles used, Q is the vehicle load limit, and a i is the arrival time of the ith customer, which is calculated by formula (10), e i is the latest service start time for the i-th customer, which is read from the database; set the solution satisfying δ=0 as a feasible solution, and the solution satisfying δ>0 as an infeasible solution, and count P t The total number of feasible solutions in ;

[0136] Set P t The set of all feasible solutions in is F t , β represents the total number of feasible solutions, which is an integer in the range of [0,50]. If β=0, execute steps 3.3-3.4 and then step 3.6. If β≥2, execute steps 3.4-3.5 and then step 3.6. If β=1, directly execute step 3.6.

[0137] Step 3.3 Differential evolution based on DE / rand / 1 mutation strategy:

[0138]

[0139]

[0140] p b,t+1 =[p b1,t+1 ,...,p bu,t+1 ,…,p b136,t+1 ] (19)

[0141] Among them, v b,tis the b-th mutation solution of the t-th generation population, p rand1,t is the rand of the t-th generation population 1 solution, p rand2,t is the rand of the t-th generation population 2 solution, rand 1 and rand 2 are two different integers randomly selected from [1,50] that are different from b; p bu,t+1 is the u-th dimension component of the b-th solution in the t+1-th generation population, rand bu is the random number taken by the b-th solution on the i-th dimension component, 0≤rand bu ≤1;

[0142] P b,t+1 Each dimension is rounded down:

[0143]

[0144] in, Indicates p b,t+1 The value of the u-th dimension component rounded down;

[0145] If 101≤p bu,t+1 ≤136,101≤p b(u+1),t+1 ≤136, and 1≤p b(u+2),t+1 ≤100, that is, p b,t+1 When there is an empty route, p b(u+1),t+1 Remove the route group, otherwise set u = u + 1 and continue to check empty routes until p b,t+1 Delete all empty routes in the .

[0146] Step 3.4 Fitness value calculation:

[0147] Calculate the dominance strength of each solution in the candidate solution set:

[0148] If p is a feasible solution and q is an infeasible solution, then p dominates q, and the dominance strength of p increases by 1;

[0149] If p and q are both feasible solutions, f 1 (p)≤f 1 (q), f 2 (p)≤f 2 (q), then p dominates q, and the dominance strength of p increases by 1;

[0150] If both p and q are infeasible solutions, δ(p)<δ(q), then p dominates q, and the dominance strength of p increases by 1;

[0151] Among them, p and q are two different solutions in the candidate solution set, f 1 (p) and f 1(q) are the total number of vehicles corresponding to p and the total number of vehicles corresponding to q, respectively, and f 2 (p) and f 2 (q) are the total vehicle travel distance corresponding to p and the total vehicle travel distance corresponding to q, respectively. δ(p) and δ(q) are the constraint violation degree corresponding to p and the constraint violation degree corresponding to q, respectively.

[0152] Calculate P t and P t+1 The rank of each solution in the candidate solution set is:

[0153] For the solution with δ = 0, i.e., a feasible solution, its rank is set to 1;

[0154] For solutions with δ>0, i.e. infeasible solutions, they are sorted in ascending order according to the size of δ, and the rank of each infeasible solution is set to the sum of the total number of feasible solutions and the corresponding δ ranking;

[0155] Add the dominance strength and the level value of each solution to obtain the fitness value of each solution;

[0156] Step 3.5: Differential mutation operation guided by feasible solution:

[0157] ① According to the fitness value, F t Sort in ascending order to get F t =(p 1 ,...,p β );

[0158] ② Calculate p 1 and p 2 Similarity:

[0159]

[0160] Among them, r p1,rand and r p2,rand Respectively represent from p 1 A delivery route randomly selected from p 2 A delivery route randomly selected from max ,|r p1,rand | and |r p2,rand |represents r p1,rand and r p2,rand The number of customers, p1,rand ∩r p2,rand | indicates that at the same time p1,rand and r p2,rand The number of customers that appear in |r p1,rand ∩r p2,rand |≥0;

[0161] ③ Adjust F t The order:

[0162] p 2,t Adjust to p β , and (p 3 ,...,p β ) is adjusted to (p 2 ,...,p β-1 );

[0163] ④If S≤0.5, execute step 3.5⑤; if S>0.5, return to step 3.5②;

[0164] ⑤Select p in the order of route number 1 A route p1,k1 , k 1 Yes 1 The corresponding route number, 1≤k 1 ≤36, set c(1,r p1,k1 ) indicates p 1 k 1 The first customer of the delivery route, in p 2 Find c(1,r p1,1 )Route p2,k2 , k 2 Yes 2 The corresponding route number, 1≤k 2 ≤36, determine r p1,k1 and r p2,k2 The difference vector pair of :

[0165]

[0166] Among them, r p1,k1 -r p1,k1 ∩r p2,k2 Indicates that p1,k1 Rather than p2,k2 The complement of the delivery routes on p2,k2 -r p1,k1 ∩r p2,k2 Indicates that p2,k1 Rather than p1,k2 The complement of the delivery routes on the p1,k1 -r p1,k1 ∩r p2,k2 and r p2,k2 -r p1,k1 ∩r p2,k2 first customer;

[0167] ⑥Let p b,t+1 =p 1 , in p 1 Find the route r where n is located p1,k3 , k 3 Yes1 Different from k 1 The route number, i.e. 1≤k 3 ≤36, k 3 ≠k 1 ; Change m and n from r p1,k1 and r p1,k3 Removed;

[0168] Under the premise of satisfying formula (3)-formula (10), find r p1,k3 All positions of m can be inserted so that δ(p b,t+1 )=0 holds, calculate the route increments corresponding to all insertion positions, set the position corresponding to the smallest route increment as the insertion position of m, and insert m into this position;

[0169] Under the premise of satisfying formula (3)-formula (10), find r p1,k3 All positions of n can be inserted so that δ(p b,t+1 )=0 holds, calculate the route increments corresponding to all insertion positions, set the position corresponding to the smallest route increment as the insertion position of n, and insert n into this position;

[0170] If there is no such b,t+1 )=0, set m or n to r p1,k3 or p1,k1 The first customer to complete p b,t+1 The construction of

[0171] ⑦Let k 1 =k 1 +1, count the population size of the t+1th generation |P t+1 |, if k 1 =K 1 , K 1 Yes 1 The corresponding number of vehicles, and |P t+1 |=50, then go to step 3.6. If k 1 =K 1 , and |P t+1 |<50, return to step 3.5②; if k 1 <K 1 , and |P t+1 |<50, return to step 3.5⑤;

[0172] Step 3.6 Large-scale neighborhood search operation:

[0173] In p b,w Randomly select a one-dimensional component p from bu,w =i, w=t or t+1, 1 <p bu,w <100, based on the distance d between customer i and customer jij Select the 14 customers closest to customer i, and find the b,w corresponding 14-dimensional component p bv,w in p, where 1 < p bv,w < N, v ≠ u, and remove p bv,w and the selected 14-dimensional component from p b,t ;

[0174] On the premise of satisfying formulas (3) - (10), find all the positions in p b,w where p bu,w can be inserted, so that δ(p b,w ) = 0 holds. Calculate the route increment of all insertion positions, and set the position corresponding to the smallest route increment as the insertion position of p bu,w , and insert p bu,w into this position;

[0175] On the premise of satisfying formulas (3) - (10), find all the positions in p b,w where p bv,w can be inserted, so that δ(p b,w ) = 0 holds. Calculate the route increment of all insertion positions, and set the position corresponding to the smallest route increment as the insertion position of p bv,w , and insert p bv,w into this position;

[0176] Step 3.7 Selection operation:

[0177] If P t+1 is empty, that is, there is no (t + 1)-th generation population, go to Step 3.8. Otherwise, sort the candidate solution set in ascending order according to the fitness value of the solutions, and select the first 50 solutions to construct P t+1 ;

[0178] Step 3.8 Judge the termination condition of the algorithm:

[0179] If t < 1000, then let t = t + 1, and repeat Steps 3.2 to 3.7. If t = T, then stop the calculation, output the vehicle dispatching plan with the minimum number of vehicles and the vehicle driving distance, and the corresponding path planning diagram, as shown in Table 2 and Figure 2 shown. Table 2 shows that the logistics center dispatched a total of 10 vehicles to serve 100 customers, and the total driving distance corresponding to the optimal vehicle dispatching plan is 828.94, Figure 2 and the corresponding vehicle path planning diagram is given.

[0180] Table 1. C101 vehicle path optimization data

[0181]

[0182]

[0183]

[0184] Table 2. Optimal vehicle dispatching solution

[0185]

Claims

1. A vehicle path optimization method based on an adaptive multi-constrained differential evolution algorithm, characterized in that: It includes the following steps: Step 1: Determine the optimization objective of the vehicle routing optimization problem: Assume that K max All transport vehicles of the same model are parked in a logistics center, ready to provide one-way delivery services to N customers. The vehicle load limit is Q, the logistics center opening time is L and the logistics center closing time is E; each customer corresponds to an order data, including customer number i, customer location coordinates (x i ,y i ), the cargo weight q of the i-th customer i , the earliest service start time for the i-th customer l i , the latest service start time e for the i-th customer i and the service time s of the i-th customer i , where i = 0 represents the logistics center; 0 <K max <50,0 <N<200,0<Q<300,0<L<200,0<E<200,i=0,…,N,0<x i <100,0 <y i <100,0 i <300,0 <l i <1000,0 <e i <1440,0 i <200;​​ The optimization objective of the vehicle routing optimization problem is to minimize the following objective function: f1 = K (1) Among them, f1 is the total number of vehicles K used in the distribution process, f2 is the total driving distance of the vehicles; <K<K max , d ij is the distance between customer i and customer j, x ijk Indicates whether the kth vehicle travels from i to j, x ijk =1 for yes, x ijk =0 means no, j=0,…,N, k=1,…,K; Each vehicle needs to meet the following conditions during the distribution process, and each customer is served by only one vehicle: Each vehicle departs from the logistics center, completes the distribution task according to the corresponding customer service order, and then returns to the logistics center: The total demand of customers of the k-th vehicle should not exceed Q: The time a that the vehicle arrives at the location of customer i i It should be within its service time window, and the vehicle should return to the logistics center before the logistics center closes: l i ≤a i ≤e i (8) L ≤ a0 ≤ E (9) The time a for vehicle to reach location i after passing j i for: Among them, w j is the time the vehicle waits for customer j, w j ≥0,t 0i is the time taken to travel from the logistics center to customer i, t ji represents the time it takes to travel from customer j to customer i, t j0 represents the time it takes to travel from customer j to the logistics center; Step 2: Obtain the distance and time data of the vehicle routing optimization problem: Calculate the distance between customer i and customer j: Among them, x i and i Respectively represent the horizontal and vertical coordinates of customer i’s location, x j and j Represent the horizontal and vertical coordinates of the customer j's position respectively; according to the average vehicle speed v, calculate the corresponding to d ij The travel time t ij : Step 3: Design an adaptive multi-constraint differential evolution algorithm: Step 3.1: Parameter setting and population initialization: Set the maximum number of generations of the multi-objective differential evolution algorithm to T, where T is an integer in the range of [250, 2000], the mutation rate is F, 0 ≤ F ≤ 1, the crossover rate is Cr, 0 ≤ Cr ≤ 1, the population size is H, H is an integer in the range of [10, 500], the neighborhood size is c, c is an integer in the range of [1, N], each solution represents a vehicle scheduling plan, and randomly generate the initial population P0: Among them, p b,0 is the b-th solution of the initial population, b=1,…,H; Step 3.2: Judgment of the number of feasible solutions: The tth generation population P t It is expressed as: Among them, p b,t is the bth solution of the tth generation population, p bu,t is the u-th dimension component of the b-th solution of the t-th generation population; t=0,…,T,p bu,t =1,…,N+K max ,u=1,…,N+K max ;p bu,t =1,…,N indicates the customer, p bu,t =N+1,…,N+K max represents a logistics center; between two logistics centers, the customer order is sorted in ascending order according to the time when the customer starts the service. The time when the customer starts the service is read from the database, and K delivery routes r are obtained. k , k is the route number, k=1,…,K,K≤K max ; Calculate the constraint violation degree δ of each solution: Among them, q i is the cargo weight of the i-th customer, K is the total number of routes, i.e. the total number of vehicles used, Q is the vehicle load limit, and a i is the arrival time of the ith customer, which is calculated by formula (10), e i is the latest service start time for the i-th customer, which is read from the database; set the solution satisfying δ=0 as a feasible solution, and the solution satisfying δ>0 as an infeasible solution, and count P t The total number of feasible solutions in ; Set P t The set of all feasible solutions in is F t , β represents the total number of feasible solutions, which is an integer in the range [0,H]. If β=0, execute steps 3.3-3.4 and then step 3.

6. If β≥2, execute steps 3.4-3.5 and then step 3.

6. If β=1, directly execute step 3.

5. Step 3.3: Differential evolution based on the DE / rand / 1 mutation strategy: Among them, v b,t is the b-th mutation solution of the t-th generation population, p rand1,t is the rand1th solution of the tth generation population, p rand2,t is the rand2th solution of the tth generation population, rand1 and rand2 are two different integers randomly selected from [1,H] that are different from b; p bu,t+1 is the u-th dimension component of the b-th solution in the t+1-th generation population, rand bu is the random number taken by the b-th solution on the i-th dimension component, 0≤rand bu ≤1; P b,t+1 Each dimension is rounded down: in, Indicates p b,t+1 The value of the u-th dimension component rounded down; If N+1≤p bu,t+1 ≤N+K max ,N+1≤p b(u+1),t+1 ≤N+K max , and satisfy 1≤p b(u+2),t+1 ≤N, that is, p b,t+1 When there is an empty route, p b(u+1),t+1 Remove the route group, otherwise set u = u + 1 and continue to check empty routes until p b,t+1 Delete all empty routes in the . Step 3.4: Fitness value calculation: Calculate the domination strength of each solution in the candidate solution set: If p is a feasible solution and q is an infeasible solution, then p dominates q and the domination strength of p is incremented by 1; If both p and q are feasible solutions, f1(p) ≤ f1(q), f2(p) ≤ f2(q), then p dominates q and the domination strength of p is incremented by 1; If both p and q are infeasible solutions and δ(p) < δ(q), then p dominates q and the domination strength of p is incremented by 1; where p and q are two different solutions in the candidate solution set, f1(p) and f1(q) are the total number of vehicles corresponding to p and q respectively, f2(p) and f2(q) are the total driving distances of the vehicles corresponding to p and q respectively, and δ(p) and δ(q) are the constraint violation degrees corresponding to p and q respectively; Calculate P t and P t+1 The rank value of each solution in the candidate solution set is: For a solution with δ = 0, that is, a feasible solution, its rank is set to 1; For a solution with δ > 0, that is, an infeasible solution, sort it in ascending order according to δ, and the rank of each infeasible solution is the sum of the total number of feasible solutions and its δ sorting value; Add the domination strength of each solution to the rank value to obtain the fitness value of each solution; Step 3.5: Differential mutation operation guided by feasible solutions: ① According to the fitness value, F t Sort in ascending order and get F t =(p1,...,p β ); ② Calculate the similarity S between p1 and p2: Among them, r p1,rand and r p2,rand They represent a delivery route randomly selected from p1 and a delivery route randomly selected from p2, 1≤rand≤K max ,|r p1,rand | and |r p2,rand |represents r p1,rand and r p2,rand The number of customers, p1,rand ∩r p2,rand | indicates that at the same time p1,rand and r p2,rand The number of customers that appear in |r p1,rand ∩r p2,rand |≥0; ③ Adjust F t The order: p2 is adjusted to p β , and (p3,...,p β ) is adjusted to (p2,...,p β-1 ); ④ If S ≤ 0.5, execute Step 3.5⑤, if S > 0.5, return to Step 3.5②; ⑤ Select a route r from p1 in the order of route numbers p1,k1 , k1 is the route number corresponding to p1, 1≤k1≤K max , set c(1,r p1,k1 ) represents the first customer of the k1th delivery route in p1, and finds c(1,r p1,1 )Route p2,k2 , k2 is the route number corresponding to p2, 1≤k2≤K max , determine r p1,k1 and r p2,k2 The difference vector pair of : Among them, r p1,k1 -r p1,k1 ∩r p2,k2 Indicates that p1,k1 Rather than p2,k2 The complement of the delivery routes on p2,k2 -r p1,k1 ∩r p2,k2 Indicates that p2,k1 Rather than p1,k2 The complement of the delivery routes on the p1,k1 -r p1,k1 ∩r p2,k2 and r p2,k2 -r p1,k1 ∩r p2,k2 first customer; ⑥Let the bth solution p of the t+1th generation population b,t+1 =p1, find the route r where n is located in p1 p1,k3 , k3 is the route number in p1 that is different from k1, 1≤k3≤K max , k3≠k1, change m and n from r p1,k1 and r p1,k3 Removed: Under the premise of satisfying formula (3)-formula (10), find r p1,k3 All positions of m can be inserted so that δ(p b,t+1 )=0 holds, calculate the route increments corresponding to all insertion positions, set the position corresponding to the smallest route increment as the insertion position of m, and insert m into this position; Under the premise of satisfying formula (3)-formula (10), find r p1,k3 All positions of n can be inserted so that δ(p b,t+1 )=0 holds, calculate the route increments corresponding to all insertion positions, set the position corresponding to the smallest route increment as the insertion position of n, and insert n into this position; If there is no such b,t+1 )=0, set m or n to r p1,k3 or p1,k1 The first customer to complete p b,t+1 The construction of ⑦ Let \(k1 = k1 + 1\), and count the population size \(|P|\) of the \((t + 1)\)-th generation. t+1 If \(k1 = K1\), where \(K1\) is the number of vehicles corresponding to \(p1\), and \(|P|\) t+1 \(= H\), where \(H\) is the population size, then go to step 3.

6. If \(k1 = K1\) and \(|P|\) t+1 \(< H\), then return to step 3.5②; if \(k1 < K1\) and \(|P|\) t+1 \(< H\), then return to step 3.5⑤. Step 3.6: Large-scale neighborhood search operation: Among p b,w randomly select a one-dimensional component p representing a customer bu,w , p bu,w = i, w = t or t + 1, 1 < p bu,w < N, where N is the total number of customers served. According to the distance d between customer i and customer j ij select c customers closest to customer i, where c is the neighborhood size, and find the corresponding c-dimensional components p b,w in bv,w , 1 < p bv,w < N, v ≠ u, and remove p bv,w and the selected c-dimensional components from p b,t ; Under the premise of satisfying formula (3)-formula (10), find p b,w You can insert p bu,w All positions of δ(p b,w )=0, calculate the route increments of all insertion positions, and set the position corresponding to the minimum route increment as p bu,w Insertion position, and p bu,w Insert at this location; Under the premise of satisfying formula (3)-formula (10), find p b,w You can insert p bv,w All positions of δ(p b,w )=0, calculate the route increments of all insertion positions, and set the position corresponding to the minimum route increment as p bv,w Insertion position, and p bv,w Insert at this location; Step 3.7: Selection operation: If P t+1 is empty, that is, there is no t+1 generation population, go to step 3.8, otherwise sort the candidate solution set in ascending order according to the fitness value of the solution, select the first H solutions to construct P t+1 , H is the population size; Step 3.8: Judge the termination condition of the algorithm: If t < T, let t = t + 1, and repeat Steps 3.2 to 3.

7. If t = T, stop the calculation and output the optimal vehicle scheduling plan.

Images