# Optimization of VRR for Cold Chain with Minimum Loss Based on Actual Traffic Conditions

Authors: Lishuan Hu; Caihong Xiang; Chengming Qi
Journal: Wireless Communications and Mobile Computing
 (2021)

---

## Abstract

Recently, fresh agricultural cold-chain logistics have been greatly developed with the increasing needs of people’s life. Reducing costs of cold-chain distribution has become the main object of loss control in logistics enterprises. The objective of this research is to find a set of optimal routes that minimize the total loss, including fuel cost, refrigeration cost, soft time window penalty cost, and cargo damage cost over transit time. In this paper, the definition and model construction of vehicle routing problem (VRP) with multiobjective minimum lost are introduced first. Then, an ant colony optimization (ACO) algorithm combined with Pareto local search (PLS) is put forward to solve the minimum loss model. In order to avoid the influence of complex road conditions during distribution, the distance matrix and the transit time matrix are both derived from the recommended navigation road based on E-map API. At last, a compare experiment between the traditional method and our proposed method is performed. The results indicate that our method has strong applicability and potential advantages in cold-chain logistic and has important practical significance and application value.


---

## Body

## 1. Introduction

According to relevant media reports, rotting fruit, vegetables, and other foods transported by truck alone are worth about 70 billion yuan each year, causing a huge economic waste. Agricultural products especially for fruit, vegetables, fish, meat, etc. require strictly limited temperature, humidity, and time in the process of transportation and storage. As a branch of the logistics industry, cold-chain logistics provides a guarantee for the safe transportation of fresh agricultural products. In the process of transportation from cold storage to customers, a complete cold-chain logistics realizes the temperature control of the whole process of refrigerated and frozen food, as well as the closed environment, storage, and transportation during loading and unloading of goods. By selecting the optimal path in the process of cold-chain logistics transmission, the circulating rate of fruits and vegetables, meat, and aquatic products is reduced; the waste of resources and the cost of logistics costs are reduced too. The vehicle routing problem (VRP) introduced by Dantzig and Ramser [1] plays a central role in the optimization of distribution networks.Recently, the loss problem of cold-chain distribution has attracted the attention of many scholars and experts. For distribution problem of fresh agricultural products, fuel cost, refrigeration cost, cargo damage cost, and soft time window penalty cost have become the important components of loss control in distribution companies [2]. Li et al. [3] proposed a low-carbon model for fresh food and a genetic simulated annealing algorithm to solve the model in the cold-chain distribution. Yao et al. [4] proposed a minimizing fuel consumption solution to time-dependent VRP with time window. Kim et al. [5] proposed a Markov decision process method to solve a dynamic VRP (DVRP) model with nonstationary transit times under actual traffic congestion. Chen et al. [6] developed a hybrid heuristic algorithm including harmony search and neighborhood descent to solve DVRP with time window. Abidi et al. [7] proposed a GA with a simple heuristic to solve a variant of RVPR (Rich VPR) with time widows and dynamically changing orders. Dongdong and Yinzhen [8] proposed a green multitype VRP with time windows to reduce the wastes of fuel consumption and carbon emission and use an improved tabu search algorithm to solve the G-MVRPTW. Fan et al. [9] pointed out that in fresh agricultural product cold-chain logistics, the total costs are composed of five kinds: fixed, transportation, damage, penalty, and energy consumption. Fang and Ai [10] proposed a mathematical model with soft time window penalty cost, refrigeration cost, cargo damage cost, and a hybrid ant colony algorithm to minimize the total costs. The DVRP mainly deal with the time-variation information of customer demands and road conditions. To share the traffic information and classify traffic conditions, Big Data and classification techniques are used in logistics and transportation. If the drivers receive the information of the traffic congestion or poor weather, they can change their way to reduce the costs and save time. Dimensionality reduction must be performed first in Big Data transmitted and stored process. Thippa et al. [11] proposed a machine learning (ML) algorithm with PCA to reduce the dimension of Big Data when the data sets are high. Gadekallu et al. [12] have studied the hybrid PCA-whale optimization algorithm to extract features and used a deep neural network to classify the diseases of tomato. Guha et al. [13] proposed an ANN-based content classification in combination with n-grame TF-IDF feature descriptor to classify the documents with accurate, sensitive information. The ant colony optimization (ACO) proposed by Colorni et al. [14] has been widely used in solving NP-hard vehicle routing problems. Many scholars used ACO algorithms to solve multiobjective combinatorial optimization problems [15–19]. These papers considered the costs of fuel, refrigeration, cargo damage, and delay simultaneously. However, their studies are not comprehensive enough in factors affecting the loss, such as the cost of energy consumption and rotting consumption when the door of compartment is opening.We construct a total cost model, including soft time window penalty cost, fuel cost, refrigeration cost in transit and during unloading, and cargo damage cost in transit and during unloading. In our model, to take into account the actual traffic condition, the distance matrix and the transit time matrix are both archived from the navigation functions based on E-Map API. The traditional heuristic algorithm ACO has some limitations, such as easy stagnation in the initial stage and slow search speed. We perform an ACO algorithm with Pareto local search (PLS) on ants to obtain the uniform Pareto-optimal frontier and keep the diversity of Pareto solution set.The paper is organized as follows. Section1 introduces the definition and structure of VRP for cold-chain logistics. Section 2 constructs a minimum cost model based on actual traffic conditions. Then, an improved ACO algorithm with PLS is presented in Section 3. In Section 4, a compare experiment is performed to indicate our proposed model and the improved method. At last, we make a conclusion of our contributions in Section 5.


## 2. Problem Description

A cold-chain supplier has a warehouse with a certain number of transportation vehicles that delivers a variety of fresh products to a certain number of customers. The capacity of the transportation vehicles is limited, and the vehicles are all the same type with a total capacity at the starting point. The object of vehicle routing problem is to distribute goods to each customer with correct goods, limited time, minimum cost, and so on. VRP with time window means that a vehicle has to visit a customer within a certain time window. So, the service time, the transit time, and the total time must be calculated when dealing with the VRP. The mathematical model of VRPTW in cold-chain distribution is defined as follows [20].There areK transportation trucks in this distribution center. The maximum load and the maximum transit time of per vehicle are Q and T. The kth distribution vehicle is responsible for path k. N is the total customers. The needs and the service time of customer point i are qi and ui; tij means the delivery time from customer points i to j. xik is a 0-1 variable, and xik=1 means the point i in the path k. xijk is a 0-1 variable, and xijk=1 means that the kth distribution vehicle travels from i to j.


