Dynamic Memory Memetic Algorithm for VRPPD With Multiple Arrival Time and Traffic Congestion Constraints

With the new distribution demands emerging continuously in the last decade, the distribution mode is changing gradually in various applications, such as e-commerce, emergency relief supplies distribution, the last-mile delivery, and so on. We formulate vehicle routing problem with pickup and delivery (VRPPD), while simultaneously considering multiple arrival time and traffic congestion constraints. Our model focuses on the clear changes of the distribution mode to meet the fast-delivery requirements in the last decade, which is characterized by the multi-batch arrival of goods in 24 hours and time-varying congestion in various time windows. Besides, we propose dynamic memory memetic algorithm, which updates its dynamic memory by whether to promote populations to find new better solutions or not. This is an effective acceleration mechanism to promote the population progress. Meanwhile, dynamic memory memetic algorithm determines the serious congestion tasks in the delivery route and transforms them into normal congestion or even non-congestion tasks. Test sets with 30 test problems are constructed by using real distribution data from Alibaba Cloud and traffic congestion data from Baidu Map in Shanghai. By comparing with four compared algorithms, the effectiveness, efficiency, and robustness of our proposed algorithm in non-congestion and congestion tests are simultaneously demonstrated.


I. INTRODUCTION
With the new distribution demands emerging continuously in the last decade, the distribution mode is changing gradually in various applications, such as e-commerce, emergency relief supplies distribution, the last-mile delivery, and so on. The modern distribution systems are characterized by their fastdelivery requirements. In the last decade, there are some clear changes of the distribution mode, such as large-scale pickup & delivery tasks, the multi-batch arrival of goods, traffic congestion, and so on. For instance, large-scale pickup & delivery tasks in Double 11 event of Alibaba is about 200 million packages. Moreover, these tasks are constrained by the multi-batch arrival of goods in 24 hours. Meanwhile, traffic The associate editor coordinating the review of this manuscript and approving it for publication was Christian Pilato . congestion in large cities is daily occurring and producing more pollutants. These requirements should be formulated to express new distribution mechanisms.
We explain our research problem as follows. Vehicles pickup goods from depots and deliver them to customers. This is defined as vehicle routing problem with pickup and delivery (VRPPD). Fig. 1 (a) is able to demonstrate VRPPD, its new variants with various constraints [1]- [5], and its research methods [6], [7]. Different from these researches, we focus on the multi-batch arrival of goods and serious traffic congestion. We formulate VRPPD with multiple arrival time and traffic congestion constraints, as shown in Fig. 1 (b). Our model is able to meet the fast-delivery requirements of the modern distribution systems. The differences are as follows.
(i) One is that goods arrive at different times, such as 6:00, 7:00, and 8:00. Thus, one vehicle must take a good from a VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ depot, after this good arrives at this depot. This constraint about the multi-batch arrival of goods is also very common in many other fields, such as disaster relief. (ii) Another is that time-varying congestion leads to the obvious differences of vehicle speed in various time windows. Although the distances among depots and customers are fixed, this results in the remarkable differences of delivery time. Inspired by memetic algorithms, we propose dynamic memory memetic algorithm (DMMA) with local individual and global population searches to address our research problem. Our contributions are as follows. (i) From a local search perspective, single-vehicle and multi-vehicle Hungarian searches are proposed to determine serious congestion tasks and address traffic congestion, which are effective to respond to the challenge of time-varying traffic congestion. Meanwhile, we introduce the individual coding and its repair strategies to deal with the challenge of the multi-batch arrival of goods. (ii) From a global search perspective, a new kind of dynamic memory is proposed to find new better solutions. DMMA updates its dynamic memory by whether to promote populations to find new better solutions or not. This is an effective acceleration mechanism to promote the population progress. (iii) Besides, we use real logistics data from Alibaba Cloud and traffic congestion data from Baidu Map in Shanghai to construct our test sets. Non-congestion results indicate that the effectiveness and efficiency of DMMA are satisfactory. Meanwhile, various-congestion-situation results mean that DMMA is effective to address VRPPD with multiple arrival time and various real-congestion constraints.
The rest of this paper is organized as follows. In Section II, we review related works, summarize our research background, and discuss the real requirements of our model. In Section III, we formulate VRPPD with multiple arrival time and traffic congestion constraints. In Section IV, we introduce DMMA. In Section V, we provide non-congestion and various-congestion-situation results. In Section VI, we present concluding remarks and future work.

