Research on the Time-Dependent Split Delivery Green Vehicle Routing Problem for Fresh Agricultural Products with Multiple Time Windows

: Due to the diversity and the different distribution conditions of agricultural products, split delivery plays an important role in the last mile distribution of agricultural products distribution. The time-dependent split delivery green vehicle routing problem with multiple time windows (TDSDGVRPMTW) is studied by considering both economic cost and customer satisfaction. A calculation method for road travel time across time periods was designed. A satisfaction measure function based on a time window and a measure function of the economic cost was employed by considering time-varying vehicle speeds, fuel consumption, carbon emissions and customers’ time windows. The object of the TDSDGVRPMTW model is to minimize the sum of the economic cost and maximize average customer satisfaction. According to the characteristics of the model, a variable neighborhood search combined with a non-dominated sorting genetic algorithm II (VNS-NSGA-II) was designed. Finally, the experimental data show that the proposed approaches effectively reduce total distribution costs and promote energy conservation and customer satisfaction.


Introduction
Fresh e-commerce has become one of the main channels for fresh agricultural products, however, logistics distribution has always been the bottleneck of the development of fresh e-commerce.As a major consumer of fresh agricultural products, China has a huge and growing e-commerce business market in fresh agricultural products.According to the market report on e-commerce businesses of fresh agricultural products in China by iResearch, a consulting company, the trade volume of China's e-commerce market of fresh agricultural products will reach 458.5 billion yuan in 2020, an increase of 64% over 2019.The overall trade volume of the fresh food business in China is estimated to reach 11.1971 billion yuan by 2023 [1].The e-commerce business of fresh agricultural products has developed rapidly, but it is difficult to achieve decent profits due to the low implementation rate of cold-chain facilities, redundant circulation links and considerable cargo damage.Moreover, the "last mile" distribution problem has become the most prominent one, impeding the growth of e-commerce businesses in fresh agricultural products.According to statistics [1], the cargo loss rate in the distribution of fresh agricultural products is around 10%, and the cost of the distribution accounts for around 40% of the overall transportation costs.Furthermore, due to delayed delivery and high cargo loss rate, the complaint rate also remains high, and 90% of the complaints are about poor delivery service.Therefore, management decisions in the distribution process have become the key to the success of delivery service provided by e-commerce platforms of fresh agricultural products, as well as a fascinating subject of academic research.

