# Optimization Model and Algorithm of Cigarette Distribution Route Based on Cluster Analysis

Authors: Hong Xiao; Huizhen Xu; Xiaokun Wang; Bo Hang
Journal: Wireless Communications and Mobile Computing
 (2022)

---

## Abstract

Cigarette distribution is an important part of the tobacco logistics supply chain, and it will affect the distribution cost, efficiency, and service quality. In this paper，we choose a two-stage optimization method to analyze the cigarette distribution route across the Administrative Region in China. First, we use a K-means clustering algorithm to generate an initial center. Then, we optimize the clustering region with the vehicle load and workload as constraints. Finally, we establish the cigarette distribution route optimization model by taking the lowest distribution cost as the objective function. In the model calculation, we not only select the adaptive genetic algorithm as the solution algorithm but also use the scanning algorithm to generate the initial population ,design adaptive crossover and mutation probability, and also design inverse operator to correct the error, and we wish to improve the solution performance of the algorithm by this way.In addition, taking the cigarette distribution of Chongqing Tobacco Logistics as an example, this paper analyzes the number of vehicles, average loading rate, and total distribution distance and distribution cost, and finds that the cross-regional joint distribution route is better than independent distribution. And we also find that using clustering analysis and the genetic algorithm could reduce the amount of computation and solve the problem effectively.


---

## Body

## 1. Introduction

China is the country most seriously harmed by tobacco, and its cigarette consumption accounts for more than 43% of the world’s total. China is also a major cigarette producer. But there has been a slight decline in cigarette production in China in recent years. Production fell to 2,364.24 billion sticks in 2019 after falling to 2,382.57 billion sticks in 2016. In 2020, there was a slight increase, with cigarette production of 2,386.37 billion sticks, an increase of 0.9 percent year-on-year. With the improvement of residents’ health consciousness and the strengthening of tobacco control, the number of tobacco enterprises will continue to decrease in the future. Since 2015, the loss of China’s tobacco industry has continued to increase, reaching 2.28 billion yuan in 2019.The cost of cigarette is mainly composed of production cost and logistics cost. In China, the cigarette production technology and equipment of the tobacco industry have been upgraded step by step, more investment has been made in informatization, and the degree of informatization is relatively high. Many enterprises have implemented ERP, MES, and various central control systems, and the production efficiency and the cost control level have also been greatly improved with the enhancement of the degree of informatization. In 2020, the industry cost-expense margin of China’s tobacco industry was 29.81 percent, the second highest among 41 industries and well above the industry average of 7.11 percent. So, at present, the production cost of cigarettes in China can be crushed. In order to reduce the cost of cigarettes, it is necessary to study the distribution cost of cigarettes.In 2020, the total logistics expenses of cigarette enterprises nationwide totaled 11.122 billion yuan, an increase of 3.05 percent year-on-year, and the average logistics expenses per case were 93.49 yuan, a decrease of 0.65 percent year-on-year. Cigarette logistics costs accounted for 0.86% of sales revenue, down 0.03% year-on-year. The per capita delivery efficiency was 990.06 boxes per person, up 8.04% year-on-year. Cigarette distribution holds the characteristics of large-scale coverage, wide-range coverage, numerous and scattered merchants, timeliness, variety dispersion, and so on. Depending on the investigation, as cigarettes are controlled commodities in China, they are distributed according to administrative regions, and the distribution route is not reasonable, resulting in high distribution costs.Since joining the World Health Organization’s “Framework Treaty on Tobacco Control” in August 2005, the annual cigarette sales have slowed down in China, and the tobacco industry should face reducing operating costs and making more profits, so the tobacco industries have to reduce logistics costs to improve their economic benefits. Cigarette distribution is a key part of the tobacco supply chain. The reasonable choice of cigarette distribution route will affect cigarette distribution cost, work efficiency, and service quality.With the rising of oil price and labor cost, the cost of tobacco distribution has exceeded inventory costs. Route optimization of tobacco distribution vehicles is the most critical part of the tobacco distribution process, which is crucial to the cost, efficiency, and speed of the whole tobacco distribution. Therefore, how to choose the best route of cigarette delivery and send the cigarette to the merchant in time has become a heated issue in the research field of cigarette [1–3].In the field of combinatorial optimization problems, the vehicle routing problem (VRP) is one of the most challenges [4]. Many researchers put forward different methods of vehicle routing optimization and opinions in their respective research fields. For the problem of vehicle routing optimization of cigarette distribution, we get some inspiration from these researches. In the research of agricultural products’ logistics distribution path optimization, Liu et al. established a distribution route optimization model which is affected by the cost and has a time limitation to improve the efficiency of logistics distribution [5]. Schmid et al. investigated mode and user-type effects in the value of travel time savings (VTTS) using a pooled RP/SP Mixed Logit modeling approach for mode, route, and destination choice data [6]. Meanwhile, Liu et al. proposed a novel distributed robust chance constrained formulation and developed the sample average approximation method and a model-based hierarchical approach to handle the stochastic inventory routing problem of modules [7]. Then, Zhen et al. developed a hybrid particle swarm optimization algorithm and a hybrid genetic algorithm to solve the multijourney vehicle routing problem [8]. In order to overcome the limitation of distribution path, Gonzalez-R et al. proposed an iterated greedy heuristic based on the iterative [9]. Liu et al. developed a deep inverse reinforcement learning (IRL) algorithm and employed Dijkstra’s algorithm instead of value iteration, to determine the current policy and compute the gradient of IRL [10]. In order to optimize the route selection of long distance vehicles in Brazil, Mayerle et al. proposed a joint solution for the routing and scheduling problem with time window constraints, determining simultaneously a routing plan [11].Researchers put forward many optimization strategies of cigarette distribution path according to the distribution characteristics of tobacco logistics. Chen presented the ant colony algorithm for VRP with objective to minimize the delivery time and created two different system models through the case investigation of Hangzhou tobacco network [4]. One is the division of distribution route of distribution center, and the other is the optimization of single vehicle route. In order to optimize the cigarette distribution route, Hu et al. built a multicriteria balanced partition model and designed an immune coevolutionary algorithm with two stages to search the optimal balanced partitions then introduced three methods including fixed districts and routes, dynamic VRP scheduling, and periodic balanced partition for comparison to show the value of the method [12]. On uncertainty and fuzziness existing in the tobacco distribution field, Li et al. researched on the application of correction C-W saving algorithm applied in the field of tobacco distribution and showed that fixed C-W saving algorithm based on the centroid and deviation sorting is an effective algorithm, and it can be well applied to cigarette distribution network optimization problem [13]. Han et al. established the mathematical model of tobacco distribution with multivehicle and proposed a novel optimization distribution strategy, which combines Fuzzy C-Means Clustering (FCM) with Ant Colony Algorithm (ACA) and proved that this strategy can shorten the distance of tobacco distribution vehicle and optimize the distribution route [1].Based on the research of tobacco supply chain, some researchers have designed some management optimization models and management information systems, which reduce the cost of logistics distribution and improve the efficiency of distribution. Li and Luo analyzed the value chain embedded in the modern cigarette logistics and used value chain method to optimize cigarette logistics distribution [14]. Zheng et al. put forward the model of optimization of fuel management of distribution vehicle and optimized the tobacco company logistics distribution costs based on the information environment [15]. Then, Wei and Tu realized the whole process tracking and visual management of tobacco logistics through the establishment of RDC-SCM system, the establishment of cigarette logistics information system based on warehouse management [16].Most of the models are established and solved for classic vehicle routing problems, and the depth and detail of the research in further are needed. This paper takes the integration and optimization of cigarette distribution resources as the goal. Under the premise of ensuring the timeliness of cigarette distribution, it breaks the restrictions of administrative divisions and implements cross-regional joint distribution to optimize the cigarette distribution network and reduce logistics costs.


## 2. Two-Stage Optimization Model of Cigarette Distribution Route

### 2.1. Problem Definition and Model Framework Design

#### 2.1.1. Problem Definition

In order to establish the cigarette distribution model and calculate it, we simplify the vehicle routing problem in view of the tobacco industry control policy and the particularity of cigarette distribution.Those are as follows:(1)
There is a time window(2)
Nonfull load distribution problem: the cargo demand of each customer is less than the vehicle loading capacity, and it can carry goods of multiple customers when the same vehicle is stowed(3)
Single objective problem: the objective is to minimize the total cost of cigarette distribution(4)
Pure delivery problem: the distribution center provides the distribution service that can meet the needs of customers to multiple customers (specification and quantity of goods)(5)
Single distribution vehicle model: the distance of distribution vehicle is within the maximum theoretical driving distance


#### 2.1.2. Model Frame Design

Two-stage algorithm is a method which is widely used and has some achievements. The two-stage algorithm mainly includes path before grouping algorithm or group before path algorithm [20]. In this paper, we use the algorithm of grouping before routing, that is, the distribution area is generated first, and then, the optimal distribution route is obtained within this area.(1) First Stage: Distribution Areas’ Division. We determine the distribution center and distribution areas, transfer station (clustering center), and covered retail merchants according to the given distribution radiation radius. It is divided into different vehicle distribution areas according to the spatial geographical coordinates of retail merchants and different traffic conditions, such as transit area and direct distribution area. Clustering algorithm is used to divide or group customer nodes [17]. This makes the retail customers with relatively concentrated geographical distribution in the distribution area of the same car.(2) Second Stage: Distribution Route Optimization. The genetic algorithm is used to solve the TSP problem and the specific cigarette delivery routes or sequences of delivery vehicles within a single delivery area.Through the two-stage model algorithm to refine and decompose the vehicle routing problem of cigarette distribution, it can effectively improve the operational efficiency of the algorithm, shorten the calculation time, and get the approximate optimal solution.


### 2.2. Distribution Area Division Clustering and Optimization

Cigarette distribution route optimization needs clustering analysis of spatial geographic data of retail households, and it needs a high efficiency algorithm due to the large amount of data. K-means algorithm is suitable for data type data and insensitive to data input sequence. So it is more appropriate.


#### 2.2.1. Initial Clustering Center of K-Means

K-means clustering algorithm needs to solve two key problems, one is to determine the number of initial clustering centers, and the other is to determine the location of initial clustering centers.(1)
The number of initial cluster centers is the number of distribution areas, that is, the number of vehicles for cigarette distribution in these distribution areas. The number of distribution vehicles is set asm=cigarettedistributionvolume/singlevehicleload+1. Considering the maximum distance between nodes, and the fact that some vehicles cannot be fully loaded for distribution, so an additional vehicle is added for flexible configuration(2)
The location of the initial cluster center is usually randomly generated. In order to improve the clustering efficiency and reduce the amount of computation, we make some improvements to the generation of the initial clustering centers. Taking the geographical coordinates of each retail customer node as the center of the circle and the average distance between the geographical data of all retail merchants as the radius, the initial cluster center is determined according to the density of customer nodes in each circle [18]. Its implementation process mainly includes the following steps
(a)
Mark the geographical coordinates of each customer node in turn, and the mathematical distance between all nodes is calculated to generate the distance matrix. In order to simplify the model, it is assumed that the customer nodes are connected by a straight line, take the mathematical distance formula between two nodes as(1)Dij=xi−xj2+yi−yj2,i,j=1,2,⋯n(b)
Calculate the average distanceDij between all nodes in the database, let the initial radius of cluster center be R1=D, R2=2D(2)D=∑i=1n∑j=1nDijn×n(c)
Take each customer node as the center in turn, takeR1 as the initial radius to make a circle, and calculate the number of nodes in each circle, that is, sample density; for example, the sample density ρi of node Ai=xi,yi(d)
Sort from large to small according to sample density, the node with the highest sample density is taken as the first cluster centerE1, and the node with the second highest sample density is named E′. If the distance between E′and the first cluster center E1is E1E′≥R2, E′ is the cluster center E2. If the mathematical distance between the node and the previously determined cluster center is greater than R2; then, the node is the new cluster center. It iterates until m cluster centers are reached(e)
From the previous step, the numbers of cluster centers are mE1,E2,E3,⋯,andEm, as the initial clustering centers of K-meansStep 1.
set upρi=0Step 2.
Judge whether nodeA1=x1,y1 falls in the circle with Ai=xi,yi as the center and R1 as the radiu by A1Ai=xi−xj2+yi−yj2≤R1, if the above inequality is satisfied, then ρi=ρi+1Step 3.
Judge whether all nodes have been calculated, when the calculation is completed,ρiis the sample density of node Ai; otherwise, judge the next node A2, repeat the second step


#### 2.2.2. Clustering Optimization considering Vehicle Load

The initial cluster centers only consider the density distribution of customer nodes and do not consider the specific demand of customer nodes. The different cigarette distribution volume of customer nodes will affect the optimization result of cigarette distribution route. Therefore, the initial clustering results should be further optimized and adjusted according to the constraints of single vehicle loading capacity. The specific optimization and adjustment plan are as follows:Step 1.
Calculate the number of retail merchants ism, and each cluster and their total cigarette distribution is Lm.Step 2.
The total cigarette distribution volumeLm of retail merchants in each cluster is compared with the maximum vehicle volume Tm of single vehicle. If Lm≤Tm, it means that the clustering meets the constraints of vehicle load. If Lm>Tm, it means that the distribution volume of cigarette in the cluster exceeds the vehicle volume. The customer node farthest from the cluster center in the cluster is excluded, and it is proposed to be classified as the cluster nearest to the adjacent cluster center. Before the cluster is classified, it is necessary to judge the updated Lm≤Tm of the cluster. If Lm≤Tm is true, the node will be classified into the cluster. Otherwise, the cluster with poor Euclidean distance will be selected to judge Lm≤Tm. If the condition of constrained single vehicle loading capacity is met, the node will be classified into the cluster. Otherwise, the node will be searched and judged againStep 3.
When the total distribution volume of all cluster retailers is less than the single vehicle loading capacityLm≤Tm, the adjustment is completed.