## 3. Model

Minimizing the loss cost is the target of this paper. Therefore, the following paper introduces the structure of the loss cost of cold chain.


### 3.1. Distance and Transit Time Analysis

In traditional methods, the distance between two points is achieved by Euclidean distance formula and the transit time is achieved by the distance divided by assumed average speed. But, the actual transit time is influenced greatly by road congestion, traffic accidents, traffic control, raining, snowing, etc. Collecting real-time traffic data will result in large amounts of data, and Big Data analysis technology should be considered.In our literature, the navigation function of E-map is used to avoid the mentioned shortcomings. The distance matrix and the transit time matrix under actual traffic conditions among customer points (including the distribution center) are derived from the E-map public platform [21].As Figure1 shows, there are two points with their latitudes and longitudes on the E-map.Figure 1
Two points on the E-map.As Figure2 shows, the distance and the transit time between the two points are shown by the navigation function of the E-map platform.Figure 2
The output of navigation result.In our study,tij represents transit time between customer i and customer j, and dij is the distance between customer i and customer j. As abovementioned, tij and dij are both achieved from the navigation result based on E-map API.


### 3.2. Loss Cost Target Construction

#### 3.2.1. Soft Time Window Penalty Cost

Each customer has different numbers of good needs and a soft time window. Penalty cost will be incurred if the vehicle arrives beyond the time window boundary. The service time of each customer point isuii=1,2,⋯,N, Eti,Lti is the range of the time window of customer point i, C11 is the penalty coefficient when vehicle ki arrives at customer i earlier than Eti, and C12 is the penalty coefficient when vehicle ki leaves from customer i later than Lti. The arrival time of vehicle k at customer i is Tki. According to the analysis above, the penalty cost of cold chain can be defined as the following:
(1)C1=C11∑k=1K∑i=1NxikmaxEti−Tki,0+C12∑k=1K∑i=1NxikmaxTki+ui−Lti,0,s.t.Tki≥Tk0+t0i,Tki+1≥Tki+tii+1+ui.


#### 3.2.2. Fuel Cost

The fuel consumption of the distribution vehicle is inevitable for completing the distribution task. We assume that the fuel cost is proportional to the load of the vehicle pert unit distance. LetQk be the load of the k vehicle and Pk be the load rate of the kth vehicle.
(2)Qk=∑i=1nxikqi,Pk=QkQ.SetsC21 and C22 represent no load and full load fuel consumption costs per unit distance of the distribution vehicle. C2k means the fuel consumption cost per unit distance of the kth vehicle.
(3)C2k=C21+pkC22−C21.The total cost of fuel consumptionC2 can be expressed as
(4)C2=∑k=1K∑i=1N∑j=1NxijkdijC21+PkC22−C21.


#### 3.2.3. Refrigeration Cost

The temperature of the cold-chain logistic must be maintained at a certain low level to keep the freshness of goods. The refrigeration function of vehicles must be operated to achieve the required temperature, which causes the refrigeration cost immediately.(1)
Refrigeration cost in transitThe refrigeration cost generated by maintaining the required low temperature per unit time in transit isC31; the total cooling cost C31 in transit can be expressed as
(5)C31=∑k=1K∑i=1N∑j=1NxijktijC31.(2)
Refrigeration cost during unloadingOpening the door of the compartment will cause the cool air inside to flow out and the hot air from outside to flow in. To maintain the temperature inside, more energy will be consumed.C32 is set to be the cooling cost per unit time during opening the door, and the total refrigeration cost of unloading C32 is
(6)C32=∑k=1K∑i=1NxikuiC32.C3 represents the total cooling cost, which consists of the total cooling cost in transit, and the total cooling cost of unloading can be expressed as
(7)C3=∑k=1K∑i=1N∑j=1NxijktijC31+xikuiC32.


#### 3.2.4. Cargo Damage Cost

The most cargo of the cold chain is fresh goods; with the increase of transit time even in low temperature, the growth of microorganisms will happen. When opening the door of compartment during unloading process, the temperature inside will be unstable which results in the damage of fresh goods more significantly. As mentioned above,p is the unit price of goods, and Qki is the load of the kth vehicle when it leaves from point i.
(1)
Cost of cargo damage in transitC41 is given as the loss coefficient of goods per unit weight per unit time in transit; the total cost of cargo damage in transit is as follows:
(8)C41=∑k=1K∑i=1N∑j=1NxijktijpQkiC41.(2)
Cost of cargo damage during unloadingDuring unloading at a customer point, the door of the compartment is open and the cost of cargo damage will increase obviously. Let the cargo damage coefficient of goods per unit weight per unit time when unloading beC42; the cost of cargo damage during unloading can be expressed as
(9)C42=∑k=1K∑i=1NxikuipQkiC42.The total cost of damage is(10)C4=pQki∑k=1K∑i=1N∑j=1NxijktijC41+xikuiC42.In conclusion, the mathematical model of minimum loss target can be expressed as(11)minZ=C1+C2+C3+C4.


### 3.3. Mathematical Model of Minimum Loss Analysis

