Research on the Time-Dependent Vehicle Routing Problem for Fresh Agricultural Products Based on Customer Value

Abstract

With continuous improvements in people’s consumption levels, consumers’ demands for safe and fresh agricultural products increase. The increase in the number of vehicles and serious congestion on roads has led to problems, such as the weak timeliness of urban cold chain logistics, high carbon emissions, low customer value and reduced customer satisfaction. In this study, carbon emissions, customer satisfaction, customer value and cost are considered, and an optimization algorithm is established to solve the time-dependent vehicle routing problem in urban cold chain logistics. For road congestion at different time periods during the cold chain distribution process, the segment function is used to express the vehicle speed. According to the characteristics of the model, considering the constraints of the time window and vehicle capacity, an improved NSGA-II algorithm with the local optimization characteristics of the greedy algorithm (G-NSGA-II) is proposed, and the sorting fitness strategy is optimized. In addition, we carry out a series of experiments on existing vehicle routing problem examples and analyze them in a real background to evaluate and prove the effectiveness of the proposed model and algorithm. The experiment results show that the proposed approach effectively reduces the total cost, enhances customer value and promotes the long-term development of logistics companies.

1. Introduction

Due to the increase in consumer health awareness, the importance of the quality of fresh food in daily life has been highlighted. The huge demand for fresh food has also driven the development of cold chain logistics. According to the prediction of the China Business Industry Research Institute, the market scale of China’s cold chain logistics business is likely to exceed RMB 550 billion by 2022 [1]. The International Energy Agency (IEA) reports that 23% of global CO₂ emissions are generated by the transportation sector, in which nearly 75% of emissions is generated by the road sector [2]. Fresh food is mainly distributed using large fuel-guzzling vehicles. The cold logistics industry has gradually become a major energy consumer and carbon emitter due to its rapid development. Therefore, the study of the vehicle routing problem is significant in reducing CO₂ emissions, achieving green logistics and promoting economic development. Through energy conservation and emission reduction, cold chain logistics actively responds to the “14th Five-Year Plan” to achieve the “double carbon goal”. It is of great social significance to study the optimization of the cold chain path from a green perspective.

The vehicle routing problem (VRP) was first proposed by Dantzig and Ramser [3]. In order to adapt to social developments and meet the needs of life, models and algorithms are gradually diversifying. The literature can be divided into three categories based on different objective functions; the three categories are as follows:

(1)

Research on green VRP (GVRP)

At present, the main low-carbon policies include emission caps, carbon tax, carbon trading and carbon offset. By conducting the study under these four low-carbon policies, Wang et al. [4] found that the carbon tax policy and carbon trading policy achieved a better emission reduction. Many scholars have studied VRP under the carbon tax policy and carbon trading policy. Yao and Zhang [5] proposed a model for minimizing the total cost under the carbon tax system, and then used an improved genetic algorithm to solve it. G. K. Liu et al. [6] constructed a cold chain path optimization model for joint distribution under the carbon tax policy. Through a quantitative analysis of carbon tax, Ning et al. [7] established a model to minimize integrated costs with the carbon tax as a decision variable. To solve this, they proposed a quantum ant colony algorithm based on adaptive rotation angles. Jiang [8] pointed out that carbon tax policies and carbon trading policies are widely used in the global carbon market. Additionally, the carbon trading policy had a more stable effect of reducing carbon emissions than the carbon tax policy. Cheng et al. [9] demonstrated that by setting different policy parameters, carbon trading policies could guide decision makers to choose transport routes with lower emissions. There are also scholars who introduced carbon emissions directly into the model. Ye et al. [10] used spatial block area planning to construct an adaptive network topology for distribution routes. They confirmed that the distribution routes obtained in this way had lower carbon emissions. Ren et al. [11] built a cold chain path optimization model with the purpose of minimizing carbon emissions. They confirmed that incorporating the upper and lower limits of pheromone concentrations into the traditional ant colony algorithm could improve the optimization efficiency. After considering the congestion of cities in different time zones, Xiao et al. [12] used a genetic algorithm with dynamic programming to formulate delivery vehicles’ routes and departure times to minimize carbon dioxide emissions. In the study of route optimization for multimodal transport, Demir et al. [13] considered the details of the transport process comprehensively and set carbon emissions as one of the objectives of the model. Sadati et al. [14] introduced a multidepot green vehicle routing problem (MDGVRP), where alternative fuel-powered vehicles were used to serve people.

(2)

Research on VRP with customer satisfaction

Due to the perishable and time-sensitive nature of fresh food, service time windows are used to measure customer satisfaction. The degree of customer satisfaction or customer benefits cannot be intuitively indicated in this way. As a result, several scholars have studied satisfaction as a model objective. The model designed by Zulvia et al. [15] for perishable products contains four objectives: minimizing operating costs, minimizing carbon emissions, minimizing cargo damage costs and maximizing service levels, using the fuzzy membership function to represent the service level. They used the multiobjective gradient evolution algorithm (MOGE) to solve the model. According to the characteristics of fresh agricultural products, Wang et al. [16] constructed the penalty cost function based on the time window and used fuzzy logic to calculate satisfaction. Additionally, an improved quantum-behaved particle swarm optimization algorithm (IQPSO) was proposed to solve the problem. Based on delivery time and quality satisfaction, Zhang et al. [17] constructed a time satisfaction function based on a mixed time window and considered the effect of the product deterioration rate on satisfaction. The model was solved using an improved partial elite single-parent genetic algorithm, revealing the changing satisfaction trend with the cost optimization process. Ganji et al. [18] proposed that each customer has a set of time windows with different priorities that represent different values for the customer. Jie et al. [19] used a soft time window to calculate penalty costs. They introduced a mixed-integer nonlinear programming model to minimize the total cost and used three multiobjective metaheuristic algorithms to solve it. Eydi et al. [20] assumed that customers have multiple needs and a time window, prioritizing services with higher demand. Additionally, a multiobjective metaheuristic algorithm was applied to minimize distribution costs and maximize customer satisfaction. Xiao et al. [21] also used a fuzzy membership function to measure customer time satisfaction and established a dual-objective model. They then found fresh food delivery schemes with a nondominated ranking genetic (NSGA-II) algorithm. To minimize energy consumption and maximize customer satisfaction, Ghannadpour et al. [22] categorized customers by priority and used NSGA-II to solve the model.

(3)

Research on time-dependent VRP (TDVRP)

