Intelligent pickup and delivery collocation for logistics models

2.1 Briefly describe the process

In MVSPDPTW, there is a given distribution center, a group of customers with different geographical locations and a group of vehicles with multiple models. The goal of MVSPDPTW is to select the right vehicle to perform all customer tasks, so that the total length of the vehicle’s travel path is minimized.

0 means the distribution center, T ={ 1, 2, …, n } means the collection of customer points, each customer point corresponds to a collection/distribution task, the customer point corresponding to the collection task is called the pickup point, the customer point corresponding to the delivery task is called the delivery point. For customer point i ∈ T, [e _i , l _i ] represents the i-th time window, e _i represents the earliest service start time, and l _i represents the latest service start time. q _i represents the quantity of goods i need to collect/deliver. Each customer point can be visited by multiple vehicles, but it is only allowed to be visited by the same vehicle at most once. K =  { 1, 2, …, m } represents the set of vehicles. For k ∈ K,  s _k , n _k ,  Q _k ,  o _k , and f _k represent the starting point, ending point, maximum capacity, earliest departure time, and latest return time of vehicle k, respectively. As shown in Fig. 1, the vehicle needs to depart from the starting point after the earliest departure time, collect the goods at the pickup point and send them back to the distribution center, load the goods at the distribution center and deliver them to the delivery point, and finally return to the destination before the latest return time. The vehicle is allowed to visit the distribution center multiple times during the mission.

OPEN IN VIEWER

Fig. 1

Schematic diagram of vehicle execution tasks.

Since the vehicle may visit the distribution center many times, if a variable is used in the model to distinguish how many times the vehicle visits the distribution center, the difficulty of modeling and solving will increase significantly. In order to reduce the difficulty, the following processing methods are adopted: T ={ 1, 2, …, n } represents the collection of customer points, n is the actual number of customer points, for the pickup point numbered i, a numbered n + i virtual delivery point is constructed to correspond to it; for the delivery point numbered j, it is renumbered as n + j, and a virtual pickup point numbered j is constructed to correspond to it. After the above processing, T =  T _P   ∪  T _D , T _P  ={ 1, 2, …, n } represents the set of pickup points, T _D  ={ n + 1,  n + 2, …, 2n } represents the set of delivery points, and the vehicle needs to collect the goods from point i (i ∈ T _P ) and deliver them to point n + j  (n + j  ∈ T _D ). In this processing method, the distribution center is replaced with n virtual customer points, all virtual customer points are located at the same location as the distribution center and there is no service time window restriction; the amount of goods collected/delivered by each virtual customer point corresponds to the equal actual customer in the points. As shown in Fig. 2, the distribution center in Fig. 1 is virtualized and replaced with virtual customer points in the box, where i^* = n + j.

OPEN IN VIEWER

Fig. 2

Virtualization schematic.

After the above processing, MVSPDPTW is defined on the upper graph of the fully directed graph G = (V,  A), V ={ S ∪ N  ∪  T ∪  0 } represents the vertex set, A ={ (i, j) : i,  j ∈ V,  i ≠ j } represents the arc set.

S: The starting point set of the vehicle, S ={ s₁, s₂, …, s _m  };

N: the set of end points of the vehicle, N ={ n₁, n₂, …, n _m  };

T: customer point collection, T ={ 1, 2, …, 2n };

d _ij : the distance between point i and point j, i, j ∈ V;

t _ij : the time between point i and point j, i, j ∈ V;

b _ik : the time when the vehicle k arrives at the point i, k ∈ K, i ∈ T ∪ 0;

a _ik : the start time of vehicle k at point i, k ∈ K, i ∈ T ∪ 0;

r _ik : the service time of vehicle k at point i, k ∈ K, i ∈ T ∪ 0;

q _ik : the amount of cargo loaded/unloaded by vehicle k at point i, k ∈ K, i ∈ T ∪ 0;

w _ik : Cargo capacity of vehicle k after point i, k ∈ K, i ∈ T ∪ 0;

x _ijk : Indicates whether the vehicle k is directly from point i to point j, if yes, it is 1, otherwise it is 0;

M: Represents a very large integer.

2.4 Linear model