of 28
The e-commerce operation mode of fresh agricultural products allows customers to select the agricultural products they need and submit orders through online platforms.These companies dispatch vehicles to deliver goods for each customer along the preplanned route from the distribution center after receiving orders from customers.Due to the perishability of fresh agricultural products, maintaining a vehicle's temperature and transport efficiency during delivery is crucial for successful delivery.Customers buying fresh agricultural products think highly of delivery efficiency, which means that the time it takes to deliver agricultural products to customers has a significant impact on customer satisfaction.For the distribution process of small batches and various categories, a variety of products can be delivered by a vehicle with multiple compartments of different temperature settings.According to Hsu and Chen's research [2], the distribution mode of multiple compartments is more appropriate to be applied to handle small but diverse customer demands for fresh agricultural products.Reed et al. [3] and Chen et al. [4] also studied the multi-compartment vehicle routing problem.In former models, the customer's demand is limited so that it does not exceed the capacity of a certain compartment.However, according to the real cases in many studies [5][6][7][8][9][10], it is possible for a customer's demand to exceed the capacity of a compartment.Therefore, this study focuses on exploring and determining whether the delivery with multi-compartment refrigerated vehicles serves the best option even when the customer's demand exceeds the capacity of a compartment.A model under the condition of split delivery is also developed to research the cold-chain distribution of agricultural products with large and diverse customer demands.
The split delivery vehicle routing problem (SDVRP), as a subtype of vehicle routing problems, has drawn considerable interest from researchers in recent years.The split delivery vehicle routing problem depicts a scenario in which many vehicles of the same weight depart from the same warehouse to deliver items to multiple customers, with each customer accepting the service of multiple vehicles [11].The benefits of split delivery in vehicle routing problems have been demonstrated [12].The SDVRP has been applied in many practical cases [13,14].This paper will study the problem of single compartment vehicle routing optimization based on split delivery, with an aim to provide a more rational distribution strategy and route optimization method for fresh food e-commerce companies, and to reduce the economic costs while improving customer satisfaction.Considering the urban traffic conditions and the need to deliver various products, a split delivery vehicle routing problem model with time-varying road network characteristics is developed, and a multi-objective optimization algorithm is designed to obtain the optimal solution set and screen out the suitable solution from the solution set.
The characteristics of diverse demand for fresh agricultural product delivery explored in this paper are quite similar to the SDVRP model.As such, this paper will expand on earlier research by applying the SDVRP to the problem of fresh agricultural product delivery.Since customers will be serviced more than once after an order is split, it is more practical to have multiple time windows for each customer.To make the model better simulate the actual conditions of reality, the road network congestion constraints, carbon emission and fresh agricultural product loss are considered in the calculation of distribution costs.Therefore, this paper studies a time-dependent split delivery green vehicle routing problem with multiple time windows (TDSDGVRPMTW) based on the aforementioned criteria and takes the total cost and customer satisfaction as optimization objectives.The NSGA-II algorithm framework [15] is improved in this study by combining it with other heuristic algorithms [4,15] to improve the model in this paper.However, the algorithm can only find the Pareto-optimal front but cannot provide a specific solution.In reality, decision-makers must follow a single, well-defined solution when making decisions.In order to solve this issue, the technique for order preference by similarity to the ideal solution (TOPSIS) approach is used to screen the Pareto-optimal front for the best solution of clear distribution that fulfills the actual needs.
Based on the study of the traditional vehicle routing problem, the split delivery vehicle routing problem (SDVRP) was first proposed by Dror and Trudeau [12].Generally, SDVRP research has been less comprehensive than studies on other VRPs.As with VRP, several SDVRP formulations have been presented while taking into account various restrictions [16][17][18].SDVRP with time windows (SDVRPTW) [19][20][21][22][23][24][25][26] is one of those that has gained increasing attention.Notably, this study investigated SDVRPTW in the distribution of fresh agricultural products and presents a heuristic algorithm for solving it.The following parts are a literature review of literature relevant to this paper.
(1) Split delivery vehicle routing problem with multiple time windows Split delivery vehicle routing problem with multiple time windows (SDVRPTW) is a variant of SDVRP that imposes restrictions on each customer by defining an interval at which a customer must begin to be served.Due to the added time window constraints, the model of SDVRPTW is more complex and more difficult to solve than SDVRP.There is not much former research on SDVRPTW.On the basis of SDVRP, Frizzell and Giffin [16] considered the time window constraint and first proposed SDVRPTW, a path construction method and two ways to improve the path.Ho and Haugland [19] proposed a tabu search-based solution to the SDVRPTW model and analyzed the experiment results from 100 customers.Desaulniers [20] proposed a new exact branch-and-price-and-cut method to solve the SDVRPTW, and the calculation results show that the method is effective.Based on the research of Desaulniers, Archetti et al. [21] proposed an improved strategy of the branch-and-price-and-cut algorithm, which further improved the performance of the algorithm to solve SDVRPTW.Salani and Vacca [22] presented a mixed-integer program for the SDVRPTW based on arc flow formulation, proposed a branch-and-price algorithm and applied an improved approach to the pricing and master problem.Based on a new relaxed compact model, Bianchessi and Irnich [23] proposed a new and tailored branch-and-cut algorithm to solve SDVRPTW.Computational experiments show that the new method can prove the optimality of several previously unsolved examples in the literature.Li et al. [26] assumed that service time is proportionate to demand, that customers have many time windows to choose from and can only be delivered in one of them.They came up with a three-indicator traffic flow model and a combined coverage model to solve this challenge.To overcome this problem, a branch-and-price-and-cut method is presented.
It is important to note, however, that the previous research on SDVRPTW focuses on the solution method of the abstract model and rarely involves the application of the model to specific fields.For example, there is no research that applies the SDVRPTW model to the distribution of fresh agricultural products.Customers often need a variety of fresh agricultural products that require different degrees of refrigeration.Thus, a split delivery of fresh agricultural products is reasonable in such a case.Fresh agricultural products are perishable and require cold-chain transportation, which will increase the complexity of the model and require more efficient solutions.As for the service time, most research adopted constant service time, and only Salani [22] and Li et al. [26] considered the variable service times that depend on the delivery quantity.Considering the service time depending on the delivery quantity represents a more pragmatic way of thinking.Previous research overlooks traffic congestion in road networks, which makes it difficult to apply the research results to reality.This study takes into account variable service time and adds time-varying road network constraints to the SDVRPTW model, so as to make the model fully represent reality.
The VNS-NSGA-II algorithm proposed in this paper belongs to the multi-point heuristic algorithm.Compared with many works of literature [20][21][22][23][24] that use precise algorithms, the method in this paper can solve large-scale problems.Compared with the research from McNabb et al. [25], the method in this study includes an adaptive genetic evolution process and is easier to jump out of the local optimum in the process of searching for a solution.Both the model of this paper and the model of Li et al. [26] set multiple time windows, however, their model only allows customers to select one window to receive products, and the model of this paper allows customers to receive products at each time window.Most research on SDVRP focused on the optimization of a single objective, and only a few of them focused on multi-objective optimization [27][28][29].This paper focuses on the application of fresh agricultural products distribution, which requires a high level of delivery efficiency.In addition to transportation costs, customer satisfaction is also a critical aspect.Therefore, establishing a multi-objective optimization model with two distinct objectives of distribution cost and customer satisfaction as optimization objectives is a reasonable choice.
(2) The Distribution of Fresh Agricultural Products According to previous research [8][9][10], customers' demand for the distribution of fresh agricultural products mostly occurs in morning rush hours, during which the speed of vehicles is very different from other time periods.Therefore, it is a reasonable choice to regard the vehicle speed as a time-dependent function.At the same time, the amount of carbon emissions produced would change in proportion to the speed of the vehicle [30], hence considering carbon emissions in the model can make the final delivery solution more environmentally friendly.Due to the requirement of high delivery efficiency in the distribution of fresh agricultural products, whether the goods can be delivered in time will have a huge impact on customer satisfaction, which is also an important indicator to measure the service level of enterprises.Therefore, we should take customer satisfaction as an optimization objective and strive to improve it.In order to make the model better, we simulated actual conditions of fresh agricultural product distribution.
(3) Time-Dependent Vehicle Routing Problem The time-dependent vehicle routing problem is another derivative of VRP.To simulate road traffic congestion, Malandraki and Daskin [31] proposed TDVRP.They treat the vehicle travel time at any two points as a step function of the vehicle departure time.Malandraki and Dial [32] proposed a restricted dynamic programming heuristic for solving the time-dependent traveling salesman problem.However, in the above research, there are situations where vehicles that depart later arrive first, which do not meet the First-In-First-Out (FIFO) criterion.Ichoua et al. [33] proposed a new time-varying model in which the vehicle's travel speed is a step function, and the corresponding travel time is a piecewise linear function.Another time-varying model is proposed by Fleishmann et al. [34], which is based on a smooth step function of travel time.The above two methods confirm the principle of FIFO, which better simulates the actual conditions of the real world.
What is more, based on the research in the previous work and the innovations in this study are summarized as follows: 1.
The SDVRPTW model has been applied rarely in the research of fresh agricultural product distribution.The basic SDVRPTW model normally takes into account load constraints, split delivery and the precondition that the consumer has just one time window.In order to make the model better simulate the actual conditions of fresh agricultural product distribution, this study will consider and evaluate a time-varying road network, carbon emissions and customer satisfaction on the basis of SDVRPTW to develop the TDSDGVRPMTW model.

2.
While multi-objective optimization problems are frequently aimed at obtaining a Pareto-optimal front, decision-makers expect a complete and feasible solution.Furthermore, the TOPSIS method is employed to select the solution that satisfies the requirements from the Pareto-optimal front.3.
It is verified by a real-world case that as a delivery strategy, the TDSDGVRPMTW model proposed in this paper not only can effectively reduce the total cost of fresh agricultural products distribution, but also improve customer satisfaction.
The description of the problem and the hypotheses are as follows.
The TDSDVRPMTW proposed in this paper can be described as follows: a group of coldchain vehicles depart from the distribution center, visit customer nodes in the order indicated in the distribution scheme, and then return to the distribution center.K = {1, 2, • • • , k} is the set of vehicles.N = {0, 1, 2, • • • , n} is the set of all nodes in the distribution network where {0} represents the distribution center, N = N\{0} denotes the customers.The distance traveled by vehicle k between node i and node j is D ij , and the travel of the vehicle incurs a travel cost TC.Fixed cost FC is incurred when the vehicle is used.Service costs SC are incurred when the vehicle serves the customer.The speed of vehicles is related to time.Vehicles travel at varying speeds during various time periods.H = {1, 2, • • • , h} is the set of all time periods.Carbon dioxide CC is emitted into the atmosphere when a vehicle is driven.Fresh agricultural products are divided into several types based on their respective temperature requirements for transportation.W = {1, 2, • • • , w} is the set of all types of fresh agricultural products.The refrigeration cost RC of fresh product w is e w .Vehicle k is only allowed to deliver one type of product w.Each customer requires a variety of fresh agricultural products, and each customer has a number of service acceptance time windows equal to the number of fresh product varieties required by the customer.Additionally, it is allowed that a single customer can be served by multiple vehicles.The service time s ik is proportional to the quantity of products delivered.The customer satisfaction ACS i of customer i is determined by the remaining delivery time after the customer received goods.Total costs include travel costs TC, fixed costs FC, service costs SC, refrigeration costs RC and carbon emissions costs CC.The purpose of this study is to produce a delivery solution that has a lower total cost and a higher level of customer satisfaction.Table 1 lists all notations used in the proposed model.In addition, the following basic assumptions are made in this paper: 1.
The customer's demand for fresh agricultural products is split and distributed according to the temperature required for distribution, and each customer has multiple time windows for receiving services.

