Real Environment-Aware Multisource Data-Associated Cold Chain Logistics Scheduling: A Multiple Population-Based Multiobjective Ant Colony System Approach

Cold chain logistics (CCL) scheduling is important for smart cities as it directly affects the service quality and operating profits of logistics companies. However, traditional CCL models seldom reflect the real transportation environment, making the solutions hardly applicable to the real CCL scenes. Hence, this paper attempts to establish a multisource data-associated CCL model oriented to the real transportation environment. This environment is considered by employing the real-captured driving duration and distance between any two places. Three scheduling objectives (namely, quality losses, personnel and vehicle costs, and transportation costs) are taken into account. To efficiently solve the proposed multisource data-associated multiobjective CCL model, a multiple population-based multiobjective ant colony system (MPMOACS) approach is proposed. Based on the multiple populations for multiple objectives framework, the MPMOACS approach can optimize multiple objectives sufficiently, and thus obtain promising solutions distributed along the entire Pareto front. To further enhance the performance of the MPMOACS, a ranking-based local search strategy is also designed. Experiments are conducted on not only the existing benchmark instances but also a real environment-aware multisource dataset that is built based on real-captured transportation data of Guangzhou and Shenzhen, China. Compared with six state-of-the-art and very recent well-performing multiobjective optimization approaches, the proposed MPMOACS approach exhibits the overall best performance.


I. INTRODUCTION
C OLD chain logistics (CCL) is concerned with the transportation of perishable goods that require lowtemperature storage [1], [2], [3]. It has attracted increasing research attention in recent years since smart logistics forms a fundamental part of smart cities [4], [5], [6]. Nevertheless, the quality of service in CCL is still limited by the quality losses of goods and expensive distribution costs. Specifically, the acute quality losses of goods that occur during the long distribution duration may lead to serious waste and even safety problems [7], [8]. Furthermore, the high distribution costs may weaken the competitiveness of the logistics companies. Hence, effective CCL scheduling should be conducted to optimize both quality losses and distribution costs.
Having developed from the vehicle routing problem (VRP), most of the existing CCL models derive driving duration by assuming a constant driving speed [9], [10], but the driving speed between two places is not a constant since it varies with time (i.e., the driving speed is time-dependent). For example, the driving speed in peak hours is much lower than that in other hours. Thus, the scheduling results of the existing CCL models are often inadequate to a real-world environment. To integrate the time-dependent speed into the VRP model, researchers have studied the time-dependent VRP (TDVRP) [11], [12], [13], [14], [15]. Nevertheless, only a simple and uniform step function is used to describe the time-dependent This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ driving speed in the TDVRP [16], [17], which is insufficiently accurate to reflect the real transportation environment. Besides, the TDVRP hypothesizes that there is a straight path between any two places and adopts the corresponding Euclidean distance to calculate the driving duration. However, such a hypothesis does not necessarily hold due to the complex urban layout. Therefore, the TDVRP still has limitations and is inefficient when the scheduling results are returned to the real CCL scenes. In other words, although the time-dependent speed is considered, the TDVRP model is not suitable for real-world CCL scheduling.
To address the problem that the solutions of the traditional CCL models may fail to work in real CCL scenes, a real environment-aware multisource data-associated CCL model is developed in this paper. Multisource data refer to the data obtained from multiple sources, consisting of order data, service data, and transportation data. The first two are gathered from customers and logistics companies, respectively, while the transportation data are captured through map API (e.g., the Baidu map API of China). In detail, the order data include customers information; the service data include information about personnel and vehicle resources; and transportation data include the driving distance, duration, and speed in real environments. By associating order data, service data, and transportation data in our model, the gap between the problem model and the real CCL scenes can be narrowed, and thus the feasibility of scheduling results can be improved.
To satisfy the requirements of different stakeholders, the proposed CCL model considers three objectives, resulting in the multisource data-associated multiobjective CCL model. The first objective is to minimize the quality losses of goods. The second objective is to minimize the personnel and vehicle costs. The third objective is to minimize the transportation costs measured by the total driving distance and driving duration.
Recently, many scheduling approaches for solving various multiobjective VRP models have been proposed [18], [19], [20]. For example, a coevolutionary constrained multiobjective optimization (CCMO) framework to solve a multiobjective VRP with time windows (VRPTW) model that aims to minimize the number of vehicles and the total traveled distance was proposed in [21]. A multi-depot VRPTW model with five objectives was solved by a two-stage multiobjective evolutionary algorithm (TS-MOEA) [22]. The first stage of the TS-MOEA finds extreme solutions using the extremized crowded nondominated sorting genetic algorithm II (EC-NSGA-II) [23], which is developed from the nondominated sorting genetic algorithm II (NSGA-II) [24]. The second stage of the TS-MOEA obtains more Pareto solutions based on extreme solutions through a multiobjective evolutionary algorithm based on decomposition (MOEA/D) [25]. However, the quality losses of goods are not considered in these multiobjective VRP models. Hence, the scheduling approaches proposed for these multiobjective VRP models may not obtain satisfying results when being used to solve the multisource data-associated multiobjective CCL model.
In addition, multiobjective CCL scheduling has seldom been studied in previous work, which indicates a lack of existing multiobjective scheduling approaches to multiobjective CCL scheduling problems. Note that although multiple objectives are considered in the CCL model proposed in [10], their work combined multiple objectives into one objective and solved the model using a single-objective optimization approach. In such case, the multiple objectives may not be sufficiently optimized. Therefore, it is necessary to develop an effective multiobjective scheduling approach to efficiently optimize the multisource data-associated multiobjective CCL model.
Owing to the promising performance of ant colony optimization (ACO) in solving the VRP [26], [27], [28], [29] and other combinatorial optimization problems [30], [31], [32], [33], a widely used ACO variant, i.e., the ant colony system (ACS) [34], [35] is regarded as an appropriate approach for solving our model. However, the simultaneous optimization of multiple objectives is difficult due to the conflict among objectives, which challenges many existing multiobjective techniques [36], [37]. Thus, the famous and popular multiple populations for multiple objectives (MPMO) framework proposed by Zhan et al. [38] is employed and combined with the ACS, resulting in a multiple population-based multiobjective ACS (MPMOACS) approach. Furthermore, a ranking-based local search (RBLS) strategy is designed to further enhance the search efficiency of MPMOACS.
Including the problem modeling aspect and the algorithm designing aspect, the contributions of the work reported in this paper are as follows.
First, a multisource data-associated multiobjective CCL model oriented to the real transportation environment is established. By using real transportation data, our model can reflect the real-world situation more accurately, and thus improve the feasibility of scheduling results. Furthermore, three objectives (i.e., quality losses, personnel and vehicle costs, and transportation costs) are considered in the proposed model, which can satisfy various real-world demands.
Second, the new MPMOACS approach is proposed to solve the proposed model. According to the suggestion of the MPMO framework [39], [40], the MPMOACS employs three populations. In this way, all the three objectives can be sufficiently searched to obtain promising solutions distributed along the entire Pareto front.
Third, the RBLS strategy is designed to strengthen the performance of the MPMOACS. By combining the ranking of all objectives and an effective local search operator, the RBLS strategy can adjust the distribution sequence of customers reasonably and thus improve the solution quality.
The remainder of this paper is organized as follows. Section II describes and formulates the multisource dataassociated multiobjective CCL model. Section III develops the MPMOACS approach. Section IV illustrates the experiments, with results discussed. Finally, Section V concludes the work.