min ∑ i ∈ V ∑ j ∈ V ∑ k ∈ K x ijk d ij(1) ∑ i ∈ T P x s k ik - ∑ j ∈ T D x jn k k = 0 , ∀ k ∈ K(2) ∑ i ∈ S ∪ T x ijk - ∑ i ∈ T ∪ N x jik = 0 , ∀ k ∈ K , j ∈ T(3) ∑ k ∈ K ∑ i ∈ S ∪ T x ijk ⩾ 1 , ∀ j ∈ T(4) ∑ i ∈ S ∪ T x ijk ⩽ 1 , ∀ k ∈ K , j ∈ T(5) ∑ j ∈ T x ijk - ∑ j ∈ T ∪ N x ( n + i ) jk = 0 , ∀ k ∈ K , i ∈ T P(6) a ik ⩽ a ( n + i ) k , ∀ k ∈ K , i ∈ T P(7) q ik = q ( n + i ) k , ∀ k ∈ K , i ∈ T P(8) w s k k = w n k k = 0 , ∀ k ∈ K(9) 0 ⩽ w ik ⩽ Q k , ∀ k ∈ K , i ∈ T(10)

w ik + q jk - ( 1 - x ijk ) M ⩽ w jk , ∀ k ∈ K , j ∈ T P , i ∈ S ∪ T(11)

w jk + q jk - ( 1 - x ijk ) M ⩽ w ik , ∀ k ∈ K , j ∈ T D , i ∈ T(12) ∑ k ∈ K q ik = q i , i ∈ T(13) a ik + t s k i - ( 1 - x s k ik ) M ⩾ o k , ∀ k ∈ K , i ∈ T P(14) a ik + r ik + t in k - ( 1 - x in k k ) ⩽ f k , ∀ k ∈ K , i ∈ T D(15)

a ik + r ik + t ij - ( 1 - x ijk ) M ⩽ b jk , ∀ k ∈ K , i ∈ T , j ∈ T(16) b ik ⩽ a ik , ∀ k ∈ K , i ∈ T(17) e i ⩽ a ik ⩽ l i , ∀ k ∈ K , i ∈ T(18) x ijk ∈ { 0 , 1 } , ∀ k ∈ K , i ∈ V , j ∈ V(19)

Wherein, formula (1) is the objective function, which minimizes the total length of the mission vehicle’s travel path; formulas (2) and (3) indicate that the vehicle will perform the task from the starting point and finally return to the end; formulas (4) and (5) indicate every customer point can be accessed by multiple vehicles, but it is only allowed to be accessed by the same vehicle at most once; formulas (6), (7), (8) indicate that the vehicle needs to collect goods from the pickup point for delivery to the corresponding delivery point; formula (9) means that the cargo carrying capacity of the vehicle from the starting point and returning to the destination is 0; formula (10) expresses the capacity constraint; formulas (11) and (12) express the change of the vehicle carrying cargo capacity at the customer point; formula (13) means the sum of the divided cargo volume of each customer point is equal to the total cargo volume of the customer point; formulas (14) and (15) indicate that the vehicle needs to depart from the starting point after the earliest departure time to provide services to customers, and return before the latest return time end point; formulas (16) and (17) express the relationship between the time when the vehicle arrives at the customer point and the time when it starts service at the customer point; formula (18) expresses that the vehicle can only provide service within the customer’s time window.

3.1 Initial solution generation

The simulated annealing algorithm is essentially an iterative solution algorithm. The initial solution determines the starting point of the iteration. Therefore, a good initial solution helps the algorithm find a better final solution. In the initial solution generation method in this article, a simple and effective dynamic greedy algorithm is used to generate the initial solution [5]. k ^ = arg k ∈ K m a x ( f k − o k )(20) j ^ = arg j ∈ T m i n { t i j + max ( e j − b j k , 0 ) }(21)

Formulas (20) and (21) respectively represent the selection criteria of vehicles and customer points, and the selection of vehicles is based on the principle of maximizing service time. The selection of customer points is based on the principle of minimizing the total travel time and waiting time between the last visit object i in the current planning path and the customer point j to be selected. Formula (21) is used to select service customer points for the current vehicle in turn, and to check whether the selected customer points meet the insertion constraints, including time window constraints, capacity constraints, customer point visit times constraints, and demand constraints. Demand constraints refer to all delivery of the goods which come from the distribution center, and all collected goods are sent back to the distribution center. If the constraints are met, the customer points are inserted into the current vehicle planning path, and this process is repeated until all customer needs are met.

The specific process of the dynamic greedy algorithm is as follows:

Input: vehicle collection K, customer point collection T, distribution center 0.

Output: Initial feasible solution S _initial .

