Route Planning of Supermarket Delivery through Third-Party Logistics Considering Carbon Emission Cost

: China’s policies about reducing carbon emissions mainly focus on saving energy and reducing pollutant emissions. However, the carbon emissions of the logistics industry remain high. Thus, logistics companies should not only reduce distribution costs, but also consider the impact of carbon emissions when planning distribution routes. This paper studies the supermarket delivery distribution route planning problem considering carbon emissions. Aiming at the lowest economic operating cost composed of fixed cost, driving cost, and carbon emission cost, we propose a grey wolf optimized genetic algorithm to generate delivery route plans, where the social rank and hunting behavior of the grey wolf optimization algorithm are integrated into the selection operation of the genetic algorithm. Finally, a real case of a large third-party logistics enterprise in Beijing is studied. The case study results verify the effectiveness and applicability of the model and algorithm.


Introduction
With the increasing consumption demand of residents, supermarket distribution becomes increasingly important.Supermarket distribution refers to the logistics activities in which the distribution center delivers the goods to the designated supermarkets.The characteristics of supermarket distribution contain three aspects.First, the locations of the supermarkets are geographically dispersed; second, goods demand varies greatly among different supermarkets; and third, the delivery service time is strict [1].As a result, the delivery cost is high, including a massive energy consumption, excessive carbon emissions, and other problems [2,3].The Communist Party of China's Central Committee and the State Council unveiled a guiding document on the country's work to achieve carbon peaking and carbon neutrality goals in 2020.However, the carbon emissions of the logistics industry remain high.
Meanwhile, there are three main delivery modes, as follows: by suppliers, by supermarkets themselves, and by third-party logistics companies.The distribution efficiency of delivery by third-party logistics companies is high since the third-party logistics companies provide a joint delivery mode.The main advantage of joint delivery is that orders from different supermarkets can be delivered together.Third-party logistics companies adopt different processing strategies (i.e., splitting or merging orders) to complete delivery services.Thus, the third-party logistics delivery mode has rapidly developed in recent years [4,5].However, the logistics companies not only need to consider reducing distribution costs, but also the impact of carbon emissions, in line with the global trend of reducing carbon emissions.It is imperative to control and reduce the carbon emissions of third-party logistics distribution under the double carbon target.Thus, this paper proposes a grey wolf optimized genetic algorithm to plan the distribution routes between several supermarkets, from the view of the third-party logistics companies.The proposed algorithm aims to minimize the economic operating cost including the fixed cost, the driving cost, and carbon emission cost.As we know, carbon emission is high when the traffic is congested, while the traffic condition is time dependent.Thus, we further consider the correlation between the driving times, the driving velocity, and road segments.In addition, the constraint of a hard time window is also related to the requirements of strict service time.We conduct a series of experiments and analyze a real case to verify the effectiveness of the proposed model and the algorithm.Compared to the three classical genetic algorithms, our proposed grey wolf optimized genetic algorithm can reduce the economic operating cost, while considering the carbon emission cost.The convergence speed of the proposed method is faster, which can generate the high-quality supermarket distribution path plan more quickly.We also verify that the fluctuation of the carbon emission price has a great impact on the distribution scheme of supermarkets.
The rest of the paper is organized as follows: Section 2 is the literature review.We define the problem in a mathematical form in Section 3. Section 4 illustrates our proposed algorithm, which is based on the genetic algorithm optimized by the grey wolf algorithm.Section 5 is the case study.Finally, the paper is concluded in Section 6.

Literature Review 2.1. VRP Incorporating Carbon Emissions
Logistics will always consume a high amount energy and emit lots of carbon [6,7].Thus, low-carbon logistics and green logistics have attracted the attention of researchers in recent years [8,9].
The existing vehicle routing problem (VRP) works to either minimize carbon emissions costs or the delivery costs, including carbon emissions.Yu et al. [10] minimized carbon emissions and used an adaptive large neighborhood search algorithm to solve gasolineelectric hybrid distribution vehicle routing.Luo et al. [11] studied vehicle route planning for agricultural products in a cold chain, which took the sum of various distribution costs including carbon emission as the goal.Meanwhile, the traffic and weather conditions are also considered.With the multi-objective optimization, Zhu et al. [12] effectively solved the transportation and distribution problems in cold chain logistics.
The carbon emissions of logistics distribution mainly arise from the fuel consumption during the distribution vehicle's driving [13].The factors that affect the carbon emissions of logistics distribution mainly include vehicle load, driving speed, delivery distance, road network topology, traffic, etc. Suzuki et al. [14] established a vehicle routing model with capacity restrictions and vehicle number constraints, considering the impact of real-time vehicle load on carbon emissions; an heuristic algorithm is employed.Bai et al. [15] studied the vehicle routing problem of cold chain logistics, which considered the complexity of the road network and changing traffic conditions, which helps to improve the evaluation accuracy of distribution costs and carbon emissions.Zhang et al. [16] studied the dynamic delivery vehicle routing problem, which considered the impact of distribution distance on carbon emissions through deploying virtual customer points.
In summary, existing work plans the distribution routes based on the influencing factors of carbon emissions.However, they all assume that the driving times and velocity of the vehicles are fixed.In fact, the driving times and velocity are dynamic and are affected by traffic lights and road congestion in urban traffic environments.In a different manner to the existing work, this paper takes the vehicle load, driving times, and velocities into account, and aims to minimize the distribution cost including carbon emission costs.