II. MULTISOURCE DATA-ASSOCIATED MULTIOBJECTIVE CCL MODEL
A. Model Description 1) Model Information: To begin with, all variables involved in this model are presented in Table I. Given one depot and n  TABLE I   TERMINOLOGY NOTATION customers, the search space of the multisource data-associated multiobjective CCL model can be represented as a directed graph, G = (V, A). V = {0, 1, 2, . . . , n, n + 1} represents the vertex set of depot and customers. Herein, 0 is the depot vertex that vehicles depart from, and n + 1 is the depot vertex that vehicles return to after finishing distribution tasks, both of which refer to the same depot but use different indexes for clarity. Besides, the remaining vertices represent customers, and the customer vertex set is denoted as N (N ⊆ V ). The set of all edges between any two vertices is denoted as A.
Three objectives are considered in this model, which are denoted as quality losses, personnel and vehicle costs, and transportation costs. Meanwhile, the hard time window constraint and capacity constraint are considered in this model. Given a customer i and its time window [b i , e i ], the service of i cannot begin before b i or after e i . Furthermore, each vehicle departs from the depot at b 0 and must return to the depot before e n+1 , which can be regarded as the hard time window constraint of the depot. To satisfy the capacity constraint, the total demand of goods distributed by a vehicle cannot exceed the maximum capacity of this vehicle.
2) Multisource Data Association: As our CCL model is oriented to a real transportation environment, three types of data are used in our CCL model. The first type of data is order data that contain all necessary information of each customer, including their location, demand, time windows, and service duration. The second type of data is the service data, including the number of vehicles that can be used, the maximum capacity of each vehicle, the fixed cost of vehicle usage, the salary cost of unit time, and the deterioration rates of goods. The third type of data is transportation data, including the driving distance, driving duration, and driving speed. These three types of data are regarded as multisource data because they are collected from different sources. Multisource data association is a promising way to build a more practical model [41]. Therefore, the associations among these three types of data are constructed when they are used as the input of our model.

3) Real Transportation Environment Consideration:
To introduce the real transportation environment information in our model, the driving duration between any two places at different times and the driving distance between any two places are captured. Then, the average captured driving duration in one hour is adopted to calculate the driving speed in that hour. In this way, the step function between the departure time and driving speed can be constructed, and the driving duration from one place to another place at any time can be derived through the step function.
Specifically, the step function based on real-captured transportation data can be formulated as where sp i j (t) is the driving speed from vertex i to vertex j at time t (t ∈ [0, 24)). C i j contains the constant speeds from i to j of each hour in a day, i.e., {C i j (0), C i j (1), …, C i j (23)}, which is calculated based on real-captured data. For clarity, an example is provided in Table S.I in the supplementary material and Fig. 1 to help explain the step function. As a part of Table S.I, Table II shows the captured driving duration between two places during the first 9 hours in a day. The driving duration is captured twice every hour in the example and the captured data are shown as the columns named num1 and num2, respectively; and the average duration is shown as the column named Average. Given the real distance, the timedependent speed in each hour can be obtained, leading to the step function shown in Fig. 1 for each hour in a day. Moreover, unlike the TDVRP in which any two places share the same step function, each pair of places has its step function in our model. Note that the step function can be established between the departure time and either the driving speed or driving duration. The driving speed is chosen in our model because it can ensure the satisfaction of the first-in-first-out (FIFO) property [17]. Herein, the FIFO property refers to the fact that given the same origin and destination, a vehicle that departs later cannot arrive at the destination before a vehicle that departs earlier.