Step 1: If the set T is not an empty set, go to step 2, otherwise, the algorithm ends to output the result.

Step 2: Formula (20) is used to select the current vehicle k from K.

Step 3: Construct a set of customer points to be selected C = T.

Step 4: If C is not an empty set, go to step 5, otherwise, go to step 9.

Step 5: Formula (21) is used to select customer j from C.

Step 6: Check whether j satisfies the insertion constraint, if so, insert j into the current vehicle planning path R _k .

Step 7: Check whether all the requirements of j are satisfied, and delete j from T if they are satisfied.

Step 8: Delete j from C, and go to step 4.

Step 9: Check whether the last visited object of R _k is distribution center 0, if yes, go to step 11, otherwise, go to step 10.

Step 10: Check whether 0 meets the time window constraint, if it is satisfied, insert 0 into R _k , go to step 3, otherwise, go to step 11.

Step 11: Delete k from K, go to step 1

In the iterative process, accepting infeasible solutions can help increase the search range of the neighborhood and improve the ability of the algorithm to jump out of the local optimum, and it helps the algorithm to obtain a better feasible solution in the iterative step [19], so a violation constraint punishment mechanism is designed in this paper. The overall constraint penalty function is designed in formula (22), B represents the fitness value of the path planning solution, β represents the parameter penalty weight, η (i, k) , φ (i, k) , λ (i, k) represent the capacity constraint evaluation value, time window constraint evaluation value and demand constraint evaluation value of point i in R _k , respectively. Their calculation methods are given by the formula (23), (24), (25). R _k  ={ o _k , …, n _k  } represents the execution task path of vehicle k, which is arranged according to the order of vehicle visits.

B = ∑ i ∈ V ∑ i ∈ V ∑ k ∈ K x ijk d ij + β ∑ k ∈ K ∑ i ∈ R k ( η ( i , k ) + φ ( i , k ) + λ ( i , k ) )(22) η ( i , k ) = max ( w ik - Q k , 0 ) , i ∈ V(23)