II. RELATED WORK
A. VRPPD BACKGROUND AND ITS NEW VARIANTS VRPPD consists of three parts: vehicles, depots, and customers. Vehicles load goods from depots and then deliver goods to customers. Generally, VRPPDs are classified as one-to-one VRPPDs (OTO-VRPPDs), many-to-many VRP-PDs (MTM-VRPPDs), and one-to-many-to-one VRPPDs (OTMTO-VRPPDs). In OTO-VRPPD, there is a one-to-one relationship between customers and depots. For instance, one depot can only supply one customer and one customer can only receive goods from one depot. In MTM-VRPPD, depots supply multiple customers and customers receive goods from multiple depots. In principle, our model belongs to MTM-VRPPD. In OTMTO-VRPPD, vehicles not only need to transport goods to customers, but also need to collect goods from customers and return these goods to depots [8]. For example, breweries provide beer to a group of customers, while simultaneously collecting empty bottles from customers. Our model is different from OTMTO-VRPPD. We review some new variants of OTO-VRPPD and MTM-VRPPD as follows.
In the OTO-VRPPD researches, Treleaven et al. propose SPLICE for a one-to-one pickup and delivery problem, which is usually referred to as stacker crane problem (SCP) [9]. Besides, their research objective is to demonstrate the existence of simple polynomial-time algorithms for SCP with probabilistic optimality guarantees. Sahin et al. consider multi-vehicle one-to-one pickup and delivery problem with split loads (MPDPSL) [10]. Meanwhile, they propose an efficient heuristic that combines tabu search with simulated annealing. Besides, Haddad et al. also focus on MPDPSL and propose an iterated local search metaheuristic and a branchand-price algorithm, which produces high-quality solutions [11]. Naccache et al. consider a multi-pickup and delivery problem with time windows, in which a set of requests is satisfied by a fleet of vehicles [12]. In each request, items are required to be picked up from different locations to be shipped and unloaded at one common delivery location. This problem is characterized by multi-pickup and one common delivery. Besides, they solve this problem by using a hybrid large neighborhood search heuristic. Bartolini et al. present traveling salesman problem with pickup, delivery, and ride-time constraints (TSPPD-RT) [13]. Besides, they analyze TSPPD-RT with and without capacity constraints. TSPPD-RT is characterized by that each request's ride-time does not exceed its maximum ride-time.
In the MTM-VRPPD researches, Ahkamiraad and Wang consider the mixed-integer linear programming model of a special type of capacitated and multiple cross-docked vehicle routing problem with pickup, delivery, and time windows [14]. This model is characterized by multiple cross docks and each cross dock has a limited maximum capacity. Besides, they propose a hybrid of genetic algorithm and particle swarm optimization. Chen et al. formulate the paired many-to-many pickup and delivery problem, which is characterized by its paired demands between customers [15]. Besides, they use library vehicle routing as an application and test their two-stage solution algorithm by using real data of San Francisco library system. Ghaffarinasab et al. formulate planar hub location-routing problem by using a continuous approximation technique [16]. The research objective is to decide on the locations of hubs & service regions and minimize the approximate total transportation cost, including local pickup and delivery costs. Besides, an iterative Weiszfeld-type algorithm and a particle swarm optimization metaheuristic are proposed. Liu et al. propose a fast decomposition and reconstruction framework to solve the pickup and delivery problem with time windows and last-in-first-out loading [17]. Moreover, they use the adaptive memory mechanism to provide a better start when populations fall into the local optima. Zhou et al. develop a twostage heuristic algorithm for green vehicle routing problem which focus on dual services with stochastic travel times [18]. Experimental results indicate that the proposed algorithm is suitable to address this problem, especially for providing sustainable and diversified delivery services. Wang investigates and proposes two heuristic algorithms based on tabu search and simulated annealing algorithm to address the heterogeneous fleet location routing problem with simultaneous pickup and delivery [19]. Focusing the generation of a priori routes for a fleet of vehicles, Zhu and Sheu propose a failure-specific cooperative recourse strategy to explore a risk pooling mechanism for routing in the context of simultaneous pickup and delivery with stochastic demands [20]. Marinaki and Marinakis propose combinatorial-neighborhoodtopology glowworm swarm optimization for successfully solving the vehicle routing problem with stochastic demands [21]. Moreover, they use capacitated vehicle routing problem and vehicle routing problem with stochastic demands to demonstrate the effectiveness of the proposed algorithm. Ruiz et al. consider the open vehicle routing problem with split deliveries, propose a new two mixed-integer formulations of the problem, and design a cutting-plane method to improve the optimization performance [22].