2.
For each service, every customer will complete a satisfaction rating, the value of satisfaction evaluation depends on the deviation degree between the remaining time when the vehicle leaves the customer and the time window.

3.
Only one type of produce can be delivered by a single vehicle.A vehicle can only serve a customer once.4.
The day is divided into several time periods, and vehicles travel at various speeds at different periods.
The rest part of this paper is organized as follows.In Section 2, the mathematical formulation of the time-dependent green vehicle routing problem with multiple time windows (TDSDGVRPMTW) is given.Then, the algorithm for obtaining the Pareto-optimal front based on variable neighborhood search combined with non-dominated sorting genetic algorithm II (VNS-NSGA-II) is presented.Besides, the technique for order preference by similarity to the ideal solution method (TOPSIS) is introduced, and the process of experiments is shown in the final part of the section.In Section 3, the results of the numerical experiments are analyzed and compared with the relevant literature, and also the reasons for the similarities and differences are discussed.Finally, conclusions and future research directions are given in Section 4. The average level of satisfaction for customer i Q K Maximum number of vehicles that can be used Distance between nodes i and j in the logistics distribution network Q i The total demand of customer i for various fresh agricultural products

L ik
The time when vehicle k completes its service and leaves customer i. q iw The demand of customer ifor fresh agricultural product w

T ik
The time when the vehicle k arrives at the customer i ET i The earliest time for node i to receive service The time when the vehicle k starts to serve customer i and also when the product is received by customer i

LT i
The latest time for node i to receive service The time that vehicle k serves customer i t ijk The driving time of vehicle k on sec tion (i, j) The driving time of vehicle k on road sec tion (i, j) in time period h The distance of vehicle k on road sec tion (i, j) in time period h st i The total number of services received by customer i.

X ijk
Equal to 1 if the vehicle k runs on the road sec tion (i, j), and 0 otherwise Equal to 1 if the vehicle k runs on the road sec tion (i, j) in time period h, and 0 otherwise y ik Equal to 1 if vehicle k visits client i, and 0 otherwise z k Equal to 1 if vehicle k is used, and 0 otherwise

TDSDGVRPMTW Model
Based on the needs of building the model, this paper uses the corresponding symbols which are listed in Table 1.
Through the comprehensive analysis above, the multi-objective optimization model of TDSDVRPMTW is given by the following, this model objectives include minimize total cost F 1 and maxmize average customer satisfaction F 2 . Minimize Maximize Subject to: Equation ( 1) is a function to minimize the total cost.It includes the vehicle travel cost (TC), fixed cost (FC), customer service cost (SC), refrigeration cost (RC) and carbon emission cost (CC).Equation (2) measures the average customer satisfaction (ACS).ACS is the proportion of fully satisfied customers to the total number of customers.Constraint (3) expresses the carrying weight limit of each vehicle.The combined demand of all customers on each vehicle's delivery route cannot exceed the maximum vehicle load limit.Constraint (4) ensures that the demand of each customer is equal to the sum of the customer's demand for each type of fresh agricultural product.Constraint (5) ensures that each customer can only be served once by each vehicle.Constraint (6) ensures that each customer must be served and the number of times served does not exceed the total number of vehicles.Constraint (7) ensures that the number of vehicles starting and arriving at each node should be balanced.That is to say, when a vehicle arrives at a customer, it will inevitably leave the customer.Constraint (8) expresses the relationship between road section driving time and vehicle driving time within a period.In other words, the total time spent on a road section is equal to the sum of the time spent on that road section in each time period.
Constraint (9) expresses the relationship between road section distance and vehicle driving distance in each time period.That is to say, the distance traveled by the vehicle on the road segment is equal to the sum of the distance traveled by the vehicle on the road section in each time period.Constraint (10) indicates that the vehicles returning to the distribution center should be restricted to the constraint of the working time window of the distribution center.Constraint (11) indicates the relationship between service start time and arrival time.Constraint (12) indicates the relationship between service start time, service duration and departure time.Constraint (13) ensures the time the customer receives the fresh products is equal to the time the service starts.Constraint (14) indicates that the number of times service received by each customer is equal to the total number of vehicles serving this customer.Constraint (15) states the binary decision variable.

•
The travel cost The travel cost (TC) refers to the variable cost incurred by each vehicle for delivery activities, mainly consisting of fuel consumption, repairs and driver pay rate.For simplicity, only the effect of distance on travel cost is considered.The cost of the vehicle's travel is generally proportional to the distance and can be calculated as follows: where X ijk is a binary variable, X ijk = 1 when vehicle k runs through the road section (i, j), otherwise X ijk = 0. D ij is the distance between node i and j, ϕ is the vehicle's travel cost per unit distance.

•
The fixed cost Vehicles' fixed costs are typically constant.It is unrelated to customer demand or delivery distance.It mostly consists of rent or use loss, and other labor costs.Calculate the cost of the vehicle's fixed (FC) use by: where z k is a binary variable, when vehicle k is used, z k = 1; otherwise, z k = 0. δ is the fixed cost per vehicle.

•
The service cost Since each customer can be serviced by multiple vehicles, the cost impact of the number of times each customer is serviced needs to be considered.The service cost (SC) is proportional to the frequency with which each customer is served.The service cost can be calculated as follows: where st i is the number of times the consumer i receives services and κ is the cost of each service.

•
The refrigeration cost Refrigeration cost (RC) refers to the energy cost of refrigeration.According to the study of Wang et al. [35], this paper treats refrigeration costs as a time-dependent function, assigns a fixed unit time cost to each temperature, and considers the refrigeration cost of the unloading service.The cost of refrigeration can be calculated as follows: X ijk e w t ijk + s ik (19) where X ijk is a binary variable, X ijk = 1 when vehicle k runs through the road section (i, j), otherwise X ijk = 0. e w is the refrigeration cost per unit time of fresh agricultural product w.t ijk is the driving time of vehicle k on section (i, j).s ik is the time that vehicle k serves customer i.

•
The carbon cost This paper uses the MEET model in reference [30] to explore the calculation of carbon emissions.The carbon emission estimation function for freight vehicles weighing between 3.5 and 40 tonnes is v .ϕ is the ratio of the vehicle's actual load to its capacity, ν is the vehicle speed (km/h), and parameters a and b are the load correction coefficient, whose value is related to the vehicle's weight range.If vehicle k's driving distance in time period h is d h ijk (km), the carbon emissions E h ijk (kg) can be computed as follows: assigns a fixed unit time cost to each temperature, and considers the refrigeration cost of the unloading service.The cost of refrigeration can be calculated as follows: ( ) where  is a binary variable,  = 1 when vehicle  runs through the road section ,  , otherwise  = 0.  is the refrigeration cost per unit time of fresh agricultural product . is the driving time of vehicle  on section ,  . is the time that vehicle  serves customer .

 The carbon cost