The previous cold chain distribution problem is mainly studied for static road networks, which, in real life, are constantly changing. Malandraki [23] first proposed TDVRP considering the different speeds of vehicles in different periods. Some scholars optimized the delivery route under time-dependent road conditions. Jie et al. [19] introduced stochastic factors (traffic congestion and weather changes) into the traditional particle swarm optimization (PSO) algorithm to improve the velocity calculating formula. A hybrid algorithm (HA) combining the sweep algorithm (SA) and improved particle swarm optimization (IPSO) was developed to minimize the total cost of distribution. Due to the different multiple traffic modes, Asgharizadeh et al. [24] proposed an accurate travel-time estimation model by plotting the speed curves of different types of roads. Xu et al. [25] described the relationship between time and velocity in terms of a trigonometric function and developed a novel multiobjective mixed-integer nonlinear programming (MINLP) model. Additionally, they used the improved NSGA-II to solve it. Similarly, Poonthalir et al. [26] set a minimum speed limit and a maximum speed limit to calculate the expected speed using a triangular distribution. They combined the greedy mutation operator with the particle swarm algorithm to reduce cost and fuel consumption. Considering different speeds in different periods, Zhou et al. [27] used the road segment division method to calculate the speed of vehicles in each road segment. Additionally, they designed an improved ant colony algorithm to find the minimum total cost. Under real traffic work, Z. Liu et al. [28] used a continuous function to describe the relation between velocity and time. In the model, they explored the impact of speed, load and road slope on fuel consumption. To solve the model, a hybrid genetic algorithm with a variational neighborhood search was used.

For the VRP models, solving algorithms can be divided into the precise algorithm, traditional heuristic algorithm and metaheuristic algorithm.

• The precise algorithm is one that can find an optimal solution to a problem. Yu et al. [29] proposed an improved branch-and-price (BAP) algorithm to precisely solve the heterogeneous fleet green vehicle routing problem with time windows (HFGVRPTW). Based on the segmentation of vehicle tours, Lee et al. [30] constructed an extended charging stations network and developed the branch-and-price method on an extended charging station network to solve the problem to optimality. Bruglieri et al. [31] introduced the green vehicle routing problem with capacitated alternative fuel stations (AFSs) and presented two mixed-integer linear programming formulations. Additionally, they described the two variants of the exact cutting planes method for solving the model. The exact solution obtained with the CPLEX solver could be used to test the performance of the heuristic algorithm [32].
• The core idea of traditional heuristic algorithms is to search locally for a better solution than the current one until no better solution is available. Li et al. [33] constructed a linear programming model with the goal of minimizing the weighted sum of the delivery time and transportation cost. The local search process based on a variety of neighborhood structures and the uppercrossing/solution (US) algorithm were adopted to design an improved iterative local search solution algorithm. Considering the travel time variability, Song et al. [34] constructed a nonlinear and concave objective function. They proposed the augmented Lagrangian relaxation method with a block coordinate descent, linearization and Lagrangian substitution techniques. Yang et al. [35] studied the vehicle routing problem with mixed backhauls and time windows (VRPMBTW), which was solved with dynamic programming in a block nonlinear Gauss–Seidel framework.
• The metaheuristic algorithm is an improvement on the heuristic algorithm. It is a combination of the stochastic algorithm and local search algorithm. Schermer et al. [36] formulated the vehicle routing problem with drones and en route operations (VRPDERO) as a mixed-integer linear program (MILP), which solved by an algorithm based on the concepts of the variable neighborhood search (VNS) and tabu search (TS). Wang et al. [37] developed a novel adaptive granularity learning-distributed particle swarm optimization (AGLDPSO) with the help of machine learning techniques, including a clustering analysis based on locality-sensitive hashing (LSH) and adaptive granularity control based on logistic regression (LR). Compared with other large-scale optimization algorithms, this algorithm had a better local search ability.

In order to improve the solving efficiency, many scholars adopted hybrid heuristic algorithms, which combine various metaheuristic algorithms together and integrate their different advantages. Wang et al. [38] combined a variant of the ACO (multiant system) with several local search heuristics to improve the solution quality. Wang et al. [39] constructed a biobjective model to minimize the total cost and delivery vehicles, and designed the tabu search algorithm (TS) and nondominated sorting genetic algorithm (NSGA-II) to solve the problem. These algorithms provided the ideas for solutions for the research in this paper.

Based on the models constructed and algorithms used in previous studies, the previous works and the innovations in this study are summarized as follows:

• From the low-carbon perspective of routing optimization, many scholars have used low-carbon policies in their models, with research mainly focusing on calculating carbon emissions and using a carbon tax to obtain the cost of carbon, while seldom considers carbon trading prices and its spatiotemporal differences. Due to carbon trading policies having a greater advantage in reducing carbon emissions, using carbon trading policies and their uncertainties in studies is conducive to both enriching related research and reducing carbon emissions.
• Most scholars use soft time windows to calculate penalty costs and customer satisfaction in their models. However, these only consider the impact of the current distribution schemes on customers. Other factors such as the influence and importance of customers also have an impact on the development of enterprises. It is more socially relevant to measure the value of customers by integrating various factors. This study calculates customer value from both current value and potential value.
• Under time-dependent road networks, there are few studies that have comprehensively considered the perishability of fresh food, customer value and carbon emissions. It is necessary to increase the research in fresh cold chains. This model can not only effectively reduce the distribution cost of fresh products, but also improve customer satisfaction.

Therefore, a green vehicle routing problem with customer value (CV-GVRP) for fresh product distribution is established. The CV-GVRP is an extension of the traditional GVRP model, which adds time-dependent road network conditions, the characteristics of fresh products and customer value. It is a mixed-integer linear programming model with objectives of minimizing the total cost and maximizing customer satisfaction. According to the NP-hard property of the model, a hybrid heuristic algorithm is used to solve the problem. NSGA-II is extensively applied as a multiobjective optimization algorithm. A greedy algorithm combined with the nondominant sorting genetic algorithm II (G-NSGA-II) was designed to improve the local search ability of the traditional algorithm. Comparisons of the results obtained using models with different objective functions for solving the arithmetic examples confirm the superiority of this model. The conclusions provide a reference for determining decisions in the distribution of fresh cold chains.

The remaining sections of this paper are summarized in sequence: Section 2 of this paper describes the research problem and introduces a mathematical model (CV-GVRP) to solve it. In Section 3, an improved algorithm (G-NSGA-II), combining the nondominated ranking genetic algorithm II with the locally greedy algorithm, is proposed. Section 4 is devoted to computing arithmetic examples via the proposed model and algorithm, with the results confirming their performance. Finally, the conclusions and suggestions for future research are discussed in Section 5.

2. Problem Description and Model Formulation

2.1. Problem Description

Fresh cold chain logistics show that distribution networks consist of multiple customer points and distribution centers. After the goods are packed and sorted in the distribution center, they are distributed using multiple refrigerated trucks. In cases where the vehicle speed changes constantly due to traffic conditions, the optimal vehicle travel path can be planned, considering both the delivery cost and customer value, and the delivery service can be provided within the time window specified by the customer.