B. VRPPD APPLICATIONS WITH VARIOUS CONSTRAINTS
In the last decade, emergency relief supplies distribution, e-commerce, and last-mile logistics develop rapidly, such as Amazon. These VRPPD applications consider many real constraints, such as load, fuel consumption, dangerous goods transportation route, delivery time, pollution emissions, time window, and so on.
We review some existing applications of VRPPDs as follows. Yanik et al. focus on capacitated vehicle routing problem with multiple pickups, single delivery, and time windows [23]. Moreover, they propose a hybrid metaheuristic approach and a modified savings algorithm. Besides, they test their method by using a real dataset in Istanbul. Guo and Liu address time-dependent vehicle routing of free pickup and delivery service to deliver a set of airline passengers, while considering traffic congestion, time window, and maximum ride time constraints [24]. Focusing on on-demand transportation systems and ride-sharing services, Mahmoudi and Zhou propose a new time-discretized multi-commodity network flow model based on the integration of vehicles' carrying states [25]. This model is characterized by threedimensional state-space-time network construction. Besides, their dynamic programming approach is tested by using the Chicago sketch and Phoenix regional transportation networks. Zheng et al. focus on the influence of time-dependent velocity and propose a realistic model to solve urban pickup and delivery problem with considering the time-dependent fuzzy velocity of vehicles [26]. Besides, they introduce a real pickup and delivery problem with 28 customers in a logistics company. Xiao and Konak propose a new mixedinteger linear programming model for green vehicle routing and scheduling problems [27]. This model considers heterogeneous vehicles, time-varying traffic congestion, customer/vehicle time window constraints, the impact of vehicle load on emissions, and vehicle capacity/mileage constraints. Sun et al. introduce a new time-varying pickup and delivery problem, which considers the time window constraints of pickup and delivery points and traffic congestion [28]. In addition, a restricted dynamic programming heuristic algorithm is proposed that finds the best solution for 32 out of 34 instances. Zhang et al. equivalently transform the collection vehicle routing problem of the garbage facilities to the vehicle routing problem with simultaneous pickup-delivery and time windows by introducing dummy customers [29]. Besides, they propose a parallel simulated annealing algorithm to address the real application problem in the Xuanwu District of Beijing. Considering the vehicle routing problem with simultaneous pickup and delivery and handling costs, Hornstra et al. propose an adaptive large neighborhood search metaheuristic and improve 39 out of 54 best known solutions for this problem [30]. Peng et al. study the multiple vehicle pickup and delivery problem that considers the constraints of vehicle route time, vehicle capacity, and vehicle last-in-firstout [31]. Besides, they propose a learning-based memetic algorithm and improve the previous best known results for 132 out of 158 problem instances. Anwar and Younas formulate the pickup and delivery problem as a many-objective pickup and delivery problem with delay time of vehicle having six criteria to be optimized [32]. Moreover, they propose a memetic improved decomposition based evolutionary algorithm and use a variety of small, medium, and large-scale problems to test the proposed algorithm.

C. SUMMARY AND DISCUSSIONS
These previous research lines are summarized in Table 1