This paper uses the MEET model in reference [30] to explore the calculation of carbon emissions.The carbon emission estimation function for freight vehicles weighing between 3.5 and 40 tonnes is   =  +   +   +   + + + , while the load correction function is  ,  =  +   +   +   +   +   +   + . is the ratio of the vehicle's actual load to its capacity,  is the vehicle speed (km/h), and parameters  and  are the load correction coefficient, whose value is related to the vehicle's weight range.If vehicle 's driving distance in time period ℎ is  (km), the carbon emissions  (kg) can be computed as follows: = ⋅ ⋅ 1000 The cost of carbon emissions (CC) is mainly generated by the energy consumed by vehicles while on the road.The carbon emission calculation function adopted in this paper is related to the vehicle speed, travel time and load capacity.The cost of carbon emissions can be calculated as follows: where  is a binary variable,  = 1 when the vehicle  runs on the road section ,  in time period ℎ, otherwise  = 0.  is the carbon emission generated by vehicle  driving on the road section ,  in time period ℎ.  is the unit price of carbon emissions.

Customer Satisfaction Measurement Method
Customer satisfaction is an important indicator used to gauge the service quality of firms in the study on the vehicle routing problem for fresh agricultural products.The measurement methods of customer satisfaction can be divided into three categories: measure satisfaction by the freshness of products, by product delivery time, and by both freshness and delivery time.Wang et al. [36] used the freshness of perishable products to measure customer satisfaction, and proposed a multi-objective VRP optimization model with mixed time windows and perishability assessment to minimize transportation costs and maximize the freshness of perishable products.Wang et al. [37] established a customer satisfaction evaluation model, in which both the timeliness of the distribution of fresh agricultural products and the loss of freshness of agricultural products are considered.We adopts the satisfaction function based on the time window.Because this paper considers that each customer has many time windows and may receive vehicle services throughout each window, the satisfaction of each customer is averaged by the number of vehicle visits.Customer satisfaction is determined by the time the vehicle completes its servicing.The satisfaction function curve for customer  in the th time window is depicted in Figure 1.In Figure 1, the horizontal coordinate represents time, and the vertical coordinate represents customer satisfaction.When the vehicle completes the service in  ,  , and customer satisfaction is 100, this is the best time for service.The service can begin The cost of carbon emissions (CC) is mainly generated by the energy consumed by vehicles while on the road.The carbon emission calculation function adopted in this paper is related to the vehicle speed, travel time and load capacity.The cost of carbon emissions can be calculated as follows: where x h ijk is a binary variable, x h ijk = 1 when the vehicle k runs on the road section (i, j) in time period h, otherwise x h ijk = 0. E h ijk is the carbon emission generated by vehicle k driving on the road section (i, j) in time period h.µ is the unit price of carbon emissions.

Customer Satisfaction Measurement Method
Customer satisfaction is an important indicator used to gauge the service quality of firms in the study on the vehicle routing problem for fresh agricultural products.The measurement methods of customer satisfaction can be divided into three categories: measure satisfaction by the freshness of products, by product delivery time, and by both freshness and delivery time.Wang et al. [36] used the freshness of perishable products to measure customer satisfaction, and proposed a multi-objective VRP optimization model with mixed time windows and perishability assessment to minimize transportation costs and maximize the freshness of perishable products.Wang et al. [37] established a customer satisfaction evaluation model, in which both the timeliness of the distribution of fresh agricultural products and the loss of freshness of agricultural products are considered.We adopts the satisfaction function based on the time window.Because this paper considers that each customer has many time windows and may receive vehicle services throughout each window, the satisfaction of each customer is averaged by the number of vehicle visits.Customer satisfaction is determined by the time the vehicle completes its servicing.The satisfaction function curve for customer i in the gth time window is depicted in Figure 1.In Figure 1, the horizontal coordinate represents time, and the vertical coordinate represents customer satisfaction.When the vehicle completes the service in ET ig , LT ig , and customer satisfaction is 100, this is the best time for service.The service can begin when the vehicle arrives at the customer within EET ig , ET ig or LT ig , ELT ig .On the other hand, customer satisfaction is low when the service is completed within this time limit, and satisfaction decreases as the deviation from the optimal service time period increases.If the vehicle provides service to the customer prior to EET ig or after ELT ig , the customer's satisfaction level is 0. EET ig = ET ig − θs i , ELT ig = LT ig + θs i are the boundaries of the customer's tolerable time points.Where θ is the customer tolerance coefficient, s i is the customer i's service time, and L ik is the time when vehicle k completes its service and departs for customer i.G denotes the set of all of customer i's time windows, and st i denotes the number of times customer i gets served.Then, the average level of satisfaction for each customer (ACS) is as follows: where the satisfaction function is: This paper adopts average customer satisfaction to measure the overall satisfaction level of the distribution schemes.The average customer satisfaction function can be expressed as Equation (2).

Time-Dependent Vehicle Speed Calculation Method
Based on existing approaches [38,39], this paper proposed a method of time division for calculating travel time.The distribution center's opening hours are divided into many time periods, and the vehicle speed varies according to the time period.Let F be the length of the time period; H = {0, 1, 2, • • • , h} is a set of all time periods, [h, h + 1] said the number of h time period.d h ijk , t h ijk , g h ijk denote the distance, time and speed of vehicle k on the road section (i, j) in time period h, respectively, and D ij denotes the road section's distance (i, j).D h ij is the distance traveled by vehicle k to complete the remaining distance (i, j) after time h; L ik is the time when vehicle k departs from customer i; and h k is the remaining drivable time of vehicle k in time period h.As a result, the following steps are used to calculate the driving time t ijk of vehicle k on road section (i, j): Step 1: Determine how long the initial phase will last.
Step 2: Step 2 should be repeated; otherwise, . The section (i, j) driving time computation is finished.

VNS-NSGA-II Algorithm
A variable neighborhood search combined with the non-dominated sorting genetic algorithm II (VNS-NSGA-II) was designed.NSGA-II has been extensively applied in research on VRP-related problems as a multi-objective combinatorial optimization algorithm [38][39][40][41].The NSGA-II algorithm decreases the complexity of the non-dominated sorting genetic algorithm and has the advantages of rapid execution and good solution set convergence.A number of heuristic strategies [4,15] are introduced into the NSGA-II algorithm in this study to boost search efficiency and avoid local optimum.Adaptive functions [42] are introduced to the crossover and mutation processes to dynamically alter the likelihood of crossover and mutation, and the variable neighborhood search algorithm (VNS) [4] is added to conduct a variable neighborhood search for good individuals in the population.It has the potential to improve the algorithm's local search capacity.Based on the advantages of the preceding techniques, the VNS-NSGA-II algorithm is constructed to solve the TDSDGVRPMTW model in this study.The VNS algorithm is divided into three stages: initialization, genetic evolution and variable neighborhood search.The initialization step generates the initial population in a random manner.In the process of genetic evolution, there exist crossover and mutation operators whose likelihood of execution is dynamically governed by adaptive functions.In the process of the variable neighborhood search, there are three types of neighborhood search operators.Non-dominated sorting and crowding distance calculation are at the end of the variable neighborhood search process.Figure 2 shows the detailed process of the algorithm.