Typical Solutions for Vehicle Routing Problems
The vehicle routing problem (VRP) is a typical combination optimization problem, which has been proven to be NP-hard.The existing works can be divided into precise algorithms, heuristic algorithms, and hybrid algorithms.The earliest methods are precise algorithms.For example, Yang et al. [6] adopted the branch-cut-pricing algorithm, which is a precise algorithm, providing a scheme for the last kilometer distribution in rural areas.Heuristic algorithms are more suitable for solving large-scale distribution vehicle routing problems.Heuristic algorithms include annealing algorithms, genetic algorithms, search algorithms, etc.For example, Liu et al. [17] employed a simulated annealing algorithm to solve the consistent vehicle routing problem.The genetic algorithm (GA) designed by Yu et al. [7] solves the routing problem of instant delivery vehicles based on customer classification.Hybrid algorithms refer to the combination of different heuristic algorithms and precise algorithms.For example, Fatemeh et al. combined the genetic algorithm and the simulated annealing algorithm to solve the green vehicle routing problem [9].
Existing research shows that heuristic algorithms are more effective in solving large distribution vehicle routing problems.Hybrid algorithms combine different heuristic algorithms, which can complement each other.Thus, hybrid algorithms can obtain a more reasonable distribution route solution.Our work follows the line of hybrid algorithms.As a new heuristic algorithm of swarm intelligence, grey wolf optimization has the characteristics of few dependent parameters and a good optimization performance, which is suitable for solving the routing problem of distribution vehicles; however, grey wolf optimization suffers from a complex coding mode.As a classical population evolution algorithm, a genetic algorithm is easy to code and to integrate with other algorithms.Therefore, we use a genetic algorithm supported by a grey wolf heuristic to solve the problem.

Problem Description
The route planning of supermarket delivery through third-party logistics considering carbon emission is formalized as follows: given the supermarkets, including the number, locations, service demands, and service hours, and a distribution center of the third-party logistics enterprise, the distribution center dispatches different types of vehicles to complete the delivery service.We aim to minimize the economic operating cost, which is composed of the fixed cost for employing a vehicle, the driving cost, and the carbon emission cost.
In order to clarify the scope of the paper, we assume the following: (1) We only consider the one-way goods delivery procedure, which is from the distribution center to supermarkets.The pick-up procedure is not involved [18].(2) The number of available vehicles is not limited.
(3) The demand of each supermarket does not exceed the maximum capacity of a vehicle.(4) All requirements of the supermarkets are met and each supermarket is served by one vehicle once.(5) All supermarkets' windows of service hours are strict.That is, if the goods are not delivered to the supermarket within the specified service hours, the service will be refused.(6) The vehicles are subject to local restrictions and traffic management measures, such as vehicle type restriction.(Vehicle number and the type of restriction are measures taken by the local government to alleviate traffic congestion, reduce vehicle exhaust emissions, and achieve air quality standards.For instance, vehicle type restriction means that vehicles with limited types cannot enter the corresponding area.Violators will be fined).
Table 1 shows the symbols used in this paper.Driving speed of the k-th vehicle of type l from supermarket i to j x ij lk Variable 0-1: If the vehicle arrives at supermarket j from i, the value is 1; otherwise, the value is 0 y i lk Variable 0-1: If the supermarket i is served by the vehicle, the value is 1; otherwise, the value is 0 r i h Variable 0-1: If the supermarket i is in the h-th ring of the city, the value is 1; otherwise, the value is 0

Objective Function
Our algorithm aims to generate a distribution plan where all requirements of the supermarkets are met, while the economic operation cost is minimized.The economic operation cost is composed of the fixed cost for employing vehicles, the driving cost, and the carbon emission cost.The objective function is shown in Equation (1).
where Z represents the economic operation cost of the distribution plan.The first term in Equation ( 1) is the fixed cost for employing vehicles, which depends on the number of vehicles used in the plan.The fixed cost is equal to the product of the number of vehicles and the fixed cost of a vehicle.The second term is the driving cost.The driving cost is the product of the petrol consumption of the vehicles and the fuel price.The third term is the cost for the carbon emission of the delivery vehicles.The carbon emission cost is the product of the carbon emissions and the carbon price.The second and third terms will be described in Section 3.3 in detail.