In order to ensure the feasibility of the model, this paper assumed the following: (1) That all vehicles depart from the same distribution center and return to the distribution center after completing the delivery. (2) The refrigerated vehicles used are the same type and are sufficient to complete the delivery task. (3) The location of the customer, the amount of demand, the time window and their importance, as well as the temperature of the refrigerated vehicle at the opening and closing of the door, are known. (4) The demand at each customer point does not exceed the maximum capacity of the refrigerated vehicle, and each customer is delivered to using only one vehicle. (5) The quality of the fresh products transported is the same and the requirements for the refrigerated temperature are the same. (6) The preservation period and deterioration rate of fresh products are the same and known. (7) The speed of the vehicle in each period is only related to the delay factor of traffic congestion in that period. (8) Positive publicity is given to the enterprise when customer satisfaction is positive.

2.2. Problem Formulation

In this paper, some corresponding parameters and variables were used to construct the model. The variables and definitions are listed in

Table 1. Symbols and definitions.

Variable	Definition	Parameter	Definition	Set	Definition
$t_i$	The time of arrival at the customer $i$	$Q$	Maximum load of refrigerated truck	$U$	Set of all notes including distribution center
$t_{ij}$	The time it takes to travel from customer $i$ to customer $j$	$v_o$	Average travel speed of refrigerated truck	$U_c$	Set of all customers
$d_{ij}$	The delivery distance between customer point $i$ and $j$	$C_0$	Departure cost of refrigerated truck	$K$	Set of all vehicles
$Q_{ij}$	Load capacity from customer $i$ to customer $j$	$C_p$	Transportation cost per unit distance of refrigerated truck	$U$	Set of all time periods
$s_i$	Time spent serving customer $i$	$\theta$	Rate of corruption
$W_i(t_i)$	Customer satisfaction function	$R_1$	Price per unit of damage
$v_{ij}^p$	The velocity from $i$ to $j$ in a time period $p$	$R_2$	Profit per unit of fresh food
$\omega$	Carbon trading price	$q_i$	Demand of customer $i$
$Z$	The total cost of travel	$n_i$	The number of potential consumers attracted by the positive publicity of the customer $i$
$Z_1$	The cost of vehicle fixed	$(E{T}_i,L{T}_i)$	Customer’s acceptable service time frame
$Z_2$	The cost of transportation	$(E{T}_i^{\alpha },L{T}_i^{\beta })$	Customer’s expected delivery time frame
$Z_3$	The cost of product loss	${\gamma }_1$	The cost per unit of time waiting for customer $i$ to arrive before $E{T}_i$
$Z_4$	The penalty cost	${\gamma }_2$	The cost per unit of time delayed to arrive at customer $i$ after $L{T}_i$
$Z_5$	The cost of carbon emissions	${\lambda }_i$	The weight of customer $i$
$V$	Total customer value	$\phi$	Carbon emissions from refrigeration equipment per unit of goods by delivery unit
$g_i$	The current value of customer $i$	${\rho }_0$	Fuel consumption per unit distance at no load
$p_i$	The potential value of customer $i$	${\rho }^{*}$	Fuel consumption per unit distance at full load
$X_{ij}$	0–1 variable: 1 means delivery performed using vehicle $k$ from customer $i$ to customer $j$, otherwise, it is 0.	$\epsilon$	Carbon emissions per unit of fuel consumption
$Y_{ik}$	0–1 variable: 1 means customer $i$ served using vehicle $k$, otherwise, it is 0.	${\rho }^p$	The delay coefficient of traffic congestion in time period $p$

.

With the previous comprehensive analysis, the CV-GVRP multiobjective optimization model with a minimum total cost and maximum customer value was formulated. The calculations for $Z_1$, $Z_2$, $Z_3$, $Z_4$ and $Z_5$ were explained in the following section. The objective function was as follows:

$$\max V={\displaystyle \sum _{i=1}^u{\lambda }_i(g_i+p_i)}$$

(1)

$$\begin{array}{l}\begin{array}{ccc}& & \end{array}\min Z=Z_1+Z_2+Z_3+Z_4+Z_5\\ \begin{array}{cc}\begin{array}{cc}& \end{array}& =C_0{\displaystyle \sum _{k=1}^k{\displaystyle \sum _{i=1}^u{X}_{0i}^k}}+{\displaystyle \sum _{k=1}^m{\displaystyle \sum _{i=0}^u{\displaystyle \sum _{j=0}^u{d}_{ij}}}}{C}_p{X}_{ij}^k\end{array}\\ \begin{array}{cc}\begin{array}{cc}& \end{array}& +{\displaystyle \sum _{i=0}^u{\displaystyle \sum _{k=1}^m{Y}_{ik}[q_i(1-K_1{e}^{-\theta t_i})+Q_{i0}(1-K_2{e}^{-\theta \frac{q_i}{xh}})]}}{R}_1\end{array}\\ \begin{array}{cc}\begin{array}{cc}& \end{array}& +{\gamma }_1{\displaystyle \sum _{i=1}^u\max \{E{T}_i-t_i,0\}}+{\gamma }_2{\displaystyle \sum _{i=1}^u\max \{t_i-L{T}_i,0\}}\end{array}\\ \begin{array}{cc}\begin{array}{cc}& \end{array}& +\omega {\displaystyle \sum _{k=1}^m{\displaystyle \sum _{i=0}^u{\displaystyle \sum _{j=0}^u{X}_{ij}^k}}}{d}_{ij}\{\epsilon [{\rho }_0+(\rho *-{\rho }_0^{})\frac{Q_{ij}}{Q}]+\phi Q_{ij}\}\end{array}\end{array}$$

(2)

subject to:

$${\displaystyle \sum _{k=1}^m{q}_i}{Y}_{ik}\le Q,i=1,2,3\cdots u,k=1,2,3\cdots m$$

(3)

$${\displaystyle \sum _{k=1}^m{\displaystyle \sum _{j=1}^u{X}_{ij}^k}}\le m,i=0$$

(4)

$${\displaystyle \sum _{k=1}^m{\displaystyle \sum _{j=1}^u{X}_{ij}^k}}=1,i\ne j,i=0,1,2\cdot \cdot \cdot u$$

(5)

$${\displaystyle \sum _{k=1}^m{\displaystyle \sum _{i=0}^u{X}_{ij}^k}}=1,i\ne j,j=1,2,3\cdots u$$

(6)

$${\displaystyle \sum _{k=1}^m{Y}_{ik}}=1,i=1,2,3\cdots m$$

(7)

$${\displaystyle \sum _{j=0}^u{X}_{ij}^k}={\displaystyle \sum _{j=0}^u{X}_{ji}^k}\le 1,i=0,k=1,2,3,\cdots m$$

(8)

$$t_j=t_i+S_i+t_{ij},i=1,2,3\cdots u,j=1,2,3\cdots u$$

(9)