Population Initialization
The chromosome in this study uses mixed natural number and letter coding, with natural numbers ranging from 1 to n representing the customer node and capital letters A, B, C, and D representing various types of fresh agricultural products.Each vehicle delivers the same type of fresh agricultural product, and each chromosome holds the delivery schedules of all vehicles and therefore contains each customer's fresh produce needs.A chromosome is depicted in Figure 3 as an example.Where "A1" means that customer 1 needs type A fresh agricultural products.The chromosome can be divided into three segments according to the type of fresh products.Each segment is delivered by a different vehicle.When decoding chromosomes into vehicle routes, the chromosomes are first divided into segments based on the type of fresh agricultural product, and then the chromosome segments with the same type of fresh agricultural product are divided into route segments based on the vehicle load capacity and the distribution center's operating time.This study employs a random method to generate the initial population.In other words, the customer numbers are randomly arranged according to the type of fresh agricultural products that suit customers' demands.

•
Selection Operator This paper combines the optimal protection strategy and roulette selection method to select chromosomes.The specific steps of the optimal protection strategy are to find the two chromosomes with the highest fitness and the lowest fitness in the current population; the fitness value of the chromosome with the highest fitness is compared with the highest fitness value of each generation in history.If the current value is higher, it will be regarded as the chromosome with the best protection; otherwise, the best protection object remains unchanged and remains the best chromosome in history.The chromosome with the worst fitness in the current population is replaced with the one with the best protection.The probability of each chromosome being selected in the roulette method is p n = f n / ∑ f n n , and the greater the fitness of the chromosome, the greater the probability of being selected for the cross-mutation operation.

• Crossover Operator
Because chromosomes contain vehicle routes for delivering various types of agricultural products, all vehicle route segments in chromosomes are classified prior to the crossover to ensure that the crossover occurs only between vehicle route segments of the same type.By randomly generating two crossover points on parent chromosomes X and Y and separating the two parent chromosomes into three pieces, we improved the crossover approach in this paper.The center sections of chromosomes X and Y were excised and inserted into the front and back sections of offspring chromosomes Y 1 and X 1 , respectively.The remaining front and back sections of parent chromosomes X and Y were placed in the same order into chromosome Y 1 's back section and chromosome X 1 's front section.In the two offspring chromosomes, leave the chromosome segment in the two crossover sections alone and delete the duplicated chromosome segment in the remaining places.The advantage of this crossover approach is that two identical parent chromosomes produce two different offspring.Figure 4a depicts the unique crossover procedure.In order to make the crossover operator work more efficiently, this study uses the improved adaptive adjustment approach presented in reference [42], which takes into account the number of iterations, the fitness values of chromosomes and populations, and the number of unmodified chromosomes in each generation population, as indicated in Equation (24).(24) where p c stands for adaptive crossover probability, p c1 and p c2 stand for adaptive adjustment parameters, and p c1 > p c2 .f l stands for the fitness value of individuals with high fitness in the chromosomes to be crossed; f avg stands for the average fitness value of each generation population, and f max stands for the maximum fitness value of each generation population.The popsize denotes the population size, whereas gen represents the current iteration number, M represents the maximum iteration number, and U represents the number of individuals with unchanged chromosomes.

• Mutation Operator
Before mutation, route segments in chromosomes should also be classified to ensure that mutation occurs between vehicle routes of the same type.The mutation method of the random chromosome point exchange is used in this study.The specific processes are as follows: initially, the chromosomes to be mutated are chosen, and then two random mutation points on the chromosomes are chosen in a random manner.Swap two points to form a new chromosome.Figure 4b depicts the mutation process.Similar to the crossover process, this study adopts an improved adaptive adjustment function to determine the probability of the mutation operator's operation and to improve the efficiency of the mutation operator [42].Equation ( 25) is the adaptive function of mutation probability: (25) where p m represents the adaptive mutation probability, p m1 and p m2 are adaptive adjustment parameters and p m1 > p m2 , f is the fitness value of chromosomes to be mutated.

Variable Neighborhood Search Operators
In order to further improve the quality of chromosomes in the population, we added a neighborhood search operator to the algorithm to perform a deep search for some excellent solutions.The excellent chromosomes in the population are the operation object of variable neighborhood search, and the chromosomes in the population are sorted according to fitness from high to low, with the chromosomes in the top half as the excellent chromosomes.In each neighborhood search operator, a node is chosen in a random manner, the distance between it and all other nodes is calculated, and the remaining nodes are organized in ascending order to generate a list of distance values.The variable neighborhood search operator used in this paper is similar to the neighborhood structure proposed by Sánchez et al. [43].The following are the variable neighborhood search operators proposed in this paper. •

2-opt operator
Node i is randomly selected from the chromosome segment.Select the first node from the list of distance values for node i as node j.If the first node in the distance value list does not exist in the current chromosome segment, select the next node in turn.The 2-opt operator disconnects node i from the node behind it and node j from the node behind it and reconnects node i with node j.This operation will be kept if the fitness improves, otherwise, the next node in the list of distance values is tried, and the process repeats until a better chromosome is found or the maximum number of searches is reached.The process of the 2-opt operator is shown in Figure 5a.

• Single node move operator
Node i is randomly selected from the chromosome segment.Select the first node in the list of distance values for node i as node j.If the first node does not exist in the current chromosome segment, the next node in the list is selected as node j.The single node move operator removes node j from its original position and inserts it behind node i, making node i adjacent to node j.After insertion, the chromosomal fitness value is calculated.The operation is kept if the fitness was increased, otherwise, the next node in the list of distance values is tried, and the process repeats until a better chromosome is found or the maximum number of searches is reached.The process of the single node move operator is shown in Figure 5b.

• Double nodes move operator
Two adjacent chromosome nodes are randomly selected.Node i is the first of the two nodes.Select the first node in the list of distance values for node i as node j.If the first node does not exist in the current chromosome segment, the next node in the list is selected as node j.Insert two adjacent nodes containing node i after node j so that node i and node j are adjacent.The operation was kept if the fitness was increased, otherwise, the next node in the list of distance values is tried, and the process repeats until a better chromosome is found or the maximum number of searches is reached.The process of the double nodes move operator is shown in Figure 5c.