Constraints
The constraints include three aspects in this paper; that is, from the view of vehicles' operation, from the view of supermarkets' requirements, and from the view of local restrictions and traffic management measures.
(1) From the view of vehicles' operation Equation (2) indicates that the total demand of the supermarkets served by each delivery vehicle cannot exceed the vehicle capacity.Equation (3) indicates that all delivery vehicles should depart from the distribution center and eventually return to the center.Equation (4) indicates that the timestamp when the vehicle leaves the supermarket i equals the time when the supermarket is served after the vehicle arrives.Equation (5) indicates that the time when the vehicle arrives at supermarket j is the time when the vehicle leaves supermarket i, as well as the driving times from i and j.
(2) From the view of supermarkets' requirements Equation ( 6) is a constrained condition, which means that the number of vehicles reaching the supermarket equals the number of vehicles leaving the supermarket, and the same can be said for the number of vehicles serving the supermarket.Equation ( 7) indicates a hard time window required by the supermarket.That is, the arrival time of the delivery vehicle at the supermarket must be within the range of the required earliest service time and the latest service time.
(3) From the view of the local restrictions traffic management measures In this paper, we only consider type restriction.Equation (8) requires that only the vehicle whose maximum capacity is not larger than the allowed maximum capacity requirement inside the h-th loop of the city can serve the supermarket.

Fuel Consumption and Carbon Emission Evaluation
We employ the widely used CMEM (Comprehensive Modal Emission Model) [19] to evaluate the fuel consumption and the carbon emissions caused by the supermarket distribution plan, which considers the effects of vehicle load, driving times, and speed.The evaluation equation is shown as follows: where β 1 , β 2 , and β 3 are the predetermined parameters, which represent the vehicle engine coefficient, the speed coefficient, and the load coefficient, respectively; with reference to [19], the three parameters are set as 0.22, 0.036, and 0.017, respectively.It has been proven that vehicle carbon emissions and fuel consumption hold a linear relationship.Therefore, the carbon emission of the delivery vehicle driving from supermarket i to j is shown in Equation ( 10): where ψ represents the fuel emission factor, which is set as 2.62 in this paper, as in Ref. [13].
Most scholars assume that the first two items in Equation ( 9), driving times t lk ij and the driving speed v lk ij , are fixed.However, in fact, these are affected by external factors such as traffic light and road congestion in the urban traffic environment; the driving times and the driving speed of the delivery vehicles vary according to different timestamps.Therefore, in this paper, we calculate the fuel consumption and the carbon emissions of the vehicle based on the diverse driving times with respect to their different locations.
Ref. [20] verifies the lognormal distribution to be the best to fit vehicle driving times within 24 h, using 36 measured data sets.Therefore, we divide every two hours in a day into a period.Assume that the driving times t lk ij , the k-th vehicle of type l from supermarket i to j, are subject to a lognormal distribution in each period.That is, ln(t lk ij ) ∼ N(µ ij , σ 2 ij ).µ ij and σ 2 ij are the mean and the variance of the vehicle driving times, respectively.We also use the case of Sinotrans Beijing Distribution Center, whose detailed information can be found in Section 5, to test whether the driving times of delivery vehicles are subject to a lognormal distribution in each period using the K-S test.First, the driving times between the distribution center and the designated 20 supermarkets are retrieved every 15 min between 6:00 and 22:00 for 15 consecutive days using Baidu map API.The unit of driving time is minutes.Then, each day from 6:00 to 22:00 is divided into eight periods, every two hours.The retrieved driving times and driving speeds are split into each period.After that, the frequency distribution histogram of the logarithmic driving times of the vehicles is drawn.Taking the ordinary period (10:00~12:00) and the peak period (18:00~20:00) as examples, the histogram of vehicle driving times after taking logarithm are shown in Figure 1.The K-S test results are shown in Table 2.
where  represents the fuel emission factor, which is set as 2.62 in this paper, as in Ref [13].
Most scholars assume that the first two items in Equation ( 9), driving times  and the driving speed  , are fixed.However, in fact, these are affected by external factors such as traffic light and road congestion in the urban traffic environment; the driving times and the driving speed of the delivery vehicles vary according to different timestamps.Therefore, in this paper, we calculate the fuel consumption and the carbon emissions of the vehicle based on the diverse driving times with respect to their different locations.
Ref [20] verifies the lognormal distribution to be the best to fit vehicle driving times within 24 h, using 36 measured data sets.Therefore, we divide every two hours in a day into a period.Assume that the driving times  , the k-th vehicle of type l from supermarket i to j, are subject to a lognormal distribution in each period.That is, ln  ~N  ,  . and  are the mean and the variance of the vehicle driving times, respectively.
We also use the case of Sinotrans Beijing Distribution Center, whose detailed information can be found in Section 5, to test whether the driving times of delivery vehicles are subject to a lognormal distribution in each period using the K-S test.First, the driving times between the distribution center and the designated 20 supermarkets are retrieved every 15 min between 6:00 and 22:00 for 15 consecutive days using Baidu map API.The unit of driving time is minutes.Then, each day from 6:00 to 22:00 is divided into eight periods, every two hours.The retrieved driving times and driving speeds are split into each period.After that, the frequency distribution histogram of the logarithmic driving times of the vehicles is drawn.Taking the ordinary period (10:00~12:00) and the peak period (18:00~20:00) as examples, the histogram of vehicle driving times after taking logarithm are shown in Figure 1.The K-S test results are shown in Table 2.The significance level is 0.05.According to the density curve of the frequency distribution histogram in Figure 1, the logarithm of the driving times approximately and intuitively follows a normal distribution.In Table 2, the p-values of the two periods were 0.716 and 0.486, respectively.The null hypothesis is accepted.Therefore, the driving times of the vehicles in the two periods approximately followed the lognormal distribution.That is, during 10:00~12:00, ln(t lk ij ) ∼ N(7.38,0.05 2 ) and during 18:00 to 20:00, ln(t lk ij ) ∼N(7.62,0.03 2 ).Table 3 shows the K-S test results of vehicle driving times' distribution at different periods between two locations.We use the mean of the driving times to compute driving speed, v ij lk , as shown in Equation ( 11): In Equation ( 11), d ij is the road network distance between supermarket i and supermarket j, and µ ij is the average vehicle driving times between the two locations at different periods.