B. Model Formulation
Three objectives are considered in the multisource dataassociated multiobjective CCL model. The first objective is the total quality losses of goods when customers receive them, which is calculated as where q i is the demand of customer i , and ql i is the quality losses of unit goods for customer i . Since the quality of goods at time t is always measured as Q(t) = Q(0) × e −d×t (Q(0) = 1, d is the deterioration rate) [42], the delivered quality losses of unit goods for customer i is measured as where the distribution duration of goods for customer i is decomposed into the duration of driving and waiting (i.e., t 1 (i )) and the duration of serving other customers that are served by the same vehicle with i and are served before customer i (i.e., t 2 (i )). In addition, the deterioration rates of t 1 (i ) and t 2 (i ) are denoted as d 1 and d 2 , respectively. Two different deterioration rates are considered because the vehicle carriage is opened when serving customers, which can result in a higher deterioration rate when compared with the situation in which the vehicle carriage is closed (i.e., the vehicle is driving or waiting). In this way, the longer distribution duration of perishable goods (including both the driving and waiting duration t 1 (i ) and the service duration t 2 (i )) will lead to a larger value calculated by (3), reflecting worse quality losses. The second objective is the personnel and vehicle costs, which are measured by the sum of the fixed costs of vehicle usage and the salary costs for the distribution duration. Thus, the formulation of the second objective is where the first term of (4) is the fixed costs of vehicle usage, the second term of (4) is the salary costs, and the distribution duration of vehicle k is computed as the difference between s (n+1)k and s 0k . Herein, c 0 is the fixed cost of a vehicle in use, and c 1 is the salary cost of unit distribution duration. The third objective is the transportation costs that depend on the driving distance and driving duration of vehicles. Hence, the third objective can be measured as where the first term is used to compute the driving distance, and the second term is used to compute the driving duration. The variable l ik refers to the time that vehicle k departs from i .
In summary, f 1 is the minimization of the quality losses, while f 2 and f 3 are related to the distribution costs. It should be noted that given enough vehicles, each customer can be served at b i (1 ≤ i ≤ n) (i.e., f 1 is minimized), but it may lead to too much waiting time, expensive vehicle costs, and redundant driving distance and duration. That is, the minimization of f 1 may lead to increases in f 2 and f 3 . Hence, f 1 conflicts with f 2 and f 3 . Furthermore, with the constraint of hard time windows, the total distribution duration of a path with the shortest driving distance may be longer than that of a path with a longer driving distance if the waiting duration of the former is much greater than that of the latter. In this case, the minimization of f 2 may result in an increase in f 3 and vice versa. Hence, f 2 conflicts with f 3 to some extent.
To optimize the three objectives mentioned above, the multisource data-associated multiobjective CCL model can be formulated as where (6) is the objective of the CCL model. Constraint (7) means that each customer can only be served once by one vehicle. Constraints (8) and (9) indicate that each vehicle in use should depart from the depot at the beginning of the distribution and return to the depot after finishing all the distribution tasks, respectively. Constraint (10) is the capacity constraint. Constraints (11), (12), and (13) provide the calculations of a i j k , l ik , and s ik . The time window constraint is specified as constraint (14). Finally, the distribution duration of goods for customer i for quality losses calculation is formulated as constraints (15) and (16). Constraint (17) shows that both decision variables can either be 0 or 1.

III. MPMOACS APPROACH
In this section, the proposed MPMOACS approach that is used to solve the CCL scheduling problem is described first. Then, as rescheduling is required when disruption events happen during the distribution process, a rescheduling strategy is introduced at the end of this section, i.e., in Section III-F.
The MPMOACS is an approach developed from the ACS and the MPMO framework. MPMO is a concise but effective framework for multiobjective optimization, which adopts M (M ≥ 2) populations for M objectives and one population focuses on one objective [38]. Besides, an information-sharing mechanism is designed to achieve the coevolution of multiple populations. Hence, three populations are required for our CCL model because there are three objectives. Herein, these three populations are denoted as three colonies since our approach is an ACS-based approach. The first colony for f 1 (i.e., quality losses) is named Qcolony, the second colony for f 2 (i.e., personnel and vehicle costs) is named Pcolony, and the third colony for f 3 (i.e., transportation costs) is named Tcolony.
The MPMOACS procedure is shown in Algorithm 1. First, in line 2, an archive named Archive is initialized as an empty set, which is used to store the Pareto solutions, and three pheromone matrixes for three colonies are constructed and initialized. Note that each colony constructs and updates its pheromone matrix individually in order to optimize its objective more effectively. Second, the solution construction for each ant of the three colonies is shown in lines 5 to 8. Moreover, after each ant (i.e., solution) is obtained, all its three objective fitness values are calculated in line 9. Third, the pheromone is updated locally in line 10. Fourth, in line 11, the solutions obtained from the three colonies are used to update the Archive through non-dominated sort. Fifth, the pheromone of each colony is updated globally in lines 12 to 14. Finally, Local update for pheromone matrix of colony; 11 Archive update; 12 foreach colony in ColonySet do 13 foreach solution in Archive do 14 Global update for pheromone matrix of colony; 15 Local search of Archive using RBLS;

End
as shown in line 15, the RBLS strategy is applied to improve the solutions in the Archive.

A. Solution Encoding
The encoding of a solution should specify the distribution route of each vehicle and the distribution sequence of each customer. Besides, the encoding format must be concise. Thus, a solution of MPMOACS is encoded as a two-dimensional array, and each row of the array represents the distribution route of a vehicle. That is, if m vehicles are used for distribution in a solution, the solution will be encoded as an array with m rows.
An example of a solution with 5 customers and 2 vehicles is shown in Fig. 2. Route 1 is the distribution route of vehicle 1, and Route 2 is the distribution route of vehicle 2. According to Fig. 2, Route 1 is 0→3→1→0, and Route 2 is 0→2→4→5→0. That is, vehicle 1 serves customer 3 first after departing from the depot, then it serves customer 1; and finally, it returns to the depot. Meanwhile, vehicle 2 serves customer 2 first, and then it serves customers 4 and 5 in order. Identically, vehicle 2 returns to the depot once the service of customer 5 is finished. Therefore, the encoding result of the solution sol in Fig. 2 is where the last element of each row represents the depot, which is denoted as n + 1. Herein, the last element is 6 because n is 5. The first row of sol is the distribution route of vehicle 1, and the second row is for vehicle 2.