Equations (1) and (2) indicate that the objective functions of the model minimized the total distribution cost and maximized the customer value. Under the constraint of Equation (3), the total order demand assigned to the K-th vehicle could not exceed the loading capacity of the distribution vehicle. Constraint (4) indicated that the number of vehicles involved in the distribution could not exceed the total number of refrigerated vehicles. Each customer point was allowed to be served using only one refrigerated vehicle for one delivery, which was realized with constraints (5) and (6). Constraint (7) indicated that each customer point was served using only one delivery vehicle. Constraint (8) indicated that the refrigerated vehicle departed from the distribution center and returned to the distribution center after completing the transportation task. The time relationship through constraint (9) indicated that the whole distribution process was continuous.

Under a time-dependent road network, the objective function of the CV-GVRP based on the variable price of carbon trading and customer value was to minimize the total cost and maximize the customer value. The total cost contained five components: the vehicle fixed cost, transportation cost, loss cost, penalty cost and carbon emission cost, which were as follows:

2.2.1. Cost Function Composition

(1)

Vehicle Fixed Cost

The fixed cost of a vehicle is the cost of purchase, maintenance, labor paid to the driver, etc. It only depends on the number of vehicles performing the task. The function was as follows:

$$Z_1=C_0{\displaystyle \sum _{k=1}^m{\displaystyle \sum _{i=1}^u{X}_{0i}^k}}$$

(10)

(2)

Transportation Cost

In general, the shorter the distance traveled, the lower the cost. Thus, the transportation cost of the vehicle is proportional to the distance traveled. $C_p$ represents the transportation cost within a unit distance. The total transportation cost $Z_2$ could be expressed as:

$$Z_2={\displaystyle \sum _{k=1}^m{\displaystyle \sum _{i=0}^u{\displaystyle \sum _{j=0}^u{d}_{ij}}}}{C}_p{X}_{ij}^k$$

(11)

(3)

Loss Cost

Both the transit period and the load of goods incurred a loss cost. The continuous lifetime function decayed according to an exponential rate. The decay function of fresh products was: $E(t)=E_0\cdot e^{-\theta t}$ [40]. This formula expresses the quality of the product at different times, where $E_0$ is the product’s initial quality and $\theta$ is the sensitivity factor of the product to time. The loss cost of transit without opening the door was as follows:

$$C_g={\displaystyle \sum _{i=0}^u{\displaystyle \sum _{k=1}^m{Y}_{ik}}}{q}_i[1-K_1{e}^{-\theta (t_i-t_0)}]R_1$$

(12)

Temperature changes caused by opening the door during loading could also affect the freshness of the product. The door opening time was determined using the customer service time $s_i$. The loss cost caused by loading goods while opening the door was:

$$C_l={\displaystyle \sum _{i=0}^u{\displaystyle \sum _{k=1}^m{Y}_{ik}}}{Q}_{i0}(1-K_2{e}^{-\theta S_i})R_1$$

(13)

We could obtain the total loss cost in the whole distribution process as follows:

$$Z_3={\displaystyle \sum _{i=0}^u{\displaystyle \sum _{k=1}^m{Y}_{ik}{R}_1[q_i(1-K_1{e}^{-\theta (t_i-t_0)})+Q_{i0}(1-K_2{e}^{-\theta S_i})]}}$$

(14)

(4)

Penalty Cost

In order to maximize the freshness of our products, we needed to arrive within the time frames set by our customers individually. Otherwise, there would be penalty costs. Suppose that the time window required by the customer was $(E{T}_i,L{T}_i)$. We could use $t_i$ to represent the time when the delivery vehicle arrives at customer $i$. If $E{T}_i\le t_i\le L{T}_i$, no additional penalty cost would be incurred. If the goods were not delivered within the time window, the penalty cost $Z_4$ would be paid. Arriving earlier than the time requested would incur a waiting cost, and arriving later would pay a delayed service fee. The relation was stated as follows:

$$Z_4={\gamma }_1{\displaystyle \sum _{i=1}^u\max \{E{T}_i-t_i,0\}}+{\gamma }_2{\displaystyle \sum _{i=1}^u\max \{t_i-L{T}_i,0\}}$$

(15)

(5)

Carbon Emissions Cost

This cost is mainly related to carbon emissions and the price of carbon trading. Carbon emissions mainly come from two sources: one is generated by energy consumption during transportation, and the other is generated by refrigeration equipment.

• Carbon emissions caused by driving

The amount of CO₂ produced by a vehicle while transporting products is not only related to the distance, but also to the load capacity of the whole process. The fuel consumption per unit distance is different at full load and at no load. The carbon emissions of the travel process from node $i$ to $j$ were given:

$$G_1=\epsilon d_{ij}[{\rho }_0+({\rho }^{*}-{\rho }_0)\frac{Q_{ij}}{Q}]$$

(16)

2.

Carbon emissions from refrigeration

In order to maintain freshness, we used refrigeration equipment in the distribution process. The cumulative carbon emissions from refrigeration between nodes $i$ and $j$ were as follows:

$$G_2=\phi d_{ij}{Q}_{ij}$$

(17)

As a result, the total amount of carbon emitted in the whole distribution process was $(G_1+G_2)$. Due to current carbon trading policies, there are regional variations and significant fluctuations in the price of carbon trading. According to spatial and temporal distribution characteristics of the “carbon K-line”, the CV-GVRP calculated the carbon cost through random variables. The price per carbon transaction was $\omega \sim (\overline{\omega }-\widehat{\omega },\overline{\omega }+\widehat{\omega })$, and the probability of taking any value in this range was equal.

The carbon emission cost could be expressed as:

$$Z_5=\omega {\displaystyle \sum _{k=1}^m{\displaystyle \sum _{i=0}^u{\displaystyle \sum _{j=0}^u{X}_{ij}^k}}}{d}_{ij}\{\epsilon [{\rho }_0+({\rho }^{*}-{\rho }_0^{})\frac{Q_{ij}}{Q}]+\phi Q_{ij}\}$$

(18)

2.2.2. Customer Value Measurement Method

The CV-GVRP considered the customer value in two parts: the current value and the potential value. The current customer value depends on the customer demand, and the two are positively correlated. The calculation method was as follows: $g_i=\frac{q_i}{{\displaystyle \sum _{i=1}^u{q}_i}}{q}_i{R}_2$.

The potential value is primarily related to the company’s reputation, business strength, technological innovation and other factors, while the quality of the distribution has a more pronounced effect on the company’s image. When the distribution is properly optimized, the freshness of the food can improve and customer satisfaction can increase. This motivates more potential consumers and, then, the company gains more potential customer value.

The potential value was calculated using $p_i=W_i{n}_i{q}_i{R}_2$. Customer satisfaction in the CV-GVRP was computed based on the time windows, and the expected time window $\left(E{T}_i^{\alpha },L{T}_i^{\beta }\right)$ was set to be within the specified time window $\left(E{T}_i,L{T}_i\right)$. Figure 1 shows the customer satisfaction function curve of the model.