III. PROBLEM STATEMENT AND FORMULATION
We formulate VRPPD with multiple arrival time and traffic congestion constraints as follows.
First, we explain the relationship among depots, customers, congestion coefficients, and vehicle speeds as follows. (i) Depot and customer set DC is The distances among depots and customers are expressed by using distance matrix DIS in (1). DIS(DC l , DC m ) is the element of the lth row and mth column, which is the distance between DC l and DC m . (iii) Congestion coefficient set J is a 3-D matrix (i.e., NDC × NDC × NTW ) to express congestion situations in various time windows. NTW is the number of time windows and J (:, :, TW k ) is a 2-D matrix to define the congestion situation in the kth time window TW k . For example, J (DC l , DC m , TW k ) is the congestion coefficient of DIS(DC l , DC m ) in TW k , which is from 1 to +∞. The greater congestion coefficient, the more serious congestion. (iv) We assume that all vehicles have the same non-congestion speed S, which is equal to 60km/h. The real speed of one vehicle from DC l to DC m in TW k is given in (2) by using the definition of congestion coefficient [33]. Note that, (2) means the more serious congestion, the lower real speed. The definition of congestion coefficient J (DC l , DC m , TW k ) reflects its average condition in TW k . Moreover, TW k is the time window when the vehicle leaves DC l .
· · · · · · · · · · · · DIS l−1 · · · DIS l−m · · · DIS l−NDC · · · · · · · · · · · · · · · DIS NDC−1 · · · DIS NDC−m · · · DIS NDC−NDC Secondly, we define the task allocation relationship by using task allocation matrix TA. Task set T is {T 1 , T 2 , . . . T j , . . . T NT } and NT = |T |. T j is the jth task, which is equal to {D p , C q , TW k }. This indicates that this good arrives at D p in TW k and should deliver to C q . Besides, Then, TA is given in (3) and each TA i−j ∈ {0, 1}. If TA i−j = 1, this means that T j is allocated to V i . Otherwise, this means that T j is not allocated to V i . In addition, (4) expresses that each task can be executed once. (4) and (5) mean that all tasks must be allocated.
Thirdly, we give the feasible task allocation solution COD as follows. The feasible task allocation solu- is the nth depot or customer that V i should arrive at. For example by using Fig. 2, COD 1 is the V 1 s task allocation solution, which is equal to This means that V 1 sequentially reaches to D 1 , C 1 , D 2 , and C 2 . Note that, each task in Fig.2 has its own arrival time, which is only able to be picked up after its arrival time. Lastly, our objective function (i.e., fitness function in DMMA) is to minimize the delivery time of all vehicles for all tasks, which is given by where DT i is the delivery time of all V i 's tasks. Besides, DT i−n is the delivery time from COD i−(n−1) to COD i−n .

IV. DYNAMIC MEMORY MEMETIC ALGORITHM A. MOTIVATION
Our research objective is to demonstrate that DMMA is effective to address VRPPD with multiple arrival time and traffic congestion constraints. Besides, another objective is to demonstrate the effectiveness and efficiency of DMMA, even in VRPPD with multiple arrival time and non-congestion constraints. VOLUME 8, 2020 Algorithm 1 Dynamic Memory Memetic Algorithm Input: population size PopuSize, dynamic memory size DMSize, the maximum distance of the memetic search K max , the minimum number of shaking S min , the maximum number of shaking S max , the mutation probability p m , the initial vitality V initial , the vitality limit V limit , the congestion-search threshold CSThr. Initialize the population by using Algorithm 2 and evaluate the initial population. Initialize dynamic memory by using the best individuals of the population. Define TempInd as the temporary-individual variable.
Define k as the current search distance. REPEAT the following steps, until the termination condition is satisfied. FOR each individual in the population TempInd ← the selected individual T1 from the population.
FOR k ← 1 to K max Select an individual from dynamic memory randomly, which is named as D1.
Shaking in Algorithm 3 for TempInd and obtain the individual T2.
Local search in Algorithm 5 for T2 and obtain the individual T3.
Update D1 in dynamic memory by using Algorithm 9.
If T3 is better than TempInd, then TempInd ← T3 and break for-loop. Otherwise, k ← k + 1. END END Except for the best individual, perform mutation in Algorithm 4 according to the mutation probability p m .
Update the number of using the fitness functions, and check the termination condition. END Output the best individual as the best solution.