#### 2.2.3. Clustering Optimization considering Workload

In order to avoid the local optimal situation, we should not only consider the constraints of vehicle capacity but also consider the cumulative workload within the distribution area. In the actual distribution work, the average daily working hours of distribution workers should be kept within the normal working hours of 8 hours.To calculate and evaluate the workload of cigarette distribution staff scientifically and reasonably has always been a difficult problem in cigarette distribution management. For one, the external and internal factors affecting the workload of distribution staff are numerous and complex. Extrinsic factors include distribution workload, distribution distance, road conditions, and the dispersion of retail nodes; intrinsic factors include staff’s work attitude and experience, the rationality of the division of labor, and the comprehensive quality of distribution staff, etc. For another, the complexity of the influencing factors and the lack of uniformity in the magnitude of the various influencing factors do not facilitate calculation and judgement.Therefore, some experts and scholars have proposed the concept of “panworkload,” which is used to calculate and evaluate the workload of each distribution route and to adjust the clustering results with the workload as the constraint condition. The concept of “panworkload” comprehensively considers the number of customer nodes, the volume of cigarettes distribution, the mileage of delivery, and specific traffic conditions. By assigning different weights to each influencing factor, each influencing factor is uniformly quantified as second, minute, and hour time period values, so as to measure the size of the panworkload in terms of the length of time worked. Then, the amount of panworkload within the cluster center is controlled within the normal working time. The panworkloadWi of cluster center Ei (the ith distribution route) can be expressed as
(3)Wi=α1Si+α2Ni+α3Qi+FT.Among them,Wiis the panworkload of cluster center Ei (the ith distribution route), Si is the total length of the distribution route, Niis the number of retail merchants on the distribution route, Qiis the cumulative volume of cigarettes distribution on the distribution route, and FT is the fixed amount of time that the delivery staff must spend on each distribution; it is mainly used for parking and warehousing, vehicle condition inspection, and procedure handover. α1, α2, and α3 are the weights before Si, Ni, and Qi, respectively.The operation are as follows: 100 customer node coordinates are randomly generated in MATLAB software to generate a scatter plot as well as the initial clustering centers. Figure1 shows the customer node scatter plot, Figure 2 shows the initial clustering center generated by the clustering algorithm, and Figure 3 shows the optimally adjusted region delineation map.Figure 1
The customer node scatter plot.Figure 2
The initial clustering center generated by the clustering algorithm.Figure 3
The optimally adjusted region delineation map.


### 2.3. Establishment of Distribution Route Optimization Model

#### 2.3.1. Problem Description

The route optimization problem in this paper can be described as a weighted graph:C=A,B, where A=0,N,M is the collection of logistics distribution nodes, 0 means transfer yard (distribution center), N=n1,n2,⋯,nn is the distribution center collection, M=1,2,⋯,m is the set of client nodes, B=bi,bj:bi,bj∈B,bi≠bj represents the set of vehicle route arcs and represents vehicle transportation routes. Set A contains multiple nonnegative attributes, such as customer’s delivery volume. K=k1,k2,⋯,kL is the collection of distribution vehicles, each car in the set K contains a nonnegative attribute. That is, the maximum carrying capacity Qmax, the delivery demand of customer i is qi. The problem can be described as that the distribution center sends vehicles of the same model with a load capacity Q to deliver the goods and then delivers the goods in turn to the customer nodes in need, it is required that all the distribution vehicles must start from the distribution center and return to the distribution center after completing the distribution task. Each customer node must be visited by a distribution vehicle and can only be visited once. The load capacity of the vehicle cannot exceed the maximum load limit of the vehicle. When a vehicle needs to be distributed across regions, the goods should be sent to the transit area for goods reorganization and then sent to the customer nodes in turn. The goal is to find the optimal distribution route to achieve the minimum total operating costs, including the distribution mileage cost, the fixed cost of vehicles and the fuel consumption cost of vehicles, and the minimum number of vehicles under the condition of meeting all customers’ demand for goods distribution and constraints.


#### 2.3.2. Penalty Function

In the actual process of cigarette distribution, the total cost of distribution refers to all the expenses incurred by the cigarette distribution center unilaterally in the distribution process and does not take into account the penalty cost of violating the customer time window constraints. If the distribution center violates the time window constraints of customers, it will inevitably bring some economic losses to cigarette retail merchants. Therefore, the service quality and market competitiveness of the distribution center will also suffer certain losses, and some potential customers may give up sales.Therefore, it is necessary for the distribution center to consider the time effect cost caused by these losses when pursuing the minimization of distribution cost and quantify the losses appropriately. The penalty cost may increase parabolically or exponentially with the quantity of goods, which varies according to the actual situation of goods and sales. This paper assumes that the penalty cost increases linearly to simplify the cost measurement problem.(4)citi=ceai−ti,ti<ai0,ai≤t≤biclti−bi,ti>bi.Among them,ti is the time when the vehicle arrives at customer i, ce represents the unit cost of delivery vehicles delivered earlier than the agreed time window (inventory cost, capital cost, etc.), cl represents the unit cost of delivery vehicles after the agreed time window (potential sales opportunities and economic profits). The opportunity cost of vehicle delivery before ai is ceai−ti and the penalty cost of vehicle delivery after biis clti−bi, and when the delivery vehicle completes the delivery in the agreed time window ai≤ti≤bi, the time utility cost is 0.


#### 2.3.3. Model Hypothesis

(1)
Suppose that all cigarette retailers are known setN, N=0,1,,2,⋯,n, M is a collection distribution centers M=1,2,⋯,m, 0 is the transition of cigarettes, and E is the distribution area obtained by cluster analysis(2)
The distribution vehicles start from the cigarette distribution center and return to the distribution center after the retail merchants of the route. (The order in which the distribution vehicles pass through the retailers is called the route)(3)
All the distribution vehicles in the cigarette distribution center are vehicles of the same model and the number of vehicles and the maximum carrying capacity of vehicles are known(4)
The geographic coordinates of all customer nodes, distribution centers, and transit sites in the network are given asGi and the latitude and longitude of the point are used to express the data(5)
The cigarette demand of each customer node is known(6)
Each customer node can only be accessed by one delivery vehicle and only once(7)
The distance between the customer node and the other customer node and between the customer node and the distribution center is the arc distance between them(8)
The distribution personnel carry out the distribution within normal working hours (8 hour working system)(9)
The delivery vehicles should be delivered within the time window agreed by customers; otherwise, they will be punished(10)
The delivery vehicles are delivered strictly according to the order quantity


#### 2.3.4. Parameter Symbol

c1 represents the cost per unit distance of the vehicle,c2 represents the fuel consumption parameter,c3 represents the fixed cost per vehicle,ti represents the actual delivery completion time of customer i,k is the number of the k cars,Qmax is the maximum load of the vehicle,qi represents the delivery demand of customer i,Qik represents the vehicle load when vehicle karrives at customer i,Sijk represents the distance of vehicle k from customer i to customer j,ω represents the vehicle fuel consumption coefficient,H represents a fixed constant,G represents the total number of distribution vehicles used,λ represents the constraint coefficient of the sub loop of the distribution route.


#### 2.3.5. Decision Variables

(5)xij=1vehiclekmovesfromitoj,kϵK0otherwiseymk1vehiclekisassignedtodistributioncenterm,mϵM0.


#### 2.3.6. Mathematical Model

The optimal distribution model of cigarette route is as follows:(6)minF=c1∑i,j∈N∑k∈KSijkxijkymk+c2∑i,j∈N∑k∈KωQijkSijkxijkymk+c3∑k∈K,j∈Nxijkymk+citis.t.(7)∑n∈Nynk=1,∀k∈K(8)∑i,j∈M∑k∈KQijk≤Qmax,i≠j(9)∑i,j∈N∑k∈Kxijk=1,i≠j(10)0≤nkm≤Nm(11)∑m∈MNm=H(12)Qmik=∑i,j∈N∑k∈K∑m∈Mqijkxijkymk,i≠j(13)∑i∈N∑k∈KQmik≤Qmax(14)0≤qi≤Qmax(15)G=∑k∈K∑i∈Nxmik(16)λj≥λi+1−n1−∑i,j∈N∑k∈Kxijk,i≠j(17)λj≥0,j∈N(18)xijk,yijk∈0,1∀i,j,k(19)citi=ce∑i=1Nmaxai−ti,0+cl∑i=1Nmaxti−bi,0In the above model, Equation (6) represents the optimal objective function of cigarette route optimization, in which the first term represents the total distance cost of vehicles, the second represents the total fuel consumption cost of the vehicles, the third represents the total fixed cost of the vehicles, and the fourth represents the penalty cost. Equation (7) indicates that each vehicle can only belong to one depot (distribution center). Equation (8) indicates that the carrying capacity of the vehicle passing through the customer node cannot be greater than the maximum carrying capacity of the vehicle. Equation (9) indicates that all customer nodes should be served, and each customer node can only be visited by one vehicle once. Equation (10) indicates that the number of customer nodes of a single route in each distribution area does not exceed the total number of customers in the distribution area. Equation (11) indicates that all customer nodes can be served by vehicles. Equation (12) shows that the carrying capacity of the vehicle departing from the distribution center is the sum of the distribution demand of the vehicle in the distribution route. Equations (13) and (14) indicate that the carrying capacity of any customer node of the vehicle cannot be greater than the maximum carrying capacity of the vehicle and the load of the vehicle cannot be negative. Equation (15) indicates the total number of vehicles used. Equations (16) and (17) indicate the subloop constraints of the distribution route, Equation (18) indicates that the decision scalars are all 0-1 variables. Equation (19) represents the penalty function expression obtained from (4).


## 2.1. Problem Definition and Model Framework Design

### 2.1.1. Problem Definition

In order to establish the cigarette distribution model and calculate it, we simplify the vehicle routing problem in view of the tobacco industry control policy and the particularity of cigarette distribution.Those are as follows:(1)
There is a time window(2)
Nonfull load distribution problem: the cargo demand of each customer is less than the vehicle loading capacity, and it can carry goods of multiple customers when the same vehicle is stowed(3)
Single objective problem: the objective is to minimize the total cost of cigarette distribution(4)
Pure delivery problem: the distribution center provides the distribution service that can meet the needs of customers to multiple customers (specification and quantity of goods)(5)
Single distribution vehicle model: the distance of distribution vehicle is within the maximum theoretical driving distance


### 2.1.2. Model Frame Design

Two-stage algorithm is a method which is widely used and has some achievements. The two-stage algorithm mainly includes path before grouping algorithm or group before path algorithm [20]. In this paper, we use the algorithm of grouping before routing, that is, the distribution area is generated first, and then, the optimal distribution route is obtained within this area.(1) First Stage: Distribution Areas’ Division. We determine the distribution center and distribution areas, transfer station (clustering center), and covered retail merchants according to the given distribution radiation radius. It is divided into different vehicle distribution areas according to the spatial geographical coordinates of retail merchants and different traffic conditions, such as transit area and direct distribution area. Clustering algorithm is used to divide or group customer nodes [17]. This makes the retail customers with relatively concentrated geographical distribution in the distribution area of the same car.(2) Second Stage: Distribution Route Optimization. The genetic algorithm is used to solve the TSP problem and the specific cigarette delivery routes or sequences of delivery vehicles within a single delivery area.Through the two-stage model algorithm to refine and decompose the vehicle routing problem of cigarette distribution, it can effectively improve the operational efficiency of the algorithm, shorten the calculation time, and get the approximate optimal solution.


## 2.1.1. Problem Definition

In order to establish the cigarette distribution model and calculate it, we simplify the vehicle routing problem in view of the tobacco industry control policy and the particularity of cigarette distribution.Those are as follows:(1)
There is a time window(2)
Nonfull load distribution problem: the cargo demand of each customer is less than the vehicle loading capacity, and it can carry goods of multiple customers when the same vehicle is stowed(3)
Single objective problem: the objective is to minimize the total cost of cigarette distribution(4)
Pure delivery problem: the distribution center provides the distribution service that can meet the needs of customers to multiple customers (specification and quantity of goods)(5)
Single distribution vehicle model: the distance of distribution vehicle is within the maximum theoretical driving distance


## 2.1.2. Model Frame Design

Two-stage algorithm is a method which is widely used and has some achievements. The two-stage algorithm mainly includes path before grouping algorithm or group before path algorithm [20]. In this paper, we use the algorithm of grouping before routing, that is, the distribution area is generated first, and then, the optimal distribution route is obtained within this area.(1) First Stage: Distribution Areas’ Division. We determine the distribution center and distribution areas, transfer station (clustering center), and covered retail merchants according to the given distribution radiation radius. It is divided into different vehicle distribution areas according to the spatial geographical coordinates of retail merchants and different traffic conditions, such as transit area and direct distribution area. Clustering algorithm is used to divide or group customer nodes [17]. This makes the retail customers with relatively concentrated geographical distribution in the distribution area of the same car.(2) Second Stage: Distribution Route Optimization. The genetic algorithm is used to solve the TSP problem and the specific cigarette delivery routes or sequences of delivery vehicles within a single delivery area.Through the two-stage model algorithm to refine and decompose the vehicle routing problem of cigarette distribution, it can effectively improve the operational efficiency of the algorithm, shorten the calculation time, and get the approximate optimal solution.