Based on this, the linear change trend of customer satisfaction over time could be expressed with the fuzzy membership function. The time satisfaction function of a customer $i$ was as follows:

$$W_i(t_i)=\{\begin{array}{l}\frac{(t_i-E{T}_i)}{(E{T}_i^{\alpha }-E{T}_i)}\begin{array}{ccc}\begin{array}{cccc}& & & t_i\in [E{T}_i,E{T}_i^{\alpha }]\end{array}& & \end{array}\\ \frac{(L{T}_i-t_i)}{(L{T}_i-L{T}_i^{\beta })}\begin{array}{cccc}& & & \begin{array}{cc}{t}_i\in [L{T}_i^{\beta },L{T}_i]& \end{array}\end{array}\\ 1\begin{array}{cccc}\begin{array}{cccc}& & & \end{array}& & & \begin{array}{cc}{t}_i\in [E{T}_i^{\alpha },L{T}_i^{\beta }]& \end{array}\end{array}\\ 0\begin{array}{ccccc}& & & & \begin{array}{ccc}& & \begin{array}{cc}{t}_i\notin [E{T}_i^{\alpha },L{T}_i^{\beta }]& \end{array}\end{array}\end{array}\end{array}$$

(19)

The calculation of the value of a single customer could be expressed as ${\lambda }_i(g_i+p_i)$. Giving service priority to customers with high weights ${\lambda }_i$ can increase customer satisfaction and customer value, which is helpful for long-term business growth. The value of total customers was as follows:

$$V={\displaystyle \sum _{i=1}^{u_c}{\lambda }_i(g_i+p_i)}$$

(20)

2.2.3. Time-Dependent Speed Calculation Method

In real life, road conditions are complex and not static. Thus, vehicle speed in each period was calculated according to formula $v_{ij}^p=v_0/{\rho }^p$. The relationship between speed and time is shown in Figure 2.

Ichoua et al. [41] believe that vehicle speed tends to be constant in each subdivided $p$ period. Therefore, the CV-GVRP divided the total time of a day into periods, i.e., $[B_1,F_1]$, $[B_2,F_2]$, $\dots$, $[B_p,F_p]$. The speed in each period remained unchanged; the length of the $p$ period is $H_p$; the driving time in the $p$ period is ; the distance to be driven after the end of $p$ period is $L_{ij}^p$; meanwhile, $d_{ij}={\displaystyle \sum _{p\in P}{d}_{ij}^p}$ and $t_{ij}^p$ represent the driving time in $p$ period. The total driving time was calculated as follows:

Step 1: First, calculate $d_{ij}^p=L_p{v}_{ij}^p$. If $d_{ij}<d_{ij}^p$, calculate $t_{ij}^p=d_{ij}/v_{ij}^p$ next to proceed to Step 3, or calculate $L_{ij}^p=d_{ij}-d_{{}_{ij}}^p$, $t_{ij}^p=L_p$.

Step 2: To calculate the vehicle’s driving time in the remaining period, let $\xi =1$ and the vehicle run in period $p+\xi$, $d_{ij}^{p+\xi }=H_{p+\xi }{v}_{{}_{ij}}^{p+\xi }$; if $d_{ij}^{p+\xi }\le L_{ij}^p$, $t_{ij}^{p+\xi }=H_{p+\xi }$, $L_{ij}^{p+\xi }=L_{ij}^p-d_{ij}^{p+\xi }$ and $\xi =\xi +1$, then repeat this step. Otherwise, calculate $t_{ij}^{p+\xi }=d_{ij}^{p+\xi -1}/v_{ij}^{p+\xi }$ to proceed to Step 3.

Step 3: The total vehicle travel time $t_{ij}={\displaystyle \sum _{p\in P}{t}_{ij}^p}$ is obtained and the section $(i,j)$ driving time computation is finished.

3. G-NSGA-Ⅱ Algorithm

3.1. Description of the G-NSGA-II Algorithm

CV-GVRP is a multiobjective model with a comprehensive consideration of enterprise cost and customer value. The nondominated ranking genetic algorithm-II with the elite genetic strategy is widely used in multiobjective optimization and solves the problem exactly [20,39]. The Pareto solution can intuitively reflect the relationship between the two objectives. Using the greedy algorithm to obtain an initial solution can expand the local search and improve the solution quality [32]. To further improve the global search ability and convergence speed of the algorithm, this paper designed an improved NSGA-II (G-NSGA-II) with the characteristics of the local optimization of the greedy algorithm to solve the model.

• Coding Mechanism

Natural number coding was used to visually describe possible solutions to the problem. A chromosome denotes a feasible distribution solution using client sites as genes of the chromosome, and 0 denotes a distribution center. The chromosome can be expressed as (0, I1, I2..., Ir, 0), which means a vehicle departs from the distribution center (point 0) and services customers I1, I2..., Ir in order. The vehicle then returns to point 0. A schematic of chromosome coding is shown in Figure 3.

2.

Population Initialization

Since the greedy algorithm could give an initial solution in a shorter time, this paper combined the greedy algorithm to solve the model through the following steps:

Step 1: First, set the initial path to empty.

Step 2: From the initial node (warehouse), randomly select customer point $i$ and calculate the distance between customer point $i$ and neighboring points. Then, add the customer point $j$ with the shortest delivery distance to the path.

Step 3: Determine whether there are still unserved customer nodes; if so, proceed to the next step, otherwise, terminate the algorithm.

Step 4: Calculate the distance between customer point $i-new$ and the remaining customer points, choose the shortest one and then return to Step3.

3.

Improve the Sorting Fitness Strategy

Figure 4 reflects the density information of individuals. The Pareto ranking of individuals 1, 2, 3, 4, 5 and 6 was one, and the ranking of individuals a, b and c was two. However, the density around individual a was significantly larger than that of b and c. Therefore, to avoid them having the same probability of entering the next generation, the ranking values were combined with the density information around the individuals to distinguish individuals in the same layer. This could improve the diversity of the population distribution without increasing the complexity of the algorithm.

4.

Crowding Calculation

After several iterations to obtain a hierarchy of noninferior solutions, the crowding distance between neighboring individuals was calculated using objective values. The individuals with a more considerable crowding distance in the same layer were retained to ensure the diversity of solutions. $D{(i+1)}^{fm}$ and $D{(i-1)}^{fm}$ represent the m-th objective function value of individuals $(i+1)$ and $(i-1)$, respectively, and $fm$ represents the m-th function in the model. $f{m}^{\max }$ and $f{m}^{\min }$ are the maximum and minimum values of the m-th objective function values, respectively. The following formula could calculate the crowding degree of individuals:

$$D(i)=D(i)+\frac{D{(i+1)}^{fm}-D{(i-1)}^{fm}}{f{m}^{\max }-f{m}^{\min }}$$

(21)

5.

Crossover and Variation

In this paper, we adopted the natural number encoding method and chose the partial-mapped crossover (PMX) for the crossover operation in order to improve the convergence speed of the algorithm. The basic procedure was as follows: 1. Two individuals were randomly identified as crossing nodes in the parent population. 2. The genes of the two crossover nodes were exchanged between the parents. 3. Relationships were mapped based on genes between two crossover nodes, replacing duplicated genes on the same parent to obtain two offspring. The diagram of the PMX is shown in Figure 5.

This paper adopted the exchange variation to increase the diversity of individuals in the population and improve the local search ability of the algorithm, in which any two gene points were selected to be swapped. The result of the mutation is shown in Figure 6.

3.2. Algorithm Process

Figure 7 represents the whole calculation process. The detailed procedure was as follows:

Step 1: To determine the population size, maximum number of iterations, crossover probability and mutation probability, use the greedy algorithm to generate the initial population.

Step 1.1–Step 1.3 are the processes of generating the initial solution:

Step 1.1: Calculate the distance between two customer points and sort them from smallest to largest.

Step 1.2: Determine whether the solution satisfies the subpath condition. If it does, add it to the current path or determine the next path.

• The path is not closed after being added.
• After adding, no node should exceed two connecting edges.
• The total demand of the customer should not exceed the load of refrigeration.
• It can meet the requirements of the customer’s time window.

Step 1.3: Execute the previous step until all subpaths are assigned while both endpoints of each path are connected to the distribution center to form a closed loop.

Step 2: Perform an improved nondominated ranking based on the biobjective function values of the individuals in the initial population. Then, apply a tournament strategy to select and generate new offspring populations through a cyclic crossover and mutation.

Step 3: Combine the offspring and parent populations to generate a new population. Through the elite strategy, compare the improved nondominated ranking value and crowding distance to obtain a better combination of individuals to generate the parent population.

Step 4: Iteratively update the newly generated parent population with genetic manipulation.

Step 5: Judge whether the current iteration number reaches a maximum and output the final result. Otherwise, repeat Step3.

4. Experiments and Results Analysis

In order to further test the performance of the algorithm, the algorithm was compared with the traditional NSGA-II. Experiments were carried out with three classical datasets (Solomon benchmark datasets r101, c101 and rc101 [42]) to test the performance of the algorithms. Customer influence, information dissemination and customer weight information were added to form the small-scale cases (c101_25, r101_25 and rc101_25), the medium-scale cases (c101_50, r101_50 and rc101_50) and the large-scale cases (c101_100, r101_100 and rc101_100), which had 100 customers.

Table 2. Related parameters and values.

Parameter	Value	Parameter	Value
The departure cost of the vehicle $C_0$$C_0$ (RMB)	100	Constant value of product deterioration $K_1,K_2$	(1,1)
The cost of running a vehicle $C_p$ (RMB/Km)	2	Customer weight $\lambda$	(1~5)
Spoilage rate of fresh produce $\theta$	0.002	Affect the size	(0~25)
Maximum load of vehicle $Q$ (kg)	500	Customer current value threshold	150
The cost of waiting to arrive early (RMB/h)	2	Customer potential value threshold	20
The penalty cost of late arrival (RMB/h)	3	The price per unit of fresh produce $R_1$ (RMB)	30
The average speed $v_0$ (km/h)	50	Profit per unit of fresh $R_2$ (RMB)	30
Fuel consumption per unit distance when fully loaded ${\rho }^{*}$ (L/Km)	0.377	Fuel consumption per unit distance when no load ${\rho }_0$ (L/Km)	0.165
CO2 emission coefficient $\epsilon$ (Kg/L)	2.63	The amount of CO2 produced by refrigeration per unit distance traveled by goods delivered $\phi$ g/(Kg·Km)	0.0066
Benchmark carbon trading price $\overline{\omega }$	30	Range of carbon trading price changes $\widehat{\omega }$	(0~30)

shows the values for the model parameter settings. The values of the parameters were obtained from [5,40].

Table 3. Traffic congestion factor.

Period of Time	[0,5]	[5,7]	[7,9]	[9,12]	[12,14]	[14,18]	[18,20]	[20,24]
Congestion Delay Coefficient	1	1.2	1.4	1.6	1.5	1.3	1.7	1.2

represents the traffic congestion coefficients for different periods. The parameter values in both tables were from the study by Wang et al. [40]. All case tests were performed on an Intel(R) Core (TM) i5-7200U, 2.50 GHz, 4 GB RAM, Windows 10 (64-bit) computer, and the algorithms used in this paper were programmed on MATLAB R2019b.

4.1. Model Solution

There were six customer information points in

Table 4. Point-related customer information.

No.	0	1	2	3	4	5
Location	(35,35)	(41,49)	(35,17)	(55,45)	(55,20)	(15,30)
Demand	0	10	7	13	19	26
Acceptable time window	[0.0,24]	[1.8,4.8]	[0.0,2.5]	[0.0,2.0]	[0.5,3.0]	[1.0,4.0]
Optimum time window	[0.0,24]	[2.8,3.8]	[0.5,1.5]	[0.5,1.5]	[1.0,2.0]	[2.0,3.0]
Service time	0	0.3	0.3	0.3	0.3	0.3
Customer importance	1	1	1	1	1.2	1.5
Weight	1	1	1	3	2	4
Scope of influence	0	13	20	15	2	5

. This study used customer point information to solve the CV-GVRP under the G-NSGA-II algorithm, and compared the obtained results with those under the uniparental genetic algorithm. Although the model and algorithm in this paper did not have an advantage in terms of time consumption, they performed better in reducing the total cost and increasing the customer value. The solution results are shown in

Table 5. SAPGA and G-NSGA-II solution results.

Algorithm	Total Cost	Customer Value	Computation Time	Distribution Path
SAPGA	375.3001	378.9227	4.62	0-3-2-4-0 0-5-1-0
G-NSGA-II	322.2635	442.8125	13.44	0-1-0 0-3-4-2-0 0-5-0

.

4.2. Numerical Benchmark

4.2.1. Comparison between MOPSO, NSGA-II and G-NSGA-II

Using MOPSO, the conventional NSGA-II and G-NSGA-II to solve the nine cases, the parameter settings of the algorithms used in the experiments were as follows: Figure 8, Figure 9 and Figure 10 show the r101_100, r101_50 and r101_100 Pareto solution results. In this study, the algorithm parameters were set as per the study by Rabiee et al. [43]. The basic parameters of the experiment are shown in

Table 6. Parameter settings of each algorithm.

MOPSO		NSGA-II		G-NSGA-II
Parameter	Value	Parameter	Value	Parameter	Value
The population number	200	The population number	30	The population number	30
Size of the warehouse	100	Crossover probability	0.9	Crossover probability	0.9
Inertial factor	0.9	Mutation probability	0.1	Mutation probability	0.1
Local velocity factor	1	The number of iterations	500	The number of iterations	500
Global velocity factor	2
The number of iterations	300

.

The algorithm comparison continued using these nine examples, and the following three indicators were used to compare the performance of the algorithms:

• The number of Pareto solutions (NPS) indicated an algorithm’s total number of nondominated solutions.
• The diversity metric (DM) of solutions was measured with $\sqrt{{(\max f_{1i}-\min f_{1i})}^2+{(\max f_{2i}-\min f_{2i})}^2}$. It determined the diversity of nondominant solutions obtained with each algorithm.
• The formula $({\displaystyle \sum _{i=1}^n\sqrt{f_{1i}{}^2+f_{2i}{}^2}})/n$ was used to calculate the average distance between the Pareto solution and the ideal point (0, 0), which was the mean ideal distance (MID).

For the first two indicators, the higher the value, the better the algorithm performance, while the lower the average ideal distance value, the better the algorithm performance [43]. The comparison results are shown in

Table 7. Algorithm performance comparison.

Examples	NPS			DM			MID
	G-NSGA-II	NSGA-II	MOPSO	G-NSGA-II	NSGA-II	MOPSO	G-NSGA-II	NSGA-II	MOPSO
c101_25	16	15	14	1116.3	889.9	857.3	1400.1	1463.8	1422.7
r101_25	14	13	13	707.4	702.0	645.6	951.4	1150.5	1168.0
rc101_25	15	14	14	455.0	445.1	412.7	1596.3	1582.4	1579.4
c101_50	15	14	14	1352.9	1284.3	1036.5	3112.3	3213.8	3203.9
r101_50	14	12	12	1349.9	1055.2	922.3	1024.6	1447.0	1583.7
rc101_50	15	14	13	1047.1	1049.1	934.9	1907.0	2131.3	2383.2
c101_100	16	14	13	1936.4	1727.5	1504.3	5401.0	6403.8	7799.2
r101_100	15	13	13	2746.8	2469.3	2610.8	4015.9	4412.5	5139.5
rc101_100	16	15	13	1898.8	1727.9	1558.5	5045.8	5144.0	6260.4

.

In Figure 11, the NPS of G-NSGA-II and NSGA-II was significantly greater than that of MOPSO. The first two were similar in terms of dispersion. However, the median NPS of G-NSGA-II was significantly higher than that of NSGA-II. Figure 12, showing the DM, shows the boxplot of G-NSGA-II, which outperformed the others in dispersion and median. The MID in Figure 13 shows that G-NSGA-II and NSGA-II performed significantly better than MOPSO, with G-NSGA-II having a slightly lower median. Based on the results as previously seen, G-NSGA-II performed better. In summary, the G-NSGA-II algorithm outperformed the traditional NSGA-II and MOPSO ones in terms of convergence and diversity. Its Pareto solutions were more uniformly distributed than those of the traditional NSGA-II. Therefore, it had better results in solving the problem.

4.2.2. Analysis of Carbon Trading Price Sensitivity

Through the analysis of different types of data, it was found that the total cost changed significantly under the influence of variable carbon trading prices, especially in the large-scale examples. The change in price was not a single trend, which was why the carbon trading prices differed from the carbon tax policies. Under its influence, carbon emissions could increase or decrease, so the total cost also fluctuated. Due to firms gaining from selling carbon allowances, the cost could be reduced. Additionally, this behavior could incentivize firms to continue cutting carbon.

The total cost was lower than the situation with a fixed carbon trading price. The changes are shown in Figure 14. Furthermore, it indicated that a proper adjustment in carbon trading prices could motivate enterprises to adjust their distribution schemes. This could reduce the total cost and contribute to declining social carbon emissions from a macroperspective.

4.2.3. Influence of Time-Dependent Network on Routing Strategy

In cases of different scales, the model in this paper had different improvements. For example, in the small-scale case, the total cost of the CV-GVRP was saved by 6.42%, and the customer value was increased by 15.43% on average. In the medium-scale case, the total cost of the CV-GVRP was saved by 5.4%, and the customer value was increased by 17.68% on average. In the large-scale case, the total cost of the CV-GVRP was saved by 3.75%, and the customer value was increased by 17.40% on average. The results of the comparison are shown in

Table 8. Comparison of static road network and CV-GVRP results.

	Static Road Network		CV-GVRP
	Total Cost	Customer Value	Total Cost	Customer Value
c101_25	1674.5	1271.9	1716.3	1423.1
r101_25	1139.1	907.1	1231.8	1091.9
rc101_25	1197.6	1495.3	1301.0	1705.2
c101_50	3415.3	2115.6	3533.7	2714.9
r101_50	2687.3	2380.7	2893.7	2736.1
rc101_50	2435.0	2512.8	2558.0	2758.3
c101_100	6926.6	3770.9	7299.9	4211.9
r101_100	5171.1	3404.9	5379.1	4106.5
rc101_100	5791.7	5389.1	5898.9	6462.3

. The cost savings of the CV-GVRP were more significant in the small-scale cases, and the increase in customer value was more obvious in the medium-scale and large-scale cases. The average increase in the total cost was 5.192%, and the average increase in customer value was 16.838%, which indicated that the CV-GVRP performed better than the model with a static road network in improving customer value. Considering traffic congestion could increase the delivery time and increase the total cost, but the CV-GVRP brought more customer value, which is conducive to the long-term development of the enterprise.

4.2.4. Influence of Customer Value on Routing Strategy

The GVRP is a model that does not consider the potential value of customers. It has the same constraints as the CV-GVRP. However, the objective function of the GVRP is customers’ total cost and current value, which does not consider the impact of potential value on corporate reputation. Then, the G-NSGA-II was used to solve the GVRP for different cases. In the small-scale case of 25 customers, the total cost of CV-GVRP increased by 7%, the customer value increased by 14% and the satisfaction increased by 26% on average. In the medium-scale case of 50 customers, the total cost increased by 8%, the customer value increased by 14% and the satisfaction increased by 15% on average. In the large-scale case of 100 customers, the total cost increased by 3%, the customer value increased by 21% and the satisfaction increased by 13% on average. Overall, the average increase in the total cost was 3%, the average increase in customer value was 21% and the average increase in satisfaction was 13%. The increase in the total costs was most likely due to the increased delivery distance and time required to serve important customers. Additionally, the improvement in customer satisfaction brought more potential value to the enterprise and increased the total customer value. The results of the comparison are shown in

Table 9. Comparison of GVRP and CV-GVRP results.

	GVRP			CV-GVRP
	Total Cost	Current Value	Satisfaction	Total Cost	Customer Value	Satisfaction
c101_25	1598.6	1236.3	0.783	1716.3	1423.1	0.963
r101_25	1182.6	939.1	0.726	1231.8	1091.9	0.931
rc101_25	1211.9	1539.1	0.744	1301.0	1705.2	0.937
c101_50	3289.4	2418.7	0.791	3533.7	2714.9	0.916
r101_50	2734.5	2417.2	0.811	2993.7	2736.1	0.927
rc101_50	2373.9	2389.9	0.803	2558.0	2758.3	0.914
c101_100	6999.8	3417.2	0.822	7299.9	4211.9	0.907
r101_100	5266.5	3396.7	0.807	5379.1	4106.5	0.912
rc101_100	5713.6	5467.5	0.793	5898.9	6462.3	0.908

.

4.3. An instance of Distribution in Shanghai City

A refrigerated logistics company in Shanghai has to deliver to 20 distribution points. The information about customers is shown in

Table 10. Shanghai Huangpu district distribution point information.

X-axis (km)	Y-axis (km)	Time Window (min)	Acceptable Time Window (min)	Demand (t)	Service Time (min)	Importance	Weights	Influences Scope
353.914	3455.623	0–500	0–500	0	0	0	0	0
353.721	3456.089	30–90	30–120	0.80	10	1.2	2	10
354.020	3456.869	60–90	60–120	3.35	20	1.3	2	3
355.047	3456.028	30–90	30–150	2.95	15	1.0	1	2
355.179	3456.943	60–90	60–120	2.40	10	1.0	1	7
355.135	3453.781	70–130	70–160	2.75	20	1.0	1	5
355.351	3453.055	80–110	80–140	3.30	20	1.1	1	8
356.423	3455.002	70–130	70–160	2.80	25	1.2	2	3
356.849	3454.232	60–90	60–120	3.25	20	1.3	3	2
355.197	3455.186	80–140	80–180	2.15	15	1.4	4	4
355.870	3455.585	85–115	85–150	3.05	15	1.4	4	6
354.804	3454.310	80–140	80–170	3.20	25	1.5	5	7
355.729	3454.642	70–100	70–130	3.50	10	1.2	4	8
354.597	3452.140	60–120	60–150	0.55	20	1.3	3	10
354.036	3452.706	80–110	80–150	2.70	20	1.0	2	12
353.423	3454.06	80–140	80–170	1.70	15	1.0	2	13
354.337	3453.427	90–120	90–160	2.25	20	1.0	1	7
354.358	3454.64	90–150	90–180	2.75	25	1.0	2	10
355.082	3453.816	70–100	70–140	1.90	10	1.5	5	15
353.345	3454.992	60–120	60–150	3.15	20	1.3	3	13
353.180	3453.517	80–110	80–140	1.00	15	1.0	2	7

. Assuming that the traffic congestion in the distribution is not considered, the vehicle is driven at a uniform speed of 40km/h and the carbon trading price is fixed at 54 RMB/ton. Moreover, the potential value of the customer is not considered. The solution can be solved with G-NSGA-II and then compared with the results of CV-GVRP. This paper’s model considered the variable speed, variable carbon trading price and customer value, resulting in an 8.18% reduction in the total cost, a 23.82% increase in customer value and a 4.05% reduction in time. The results are shown in

Table 11. Example solution result comparison.

	Pre-Optimization	Optimized
Total cost	1858.5	1706.5
Customer value	1730.0	2142.1
Calculating time	19.6442	18.8483
Delivery route	0→3→1→0 0→9→10→17→0 0→4→2→15→0 0→8→12→14→16→0 0→19→11→18→6→13→0 0→7→5→20→0	0→9→18→19→5→17→0 0→13→1→20→14→15→0 0→2→3→12→6→16→11→0 0→4→8→10→7→0

. In addition, the number of vehicles used was two less. The two solutions are shown in Figure 15 and Figure 16.

4.4. Discussion on the Model and the Algorithm

In general, the experiments focused on two aspects: the algorithm and the model.

To analyze the algorithm performance, six customer points were used as examples to verify the feasibility of the G-NSGA-II algorithm. To further explore the superiority of the G-NSGA-II algorithm, the small–medium–large-scale test sets were constructed using Solomon data. The solution results of the MOPSO, NSGA-II and G-NSGA-II algorithms were compared in convergence and diversity. By plotting boxplots in NPS, DM and MID and comparing their median and dispersion, the G-NSGA-II algorithm outperformed the others.

The CV model was used to analyze the influence of carbon trading uncertainty, time-varying network and customer potential value on the distribution schemes.

• Variable carbon trading prices were used to reflect the temporal and spatial aspects of carbon trading policies. The experimental results showed that there was no linear relationship between the range of carbon trading price and the total cost, and the reasonable setting of the carbon trading price helped to reduce the cost and achieve a green distribution.
• The experiments were designed to explore the influence of time-dependent road networks on distribution strategies, and it was found that the total costs, customer values and customer satisfaction would be affected on different scales. For large customer groups, the increase in customer value was more significant.
• This paper studied the impact of potential value from customers on the distribution strategy. The experimental results showed that considering the potential value of customers could greatly improve customer satisfaction, especially in small groups.

Finally, the CV-GVRP model and the G-NSGA-II algorithm were applied to practical problems, and a planning scheme was proposed for real-life cold chain distributions. Constructing a cold chain distribution model under a carbon trading policy was conducive to reducing carbon emissions and achieving sustainable developments in enterprises. Taking into account the time window and potential value of the customer could help improve the reputation and benefits of the business by providing timely services to important customers.

5. Conclusions

From the perspective of variable carbon trading prices, this paper built a multiobjective optimization model, the CV-GVRP, which aimed to minimize the total cost and maximize customer value, considering both customer value and the time-dependent network. This model was an improvement of the traditional cold chain one, which was mainly applicable to fresh agricultural products with a short shelf life, low temperature and simple requirements for distribution conditions, such as leafy vegetables and fruits. By setting different ranges of carbon trading prices, it was found that flexible carbon trading prices positively affected companies in reducing their carbon emissions. Compared to static networks, the CV-GVRP was found to be advantageous in the total cost savings and customer value growth. Moreover, the CV-GVRP could improve customer satisfaction at a lower cost than models that do not take customer value into account. Small, medium and large-scale examples of c, r and rc, respectively, were constructed using Solomon data. Using the nine examples, G-NSGA-II, NSGA-II and MOPSO were used to solve this model. The values of NPS, DM and MID showed that the G-NSGA-II algorithm performed better than the others.

Future research can incorporate additional factors into the consideration of customer value, such as service quality and product quality. For the calculation of variable speed, only the congestion case was considered. However, the speed limits, emergencies and traffic control also affect vehicle speed in real life. This study considered only the carbon trading policy; in the future, the carbon tax policy and carbon trading policy can be used in coordination to reduce carbon emissions. To better address practical problems, there is still a lot of space to enrich the details of this model.