Grey Wolf Optimized Genetic Algorithm
As we know, the number of supermarkets served by third-party logistics enterprises is large.Their locations are sparse and their demands are quite diverse.Therefore, we propose a grey wolf optimized genetic algorithm to solve the supermarket distribution route plan problem of third-party logistics enterprises.The grey wolf optimized genetic algorithm is a hybrid heuristic algorithm combining the genetic algorithm and the grey wolf optimization idea.The genetic algorithm is simple and easy to combine with other algorithms.At the same time, the grey wolf optimization algorithm is less dependent on parameters and has a good optimization performance.Thus, we propose to integrate the social rank and hunting behavior in the grey wolf optimization algorithm into the selection operation of the genetic algorithm, which is suitable for solving complex distribution path problems.For ease of understanding, Table 4 gives the common terms and the correspondence meanings in our proposed new algorithm.
Figure 2 shows the flow chart of the proposed algorithm.At the beginning of the algorithm, the distribution center and the supermarkets are coded based on the chromosome coding.As a result, the distribution centers and supermarkets are mapped to genes.This is the basis for the following steps.The chromosome, i.e., a grey Wolf, is a supermarket distribution path plan.The chromosome is composed of genes (which are the distribution center and supermarkets).This paper uses integer coding to encode genes in chromosomes, where "0" represents distribution centers and "1, 2, . .., N" are supermarkets.After encoding, the initial grey wolf population is generated.That is, several supermarket distribution path plans are generated randomly.The specific procedure will be described in detail in Section 4.1.
Secondly, the fitness of each grey wolf in the population was calculated.Fitness is computed using the objective function; that is, the economic operating cost of the plan.The distribution route plans are checked based on the fitness, such that the social classifications of the grey wolves are classified.Then, the grey wolves start to hunt.In the first generation, all grey wolves (i.e., supermarket distribution path plans) regard the first three wolves with the highest social rank (i.e., the first three distribution route plans with the lowest economic operating cost) as the target to update locations.Through hunting behavior, it is possible to renew the location of the target for grey wolves.The two grey wolves closest to the target are found and the subsequent selection operations are conducted.The specific social hierarchy criteria and hunting behavior rules of the grey wolf will be described in detail in Section 4.2.Secondly, the fitness of each grey wolf in the population was calculated.Fitness is computed using the objective function; that is, the economic operating cost of the plan.The distribution route plans are checked based on the fitness, such that the social classifications of the grey wolves are classified.Then, the grey wolves start to hunt.In the first generation, all grey wolves (i.e., supermarket distribution path plans) regard the first three wolves with the highest social rank (i.e., the first three distribution route plans with the lowest economic operating cost) as the target to update locations.Through hunting behavior, it is possible to renew the location of the target for grey wolves.The two grey wolves closest to the target are found and the subsequent selection operations are conducted.The specific social hierarchy criteria and hunting behavior rules of the grey wolf will be described in detail in Section 4.2.
After that, through the genetic operators such as selection, crossover, and mutation, the population are evolved.That is, a new supermarket distribution route plan is generated.In order to improve the local search ability, the genetic operator operation of the proposed algorithm is based on the hunting behavior rule of the grey wolf.The specific operation will be described in detail in Section 4.3.
Finally, the termination condition is tested.The termination condition is the number of evolutionary generations (i.e., number of iterations).If the termination condition is not met, the above steps are repeated.Otherwise, the best individual in the last generation of the grey wolf population is output; that is, the best supermarket distribution route plan.

Term Meaning Meaning Description
A grey wolf/chromosome A candidate solution.
A supermarket distribution route plan.

Gene
Components of a grey wolf/chromosome, component in the candidate solution.
Distribution centers and supermarket'stores.
Grey wolf population A group of grey wolves, a collection of candidate solutions.
A collection of supermarket distribution route plans.

Population size
The number of grey wolves in the population and the number of candidate solutions in the set.
The number of supermarket distribution route plans generated each time.

Fitness
Measurement for the superiority and inferiority of the grey wolf population to carry out social classification.
The index to judge the best distribution route plan of supermarkets according to the objective function