## 2.2. Distribution Area Division Clustering and Optimization

Cigarette distribution route optimization needs clustering analysis of spatial geographic data of retail households, and it needs a high efficiency algorithm due to the large amount of data. K-means algorithm is suitable for data type data and insensitive to data input sequence. So it is more appropriate.


### 2.2.1. Initial Clustering Center of K-Means

K-means clustering algorithm needs to solve two key problems, one is to determine the number of initial clustering centers, and the other is to determine the location of initial clustering centers.(1)
The number of initial cluster centers is the number of distribution areas, that is, the number of vehicles for cigarette distribution in these distribution areas. The number of distribution vehicles is set asm=cigarettedistributionvolume/singlevehicleload+1. Considering the maximum distance between nodes, and the fact that some vehicles cannot be fully loaded for distribution, so an additional vehicle is added for flexible configuration(2)
The location of the initial cluster center is usually randomly generated. In order to improve the clustering efficiency and reduce the amount of computation, we make some improvements to the generation of the initial clustering centers. Taking the geographical coordinates of each retail customer node as the center of the circle and the average distance between the geographical data of all retail merchants as the radius, the initial cluster center is determined according to the density of customer nodes in each circle [18]. Its implementation process mainly includes the following steps
(a)
Mark the geographical coordinates of each customer node in turn, and the mathematical distance between all nodes is calculated to generate the distance matrix. In order to simplify the model, it is assumed that the customer nodes are connected by a straight line, take the mathematical distance formula between two nodes as(1)Dij=xi−xj2+yi−yj2,i,j=1,2,⋯n(b)
Calculate the average distanceDij between all nodes in the database, let the initial radius of cluster center be R1=D, R2=2D(2)D=∑i=1n∑j=1nDijn×n(c)
Take each customer node as the center in turn, takeR1 as the initial radius to make a circle, and calculate the number of nodes in each circle, that is, sample density; for example, the sample density ρi of node Ai=xi,yi(d)
Sort from large to small according to sample density, the node with the highest sample density is taken as the first cluster centerE1, and the node with the second highest sample density is named E′. If the distance between E′and the first cluster center E1is E1E′≥R2, E′ is the cluster center E2. If the mathematical distance between the node and the previously determined cluster center is greater than R2; then, the node is the new cluster center. It iterates until m cluster centers are reached(e)
From the previous step, the numbers of cluster centers are mE1,E2,E3,⋯,andEm, as the initial clustering centers of K-meansStep 1.
set upρi=0Step 2.
Judge whether nodeA1=x1,y1 falls in the circle with Ai=xi,yi as the center and R1 as the radiu by A1Ai=xi−xj2+yi−yj2≤R1, if the above inequality is satisfied, then ρi=ρi+1Step 3.
Judge whether all nodes have been calculated, when the calculation is completed,ρiis the sample density of node Ai; otherwise, judge the next node A2, repeat the second step


### 2.2.2. Clustering Optimization considering Vehicle Load

The initial cluster centers only consider the density distribution of customer nodes and do not consider the specific demand of customer nodes. The different cigarette distribution volume of customer nodes will affect the optimization result of cigarette distribution route. Therefore, the initial clustering results should be further optimized and adjusted according to the constraints of single vehicle loading capacity. The specific optimization and adjustment plan are as follows:Step 1.
Calculate the number of retail merchants ism, and each cluster and their total cigarette distribution is Lm.Step 2.
The total cigarette distribution volumeLm of retail merchants in each cluster is compared with the maximum vehicle volume Tm of single vehicle. If Lm≤Tm, it means that the clustering meets the constraints of vehicle load. If Lm>Tm, it means that the distribution volume of cigarette in the cluster exceeds the vehicle volume. The customer node farthest from the cluster center in the cluster is excluded, and it is proposed to be classified as the cluster nearest to the adjacent cluster center. Before the cluster is classified, it is necessary to judge the updated Lm≤Tm of the cluster. If Lm≤Tm is true, the node will be classified into the cluster. Otherwise, the cluster with poor Euclidean distance will be selected to judge Lm≤Tm. If the condition of constrained single vehicle loading capacity is met, the node will be classified into the cluster. Otherwise, the node will be searched and judged againStep 3.
When the total distribution volume of all cluster retailers is less than the single vehicle loading capacityLm≤Tm, the adjustment is completed.


### 2.2.3. Clustering Optimization considering Workload

In order to avoid the local optimal situation, we should not only consider the constraints of vehicle capacity but also consider the cumulative workload within the distribution area. In the actual distribution work, the average daily working hours of distribution workers should be kept within the normal working hours of 8 hours.To calculate and evaluate the workload of cigarette distribution staff scientifically and reasonably has always been a difficult problem in cigarette distribution management. For one, the external and internal factors affecting the workload of distribution staff are numerous and complex. Extrinsic factors include distribution workload, distribution distance, road conditions, and the dispersion of retail nodes; intrinsic factors include staff’s work attitude and experience, the rationality of the division of labor, and the comprehensive quality of distribution staff, etc. For another, the complexity of the influencing factors and the lack of uniformity in the magnitude of the various influencing factors do not facilitate calculation and judgement.Therefore, some experts and scholars have proposed the concept of “panworkload,” which is used to calculate and evaluate the workload of each distribution route and to adjust the clustering results with the workload as the constraint condition. The concept of “panworkload” comprehensively considers the number of customer nodes, the volume of cigarettes distribution, the mileage of delivery, and specific traffic conditions. By assigning different weights to each influencing factor, each influencing factor is uniformly quantified as second, minute, and hour time period values, so as to measure the size of the panworkload in terms of the length of time worked. Then, the amount of panworkload within the cluster center is controlled within the normal working time. The panworkloadWi of cluster center Ei (the ith distribution route) can be expressed as
(3)Wi=α1Si+α2Ni+α3Qi+FT.Among them,Wiis the panworkload of cluster center Ei (the ith distribution route), Si is the total length of the distribution route, Niis the number of retail merchants on the distribution route, Qiis the cumulative volume of cigarettes distribution on the distribution route, and FT is the fixed amount of time that the delivery staff must spend on each distribution; it is mainly used for parking and warehousing, vehicle condition inspection, and procedure handover. α1, α2, and α3 are the weights before Si, Ni, and Qi, respectively.The operation are as follows: 100 customer node coordinates are randomly generated in MATLAB software to generate a scatter plot as well as the initial clustering centers. Figure1 shows the customer node scatter plot, Figure 2 shows the initial clustering center generated by the clustering algorithm, and Figure 3 shows the optimally adjusted region delineation map.Figure 1
The customer node scatter plot.Figure 2
The initial clustering center generated by the clustering algorithm.Figure 3
The optimally adjusted region delineation map.


## 2.2.1. Initial Clustering Center of K-Means

K-means clustering algorithm needs to solve two key problems, one is to determine the number of initial clustering centers, and the other is to determine the location of initial clustering centers.(1)
The number of initial cluster centers is the number of distribution areas, that is, the number of vehicles for cigarette distribution in these distribution areas. The number of distribution vehicles is set asm=cigarettedistributionvolume/singlevehicleload+1. Considering the maximum distance between nodes, and the fact that some vehicles cannot be fully loaded for distribution, so an additional vehicle is added for flexible configuration(2)
The location of the initial cluster center is usually randomly generated. In order to improve the clustering efficiency and reduce the amount of computation, we make some improvements to the generation of the initial clustering centers. Taking the geographical coordinates of each retail customer node as the center of the circle and the average distance between the geographical data of all retail merchants as the radius, the initial cluster center is determined according to the density of customer nodes in each circle [18]. Its implementation process mainly includes the following steps
(a)
Mark the geographical coordinates of each customer node in turn, and the mathematical distance between all nodes is calculated to generate the distance matrix. In order to simplify the model, it is assumed that the customer nodes are connected by a straight line, take the mathematical distance formula between two nodes as(1)Dij=xi−xj2+yi−yj2,i,j=1,2,⋯n(b)
Calculate the average distanceDij between all nodes in the database, let the initial radius of cluster center be R1=D, R2=2D(2)D=∑i=1n∑j=1nDijn×n(c)
Take each customer node as the center in turn, takeR1 as the initial radius to make a circle, and calculate the number of nodes in each circle, that is, sample density; for example, the sample density ρi of node Ai=xi,yi(d)
Sort from large to small according to sample density, the node with the highest sample density is taken as the first cluster centerE1, and the node with the second highest sample density is named E′. If the distance between E′and the first cluster center E1is E1E′≥R2, E′ is the cluster center E2. If the mathematical distance between the node and the previously determined cluster center is greater than R2; then, the node is the new cluster center. It iterates until m cluster centers are reached(e)
From the previous step, the numbers of cluster centers are mE1,E2,E3,⋯,andEm, as the initial clustering centers of K-meansStep 1.
set upρi=0Step 2.
Judge whether nodeA1=x1,y1 falls in the circle with Ai=xi,yi as the center and R1 as the radiu by A1Ai=xi−xj2+yi−yj2≤R1, if the above inequality is satisfied, then ρi=ρi+1Step 3.
Judge whether all nodes have been calculated, when the calculation is completed,ρiis the sample density of node Ai; otherwise, judge the next node A2, repeat the second step


## 2.2.2. Clustering Optimization considering Vehicle Load

The initial cluster centers only consider the density distribution of customer nodes and do not consider the specific demand of customer nodes. The different cigarette distribution volume of customer nodes will affect the optimization result of cigarette distribution route. Therefore, the initial clustering results should be further optimized and adjusted according to the constraints of single vehicle loading capacity. The specific optimization and adjustment plan are as follows:Step 1.
Calculate the number of retail merchants ism, and each cluster and their total cigarette distribution is Lm.Step 2.
The total cigarette distribution volumeLm of retail merchants in each cluster is compared with the maximum vehicle volume Tm of single vehicle. If Lm≤Tm, it means that the clustering meets the constraints of vehicle load. If Lm>Tm, it means that the distribution volume of cigarette in the cluster exceeds the vehicle volume. The customer node farthest from the cluster center in the cluster is excluded, and it is proposed to be classified as the cluster nearest to the adjacent cluster center. Before the cluster is classified, it is necessary to judge the updated Lm≤Tm of the cluster. If Lm≤Tm is true, the node will be classified into the cluster. Otherwise, the cluster with poor Euclidean distance will be selected to judge Lm≤Tm. If the condition of constrained single vehicle loading capacity is met, the node will be classified into the cluster. Otherwise, the node will be searched and judged againStep 3.
When the total distribution volume of all cluster retailers is less than the single vehicle loading capacityLm≤Tm, the adjustment is completed.


## 2.2.3. Clustering Optimization considering Workload

In order to avoid the local optimal situation, we should not only consider the constraints of vehicle capacity but also consider the cumulative workload within the distribution area. In the actual distribution work, the average daily working hours of distribution workers should be kept within the normal working hours of 8 hours.To calculate and evaluate the workload of cigarette distribution staff scientifically and reasonably has always been a difficult problem in cigarette distribution management. For one, the external and internal factors affecting the workload of distribution staff are numerous and complex. Extrinsic factors include distribution workload, distribution distance, road conditions, and the dispersion of retail nodes; intrinsic factors include staff’s work attitude and experience, the rationality of the division of labor, and the comprehensive quality of distribution staff, etc. For another, the complexity of the influencing factors and the lack of uniformity in the magnitude of the various influencing factors do not facilitate calculation and judgement.Therefore, some experts and scholars have proposed the concept of “panworkload,” which is used to calculate and evaluate the workload of each distribution route and to adjust the clustering results with the workload as the constraint condition. The concept of “panworkload” comprehensively considers the number of customer nodes, the volume of cigarettes distribution, the mileage of delivery, and specific traffic conditions. By assigning different weights to each influencing factor, each influencing factor is uniformly quantified as second, minute, and hour time period values, so as to measure the size of the panworkload in terms of the length of time worked. Then, the amount of panworkload within the cluster center is controlled within the normal working time. The panworkloadWi of cluster center Ei (the ith distribution route) can be expressed as
(3)Wi=α1Si+α2Ni+α3Qi+FT.Among them,Wiis the panworkload of cluster center Ei (the ith distribution route), Si is the total length of the distribution route, Niis the number of retail merchants on the distribution route, Qiis the cumulative volume of cigarettes distribution on the distribution route, and FT is the fixed amount of time that the delivery staff must spend on each distribution; it is mainly used for parking and warehousing, vehicle condition inspection, and procedure handover. α1, α2, and α3 are the weights before Si, Ni, and Qi, respectively.The operation are as follows: 100 customer node coordinates are randomly generated in MATLAB software to generate a scatter plot as well as the initial clustering centers. Figure1 shows the customer node scatter plot, Figure 2 shows the initial clustering center generated by the clustering algorithm, and Figure 3 shows the optimally adjusted region delineation map.Figure 1
The customer node scatter plot.Figure 2
The initial clustering center generated by the clustering algorithm.Figure 3
The optimally adjusted region delineation map.