According to the above analysis, the minimum loss model of the cold-chain distribution problem is as follows:(12)minZ=C11∑k=1K∑i=1NxikmaxEti−Tki,0+C12∑k=1K∑i=1NxikmaxTki+ui−Lti,0+∑k=1K∑i=1N∑j=1NxijkC21+PkC22−C21+∑k=1K∑i=1N∑j=1NxijktijC31+xikuiC32pQki∑k=1K∑i=1N∑j=1NxijktijC41+xikuiC42(13)s.t.xik=1,pointiisservicedbycarK,0,other,(14)xijk=1,vehicleKtravelsfrompointitopointj,0,other,(15)∑k=1Kxik=1,∀i∈N,(16)∑i=1Nxikqi≤Q,∀k∈K,(17)∑j=1Nx0jk=∑i=1Nxi0k≤1,∀k∈K.Equation (12) is the objective optimization function which is aimed at finding the minimum total cost during the whole distribution process. Equation (13)–Equation (17) are the constraint conditions, where Equation (13) is a 0-1 variable; the value 1 means the given customer is serviced by the given vehicle. Equation (14) is a 0-1 variable too; the value 1 means the given route is travelled by the given vehicle. Equation (15) means each customer point is serviced by one and only one vehicle for distribution. Equation (16) means that the total demand served by the vehicle cannot exceed its maximum load. Equation (17) means that the vehicle can only start or end once.


## 3.1. Distance and Transit Time Analysis

In traditional methods, the distance between two points is achieved by Euclidean distance formula and the transit time is achieved by the distance divided by assumed average speed. But, the actual transit time is influenced greatly by road congestion, traffic accidents, traffic control, raining, snowing, etc. Collecting real-time traffic data will result in large amounts of data, and Big Data analysis technology should be considered.In our literature, the navigation function of E-map is used to avoid the mentioned shortcomings. The distance matrix and the transit time matrix under actual traffic conditions among customer points (including the distribution center) are derived from the E-map public platform [21].As Figure1 shows, there are two points with their latitudes and longitudes on the E-map.Figure 1
Two points on the E-map.As Figure2 shows, the distance and the transit time between the two points are shown by the navigation function of the E-map platform.Figure 2
The output of navigation result.In our study,tij represents transit time between customer i and customer j, and dij is the distance between customer i and customer j. As abovementioned, tij and dij are both achieved from the navigation result based on E-map API.


## 3.2. Loss Cost Target Construction

### 3.2.1. Soft Time Window Penalty Cost

Each customer has different numbers of good needs and a soft time window. Penalty cost will be incurred if the vehicle arrives beyond the time window boundary. The service time of each customer point isuii=1,2,⋯,N, Eti,Lti is the range of the time window of customer point i, C11 is the penalty coefficient when vehicle ki arrives at customer i earlier than Eti, and C12 is the penalty coefficient when vehicle ki leaves from customer i later than Lti. The arrival time of vehicle k at customer i is Tki. According to the analysis above, the penalty cost of cold chain can be defined as the following:
(1)C1=C11∑k=1K∑i=1NxikmaxEti−Tki,0+C12∑k=1K∑i=1NxikmaxTki+ui−Lti,0,s.t.Tki≥Tk0+t0i,Tki+1≥Tki+tii+1+ui.


### 3.2.2. Fuel Cost

The fuel consumption of the distribution vehicle is inevitable for completing the distribution task. We assume that the fuel cost is proportional to the load of the vehicle pert unit distance. LetQk be the load of the k vehicle and Pk be the load rate of the kth vehicle.
(2)Qk=∑i=1nxikqi,Pk=QkQ.SetsC21 and C22 represent no load and full load fuel consumption costs per unit distance of the distribution vehicle. C2k means the fuel consumption cost per unit distance of the kth vehicle.
(3)C2k=C21+pkC22−C21.The total cost of fuel consumptionC2 can be expressed as
(4)C2=∑k=1K∑i=1N∑j=1NxijkdijC21+PkC22−C21.


### 3.2.3. Refrigeration Cost

The temperature of the cold-chain logistic must be maintained at a certain low level to keep the freshness of goods. The refrigeration function of vehicles must be operated to achieve the required temperature, which causes the refrigeration cost immediately.(1)
Refrigeration cost in transitThe refrigeration cost generated by maintaining the required low temperature per unit time in transit isC31; the total cooling cost C31 in transit can be expressed as
(5)C31=∑k=1K∑i=1N∑j=1NxijktijC31.(2)
Refrigeration cost during unloadingOpening the door of the compartment will cause the cool air inside to flow out and the hot air from outside to flow in. To maintain the temperature inside, more energy will be consumed.C32 is set to be the cooling cost per unit time during opening the door, and the total refrigeration cost of unloading C32 is
(6)C32=∑k=1K∑i=1NxikuiC32.C3 represents the total cooling cost, which consists of the total cooling cost in transit, and the total cooling cost of unloading can be expressed as
(7)C3=∑k=1K∑i=1N∑j=1NxijktijC31+xikuiC32.


### 3.2.4. Cargo Damage Cost

The most cargo of the cold chain is fresh goods; with the increase of transit time even in low temperature, the growth of microorganisms will happen. When opening the door of compartment during unloading process, the temperature inside will be unstable which results in the damage of fresh goods more significantly. As mentioned above,p is the unit price of goods, and Qki is the load of the kth vehicle when it leaves from point i.
(1)
Cost of cargo damage in transitC41 is given as the loss coefficient of goods per unit weight per unit time in transit; the total cost of cargo damage in transit is as follows:
(8)C41=∑k=1K∑i=1N∑j=1NxijktijpQkiC41.(2)
Cost of cargo damage during unloadingDuring unloading at a customer point, the door of the compartment is open and the cost of cargo damage will increase obviously. Let the cargo damage coefficient of goods per unit weight per unit time when unloading beC42; the cost of cargo damage during unloading can be expressed as
(9)C42=∑k=1K∑i=1NxikuipQkiC42.The total cost of damage is(10)C4=pQki∑k=1K∑i=1N∑j=1NxijktijC41+xikuiC42.In conclusion, the mathematical model of minimum loss target can be expressed as(11)minZ=C1+C2+C3+C4.


## 3.2.1. Soft Time Window Penalty Cost