Non-Dominated Sorting and Crowding Distance Calculation
The objective function values of all chromosomes in the population are calculated, and the crowding distance is calculated.Finally, the Pareto-optimal front is obtained.Let F max m and F min m be the maximum and minimum values of the mth objective function, respectively, while F i−1 m and F i+1 m are the mth objective function values of the two solutions adjacent to the ith solution.Then, the crowding distance CD i of the ith solution can be calculated by Equation (26).
The new population is selected according to the Pareto-optimal front and crowding distance.Then, it is judged whether the maximum number of iterations is reached.If the maximum number of iterations is not reached, return to the genetic evolution process and continue to iterate.When the maximum number of iterations is reached, the algorithm is terminated and the Pareto-optimal front is printed.The Pareto-optimal front can be regarded and thought of as a solution set, including numerous potential delivery schemes.To determine which of these delivery schemes best fits the needs, the TOPSIS approach is employed in this study.

Select the Optimal Solution Strategy
When the Pareto-optimal front is created using the VNS-NSGA-II algorithm, we used the TOPSIS method to analyze all solutions and select the optimal strategy.The following are the specific steps: Step 1: Complete the index homogenization process.The objective function F 2 is adjusted to F 2 = 100 − F 2 in this study so that it can be minimized alongside F 1 .
Step 2: Construct the original data matrix.Assume that there are n solutions and m objective functions in the solution set, then matrix B is shown in Equation (27).
Step 3: Vector normalization of indicators: Obtain the normalized matrix Z.
Step 4: Determine the optimal scheme Z + and the worst scheme Z − for each indicator: Step 5: Calculate the proximity of each solution to the optimal scheme and the worst scheme: where ι j is the weight of the jth indicator, which is determined based on actual need.
Step 6: Calculate the closeness of each solution to the optimal solution.
where 0 ≤ C i ≤ 1, the closer C i is to 0, the better the evaluated solution.After selecting the solution that best fits the requirements using the TOPSIS technique, the decision-making process is complete.

Validation of the Simulation Model
In order to verify the validity of the method proposed in this study, we conducted three experiments, respectively is algorithm comparison experiment, solution selection experiment and real case experiment.The algorithm comparison experiment and solution selection experiment used R201 dataset from Solomon [44] benchmark.The customer demands in the dataset are randomly divided into several parts, whose numbers are limited up to four.The time window is also randomly split into several periods.Python 3.8 programming is used to carry out the experiment.According to the work of Fan et al. [42], the algorithm's parameter settings are connected to the size of the dataset utilized in the experiment as follows: p c1 = 0.7, p c2 = 0.5, p m1 = 0.01, p m2 = 0.008, maxit = 100 ∼ 300 iterations, population size popsize = 100 ∼ 200, and maximum field search times S t = 15 ∼ 30.The related parameters for calculating the vehicle travel time and carbon emissions are identical to those in the work of Liu et al. [45], which are shown as follows: the distribution center's time 0 is 7:00 a.m., the traffic congestion periods are 8:00~9:00, 18:00~19:00, and the vehicle speed is 20 km per hour.For the time period h, according to the remainder function η = h mod 3, η is (1, 2, 0) corresponding to (54, 72, 42) km/h, respectively, with three time-varying velocities.The correlation coefficients for the carbon emission model are as follows: 33.The carbon emission price µ = 0.0528 yuan per kilogram, the time window correlation tolerance coefficient is θ = 0.5, the driving cost per kilometer is ϕ = 1.3, the fixed usage fee is set at δ = 20, and the single service cost is set at κ = 3.All the experiments were carried out ten times, with the best result being chosen.The results of all experiments will be analyzed in detail in the Section 3.

Comparison with Other Efficient Algorithms
This paper compares the VNS-NSGA-II algorithm with other algorithms designed to solve multi-objective VRP problems.The VNS-NSGA-II algorithm is also compared to other algorithms designed to solve multi-objective VRP problems.As the VNS-NSGA-II algorithm is improved on the basis of the NSGA-II algorithm, the NSGA-II algorithm [38] is selected for comparison in order to verify the improvement.Three neighborhood search operators are added to the VNS-NSGA-II algorithm, while the many-objective gradient evolution (MOGE) algorithm [46] also has three search operators.In order to test the search ability of the variable neighborhood operators, the MOGE algorithm is selected for comparison.
Both the MOGE algorithm and NSGA-II algorithm are designed to solve multiobjective VRP.These two algorithms can be applied to the model in this paper.There are some differences between the two algorithms.The MOGE algorithm is designed on the basis of the gradient approximation, while the NSGA-II algorithm is designed based on the laws of biological evolution.The MOGE algorithm uses three operators to explore search space, improve the quality of the solution, avoid local optima, and promote population diversity.The NSGA-II algorithm uses two operators to explore search space, which are used for global search and local search, respectively.
This experiment used datasets of various sizes to verify the efficiency of the VNS-NSGA-II algorithm provided in this paper.The Pareto-optimal front and convergence of the VNS-NSGA-II, NSGA-II [38] and MOGE algorithms [46] were compared.For the experiment, 30, 50, 70, and 100 customers were chosen from the datasets.According to Wang et al. [35], the fresh agricultural products in this experiment are classified into four types based on their required temperatures.For simplicity, we use A, B, C and D for different types of products.For each type of agricultural product, a time-sensitive regulatory factor r was assigned [26].Table 2 shows the features of each scale dataset, including the number of customers, the number of vehicles, the total demand, the type of agricultural products, the value of time-sensitive adjustment factors, and the unit price of agricultural products.Figure 6 shows the Pareto-optimal front of each algorithm obtained in the experiment.Table 3 shows the optimal value of each objective function in the Pareto-optimal front of each algorithm.It should be noticed, however, that TTC in Table 3 represents the total cost.The optimal total cost and satisfaction values, as well as the number of iterations, is shown in Figure 7. Table 3 and Figure 7 show the total costs and satisfaction for two different solutions.Next, we will analyze these figures and tables in detail.