## 2.3. Establishment of Distribution Route Optimization Model

### 2.3.1. Problem Description

The route optimization problem in this paper can be described as a weighted graph:C=A,B, where A=0,N,M is the collection of logistics distribution nodes, 0 means transfer yard (distribution center), N=n1,n2,⋯,nn is the distribution center collection, M=1,2,⋯,m is the set of client nodes, B=bi,bj:bi,bj∈B,bi≠bj represents the set of vehicle route arcs and represents vehicle transportation routes. Set A contains multiple nonnegative attributes, such as customer’s delivery volume. K=k1,k2,⋯,kL is the collection of distribution vehicles, each car in the set K contains a nonnegative attribute. That is, the maximum carrying capacity Qmax, the delivery demand of customer i is qi. The problem can be described as that the distribution center sends vehicles of the same model with a load capacity Q to deliver the goods and then delivers the goods in turn to the customer nodes in need, it is required that all the distribution vehicles must start from the distribution center and return to the distribution center after completing the distribution task. Each customer node must be visited by a distribution vehicle and can only be visited once. The load capacity of the vehicle cannot exceed the maximum load limit of the vehicle. When a vehicle needs to be distributed across regions, the goods should be sent to the transit area for goods reorganization and then sent to the customer nodes in turn. The goal is to find the optimal distribution route to achieve the minimum total operating costs, including the distribution mileage cost, the fixed cost of vehicles and the fuel consumption cost of vehicles, and the minimum number of vehicles under the condition of meeting all customers’ demand for goods distribution and constraints.


### 2.3.2. Penalty Function

In the actual process of cigarette distribution, the total cost of distribution refers to all the expenses incurred by the cigarette distribution center unilaterally in the distribution process and does not take into account the penalty cost of violating the customer time window constraints. If the distribution center violates the time window constraints of customers, it will inevitably bring some economic losses to cigarette retail merchants. Therefore, the service quality and market competitiveness of the distribution center will also suffer certain losses, and some potential customers may give up sales.Therefore, it is necessary for the distribution center to consider the time effect cost caused by these losses when pursuing the minimization of distribution cost and quantify the losses appropriately. The penalty cost may increase parabolically or exponentially with the quantity of goods, which varies according to the actual situation of goods and sales. This paper assumes that the penalty cost increases linearly to simplify the cost measurement problem.(4)citi=ceai−ti,ti<ai0,ai≤t≤biclti−bi,ti>bi.Among them,ti is the time when the vehicle arrives at customer i, ce represents the unit cost of delivery vehicles delivered earlier than the agreed time window (inventory cost, capital cost, etc.), cl represents the unit cost of delivery vehicles after the agreed time window (potential sales opportunities and economic profits). The opportunity cost of vehicle delivery before ai is ceai−ti and the penalty cost of vehicle delivery after biis clti−bi, and when the delivery vehicle completes the delivery in the agreed time window ai≤ti≤bi, the time utility cost is 0.


### 2.3.3. Model Hypothesis

(1)
Suppose that all cigarette retailers are known setN, N=0,1,,2,⋯,n, M is a collection distribution centers M=1,2,⋯,m, 0 is the transition of cigarettes, and E is the distribution area obtained by cluster analysis(2)
The distribution vehicles start from the cigarette distribution center and return to the distribution center after the retail merchants of the route. (The order in which the distribution vehicles pass through the retailers is called the route)(3)
All the distribution vehicles in the cigarette distribution center are vehicles of the same model and the number of vehicles and the maximum carrying capacity of vehicles are known(4)
The geographic coordinates of all customer nodes, distribution centers, and transit sites in the network are given asGi and the latitude and longitude of the point are used to express the data(5)
The cigarette demand of each customer node is known(6)
Each customer node can only be accessed by one delivery vehicle and only once(7)
The distance between the customer node and the other customer node and between the customer node and the distribution center is the arc distance between them(8)
The distribution personnel carry out the distribution within normal working hours (8 hour working system)(9)
The delivery vehicles should be delivered within the time window agreed by customers; otherwise, they will be punished(10)
The delivery vehicles are delivered strictly according to the order quantity


### 2.3.4. Parameter Symbol

c1 represents the cost per unit distance of the vehicle,c2 represents the fuel consumption parameter,c3 represents the fixed cost per vehicle,ti represents the actual delivery completion time of customer i,k is the number of the k cars,Qmax is the maximum load of the vehicle,qi represents the delivery demand of customer i,Qik represents the vehicle load when vehicle karrives at customer i,Sijk represents the distance of vehicle k from customer i to customer j,ω represents the vehicle fuel consumption coefficient,H represents a fixed constant,G represents the total number of distribution vehicles used,λ represents the constraint coefficient of the sub loop of the distribution route.


### 2.3.5. Decision Variables

(5)xij=1vehiclekmovesfromitoj,kϵK0otherwiseymk1vehiclekisassignedtodistributioncenterm,mϵM0.


### 2.3.6. Mathematical Model

The optimal distribution model of cigarette route is as follows:(6)minF=c1∑i,j∈N∑k∈KSijkxijkymk+c2∑i,j∈N∑k∈KωQijkSijkxijkymk+c3∑k∈K,j∈Nxijkymk+citis.t.(7)∑n∈Nynk=1,∀k∈K(8)∑i,j∈M∑k∈KQijk≤Qmax,i≠j(9)∑i,j∈N∑k∈Kxijk=1,i≠j(10)0≤nkm≤Nm(11)∑m∈MNm=H(12)Qmik=∑i,j∈N∑k∈K∑m∈Mqijkxijkymk,i≠j(13)∑i∈N∑k∈KQmik≤Qmax(14)0≤qi≤Qmax(15)G=∑k∈K∑i∈Nxmik(16)λj≥λi+1−n1−∑i,j∈N∑k∈Kxijk,i≠j(17)λj≥0,j∈N(18)xijk,yijk∈0,1∀i,j,k(19)citi=ce∑i=1Nmaxai−ti,0+cl∑i=1Nmaxti−bi,0In the above model, Equation (6) represents the optimal objective function of cigarette route optimization, in which the first term represents the total distance cost of vehicles, the second represents the total fuel consumption cost of the vehicles, the third represents the total fixed cost of the vehicles, and the fourth represents the penalty cost. Equation (7) indicates that each vehicle can only belong to one depot (distribution center). Equation (8) indicates that the carrying capacity of the vehicle passing through the customer node cannot be greater than the maximum carrying capacity of the vehicle. Equation (9) indicates that all customer nodes should be served, and each customer node can only be visited by one vehicle once. Equation (10) indicates that the number of customer nodes of a single route in each distribution area does not exceed the total number of customers in the distribution area. Equation (11) indicates that all customer nodes can be served by vehicles. Equation (12) shows that the carrying capacity of the vehicle departing from the distribution center is the sum of the distribution demand of the vehicle in the distribution route. Equations (13) and (14) indicate that the carrying capacity of any customer node of the vehicle cannot be greater than the maximum carrying capacity of the vehicle and the load of the vehicle cannot be negative. Equation (15) indicates the total number of vehicles used. Equations (16) and (17) indicate the subloop constraints of the distribution route, Equation (18) indicates that the decision scalars are all 0-1 variables. Equation (19) represents the penalty function expression obtained from (4).


## 2.3.1. Problem Description

The route optimization problem in this paper can be described as a weighted graph:C=A,B, where A=0,N,M is the collection of logistics distribution nodes, 0 means transfer yard (distribution center), N=n1,n2,⋯,nn is the distribution center collection, M=1,2,⋯,m is the set of client nodes, B=bi,bj:bi,bj∈B,bi≠bj represents the set of vehicle route arcs and represents vehicle transportation routes. Set A contains multiple nonnegative attributes, such as customer’s delivery volume. K=k1,k2,⋯,kL is the collection of distribution vehicles, each car in the set K contains a nonnegative attribute. That is, the maximum carrying capacity Qmax, the delivery demand of customer i is qi. The problem can be described as that the distribution center sends vehicles of the same model with a load capacity Q to deliver the goods and then delivers the goods in turn to the customer nodes in need, it is required that all the distribution vehicles must start from the distribution center and return to the distribution center after completing the distribution task. Each customer node must be visited by a distribution vehicle and can only be visited once. The load capacity of the vehicle cannot exceed the maximum load limit of the vehicle. When a vehicle needs to be distributed across regions, the goods should be sent to the transit area for goods reorganization and then sent to the customer nodes in turn. The goal is to find the optimal distribution route to achieve the minimum total operating costs, including the distribution mileage cost, the fixed cost of vehicles and the fuel consumption cost of vehicles, and the minimum number of vehicles under the condition of meeting all customers’ demand for goods distribution and constraints.


## 2.3.2. Penalty Function

In the actual process of cigarette distribution, the total cost of distribution refers to all the expenses incurred by the cigarette distribution center unilaterally in the distribution process and does not take into account the penalty cost of violating the customer time window constraints. If the distribution center violates the time window constraints of customers, it will inevitably bring some economic losses to cigarette retail merchants. Therefore, the service quality and market competitiveness of the distribution center will also suffer certain losses, and some potential customers may give up sales.Therefore, it is necessary for the distribution center to consider the time effect cost caused by these losses when pursuing the minimization of distribution cost and quantify the losses appropriately. The penalty cost may increase parabolically or exponentially with the quantity of goods, which varies according to the actual situation of goods and sales. This paper assumes that the penalty cost increases linearly to simplify the cost measurement problem.(4)citi=ceai−ti,ti<ai0,ai≤t≤biclti−bi,ti>bi.Among them,ti is the time when the vehicle arrives at customer i, ce represents the unit cost of delivery vehicles delivered earlier than the agreed time window (inventory cost, capital cost, etc.), cl represents the unit cost of delivery vehicles after the agreed time window (potential sales opportunities and economic profits). The opportunity cost of vehicle delivery before ai is ceai−ti and the penalty cost of vehicle delivery after biis clti−bi, and when the delivery vehicle completes the delivery in the agreed time window ai≤ti≤bi, the time utility cost is 0.


## 2.3.3. Model Hypothesis

(1)
Suppose that all cigarette retailers are known setN, N=0,1,,2,⋯,n, M is a collection distribution centers M=1,2,⋯,m, 0 is the transition of cigarettes, and E is the distribution area obtained by cluster analysis(2)
The distribution vehicles start from the cigarette distribution center and return to the distribution center after the retail merchants of the route. (The order in which the distribution vehicles pass through the retailers is called the route)(3)
All the distribution vehicles in the cigarette distribution center are vehicles of the same model and the number of vehicles and the maximum carrying capacity of vehicles are known(4)
The geographic coordinates of all customer nodes, distribution centers, and transit sites in the network are given asGi and the latitude and longitude of the point are used to express the data(5)
The cigarette demand of each customer node is known(6)
Each customer node can only be accessed by one delivery vehicle and only once(7)
The distance between the customer node and the other customer node and between the customer node and the distribution center is the arc distance between them(8)
The distribution personnel carry out the distribution within normal working hours (8 hour working system)(9)
The delivery vehicles should be delivered within the time window agreed by customers; otherwise, they will be punished(10)
The delivery vehicles are delivered strictly according to the order quantity


## 2.3.4. Parameter Symbol

c1 represents the cost per unit distance of the vehicle,c2 represents the fuel consumption parameter,c3 represents the fixed cost per vehicle,ti represents the actual delivery completion time of customer i,k is the number of the k cars,Qmax is the maximum load of the vehicle,qi represents the delivery demand of customer i,Qik represents the vehicle load when vehicle karrives at customer i,Sijk represents the distance of vehicle k from customer i to customer j,ω represents the vehicle fuel consumption coefficient,H represents a fixed constant,G represents the total number of distribution vehicles used,λ represents the constraint coefficient of the sub loop of the distribution route.


## 2.3.5. Decision Variables

(5)xij=1vehiclekmovesfromitoj,kϵK0otherwiseymk1vehiclekisassignedtodistributioncenterm,mϵM0.


## 2.3.6. Mathematical Model

The optimal distribution model of cigarette route is as follows:(6)minF=c1∑i,j∈N∑k∈KSijkxijkymk+c2∑i,j∈N∑k∈KωQijkSijkxijkymk+c3∑k∈K,j∈Nxijkymk+citis.t.(7)∑n∈Nynk=1,∀k∈K(8)∑i,j∈M∑k∈KQijk≤Qmax,i≠j(9)∑i,j∈N∑k∈Kxijk=1,i≠j(10)0≤nkm≤Nm(11)∑m∈MNm=H(12)Qmik=∑i,j∈N∑k∈K∑m∈Mqijkxijkymk,i≠j(13)∑i∈N∑k∈KQmik≤Qmax(14)0≤qi≤Qmax(15)G=∑k∈K∑i∈Nxmik(16)λj≥λi+1−n1−∑i,j∈N∑k∈Kxijk,i≠j(17)λj≥0,j∈N(18)xijk,yijk∈0,1∀i,j,k(19)citi=ce∑i=1Nmaxai−ti,0+cl∑i=1Nmaxti−bi,0In the above model, Equation (6) represents the optimal objective function of cigarette route optimization, in which the first term represents the total distance cost of vehicles, the second represents the total fuel consumption cost of the vehicles, the third represents the total fixed cost of the vehicles, and the fourth represents the penalty cost. Equation (7) indicates that each vehicle can only belong to one depot (distribution center). Equation (8) indicates that the carrying capacity of the vehicle passing through the customer node cannot be greater than the maximum carrying capacity of the vehicle. Equation (9) indicates that all customer nodes should be served, and each customer node can only be visited by one vehicle once. Equation (10) indicates that the number of customer nodes of a single route in each distribution area does not exceed the total number of customers in the distribution area. Equation (11) indicates that all customer nodes can be served by vehicles. Equation (12) shows that the carrying capacity of the vehicle departing from the distribution center is the sum of the distribution demand of the vehicle in the distribution route. Equations (13) and (14) indicate that the carrying capacity of any customer node of the vehicle cannot be greater than the maximum carrying capacity of the vehicle and the load of the vehicle cannot be negative. Equation (15) indicates the total number of vehicles used. Equations (16) and (17) indicate the subloop constraints of the distribution route, Equation (18) indicates that the decision scalars are all 0-1 variables. Equation (19) represents the penalty function expression obtained from (4).