Each customer has different numbers of good needs and a soft time window. Penalty cost will be incurred if the vehicle arrives beyond the time window boundary. The service time of each customer point isuii=1,2,⋯,N, Eti,Lti is the range of the time window of customer point i, C11 is the penalty coefficient when vehicle ki arrives at customer i earlier than Eti, and C12 is the penalty coefficient when vehicle ki leaves from customer i later than Lti. The arrival time of vehicle k at customer i is Tki. According to the analysis above, the penalty cost of cold chain can be defined as the following:
(1)C1=C11∑k=1K∑i=1NxikmaxEti−Tki,0+C12∑k=1K∑i=1NxikmaxTki+ui−Lti,0,s.t.Tki≥Tk0+t0i,Tki+1≥Tki+tii+1+ui.


## 3.2.2. Fuel Cost

The fuel consumption of the distribution vehicle is inevitable for completing the distribution task. We assume that the fuel cost is proportional to the load of the vehicle pert unit distance. LetQk be the load of the k vehicle and Pk be the load rate of the kth vehicle.
(2)Qk=∑i=1nxikqi,Pk=QkQ.SetsC21 and C22 represent no load and full load fuel consumption costs per unit distance of the distribution vehicle. C2k means the fuel consumption cost per unit distance of the kth vehicle.
(3)C2k=C21+pkC22−C21.The total cost of fuel consumptionC2 can be expressed as
(4)C2=∑k=1K∑i=1N∑j=1NxijkdijC21+PkC22−C21.


## 3.2.3. Refrigeration Cost

The temperature of the cold-chain logistic must be maintained at a certain low level to keep the freshness of goods. The refrigeration function of vehicles must be operated to achieve the required temperature, which causes the refrigeration cost immediately.(1)
Refrigeration cost in transitThe refrigeration cost generated by maintaining the required low temperature per unit time in transit isC31; the total cooling cost C31 in transit can be expressed as
(5)C31=∑k=1K∑i=1N∑j=1NxijktijC31.(2)
Refrigeration cost during unloadingOpening the door of the compartment will cause the cool air inside to flow out and the hot air from outside to flow in. To maintain the temperature inside, more energy will be consumed.C32 is set to be the cooling cost per unit time during opening the door, and the total refrigeration cost of unloading C32 is
(6)C32=∑k=1K∑i=1NxikuiC32.C3 represents the total cooling cost, which consists of the total cooling cost in transit, and the total cooling cost of unloading can be expressed as
(7)C3=∑k=1K∑i=1N∑j=1NxijktijC31+xikuiC32.


## 3.2.4. Cargo Damage Cost

The most cargo of the cold chain is fresh goods; with the increase of transit time even in low temperature, the growth of microorganisms will happen. When opening the door of compartment during unloading process, the temperature inside will be unstable which results in the damage of fresh goods more significantly. As mentioned above,p is the unit price of goods, and Qki is the load of the kth vehicle when it leaves from point i.
(1)
Cost of cargo damage in transitC41 is given as the loss coefficient of goods per unit weight per unit time in transit; the total cost of cargo damage in transit is as follows:
(8)C41=∑k=1K∑i=1N∑j=1NxijktijpQkiC41.(2)
Cost of cargo damage during unloadingDuring unloading at a customer point, the door of the compartment is open and the cost of cargo damage will increase obviously. Let the cargo damage coefficient of goods per unit weight per unit time when unloading beC42; the cost of cargo damage during unloading can be expressed as
(9)C42=∑k=1K∑i=1NxikuipQkiC42.The total cost of damage is(10)C4=pQki∑k=1K∑i=1N∑j=1NxijktijC41+xikuiC42.In conclusion, the mathematical model of minimum loss target can be expressed as(11)minZ=C1+C2+C3+C4.


## 3.3. Mathematical Model of Minimum Loss Analysis

According to the above analysis, the minimum loss model of the cold-chain distribution problem is as follows:(12)minZ=C11∑k=1K∑i=1NxikmaxEti−Tki,0+C12∑k=1K∑i=1NxikmaxTki+ui−Lti,0+∑k=1K∑i=1N∑j=1NxijkC21+PkC22−C21+∑k=1K∑i=1N∑j=1NxijktijC31+xikuiC32pQki∑k=1K∑i=1N∑j=1NxijktijC41+xikuiC42(13)s.t.xik=1,pointiisservicedbycarK,0,other,(14)xijk=1,vehicleKtravelsfrompointitopointj,0,other,(15)∑k=1Kxik=1,∀i∈N,(16)∑i=1Nxikqi≤Q,∀k∈K,(17)∑j=1Nx0jk=∑i=1Nxi0k≤1,∀k∈K.Equation (12) is the objective optimization function which is aimed at finding the minimum total cost during the whole distribution process. Equation (13)–Equation (17) are the constraint conditions, where Equation (13) is a 0-1 variable; the value 1 means the given customer is serviced by the given vehicle. Equation (14) is a 0-1 variable too; the value 1 means the given route is travelled by the given vehicle. Equation (15) means each customer point is serviced by one and only one vehicle for distribution. Equation (16) means that the total demand served by the vehicle cannot exceed its maximum load. Equation (17) means that the vehicle can only start or end once.


## 4. Ant Colony Optimization Algorithm

### 4.1. Ant Colony Optimization Algorithm

The ant colony optimization (ACO) algorithm with parallelism, positive feedback, and strong robustness is used to solve the NP-hard and highly constrained problem. In the 1990s, Italian scholars Dorigo and Maniezzo found that ants will release a certain amount of substance called pheromone in the path when they are searching for food. When they need to choose a path, they prefer to choose the path with high concentration of pheromone. At last, the route with the highest concentration will be selected as an optimist route between their home and the food source.The ACO algorithm is defined as follows:(18)pijkτijαt⋅ηijβt∑s∈allowedkτisαt⋅ηisβt,ifj∈allowedk,0,otherwise,τijt+n=ρ⋅τijt+Δτij,Δτij=∑k=1mΔτij,where pijk is probability of ant k transferring from point i to point j at time t; ηijt is the heuristic function, and β is the heuristic function important factor; τijt means the pheromone concentration on the path at time t, and α is the pheromone important factor. allowedk represents points that ant k is allowed to select in the next step; τijt+n is the pheromone update over time; ρ is the pheromone factor; and Δτijk represents the pheromone enhancer.The rough process of the algorithm is as follows:Algorithm 1
Initialize the coefficient variables and pheromone;for iteration 1,…,M{Initialize a taboo table and a start position for each ant.for ant, 1,…,N{Select the next node according to rulesThe selected node is stored in the taboo tableIf all nodes are stored in taboo tables, the iteration is completed, break;The total cost is calculatedThe local pheromone is updated}The global pheromone is updatedAll ants’ total costs are comparedThe current optimal solution is stored}The process is stopped and the current path with the shortest cost is output.