B. FRAMEWORK
The framework of DMMA is given in Algorithm 1. Comparing with memetic algorithms [31], [32], [34], the main differences are as follows. (i) We design dynamic memory to update the population and this is effective to prompt the global-population-search progress. (ii) Our local search is to address congestion by using the local-individual-search method.
First, one feasible task allocation solution is also the coding of one individual in the population, as shown in Fig. 2.
In this paper, one coding corresponds to one individual in the population. Moreover, various operator processes may lead to infeasible solutions. There are three kinds of infeasible solutions and their repair strategies, as shown in Fig. 3. These cases are only related to the single-vehicle task sequences and details are as follows. Fig. 3 (a) means that V i must Algorithm 2 The Population Initialization Input: population size PopuSize Step 1: Obtain all tasks in the first batch and randomly allocate these tasks to each vehicle.
Step 2: Repeat Step 1 for the following batches, sequentially.
Then, obtain the coding of one individual.
Step 3: Execute 10 times of swap operators in Fig. 4. If found a better one, update this individual immediately.
Step 4: Repeat Steps 1-4 until all individuals are generated.
firstly pick up and then deliver. Fig. 3 (b) indicates that D 2 is not reasonable, because V i is not allocated to accomplish T 2 . Fig. 3 (c) indicates that V i is not necessary to arrive at D 1 twice only for delivering C 1 . Besides, if one infeasible solution is found, we should use the corresponding repair strategies to obtain a feasible solution. Secondly, the population initialization is given in Algorithm 2, which randomly allocates tasks to vehicles and uses swap operator in Fig. 4 to optimize a single-vehicle task sequence.
Thirdly, shaking in Algorithm 3 is used to inject the continuously-updated genes of dynamic memory into the population. In each time of shaking, a k-length slice is chosen from a randomly selected individual of dynamic memory. Besides, the shaking illustrations are given in Fig. 5. Note that, k is equal to the current search distance.
Lastly, mutation in Algorithm 4 uses swap operator to improve the population diversity.  . The process of swap operator, which randomly selects a segment (i.e., from DC 2 to DC 5 ) from a single-vehicle sequence. If the swapped segment (i.e., from DC 5 to DC 2 ) is better than before, then replace the randomly selected segment with the swapped segment.

Algorithm 3 Shaking
Input: the randomly selected individual D1 from dynamic memory, the current search distance k, the temporary-individual variable TempInd, the minimum number of shaking S min , the maximum number of shaking S max .
Step 1: Define a temporary vector as TV and TV ← TempInd.
Then, randomly generate an integer number S n , which is from S min to S max .
Step 2: Repeat the next steps, until the S n times of shaking are accomplished.
Step 3: Randomly selected a single-vehicle sequence with more than k genes from D1, as shown from Fig. 5 (a) to Fig. 5 (b). Then, randomly selected a segment with k genes as Seg, as shown from Fig. 5  multi-vehicle Hungarian searches are used to optimize the single-vehicle task sequences and multi-vehicle task sequences, respectively. Besides, Fig. 6 explains how to determine a serious congestion task. This is frequently used in local search. As shown in Fig. 6, we select the most serious congestion paths from a single-vehicle task sequence according to congestion coefficient set J , obtain the corresponding serious congestion tasks, and randomly choose one as the selected serious congestion task.

Algorithm 5 Local Search
Input: the individual T2 Step 1: Perform single-vehicle search in Algorithm 6 for T2 and obtain the individual T2-A.
Step 2: Define a temporary individual as TI. If T2-A is better than T2, then TI ← T2-A. Otherwise, TI ← T2.
Step 3: Perform multi-vehicle Hungarian search in Algorithm 7 for TI and then obtain the individual T2-B.
Step 4: If T2-B is better than TI, use T2-B as output. Otherwise, use TI as output.