Social hierarchy of the grey wolf α Wolf
The best solution in the current candidate solutions.
The best solution in the set of supermarket distribution route plans.
β Wolf The second best solution in the current candidate solutions.
The second best solution in the set of supermarket distribution route plans.
δ Wolf The third best solution in the current candidate solutions.
The third best solution in the set of supermarket distribution route plans.
ω Wolf The remaining solutions in the current candidate solution set.
Other solutions in the set of supermarket distribution route plans.

Prey
The best solution.
The best supermarket distribution route plan.

Surround
The distance between individual grey wolves and prey is calculated for grey wolf location update.
Calculate the difference between the fitness value of each supermarket distribution path scheme and the fitness value of the global best scheme.
Hunt ω wolves update the locations based on α, β, and δ wolves' locations.
The top three solutions are assumed to be the current best solutions and the target fitness of the other solutions is calculated.The process of obtaining a new set of distribution solutions for supermarkets.

Select
Selection acts on the grey wolf population to achieve genetic purposes and to prepare for crossover and mutation.
The selection operation is carried out according to the law of grey wolf hunting behavior.

Crossover
Crossover acts on a grey wolf population to produce new individuals in the population.
The crossover operation is carried out using the sequential crossover method.

Mutation
Mutation acts on a grey wolf population to produce new individuals in the population.
The mutation operation was carried out using the binary mutation method.
After that, through the genetic operators such as selection, crossover, and mutation, the population are evolved.That is, a new supermarket distribution route plan is generated.In order to improve the local search ability, the genetic operator operation of the proposed algorithm is based on the hunting behavior rule of the grey wolf.The specific operation will be described in detail in Section 4.3.
Finally, the termination condition is tested.The termination condition is the number of evolutionary generations (i.e., number of iterations).If the termination condition is not met, the above steps are repeated.Otherwise, the best individual in the last generation of the grey wolf population is output; that is, the best supermarket distribution route plan.

Population Generation Initialization
The initial population must meet the maximum vehicle capacity constraint and the service time window requirement of the supermarket.
First, supermarket visiting sequences are generated randomly.Then, the time of arrival and departure of delivery vehicles at each supermarket in the sequence and the load of delivery vehicles during the service procedure are calculated.After that, route segmentation is carried out for the supermarkets.The constraints of the route segmentation are the maximum capacity constraint of delivery vehicles and the service time window requirements of supermarkets.For each supermarket in the sequence, if the two conditions are met, the supermarket will be assigned a distribution vehicle.The route segmentation procedure is stopped until all supermarkets are served.The initial grey wolf population is obtained (i.e., sets of supermarket distribution route plans).

Grey Wolf Social Hierarchy and Hunting Behavior Rules
The social classification of grey wolves includes α, β, δ, and ω levels, according to the fitness.Each grey wolf is classified as one kind of wolf.There is only one α-wolf, one β-wolf, and one δ-wolf.The α-wolf, the β-wolf, and the δ-wolf correspond to the first, second, and third supermarket distribution route plans ranked by the fitness.ω-wolves correspond to the other supermarket distribution route plans.
The traditional hunting behavior of grey wolves includes surrounding, hunting, attacking, and searching.We adapt two activities, i.e., surrounding and hunting.The hunting behavior of the grey wolves depends on the social hierarchy.ω wolves take α, β, and δ wolves as targets to update locations.The specific behavior is as follows: (1) Surrounding Surrounding is the way in which grey wolves approach their prey when the prey's location is known.In our problem, surrounding means that in a generation population (i.e., sets of supermarket distribution route plans), a route plan is selected as the prey (i.e., the best supermarket distribution route plan) and other plans are adjusted as the best one.The adjustment method is as follows: Z(t) is fitness of the current grey wolf; Z p (t) is the fitness of the prey; and t is the iteration number.D is the fitness gap between grey wolves and their prey.Formula ( 13) is the fitness after updating locations, while 2ar 1 − a and 2r 2 are the coefficient vectors.The coefficient vectors help to avoid falling into the local optimal solution [21].In Equation ( 13), a is the convergence factor, which decreases linearly from 2 to 0; r 1 and r 2 are random numbers between [0, 1].
(2) Hunting We assume that α, β, and δ wolves know more about their prey, and the ω wolf pursues prey based on the fitness of these three wolves.The hunting method is as follows: A 1 (C 1 ), A 2 (C 2 ), and A 3 (C 3 ) are the coefficient vectors for α, β, and δ wolves, respectively.As in Equations ( 12) and ( 13), A = 2ar 1 − a, and C = 2r 2 , as in Equation (13).In Equations ( 14)-( 16), the current best top three supermarket distribution route plans are regarded as prey.The plan is adjusted according to the adjustment method in the surrounding.Equation ( 16) means that the average fitness value is taken as the adjustment target for the fitness of the remaining route plans.

Selection
The purpose of selection is to select the distribution route plan with the highest fitness during each evolution, which is the basis for crossover and mutation.We selected the plan based on the social classification and hunting behavior rules of grey wolves.It mainly consists of two steps.First, the grey wolf, which pass down to the next generation (a supermarket distribution route plan) directly, is selected.Second, the grey wolves that are crossed and mutated are selected.
The specific steps are as follows: First, the α, β, and δ wolves are inherited by the next generation directly.That is, the top three route plans with the best fitness in the current generation are retained.Then, based on the fitness adjustment target of grey wolf hunting behaviors, the two grey wolves closest to the target fitness, i.e., the two supermarket distribution route plans with the closest economic cost, were selected for local adjustment.