## 3. Solution and Algorithm Design of Cigarette Distribution Route Optimization Model

### 3.1. Design of Clustering Algorithm for Distribution Area Division

Taking the highest single vehicle loading rate (or the least distribution vehicles) as the optimization objective and the vehicle load and daily workload as the constraints, the clustering algorithm is used to realize the solution and divide the distribution area.In order to calculate and evaluate the workload of distribution personnel in a scientific and reasonable way, we can follow the distribution vehicles, and by extensive contact and investigation with the distribution center management personnel and distribution personnel, we can understand and be familiar with the actual work of distribution personnel. In order to calculate and evaluate the workload of distribution personnel more comprehensively and reasonably, it is necessary to modify and adjust the specific connotation of the workload. Through the statistics and analysis of the historical data of cigarette distribution and following the distribution vehicles, we can take approximately exact values forα1, α2, and α3.
(1)
To determine the coefficientα1, we have to consider two factors: the speed of distribution vehicles v, the other is the distribution road condition, and β is used to express the traffic condition coefficient. The worse the road condition is, the smaller the β, and the longer the delivery time. In the actual distribution process, the vehicle speed is uncertain, so it can only take the approximate average value. At the same time, according to the distribution of urban network and rural network, the road conditions are divided into two levels (or other categories). And according to the difference between urban network and rural network, different α1 values can be obtained. For example, the driving speed of IVECO is v=50 km/h in rural distribution area and v=30 km/h in urban distribution area. The specific value of α1 coefficient is shown in Table 1(20)α1=1βv(2)
Determination ofFTvalue,FTis the fixed time that the delivery personnel must spend in each delivery. It is mainly used for parking and warehousing, vehicle condition inspection, and procedure handover. After field investigation and statistics, FT=20minutesTable 1
α1 coefficient value table.
Distribution routeRural networkUrban networkTraffic condition coefficientβ0.750.5Iveco driving speedv5030α1 coefficient value0.0270.067After determining the above weight coefficients and fixing the time, the general workload of each distribution route is counted to ensure that the workload of each cluster center is roughly balanced and the average daily working time is within the range of normal working hours.


### 3.2. Design and Improvement of Adaptive Genetic Algorithm

Although the accurate algorithm can get a satisfactory solution, the calculation is too large, so it can only be applied to small-scale and simple vehicle routing problems. For the current vehicle routing problems, it has almost no practical value. Compared with the precise algorithm, the classical heuristic algorithm has short computation time and can get the initial feasible solution quickly. However, it is limited to small and medium-sized vehicle routing problems and is often used for local optimization of existing paths. Modern heuristic algorithm not only has stable performance, fast calculation speed, but also has flexible and powerful search ability, so it can solve large-scale complex vehicle routing problems. Modern heuristic algorithm becomes the main research direction for solving vehicle routing problem, such as genetic algorithm is used in many literature, due to its good solving performance, and is widely used by experts and scholars at home and abroad to solve vehicle path problems. However, with the deepening of research, the single modern heuristic algorithm cannot meet the needs of research gradually, and the hybrid algorithm arises at the historic moment. The hybrid algorithm using multiple search mechanisms improves the shortcoming of single modern heuristic algorithm, overcomes the limitation of single heuristic algorithm effectively, and improves the solving speed and quality of the algorithm.Compared with other algorithms, genetic algorithm indeed has many excellent properties, which provide a good way to solve vehicle routing problems. With the deepening of the research, some defects are also exposed, such as premature phenomenon. But genetic algorithm is parallel, so it can be well combined with other algorithms, such as scan algorithm, tabu search algorithm, simulated annealing, and immune algorithm, so as to make up for these defects and form a better hybrid algorithm.We design an improved adaptive genetic algorithm that in view of the above defects of genetic algorithm. Genetic algorithm can be improved in coding structure, crossover, and mutation operation and population. We use natural number coding, design adaptive crossover and mutation probabilities, design reversal operator to correct errors, and optimize the initial population with scanning algorithm, so as to reduce the premature defects of genetic algorithm and improve the performance of the algorithm. The specific design is discussed in the following sections of this chapter.


#### 3.2.1. Chromosome Coding Design

Due to the special structure of the solution of the vehicle routing problem in cigarette distribution, a more intuitive natural number coding method is adopted in the chromosome coding design. Each chromosome is composed of the serial number of the depot and the customer point, and the customers and the depot are arranged together. If an existing distribution network composed ofn customers’ needs m vehicles for distribution, and the vehicles start from the depot and do not return to the depot after completing the task, the chromosome code can be expressed as follows:
(21)0,s11,s12⋯,s1p,0,s21,s22,⋯,s2q,0,⋯,0,sm1,sm1,⋯,smw.The length of chromosome isn+m, where skj is the jth customer served by the kth car, and m zero numbers represent the depot, which divides the chromosome into m segments, that is, there are m paths in total.


#### 3.2.2. Population Initialization

Starting from a better initial population, genetic algorithm can find a better solution faster, so it is necessary to find a high-quality initial population as the initial solution after completing the chromosome coding. The scanning algorithm is used to generate the initial population. The steps are as follows:Step 1.
Taking the depot as the pole and the connecting line between the depot and any customer as the polar axis, the polar coordinate system is established and the polar coordinates of each customer point are calculatedStep 2.
Scan the customer points’ clockwise. Check the constraint condition once for each additional customer. When the constraint condition is violated, eliminate the customer and insert 0 at the beginning of the path. Then continue scanning until all customers are scanned. The customer sequence obtained after scanning is a chromosome. For multiple customers at the same angle, the order is from near to far, until the initial population with n chromosomes is generated


#### 3.2.3. Fitness Function

In general, the higher the fitness value, the better the performance of the chromosome, and the greater the probability of entering the next generation, and vice versa. Therefore, the fitness function can be set as the reciprocal of the objective function. If the objective function value of chromosomei is zi, then the fitness value of chromosome i is fi=1/zi.


#### 3.2.4. Genetic Operator

(1) Selection Operator. In order to complete the genetic evolution operation, we need to select a certain number of individuals from the initial population to form a new population. In this paper, we use the improved roulette selection method as the selection operator. The specific steps are as follows:Step 1.
The fitness values of individuals were ranked to ensure that the top m individuals with high fitness entered the next generationStep 2.
Roulette is used to select and copy the remaining individuals. The probability is calculated by formulaPi=c1−ci−1, where c is the selection probability of the ordered individuals, and the cumulative probability is calculated by formula Qi=∑i=1nPiStep 3.
Randomly generated numberr between [0, 1] is used as the selection pointer to determine the selected individual. If r≤Qi, individual 1 is selected; if Qi−1≤r≤Qi, individual i is selected;Step 4.
Repeat step 3 until (n-m) individuals are generatedIn order to ensure the genetic evolution operation and the diversity of the population, the crossover and mutation operation adopt the adaptive adjustment strategy. The functions of adaptive crossover probabilitypc and mutation probability pm are as follows:
(22)pc=pc1f−favgfmax−favg,f>favgpc2,f≤favgpm=pm1f′−favgfmax−favgf′>favgpm2f′≤favgAmong them,f is the larger fitness function value of the two crossed chromosomes; fmax is the maximum fitness function value of the population; favg is the average fitness function value of the population; f′ is the fitness function value of the chromosome to be mutated; Pc1 and Pc2 are in the [0.8, 1] interval; and pm1 and pm2 are in the [0.001, 0.1] interval.(2) Crossover Operator. Crossover operation is very important in genetic algorithm, which is the process of exchanging some genes of paired chromosomes in some way to form two new chromosomes. In order to avoid destroying excellent genes, the substring is moved to the first place when crossing, instead of directly copying in the past. This can not only protect excellent substrings but also avoid the generation of new individuals when the parents are the same, which improves the optimization ability and convergence speed of the algorithm. Crossover operation is the process of exchanging some genes of paired chromosomes in some way to form two new chromosomes, which is very important in genetic algorithms. In order to avoid destroying excellent genes, it is necessary to move the substrings to the first place during the crossover operation instead of copying the substrings directly over. This can protect the good substrings and avoids producing new individuals when the parents are identical, which improves the algorithm’s optimization ability and convergence speed. The specific operations are as follows:Step 1.
The adaptive crossover probabilitypc is calculated and the parent chromosomes A and B to be crossed are selected.Step 2.
In the interval1,k, integers r1 and r2 are generated randomly, and the sub path r1 in A is copied to the front end of the child B*, and the sub path r2 in B is copied to the front end of the child A*.Step 3.
The remaining k-1 zeros are randomly inserted into the remaining spaces, and the last space is not zero, and any two zeros cannot be adjacent.Step 4.
The remaining genes in the parent chromosomeA were inserted into the offspring A* to generate the complete offspring chromosome A*. Similarly, the progeny chromosome B* was generated.Step 5.
The decoding operation is performed to check whether the constraint conditions are met, and if not, discard.(3) Mutation Operator. In order to avoid losing some useful genes and prevent some chromosomes from being constant, mutation operation changes the position of some genes in the operation process, which can effectively prevent the algorithm from premature. In this adaptive mutation operation, the exchange mutation mode is adopted, and the specific operation is as follows:Step 1.
The adaptive mutation probabilitypm was calculated to select chromosomes.Step 2.
Two unequal integerst1 and t2 are randomly generated in the range of 1,n+k.Step 3.
Ift1 and t2 are both non-zero genes, the position will be exchanged, otherwise it will be generated again.Step 4.
Check whether the new offspring chromosome meets the constraints, if not, it will be generated again.Take chromosome 056104780392 as an example, ift1=2 and t2=6, and then, the chromosome generated by the above operator is 046105780392, as shown in Figure 1.(4) Evolution Reversal Operation. In order to correct the error in the process of evolution and improve the local search ability of the genetic algorithm, a gene segment is randomly selected from the parent generation and reversed in one direction. The specific operation is as follows:Step 1.
In interval [1,n+k], s1 and s2 are generated randomly.Step 2.
Ifs1 and s2 are both nonzero genes, the client genes between them will be reversed, otherwise they will be generated again.Step 3.
Check whether the offspring chromosomes meet the constraints, and whether the fitness value is improved. If not, discard them.Take chromosome 056104780392 as an example, ifs1=2, s2=6, and then using the above operator, the generated chromosome is 036105870429, as shown in Figure 2.


#### 3.2.5. Termination Conditions

We set appropriate termination conditions to avoid the endless running of the algorithm or a dead loop without knowing. When an optimization result appears several times or reaches the maximum number of iterations, the operation is terminated and we output the satisfactory solution or approximate optimal solution. In this paper, the maximum number of iterations is 1000 as the termination condition.


## 3.1. Design of Clustering Algorithm for Distribution Area Division

Taking the highest single vehicle loading rate (or the least distribution vehicles) as the optimization objective and the vehicle load and daily workload as the constraints, the clustering algorithm is used to realize the solution and divide the distribution area.In order to calculate and evaluate the workload of distribution personnel in a scientific and reasonable way, we can follow the distribution vehicles, and by extensive contact and investigation with the distribution center management personnel and distribution personnel, we can understand and be familiar with the actual work of distribution personnel. In order to calculate and evaluate the workload of distribution personnel more comprehensively and reasonably, it is necessary to modify and adjust the specific connotation of the workload. Through the statistics and analysis of the historical data of cigarette distribution and following the distribution vehicles, we can take approximately exact values forα1, α2, and α3.
(1)
To determine the coefficientα1, we have to consider two factors: the speed of distribution vehicles v, the other is the distribution road condition, and β is used to express the traffic condition coefficient. The worse the road condition is, the smaller the β, and the longer the delivery time. In the actual distribution process, the vehicle speed is uncertain, so it can only take the approximate average value. At the same time, according to the distribution of urban network and rural network, the road conditions are divided into two levels (or other categories). And according to the difference between urban network and rural network, different α1 values can be obtained. For example, the driving speed of IVECO is v=50 km/h in rural distribution area and v=30 km/h in urban distribution area. The specific value of α1 coefficient is shown in Table 1(20)α1=1βv(2)
Determination ofFTvalue,FTis the fixed time that the delivery personnel must spend in each delivery. It is mainly used for parking and warehousing, vehicle condition inspection, and procedure handover. After field investigation and statistics, FT=20minutesTable 1
α1 coefficient value table.
Distribution routeRural networkUrban networkTraffic condition coefficientβ0.750.5Iveco driving speedv5030α1 coefficient value0.0270.067After determining the above weight coefficients and fixing the time, the general workload of each distribution route is counted to ensure that the workload of each cluster center is roughly balanced and the average daily working time is within the range of normal working hours.