1) SINGLE-VEHICLE SEARCH
Single-vehicle search in Algorithm 6 uses swap and congestion shift operators. In Algorithm 6, SVN and PON are defined only to adjust the search efforts, which are equal to 10 and 4, respectively. Besides, congestion shift operator is used to transform the serious congestion tasks into normalcongestion or even non-congestion tasks, as shown in Fig. 7. VOLUME 8, 2020 FIGURE 7. The process of congestion shift operator, which randomly selects a serious congestion task (i.e., DC 3 ) and randomly moves this task. If found a better single-vehicle task sequence, then update.

2) MULTI-VEHICLE HUNGARIAN SEARCH
Multi-vehicle Hungarian search in Algorithm 7 uses Hungarian algorithm [35] in Algorithm 8. Besides, Fig. 8  Step 2: According to CMSize, randomly select the singlevehicle task sequences from TI and obtain the selected vehicle set SelVehicleSet. Then, use the following process to select the tasks from these sequences and obtain the selected task set SelTaskSet. If rand() < CSThr, select a serious congestion task. Otherwise, select a task randomly.
Step 3: According to these selected sequences and tasks, obtain the consumption matrix CM.
If found a better solution, update TI. Otherwise, don't update.
to find new better solutions or not. As shown in (8) and (9), the individual vitality value varies from (V initial -V limit ) to (V initial + V limit ). V initial is the initial vitality, V limit is the vitality limit, and they are integers. When an individual in dynamic memory (like D1 in Algorithm 1) promotes one individual in the population (like TempInd and T3 in Algorithm 1) to find new better solutions, then we increase the individual   vitality V d by +1. Otherwise, we decrease this individual vitality by −1.
Vitality selection is used to update dynamic memory in Algorithm 9. Fig. 9 is used to explain the meaning of the individual-vitality trend as follows. (i) In Fig. 9 (a), one individual vitality continuously increases. This kind of TABLE 3. Test sets used real distribution data from Alibaba Cloud in Shanghai. Besides, ND is the number of depots, NV is the number of vehicles, and NT is the number of tasks. Moreover, the arrival time of the multi-batch goods are 7:00, 9:00, 11:00, 13:00, 15:00, 17:00 and 19:00, respectively.

Algorithm 8 Hungarian Algorithm
Input: the consumption matrix CM, the selected task set SelTaskSet, the selected vehicle set SelVehicleSet. Output: the output result that assigns tasks (row) to vehicles (column).
Step 1: In each row of CM, subtract the smallest element of the row.
Step 2: In each column of CM, subtract the smallest element of the column.
Step 3: Draw lines through the rows and columns that have the 0 elements by using the fewest lines.
Step 4: If there are CMSize lines drawn, then the optimal assignment of 0 elements is the assignment scheme as the output result. The 0 element in column h of line g represents we should allocate the SelTaskSet h to the SelVehicleSet g . Otherwise, go to Step 5.
Step 5: Find the smallest element in CM, which is not covered by any line. Subtract this element from each row that isn't covered and add it to each column that is covered. Then, go to Step 3.
individual-vitality trend lines is the evidence to prove that this individual is effective to promote the population progress. Or to say, this kind of task segments that are inserted into the population is effective to prevent populations from trapping in local optima. Note that, these task segments are inserted as a whole (like Fig. 5). Thus, this changes the arrangement layout of the inserted individuals in the population. (ii) Otherwise, when one individual vitality continuously decreases (like

Algorithm 9 Vitality Selection to Update Dynamic Memory
Input: the individual D1 in dynamic memory, the individual TempInd and the individual T3, the initial vitality V initial , the vitality limit V limit Step 1: According to TempInd and T3, update the vitality value of the individual D1 by using (8) or (9).
Step 2: If the vitality value of the individual D1 is equal to (V initial -V limit ), then replace D1 with the selected individual from the population. Note that, there are two requirements for this selected individual as follows. One is that dynamic memory does not include this selected individual before replacing D1. The other is that, if there are more than one individual that satisfy the previous requirement, then select the best fitness individual from these individuals. Fig. 9 (b)), this kind of individual-vitality trend lines indicates that this individual is not effective to promote the population progress. Then, when its vitality is equal to (V initial -V limit ), this individual is discarded.
2) PRINCIPLE ANALYSIS Fig. 10 demonstrates the function of dynamic memory and explains how to establish an effective acceleration mechanism for promoting the population search. This acceleration mechanism is goal-oriented by using its two kinds of updating criteria as follows. (i) According to fitness, one search effort updates the population. (ii) According to whether to promote populations to find new better solutions or not, another search effort updates dynamic memory. VOLUME 8, 2020 Explanations are as follows. (i) The premature phenomena of populations are one of the major obstacles in evolutionary computation. Different schemes of protecting population diversity are widely used that do not be effective to find new better solutions in all cases. For example, injecting new genes into populations is not able to ensure the population to find new better solutions, definitely. (ii) However, from an updating-criterion perspective, better individuals in populations usually enter dynamic memory in the early stage of population evolution. This is effective to promote individuals to make progress. Then, in the later stage, better individuals do not enter dynamic memory, because dynamic memory has included them. Thus, many individuals with gene diversity and poor fitness could enter dynamic memory. Different from better individuals, they may include various task segments. These task segments could improve the probability of finding new better solutions.