Analysis of Optimal Solution Selection
The primary purpose of the method suggested in this paper is to develop a clear and effective distribution scheme for fresh agricultural products.Thus, after generating the Pareto-optimal solution set using the VNS-NSGA-II algorithm, it is important to select a solution and decode it as a distribution plan using the TOPSIS method.In this experiment, the dataset of 100 customers is used as an example, and the TOPSIS method is used to select the most suitable solution from the Pareto-optimal front.Figure 8 illustrates the Pareto-optimal front of 100 customers.The horizontal coordinate of Figure 8 represents the total cost, and the vertical coordinate represents customer satisfaction.Each point in the graph represents a solution in the Pareto-optimal front.Table 4 shows the process of evaluating each solution in Figure 8 using the TOPSIS method.As indicated in Table 4, the TOPSIS method is utilized to determine the S + i , S − i and C i values for each solution.The filtered solution is compared to the solution with the ideal index values, as shown in Table 5. Note, however, that CC in Table 5 is the cost of carbon emissions and RC is the cost of refrigeration.The data in Tables 4 and 5 will be analyzed in detail next.We need to choose an appropriate solution based on the data in Table 4.Given that the objective function in this paper is to reduce the total cost and TOPSIS minimizes satisfaction, the smallest value should be chosen from the final closeness degree C i .As a result, solution No. 5 should be chosen as the final delivery strategy.Figure 8  In Table 5, the total cost of solution No. 5 is compared with the total cost of solution No. 1 and solution No. 20, which include the proportion of carbon emission cost to the total cost, the proportion of refrigeration cost to the total cost, the number of vehicles used, and the average load rate and satisfaction.Solution No. 1 has the highest customer satisfaction but also the highest total cost, and solution No. 20 has the lowest total cost but the worst customer satisfaction.The total cost and customer satisfaction of solution No. 5 are between solutions No. 1 and No. 20.It can be determined that the number of vehicles used must grow, while the average loading rate must drop in order to boost customer satisfaction.Additionally, more vehicles will lead to increased carbon emissions and cooling costs, and refrigeration costs will increase as well; on the other hand, to reduce total costs, the number of vehicles will be reduced, the average loading rate will be increased, and consequently, fewer delivery vehicles will delay the delivery of produce from the distribution center, lowering carbon emissions and cooling costs.
From the foregoing study, it is clear that the TOPSIS method can well screen out the suitable solutions from the Pareto-optimal solution set.

Optimisation of the Fresh Agricultural Products Distribution Routes for the e-Commerce Business in the Sample
In order to verify the effectiveness of the TDSDGVRPMTW model in reality, this study uses a case in Shanghai to compare and analyze the differences between the three delivery strategies.The current common distribution strategy for fresh agricultural products is to distribute all types of products in one vehicle at one temperature without split delivery or to use multi-compartment vehicles for multi-temperature distribution.According to the two delivery strategies and the mathematical model proposed in this paper, the two delivery strategies can be modeled as a time-dependent green vehicle routing problem with a time windows model (TDGVRPTW) and a time-dependent multi-compartment green vehicle routing problem with a time windows model (TDMCGVRPTW).Among them, the TDMCGVRPTW model adopts the method from Reed et al. [3].Therefore, the three strategies involved in the comparison are TDGVRPTW, TDMCGVRPTW and TDSDGVRPMTW.
We used the distribution data from an e-commerce business of fresh agricultural products in Shanghai.The data of one distribution center and 24 customers are shown in Table 6, including the locations of the distribution center, individual customers, demand, types of demand and time windows.Note, however, that the latitude and longitude of the customer's location are converted to X/Y coordinates for convenience.Types represent the types of fresh agricultural products that customers need.All types of fresh agricultural products in a refrigerated vehicle are loaded, and a certain temperature for delivery is set.The storage temperature of the vehicle is determined by the lowest temperature required among all types of fresh products.This is the company's current distribution strategy.The corresponding cost is calculated based on the company's actual distribution schemes and the constraints of the mathematical model established in this paper.According to the real vehicle data from the company, the total number of available vehicles is adjusted to Q L = 15 and the maximum capacity is adjusted to Q K = 1500.Other than that, all other parameters remain the same as in Section 2.6.Table 7 shows the relevant data for each model.Note that TC in Table 7 is the travel cost, FC is the fixed cost, and SC is the service cost.The relevant data of the TDGVRPTW model in Table 7 is the calculation result.The data corresponding to the TDMCGVRPTW model and the TDSDGVRPMTW model in Table 7 are obtained by using the VNS-NSGA-II algorithm and TOPSIS method proposed in this paper.Note, however, that Gap_mc in Table 7 refers to the difference between the related data of TDMCGVRPTW and TDGVRPTW, and Gap_sd refers to the difference between the related data of TDSDGVRPMTW and TDGVRPTW.Figure 9 is a comparison of all data of the three models.Each column in the figure represents the total cost of a model, where different colors represent different itemized costs.Broken lines represent customer satisfaction for different models.The scale on the left of the vertical coordinate corresponds to the value of total cost, while the scale on the right of the vertical coordinate corresponds to the value of satisfaction.The data presented in Table 8 can be used to analyze the different choices made by companies on the basis of different strategies when they have to give up part of the customer orders or the number of vehicles is limited.Table 8 shows the quantity and on-time rate of delivery of each type of agricultural product under three different strategies.On-time delivery means that the vehicle completes the delivery service for the customer within the customer's time window, under the circumstance of which the customer satisfaction is greater than 0. Note, that PD refers to the number of items delivered punctually.Proportion means the percentage of the quantity of products delivered punctually in the total demand for that type of product.Refrigeration cost per minute represents the refrigeration cost per minute for preserving the corresponding type of product.Total demand refers to the total demand of all customers ordering the corresponding type of product.As shown in Table 7 and Figure 9, the proportion of each cost in the total cost is roughly the same despite the different strategies taken, which shows that the delivery schemes produced by the method with two different models in this study are feasible.TDGVRPTW and TDMCGVRPTW have the same fixed cost, which means that the two strategies use the same number of vehicles.Compared with the RC of TDGVRPTW, the RC of TDMCGVRPTW and TDSDGVRPMTW has been greatly reduced, which shows that both the multi-compartment distribution strategy and the split delivery with a single compartment distribution strategy can greatly reduce the refrigeration cost.
However, it can be found that there are many differences between the three distribution strategies.As shown in Table 7 and Figure 9, the total cost and satisfaction of TDSDGVRPMTW are optimal among the three delivery strategies.The total cost of TDSDGVRPMTW is 5.81% lower than that in TDGVRPTW, and the satisfaction of TDS-DGVRPMTW is 10.68 points higher than that in TDGVRPTW.In other words, the distribution strategy of split delivery with a single compartment proposed in this paper can reduce the total cost and improve satisfaction.For the total cost, FC, RC and CC of TDSDGVRPMTW are all lower than those in TDGVRPTW.Although TDSDGVRPMTW is higher than TDGVRPTW in terms of TC and SC, the small difference will not trigger a large increase in TTC.Comparing TDMCGVRPTW with TDGVRPTW, it can be found that although the TTC of TDMCGVRPTW is 3.49% lower than that of TDGVRPTW, the customer satisfaction of TDMCGVRPTW is 8.21 points lower than that of TDGVRPTW.This means that the strategy reduced total costs at the expense of satisfaction.The RC of TDMCGVRPTW is lower than those of TDGVRPTW, but the TC, SC and CC are higher than those of TDGVRPTW.A high CC means that more carbon emissions are generated, which is not conducive to environmental protection.
As shown in Table 8, customers have the largest demand for A and the smallest demand for D. A has the lowest refrigeration cost, while D has the highest refrigeration cost.Three different delivery strategies show three different ways of allocating shipping capacity.The on-time rates of delivery of TDGVRPTW to A, B, and C are similar, which shows that the strategy will choose not to deliver the products with the lowest demand and the highest refrigeration cost when the number of vehicles is limited and some customer orders must be abandoned.TDMCGVRPTW has lower on-time rates of delivery for A and B than it does for C and D. This suggests that the strategy will choose to not deliver products that are in greater demand and have lower refrigeration costs when the number of vehicles is limited.This result is consistent with the findings of Hsu and Chen [2].However, this strategy produces very low satisfaction.TDSDGVRPMTW has lower on-time rates of delivery for C and D than for A and B, which suggests that the strategy will not deliver products with lower demand and higher refrigeration costs when the number of vehicles is limited.The on-time rate of delivery of this strategy for each product is much higher than that of TDGVRPTW, which is why the satisfaction of this strategy is higher than that of TDGVRPTW.
From the above analyses of results shown in this experiment, it can be concluded that the distribution strategy represented by TDSDGVRPMTW is most suitable for the realworld case presented in this study.There are many differences between the experimental results of TDSDGVRPMTW and TDMCGVRPTW.Comparing the results of this experiment with research on multi-compartment vehicle routing problems [2][3][4], it can be found that the reasons for these differences are as follows: 1.
Split delivery would allow customers ordering multiple agricultural products to be served by multiple vehicles, meaning that each vehicle would serve more customers and travel longer routes, which leads to higher travel costs and service costs.However, split delivery keeps each product at the optimum temperature for transport, which greatly reduces refrigeration costs and lowers the total cost.2.
Split delivery allows each vehicle to deliver a smaller number of products to customers, which leads to a shorter service time that makes the vehicle more likely to finish each delivery and leave the customer within the optimal service time window, resulting in higher customer satisfaction.Another reason for the high level of satisfaction is that TDSDGVRPMTW chooses to refuse orders with products of small quantity and prioritizes serving customers ordering products of large quantity.