## 3.2. Design and Improvement of Adaptive Genetic Algorithm

Although the accurate algorithm can get a satisfactory solution, the calculation is too large, so it can only be applied to small-scale and simple vehicle routing problems. For the current vehicle routing problems, it has almost no practical value. Compared with the precise algorithm, the classical heuristic algorithm has short computation time and can get the initial feasible solution quickly. However, it is limited to small and medium-sized vehicle routing problems and is often used for local optimization of existing paths. Modern heuristic algorithm not only has stable performance, fast calculation speed, but also has flexible and powerful search ability, so it can solve large-scale complex vehicle routing problems. Modern heuristic algorithm becomes the main research direction for solving vehicle routing problem, such as genetic algorithm is used in many literature, due to its good solving performance, and is widely used by experts and scholars at home and abroad to solve vehicle path problems. However, with the deepening of research, the single modern heuristic algorithm cannot meet the needs of research gradually, and the hybrid algorithm arises at the historic moment. The hybrid algorithm using multiple search mechanisms improves the shortcoming of single modern heuristic algorithm, overcomes the limitation of single heuristic algorithm effectively, and improves the solving speed and quality of the algorithm.Compared with other algorithms, genetic algorithm indeed has many excellent properties, which provide a good way to solve vehicle routing problems. With the deepening of the research, some defects are also exposed, such as premature phenomenon. But genetic algorithm is parallel, so it can be well combined with other algorithms, such as scan algorithm, tabu search algorithm, simulated annealing, and immune algorithm, so as to make up for these defects and form a better hybrid algorithm.We design an improved adaptive genetic algorithm that in view of the above defects of genetic algorithm. Genetic algorithm can be improved in coding structure, crossover, and mutation operation and population. We use natural number coding, design adaptive crossover and mutation probabilities, design reversal operator to correct errors, and optimize the initial population with scanning algorithm, so as to reduce the premature defects of genetic algorithm and improve the performance of the algorithm. The specific design is discussed in the following sections of this chapter.


### 3.2.1. Chromosome Coding Design

Due to the special structure of the solution of the vehicle routing problem in cigarette distribution, a more intuitive natural number coding method is adopted in the chromosome coding design. Each chromosome is composed of the serial number of the depot and the customer point, and the customers and the depot are arranged together. If an existing distribution network composed ofn customers’ needs m vehicles for distribution, and the vehicles start from the depot and do not return to the depot after completing the task, the chromosome code can be expressed as follows:
(21)0,s11,s12⋯,s1p,0,s21,s22,⋯,s2q,0,⋯,0,sm1,sm1,⋯,smw.The length of chromosome isn+m, where skj is the jth customer served by the kth car, and m zero numbers represent the depot, which divides the chromosome into m segments, that is, there are m paths in total.


### 3.2.2. Population Initialization

Starting from a better initial population, genetic algorithm can find a better solution faster, so it is necessary to find a high-quality initial population as the initial solution after completing the chromosome coding. The scanning algorithm is used to generate the initial population. The steps are as follows:Step 1.
Taking the depot as the pole and the connecting line between the depot and any customer as the polar axis, the polar coordinate system is established and the polar coordinates of each customer point are calculatedStep 2.
Scan the customer points’ clockwise. Check the constraint condition once for each additional customer. When the constraint condition is violated, eliminate the customer and insert 0 at the beginning of the path. Then continue scanning until all customers are scanned. The customer sequence obtained after scanning is a chromosome. For multiple customers at the same angle, the order is from near to far, until the initial population with n chromosomes is generated


### 3.2.3. Fitness Function

In general, the higher the fitness value, the better the performance of the chromosome, and the greater the probability of entering the next generation, and vice versa. Therefore, the fitness function can be set as the reciprocal of the objective function. If the objective function value of chromosomei is zi, then the fitness value of chromosome i is fi=1/zi.


### 3.2.4. Genetic Operator

(1) Selection Operator. In order to complete the genetic evolution operation, we need to select a certain number of individuals from the initial population to form a new population. In this paper, we use the improved roulette selection method as the selection operator. The specific steps are as follows:Step 1.
The fitness values of individuals were ranked to ensure that the top m individuals with high fitness entered the next generationStep 2.
Roulette is used to select and copy the remaining individuals. The probability is calculated by formulaPi=c1−ci−1, where c is the selection probability of the ordered individuals, and the cumulative probability is calculated by formula Qi=∑i=1nPiStep 3.
Randomly generated numberr between [0, 1] is used as the selection pointer to determine the selected individual. If r≤Qi, individual 1 is selected; if Qi−1≤r≤Qi, individual i is selected;Step 4.
Repeat step 3 until (n-m) individuals are generatedIn order to ensure the genetic evolution operation and the diversity of the population, the crossover and mutation operation adopt the adaptive adjustment strategy. The functions of adaptive crossover probabilitypc and mutation probability pm are as follows:
(22)pc=pc1f−favgfmax−favg,f>favgpc2,f≤favgpm=pm1f′−favgfmax−favgf′>favgpm2f′≤favgAmong them,f is the larger fitness function value of the two crossed chromosomes; fmax is the maximum fitness function value of the population; favg is the average fitness function value of the population; f′ is the fitness function value of the chromosome to be mutated; Pc1 and Pc2 are in the [0.8, 1] interval; and pm1 and pm2 are in the [0.001, 0.1] interval.(2) Crossover Operator. Crossover operation is very important in genetic algorithm, which is the process of exchanging some genes of paired chromosomes in some way to form two new chromosomes. In order to avoid destroying excellent genes, the substring is moved to the first place when crossing, instead of directly copying in the past. This can not only protect excellent substrings but also avoid the generation of new individuals when the parents are the same, which improves the optimization ability and convergence speed of the algorithm. Crossover operation is the process of exchanging some genes of paired chromosomes in some way to form two new chromosomes, which is very important in genetic algorithms. In order to avoid destroying excellent genes, it is necessary to move the substrings to the first place during the crossover operation instead of copying the substrings directly over. This can protect the good substrings and avoids producing new individuals when the parents are identical, which improves the algorithm’s optimization ability and convergence speed. The specific operations are as follows:Step 1.
The adaptive crossover probabilitypc is calculated and the parent chromosomes A and B to be crossed are selected.Step 2.
In the interval1,k, integers r1 and r2 are generated randomly, and the sub path r1 in A is copied to the front end of the child B*, and the sub path r2 in B is copied to the front end of the child A*.Step 3.
The remaining k-1 zeros are randomly inserted into the remaining spaces, and the last space is not zero, and any two zeros cannot be adjacent.Step 4.
The remaining genes in the parent chromosomeA were inserted into the offspring A* to generate the complete offspring chromosome A*. Similarly, the progeny chromosome B* was generated.Step 5.
The decoding operation is performed to check whether the constraint conditions are met, and if not, discard.(3) Mutation Operator. In order to avoid losing some useful genes and prevent some chromosomes from being constant, mutation operation changes the position of some genes in the operation process, which can effectively prevent the algorithm from premature. In this adaptive mutation operation, the exchange mutation mode is adopted, and the specific operation is as follows:Step 1.
The adaptive mutation probabilitypm was calculated to select chromosomes.Step 2.
Two unequal integerst1 and t2 are randomly generated in the range of 1,n+k.Step 3.
Ift1 and t2 are both non-zero genes, the position will be exchanged, otherwise it will be generated again.Step 4.
Check whether the new offspring chromosome meets the constraints, if not, it will be generated again.Take chromosome 056104780392 as an example, ift1=2 and t2=6, and then, the chromosome generated by the above operator is 046105780392, as shown in Figure 1.(4) Evolution Reversal Operation. In order to correct the error in the process of evolution and improve the local search ability of the genetic algorithm, a gene segment is randomly selected from the parent generation and reversed in one direction. The specific operation is as follows:Step 1.
In interval [1,n+k], s1 and s2 are generated randomly.Step 2.
Ifs1 and s2 are both nonzero genes, the client genes between them will be reversed, otherwise they will be generated again.Step 3.
Check whether the offspring chromosomes meet the constraints, and whether the fitness value is improved. If not, discard them.Take chromosome 056104780392 as an example, ifs1=2, s2=6, and then using the above operator, the generated chromosome is 036105870429, as shown in Figure 2.


### 3.2.5. Termination Conditions

We set appropriate termination conditions to avoid the endless running of the algorithm or a dead loop without knowing. When an optimization result appears several times or reaches the maximum number of iterations, the operation is terminated and we output the satisfactory solution or approximate optimal solution. In this paper, the maximum number of iterations is 1000 as the termination condition.


## 3.2.1. Chromosome Coding Design

Due to the special structure of the solution of the vehicle routing problem in cigarette distribution, a more intuitive natural number coding method is adopted in the chromosome coding design. Each chromosome is composed of the serial number of the depot and the customer point, and the customers and the depot are arranged together. If an existing distribution network composed ofn customers’ needs m vehicles for distribution, and the vehicles start from the depot and do not return to the depot after completing the task, the chromosome code can be expressed as follows:
(21)0,s11,s12⋯,s1p,0,s21,s22,⋯,s2q,0,⋯,0,sm1,sm1,⋯,smw.The length of chromosome isn+m, where skj is the jth customer served by the kth car, and m zero numbers represent the depot, which divides the chromosome into m segments, that is, there are m paths in total.


## 3.2.2. Population Initialization

Starting from a better initial population, genetic algorithm can find a better solution faster, so it is necessary to find a high-quality initial population as the initial solution after completing the chromosome coding. The scanning algorithm is used to generate the initial population. The steps are as follows:Step 1.
Taking the depot as the pole and the connecting line between the depot and any customer as the polar axis, the polar coordinate system is established and the polar coordinates of each customer point are calculatedStep 2.
Scan the customer points’ clockwise. Check the constraint condition once for each additional customer. When the constraint condition is violated, eliminate the customer and insert 0 at the beginning of the path. Then continue scanning until all customers are scanned. The customer sequence obtained after scanning is a chromosome. For multiple customers at the same angle, the order is from near to far, until the initial population with n chromosomes is generated


## 3.2.3. Fitness Function

In general, the higher the fitness value, the better the performance of the chromosome, and the greater the probability of entering the next generation, and vice versa. Therefore, the fitness function can be set as the reciprocal of the objective function. If the objective function value of chromosomei is zi, then the fitness value of chromosome i is fi=1/zi.


## 3.2.4. Genetic Operator

(1) Selection Operator. In order to complete the genetic evolution operation, we need to select a certain number of individuals from the initial population to form a new population. In this paper, we use the improved roulette selection method as the selection operator. The specific steps are as follows:Step 1.
The fitness values of individuals were ranked to ensure that the top m individuals with high fitness entered the next generationStep 2.
Roulette is used to select and copy the remaining individuals. The probability is calculated by formulaPi=c1−ci−1, where c is the selection probability of the ordered individuals, and the cumulative probability is calculated by formula Qi=∑i=1nPiStep 3.
Randomly generated numberr between [0, 1] is used as the selection pointer to determine the selected individual. If r≤Qi, individual 1 is selected; if Qi−1≤r≤Qi, individual i is selected;Step 4.
Repeat step 3 until (n-m) individuals are generatedIn order to ensure the genetic evolution operation and the diversity of the population, the crossover and mutation operation adopt the adaptive adjustment strategy. The functions of adaptive crossover probabilitypc and mutation probability pm are as follows:
(22)pc=pc1f−favgfmax−favg,f>favgpc2,f≤favgpm=pm1f′−favgfmax−favgf′>favgpm2f′≤favgAmong them,f is the larger fitness function value of the two crossed chromosomes; fmax is the maximum fitness function value of the population; favg is the average fitness function value of the population; f′ is the fitness function value of the chromosome to be mutated; Pc1 and Pc2 are in the [0.8, 1] interval; and pm1 and pm2 are in the [0.001, 0.1] interval.(2) Crossover Operator. Crossover operation is very important in genetic algorithm, which is the process of exchanging some genes of paired chromosomes in some way to form two new chromosomes. In order to avoid destroying excellent genes, the substring is moved to the first place when crossing, instead of directly copying in the past. This can not only protect excellent substrings but also avoid the generation of new individuals when the parents are the same, which improves the optimization ability and convergence speed of the algorithm. Crossover operation is the process of exchanging some genes of paired chromosomes in some way to form two new chromosomes, which is very important in genetic algorithms. In order to avoid destroying excellent genes, it is necessary to move the substrings to the first place during the crossover operation instead of copying the substrings directly over. This can protect the good substrings and avoids producing new individuals when the parents are identical, which improves the algorithm’s optimization ability and convergence speed. The specific operations are as follows:Step 1.
The adaptive crossover probabilitypc is calculated and the parent chromosomes A and B to be crossed are selected.Step 2.
In the interval1,k, integers r1 and r2 are generated randomly, and the sub path r1 in A is copied to the front end of the child B*, and the sub path r2 in B is copied to the front end of the child A*.Step 3.
The remaining k-1 zeros are randomly inserted into the remaining spaces, and the last space is not zero, and any two zeros cannot be adjacent.Step 4.
The remaining genes in the parent chromosomeA were inserted into the offspring A* to generate the complete offspring chromosome A*. Similarly, the progeny chromosome B* was generated.Step 5.
The decoding operation is performed to check whether the constraint conditions are met, and if not, discard.(3) Mutation Operator. In order to avoid losing some useful genes and prevent some chromosomes from being constant, mutation operation changes the position of some genes in the operation process, which can effectively prevent the algorithm from premature. In this adaptive mutation operation, the exchange mutation mode is adopted, and the specific operation is as follows:Step 1.
The adaptive mutation probabilitypm was calculated to select chromosomes.Step 2.
Two unequal integerst1 and t2 are randomly generated in the range of 1,n+k.Step 3.
Ift1 and t2 are both non-zero genes, the position will be exchanged, otherwise it will be generated again.Step 4.
Check whether the new offspring chromosome meets the constraints, if not, it will be generated again.Take chromosome 056104780392 as an example, ift1=2 and t2=6, and then, the chromosome generated by the above operator is 046105780392, as shown in Figure 1.(4) Evolution Reversal Operation. In order to correct the error in the process of evolution and improve the local search ability of the genetic algorithm, a gene segment is randomly selected from the parent generation and reversed in one direction. The specific operation is as follows:Step 1.
In interval [1,n+k], s1 and s2 are generated randomly.Step 2.
Ifs1 and s2 are both nonzero genes, the client genes between them will be reversed, otherwise they will be generated again.Step 3.
Check whether the offspring chromosomes meet the constraints, and whether the fitness value is improved. If not, discard them.Take chromosome 056104780392 as an example, ifs1=2, s2=6, and then using the above operator, the generated chromosome is 036105870429, as shown in Figure 2.