φ ( i , k ) = { q ik , a ik - l i > 0 , i ∈ T Q k , a ik + r ik + t in k > f k , i is the penultimate access object in R k 0 , other(24)

λ ( i , k ) = { q ik , i is the pickup point , and there is no 0 behind i in R k q ik , i is the delivery point , and there is no 0 in front of i in R k 0 , other(25)

Tabu search is a meta-heuristic algorithm that simulates the human thought process. Its basic idea is to avoid repeated searches in the searched space and force the algorithm to explore other new solutions.

The current solution of each iteration is taken as the taboo object and stored in the tabu list [20]. In the subsequent iterations of l (the length of the tabu list), the object is set to the taboo state, which is different from the side movement taboo in taboo structures [21, 22]. The length of the tabu list refers to the number of objects which are stored in the tabu list. When the number of the stored objects reaches the upper limit and another object is added, the first in, first out principle is adopted to remove the previous objects from the tabu list and then add them.

3.4 Multi neighborhood structure

The multi-neighbor structure can increase the diversity of solutions and improve the possibility of the algorithm to find the best global solution. In this paper, five kinds of neighborhood search operators are used, the migration operator, the migration division operator, and the exchange operator are derived from the classic neighborhood search operator [23]. In the iterative process of this article, the repair operator and the car exchange operator are designed to accept the infeasible solution and the characteristics of multiple models, and they are used to repair the violation of capacity constraints and the operation of changing vehicles respectively.

In order to better explain the operator operation, the following symbols are defined. R _k  (k ∈ K) and R _h  (h ∈ K) represent the execution task path of vehicle k and vehicle h, respectively. For client points i ∈ R k , z ik 0 and z ik l respectively represents the amount of goods which are loaded/unloaded by vehicle k on customer i before and after the operator operation. The calculation method of B (R _k ) is given in Equation (26), and B (R _k ) represents the fitness value of path R _k .

B ( R k ) = ∑ i ∈ V ∑ i ∈ V x ijk d ij + β ∑ i ∈ R k ( η ( i , k ) + φ ( i , k ) + λ ( i , k ) ) , ∀ k ∈ K(26)

Task set: In the planned vehicle path, the starting point, ending point of the vehicle and non-customer points such as the distribution center are used as the boundary, consecutive customer points between two adjacent non-customer points form a task set. If there is no client point between two adjacent non-client points, it does not constitute a task set. As shown in Fig. 3, customer points 1 and 2 form a task set. If there is no customer point between the destination and the distribution center, it does not constitute a task set.

OPEN IN VIEWER

Fig. 3

Schematic diagram of the task set.

Migration operator: The migration operator is used to delete the customer point from the original path and insert it into other paths. As shown in Fig. 4. For the customer point i ∈  R _k , i is deleted from R _k and i is inserted into R _h with the principle to minimize B (R _h ).

OPEN IN VIEWER

Fig. 4

Schematic diagram of migration operator operation.

Migration segmentation operator: As shown in Fig. 5, customer i is served by vehicles k and h, customer j is served by vehicle h, and i and j belong to the same task set in the path R _h . Delete i from the path R _k , and then insert j into R _k according to the principle of minimizing B (R _k ). At this time, i is served by vehicle h, j is served by vehicles k and h together. When vehicles k and h are loaded at i and j. the amount of unloaded cargo is adjusted to z ih l = z ik 0 + z ih 0 , z ik l = 0 , z jk l = z ik , 0 z jh l = z jh 0 - z jk l .

OPEN IN VIEWER

Fig. 5

Schematic diagram of migration and division operator operation.

Exchange operator: Exchange operator is used to exchange customer points in two paths. As shown in Fig. 6, customer point i ∈ R _k and customer point j ∈ R _h . Exchange i and j. After the exchange, i ∈ R _h , j ∈ R _k .

OPEN IN VIEWER

Fig. 6

Schematic diagram of exchange operator operation.

Repair operator: The repair operator is used to repair violations of vehicle capacity constraints. As shown in Fig. 7, if the cargo capacity of vehicle k exceeds the maximum capacity Q _k of the vehicle when serving customer i, choose a customer point j of the tasks set to which i belongs, j can be i, choose a non-empty path R _h , insert j into R _h with the principle of minimizing B (R _h ), and adjust the amount of goods which are loaded/unloaded by vehicles k and h at j to: z jk l = z jk 0 - η ( i , k ) , z jh l = η ( i , k ) , the calculation method of η(i,k) is given by formula (23).

OPEN IN VIEWER

Fig. 7

Schematic diagram of repair operator operation.

Car exchange operator: As shown in Fig. 8, for non-empty paths R _k and empty paths R _h , all customer points in R _k form a customer point set C, and then according to the path generation method in the initial solution, the customer points in C are inserted into R _h . If all the customer points in C are successfully inserted into R _h , then the path R _k is cleared, otherwise, the customer points which are inserted into R _h will be deleted from R _k .

OPEN IN VIEWER

Fig. 8

Schematic diagram of car change operator operation.

Since a single customer point is only allowed to be accessed by the same vehicle at most once, after the operator operation, if the situation shown in Fig. 9 occurs, R _k contains two customer points i, then one i is selected to be deleted from R _k based on the principle of minimizing B (R _k ), and to combine the quantity of goods which are loaded/unloaded by two services i of k.

OPEN IN VIEWER

Fig. 9

Schematic diagram of merged customer points.

The dynamic greedy algorithm is used to generate the initial feasible solution, and it is set to the current solution S _now of the 0-th generation and the global best solution S _best , and then the iteration starts. In each iteration, the design criteria is used to update S _now and S _best . Finally, when the current temperature is less than the set minimum temperature, the algorithm terminates and outputs S _best . The specific algorithm flow is as follows:

Input: Initial temperature T₀, termination temperature T _end , cooling rate R _T , iteration step size L, number of neighborhood solutions N _u , initial feasible solution S _initial .

Output: S _best is the best solution globally.

Step 1: Set S _initial to S _now and S _best .

Step 2: Initialize the current temperature T′ = T₀ and initialize the number of iterations L′ = 0 at the current temperature.

Step 3: Initialize the tabu list and add S _now to the tabu list.

Step 4: If T> T _end , go to step 5, otherwise, the algorithm ends to output the result.

Step 5: If L = L+ 1 go to step 6, otherwise, , go to step 4.

Step 6: A neighborhood operator is randomly selected from the five neighborhood operators to convert the current solution and to generate a neighborhood solution.

Step 7: Repeat step 6 until the number of neighborhood solutions reaches N _u .

Step 8: Select the best neighborhood solution S _can which is not taboo from N _u neighborhood solutions.

Step 9: If S _can  < S _now , S _can will become the new S _now , go to step 11, otherwise, go to step 10.

Step 10:

If e^{(S_can-S_now)/T} > Random (0, 1), S _can will become the new S _now , Random(0,1) is used to generate random data uniformly distributed on (0,1), otherwise keep S _now .

Step 11: If S _can is a feasible solution and satisfies S _can  < S _best , S _can will become the new S _best , otherwise keep S _best .

Step 12: If S _now changes, add S _now to the taboo list.

Step 13: Go to step 5.

4.1 Parameter setting

Setting reasonable parameters can effectively improve the performance of the algorithm [26]. First, the parameter settings in references [27, 28] are used to determine the test interval of each parameter in this article, and then different scale examples are selected to perform multiple tests by changing the parameter combination. Finally, the best combination of parameters is determined.

The value range and step length of each parameter are shown in Table 1. In order to simplify the calculation complexity, the controlled variable method is used to determine the optimal value of each parameter in turn. The final result of TSA parameter setting is as follows: initial temperature T₀ = 100, termination temperature T _end  = 0.1, iteration length L = 30, cooling rate R _T  = 0.93, penalty weight β= 4, number of neighborhood solutions N _u  = 2n + 150 (n is the number of customers), the length of the tabu list is l = 8.

As there is currently no recognized data set available in the field of MVSPDPTW, there is no way to directly test the algorithm in this article. MVSPDPTW is an extension of the classic VRPTW and SDVRPTW. Many of the constraints in the problem are similar. Therefore, if the algorithm in this article can solve VRPTW and SDVRPTW well, then there is reason to believe that the algorithm in this paper can also effectively solve MVSPDPTW.

12 MVSPSPTW simulation examples were constructed based on 12 SDVRPTW examples. Without changing the distribution center location, customer location, customer cargo volume, and customer time window in the SDVRPTW example, only half of the customers were randomly selected to change their distribution tasks to collection tasks. The starting point and ending point of the vehicle are no longer set as the distribution center, but they are randomly generated within the 100×100 square area that the customer belongs to. The starting service time of a vehicle is a randomly generated integer between 0 and 800, and the end service time is equal to a randomly generated integer between 500 and 900 plus the starting service time of the vehicle; the maximum vehicle capacity is a randomly generated integer between 70 and 120. In addition, when the number of customers is 25, 15 vehicles are randomly generated, and when the number of customers is 50, 30 vehicles are randomly generated to ensure sufficient vehicle resources.

In Table 4, the TSA algorithm, the simulation results’ comparison of the traditional simulated annealing (SA) algorithm, and particle swarm optimization (PSO) algorithm to solve the MVSPDPTW are recorded. The design and processing method of the PSO algorithm comes from the reference [29, 30]. The parameters of the traditional SA algorithm and the PSO algorithm are determined according to the parameter setting method in Table 1. Gap1 and Gap2 respectively represent the gap between TSA algorithm solving each example result (R _best ) and SA algorithm solving each example result (T _best ), PSO algorithm solving each example result (P _best ), Gap1 = (T _best - R _best )/R _best , Gap2 = (P _best  - R _best )/R _best .

According to the data in Table 4, the TSA algorithm in this paper has better solution performance than the traditional SA algorithm and PSO algorithm. The main factors are: (i) The TSA algorithm has the best solution result, which is up to 2.40% higher than the traditional SA algorithm, and there is an average increase of 1.14%. Compared to the PSO algorithm, the highest increase is 5.19%, and an average increase is 3.52%. (ii) In terms of solving time, the SA algorithm is the longest, the TSA algorithm is second, and the PSO algorithm is the shortest. However, the overall time-consuming gap between the TSA algorithm and the PSO algorithm is not large.

In addition, the required time for the TSA algorithm to solve the same number of customers is basically the same, and it has good stability. The main reasons for the better solution performance of the TSA algorithm are: (i) Due to the characteristics of the multi-models in this paper, the result of the car selection will have a vital impact on the pros and cons of the final solution. The car change operator is helpful to find the most suitable vehicle to perform the task. In addition, a variety of operator coordination methods are beneficial to expand the neighborhood search range and improve the algorithm’s global optimization capability. (ii) A taboo mechanism and a penalty mechanism for violation of constraints are added to the TSA algorithm, the effective tailoring is achieved in the search space, avoiding the algorithm from falling into a local optimum and saving solution time. In contrast, the traditional SA algorithm and the PSO algorithm are easy to fall into a local optimum.