B. Solution Construction
To construct a solution, an ant first initializes a solution as a two-dimensional array with an element 0 at the beginning of the first row, which means that the first vehicle is scheduled to depart from the depot. Then, the ant constructs a complete solution by iteratively determining the next customer to be served following the state transition rule. The state transition rule is formulated as where i is the current customer vertex visited by the current vehicle, and β is a parameter. J is the set of customers that have not been assigned to any vehicle until now and can be served by the current vehicle without constraint violation. If J is empty, another empty vehicle is put into use. In this case, the solution (i.e., the array) sets the next row to present that a new vehicle is used. The τ * (i, u) and η * (i, u) represent the pheromone information and heuristic information between the customers i and u, respectively. According to the MPMO framework, each colony has its pheromone matrix and the heuristic information tailored for its objective. To distinguish the pheromone matrix and heuristic information of a colony from the other colonies, a variable denoted as * is used. If the ant comes from Qcolony, * is Qcolony. Similarly, * is Pcolony when the ant comes from Pcolony, and * is Tcolony when the ant comes from Tcolony.
Moreover, a parameter q 0 (0 ≤ q 0 ≤ 1) is introduced to help determine the next customer to be served. If a random number q ranging from 0 to 1 is less than or equal to q 0 , customer u with the greatest value of τ * (i , u) × [η * (i , u)] β is chosen as the next customer. Otherwise, the next customer is determined through the roulette wheel selection, and the selection probability of customer j is calculated by

C. Heuristic Information Design
To determine Pareto solutions with high quality, different heuristic information should be designed for different colonies according to their objectives. For the ants of Qcolony, the formulation of heuristic information is First, to minimize the quality losses, the goods should be delivered to customers as soon as customers can accept the services. Hence, the ants of Qcolony should prefer to select customer j if j can be served as early as possible once j is ready to accept service. In this case, the first term of where the maximum value of h 1 ( j ) is 1, which occurs when s j k equals b j . The value of h 1 ( j ) decreases as the difference between s j k and b j increases. Meanwhile, considering that the deterioration rate is identical for goods of different customers, ants may prefer to serve customers with greater demand first because the distribution duration of more goods can be reduced, which is conducive to reducing the total quality losses of goods. Hence, the second term of η Qcolony (i , j ) can be represented as where the denominator term is the remaining capacity of vehicle k, and V k represents the set of customers that have been assigned to vehicle k. In this way, h 2 ( j ) is proportionate to q j . Nevertheless, if only h 1 ( j ) and h 2 ( j ) are considered, few feasible solutions can be found by Qcolony since the number of vehicles required by the ants exceeds the number of vehicles that are provided. This case occurs because the driving and waiting duration is ignored by h 1 ( j ) and h 2 ( j ), which may mislead ants to make vehicle k wait at the destination of a customer and assign other customers to new vehicles even though the waiting duration is sufficient for vehicle k to serve other customers. Therefore, the duration of driving and waiting should be considered to reduce the wastage of time so that the number of vehicles used can be reduced. Based on the analysis mentioned above, the last term of η Qcolony (i , j ) is formulated as which is the reciprocal of the duration of driving and waiting required for vehicle k from leaving customer i at l ik to start serving customer j . A greater value of h 3 ( j ) indicates less driving duration or less waiting duration. To minimize f 2 (i.e., the personnel and vehicle costs), the number of vehicles being used and the total distribution duration should be optimized. In other words, customers with less waiting or driving duration should be served first. Meanwhile, customers with higher urgency should have higher priority. The motivation of this heuristic information comes from the previous work [43] on the VRPTW. Thus, the formulation of heuristic information for f 2 is where the first term of (25) is the same as (24), and the second term is the measurement of urgency. As for f 3 , the objective of minimizing the transportation costs, the heuristic information should focus on the minimization of the driving distance and driving duration. One note is that both the distance and duration should be considered since the driving distance is fixed while the driving duration is time-dependent. Thus, the heuristic information of f 3 is simply determined as D. Pheromone Setting 1) Pheromone Initialization: The pheromone initialization of each colony is conducted according to where Num refers to the number of both customers and the depot, and sol greed y is an initial solution constructed through a greedy heuristic. In detail, sol greed y is constructed by iteratively choosing the customer nearest to the current position of the vehicle as the next customer to be served. If no customers can be served by the vehicle being used, a new empty vehicle is put into use. Note that although the τ 0 of the three colonies are calculated based on the same solution (i.e., sol greed y ), they are different since different objective fitness values of sol greed y are used.
2) Pheromone Update: Once an ant finishes the construction of a solution, the local update of the pheromone is carried out. For each edge traversed by the ant, the pheromone on the edge is updated according to where * represents Qcolony, Pcolony, or Tcolony, which depends on the colony to which the ant belongs. Parameter ρ is the degree of update of the pheromone.
The global update is executed when all ants in three colonies complete the solution construction. Since the CCL model is multiobjective, a set of solutions (i.e., the Pareto solutions) can be regarded as the global best solutions due to the nondominated relationship. Therefore, an archive with limited capacity, termed Archive, is introduced to store the Pareto solutions obtained by the three colonies during the search process. The update of Archive can be described as follows. In every generation, the solutions obtained by three colonies are added to Archive. Then, the solutions dominated by any other solutions are deleted from Archive. If the number of solutions in Archive still exceeds the maximum capacity of Archive, the crowding distance [24] of each solution is calculated. Then, solutions with the minimum crowding distance are removed from Archive until the number of solutions in Archive is reduced to the maximum capacity of Archive. Finally, the global update of pheromone is conducted on all solutions in Archive, which is shown as where sol is a solution in Archive. Parameter ε represents the pheromone enhancement degree.