## 3.2.5. Termination Conditions

We set appropriate termination conditions to avoid the endless running of the algorithm or a dead loop without knowing. When an optimization result appears several times or reaches the maximum number of iterations, the operation is terminated and we output the satisfactory solution or approximate optimal solution. In this paper, the maximum number of iterations is 1000 as the termination condition.


## 4. Numerical Example Analysis

### 4.1. Case Introduction

#### 4.1.1. Introduction of Chongqing Tobacco Logistics Background

(1) Introduction to Chongqing Tobacco Logistics Background. Chongqing Tobacco Company owned branches in all districts and counties of Chongqing, and its chain stores basically covered the county-level areas in 2005. Chongqing Tobacco Logistics Company is determined to spread Chongqing Tobacco all over Chongqing. Nowadays, Chongqing tobacco logistics has distribution centers and subwarehouses in middle, northern, and western Chongqing. It has realized rapid distribution of cigarette logistics in the whole Chongqing.After the further integration of cigarette logistics resources in 2015, Chongqing tobacco logistics now has four distribution centers, Jiangbei, Fuling, Qianjiang, and Wanzhou, its radiation range is shown in Table2. Chongqing tobacco logistics adopts the distribution mode of combining direct distribution and transit distribution, which makes the management more simpler and flexible. It not only ensures the distribution of cigarettes on time and in quantity but also effectively reduces the distribution cost and improves the distribution speed.Table 2
Radiation scope of Chongqing Tobacco Logistics’ distribution centers.
Distribution centerRadiation scopeJiangbeiJiangbei District, Nanan District, Yuzhong District, Yubei District, Beibei District, Hechuan District, Dadukou District, Shapingba District, Jiulongpo District, Yongchuan District, Banan District, Tongliang County, Bishan County, Qijiang County, Dazu County, Rongchang County, Tongnan CountyFulingWansheng District, Changshou District, Nanchuan District, Fuling District, Fengdu County, Dianjiang County, Wulong County, Shizhu CountyQianjiangQianjiang District, Pengshui County Youyang County, Xiushan CountyWanzhouWanzhou District, Yunyang County, Fengjie County, Kaixian County, Wuxi County, Liangping County, Wushan County, Zhongxian County, Chengkou County


#### 4.1.2. Case Data Information

This case uses 50 customer nodes of distribution, 5 distribution centers, and 1 transit as an example. The same type of IVECO transport vehicle is used to carry out the distribution task. The dimensions of IVECO box are 5.9 meters long, 2.7 meters wide, and 2 meters high, and the carrying capacity is 10 tons. It can carry about 30 boxes of standard box volume of tobacco.The location coordinates of customer nodes are subject to longitude and latitude. The geographical coordinates of distribution center and transfer station are as follows: (a) distribution center (106.71, 26.50), (b) transfer field (106.46, 25.32), (c) transfer field (106.72, 25.34), (d) transfer field (106.82, 25.28), (e) transfer field (106.79, 25.21), and (f) transfer field (106.27, 25.29). The specific customer information data are shown in Table3.Table 3
Geographical location and demand of customer nodes.
Branch numberPosition coordinatesDelivery volume/caseBranch numberPosition coordinatesDelivery volume/case1(105.21,25.83)426(106.62,27.80)32(105.49,25.11)327(105.79,25.00)33(106.09,25.17)328(105.64,25.81)34(105.82,26.58)529(105.76,26.66)45(106.59,26.84)430(105.92,26.25)56(106.66,26.14)331(106.46,25.16)37(106.55,25.83)432(106.51,25.32)38(106.54,25.84)433(106.67,25.53)69(106.97,26.59)534(106.66,25.43)210(106.88,27.22)435(106.72,25.33)511(105.61,27.16)436(106.64,25.03)512(105.38,26.77)437(106.75,25.05)513(106.22,27.46)238(106.80,25.24)314(106.41,27.81)439(106.88,25.31)515(106.71,27.13)440(106.95,25.16)316(106.58,26.49)541(106.83,25.11)417(106.50,25.94)542(106.81,25.08)218(106.48,25.92)343(106.53,25.05)419(106.68,26.98)644(106.36,25.12)420(106.82,27.71)545(106.33,25.34)421(106.72,27.21)446(106.28,25.03)322(106.85,27.73)347(106.25,25.01)423(107.08,27.71)448(106.12,25.37)324(106.50,27.76)349(106.55,25.39)625(106.24,27.52)550(106.70,25.56)4(Note: The data in this table are collected from Baidu Map, and the driving distance in this paper adopts the mathematical coordinate distance between two points).


### 4.2. Case Analysis

MATLAB is used to write the genetic algorithm program. In the parameters part, the unit distance cost ofc1 vehicle is c1=2 yuan/km, the unit fuel consumption cost is c2= 5.6 yuan/L, and the fixed loss cost of each vehicle is c3= 1000 yuan/vehicle. Query the international fuel consumption statistics literature, this paper adopts the fuel consumption coefficientω=6.20×10−3. Firstly, the clustering algorithm is used to divide the customer nodes into regions, and then, the genetic algorithm is used to solve the problem. In this paper, the population size of the algorithm is set as 200, the number of iterations is 300, the probability of crossover and mutation is Pc=0.9 and Pm=0.1, respectively, and the optimal scheme is selected as the final scheme after running the algorithm for ten times. The distribution optimization path diagram is shown in Figure 4.Figure 4
Distribution route optimization diagram.As can be seen from Figure4, according to the clustering algorithm, the whole distribution area is divided into six areas, and two distribution methods, direct distribution and transfer distribution, are adopted. Distribution center a will directly distribute the goods in the distribution area after sorting, while distribution center a will sort and package the goods in b,c,d,e, and f and then transport them to the transfer station in the area, and the transfer station will organize the next distribution. Customer nodes 4 and 25 belong to the areas covered by transit yard c and transit yard e, respectively. However, due to the fact that customer node 4 is close to the area covered by transit yard b in the process of distribution, the principle of nearby distribution is adopted. Customer node 4 is delivered by transit vehicles in transit yard b, and customer node 25 is delivered by transit vehicles in transit yard f in the same way, thus realizing cross regional joint distribution and saving distribution costs. The distribution route results are shown in Table 4.Table 4
Solution results of cross regional joint distribution routes.
Vehicle numberDistribution routeLoading rateDistance/kmCost/yuan1a-37-10-19-9-17-16-a100%54.621166.132a-5-29-30-6-8-7-36-a96.7%73.561221.193a-b-32-4-28-1-2-27-3-31-b90%185.711545.484a-c-12-11-35-50-15-34-33-c96.7%135.821408.335a-d-22-38-39-13-14-23-40-24-d90%193.401568.106a-e-49-41-42-26-43-20-21-e93.3%145.311431.857a-f-45-25-44-46-47-49-48-18-f90%196.621577.50Total985.049918.59It can be seen from Table4 that there are 7 vehicles for this distribution, the total driving distance is 985.04 km, and the total cost is 9918.59 yuan.In order to prove the effectiveness of cross-regional joint distribution route optimization, this paper makes a comparative analysis with the independent distribution scheme of each region without cross-regional joint distribution. The distribution vehicles are all distributed in their respective districts and shall not exceed their respective districts. The distribution routes without cross-regional joint distribution are shown in Figure5.Figure 5
Independent distribution route map.Genetic algorithm was adopted to solve the scheme of independent distribution in each region by setting the same parameters as the above experiment. The results are shown in Table5.Table 5
Solution results of independent distribution scheme.
Vehicle numberDistribution routeLoading rateDistance/kmCost/yuan1a-37-10-19-9-17-16-a100%54.621166.132a-5-29-30-6-8-7-36-a96.7%73.561221.193a-b-32-28-1-2-27-3-31-b73.3%168.621466.044a-c-4-12-11-35-c60%98.731259.165c-50-15-34-33-c53.3%76.191194.716a-d-22-38-39-13-14-23-40-24-d90%193.401568.107a-e-49-41-26-42-e50%91.851231.548e-21-25-43-20-e63.3%68.671187.419a-f-45-44-46-47-49-48-18-f70%183.521500.85Total1009.1611795.12It can be seen from Table6 that there are 7 vehicles in cross-regional joint distribution and 9 vehicles in independent regional distribution. The average loading rate of vehicles in cross-regional joint distribution is 93.4%, and that in independent regional distribution is 73.1%. The utilization rate of vehicles in cross-regional joint distribution is better than that in independent regional distribution, because cross-regional joint distribution can realize resource sharing and vehicle sharing between regions and then make full use of resources and vehicles. The total driving distance of cross-regional joint distribution is 985.04 km, and the total cost is 9918.59 yuan, while the total driving distance of independent regional distribution is 1009.16 km, and the total cost is 11795.12 yuan. The total driving distance and assembly cost of cross-regional joint distribution are better than those of independent regional distribution. This is because the cross-regional joint distribution adopts the principle of the nearest distribution, which realizes resource sharing and information sharing, so that the vehicle can maximize the utilization rate of the vehicle, avoiding the situation of redeparture due to the full load during the delivery process, thereby reducing the cost of distribution, and proving the superiority of cross-regional joint distribution.Table 6
Comparison of solution results of two schemes.
CategoryNumber of vehiclesAverage loading rateTotal distribution distance/kmTotal cost/yuanCross-regional joint distribution793.4%985.049918.59Independent regional distribution973.1%1009.1611795.12It can be seen from Figure6 that the total cost obtained by the genetic algorithm is about 9600 yuan, and the iteration tends to converge in about 110 times, which shows that the algorithm used in this paper is effective and feasible, and the optimization result is ideal.Figure 6
Iterative graph of adaptive genetic algorithm.


## 4.1. Case Introduction

### 4.1.1. Introduction of Chongqing Tobacco Logistics Background

(1) Introduction to Chongqing Tobacco Logistics Background. Chongqing Tobacco Company owned branches in all districts and counties of Chongqing, and its chain stores basically covered the county-level areas in 2005. Chongqing Tobacco Logistics Company is determined to spread Chongqing Tobacco all over Chongqing. Nowadays, Chongqing tobacco logistics has distribution centers and subwarehouses in middle, northern, and western Chongqing. It has realized rapid distribution of cigarette logistics in the whole Chongqing.After the further integration of cigarette logistics resources in 2015, Chongqing tobacco logistics now has four distribution centers, Jiangbei, Fuling, Qianjiang, and Wanzhou, its radiation range is shown in Table2. Chongqing tobacco logistics adopts the distribution mode of combining direct distribution and transit distribution, which makes the management more simpler and flexible. It not only ensures the distribution of cigarettes on time and in quantity but also effectively reduces the distribution cost and improves the distribution speed.Table 2
Radiation scope of Chongqing Tobacco Logistics’ distribution centers.
Distribution centerRadiation scopeJiangbeiJiangbei District, Nanan District, Yuzhong District, Yubei District, Beibei District, Hechuan District, Dadukou District, Shapingba District, Jiulongpo District, Yongchuan District, Banan District, Tongliang County, Bishan County, Qijiang County, Dazu County, Rongchang County, Tongnan CountyFulingWansheng District, Changshou District, Nanchuan District, Fuling District, Fengdu County, Dianjiang County, Wulong County, Shizhu CountyQianjiangQianjiang District, Pengshui County Youyang County, Xiushan CountyWanzhouWanzhou District, Yunyang County, Fengjie County, Kaixian County, Wuxi County, Liangping County, Wushan County, Zhongxian County, Chengkou County


### 4.1.2. Case Data Information