E. COMPUTATIONAL COMPLEXITY
The computational complexity of the proposed model with DMMA in (10) is dependent on the test set size and DMMA's parameters [37], such as the number of depot and customer set NDC, the number of vehicle set NV, termination generation TG, population size PopuSize, the maximum distance of the memetic search K max , the maximum number of shaking S max , and the maximum size of consumption matrix CMMaxSize. In (10), T () returns the computational complexity of one algorithm's item. Besides, Table 2 gives the computational complexity of each item in (10). Thus, the total time complexity of the proposed model with DMMA could be obtained by using (11).

V. EXPERIMENTS AND DISCUSSIONS A. EXPERIMENTAL DESIGN 1) TEST SETS
We used distribution data from Alibaba Cloud [38] in Shanghai and constructed test sets with 30 test problems, as shown in Table 3. Congestion coefficient sets were from Baidu Map [39] on Monday, Wednesday, and Sunday. Test and congestion coefficient sets could be downloaded from [40].

2) COMPARED ALGORITHMS
We used MAPSO [41], LSVNS [42], MTWPS [43], and BRKGA [44] as compared algorithms. Because our model is a new variant of VRPPD, the coding method and repair strategies of infeasible solutions are different from other researches. Therefore, to compare with DMMA objectively, these compared algorithms used the same coding method in Fig. 2, the population initialization in Algorithm 2, and the repair strategies in Fig. 3.

3) PARAMETER SETTING B. NON-CONGESTION RESULTS
We tested the effectiveness and efficiency of DMMA under the non-congestion condition. Note that, the non-congestion condition means that all elements in congestion coefficient set J were equal to 1.  BEST is the best solution for 20 runs. Besides, the percentage deviations of the average solution (PDav) and the best solution (PDbest) are given in (12) and (13). The best known To further demonstrate the effectiveness of DMMA, we used Wilcoxon test at 0.05 level of significance for the best-fitness data of 20 runs between one compared algorithm and DMMA, as shown in Table 5. The symbols of (+), (≈), and (−) mean that DMMA is better than, approximately equal to, and worse than the compared algorithm, respectively. DMMA obtained the best performance in all cases. The reason is that DMMA uses the effective acceleration mechanism for promoting the population search, as shown in Fig. 10. Fig. 13 presents the efficiency results of the average computation times (ACTs) for 20 runs. The smaller ACT, the better algorithm efficiency. These results in Fig. 13 indicate that DMMA has lower computation costs, especially for large scale test set.