E. Ranking-Based Local Search Strategy
The RBLS strategy is applied to the solutions in Archive, and the neighbor operator adopted in the RBLS is the best cost route crossover (BCRC) proposed in [44]. The BCRC is a crossover operator originally designed for the genetic algorithm, which can also be applied in the MPMOACS. The main idea of the BCRC can be described as follows. First, two routes, r 1 and r 2 , are randomly selected from two different parent solutions, p 1 and p 2 , respectively. Second, customers served in r 1 are removed from p 2 , and customers served in r 2 are removed from p 1 . Finally, customers in r 1 are reinserted into the most feasible position of p 2 in random order. A similar operator is conducted on r 2 and p 1 . Note that a position is considered feasible for a customer as long as the customer can be served without violating any constraints.
A schematic diagram of the BCRC is shown in Fig. 3. In the first step, the first route of p 1 is selected as r 1 , and the second route of p 2 is selected as r 2 . Then, customer 3 and customer 1 are removed from p 2 ; and customers 2, 3, and 5 are removed from p 1 , resulting in c 1 and c 2 shown in Fig. 3(b). Finally, the customers of r 1 and r 2 are reinserted into the most feasible positions of c 2 and c 1 , respectively. If there is no feasible position for a customer, a new empty vehicle is added.
As a problem-specific operator for the VRPTW, the BCRC measures the most feasible position of a customer through the total costs of routes. However, the truly optimal position of a customer cannot be determined directly because there are three objectives in our model. To address the measurement of the feasible degree of a position, ranking-based fitness (RF) is introduced. The position with the maximum RF is regarded as the most feasible position for a customer to insert. The major steps of the calculation of RF are as follows.
First, the ranking vectors of p 1 and p 2 , which are denoted as rk 1 and rk 2 , respectively, are obtained. Specifically, both rk 1 w io = rk io rk i1 + rk i2 + rk i3 (31) where w io is the weight depending on the ranking of objective o of p i . A higher ranking indicates a larger weight. It is designed in this way because we think a higher ranking may imply more room for improvement. However, due to the conflict among objectives, the improvement of an objective may result in recessions of other objectives. Herein, c i is regarded as better than p i if RF(c i ) is larger than 1. This is reasonable since at least one objective is improved when RF(c i ) exceeds 1. Finally, with the combination of the BCRC, the complete RBLS strategy is shown in Algorithm 2. To begin with, two different solutions p 1 and p 2 are selected randomly from the Archive, and then the BCRC operator is conducted on p 1 and p 2 , resulting in two new solutions (i.e., c 1 and c 2 ). As described in lines 8 to 12, a local search is considered successful only when both c 1 and c 2 are feasible and at least one of them is better than their parents (i.e., RF(c i ) > 1), and the variable cnt is used to record the number of failed local search operations. Once the number of failed local search operations reaches maxCnt or the size of Archive is less than 2,

Algorithm 3: Rescheduling Strategy
Input: an old distribution scheme, oldSol Output: the rescheduling distribution scheme, newSol 1 Begin 2 rnum ← the number of routes in oldSol; 3 for i = 1 : rnum do 4 r ← the i th route in oldSol; 5 ← the set of all undistributed orders in r ;

F. Rescheduling Strategy
During the distribution process, the time-dependent driving speed in the real traffic network may be changed (e.g., reduced) unexpectedly due to disruption events like traffic accidents and/or bad weather. In this case, the distribution scheme obtained based on the normal situation without disruption events may become not as efficient as expected, rescheduling of the remaining undistributed orders is thus required. Therefore, a rescheduling strategy is proposed to deal with the rescheduling problem.
Algorithm 3 describes the procedure of the rescheduling strategy. The main idea is to reschedule the routes of an old solution (i.e., the old distribution scheme) one by one based on the 2-opt [22] operator. As shown in lines 5 to 8, when rescheduling an old route, the distribution sequence of any two undistributed orders and all orders between them are reversed, leading to a new route. For example, for a route A→B→C→D→E→F with two selected undistributed orders as B and E, the resulting route is A→E→D→C→B→F. In this way, every two undistributed orders in an old route can result in a new route, and a series of new routes for an old route are thus obtained. Then, the best new route among all the new routes is used to replace the old route if it is better than the old route. When all the old routes have been rescheduled, the resulting new solution is adopted as the rescheduling result.

A. Dataset Generation
Since the CCL model in this paper is oriented to the real transportation environment, a CCL dataset that considers the real transportation environment, which is not provided in previous work, is needed. Therefore, we establish a multisource dataset based on real-captured transportation data of Guangzhou, China and Shenzhen, China.
First, to construct a test instance with n customers and 1 depot, n + 1 places are randomly selected to be the positions of customers and depot. Meanwhile, the road network that contains the driving paths between any two selected places is constructed. A schematic diagram of the road network based on the four selected places of Shenzhen, China is shown in Fig. 4. Second, the driving distance and time-dependent driving duration between any two places selected in the previous step are captured according to the road network. Note that the driving distance and driving duration are captured multiple times a day, and the capture operation is repeated for several days. Third, the average driving duration between any two places for each hour in a day can be obtained. Then, the average driving speed for each hour in a day that depends on both the driving distance and the average driving duration is calculated, resulting in the step function between the departure time and driving speed. In this case, the driving duration between any two places at any time can be calculated according to the step function.
Besides transportation data, order data and service data are indispensable for a complete test instance, which are generated randomly. To increase the diversity of instances, two distributions of time windows are adopted for each instance: a Gaussian distribution and a random distribution. In this way, 28 instances (14 for Guangzhou and 14 for Shenzhen) can be designed, and n is 30 and 40 for the instances of Guangzhou and Shenzhen, respectively. The 7 instances of Guangzhou with Gaussian distribution are denoted as GZ-G1, GZ-G2, …, and GZ-G7, while the 7 instances of Guangzhou with random distribution are denoted as GZ-R1, GZ-R2, …, and GZ-R7. Similarly, the 14 instances of Shenzhen are denoted as SZ-G1, SZ-G2, …, SZ-G7, SZ-R1, SZ-R2, …, and SZ-R7. The multisource dataset can be downloaded from [45].