This case uses 50 customer nodes of distribution, 5 distribution centers, and 1 transit as an example. The same type of IVECO transport vehicle is used to carry out the distribution task. The dimensions of IVECO box are 5.9 meters long, 2.7 meters wide, and 2 meters high, and the carrying capacity is 10 tons. It can carry about 30 boxes of standard box volume of tobacco.The location coordinates of customer nodes are subject to longitude and latitude. The geographical coordinates of distribution center and transfer station are as follows: (a) distribution center (106.71, 26.50), (b) transfer field (106.46, 25.32), (c) transfer field (106.72, 25.34), (d) transfer field (106.82, 25.28), (e) transfer field (106.79, 25.21), and (f) transfer field (106.27, 25.29). The specific customer information data are shown in Table3.Table 3
Geographical location and demand of customer nodes.
Branch numberPosition coordinatesDelivery volume/caseBranch numberPosition coordinatesDelivery volume/case1(105.21,25.83)426(106.62,27.80)32(105.49,25.11)327(105.79,25.00)33(106.09,25.17)328(105.64,25.81)34(105.82,26.58)529(105.76,26.66)45(106.59,26.84)430(105.92,26.25)56(106.66,26.14)331(106.46,25.16)37(106.55,25.83)432(106.51,25.32)38(106.54,25.84)433(106.67,25.53)69(106.97,26.59)534(106.66,25.43)210(106.88,27.22)435(106.72,25.33)511(105.61,27.16)436(106.64,25.03)512(105.38,26.77)437(106.75,25.05)513(106.22,27.46)238(106.80,25.24)314(106.41,27.81)439(106.88,25.31)515(106.71,27.13)440(106.95,25.16)316(106.58,26.49)541(106.83,25.11)417(106.50,25.94)542(106.81,25.08)218(106.48,25.92)343(106.53,25.05)419(106.68,26.98)644(106.36,25.12)420(106.82,27.71)545(106.33,25.34)421(106.72,27.21)446(106.28,25.03)322(106.85,27.73)347(106.25,25.01)423(107.08,27.71)448(106.12,25.37)324(106.50,27.76)349(106.55,25.39)625(106.24,27.52)550(106.70,25.56)4(Note: The data in this table are collected from Baidu Map, and the driving distance in this paper adopts the mathematical coordinate distance between two points).


## 4.1.1. Introduction of Chongqing Tobacco Logistics Background

(1) Introduction to Chongqing Tobacco Logistics Background. Chongqing Tobacco Company owned branches in all districts and counties of Chongqing, and its chain stores basically covered the county-level areas in 2005. Chongqing Tobacco Logistics Company is determined to spread Chongqing Tobacco all over Chongqing. Nowadays, Chongqing tobacco logistics has distribution centers and subwarehouses in middle, northern, and western Chongqing. It has realized rapid distribution of cigarette logistics in the whole Chongqing.After the further integration of cigarette logistics resources in 2015, Chongqing tobacco logistics now has four distribution centers, Jiangbei, Fuling, Qianjiang, and Wanzhou, its radiation range is shown in Table2. Chongqing tobacco logistics adopts the distribution mode of combining direct distribution and transit distribution, which makes the management more simpler and flexible. It not only ensures the distribution of cigarettes on time and in quantity but also effectively reduces the distribution cost and improves the distribution speed.Table 2
Radiation scope of Chongqing Tobacco Logistics’ distribution centers.
Distribution centerRadiation scopeJiangbeiJiangbei District, Nanan District, Yuzhong District, Yubei District, Beibei District, Hechuan District, Dadukou District, Shapingba District, Jiulongpo District, Yongchuan District, Banan District, Tongliang County, Bishan County, Qijiang County, Dazu County, Rongchang County, Tongnan CountyFulingWansheng District, Changshou District, Nanchuan District, Fuling District, Fengdu County, Dianjiang County, Wulong County, Shizhu CountyQianjiangQianjiang District, Pengshui County Youyang County, Xiushan CountyWanzhouWanzhou District, Yunyang County, Fengjie County, Kaixian County, Wuxi County, Liangping County, Wushan County, Zhongxian County, Chengkou County


## 4.1.2. Case Data Information

This case uses 50 customer nodes of distribution, 5 distribution centers, and 1 transit as an example. The same type of IVECO transport vehicle is used to carry out the distribution task. The dimensions of IVECO box are 5.9 meters long, 2.7 meters wide, and 2 meters high, and the carrying capacity is 10 tons. It can carry about 30 boxes of standard box volume of tobacco.The location coordinates of customer nodes are subject to longitude and latitude. The geographical coordinates of distribution center and transfer station are as follows: (a) distribution center (106.71, 26.50), (b) transfer field (106.46, 25.32), (c) transfer field (106.72, 25.34), (d) transfer field (106.82, 25.28), (e) transfer field (106.79, 25.21), and (f) transfer field (106.27, 25.29). The specific customer information data are shown in Table3.Table 3
Geographical location and demand of customer nodes.
Branch numberPosition coordinatesDelivery volume/caseBranch numberPosition coordinatesDelivery volume/case1(105.21,25.83)426(106.62,27.80)32(105.49,25.11)327(105.79,25.00)33(106.09,25.17)328(105.64,25.81)34(105.82,26.58)529(105.76,26.66)45(106.59,26.84)430(105.92,26.25)56(106.66,26.14)331(106.46,25.16)37(106.55,25.83)432(106.51,25.32)38(106.54,25.84)433(106.67,25.53)69(106.97,26.59)534(106.66,25.43)210(106.88,27.22)435(106.72,25.33)511(105.61,27.16)436(106.64,25.03)512(105.38,26.77)437(106.75,25.05)513(106.22,27.46)238(106.80,25.24)314(106.41,27.81)439(106.88,25.31)515(106.71,27.13)440(106.95,25.16)316(106.58,26.49)541(106.83,25.11)417(106.50,25.94)542(106.81,25.08)218(106.48,25.92)343(106.53,25.05)419(106.68,26.98)644(106.36,25.12)420(106.82,27.71)545(106.33,25.34)421(106.72,27.21)446(106.28,25.03)322(106.85,27.73)347(106.25,25.01)423(107.08,27.71)448(106.12,25.37)324(106.50,27.76)349(106.55,25.39)625(106.24,27.52)550(106.70,25.56)4(Note: The data in this table are collected from Baidu Map, and the driving distance in this paper adopts the mathematical coordinate distance between two points).


## 4.2. Case Analysis

MATLAB is used to write the genetic algorithm program. In the parameters part, the unit distance cost ofc1 vehicle is c1=2 yuan/km, the unit fuel consumption cost is c2= 5.6 yuan/L, and the fixed loss cost of each vehicle is c3= 1000 yuan/vehicle. Query the international fuel consumption statistics literature, this paper adopts the fuel consumption coefficientω=6.20×10−3. Firstly, the clustering algorithm is used to divide the customer nodes into regions, and then, the genetic algorithm is used to solve the problem. In this paper, the population size of the algorithm is set as 200, the number of iterations is 300, the probability of crossover and mutation is Pc=0.9 and Pm=0.1, respectively, and the optimal scheme is selected as the final scheme after running the algorithm for ten times. The distribution optimization path diagram is shown in Figure 4.Figure 4
Distribution route optimization diagram.As can be seen from Figure4, according to the clustering algorithm, the whole distribution area is divided into six areas, and two distribution methods, direct distribution and transfer distribution, are adopted. Distribution center a will directly distribute the goods in the distribution area after sorting, while distribution center a will sort and package the goods in b,c,d,e, and f and then transport them to the transfer station in the area, and the transfer station will organize the next distribution. Customer nodes 4 and 25 belong to the areas covered by transit yard c and transit yard e, respectively. However, due to the fact that customer node 4 is close to the area covered by transit yard b in the process of distribution, the principle of nearby distribution is adopted. Customer node 4 is delivered by transit vehicles in transit yard b, and customer node 25 is delivered by transit vehicles in transit yard f in the same way, thus realizing cross regional joint distribution and saving distribution costs. The distribution route results are shown in Table 4.Table 4
Solution results of cross regional joint distribution routes.
Vehicle numberDistribution routeLoading rateDistance/kmCost/yuan1a-37-10-19-9-17-16-a100%54.621166.132a-5-29-30-6-8-7-36-a96.7%73.561221.193a-b-32-4-28-1-2-27-3-31-b90%185.711545.484a-c-12-11-35-50-15-34-33-c96.7%135.821408.335a-d-22-38-39-13-14-23-40-24-d90%193.401568.106a-e-49-41-42-26-43-20-21-e93.3%145.311431.857a-f-45-25-44-46-47-49-48-18-f90%196.621577.50Total985.049918.59It can be seen from Table4 that there are 7 vehicles for this distribution, the total driving distance is 985.04 km, and the total cost is 9918.59 yuan.In order to prove the effectiveness of cross-regional joint distribution route optimization, this paper makes a comparative analysis with the independent distribution scheme of each region without cross-regional joint distribution. The distribution vehicles are all distributed in their respective districts and shall not exceed their respective districts. The distribution routes without cross-regional joint distribution are shown in Figure5.Figure 5
Independent distribution route map.Genetic algorithm was adopted to solve the scheme of independent distribution in each region by setting the same parameters as the above experiment. The results are shown in Table5.Table 5
Solution results of independent distribution scheme.
Vehicle numberDistribution routeLoading rateDistance/kmCost/yuan1a-37-10-19-9-17-16-a100%54.621166.132a-5-29-30-6-8-7-36-a96.7%73.561221.193a-b-32-28-1-2-27-3-31-b73.3%168.621466.044a-c-4-12-11-35-c60%98.731259.165c-50-15-34-33-c53.3%76.191194.716a-d-22-38-39-13-14-23-40-24-d90%193.401568.107a-e-49-41-26-42-e50%91.851231.548e-21-25-43-20-e63.3%68.671187.419a-f-45-44-46-47-49-48-18-f70%183.521500.85Total1009.1611795.12It can be seen from Table6 that there are 7 vehicles in cross-regional joint distribution and 9 vehicles in independent regional distribution. The average loading rate of vehicles in cross-regional joint distribution is 93.4%, and that in independent regional distribution is 73.1%. The utilization rate of vehicles in cross-regional joint distribution is better than that in independent regional distribution, because cross-regional joint distribution can realize resource sharing and vehicle sharing between regions and then make full use of resources and vehicles. The total driving distance of cross-regional joint distribution is 985.04 km, and the total cost is 9918.59 yuan, while the total driving distance of independent regional distribution is 1009.16 km, and the total cost is 11795.12 yuan. The total driving distance and assembly cost of cross-regional joint distribution are better than those of independent regional distribution. This is because the cross-regional joint distribution adopts the principle of the nearest distribution, which realizes resource sharing and information sharing, so that the vehicle can maximize the utilization rate of the vehicle, avoiding the situation of redeparture due to the full load during the delivery process, thereby reducing the cost of distribution, and proving the superiority of cross-regional joint distribution.Table 6
Comparison of solution results of two schemes.
CategoryNumber of vehiclesAverage loading rateTotal distribution distance/kmTotal cost/yuanCross-regional joint distribution793.4%985.049918.59Independent regional distribution973.1%1009.1611795.12It can be seen from Figure6 that the total cost obtained by the genetic algorithm is about 9600 yuan, and the iteration tends to converge in about 110 times, which shows that the algorithm used in this paper is effective and feasible, and the optimization result is ideal.Figure 6
Iterative graph of adaptive genetic algorithm.


## 5. Conclusion

The developments presented in this paper was aimed at providing a theoretical framework for formulating and solving the cigarette distribution route optimization problem to reducing logistics cost of cigarettes in China. The conclusions are as follows:(1)
We design an improved adaptive genetic algorithm that in view of the above defects of genetic algorithm the coding structure, crossover and mutation operation, and population were improved. We use natural number coding, design adaptive crossover and mutation probabilities, design reversal operator to correct errors, and optimize the initial population with scanning algorithm, so as to reduce the premature defects of genetic algorithm and improve the performance of the algorithmThis paper makes a comparative analysis of the cross-regional joint distribution of logistics distribution scheme and the independent distribution scheme of each region without cross-regional joint distribution. There are 7 vehicles in cross-regional joint distribution and 9 vehicles in independent regional distribution. The average loading rate of vehicles in cross-regional joint distribution is 93.4%, and that in independent regional distribution is 73.1%. The utilization rate of vehicles in cross-regional joint distribution is better than that in independent regional distribution.The innovative points of this paper are as follows:(1)
In this paper, a two-stage method is used to establish the model according to the characteristics of cigarette distribution(a) The first stage: distribution areas are dividedWe determine the distribution center and distribution areas, transfer station (clustering center), and covered retail merchants according to the given distribution radiation radius. Clustering algorithm is used to divide or group customer nodes. This makes the retail customers with relatively concentrated geographical distribution in the distribution area of the same car(b) The second stage: distribution routes are optimizedThe genetic algorithm is used to solve the TSP problem and the specific cigarette delivery routes or sequences of delivery vehicles within a single delivery area(2)
We design an improved adaptive genetic algorithm aimed at the above defects of genetic algorithm. There are some improvements of the coding structure, crossover, and mutation operation and populationWe hope that these areas will be improved in the future with possible extensions(a)
In terms of model: In order to realize the mathematical modelling, we make a series of assumptions about the VRP problem to simplify the real problem, so there may be biases between the resulting data and the actual data in this paper, which need to be improved in the future research(b)
In terms of algorithm: Due to the poor understanding of the hybrid heuristic algorithm, we only used the scanning algorithm to optimize the initial population and improve some operations such as crossover and variation but did not use other heuristics to improve the deeper content. This study still have many limitations but also have some issues for further in-depth study



---

*Source: 3830137-2022-04-05.xml*