2) ON ALGORITHM EFFICIENCY
To further demonstrate the efficiency of DMMA, we also used Wilcoxon test at 0.05 level of significance, as shown in Table 6. DMMA obtained the best performance in most cases. The reason is that the computation costs of singlevehicle search, multi-vehicle Hungarian search, and vitality selection are lower, such as the shift operator in Algorithm 6, the swap operator in Algorithm 6, and the Hungarian algorithm only with a short task sequence (i.e., CMSize) in Algorithm 8, updating the individual vitality (i.e., (8) and (9)).

1) DELIVERY TIME AND SENSITIVITY ANALYSIS
The best fitness in (6) is the delivery time, which is the most important index to express the pickup & delivery efficiency of real distribution systems. The shorter delivery time, the better pickup & delivery efficiency. Fig. 14 provides the delivery time results. In Fig.14, the sum of MEANs for various conges-tion situations in 20 runs is used as the evaluation criterion, while each kind of MEANs for each congestion coefficient set is also given. DMMA obtained the best performance in all cases. The reason is that DMMA uses effective schemes to determine serious congestion tasks and address traffic congestion, as shown in Figs. 6, 7, and 8. These results in Fig. 14 indicate that DMMA is robust and effective for VRPPD with multiple arrival time and various congestion constraints.
To demonstrate the robustness of DMMA for various distribution scenarios and various congestion situations, the sensitivity-analysis explanations were as follows. (i) From a distribution perspective, there were various kinds of depot's numbers (i.e., 3,7,12,18,24,30), vehicle's numbers (i.e., 3,4,5,8,9,10,20,25,30,45), and task's numbers (i.e., 35, 45, 55, . . . , 484) in 30 test problems, as shown in Table 3. (ii) From a congestion perspective, we used various congestion coefficient sets, as shown in Fig. 15. Three kinds of congestion situations were the representatives of our daily traffic congestion in one week, because they corresponded to the traffic congestion situations of Monday, Wednesday, and Sunday. These indicate that DMMA is robust, even for the changes of distribution scenarios and congestion situations. Fig. 16 gives the sum of ACTs for various congestion situations in 20 runs, while each kind of ACTs for each congestion coefficient set was also given. DMMA obtained the best results in most cases. Generally, ACTs of DMMA for various congestion situations were similar to non-congestion ACTs. For example, ACTs of DMMA for JSet1 were roughly the same with non-congestion ACTs in Fig. 13.

2) COMPUTATION TIME AND ALGORITHM USABILITY
On the other hand, we discussed the usability of DMMA as follows. The maximum value of all ACTs in Fig. 16 was 158.88min for P25. In P25, there are 24 depots, 30 vehicles, and 472 tasks. From a computation efficiency perspective, this ACT's maximum value is small enough for the size of this test problem. From an application perspective, this means that DMMA is able to complete task allocations within acceptable time, because real distribution scenarios are able to obtain all pickup & delivery tasks' information at least one day in advance.
To further demonstrate the usability of DMMA in real applications, Fig. 17 shows the real cases by using map route and their convergence curves. These maps were all based on real pickup & delivery data and real congestion data. In Fig. 17 (a) and (c), each kind of the color paths corresponded to one vehicle and each kind of line types corresponded to one kind of congestion degree. Meanwhile, Fig. 17 (b) and (d) indicate that the convergence performance of DMMA is satisfactory. Generally speaking, Fig. 17 displays a fleet pickup & delivery task allocation and this implies that DMMA is able to be used in the real distribution cases.

VI. CONCLUSION
We formulate VRPPD with multiple arrival time and traffic congestion constraints. Moreover, our model is able to express the fast-delivery requirements of the modern distribution systems in our production and life, such as largescale emergency relief supplies distribution. As the extensive results shown, DMMA is able to effectively address the real pickup & delivery application problems with large-scale tasks, the multi-batch arrival of goods, and time-dependent traffic congestion. The robustness and usability of DMMA are proved by non-congestion and congestion tests, which are based on the real distribution data and traffic congestion data in Shanghai.
Recently, the deployment of unmanned pickup & delivery platforms gradually increases in the VRPPD applications. Our future researches focus on online pickup-delivery task allocations, which are able to respond to the newly added tasks and real-time traffic congestion.