### 4.2. Combining ACO Algorithm with PLS

We applied Pareto local search (PLS) in our paper. After finding a path, PLS is used to further optimize the path. PLS proposed by Paquete et al. is a heuristic algorithm for tackling NP-hard multiobjective combinatorial optimization problems in the Pareto sense [22]. Pareto dominance defines a partial order on the set of feasible solutions. The goal when tackling the multiobjective problems in the Pareto sense is to find the set of Pareto-optimal solutions. The weak component-wise ordering is used as a mutually nondominated criterion in PLS. Let a solution s∈A and its neighborhood s′∈Ns. If each s′ in Ns is not dominated by any solution in A, s is marked visited and is added into A. When A contains only solutions that have been visited, a Pareto local optimum is achieved [23].An improved ACO with PLS is proposed in our research to deal with the mode. In our approach, ants are used to optimize solutions and generate new Pareto solutions, and the new pheromone updating strategy is used to control the path selection. The procedure is described as follows:Step 1.
Computepijk of n customers and update the trail level τij_new for each k.Step 2.
Run PLS, and update the current solution according to calculating results.Step 3.
Compare (Cost, CostBest), and record the costBest and BestSolution.Step 4.
For each move (i,j) in BestSolution to update the trail level τij_new.Step 5.
Repeat step 1 until the maximum computation time reaches.The proportion of the pheromone update model of the ACO algorithm with PLS is as follows:(19)ω=1−k1Smax−StSmax,(20)St=−K∑j∈Dpijtlogpijt,where St means information entropy value and pijt is transition probability. Equation (19) means that the proportion of the pheromone update will decreases with the increase of the pheromone value. k1 is the positive constant, and its value range is [0,1]. The value 1 means the proportion of the pheromone update is high influenced by the information entropy value, and 0.75 is usually used as the value of k1. This definition combined a sequence of arithmetic with information entropy to adjust the pheromone adaptively.


## 4.1. Ant Colony Optimization Algorithm

The ant colony optimization (ACO) algorithm with parallelism, positive feedback, and strong robustness is used to solve the NP-hard and highly constrained problem. In the 1990s, Italian scholars Dorigo and Maniezzo found that ants will release a certain amount of substance called pheromone in the path when they are searching for food. When they need to choose a path, they prefer to choose the path with high concentration of pheromone. At last, the route with the highest concentration will be selected as an optimist route between their home and the food source.The ACO algorithm is defined as follows:(18)pijkτijαt⋅ηijβt∑s∈allowedkτisαt⋅ηisβt,ifj∈allowedk,0,otherwise,τijt+n=ρ⋅τijt+Δτij,Δτij=∑k=1mΔτij,where pijk is probability of ant k transferring from point i to point j at time t; ηijt is the heuristic function, and β is the heuristic function important factor; τijt means the pheromone concentration on the path at time t, and α is the pheromone important factor. allowedk represents points that ant k is allowed to select in the next step; τijt+n is the pheromone update over time; ρ is the pheromone factor; and Δτijk represents the pheromone enhancer.The rough process of the algorithm is as follows:Algorithm 1
Initialize the coefficient variables and pheromone;for iteration 1,…,M{Initialize a taboo table and a start position for each ant.for ant, 1,…,N{Select the next node according to rulesThe selected node is stored in the taboo tableIf all nodes are stored in taboo tables, the iteration is completed, break;The total cost is calculatedThe local pheromone is updated}The global pheromone is updatedAll ants’ total costs are comparedThe current optimal solution is stored}The process is stopped and the current path with the shortest cost is output.


## 4.2. Combining ACO Algorithm with PLS

We applied Pareto local search (PLS) in our paper. After finding a path, PLS is used to further optimize the path. PLS proposed by Paquete et al. is a heuristic algorithm for tackling NP-hard multiobjective combinatorial optimization problems in the Pareto sense [22]. Pareto dominance defines a partial order on the set of feasible solutions. The goal when tackling the multiobjective problems in the Pareto sense is to find the set of Pareto-optimal solutions. The weak component-wise ordering is used as a mutually nondominated criterion in PLS. Let a solution s∈A and its neighborhood s′∈Ns. If each s′ in Ns is not dominated by any solution in A, s is marked visited and is added into A. When A contains only solutions that have been visited, a Pareto local optimum is achieved [23].An improved ACO with PLS is proposed in our research to deal with the mode. In our approach, ants are used to optimize solutions and generate new Pareto solutions, and the new pheromone updating strategy is used to control the path selection. The procedure is described as follows:Step 1.
Computepijk of n customers and update the trail level τij_new for each k.Step 2.
Run PLS, and update the current solution according to calculating results.Step 3.
Compare (Cost, CostBest), and record the costBest and BestSolution.Step 4.
For each move (i,j) in BestSolution to update the trail level τij_new.Step 5.
Repeat step 1 until the maximum computation time reaches.The proportion of the pheromone update model of the ACO algorithm with PLS is as follows:(19)ω=1−k1Smax−StSmax,(20)St=−K∑j∈Dpijtlogpijt,where St means information entropy value and pijt is transition probability. Equation (19) means that the proportion of the pheromone update will decreases with the increase of the pheromone value. k1 is the positive constant, and its value range is [0,1]. The value 1 means the proportion of the pheromone update is high influenced by the information entropy value, and 0.75 is usually used as the value of k1. This definition combined a sequence of arithmetic with information entropy to adjust the pheromone adaptively.


## 5. Experimental Results and Discussion

### 5.1. Data Set Description