3.
Divide a vehicle's compartment into multiple sections.Each of them transports one type of agricultural product with different optimum temperatures as needed.Although this can reduce refrigeration costs, if a customer's demand for a certain type of agricultural product exceeds the capacity of the divided compartment, multiple deliveries are required to meet the customer's demand for this type of agricultural product and more vehicles are needed, leading to high travel costs and carbon emissions.4.
Customers have multiple time windows to choose and therefore vehicles have more opportunities to arrive at locations and complete services during a certain time window.The advantages of multiple time windows over a single time window will be leveraged, especially when a customer needs to be served by vehicles multiple times.One customer in TDMCGVRPTW needs the service of multiple vehicles, however, each customer has only one time window, and consequently many vehicles arrive at the location of the customer out of the time window, leading to very low customer satisfaction.Another reason for the low satisfaction is that TDMCGVRPTW chooses to refuse many orders for items that are in high demand, leaving many customers unserved.5.
Therefore, the number of vehicle compartments, the capacity of each compartment, and the number of product types demanded by customers are the factors that determine why TDMCGVRPTW chooses not to deliver products that are in high demand and TDSDGVRPMTW chooses the opposite.The compartment of the vehicle in TDMCGVRPTW is divided into four parts.When a customer's demand for a certain product exceeds the compartment capacity, multiple vehicles are required to deliver the same product to that customer.When the number of vehicles needed to serve the customer exceeds the number of types of products the customer orders, the service cost is too high and the order will be refused by TDMCGVRPTW.Therefore, under this circumstance, the advantages of TDSDGVRPMTW over TDMCGVRPTW can be observed.In terms of product types, TDSDGVRPMTW needs fewer vehicles to complete the distribution to the customer with lower service costs, thus, this customer order will not be refused.When the number of such customers is large, TDSDGVRPMTW would naturally become the best strategy.
in this paper to model and select from different distribution strategies, so as to find the optimal one that suits their needs.
Comparing the experimental results in this study with those from Chen et al. [4], Reed et al. [3] and Hsu and Chen [2] can provide some insights into distribution management.The customer's demand, the number of types of products that the customer orders, and the capacity of each compartment in the multi-compartment delivery system, are the major factors that determine whether to use delivery by use of a multi-compartment vehicle or split delivery by a single compartment vehicle.As shown in the experiment of Reed et al. [3], it is reasonable to use multi-compartment vehicles for delivery when the customer requires several types of products and the customer's demand for a certain product does not exceed the capacity of each compartment.The experiment of a real case in this study shows that it is reasonable to split delivery with a single compartment vehicle when a customer requires a few types of products and the customer demand for a certain product exceeds the capacity of each compartment.The experiments in this paper also found that loading all types of fresh agricultural products in a vehicle and setting one single temperature for distribution is worse than the above two strategies.
On the other hand, this study has some limitations.In this paper, the calculation of carbon emissions does not consider the road slope, and it is not accurate enough since it only considers the vehicle speed, load weight and travel distance.Although the VNS-NSGA-II algorithm proposed has a good performance, there is still room for further optimization.For example, some heuristics are added in the initial population generation stage to improve the quality of the initial population.In this paper, the TOPSIS method is selected to screen suitable solutions from the solution set and to represent the effectiveness of this method.However, there are many multi-attribute decision-making methods, such as the elimination and choice expressing reality (ELECTRE) method, the preference ranking organization method for enrichment evaluations (PROMETHEE) method, and the analytical hierarchal process (AHP) method, and so on.In future, these methods can be applied to compare with the TOPSIS methods and to find more suitable ones for future research on TDSDGVRPMTW.In addition, it was found that the number of types of products needed by customers, the quantity of each ordered product, and the capacity of each compartment in a multi-compartment vehicle were factors that have an impact on what delivery strategy should be selected.Furthermore, how each factor affects the selection of delivery strategies is worth further study.

Figure 3 .
Figure 3.One example of a chromosome.

Figure 4 .
Figure 4. Examples of crossover and mutation.

Figure 5 .
Figure 5. Examples of variable neighborhood search.

Figure 7 .
Figure 7.Comparison of convergence of algorithms.
illustrates the location of solution No. 5 in the distribution of 100 customers' Pareto-optimal front.Solution No. 1 has the highest satisfaction value, while solution No. 20 has the lowest total cost value.While solution No. 5 provides 4.02 points less satisfaction than solution No. 1, the total cost is lowered by 17.4%.Although the total cost of solution No. 5 is 20.2% greater than that of solution No. 20, the satisfaction value of solution No. 5 is 30.54 points greater than that of solution No. 20.Taken together, solution No. 5 is an equilibrium solution in the middle of the weight setting's extreme values.

Figure 9 .
Figure 9.Total cost and satisfaction comparison of different models.

Table 1 .
Symbol definitions in the TDSDGVRPMTW optimization.

Table 2 .
Characteristics of datasets of different sizes.

Table 3 .
Experimental results of algorithm comparison.

Table 4 .
Evaluation of results using the TOPSIS approach.

Table 5 .
Comparison of selected solution and optimal value.

Table 6 .
Relevant data of distribution center and 24 customers.

Table 7 .
Comparison of three delivery strategies.

Table 8 .
The quantity of each type of product delivered punctually under three strategies.