1) Experimental Purposes:
The experiments are conducted to achieve three purposes. The first purpose is to investigate the importance of real transportation environment consideration, and the second purpose is to evaluate the performance of the MPMOACS on the multisource data-associated multiobjective CCL model. The last purpose is to investigate the effectiveness of the RBLS strategy and the rescheduling strategy.
To achieve the first purpose, four modes of the CCL model are designed according to the adoption of different transportation data, and the MPMOACS approach is applied to each mode individually. The four modes correspond to four real transportation environments with different accuracies. The first mode is the CCL model with linear distance and constant speed, denoted as the LDCS mode. The linear distance between two places is similar to the Euclidean distance but considers the spherical property of the earth. The constant speed is determined as 60km/h. The second mode is the CCL model with linear distance and r eal speed, and the real speed is the time-dependent speed. Herein, the second mode is denoted as the LDRS mode. The third mode named RDCS is the CCL model with r eal distance and constant speed, and the real distance is the driving distance. The last mode is the CCL model with r eal distance and r eal speed, that is, with real driving distance and time-dependent speed, which is also the mode adopted in our model. Similarly, the last mode is denoted as the RDRS mode.
The experiments of the second purpose can be divided into two parts. The first part is experiments conducted on the multisource dataset, and the MPMOACS is compared with the other six approaches, including the NSGA-II [24], MOEA/D [25], CCMO [21], TS-MOEA [22], the multi-objective local search (MOLS) approach [18], and the hybrid ant colony optimization (HACO) approach [20]. The NSGA-II and MOEA/D are the classic approaches for multiobjective optimization. CCMO is a coevolutionary framework developed recently and has been used for the optimization of the multiobjective VRP. The TS-MOEA is also a recently proposed optimization approach for the multiobjective VRP. The MOLS is proposed to solve the multiobjective VRP through the cooperation of multiple local search operations, and the HACO is proposed to solve the multiobjective VRP by integrating mutation operators into the ant colony optimization approach. Since the multisource data-associated multiobjective CCL model is also developed from the VRP, all approaches mentioned above can be adapted to solve this model. Therefore, it is reasonable to compare the MPMOACS with these approaches. The second part is the comparison experiments of the above approaches on six benchmark instances [43] for the TDVRP. The benchmark instances for the TDVRP are adopted due to the time-dependent speed used in the instances, which is similar to the real driving speed in the multisource data-associated CCL model. With the suggestion of [13], the working duration of the depot is evenly divided into 5 intervals, and the speed of each interval is set as [1.00, 1.60, 1.05, 1.60, 1.00] in sequence.
The third purpose is to investigate the effectiveness of the RBLS strategy and the rescheduling strategy. Herein, for the investigation of RBLS strategy, the compared approach named MPMOACS-noLS is the MPMOACS without the RBLS, which is designed by removing the RBLS from the MPMOACS.
2) Experimental Setting: The parameter configurations of all approaches are given as follows. Unless otherwise specified, parameter configurations of all the compared approaches follow the suggestions of their original papers since these parameter configurations have been well-tuned in their original papers to solve the problems that are related to our proposed multiobjective CCL model. For the MPMOACS, the size of each colony is set as 10. The parameter q 0 for state transition is 0.9, and β is 2. The values of both ρ and ε are set as 0.1. The maximum size of Archive is 300. The population size of the NSGA-II is 300, the crossover rate p c is 0.9, and the mutation rate p m is 1.0/n. As introduced in [25], the population size of the MOEA/D is 351 since there are three objectives. The size of the neighbor, T , is 20. The crossover rate p c is 1.0, and the mutation rate p m is 1.0/n. Regarding CCMO, the population size is 50, p c is 1.0, and p m is also 1.0/n. The TS-MOEA is a two-stage approach. The first stage is based on the EC-NSGA-II with a population of 100, and the second stage is developed from the MOEA/D with a population of 351 and T of 5. For MOLS, the population size is 300, the coefficient of tournament size is 0.5, and the population is updated in every iteration. In HACO, the population size is 35, the mutation rate is 0.6, q 0 is 0.6, ρ is 0.7, and the values of α, β, γ are 1, 4, 2 in order.
The execution time of all approaches on each instance in the multisource dataset is set as 200 seconds. For benchmark instances, the execution time is set as 300 seconds since the instance scales are larger. Considering the stochastic factors, each approach is conducted on each instance 10 times independently and the average results are compared. The experimental environment is a PC with an Intel i5-7400 CPU and 8.00 GB RAM.
3) Performance Metrics: Two performance metrics are employed for the comparisons among various modes or approaches.
The first metric is the hypervolume (HV) performance metric. The HV is a performance metric that can measure both the convergence and diversity of a mode/approach through the hypervolume constructed by the solutions obtained by the mode/approach and a reference point. A better mode/approach always has a larger HV. To calculate the HV, we first normalize the objective fitness values of all solutions obtained by all modes/approaches in all runs and set the reference point as [1.00, 1.00, 1.00]. Then, the HV is computed for each mode/approach in each run, and the average HV is adopted for comparison. Besides, the Wilcoxon rank-sum test at 0.05 significance level is used to compare the HV of RDRS mode/MPMOACS with the HV of the compared modes/approaches, and the comparison results are represented as '+', '-', or '≈', meaning that the performance of RDRS mode/MPMOACS is significantly superior, significantly inferior, or similar to that of the compared modes/approaches.
The second metric is C( A, B) [46], which is used to compare the dominant relationship between solutions of the MPMOACS and other approaches. For approaches A and B, C ( A, B) is defined as C(A, B) = |{b ∈ B|∃a ∈ A : a dominates or equals b}| |B| Approach A is regarded as better than approach B if C ( A, B) is greater than C(B, A), meaning that more solutions of B are dominated by or equal to solutions of A. Herein, the C (A, B) is calculated based on all non-dominated solutions found in all runs.