Crossover
The aim of crossover is to generate a new supermarket distribution route plan.Typical crossover operators include partial matching crossover, sequential crossover, cyclic crossover, etc.This paper adopts a sequential crossover strategy to carry out crossover operations.
The basic steps are as follows: First, two grey wolves (i.e., supermarket distribution route plans) generated from the previous selection operation are chosen as parents.Then, several genes in one parent are left unchanged and the remaining genes are crossed with the genes in the other parent, in order.The filial generation is obtained.In other words, the service order of some supermarkets in a route plan are kept unchanged and the service order of other supermarkets changes according to the service order in another route plan.The process diagram of sequential crossing is shown in Figure 3.Each crossover will produce two children (i.e., the new supermarket distribution route plans).
The aim of crossover is to generate a new supermarket distribution route plan.Typical crossover operators include partial matching crossover, sequential crossover, cyclic crossover, etc.This paper adopts a sequential crossover strategy to carry out crossover operations.
The basic steps are as follows: First, two grey wolves (i.e., supermarket distribution route plans) generated from the previous selection operation are chosen as parents.Then, several genes in one parent are left unchanged and the remaining genes are crossed with the genes in the other parent, in order.The filial generation is obtained.In other words, the service order of some supermarkets in a route plan are kept unchanged and the service order of other supermarkets changes according to the service order in another route plan.The process diagram of sequential crossing is shown in Figure 3.Each crossover will produce two children (i.e., the new supermarket distribution route plans).

Mutation
Mutation generates new supermarket distribution route plans as well, whose function is to expand the search range and avoid falling into the local best solution.In a different manner to the traditional mutation only requiring one operator, we adopt a binary mutation strategy to carry out mutation operations.That is, we require two supermarket distribution route plans to participate.
The basic steps are as follows: First, two grey wolves generated from the previous selection operation are chosen as parents.Second, a different gene (i.e., a supermarket) from the two grey wolves was selected randomly.Then, the locations of the two selected

Mutation
Mutation generates new supermarket distribution route plans as well, whose function is to expand the search range and avoid falling into the local best solution.In a different manner to the traditional mutation only requiring one operator, we adopt a binary mutation strategy to carry out mutation operations.That is, we require two supermarket distribution route plans to participate.
The basic steps are as follows: First, two grey wolves generated from the previous selection operation are chosen as parents.Second, a different gene (i.e., a supermarket) from the two grey wolves was selected randomly.Then, the locations of the two selected genes exchange.That is, the service order of the different selected supermarkets is exchanged, such that two new supermarket distribution route plans are generated.The process of a mutation operation is shown in Figure 4. Similarly to crossover, the process of each binary mutation will produce new supermarket distribution route plans.genes exchange.That is, the service order of the different selected supermarkets is exchanged, such that two new supermarket distribution route plans are generated.The process of a mutation operation is shown in Figure 4. Similarly to crossover, the process of each binary mutation will produce new supermarket distribution route plans.