In our experiments, the distribution situation of JingKeLong in Beijing City with one distribution center and 20 customer points is selected to verify our proposed mathematical model and algorithm. The distribution center has 4 vehicles of the same type with an average speed of 60 km/h, a maximum driving time of 1200 minutes, and a maximum load of 2000 kg. Table1 shows the settings of parameters in our given minimum loss model and ACO+PLS algorithm.Table 1
Parameter setting.
Parameter of the symbolMeaningValueParameter of the symbolMeaningValueKThe total vehicles4pUnit price of goods5NThe total customer points20mNumber of ants50QFull load of one vehicle2000nNumber of iterations60C11Penalty coefficients earlier than TW2.0AlphaThe pheromone important factor1C12Penalty coefficients later than TW3.0BetaHeuristic function important factor3C21Unit fuel cost of no load2.5RhoThe pheromone factor0.85C22Unit fuel cost of full load3.0k1Positive constant0.75C31Unit refrigeration cost in transit2.0C32Unit refrigeration cost in unloading2.5C41Coefficient of damage in transit0.01C42Coefficient of damage in unloading0.015The data of the distribution center and customers are shown in Table2. The distribution center is defined as 0, and customer points are defined as 1 to 20. Figure 3 shows the latitude and longitude of each point in the XY system, and Figure 4 shows the actual locations of them in E-map.Table 2
The data of the distribution center and customers.
Number of customersLatitude and longitude on E-mapThe demand of customer (kg)Starting time (min)Ending time (min)Server time (min)Latitude (easting)Longitude (northing)0116.484339.876800001116.48608139.801535300912967152116.4328739.8516571100825870303116.57147439.85701112565146104116.5177439.86957100727782105116.25816539.8962872001567106116.59800139.920533150621702107116.47589239.92811150170225108116.65939339.9282450255324209116.11054839.9434143005346052010116.5012139.9677271003574101011116.45043339.9711269504485053012116.40820639.9737341256527211013116.33998839.97835415030901014116.46883240.0061831505676201015116.58448540.0079855503844292016116.60142140.0546171504755281017116.43929740.057803100991481018116.67036340.1406511501792541019116.22662640.2281014002783452020116.65328440.332406300107320Figure 3
Distribution of warehouse and customers in theXY system.Figure 4
Geographic locations of warehouse and customers in E-map.Our proposed method needs an actual transit time matrix and distance matrix both derived from E-map API. We supposed that the coordinates of points are stored into vertex [N] [2]; the core codes with E-map API to achieve the transit time matrix and distance matrix are described as follows.Algorithm 2
......var output_t = “The transit time:”;var output_d = “The traveling distance:”;var distance_time_search = function (customerPoints){.if (vehicle_transit.getStatus() != BMAP_STATUS_SUCCESS){.return ;}var t_d_result = customerPoints.getPlan(0);travel_time += t_d_result.getDuration(true) + “\n”; //get transit time of two points.trave_distance+= t_d_result.getDistance(true)+”\n”; //get travel distance of two points.}var vehicle_transit = new BMapGL.DrivingRoute(map, {renderOptions: {map: map},onSearchComplete: distance_time_search;}});for(var i=0;i<=NumberOfCustomers;i++)for( var j=i+1;j<= NumberOfCustomers;j++){var start=new BMapGL.Point(vertex[i][0], vertex[i][1]);var end=new BMapGL.Point(vertex[j][0], vertex[j][1]);vehicle_transit.search(start, end);}......Table3 shows the transit time matrix among customer points and distribution center obtained by E-map API.Table 3
Transit time matrix under actual traffic conditions (minutes).
0123456789101112131415161718192000222525263513162443251318232521241543464012201719282721263817283451343732434059585522517025202632301843332835343837456342616132519250162424142636342843573427515553413242628201609232421222838523529393629575153535272624902628352750423351365065495449716132132242326013272026253331384119595335507162630142428130923171424252027284847424982438182621352790242732423524633237415857943174336222720232403940464651606342486978102528333428502617273901829232534263351484911133428283842251432401808161023402026343812185135435233332442462980161322342643244013233434573551312535462316160223436363252501425373834293638202451251013220182435392637152132372739504127636034232234180193249364916244345513665192832632640343624190253145461715406355294959483742332026363532250433841184359425357545347414851264332394931430334519465861415149354258694834245226364538330422040556132537150495778493840503749464145420Table4 shows the distance matrix among customer points and distribution center obtained by E-map API.Table 4
The distance matrix derived from the E-map navigation function (KM).
01234567891011121314151617181920002213.620.717.66.519.25.611.321.231.616.811.47.217.441.62214.253.441.635.112209.112.711.32119.526.613.129.545.423.527.832.737.927.339.535.462.286.464.7213.69.1012.319.111.722.422.96.529.316.24036.3272715.437.255.667.33469.1320.712.712.309.88.815.218.210.324.639.627.421.524.629.831.555.634.152.574.259.1417.611.319.19.803.816.910.517.91928.918.733.626.149.621.925.728.353.179.855.356.52111.78.83.8022.624.126.832.217.837.533.91929.84557.83374.844.660.5619.219.522.415.216.922.60617.111.521.932.324.420.612.317.664.455.630.431.364.175.626.622.918.210.524.16018.38.211.916.9613.52121.457.240.918.54438811.313.16.510.317.926.817.118.3021.426.228.714.936.626.518.725.459.834.452.465.2921.229.529.324.61932.211.58.221.4029.432.655.442.136.940.285.544.641.986.9751031.645.416.239.628.917.821.911.926.229.4017.712.527.826.91724.681.227.25259.41116.823.54027.418.737.532.316.928.732.617.705.86.111.624.213.226.736.851.1381211.427.836.321.533.633.924.4614.955.412.55.8011.722.19.37.726.638.435.752.7137.232.72724.626.11920.613.536.642.127.86.111.701428.437.15944.72132.91417.437.92729.849.629.812.32126.536.926.911.622.11408.721.216.637.64833.71541.627.315.431.521.94517.621.418.740.21724.29.328.48.708.516.423.357.350.5162239.537.255.625.757.864.457.225.485.524.613.27.737.121.28.5018.851.917.650.31714.235.455.634.128.33355.640.959.844.681.226.726.65916.616.418.803644.129.71853.462.267.352.553.174.830.418.534.441.927.236.838.444.737.623.351.93602547.51941.686.43474.279.844.631.34452.486.95251.135.7214857.317.644.125056.32035.164.769.159.155.360.564.13865.27559.43852.732.933.750.550.329.747.556.30


### 5.2. Experimental Result and Discussion

In the traditional method, the distance matrix is achieved by the Euclidean distance formula. The Euclidean distance formula of two points on earth is defined as the following:(21)DAB=R∗arccoscoslonA∗coslonB∗coslatA−latB+sinlonA∗sinlonB∗π180,where A and B are two points on the surface of earth and R is the radius of the earth. LatA and lonA represent the longitude and latitude of point A.First, the traditional method ACO with the Euclidean distance matrix is performed to solve the VRP with minimum cost. Figures5 and 6 are the optimal distribution routes and the trend of cost during iteration processed. The optimal distribution routing scheme, the total transit time, and the total cost based on the traditional method are shown in Table 5.Figure 5
The optimal path plan based on the traditional method.Figure 6
The trend of cost during the iteration.Table 5
Distribution paths of 4 vehicles based on the traditional method.
VehiclesThe route of each vehicleThe realistic delivery on E-map (minutes)Fist vehicle0−>20−>3−>7−>18−>8−>19−>10−>4−>0334Second vehicle0−>6−>12−>2−>1−>0120Third vehicle0−>15−>16−>11−>9−>0163Fourth vehicle0−>5−>13−>17−>14−>0182The total time799The total cost4504.5Second, our proposed method ACO+PLS with the actual transit time matrix and travel distance matrix is performed to solve the VRP with minimum cost. Figures7 and 8 are the calculated distribution routes on map and the trend of cost during iteration processed. The optimal distribution routing scheme, the total transit time, and the total cost based on the proposed method are shown in Table 6.Figure 7
The optimal path plan based on the proposed method.Figure 8
The trend of cost during the iteration.Table 6
Distribution paths of 3 vehicles based on the proposed method.
VehiclesThe route of each vehicleThe realistic delivery on E-map (minutes)Fist vehicle0−>20−>17−>18−>19−>10−>15−>16−>14−>4−>0331Second vehicle0−>5−>3−>7−>8−>11−>12−>0113Third vehicle0−>13−>9−>6−>2−>1−>0160The total time593The total cost3673From the above compare experimental results, we can see that ACO+PLS with actual traffic conditions can get smaller total cost and save more transit time. Our proposed method can be used to get the optimal solution effectively for VRP of cold-chain logistics.


## 5.1. Data Set Description

In our experiments, the distribution situation of JingKeLong in Beijing City with one distribution center and 20 customer points is selected to verify our proposed mathematical model and algorithm. The distribution center has 4 vehicles of the same type with an average speed of 60 km/h, a maximum driving time of 1200 minutes, and a maximum load of 2000 kg. Table1 shows the settings of parameters in our given minimum loss model and ACO+PLS algorithm.Table 1
Parameter setting.
Parameter of the symbolMeaningValueParameter of the symbolMeaningValueKThe total vehicles4pUnit price of goods5NThe total customer points20mNumber of ants50QFull load of one vehicle2000nNumber of iterations60C11Penalty coefficients earlier than TW2.0AlphaThe pheromone important factor1C12Penalty coefficients later than TW3.0BetaHeuristic function important factor3C21Unit fuel cost of no load2.5RhoThe pheromone factor0.85C22Unit fuel cost of full load3.0k1Positive constant0.75C31Unit refrigeration cost in transit2.0C32Unit refrigeration cost in unloading2.5C41Coefficient of damage in transit0.01C42Coefficient of damage in unloading0.015The data of the distribution center and customers are shown in Table2. The distribution center is defined as 0, and customer points are defined as 1 to 20. Figure 3 shows the latitude and longitude of each point in the XY system, and Figure 4 shows the actual locations of them in E-map.Table 2
The data of the distribution center and customers.
Number of customersLatitude and longitude on E-mapThe demand of customer (kg)Starting time (min)Ending time (min)Server time (min)Latitude (easting)Longitude (northing)0116.484339.876800001116.48608139.801535300912967152116.4328739.8516571100825870303116.57147439.85701112565146104116.5177439.86957100727782105116.25816539.8962872001567106116.59800139.920533150621702107116.47589239.92811150170225108116.65939339.9282450255324209116.11054839.9434143005346052010116.5012139.9677271003574101011116.45043339.9711269504485053012116.40820639.9737341256527211013116.33998839.97835415030901014116.46883240.0061831505676201015116.58448540.0079855503844292016116.60142140.0546171504755281017116.43929740.057803100991481018116.67036340.1406511501792541019116.22662640.2281014002783452020116.65328440.332406300107320Figure 3
Distribution of warehouse and customers in theXY system.Figure 4
Geographic locations of warehouse and customers in E-map.Our proposed method needs an actual transit time matrix and distance matrix both derived from E-map API. We supposed that the coordinates of points are stored into vertex [N] [2]; the core codes with E-map API to achieve the transit time matrix and distance matrix are described as follows.Algorithm 2
......var output_t = “The transit time:”;var output_d = “The traveling distance:”;var distance_time_search = function (customerPoints){.if (vehicle_transit.getStatus() != BMAP_STATUS_SUCCESS){.return ;}var t_d_result = customerPoints.getPlan(0);travel_time += t_d_result.getDuration(true) + “\n”; //get transit time of two points.trave_distance+= t_d_result.getDistance(true)+”\n”; //get travel distance of two points.}var vehicle_transit = new BMapGL.DrivingRoute(map, {renderOptions: {map: map},onSearchComplete: distance_time_search;}});for(var i=0;i<=NumberOfCustomers;i++)for( var j=i+1;j<= NumberOfCustomers;j++){var start=new BMapGL.Point(vertex[i][0], vertex[i][1]);var end=new BMapGL.Point(vertex[j][0], vertex[j][1]);vehicle_transit.search(start, end);}......Table3 shows the transit time matrix among customer points and distribution center obtained by E-map API.Table 3
Transit time matrix under actual traffic conditions (minutes).
0123456789101112131415161718192000222525263513162443251318232521241543464012201719282721263817283451343732434059585522517025202632301843332835343837456342616132519250162424142636342843573427515553413242628201609232421222838523529393629575153535272624902628352750423351365065495449716132132242326013272026253331384119595335507162630142428130923171424252027284847424982438182621352790242732423524633237415857943174336222720232403940464651606342486978102528333428502617273901829232534263351484911133428283842251432401808161023402026343812185135435233332442462980161322342643244013233434573551312535462316160223436363252501425373834293638202451251013220182435392637152132372739504127636034232234180193249364916244345513665192832632640343624190253145461715406355294959483742332026363532250433841184359425357545347414851264332394931430334519465861415149354258694834245226364538330422040556132537150495778493840503749464145420Table4 shows the distance matrix among customer points and distribution center obtained by E-map API.Table 4
The distance matrix derived from the E-map navigation function (KM).
01234567891011121314151617181920002213.620.717.66.519.25.611.321.231.616.811.47.217.441.62214.253.441.635.112209.112.711.32119.526.613.129.545.423.527.832.737.927.339.535.462.286.464.7213.69.1012.319.111.722.422.96.529.316.24036.3272715.437.255.667.33469.1320.712.712.309.88.815.218.210.324.639.627.421.524.629.831.555.634.152.574.259.1417.611.319.19.803.816.910.517.91928.918.733.626.149.621.925.728.353.179.855.356.52111.78.83.8022.624.126.832.217.837.533.91929.84557.83374.844.660.5619.219.522.415.216.922.60617.111.521.932.324.420.612.317.664.455.630.431.364.175.626.622.918.210.524.16018.38.211.916.9613.52121.457.240.918.54438811.313.16.510.317.926.817.118.3021.426.228.714.936.626.518.725.459.834.452.465.2921.229.529.324.61932.211.58.221.4029.432.655.442.136.940.285.544.641.986.9751031.645.416.239.628.917.821.911.926.229.4017.712.527.826.91724.681.227.25259.41116.823.54027.418.737.532.316.928.732.617.705.86.111.624.213.226.736.851.1381211.427.836.321.533.633.924.4614.955.412.55.8011.722.19.37.726.638.435.752.7137.232.72724.626.11920.613.536.642.127.86.111.701428.437.15944.72132.91417.437.92729.849.629.812.32126.536.926.911.622.11408.721.216.637.64833.71541.627.315.431.521.94517.621.418.740.21724.29.328.48.708.516.423.357.350.5162239.537.255.625.757.864.457.225.485.524.613.27.737.121.28.5018.851.917.650.31714.235.455.634.128.33355.640.959.844.681.226.726.65916.616.418.803644.129.71853.462.267.352.553.174.830.418.534.441.927.236.838.444.737.623.351.93602547.51941.686.43474.279.844.631.34452.486.95251.135.7214857.317.644.125056.32035.164.769.159.155.360.564.13865.27559.43852.732.933.750.550.329.747.556.30


## 5.2. Experimental Result and Discussion

In the traditional method, the distance matrix is achieved by the Euclidean distance formula. The Euclidean distance formula of two points on earth is defined as the following:(21)DAB=R∗arccoscoslonA∗coslonB∗coslatA−latB+sinlonA∗sinlonB∗π180,where A and B are two points on the surface of earth and R is the radius of the earth. LatA and lonA represent the longitude and latitude of point A.First, the traditional method ACO with the Euclidean distance matrix is performed to solve the VRP with minimum cost. Figures5 and 6 are the optimal distribution routes and the trend of cost during iteration processed. The optimal distribution routing scheme, the total transit time, and the total cost based on the traditional method are shown in Table 5.Figure 5
The optimal path plan based on the traditional method.Figure 6
The trend of cost during the iteration.Table 5
Distribution paths of 4 vehicles based on the traditional method.
VehiclesThe route of each vehicleThe realistic delivery on E-map (minutes)Fist vehicle0−>20−>3−>7−>18−>8−>19−>10−>4−>0334Second vehicle0−>6−>12−>2−>1−>0120Third vehicle0−>15−>16−>11−>9−>0163Fourth vehicle0−>5−>13−>17−>14−>0182The total time799The total cost4504.5Second, our proposed method ACO+PLS with the actual transit time matrix and travel distance matrix is performed to solve the VRP with minimum cost. Figures7 and 8 are the calculated distribution routes on map and the trend of cost during iteration processed. The optimal distribution routing scheme, the total transit time, and the total cost based on the proposed method are shown in Table 6.Figure 7
The optimal path plan based on the proposed method.Figure 8
The trend of cost during the iteration.Table 6
Distribution paths of 3 vehicles based on the proposed method.
VehiclesThe route of each vehicleThe realistic delivery on E-map (minutes)Fist vehicle0−>20−>17−>18−>19−>10−>15−>16−>14−>4−>0331Second vehicle0−>5−>3−>7−>8−>11−>12−>0113Third vehicle0−>13−>9−>6−>2−>1−>0160The total time593The total cost3673From the above compare experimental results, we can see that ACO+PLS with actual traffic conditions can get smaller total cost and save more transit time. Our proposed method can be used to get the optimal solution effectively for VRP of cold-chain logistics.


## 6. Conclusion

This paper studies the model of VPR with minimum cost of cold-chain distribution. Minimizing the total cost of distribution can maximize the economic benefit of distribution enterprises. To consider the dynamic change of transit speed under actual traffic conditions, the transit time matrix and distance matrix are both derived from the navigation function based on E-map API. We proposed a heuristic approach ACO+PLS to solve the minimum loss model of cold-chain logistics. The experiments show that our proposed method has strong applicability and potential advantages in cold-chain distribution. In future studies, a combination of multiple spatial information technologies such as geographic information technology and remote sensing technology can be realized to make the problem more practical.



---

*Source: 2930366-2021-12-10.xml*