C. Importance of Real Transportation Environment Consideration
As introduced in Section IV-B, four modes of the CCL scheduling model, which are denoted as LDCS, LDRS, RDCS, and RDRS, are proposed to investigate the importance of considering the real transportation environment. Herein, the solutions obtained in the LDCS, LDRS, and RDCS are directly applied to the RDRS mode; and the average feasibility rates are shown in Table S.II in the supplementary material. The best results are shown in boldface for clarity.
It can be concluded from Table S.II in the supplementary material that the feasibility rates of all three modes are less than 40% and some of them are even 0, meaning that most solutions of traditional CCL models cannot be directly applied in the real transportation environment, which also indicates the importance of considering the real transportation environment. Besides, the LDRS dominates the other two modes on 19 of 28 instances, the RDCS surpasses the LDRS and LDCS on 8 instances, and the LDCS outperforms the other two modes on no instances. That is, both the LDRS and RDCS are better than the LDCS, and LDRS is even better than RDCS. This is because both LDRS and RDCS adopt more real information (i.e., real speed and real distance, respectively) than LDCS, leading to a more realistic transportation environment. In other words, the adoption of both real speed and real distance can improve the feasibility of scheduling results. Moreover, as LDRS is better than both LDCS and RDCS, the consideration of the real speed seems to be more important than that of the real distance. Table S.III in the supplementary material presents the HV comparison results of the LDCS, LDRS, RDCS, and RDRS. It should be noted that only feasible solutions of the LDCS, LDRS, and RDCS are considered in the comparison. Besides, the HV of an instance with no feasible solutions is denoted as N/A, and the best results are shown in boldface. As can be observed, the HV of the RDRS is the greatest among the four modes on all instances in the multisource dataset. That is, the hypervolume constructed by solutions of the RDRS is larger than those of the other modes on the multisource dataset. Furthermore, the RDRS dominates the other three modes on all instances in terms of the statistical test, which shows the superiority of the RDRS. Thus, the necessity and significance of considering the real transportation environment for CCL scheduling are verified.  boldface. Then, the Pareto solutions of each approach on 2 instances are provided for a more intuitive comparison. Table III provides the comparison results on the HV for all compared approaches and the MPMOACS. First, the HV of the MPMOACS is much better than that of both the NSGA-II and MOEA/D on all instances, which means that the MPMOACS is better than the NSGA-II and MOEA/D for multisource data-associated CCL model optimization on the multisource dataset. Moreover, the comparison results based on the statistical test also support the conclusion that the MPMOACS outperforms the NSGA-II and MOEA/D. This is reasonable since both the NSGA-II and MOEA/D are not problem-specific approaches designed for multisource data-associated CCL scheduling. Besides, the use of effective heuristic information also accounts for the better performance of the MPMOACS.
Comparison results of the CCMO, the TS-MOEA, and the MPMOACS on the multisource dataset also show the advantages of the MPMOACS. Specifically, the MPMOACS surpasses the CCMO and the TS-MOEA on all instances when the average HV is considered. Meanwhile, the last row of Table III shows that the MPMOACS significantly dominates the CCMO on all instances. Besides, the MPMOACS has significantly better performance on 26 instances when compared with the TS-MOEA while the TS-MOEA only obtains comparable results on 2 instances. Thus, it is reasonable to conclude that the MPMOACS has better performance on CCL scheduling when compared with the CCMO and the TS-MOEA.
Meanwhile, Table III shows that MPMOACS outperforms MOLS and HACO on all instances in the multisource dataset in terms of both average HV and statistical test, while MOLS and HACO cannot surpass MPMOACS on any instance.
The C( A, B) comparison results for all solutions obtained by each approach in 10 runs are shown in Table S Furthermore, the Pareto solutions of the MPMOACS and compared approaches on 2 instances are shown in Fig. 5, including one instance with Gaussian distribution and one instance with random distribution. The figure also shows that the performance of the MPMOACS is better.
In summary, the comparisons on two performance metrics (i.e., HV and C ( A, B)) and the presentation of Pareto solutions show that MPMOACS has overall better performance on the multisource dataset when compared with the NSGA-II, MOEA/D, CCMO, TS-MOEA, MOLS, and HCAO.
2) Experimental Results on the Benchmark Instances: Without loss of generality, experiments are conducted on six benchmark instances chosen from Solomon's benchmark [43], and the time-dependent speed was provided by Figliozzi [13]. The time-dependent speed of [13] is used because it is suggested by Figliozzi that the time-dependent speed can be directly applied to Solomon's benchmark, leading to a time-dependent benchmark of the VRPTW. The six instances chosen in our experiment are C101, C201, R101, R201, RC101, and RC201.
The experimental results on HV of the compared approaches and the MPMOACS are shown in Table IV, where N/A means that the approach cannot obtain any feasible solution within the maximum execution time. Besides, the best results are highlighted in boldface. The table shows that the performance of the MPMOACS is much better than those of the NSGA-II and MOEA/D. For example, the HV of NSGA-II and MOEA/D on RC201 are 0.233 and 0.057, respectively; while it is 0.729 for the MPMOACS. Furthermore, the MOEA/D and the NSGA-II cannot obtain feasible solutions on R101 within the execution time provided.
Different from the NSGA-II and MOEA/D, both the CCMO and the TS-MOEA can obtain feasible solutions on all benchmark instances. However, the MPMOACS is superior to the CCMO on all instances, which is similar to the comparison results on the multisource dataset. The statistical test results also show that the MPMOACS has significantly better performance than the CCMO on all benchmark instances. Besides, the MPMOACS also dominates the TS-MOEA on all instances.
The average HV of MOLS and HACO are significantly inferior to that of MPMOACS, meaning that the performance of MPMOACS surpasses those of MOLS and HACO on benchmark instances.
Similar to the multisource dataset, the performance metric C ( A, B) is applied on all Pareto solutions of the MPMOACS and other approaches to measuring the dominant relationship between the MPMOACS and compared approaches. The comparison results are shown in Table S The Pareto solutions of the MPMOACS and compared approaches on three randomly selected benchmark instances are plotted in Fig. S.1 in the supplementary material. The figure shows that the Pareto front of the MPMOACS is rather different from those of the compared approaches since the former shows better convergence and diversity. Thus, the comparison of the Pareto solution distribution of the MPMOACS and those of other approaches further shows that the MPMOACS has better performance on the optimization of benchmark instances.