Data Analysis
We use the data from Sinotrans (https://www.sinotrans.com/), a pioneer in China's third-party logistics industry, as the studying case.The delivery demands of the supermarkets on one day in Beijing were used to verify the model and the algorithm.The distribution center is marked as a star and 20 supermarkets are marked as small red and blue balls.The distances between the distribution center and the supermarkets are retrieved using Baidu map API.The distances are drawn as a hot map, as Figure 5 shows.From Figure 5, the distances between the distribution center and each supermarket varies from 0 to 80 km.More than half of the supermarkets are within 50 km of each other.That is because most of the supermarkets served by Sinotrans are in the sixth loops, which are near the shopping malls with large passenger flow in Beijing.A small number of supermarkets are located in the suburbs outside the sixth loop.

Case Study 5.1. Data Analysis
We use the data from Sinotrans (https://www.sinotrans.com/), a pioneer in China's third-party logistics industry, as the studying case.The delivery demands of the supermarkets on one day in Beijing were used to verify the model and the algorithm.The distribution center is marked as a star and 20 supermarkets are marked as small red and blue balls.The distances between the distribution center and the supermarkets are retrieved using Baidu map API.The distances are drawn as a hot map, as Figure 5 shows.From Figure 5, the distances between the distribution center and each supermarket varies from 0 to 80 km.More than half of the supermarkets are within 50 km of each other.That is because most of the supermarkets served by Sinotrans are in the sixth loops, which are near the shopping malls with large passenger flow in Beijing.A small number of supermarkets are located in the suburbs outside the sixth loop.
using Baidu map API.The distances are drawn as a hot map, as Figure 5 shows.From Figure 5, the distances between the distribution center and each supermarket varies from 0 to 80 km.More than half of the supermarkets are within 50 km of each other.That is because most of the supermarkets served by Sinotrans are in the sixth loops, which are near the shopping malls with large passenger flow in Beijing.A small number of supermarkets are located in the suburbs outside the sixth loop.The specific locations of the distribution and the supermarkets and their orders' demands are shown in Table 5.The demand of each supermarket varies from [13,1604] pieces.The demands of all supermarkets are quite different.The available vehicle types are listed in Table 6.The traffic restriction policy claims that trucks larger than 8 t are prohibited within the Fifth Ring Road of Beijing.Thus, we divide the 20 supermarkets into two groups, as follows: Group 1, i.e., the supermarkets within the Fifth Ring Road, includes supermarkets 4, 5, 6, 7, 12, 14, and 15; Group 2, i.e., the supermarkets outside the Fifth Ring Road, includes the rest of the supermarkets.Considering the restrictions on vehicle types, the first group of supermarkets can only have a delivery from 6.8 t vehicles.For the second group, the 9.6 t vehicle is guaranteed to be fully loaded and the 9.6 t or 6.8 t model will be randomly selected for the remaining demand.The specific locations of the distribution and the supermarkets and their orders' demands are shown in Table 5.The demand of each supermarket varies from [13,1604] pieces.The demands of all supermarkets are quite different.The available vehicle types are listed in Table 6.The traffic restriction policy claims that trucks larger than 8 t are prohibited within the Fifth Ring Road of Beijing.Thus, we divide the 20 supermarkets into two groups, as follows: Group 1, i.e., the supermarkets within the Fifth Ring Road, includes supermarkets 4, 5, 6, 7, 12, 14, and 15; Group 2, i.e., the supermarkets outside the Fifth Ring Road, includes the rest of the supermarkets.Considering the restrictions on vehicle types, the first group of supermarkets can only have a delivery from 6.8 t vehicles.For the second group, the 9.6 t vehicle is guaranteed to be fully loaded and the 9.6 t or 6.8 t model will be randomly selected for the remaining demand.

Analysis of Results
We use the proposed grey wolf optimized genetic algorithm to solve the above practical distribution route problem.Python 3.8 was used to implement the algorithm.This paper verifies the effectiveness of the model and algorithm from the following three aspects: the influence of different parameters, the comparison of different algorithms, and the change in carbon emission price.

Impact Analysis of Different Algorithm Parameters
The parameter setting has an important effect on the overall performance.With reference to the parameter setting experience of previous genetic algorithms, the crossover probability was set to 0.8 and the mutation probability was set to 0.2.The population size was 50, 75, 100, 125, and 150, and the number of iterations was 50, 100, 150, 200, and 250, respectively.The economic operation cost of the supermarket distribution scheme under different parameter combinations is shown in Table 7. From Table 7, when the population size is 100 and the iteration number is 200, the economic operation cost is minimized.Through our proposed grey wolf optimized genetic algorithm, a distribution route plan is obtained, as shown in Table 8.The two numbers in the brackets under the column "delivery time" represent the time of arrival at or departure from the distribution center or supermarkets, respectively.In order to visualize the delivery routes, we transformed the latitudes and longitudes into the coordinates in the Cartesian coordinate system.Then, the distribution vehicle routes are shown in Figure 6.Each line with one color is a vehicle route.ber of supermarkets a delivery vehicle serves for varies greatly.The difference mainly comes from the strict service window requirements.In order to ensure the provision of a distribution service within the required time window, the distribution vehicles cannot be fully loaded, resulting in the waste of distribution resources.This shows that the time window requirement of supermarket stores is an important factor affecting the distribution route plan of supermarket stores.From Table 8 and Figure 6, a total of 6 vehicles are needed to complete the distribution demand of 20 supermarkets, including 3 vehicles of 9.6 t and 6.8 t, respectively.The number of supermarkets a delivery vehicle serves for varies greatly.The difference mainly comes from the strict service window requirements.In order to ensure the provision of a distribution service within the required time window, the distribution vehicles cannot be fully loaded, resulting in the waste of distribution resources.This shows that the time window requirement of supermarket stores is an important factor affecting the distribution route plan of supermarket stores.

Comparative Analysis with the Genetic Algorithm
In order to verify its effectiveness, the proposed method is compared with the classical genetic algorithm.The difference between the proposed method and the classical genetic algorithm mainly lies in the selection operation.The classical genetic selection operators include roulette selection, binary tournament selection, linear sorting selection, etc.Therefore, we compare the genetic algorithm with the above classical genetic selection operators with our proposed methods under the same parameter settings.The comparison results are shown in Table 9 and the population evolution trend of the four algorithms is shown in Figure 7.
Obviously, the cost of carbon emissions increases linearly with the increase in carbon prices in Figure 8e.With the increase in carbon emission price, fuel consumption, carbon emission, and driving cost show the same trends as shown in Figure 8a,b,d.The three variations decrease first and then stabilize as the carbon emission price is larger than 5. From Table 9, we observe that the vehicle number increases by 1 with the price of carbon emissions increasing in order to serve all customers.As a result, fuel consumption is reduced since more vehicles can improve the service efficiency, while the fixed cost increases.Thus, the helpfulness of increasing vehicle numbers is limited.
From Figure 8c, with the price of carbon emissions increases, the fixed cost of vehicle use shows a "step-like" increasing trend.This is because the fixed cost is related the number of vehicles.When the carbon emission price changes at [0, 5.0] CNY/kg, the number of vehicles used remains unchanged, so the fixed cost remains unchanged.When the carbon emission price continues to increase, in order to reduce carbon emission costs, companies increase the number of vehicles used, resulting in a sharp increase in the fixed costs of vehicle use.From Figure 8c, with the price of carbon emissions increases, the fixed cost of vehicle use shows a "step-like" increasing trend.This is because the fixed cost is related to the number of vehicles.When the carbon emission price changes at [0, 5.0] CNY/kg, the number of vehicles used remains unchanged, so the fixed cost remains unchanged.When the carbon emission price continues to increase, in order to reduce carbon emission costs, With the increase in carbon emission price, the economic operating cost of the supermarket distribution scheme shows an increasing trend, as shown in Figure 8f.When the carbon emission price changes in the range of [0, 5.0] and [5.5, 7.0] CNY/kg, the change in trend of economic operating cost is only slight, which is mainly affected by the driving cost and carbon emission cost.The trend changes sharply when the carbon emission price is in the range of [5.0, 5.5] CNY/kg, due to the fixed cost of vehicle use sharply increasing.
From the above analysis, carbon emission price has a great impact on the supermarket distribution route plan.Third-party logistics enterprises should formulate reasonable supermarket distribution route schemes according to different carbon emission prices.When the carbon emission price is in the range of [0, 2.0] CNY/kg, the enterprises should give priority to reducing the fixed cost and driving cost, because the carbon emission cost accounts for a relatively small proportion of the economic operation cost, about 1% to 4%.When the carbon emission price is in the range of [2.0, 5.0] CNY/kg, the carbon emission cost accounts for a relatively large part of the economic operation cost.As a result, fuel consumption and carbon emission have significant changes.Therefore, the enterprises should focus on reducing the carbon emission cost.When the carbon emission price is in the range of [5.0, 7.0] CNY/kg, the enterprises should balance the fixed cost, driving cost, and carbon emission cost to ensure the lowest economic operating cost.

Conclusions
Through analyzing the characteristics of the route planning of supermarket deliveries through third-party logistics, this paper proposes a grey wolf optimized genetic algorithm for supermarkets distribution route planning, considering carbon emission cost.The main conclusions are as follows: (1) When building the supermarket distribution route planning model, carbon emissions are calculated based on dynamic driving times on the road networks.The route planning can reduce the fuel consumption and carbon emissions of the supermarket distribution service.(2) Compared with the three classical genetic algorithms, the proposed grey wolf optimized genetic algorithm reduces the economic operating cost of CNY 98.4, CNY 32.4, and CNY 82.1, respectively.The convergence speed is faster, which can reach the high-quality supermarket distribution path plan more quickly.(3) The fluctuation of carbon emission price has a great impact on the distribution scheme of supermarkets.Therefore, it is suggested that third-party logistics enterprises develop supermarket distribution route schemes applicable to different carbon emission prices, so as to reduce carbon emissions and improve profit margins.
This paper only considers the delivery problem with one single distribution center, without time constraints.In the future, we will expand the study on the distribution path of supermarkets with multiple distribution centers for simultaneous pickup and delivery with time constraints and the order of vehicles selection on the different route planning of supermarket delivery.

Figure 1 .
Figure 1.Frequency distribution histogram of distribution vehicle travel time logarithm.(a) Ordinary period; (b) Peak hours.

Figure 1 .
Figure 1.Frequency distribution histogram of distribution vehicle travel time logarithm.(a) Ordinary period; (b) Peak hours.

Figure 2 .
Figure 2. Flow chart of grey wolf optimized genetic algorithm.

Figure 2 .
Figure 2. Flow chart of grey wolf optimized genetic algorithm.

Figure 5 .
Figure 5.A distances hot-map between the distribution center and the supermarkets.

Figure 5 .
Figure 5.A distances hot-map between the distribution center and the supermarkets.

Figure 8 .
Figure 8. Change in trend of various indicators under different carbon emission prices.(a) Fuel consumption; (b) Carbon emissions; (c) Fixed cost; (d) Vehicle running cost; (e) Carbon emission cost; (f) Economic operation cost.

Figure 8 .
Figure 8. Change in trend of various indicators under different carbon emission prices.(a) Fuel consumption; (b) Carbon emissions; (c) Fixed cost; (d) Vehicle running cost; (e) Carbon emission cost; (f) Economic operation cost.

Table 2 .
Hypothesis test summary using K-S test.

Table 2 .
Hypothesis test summary using K-S test.
Decision Maker 1-The distribution of travel time log from 10 o'clock to 12 o'clock is normal, the mean value is 7.38, and the standard deviation is 0.05

Table 3 .
Distribution of vehicle driving times between two locations at different times.

Table 5 .
Distribution center and supermarket stores information.

Table 6 .
Distribution center vehicle types.

Table 7 .
Economic operation cost under different parameter combinations.