E. Investigation of the RBLS Strategy
To investigate the effectiveness of the RBLS strategy, the MPMOACS-noLS approach is developed from the MPMOACS by removing the RBLS strategy. Comparison results of the MPMOACS-noLS and the MPMOACS are shown in Table S.VI in the supplementary material, with the best results presented in boldface. The MPMOACS is regarded better than the MPMOACS-noLS since the HV of the MPMOACS is greater than that of the MPMOACS-noLS on all instances. For example, the HV of the MPMOACS on GZ-G3 is 0.679, while it is 0.427 for the MPMOACS-noLS. Thus, with the assistance of the RBLS strategy, the performance of the MPMOACS can be improved.
The box plot for HV of the MPMOACS-noLS and the MPMOACS is given as Fig. S.2 in the supplementary material. Box plot is a kind of graph that can illustrate the convergence and stability of an approach simultaneously. The convergence of an approach is reflected through the height of a box, and the stability is demonstrated by the thickness of a box. An approach with a thicker box is considered less stable. As presented in Fig. S.2, boxes of the MPMOACS are much higher than those of the MPMOACS-noLS, which indicates that the MPMOACS has an advantage in terms of convergence when compared with the MPMOACS-noLS. Moreover, it can be observed from Fig. S.2 that boxes of the MPMOACS-noLS are thicker than those of the MPMOACS on most instances, meaning that the MPMOACS is more stable than the MPMOACS-noLS.

F. Investigation of the Rescheduling Strategy 1) Disruption Events Simulation:
To simulate the disruption events, two model parameters denoted as θ and sf are introduced.
Parameter θ is the number of edges that are affected by the disruption events. When the disruption events happen, θ untraveled edges from an old distribution scheme are selected randomly, and the driving speed of these edges is changed.
Parameter sf (0 < sf < 1) is the speed f actor that indicates the reduction degree of driving speed. In detail, the new driving speed is calculated by multiplying sf with the original driving speed. Meanwhile, for a distribution scheme, the time that the disruption events happen is randomly generated, and the change of the driving speed is assumed to keep until the end of the distribution process.
2) Objective Evaluation Criterion: As the change of driving speed due to the disruption events may cause the violation of the hard time windows constraints of the CCL model, the rescheduling is expected to minimize the total late delivery duration at first. The total late delivery duration of a new solution named sol is denoted by ld(sol), and is calculated as where s ik is the time that vehicle k begins to serve customer i , and e i is the upper bound of customer i 's time window. Besides, it is also preferred if the rescheduling strategy can obtain solutions with less objective fitness values (i.e., f 1 , f 2 , and f 3 ). Therefore, a two-level objective evaluation criterion is proposed to evaluate the new solutions obtained by the rescheduling strategy. Specifically, a new solution sol 1 is considered better than another new solution sol 2 if one of the following two conditions is satisfied: i) ld(sol 1 ) is less than ld(sol 2 ); ii) ld(sol 1 ) equals ld(sol 2 ) and the product of all the objective fitness values (i.e., f 1 , f 2 , and f 3 ) of sol 1 is less than that of sol 2 .
3) Results Discussion: Without loss of generality, we simulate the disruption events on seven typical selected instances (i.e., GZ-G1, GZ-G2, GZ-R4, SZ-G4, SZ-G5, SZ-R6, and SZ-R7). Moreover, three old solutions obtained by MPMOACS to each of these instances are randomly selected to conduct the rescheduling. The disruption event is simulated on each of the 7 × 3=21 solutions, with the number of edges that are affected by the disruption event (i.e., θ) set as 5 and the speed factor sf set as 0.2, 0.4, 0.6, and 0.8, resulting in totally 21 × 4 = 84 scenarios.
The total late delivery duration and the product of all the objective fitness values for both old solutions and new solutions in all the 84 scenarios are presented in Table S.VII in the supplementary material. Note that the new solutions with higher qualities (smaller total late delivery duration or/and smaller product of all the objective fitness values) are shown in boldface for clarity. The results show that the rescheduling strategy can obtain new solutions that greatly reduce the total late delivery duration or at least reduce the objective values in 42 scenarios. Moreover, the results also show that the rescheduling strategy can find new solutions with better performance in more scenarios when sf decreases. This indicates that the rescheduling strategy has more advantages when disruption events become more serious. Thus, the proposed rescheduling strategy is regarded as effective for the rescheduling in the CCL model.

V. CONCLUSION
In this paper, we first establish the multisource dataassociated multiobjective CCL model, which is a multiobjective optimization problem with the three objectives of minimizing the quality losses, personnel and vehicle costs, and transportation costs. By adopting real transportation data, our model can alleviate the problem that most solutions of traditional CCL models are inapplicable in reality. Then, the MPMOACS approach developed from the ACS and MPMO framework is proposed to optimize the model. Following the suggestion of MPMO, three ant colonies are employed in the MPMOACS. The first colony aims to optimize the first objective, and the second and third colonies are responsible for the optimization of the second and third objectives, respectively. Meanwhile, the RBLS strategy designed for multiobjective local search also improves the performance of MPMOACS. In addition, a rescheduling strategy is proposed to deal with the rescheduling caused by disruption events.
Comprehensive experiments on a multisource dataset and six benchmark instances are conducted to investigate the necessity of considering the real transportation environment and to evaluate the performance of the MPMOACS. As shown in the experimental results, the feasibility rates of modes without real transportation data are relatively low. Hence, the consideration of the real transportation environment in CCL scheduling is important. By comparing the MPMOACS with six state-of-the-art and very recent well-performing multiobjective approaches (i.e., NSGA-II, MOEA/D, CCMO, TS-MOEA, MOLS, and HACO), it can be concluded that MPMOACS is an effective approach for multisource dataassociated multiobjective CCL scheduling. Moreover, the employment of the RBLS strategy is shown to be helpful to the MPMOACS. The rescheduling strategy is seen as effective in dealing with the rescheduling problem.