# Green Logistics Location-Routing Optimization Solution Based on Improved GA A1gorithm considering Low-Carbon and Environmental Protection

Authors: Anqing Zhu; Youyun Wen
Journal: Journal of Mathematics
 (2021)

---

## Abstract

This paper proposes a green logistics location-routing optimization problem based on improved genetic algorithm (GA) from the perspective of low-carbon and environmental protection. First, considering the cost factor, time window, deterioration rate of agricultural products, inventory and distribution capacity, carbon trading mechanism, and other factors, and with the total cost minimization as the optimization goal, a low-carbon and environmental protection logistics location-routing optimization model is constructed. Then, the adaptive operator and cataclysm operator are introduced to improve the GA algorithm, which can adjust crossover and mutation probability according to the needs, reducing the influence of parameters and running time. Furthermore, the improved GA algorithm is used to solve the location-routing optimization problem in green logistics, so as to obtain a low-carbon, economical, and efficient distribution path. Finally, perform experimental analysis of the proposed method using the relevant data of U company. The results show that the total distribution cost is 6771.3 yuan, which meets the design requirements of economy and environmental protection.


---

## Body

## 1. Introduction

With the development of modern logistics to informatization, globalization, and integration, distribution plays an important role in the entire logistics system [1]. In particular, the development of cold chain logistics not only needs to meet safety and freshness but also needs to meet relevant policies, laws, and regulations of cold chain logistics in terms of energy consumption and carbon emissions [2, 3]. Since the energy consumption and carbon emissions of cold chain logistics are usually much higher than those of traditional room temperature logistics, the location of distribution points and optimization of transportation routes show a strong correlation with energy consumption and emissions. Therefore, the location for cold chain green logistics-distribution route optimization has become a key research issue.Taking green logistics as the research object, consider the location-routing optimization problem for cold chain distribution of fresh agricultural products based on carbon trading. This not only has a certain theoretical significance but also has great practical significance [4]. From a theoretical point of view, the research on cold chain distribution of fresh agricultural products at home and abroad is mostly about the location of distribution centers or optimization of vehicle routing. The location problem and path optimization problem are actually related to each other and restrict each other, and joint optimization is of great significance [5]. At the same time, it is necessary to consider factors such as distance, vehicle speed, load capacity, and the carbon emissions of delivery vehicle itself and carbon emissions of refrigeration devices. From a practical point of view, the joint optimization of location and vehicle routing is more conducive to the effective use of resources by enterprises and lower storage and transportation costs. This can improve the level of customer service, and at the same time, the carbon trading mechanism may bring new sources of profit for enterprises. Through low-carbon and environmentally friendly distribution, it can reflect the company’s sustainable development concept and bring new competition points for companies [6]. Secondly, the entire cold chain in warehousing and distribution process can ensure food safety. Thus, the green logistics distribution location-routing optimization problem from the perspective of low-carbon and environmental protection is a comprehensive and complex comprehensive problem [7].For traditional logistics research, whether it is basic location or route planning, most of them only take the minimization of operating costs as optimization goal, and they rarely consider time windows and carbon emissions [8]. As the concept of low-carbon and environmental protection gradually gains popularity, the carbon emission problem in cold chain distribution has become a key issue for logistics companies to solve [9]. In addition, part of research in the cold chain distribution vehicle routing optimization problem that considers carbon emissions is mostly considered in the form of carbon tax, and carbon trading mechanisms are rarely considered [10, 11].However, with the improvement of environmental protection policies and carbon emission standards, the location of cold chain logistics is more than just optimization in terms of cost and timeliness. It needs to consider more factors in terms of carbon emissions. To this end, a green logistics location-routing optimization model from the perspective of low-carbon and environmental protection has been constructed. And, it used the improved GA algorithm to solve in order to get a reasonable and efficient logistics address and delivery route. Combined with carbon trading mechanism, it formulates the optimal distribution path based on full consideration of carbon emissions and time costs. The contribution of the proposed method is as follows:(1)
In view of the fact that most of existing logistics distribution path planning does not consider the low-carbon and environmental protection problems, a location-routing optimization problem of green logistics using improved the genetic algorithm (GA) from the perspective of low-carbon and environmental protection is proposed.(2)
The standard GA algorithm is susceptible to the influence of parameters and falls into local optimum. To this end, the adaptive operator and cataclysm operator are introduced. The improved catastrophic adaptive GA can adjust the probability of crossover and mutation in GA as needed, which improves optimization performance and reduces the running time.(3)
In order to realize the low-carbon and environmental protection of distribution route, the proposed method constructs a green logistics location-routing optimization problem model. It considers cost factors, time window, deterioration rate of agricultural products, inventory capacity and distribution capacity, carbon trading mechanism, and other influencing factors and takes the minimization of total cost as optimization goal. Besides, the improved GA algorithm is used to solve the problem, which improves the economy and environmental protection of logistics.


## 2. Related Work

The location-routing optimization problem of green logistics aims to improve service response speed, save costs, and improve customer satisfaction by reasonable spatial layout and planning [12]. Scholars at home and abroad can review research on the location-routing problem from the following three parts, namely, location-routing problem (LRP), vehicle-routing problem (VRP), and location-routing problem [13].First of all, there have been many scientific research results for LRP research. Reference [14] proposed a low-carbon cold chain logistics path optimization algorithm based on ribonucleic acid-ant colony. Based on the established cold chain logistics path optimization model including carbon emission cost, it used ribonucleic acid calculation and ant colony optimization to prevent the impact of unreasonable parameter selection on algorithm performance. Reference [15] proposed a rough set-based method of scoring advantages and disadvantages for the location of third-party cold chain drug logistics. They assigned different weights to different indicators and provided a new solution to the selection of third-party logistics providers. However, the above two methods only optimize the supplier location problem, and the overall location-routing optimization problem still needs to be studied in depth.Secondly, in terms of VRP, only the single route optimization problem of distribution cost was initially considered. Then, consider the path optimization problem of cost and service quality. Finally, the multiobjective path optimization problems such as cost, service quality, and environmental impact are considered, and the influence of multiple factors such as time window, real-time road conditions, and deterioration rate on the optimization model was considered [16]. In terms of model solving, VRP mainly has two types of solving methods: precise algorithm and heuristic algorithm. Accurate algorithms include branch and bound method, cutting plane method, and dynamic programming. Heuristic algorithms mainly include GA, ant colony algorithm, tabu search, and simulated annealing algorithm [17]. For example, reference [18] proposed a GA-based cold chain logistics network optimization model, with the total operating cost of logistics network as the goal. It emended a nonlinear mixed integer programming involving two allocation problems, which effectively determines the optimal layout of logistics centers and transportation routes. But it did not involve the issue of carbon emissions. Reference [19] proposed a single-objective optimization scheme based on the particle swarm optimization algorithm for the VRP problem. It effectively minimized the cost of logistics and distribution, but did not involve timeliness and carbon emission factors.Finally, scholars at home and abroad have done some research on the location-routing optimization of green logistics. However, there are few studies on such comprehensive multiobjectives. Reference [20] proposed an optimization scheme based on improved ant colony algorithm for multiobjective solution of cost, carbon emission, and customer satisfaction in the logistics distribution process. On the basis of ensuring customer satisfaction, this solution effectively reduced distribution costs and carbon emissions.


## 3. Model Construction and Solution

### 3.1. Model Construction of Green Logistics Location-Routing Optimization Problem considering Carbon Trading

#### 3.1.1. Model Assumption

The following assumptions are made for green logistics location-routing optimization problem considering carbon trading:(1)
Consider the capacity limitation of distribution centers(2)
During offloading time, the deterioration rate will be higher than deterioration rate during transportation.(3)
A demand point can only be delivered by one vehicle, and other vehicles are not allowed to deliver, that is, split delivery is not allowed.(4)
The vehicle runs at a constant speed, and temperature is constant during storage and transportation. After completing the distribution task, it will return to distribution centers empty. In this process, refrigeration is not required.(5)
It is assumed that, as the quality of fresh agricultural products decreases, the price will also decrease correspondingly, resulting in a certain loss of goods cost.(6)
For a given amount of carbon emission allowance, if the cold chain enterprise’s emissions do not exceed allowance, the part less than allowance can be converted into revenue at carbon trading price to offset part of the cost. If the cold chain company’s emissions exceed the quota, it needs to purchase additional quotas from carbon trading market, incurring additional carbon trading costs.


#### 3.1.2. Function Analysis

(1) Metamorphosis Rate Function. In the process of storage and distribution, the quality of fresh agricultural products will change with time and temperature. People’s requirements for food quality are getting higher and higher. When the quality of fresh agricultural products declines, sales price will also decline accordingly. Assuming that the temperature during storage and distribution is constant and strictly controlled, the spoilage rate of fresh agricultural products at this time is only related to time. Due to short loading and offloading time, the temperature change when opening the door is negligible. Therefore, at this constant temperature, the deterioration rate can be regarded as a constant. The quality of fresh produce can be regarded as exponentially changing over time.Assume thatQ0 is the quality of fresh agricultural products at time 0, that is, the quality of fresh agricultural products when the refrigerated vehicle departs from selected distribution centers. Use θ to represent the spoilage rate of fresh agricultural products. At a constant temperature, the spoilage rate of a fresh agricultural product is only related to time. If Qt is used to represent the quality of fresh agricultural products at time t, the mathematical expression is(1)Qt=Q0e−θt.(2) Carbon Trading Function. Considering the location-routing optimization problem of cold chain distribution of fresh agricultural products under carbon trading mechanism, the total amount of carbon emissions is obtained. In order to study the impact of carbon trading mechanism on the location-routing problem of cold chain distribution of fresh agricultural products, the total carbon emissions need to be monetized and converted into carbon trading costs. The conversion of total carbon emissions to carbon trading costs is achieved through carbon trading functions. The size of carbon trading costs depends on the total carbon emissions, carbon trading costs, and initial carbon allowances [21–25].When the company’s actual total carbon emissions are less than initial quota, the part that is less than the quota can be sold to carbon trading agency. This can yield benefits and offset some of the costs. When the company’s actual total carbon emissions are higher than initial allowance, it will need to spend a certain cost to purchase the carbon emissions allowance that is higher than initial allowance. Thus, the carbon trading function can be expressed as(2)C1=pc×ψC−ψ0=pc×ψCT,where pc represents the carbon transaction price, that is, the transaction price of unit carbon emission allowance; ψC represents the actual total carbon emission of enterprise; ψ0 represents the initial carbon emission allowance of enterprise; and ψCT is used to denote the company’s carbon trading volume, which is the difference between company’s actual total carbon emissions minus initial carbon emission allowances. The value of ψCT may be positive or negative. When ψCT≥0, it means that the company’s actual total carbon emissions are higher than initial carbon emissions quota. At this time, ψCT represents the carbon trading quota that the company has purchased from carbon trading market. When ψCT<0, it means that the company’s actual total carbon emission is lower than initial carbon emission quota. At this time, ψCT represents the carbon trading quota sold by the company to carbon trading market.(3) Objective Function. For the objective function of location-routing problem, the most common ones are cost minimization, time minimization, and profit maximization. For different situations, the selected objective function is also quite different. For the location-routing problem of distribution center in green logistics, considering the weak cost of green logistics, the time is generally chosen to minimize objective function. In cold chain logistics, from the perspective of cold chain logistics companies, the minimum cost is particularly important.The proposed method establishes a model with the goal of minimizing the total cost, and the total cost is equal to the sum of various costs. In the study of low-carbon and environmentally friendly fresh agricultural products cold chain distribution location-routing problem, cold chain companies choose a suitable cold storage from the existing alternative cold storage as distribution centers and rationally arrange the distribution path for each demand point. The total cost considers the nature of fresh agricultural products and carbon trading mechanism in addition to fixed usage cost of cold storage, refrigeration cost of cold storage, fixed usage cost of vehicles, and vehicle distribution costs. The cost of damage to fresh agricultural products and carbon transaction cost of enterprises should also be considered.(1)
Fixed usage costC2 of cold storageChoose a suitable cold storage from available cold storage options as distribution centers for fresh agricultural products. The distribution center mainly completes simple processing and storage operations of fresh agricultural products. When choosing a cold storage, an enterprise should focus on the geographic location of cold storage. A reasonable geographical location can save delivery time and reduce the quality loss of fresh agricultural products. In addition, the temperature and humidity and whether there is a complete operating system should also be considered.Factors such as the location, temperature, and humidity, and the integrity of supporting devices of alternative cold storage will also affect the rental price of cold storage. Generally speaking, fixed usage cost of cold storage is a constant, and it is only related to the unit fixed usage cost and storage capacity of cold storage. Assuming that the unit fixed usage cost of alternative cold storagem is Fm and the storage capacity is Am, the fixed usage cost C2 of cold storage is calculated as follows:(3)C2=∑m∈MxmFmAm,xm=1,select candidate distribution centerm,0,do not select candidate distribution centerm,m⊂M,whereM is the set of candidate distribution centers and m represents candidate distribution centers, m=1,2,…,M.(2)
Cooling costC3 of cold storageIn order to ensure the quality of fresh agricultural products, it is necessary to maintain a low temperature environment during simple processing and storage. In order to control the temperature, refrigeration devices are needed for cooling. The cost of refrigeration in cold storage mainly comes from electricity consumed by refrigeration devices. The power consumed by refrigeration devices is related to the temperature and capacity of cold storage. Assuming that refrigeration devices used in cold storage are all refrigerators of the same model with powerP, the power unit is W, the number of refrigerators is qm, and the electricity price is set to p∗. The average working time is expressed by tm, and the unit is h. Then, the cooling cost C3 of cold storage is calculated as follows: (4)C3=p∗×∑m∈MxmPqmtm.(3)
Fixed usage costC4 of vehiclesThe fixed usage cost of vehicles refers to related expenses such as daily maintenance, maintenance, and depreciation of vehicles, as well as the driver’s salary and labor cost of other related personnel. The fixed cost of vehicles is generally a constant, and this cost has nothing to do with factors such as the delivery distance and load capacity of vehicles. Regardless of multiple models, it is assumed that the same model of refrigerated vehicle is used for the distribution of fresh agricultural products. The fixed usage cost of each refrigerated vehicle isfk, and then the fixed usage cost C4 of vehicles is calculated as follows: (5)C4=∑i∈M∑j∈C∑k∈Vfkyijk,yijk=1,vehiclekreachesjfrom demand pointi,0,otherwise,j⊂M,i⊂N,whereV is the collection of vehicles, and vehicles are represented by k, k=1,2,…,K; C=M∪N. N is the set of demand points, and n represents demand points, n=1,2,…,N.(4)
Distribution costC5 of vehiclesThe cost of vehicle distribution is mainly variable cost of refrigerated vehicles in the distribution process. This cost is closely related to factors such as delivery distance and cargo capacity. Among them, the vehicle distribution cost includes two parts, the transportation cost of vehicles and cooling cost of vehicles. The transportation cost of vehicles is mainly caused by the fuel consumption in distribution process. This part of cost is positively related to the transportation distance. UseC5′ to represent the transportation cost of vehicles, which is calculated as follows: (6)C5′=∑k∈V∑i∈C∑j∈CCkyijkdij,wheredij represents the distance between node i and node j and Ck represents the fixed usage cost of vehicle k.Fresh agricultural products are perishable. In the entire distribution process of refrigerated vehicles, in order to ensure the quality of fresh agricultural products, the distribution environment needs to be maintained at a suitable temperature. The refrigeration cost of a refrigerated vehicle is the cost of refrigerant consumed to maintain the temperature inside vehicles. The refrigerant consumption of refrigerated vehicles in distribution process can be determined by the heat load. In general, the temperature inside the refrigerated vehicle needs to be higher than the outdoor temperature. The heat load is mainly the heat exchange by compartment wall, the heat exchange caused by the air exchange in the compartment, and the heat load caused by solar radiation. During the distribution process, the heat load coefficientZk1 generated by refrigerated vehicle k is calculated as follows:(7)Zk1=1+α×β×S×tw−tn,whereα represents the deterioration degree coefficient of carriage. The greater the degree of aging and leakage of carriage, the greater the heat conduction caused. β represents the heat transfer coefficient of carriage, which refers to the heat conduction per unit time per unit temperature difference per unit area, and the unit is kcal/h∗m2∗C°. S represents the heat transfer area of carriage, which is equal to the arithmetic square root of the inner and outer surface areas of carriage, that is, S=swsn, where sw and sn represent the surface area of outer compartment and surface area of inner compartment, respectively. tw and tn represent the outdoor temperature and temperature inside the cabin, respectively.In the distribution process of refrigerated vehicles, take one of routes as an example, as shown in Figure1. Suppose that the refrigerated vehicle starts at 0 point of distribution centers, serves demand point 1, demand point 2, and demand point 3 in turn, and finally returns to distribution center 0. In this process, the distribution center 0-demand point 1, demand point 1-demand point 2, and demand point 2-demand point 3 contain products and require refrigeration. The demand point 3-distribution center 0 section is offloaded, and this section does not require refrigeration.UseC5″ to represent the cooling cost of vehicles due to heat exchange, which is calculated as follows:(8)C5″=p5∑k∈VZkl×∑i∈C∑j∈Nyijktij,tijk=dijv,wherep5 is the unit refrigeration cost; ∑i∈C∑j∈Nyijktij represents the total time taken by refrigerated vehicle k from distribution center to the last demand point of its service. The average speed of vehicle k is V, and tijk represents the time for vehicle k from node i to node j.When the vehicle reaches demand points, the door needs to be opened for offloading. At this time, the convection of air outside the vehicle and air inside the vehicle generates a heat load. The heat load of vehiclek when the door is opened is calculated as follows:(9)Zk2=γ×0.55Vk+3.12×tw−tn,whereγ is the door opening frequency coefficient, taking 0.54, and Vk is the volume of the refrigerated vehicle.UseC5‴ to represent the cooling cost of vehicles due to opening the door, which is calculated as follows:(10)C5‴=p5∑k∈VZk2×∑i∈N∑j∈Cyijktsi,wheretsi represents the service time of vehicle k at node i.Then, the vehicle distribution costC5 is calculated as follows:(11)C5=C5′+C5″+C5‴=∑k∈V∑i∈C∑j∈CCkyijkdij+p5∑k∈VZk1×∑i∈C∑j∈Nyijktij+p5∑k∈VZk2×∑i∈C∑j∈Nyijktsi.(5)
Total cargo damage costC6It is assumed that the temperature is constant during the distribution process. Therefore, the cargo damage cost of fresh agricultural products is mainly related to transportation distance. With the accumulation of transportation time, the quality of fresh agricultural products gradually declines. This has led to a corresponding reduction in its transaction price. According to the analysis of deterioration rate function, the cost of cargo damage during transportation of refrigerated vehicles is calculated as follows:(12)C6′=∑k∈V∑i∈C∑j∈Nyijkpbi1−e−θTik−T0k,wherep represents the price of fresh agricultural products without corruption; bi represents the demand at demand point i; Tik represents the time when vehicle k arrives at node i, assuming that the vehicle departs from distribution centers at time 0. For the convenience of calculation, for i∈M, the value of bi needs to be assumed to be 0.The cost of cargo damage caused by opening the door is calculated as follows:(13)C6″=∑k∈V∑i∈C∑j∈NyijkpWijk1−e−θtsi,whereWijk represents the load capacity of vehicle k when it is transported from demand point i to demand point i.Then, the total cargo damage costC6 is calculated as follows:(14)C6=C6′+C6″=∑k∈V∑i∈C∑j∈Nyijkpbi1−e−θTik−T0k+∑k∈V∑i∈C∑j∈NyijkpWijk1−e−θtsi.(6)
Carbon trading costC7The carbon emission of refrigerated vehicles between nodesi,j is ψCij=ePqijdij, and carbon emission in the whole process can be expressed as(15)C6=pc×ψC−ψ0=pc×∑m∈Mxmqpt∗g1000+∑k∈V∑i∈C∑j∈CeP0+qP∗−P0ψ∗dij−ψ0.Figure 1
Schematic diagram of refrigerated vehicle distribution path.


### 3.2. Model Solution

#### 3.2.1. Algorithm Selection

In the established green logistics location-routing optimization model, there is a phenomenon that cannot be compared due to the conflict between targets and different dimensions. The optimization of one objective may lead to a poor solution for another objective [26]. Since the 1990s, many scholars have begun to focus on researching intelligent optimization algorithms and introducing them into the routing problem of distribution vehicles. Intelligent optimization algorithms can be divided into two categories based on the individual behavior and based on group intelligence according to the group size. The first category includes simulated annealing algorithm and tabu search algorithm. The second category includes GA, ant colony algorithm, and particle swarm algorithm. Algorithms based on individual behavior have memory characteristics, which is deterministic, but also random. Such characteristics are often sporadic, random, and unpredictable. The algorithm based on swarm intelligence is to use the method of simulating learning mechanism and information exchange and sharing between group creatures and gradually select high-quality solutions by iteration. Among them, GA can handle constraints well, obtain global optimal solution, and has strong global search ability. But its convergence speed is slow, the local search ability is weak, the running time is longer, and it is easily affected by parameters [27].Therefore, the adaptive operator and cataclysm operator are introduced into GA. Among them, the catastrophe adaptive GA can adjust the probability of crossover and mutation in GA according to the needs, which reduces the shortcomings of standard GA that is susceptible to parameters, and also reduces the running time [28, 29]. The flow of adaptive catastrophe algorithm is shown in Figure 2.Figure 2
Adaptive catastrophe algorithm flow.


#### 3.2.2. Algorithm Design

GA has the advantages of fast calculation speed, strong optimization ability, and good robustness. However, there are defects in the solution process, such as premature maturity and local trapping. Based on the characteristics of green logistics location-distribution route studied, improvements are made on the basis of traditional GA. By improving the coding, selection, crossover, and mutation methods of genetic operators, the generation of errors is reduced and the disadvantage of premature convergence is improved.(1) Coding Method. The primary problem of using GA is to solve chromosome coding, and an important basis of it is to choose the appropriate way of expression of solution. At present, many scholars have launched research on this. According to the actual situation, different encoding methods are summarized, which can be summarized as binary encoding, real number encoding, integer permutation encoding, and data structure encoding. The most used encoding method is binary encoding. 1 and 0 constitute its symbol set, which, respectively, represents the two states of on and off. This makes subsequent operations easier and easier. However, when solving sequence-related problems such as path planning, the binary encoding method has poor operability. It is very likely that infeasible solutions will often appear, which will bring greater complexity to the future. Thus, in order to improve the genetic efficiency and avoid generating infeasible solutions, the intuitive natural number and integer encoding methods for the path planning problem can well express the order of solutions obtained.In summary, the proposed method is more appropriate to use natural number coding. That is, each individuali is coded as a full array of natural numbers from 1 to n, and each natural number corresponds to the number of customer nodes that need to be served in the vehicle path. Its expression sequence means the actual driving sequence of each logistics vehicle service customer point. Because the starting point and return point of vehicles’ entire mission are distribution centers, that is, the logistics vehicles depart from distribution centers and eventually return to distribution centers after performing tasks at each customer point along the way. Based on this coding method, it can be ensured that each client node will be accessed in a certain order without duplication. It can make the processing of model constraints more convenient and quicker and facilitate the reading and modification of program. If there are n customer demand points, there are k vehicles to provide services. In this way, the chromosome is composed of distribution centers and customer point together, and they are arranged in a row in the chromosome. And, each node corresponds to a locus, and finally a n + k + 1 chromosome is formed. The code of chromosome can be expressed as(16)0,i11,i12,i1s;0,i2s,…,i2r;0,ik1,…,ikn,0,where ikm represents the n customer served by vehicle k. 0 represents a distribution center, and each 0 represents a vehicle distribution path. That is to say, the first vehicle departs from distribution center 0, goes through customers i11, i12, and i1s, and then returns to distribution center 0. The same goes for vehicle k, and the end of coding means that delivery is complete. For example, chromosome 0137029805670, which means that 3 vehicles serve 9 customers, and the corresponding path is as follows:Subpath 1: distribution center—customer 1—customer 3—customer 7—distribution centerSubpath 2: distribution center—customer 2—customer 9—customer 8—distribution centerSubpath 3: distribution center—customer 5—customer 6—customer 7—distribution center(2) Generate Initial Population. GA is different from other algorithms in that it searches for populations. Therefore, the generation of initial population is very important in the entire solution process. The pros and cons of initial population have a greater impact on the operational efficiency and computational complexity of algorithm during the calculation process. The key to determining the quality of population is the generation method of initial population and population size. If the population size is small, the algorithm’s solution process will be affected to a certain extent, and generally speaking, it will not be too good. However, choosing a larger group can avoid the situation where solution obtained is a local optimum to a certain extent. But this time will increase the space complexity, and the calculation amount and calculation time will also increase. Generally, for the case where the chromosome length is not very large, the appropriate value range for the population size is 20–150. But for large-scale solving problems, the size of population should be increased to a certain extent based on actual situation.The chromosome chooses real-valued natural number coding, and the expression of solution is related to the order of arrangement. Thus, when initializing the population, a full array ofPopsize client nodes can be generated to form a two-dimensional matrix of n×Popsize, where n represents the number of client nodes and Popsize represents the number of populations.(3) Determine the Fitness Evaluation Function. Fitness is not only a criterion for judging the quality of an individual but also a key indicator for subsequent operations. The greater the adaptability of solution is, the better it is. Therefore, the possibility of being selected will be greater. Because the proposed model needs to meet the goal of minimizing total distribution cost on the one hand, on the other hand, it is also necessary to consider time window, load, and other constraints of vehicles at any customer node. So, the feasibility of solution in the population is very important. The number of vehicles in route optimization problem is undetermined, and the number of delivery vehicles can be efficiently configured in combination with realistic conditions and related restrictions. Thus, the feasibility of solution can be determined by setting an appropriate fitness function value. At the same time, the feasible solution of a single individual can be calculated, which also represents the driving route and the number of vehicles arranged for logistics vehicle. The specific process is as follows:Define a setV to store the set of vehicles in path and use “0” element to separate different vehicles. Its purpose is to show the distribution of vehicle paths by collection containing “0” elements. This means that there is only one “0” element when the collection is initialized. Then, for single individual i, starting from the first gene that has never participated in accumulation, the delivery volume, pickup volume, and time of customer point corresponding to each gene are accumulated one by one. When the accumulation to a certain customer point does not meet constraint conditions, record the genes of cycle in set V, and add a “0” element at the end to represent the end of the cycle. If there are genes in the individual that are not involved in accumulation, they will be assigned another vehicle for the next cycle operation. Otherwise, the loop calculation ends.Therefore, the fitness function used is expressed as follows:(17)fitnessi=1∑Ci,i=1,2,…,Popsize,where ∑Ci is the comprehensive total cost of individual i, that is, the sum of fixed cost, transportation cost, carbon emission cost, and time penalty cost. The total cost changes with the distance traveled by vehicles, so it is necessary to consider whether a certain customer point reached is the starting point of the vehicles’ route. If so, you need to calculate the distance from distribution centers and the distance between this point and next point. If it is not, you only need to calculate the distance between the node and node behind it. The vehicle departs from the starting point (distribution center) and serves customers on complete path along the way. After serving the last customer, it needs to return to distribution centers. Thus, the total distance of the complete path is calculated by adding the distance from the last customer point to distribution centers to the sum of distances during cycle before the end of each cycle.(4) Selection Process. The selection operation is evolutionary process of GA. It is to select individuals with large fitness values from population to perform crossover, mutation operations, or directly as the next generation of new individuals to form a new population. This process embodies the population evolution law of “natural selection, survival of the fittest.” The selection mechanism adopted by the proposed method is a mixed selection mechanism. That is, first sort the fitness of individuals in population by size. And, according to a certain proportion, the population with higher adaptability is selected to enter the next stage to prevent chromosomes with better adaptability from being mistakenly eliminated. Then, the remaining individuals are selected in a form similar to that of betting market selection. Throughout the process, the probability of an individual being selected will increase with the increase in fitness.The traditional roulette method is simple and easy to operate. However, the result of random operation is that the selection error is relatively large, and sometimes even degradation occurs. This means that there may be individuals with higher fitness levels that are not selected but eliminated, which results in the algorithm’s difficulty in obtaining high-quality solutions and large fluctuations in the subsequent process. Therefore, the proposed method operates by sorted multiround roulette method without replacement. The specific process is shown in Algorithm1.Algorithm 1: Pseudo code of multiround roulette wheel selection method with replacement based on sorting.
Assuming the population size isH, the parent population H′ is generated, which is recorded as H′=a1,a2,…,ai,…,aH, where the value of fai represents the fitness of each individual. The initial state of offspring population is X=.Begin(1)
Calculate the fitness value of all individuals and arrange them in order of size. The ranked populationH^′ is denoted as H^′=b1,b2,…,bi,…,bH, where fbi−1>fbi>fbi+1. At the same time, calculate the individual h with the greatest fitness in the parent population, namely, fah=maxfa1,fa2,…,faH(2)
Calculate the total fitness of all chromosomes inH^′, ∑Hfbii=1,2,…,H(3)
Calculate the probability of single chromosome being selected,Pbi=fbi/∑i=1Hfbi(4)
Calculate the cumulative probability of each chromosome being selectedP˜K, P˜K=∑i=1hPi, respectively(5)
Rotate the wheel to generateH uniformly distributed random numbers r between [0,1](6)
Ifr≤P˜1then(7)
Select chromosomeb1;(8)
Otherwise select the h chromosome bh2≤h≤H so that P˜h−1≤r≤P˜h;(9)
Count theζ value of each interval, ζ1,ζ2,…,ζH, ζi is the number of random numbers representing the i area, i=1,2,…,H;(10)
Take the largestζ value and individual ζj=maxζ1,ζ2,…,ζH corresponding to the interval j where bj is located is the individual Hj selected for this round, there is; Hj=bjotherwiseHj=bmini1,i2,…,ijthere are the sameζvalue interval(11)
IncorporateHj into X, that is, X0=∅Xt=Xt−1∪Hj, to get X=c1,c2,…,cH;(12)
Calculate the individual with the smallest fitness value inX, fci=minfc1,fc2,…,fcH and substitute parent individual h for offspring individual j;(13)
Record all newly selected individuals.EndThe proposed method uses a sort-based multiround roulette method without replacement, which greatly increases the number of random numbers generated. It can prevent better chromosomes from being eliminated and effectively reduce selection errors. And, the higher the fitness of individual, the higher the order of ranking, which increases the probability of fitness individual being selected. Although the traditional selection operation can reduce the difference in individual selection to a certain extent, this does not ensure that the best individual in the current generation of population can enter the next step. Therefore, in the improved operator selection method, after each operation, the solution in child with the worst fitness will be replaced. This allows high-quality solution to enter the next step, so that the convergence efficiency is improved.(5) Crossover Operation. Compared with the traditional crossover method, the partial matching crossover method does not directly exchange the crossover segments of chromosomes, but the first adds the gene segments that need to be crossed before the first gene of the opposite chromosome. Then, one by one and the original individual based on comparison, if it is the same gene, it will be retained, if it is not the same, it will be removed, so a new individual is generated. This method can not only avoid the phenomenon of premature maturity but also does not require the diversity of population. The schematic diagram of partial matching crossover operation is shown in Figure 3.Figure 3
Partial matching crossover operation.Assuming parent A and B, the chromosome codes are “235|4891|67” and “584|3192|67,” respectively, and “|” represents the cross position of chromosomes. When performing the crossover operation, first put the gene segment “4891” of the selected crossover of parent A in front. Then, arrange the chromosome of other parent B in order from the position of the fifth locus. Then, compare the elements in parent B with gene segment “4891” one by one. In this process, the same gene is removed; otherwise, the gene remains unchanged. From this, we can get offspring A as “489143267.” Similarly, child B can be obtained as “319254867.”(6) Mutation Operation. In order to maintain the diversity of individuals in population and avoid the convergence of individuals in later stages of evolution, mutation operations can repair this. The mutation operation can expand the range of local optimization comprehensive target value, thereby enhancing the global search ability of algorithm. The proposed method adopts the mutation method of reversal operation to operate on an individual. First, select individuals in the population according to a certain probability and randomly generate two mutation points in the individual’s genes, resulting in mutation. Then, perform reversal operation on variant segment. The mutation operation process is shown in Figure 4.Figure 4
Mutation operation.Randomly select a certain body in populationA = “156827394” and then randomly generate two mutations “3” and “8,” namely, A = 156|8273|94. Finally, the variant segment was reversed to obtain a new individual A “156372894.”(7) Catastrophe Operation. When there are no better individuals for N consecutive generations or the population diversity change is very low, the algorithm is very likely to fall into a local optimal solution. By using the catastrophic method to kill the optimal individual of population in the next iteration, it will jump out of the local extreme value.The catastrophe judgment process is as follows: letϕ0 be the initial setting value, and the catastrophe counter starts counting from 0. When the catastrophe counter value exceeds ϕ0, it is considered that local search has produced the optimal solution, and the catastrophe is carried out. If the new optimal solution continues to be produced, the catastrophe count value is cleared. In order to prevent insufficient local search during catastrophe, the set value is enlarged and catastrophe coefficient is introduced, and then the catastrophe judgment value ϕ is as follows.Whenϕ is less than the initial setting value, count again. When ϕ is greater than the initial set value, it indicates that the population diversity change is very low. Use this value to update the catastrophe judgment value according to equation (18) and then start the countdown.(8) Termination Conditions. The termination criterion is used to judge whether the genetic operation is over. The termination conditions of algorithm generally have the following forms: the algorithm is executed to the maximum number of iterations, when the adjacent generations are unable to evolve, and the absolute value of the difference between the two target values is less than a precision. The proposed method adopts setting the maximum number of iterations N, when the evolutionary algebra reaches the maximum number of iterations, the algorithm ends. Otherwise, continue to repeat the above process.


## 3.1. Model Construction of Green Logistics Location-Routing Optimization Problem considering Carbon Trading

### 3.1.1. Model Assumption

The following assumptions are made for green logistics location-routing optimization problem considering carbon trading:(1)
Consider the capacity limitation of distribution centers(2)
During offloading time, the deterioration rate will be higher than deterioration rate during transportation.(3)
A demand point can only be delivered by one vehicle, and other vehicles are not allowed to deliver, that is, split delivery is not allowed.(4)
The vehicle runs at a constant speed, and temperature is constant during storage and transportation. After completing the distribution task, it will return to distribution centers empty. In this process, refrigeration is not required.(5)
It is assumed that, as the quality of fresh agricultural products decreases, the price will also decrease correspondingly, resulting in a certain loss of goods cost.(6)
For a given amount of carbon emission allowance, if the cold chain enterprise’s emissions do not exceed allowance, the part less than allowance can be converted into revenue at carbon trading price to offset part of the cost. If the cold chain company’s emissions exceed the quota, it needs to purchase additional quotas from carbon trading market, incurring additional carbon trading costs.


### 3.1.2. Function Analysis

(1) Metamorphosis Rate Function. In the process of storage and distribution, the quality of fresh agricultural products will change with time and temperature. People’s requirements for food quality are getting higher and higher. When the quality of fresh agricultural products declines, sales price will also decline accordingly. Assuming that the temperature during storage and distribution is constant and strictly controlled, the spoilage rate of fresh agricultural products at this time is only related to time. Due to short loading and offloading time, the temperature change when opening the door is negligible. Therefore, at this constant temperature, the deterioration rate can be regarded as a constant. The quality of fresh produce can be regarded as exponentially changing over time.Assume thatQ0 is the quality of fresh agricultural products at time 0, that is, the quality of fresh agricultural products when the refrigerated vehicle departs from selected distribution centers. Use θ to represent the spoilage rate of fresh agricultural products. At a constant temperature, the spoilage rate of a fresh agricultural product is only related to time. If Qt is used to represent the quality of fresh agricultural products at time t, the mathematical expression is(1)Qt=Q0e−θt.(2) Carbon Trading Function. Considering the location-routing optimization problem of cold chain distribution of fresh agricultural products under carbon trading mechanism, the total amount of carbon emissions is obtained. In order to study the impact of carbon trading mechanism on the location-routing problem of cold chain distribution of fresh agricultural products, the total carbon emissions need to be monetized and converted into carbon trading costs. The conversion of total carbon emissions to carbon trading costs is achieved through carbon trading functions. The size of carbon trading costs depends on the total carbon emissions, carbon trading costs, and initial carbon allowances [21–25].When the company’s actual total carbon emissions are less than initial quota, the part that is less than the quota can be sold to carbon trading agency. This can yield benefits and offset some of the costs. When the company’s actual total carbon emissions are higher than initial allowance, it will need to spend a certain cost to purchase the carbon emissions allowance that is higher than initial allowance. Thus, the carbon trading function can be expressed as(2)C1=pc×ψC−ψ0=pc×ψCT,where pc represents the carbon transaction price, that is, the transaction price of unit carbon emission allowance; ψC represents the actual total carbon emission of enterprise; ψ0 represents the initial carbon emission allowance of enterprise; and ψCT is used to denote the company’s carbon trading volume, which is the difference between company’s actual total carbon emissions minus initial carbon emission allowances. The value of ψCT may be positive or negative. When ψCT≥0, it means that the company’s actual total carbon emissions are higher than initial carbon emissions quota. At this time, ψCT represents the carbon trading quota that the company has purchased from carbon trading market. When ψCT<0, it means that the company’s actual total carbon emission is lower than initial carbon emission quota. At this time, ψCT represents the carbon trading quota sold by the company to carbon trading market.(3) Objective Function. For the objective function of location-routing problem, the most common ones are cost minimization, time minimization, and profit maximization. For different situations, the selected objective function is also quite different. For the location-routing problem of distribution center in green logistics, considering the weak cost of green logistics, the time is generally chosen to minimize objective function. In cold chain logistics, from the perspective of cold chain logistics companies, the minimum cost is particularly important.The proposed method establishes a model with the goal of minimizing the total cost, and the total cost is equal to the sum of various costs. In the study of low-carbon and environmentally friendly fresh agricultural products cold chain distribution location-routing problem, cold chain companies choose a suitable cold storage from the existing alternative cold storage as distribution centers and rationally arrange the distribution path for each demand point. The total cost considers the nature of fresh agricultural products and carbon trading mechanism in addition to fixed usage cost of cold storage, refrigeration cost of cold storage, fixed usage cost of vehicles, and vehicle distribution costs. The cost of damage to fresh agricultural products and carbon transaction cost of enterprises should also be considered.(1)
Fixed usage costC2 of cold storageChoose a suitable cold storage from available cold storage options as distribution centers for fresh agricultural products. The distribution center mainly completes simple processing and storage operations of fresh agricultural products. When choosing a cold storage, an enterprise should focus on the geographic location of cold storage. A reasonable geographical location can save delivery time and reduce the quality loss of fresh agricultural products. In addition, the temperature and humidity and whether there is a complete operating system should also be considered.Factors such as the location, temperature, and humidity, and the integrity of supporting devices of alternative cold storage will also affect the rental price of cold storage. Generally speaking, fixed usage cost of cold storage is a constant, and it is only related to the unit fixed usage cost and storage capacity of cold storage. Assuming that the unit fixed usage cost of alternative cold storagem is Fm and the storage capacity is Am, the fixed usage cost C2 of cold storage is calculated as follows:(3)C2=∑m∈MxmFmAm,xm=1,select candidate distribution centerm,0,do not select candidate distribution centerm,m⊂M,whereM is the set of candidate distribution centers and m represents candidate distribution centers, m=1,2,…,M.(2)
Cooling costC3 of cold storageIn order to ensure the quality of fresh agricultural products, it is necessary to maintain a low temperature environment during simple processing and storage. In order to control the temperature, refrigeration devices are needed for cooling. The cost of refrigeration in cold storage mainly comes from electricity consumed by refrigeration devices. The power consumed by refrigeration devices is related to the temperature and capacity of cold storage. Assuming that refrigeration devices used in cold storage are all refrigerators of the same model with powerP, the power unit is W, the number of refrigerators is qm, and the electricity price is set to p∗. The average working time is expressed by tm, and the unit is h. Then, the cooling cost C3 of cold storage is calculated as follows: (4)C3=p∗×∑m∈MxmPqmtm.(3)
Fixed usage costC4 of vehiclesThe fixed usage cost of vehicles refers to related expenses such as daily maintenance, maintenance, and depreciation of vehicles, as well as the driver’s salary and labor cost of other related personnel. The fixed cost of vehicles is generally a constant, and this cost has nothing to do with factors such as the delivery distance and load capacity of vehicles. Regardless of multiple models, it is assumed that the same model of refrigerated vehicle is used for the distribution of fresh agricultural products. The fixed usage cost of each refrigerated vehicle isfk, and then the fixed usage cost C4 of vehicles is calculated as follows: (5)C4=∑i∈M∑j∈C∑k∈Vfkyijk,yijk=1,vehiclekreachesjfrom demand pointi,0,otherwise,j⊂M,i⊂N,whereV is the collection of vehicles, and vehicles are represented by k, k=1,2,…,K; C=M∪N. N is the set of demand points, and n represents demand points, n=1,2,…,N.(4)
Distribution costC5 of vehiclesThe cost of vehicle distribution is mainly variable cost of refrigerated vehicles in the distribution process. This cost is closely related to factors such as delivery distance and cargo capacity. Among them, the vehicle distribution cost includes two parts, the transportation cost of vehicles and cooling cost of vehicles. The transportation cost of vehicles is mainly caused by the fuel consumption in distribution process. This part of cost is positively related to the transportation distance. UseC5′ to represent the transportation cost of vehicles, which is calculated as follows: (6)C5′=∑k∈V∑i∈C∑j∈CCkyijkdij,wheredij represents the distance between node i and node j and Ck represents the fixed usage cost of vehicle k.Fresh agricultural products are perishable. In the entire distribution process of refrigerated vehicles, in order to ensure the quality of fresh agricultural products, the distribution environment needs to be maintained at a suitable temperature. The refrigeration cost of a refrigerated vehicle is the cost of refrigerant consumed to maintain the temperature inside vehicles. The refrigerant consumption of refrigerated vehicles in distribution process can be determined by the heat load. In general, the temperature inside the refrigerated vehicle needs to be higher than the outdoor temperature. The heat load is mainly the heat exchange by compartment wall, the heat exchange caused by the air exchange in the compartment, and the heat load caused by solar radiation. During the distribution process, the heat load coefficientZk1 generated by refrigerated vehicle k is calculated as follows:(7)Zk1=1+α×β×S×tw−tn,whereα represents the deterioration degree coefficient of carriage. The greater the degree of aging and leakage of carriage, the greater the heat conduction caused. β represents the heat transfer coefficient of carriage, which refers to the heat conduction per unit time per unit temperature difference per unit area, and the unit is kcal/h∗m2∗C°. S represents the heat transfer area of carriage, which is equal to the arithmetic square root of the inner and outer surface areas of carriage, that is, S=swsn, where sw and sn represent the surface area of outer compartment and surface area of inner compartment, respectively. tw and tn represent the outdoor temperature and temperature inside the cabin, respectively.In the distribution process of refrigerated vehicles, take one of routes as an example, as shown in Figure1. Suppose that the refrigerated vehicle starts at 0 point of distribution centers, serves demand point 1, demand point 2, and demand point 3 in turn, and finally returns to distribution center 0. In this process, the distribution center 0-demand point 1, demand point 1-demand point 2, and demand point 2-demand point 3 contain products and require refrigeration. The demand point 3-distribution center 0 section is offloaded, and this section does not require refrigeration.UseC5″ to represent the cooling cost of vehicles due to heat exchange, which is calculated as follows:(8)C5″=p5∑k∈VZkl×∑i∈C∑j∈Nyijktij,tijk=dijv,wherep5 is the unit refrigeration cost; ∑i∈C∑j∈Nyijktij represents the total time taken by refrigerated vehicle k from distribution center to the last demand point of its service. The average speed of vehicle k is V, and tijk represents the time for vehicle k from node i to node j.When the vehicle reaches demand points, the door needs to be opened for offloading. At this time, the convection of air outside the vehicle and air inside the vehicle generates a heat load. The heat load of vehiclek when the door is opened is calculated as follows:(9)Zk2=γ×0.55Vk+3.12×tw−tn,whereγ is the door opening frequency coefficient, taking 0.54, and Vk is the volume of the refrigerated vehicle.UseC5‴ to represent the cooling cost of vehicles due to opening the door, which is calculated as follows:(10)C5‴=p5∑k∈VZk2×∑i∈N∑j∈Cyijktsi,wheretsi represents the service time of vehicle k at node i.Then, the vehicle distribution costC5 is calculated as follows:(11)C5=C5′+C5″+C5‴=∑k∈V∑i∈C∑j∈CCkyijkdij+p5∑k∈VZk1×∑i∈C∑j∈Nyijktij+p5∑k∈VZk2×∑i∈C∑j∈Nyijktsi.(5)
Total cargo damage costC6It is assumed that the temperature is constant during the distribution process. Therefore, the cargo damage cost of fresh agricultural products is mainly related to transportation distance. With the accumulation of transportation time, the quality of fresh agricultural products gradually declines. This has led to a corresponding reduction in its transaction price. According to the analysis of deterioration rate function, the cost of cargo damage during transportation of refrigerated vehicles is calculated as follows:(12)C6′=∑k∈V∑i∈C∑j∈Nyijkpbi1−e−θTik−T0k,wherep represents the price of fresh agricultural products without corruption; bi represents the demand at demand point i; Tik represents the time when vehicle k arrives at node i, assuming that the vehicle departs from distribution centers at time 0. For the convenience of calculation, for i∈M, the value of bi needs to be assumed to be 0.The cost of cargo damage caused by opening the door is calculated as follows:(13)C6″=∑k∈V∑i∈C∑j∈NyijkpWijk1−e−θtsi,whereWijk represents the load capacity of vehicle k when it is transported from demand point i to demand point i.Then, the total cargo damage costC6 is calculated as follows:(14)C6=C6′+C6″=∑k∈V∑i∈C∑j∈Nyijkpbi1−e−θTik−T0k+∑k∈V∑i∈C∑j∈NyijkpWijk1−e−θtsi.(6)
Carbon trading costC7The carbon emission of refrigerated vehicles between nodesi,j is ψCij=ePqijdij, and carbon emission in the whole process can be expressed as(15)C6=pc×ψC−ψ0=pc×∑m∈Mxmqpt∗g1000+∑k∈V∑i∈C∑j∈CeP0+qP∗−P0ψ∗dij−ψ0.Figure 1
Schematic diagram of refrigerated vehicle distribution path.


## 3.1.1. Model Assumption

The following assumptions are made for green logistics location-routing optimization problem considering carbon trading:(1)
Consider the capacity limitation of distribution centers(2)
During offloading time, the deterioration rate will be higher than deterioration rate during transportation.(3)
A demand point can only be delivered by one vehicle, and other vehicles are not allowed to deliver, that is, split delivery is not allowed.(4)
The vehicle runs at a constant speed, and temperature is constant during storage and transportation. After completing the distribution task, it will return to distribution centers empty. In this process, refrigeration is not required.(5)
It is assumed that, as the quality of fresh agricultural products decreases, the price will also decrease correspondingly, resulting in a certain loss of goods cost.(6)
For a given amount of carbon emission allowance, if the cold chain enterprise’s emissions do not exceed allowance, the part less than allowance can be converted into revenue at carbon trading price to offset part of the cost. If the cold chain company’s emissions exceed the quota, it needs to purchase additional quotas from carbon trading market, incurring additional carbon trading costs.


## 3.1.2. Function Analysis

(1) Metamorphosis Rate Function. In the process of storage and distribution, the quality of fresh agricultural products will change with time and temperature. People’s requirements for food quality are getting higher and higher. When the quality of fresh agricultural products declines, sales price will also decline accordingly. Assuming that the temperature during storage and distribution is constant and strictly controlled, the spoilage rate of fresh agricultural products at this time is only related to time. Due to short loading and offloading time, the temperature change when opening the door is negligible. Therefore, at this constant temperature, the deterioration rate can be regarded as a constant. The quality of fresh produce can be regarded as exponentially changing over time.Assume thatQ0 is the quality of fresh agricultural products at time 0, that is, the quality of fresh agricultural products when the refrigerated vehicle departs from selected distribution centers. Use θ to represent the spoilage rate of fresh agricultural products. At a constant temperature, the spoilage rate of a fresh agricultural product is only related to time. If Qt is used to represent the quality of fresh agricultural products at time t, the mathematical expression is(1)Qt=Q0e−θt.(2) Carbon Trading Function. Considering the location-routing optimization problem of cold chain distribution of fresh agricultural products under carbon trading mechanism, the total amount of carbon emissions is obtained. In order to study the impact of carbon trading mechanism on the location-routing problem of cold chain distribution of fresh agricultural products, the total carbon emissions need to be monetized and converted into carbon trading costs. The conversion of total carbon emissions to carbon trading costs is achieved through carbon trading functions. The size of carbon trading costs depends on the total carbon emissions, carbon trading costs, and initial carbon allowances [21–25].When the company’s actual total carbon emissions are less than initial quota, the part that is less than the quota can be sold to carbon trading agency. This can yield benefits and offset some of the costs. When the company’s actual total carbon emissions are higher than initial allowance, it will need to spend a certain cost to purchase the carbon emissions allowance that is higher than initial allowance. Thus, the carbon trading function can be expressed as(2)C1=pc×ψC−ψ0=pc×ψCT,where pc represents the carbon transaction price, that is, the transaction price of unit carbon emission allowance; ψC represents the actual total carbon emission of enterprise; ψ0 represents the initial carbon emission allowance of enterprise; and ψCT is used to denote the company’s carbon trading volume, which is the difference between company’s actual total carbon emissions minus initial carbon emission allowances. The value of ψCT may be positive or negative. When ψCT≥0, it means that the company’s actual total carbon emissions are higher than initial carbon emissions quota. At this time, ψCT represents the carbon trading quota that the company has purchased from carbon trading market. When ψCT<0, it means that the company’s actual total carbon emission is lower than initial carbon emission quota. At this time, ψCT represents the carbon trading quota sold by the company to carbon trading market.(3) Objective Function. For the objective function of location-routing problem, the most common ones are cost minimization, time minimization, and profit maximization. For different situations, the selected objective function is also quite different. For the location-routing problem of distribution center in green logistics, considering the weak cost of green logistics, the time is generally chosen to minimize objective function. In cold chain logistics, from the perspective of cold chain logistics companies, the minimum cost is particularly important.The proposed method establishes a model with the goal of minimizing the total cost, and the total cost is equal to the sum of various costs. In the study of low-carbon and environmentally friendly fresh agricultural products cold chain distribution location-routing problem, cold chain companies choose a suitable cold storage from the existing alternative cold storage as distribution centers and rationally arrange the distribution path for each demand point. The total cost considers the nature of fresh agricultural products and carbon trading mechanism in addition to fixed usage cost of cold storage, refrigeration cost of cold storage, fixed usage cost of vehicles, and vehicle distribution costs. The cost of damage to fresh agricultural products and carbon transaction cost of enterprises should also be considered.(1)
Fixed usage costC2 of cold storageChoose a suitable cold storage from available cold storage options as distribution centers for fresh agricultural products. The distribution center mainly completes simple processing and storage operations of fresh agricultural products. When choosing a cold storage, an enterprise should focus on the geographic location of cold storage. A reasonable geographical location can save delivery time and reduce the quality loss of fresh agricultural products. In addition, the temperature and humidity and whether there is a complete operating system should also be considered.Factors such as the location, temperature, and humidity, and the integrity of supporting devices of alternative cold storage will also affect the rental price of cold storage. Generally speaking, fixed usage cost of cold storage is a constant, and it is only related to the unit fixed usage cost and storage capacity of cold storage. Assuming that the unit fixed usage cost of alternative cold storagem is Fm and the storage capacity is Am, the fixed usage cost C2 of cold storage is calculated as follows:(3)C2=∑m∈MxmFmAm,xm=1,select candidate distribution centerm,0,do not select candidate distribution centerm,m⊂M,whereM is the set of candidate distribution centers and m represents candidate distribution centers, m=1,2,…,M.(2)
Cooling costC3 of cold storageIn order to ensure the quality of fresh agricultural products, it is necessary to maintain a low temperature environment during simple processing and storage. In order to control the temperature, refrigeration devices are needed for cooling. The cost of refrigeration in cold storage mainly comes from electricity consumed by refrigeration devices. The power consumed by refrigeration devices is related to the temperature and capacity of cold storage. Assuming that refrigeration devices used in cold storage are all refrigerators of the same model with powerP, the power unit is W, the number of refrigerators is qm, and the electricity price is set to p∗. The average working time is expressed by tm, and the unit is h. Then, the cooling cost C3 of cold storage is calculated as follows: (4)C3=p∗×∑m∈MxmPqmtm.(3)
Fixed usage costC4 of vehiclesThe fixed usage cost of vehicles refers to related expenses such as daily maintenance, maintenance, and depreciation of vehicles, as well as the driver’s salary and labor cost of other related personnel. The fixed cost of vehicles is generally a constant, and this cost has nothing to do with factors such as the delivery distance and load capacity of vehicles. Regardless of multiple models, it is assumed that the same model of refrigerated vehicle is used for the distribution of fresh agricultural products. The fixed usage cost of each refrigerated vehicle isfk, and then the fixed usage cost C4 of vehicles is calculated as follows: (5)C4=∑i∈M∑j∈C∑k∈Vfkyijk,yijk=1,vehiclekreachesjfrom demand pointi,0,otherwise,j⊂M,i⊂N,whereV is the collection of vehicles, and vehicles are represented by k, k=1,2,…,K; C=M∪N. N is the set of demand points, and n represents demand points, n=1,2,…,N.(4)
Distribution costC5 of vehiclesThe cost of vehicle distribution is mainly variable cost of refrigerated vehicles in the distribution process. This cost is closely related to factors such as delivery distance and cargo capacity. Among them, the vehicle distribution cost includes two parts, the transportation cost of vehicles and cooling cost of vehicles. The transportation cost of vehicles is mainly caused by the fuel consumption in distribution process. This part of cost is positively related to the transportation distance. UseC5′ to represent the transportation cost of vehicles, which is calculated as follows: (6)C5′=∑k∈V∑i∈C∑j∈CCkyijkdij,wheredij represents the distance between node i and node j and Ck represents the fixed usage cost of vehicle k.Fresh agricultural products are perishable. In the entire distribution process of refrigerated vehicles, in order to ensure the quality of fresh agricultural products, the distribution environment needs to be maintained at a suitable temperature. The refrigeration cost of a refrigerated vehicle is the cost of refrigerant consumed to maintain the temperature inside vehicles. The refrigerant consumption of refrigerated vehicles in distribution process can be determined by the heat load. In general, the temperature inside the refrigerated vehicle needs to be higher than the outdoor temperature. The heat load is mainly the heat exchange by compartment wall, the heat exchange caused by the air exchange in the compartment, and the heat load caused by solar radiation. During the distribution process, the heat load coefficientZk1 generated by refrigerated vehicle k is calculated as follows:(7)Zk1=1+α×β×S×tw−tn,whereα represents the deterioration degree coefficient of carriage. The greater the degree of aging and leakage of carriage, the greater the heat conduction caused. β represents the heat transfer coefficient of carriage, which refers to the heat conduction per unit time per unit temperature difference per unit area, and the unit is kcal/h∗m2∗C°. S represents the heat transfer area of carriage, which is equal to the arithmetic square root of the inner and outer surface areas of carriage, that is, S=swsn, where sw and sn represent the surface area of outer compartment and surface area of inner compartment, respectively. tw and tn represent the outdoor temperature and temperature inside the cabin, respectively.In the distribution process of refrigerated vehicles, take one of routes as an example, as shown in Figure1. Suppose that the refrigerated vehicle starts at 0 point of distribution centers, serves demand point 1, demand point 2, and demand point 3 in turn, and finally returns to distribution center 0. In this process, the distribution center 0-demand point 1, demand point 1-demand point 2, and demand point 2-demand point 3 contain products and require refrigeration. The demand point 3-distribution center 0 section is offloaded, and this section does not require refrigeration.UseC5″ to represent the cooling cost of vehicles due to heat exchange, which is calculated as follows:(8)C5″=p5∑k∈VZkl×∑i∈C∑j∈Nyijktij,tijk=dijv,wherep5 is the unit refrigeration cost; ∑i∈C∑j∈Nyijktij represents the total time taken by refrigerated vehicle k from distribution center to the last demand point of its service. The average speed of vehicle k is V, and tijk represents the time for vehicle k from node i to node j.When the vehicle reaches demand points, the door needs to be opened for offloading. At this time, the convection of air outside the vehicle and air inside the vehicle generates a heat load. The heat load of vehiclek when the door is opened is calculated as follows:(9)Zk2=γ×0.55Vk+3.12×tw−tn,whereγ is the door opening frequency coefficient, taking 0.54, and Vk is the volume of the refrigerated vehicle.UseC5‴ to represent the cooling cost of vehicles due to opening the door, which is calculated as follows:(10)C5‴=p5∑k∈VZk2×∑i∈N∑j∈Cyijktsi,wheretsi represents the service time of vehicle k at node i.Then, the vehicle distribution costC5 is calculated as follows:(11)C5=C5′+C5″+C5‴=∑k∈V∑i∈C∑j∈CCkyijkdij+p5∑k∈VZk1×∑i∈C∑j∈Nyijktij+p5∑k∈VZk2×∑i∈C∑j∈Nyijktsi.(5)
Total cargo damage costC6It is assumed that the temperature is constant during the distribution process. Therefore, the cargo damage cost of fresh agricultural products is mainly related to transportation distance. With the accumulation of transportation time, the quality of fresh agricultural products gradually declines. This has led to a corresponding reduction in its transaction price. According to the analysis of deterioration rate function, the cost of cargo damage during transportation of refrigerated vehicles is calculated as follows:(12)C6′=∑k∈V∑i∈C∑j∈Nyijkpbi1−e−θTik−T0k,wherep represents the price of fresh agricultural products without corruption; bi represents the demand at demand point i; Tik represents the time when vehicle k arrives at node i, assuming that the vehicle departs from distribution centers at time 0. For the convenience of calculation, for i∈M, the value of bi needs to be assumed to be 0.The cost of cargo damage caused by opening the door is calculated as follows:(13)C6″=∑k∈V∑i∈C∑j∈NyijkpWijk1−e−θtsi,whereWijk represents the load capacity of vehicle k when it is transported from demand point i to demand point i.Then, the total cargo damage costC6 is calculated as follows:(14)C6=C6′+C6″=∑k∈V∑i∈C∑j∈Nyijkpbi1−e−θTik−T0k+∑k∈V∑i∈C∑j∈NyijkpWijk1−e−θtsi.(6)
Carbon trading costC7The carbon emission of refrigerated vehicles between nodesi,j is ψCij=ePqijdij, and carbon emission in the whole process can be expressed as(15)C6=pc×ψC−ψ0=pc×∑m∈Mxmqpt∗g1000+∑k∈V∑i∈C∑j∈CeP0+qP∗−P0ψ∗dij−ψ0.Figure 1
Schematic diagram of refrigerated vehicle distribution path.


## 3.2. Model Solution

### 3.2.1. Algorithm Selection

In the established green logistics location-routing optimization model, there is a phenomenon that cannot be compared due to the conflict between targets and different dimensions. The optimization of one objective may lead to a poor solution for another objective [26]. Since the 1990s, many scholars have begun to focus on researching intelligent optimization algorithms and introducing them into the routing problem of distribution vehicles. Intelligent optimization algorithms can be divided into two categories based on the individual behavior and based on group intelligence according to the group size. The first category includes simulated annealing algorithm and tabu search algorithm. The second category includes GA, ant colony algorithm, and particle swarm algorithm. Algorithms based on individual behavior have memory characteristics, which is deterministic, but also random. Such characteristics are often sporadic, random, and unpredictable. The algorithm based on swarm intelligence is to use the method of simulating learning mechanism and information exchange and sharing between group creatures and gradually select high-quality solutions by iteration. Among them, GA can handle constraints well, obtain global optimal solution, and has strong global search ability. But its convergence speed is slow, the local search ability is weak, the running time is longer, and it is easily affected by parameters [27].Therefore, the adaptive operator and cataclysm operator are introduced into GA. Among them, the catastrophe adaptive GA can adjust the probability of crossover and mutation in GA according to the needs, which reduces the shortcomings of standard GA that is susceptible to parameters, and also reduces the running time [28, 29]. The flow of adaptive catastrophe algorithm is shown in Figure 2.Figure 2
Adaptive catastrophe algorithm flow.


### 3.2.2. Algorithm Design

GA has the advantages of fast calculation speed, strong optimization ability, and good robustness. However, there are defects in the solution process, such as premature maturity and local trapping. Based on the characteristics of green logistics location-distribution route studied, improvements are made on the basis of traditional GA. By improving the coding, selection, crossover, and mutation methods of genetic operators, the generation of errors is reduced and the disadvantage of premature convergence is improved.(1) Coding Method. The primary problem of using GA is to solve chromosome coding, and an important basis of it is to choose the appropriate way of expression of solution. At present, many scholars have launched research on this. According to the actual situation, different encoding methods are summarized, which can be summarized as binary encoding, real number encoding, integer permutation encoding, and data structure encoding. The most used encoding method is binary encoding. 1 and 0 constitute its symbol set, which, respectively, represents the two states of on and off. This makes subsequent operations easier and easier. However, when solving sequence-related problems such as path planning, the binary encoding method has poor operability. It is very likely that infeasible solutions will often appear, which will bring greater complexity to the future. Thus, in order to improve the genetic efficiency and avoid generating infeasible solutions, the intuitive natural number and integer encoding methods for the path planning problem can well express the order of solutions obtained.In summary, the proposed method is more appropriate to use natural number coding. That is, each individuali is coded as a full array of natural numbers from 1 to n, and each natural number corresponds to the number of customer nodes that need to be served in the vehicle path. Its expression sequence means the actual driving sequence of each logistics vehicle service customer point. Because the starting point and return point of vehicles’ entire mission are distribution centers, that is, the logistics vehicles depart from distribution centers and eventually return to distribution centers after performing tasks at each customer point along the way. Based on this coding method, it can be ensured that each client node will be accessed in a certain order without duplication. It can make the processing of model constraints more convenient and quicker and facilitate the reading and modification of program. If there are n customer demand points, there are k vehicles to provide services. In this way, the chromosome is composed of distribution centers and customer point together, and they are arranged in a row in the chromosome. And, each node corresponds to a locus, and finally a n + k + 1 chromosome is formed. The code of chromosome can be expressed as(16)0,i11,i12,i1s;0,i2s,…,i2r;0,ik1,…,ikn,0,where ikm represents the n customer served by vehicle k. 0 represents a distribution center, and each 0 represents a vehicle distribution path. That is to say, the first vehicle departs from distribution center 0, goes through customers i11, i12, and i1s, and then returns to distribution center 0. The same goes for vehicle k, and the end of coding means that delivery is complete. For example, chromosome 0137029805670, which means that 3 vehicles serve 9 customers, and the corresponding path is as follows:Subpath 1: distribution center—customer 1—customer 3—customer 7—distribution centerSubpath 2: distribution center—customer 2—customer 9—customer 8—distribution centerSubpath 3: distribution center—customer 5—customer 6—customer 7—distribution center(2) Generate Initial Population. GA is different from other algorithms in that it searches for populations. Therefore, the generation of initial population is very important in the entire solution process. The pros and cons of initial population have a greater impact on the operational efficiency and computational complexity of algorithm during the calculation process. The key to determining the quality of population is the generation method of initial population and population size. If the population size is small, the algorithm’s solution process will be affected to a certain extent, and generally speaking, it will not be too good. However, choosing a larger group can avoid the situation where solution obtained is a local optimum to a certain extent. But this time will increase the space complexity, and the calculation amount and calculation time will also increase. Generally, for the case where the chromosome length is not very large, the appropriate value range for the population size is 20–150. But for large-scale solving problems, the size of population should be increased to a certain extent based on actual situation.The chromosome chooses real-valued natural number coding, and the expression of solution is related to the order of arrangement. Thus, when initializing the population, a full array ofPopsize client nodes can be generated to form a two-dimensional matrix of n×Popsize, where n represents the number of client nodes and Popsize represents the number of populations.(3) Determine the Fitness Evaluation Function. Fitness is not only a criterion for judging the quality of an individual but also a key indicator for subsequent operations. The greater the adaptability of solution is, the better it is. Therefore, the possibility of being selected will be greater. Because the proposed model needs to meet the goal of minimizing total distribution cost on the one hand, on the other hand, it is also necessary to consider time window, load, and other constraints of vehicles at any customer node. So, the feasibility of solution in the population is very important. The number of vehicles in route optimization problem is undetermined, and the number of delivery vehicles can be efficiently configured in combination with realistic conditions and related restrictions. Thus, the feasibility of solution can be determined by setting an appropriate fitness function value. At the same time, the feasible solution of a single individual can be calculated, which also represents the driving route and the number of vehicles arranged for logistics vehicle. The specific process is as follows:Define a setV to store the set of vehicles in path and use “0” element to separate different vehicles. Its purpose is to show the distribution of vehicle paths by collection containing “0” elements. This means that there is only one “0” element when the collection is initialized. Then, for single individual i, starting from the first gene that has never participated in accumulation, the delivery volume, pickup volume, and time of customer point corresponding to each gene are accumulated one by one. When the accumulation to a certain customer point does not meet constraint conditions, record the genes of cycle in set V, and add a “0” element at the end to represent the end of the cycle. If there are genes in the individual that are not involved in accumulation, they will be assigned another vehicle for the next cycle operation. Otherwise, the loop calculation ends.Therefore, the fitness function used is expressed as follows:(17)fitnessi=1∑Ci,i=1,2,…,Popsize,where ∑Ci is the comprehensive total cost of individual i, that is, the sum of fixed cost, transportation cost, carbon emission cost, and time penalty cost. The total cost changes with the distance traveled by vehicles, so it is necessary to consider whether a certain customer point reached is the starting point of the vehicles’ route. If so, you need to calculate the distance from distribution centers and the distance between this point and next point. If it is not, you only need to calculate the distance between the node and node behind it. The vehicle departs from the starting point (distribution center) and serves customers on complete path along the way. After serving the last customer, it needs to return to distribution centers. Thus, the total distance of the complete path is calculated by adding the distance from the last customer point to distribution centers to the sum of distances during cycle before the end of each cycle.(4) Selection Process. The selection operation is evolutionary process of GA. It is to select individuals with large fitness values from population to perform crossover, mutation operations, or directly as the next generation of new individuals to form a new population. This process embodies the population evolution law of “natural selection, survival of the fittest.” The selection mechanism adopted by the proposed method is a mixed selection mechanism. That is, first sort the fitness of individuals in population by size. And, according to a certain proportion, the population with higher adaptability is selected to enter the next stage to prevent chromosomes with better adaptability from being mistakenly eliminated. Then, the remaining individuals are selected in a form similar to that of betting market selection. Throughout the process, the probability of an individual being selected will increase with the increase in fitness.The traditional roulette method is simple and easy to operate. However, the result of random operation is that the selection error is relatively large, and sometimes even degradation occurs. This means that there may be individuals with higher fitness levels that are not selected but eliminated, which results in the algorithm’s difficulty in obtaining high-quality solutions and large fluctuations in the subsequent process. Therefore, the proposed method operates by sorted multiround roulette method without replacement. The specific process is shown in Algorithm1.Algorithm 1: Pseudo code of multiround roulette wheel selection method with replacement based on sorting.
Assuming the population size isH, the parent population H′ is generated, which is recorded as H′=a1,a2,…,ai,…,aH, where the value of fai represents the fitness of each individual. The initial state of offspring population is X=.Begin(1)
Calculate the fitness value of all individuals and arrange them in order of size. The ranked populationH^′ is denoted as H^′=b1,b2,…,bi,…,bH, where fbi−1>fbi>fbi+1. At the same time, calculate the individual h with the greatest fitness in the parent population, namely, fah=maxfa1,fa2,…,faH(2)
Calculate the total fitness of all chromosomes inH^′, ∑Hfbii=1,2,…,H(3)
Calculate the probability of single chromosome being selected,Pbi=fbi/∑i=1Hfbi(4)
Calculate the cumulative probability of each chromosome being selectedP˜K, P˜K=∑i=1hPi, respectively(5)
Rotate the wheel to generateH uniformly distributed random numbers r between [0,1](6)
Ifr≤P˜1then(7)
Select chromosomeb1;(8)
Otherwise select the h chromosome bh2≤h≤H so that P˜h−1≤r≤P˜h;(9)
Count theζ value of each interval, ζ1,ζ2,…,ζH, ζi is the number of random numbers representing the i area, i=1,2,…,H;(10)
Take the largestζ value and individual ζj=maxζ1,ζ2,…,ζH corresponding to the interval j where bj is located is the individual Hj selected for this round, there is; Hj=bjotherwiseHj=bmini1,i2,…,ijthere are the sameζvalue interval(11)
IncorporateHj into X, that is, X0=∅Xt=Xt−1∪Hj, to get X=c1,c2,…,cH;(12)
Calculate the individual with the smallest fitness value inX, fci=minfc1,fc2,…,fcH and substitute parent individual h for offspring individual j;(13)
Record all newly selected individuals.EndThe proposed method uses a sort-based multiround roulette method without replacement, which greatly increases the number of random numbers generated. It can prevent better chromosomes from being eliminated and effectively reduce selection errors. And, the higher the fitness of individual, the higher the order of ranking, which increases the probability of fitness individual being selected. Although the traditional selection operation can reduce the difference in individual selection to a certain extent, this does not ensure that the best individual in the current generation of population can enter the next step. Therefore, in the improved operator selection method, after each operation, the solution in child with the worst fitness will be replaced. This allows high-quality solution to enter the next step, so that the convergence efficiency is improved.(5) Crossover Operation. Compared with the traditional crossover method, the partial matching crossover method does not directly exchange the crossover segments of chromosomes, but the first adds the gene segments that need to be crossed before the first gene of the opposite chromosome. Then, one by one and the original individual based on comparison, if it is the same gene, it will be retained, if it is not the same, it will be removed, so a new individual is generated. This method can not only avoid the phenomenon of premature maturity but also does not require the diversity of population. The schematic diagram of partial matching crossover operation is shown in Figure 3.Figure 3
Partial matching crossover operation.Assuming parent A and B, the chromosome codes are “235|4891|67” and “584|3192|67,” respectively, and “|” represents the cross position of chromosomes. When performing the crossover operation, first put the gene segment “4891” of the selected crossover of parent A in front. Then, arrange the chromosome of other parent B in order from the position of the fifth locus. Then, compare the elements in parent B with gene segment “4891” one by one. In this process, the same gene is removed; otherwise, the gene remains unchanged. From this, we can get offspring A as “489143267.” Similarly, child B can be obtained as “319254867.”(6) Mutation Operation. In order to maintain the diversity of individuals in population and avoid the convergence of individuals in later stages of evolution, mutation operations can repair this. The mutation operation can expand the range of local optimization comprehensive target value, thereby enhancing the global search ability of algorithm. The proposed method adopts the mutation method of reversal operation to operate on an individual. First, select individuals in the population according to a certain probability and randomly generate two mutation points in the individual’s genes, resulting in mutation. Then, perform reversal operation on variant segment. The mutation operation process is shown in Figure 4.Figure 4
Mutation operation.Randomly select a certain body in populationA = “156827394” and then randomly generate two mutations “3” and “8,” namely, A = 156|8273|94. Finally, the variant segment was reversed to obtain a new individual A “156372894.”(7) Catastrophe Operation. When there are no better individuals for N consecutive generations or the population diversity change is very low, the algorithm is very likely to fall into a local optimal solution. By using the catastrophic method to kill the optimal individual of population in the next iteration, it will jump out of the local extreme value.The catastrophe judgment process is as follows: letϕ0 be the initial setting value, and the catastrophe counter starts counting from 0. When the catastrophe counter value exceeds ϕ0, it is considered that local search has produced the optimal solution, and the catastrophe is carried out. If the new optimal solution continues to be produced, the catastrophe count value is cleared. In order to prevent insufficient local search during catastrophe, the set value is enlarged and catastrophe coefficient is introduced, and then the catastrophe judgment value ϕ is as follows.Whenϕ is less than the initial setting value, count again. When ϕ is greater than the initial set value, it indicates that the population diversity change is very low. Use this value to update the catastrophe judgment value according to equation (18) and then start the countdown.(8) Termination Conditions. The termination criterion is used to judge whether the genetic operation is over. The termination conditions of algorithm generally have the following forms: the algorithm is executed to the maximum number of iterations, when the adjacent generations are unable to evolve, and the absolute value of the difference between the two target values is less than a precision. The proposed method adopts setting the maximum number of iterations N, when the evolutionary algebra reaches the maximum number of iterations, the algorithm ends. Otherwise, continue to repeat the above process.


## 3.2.1. Algorithm Selection

In the established green logistics location-routing optimization model, there is a phenomenon that cannot be compared due to the conflict between targets and different dimensions. The optimization of one objective may lead to a poor solution for another objective [26]. Since the 1990s, many scholars have begun to focus on researching intelligent optimization algorithms and introducing them into the routing problem of distribution vehicles. Intelligent optimization algorithms can be divided into two categories based on the individual behavior and based on group intelligence according to the group size. The first category includes simulated annealing algorithm and tabu search algorithm. The second category includes GA, ant colony algorithm, and particle swarm algorithm. Algorithms based on individual behavior have memory characteristics, which is deterministic, but also random. Such characteristics are often sporadic, random, and unpredictable. The algorithm based on swarm intelligence is to use the method of simulating learning mechanism and information exchange and sharing between group creatures and gradually select high-quality solutions by iteration. Among them, GA can handle constraints well, obtain global optimal solution, and has strong global search ability. But its convergence speed is slow, the local search ability is weak, the running time is longer, and it is easily affected by parameters [27].Therefore, the adaptive operator and cataclysm operator are introduced into GA. Among them, the catastrophe adaptive GA can adjust the probability of crossover and mutation in GA according to the needs, which reduces the shortcomings of standard GA that is susceptible to parameters, and also reduces the running time [28, 29]. The flow of adaptive catastrophe algorithm is shown in Figure 2.Figure 2
Adaptive catastrophe algorithm flow.


## 3.2.2. Algorithm Design

GA has the advantages of fast calculation speed, strong optimization ability, and good robustness. However, there are defects in the solution process, such as premature maturity and local trapping. Based on the characteristics of green logistics location-distribution route studied, improvements are made on the basis of traditional GA. By improving the coding, selection, crossover, and mutation methods of genetic operators, the generation of errors is reduced and the disadvantage of premature convergence is improved.(1) Coding Method. The primary problem of using GA is to solve chromosome coding, and an important basis of it is to choose the appropriate way of expression of solution. At present, many scholars have launched research on this. According to the actual situation, different encoding methods are summarized, which can be summarized as binary encoding, real number encoding, integer permutation encoding, and data structure encoding. The most used encoding method is binary encoding. 1 and 0 constitute its symbol set, which, respectively, represents the two states of on and off. This makes subsequent operations easier and easier. However, when solving sequence-related problems such as path planning, the binary encoding method has poor operability. It is very likely that infeasible solutions will often appear, which will bring greater complexity to the future. Thus, in order to improve the genetic efficiency and avoid generating infeasible solutions, the intuitive natural number and integer encoding methods for the path planning problem can well express the order of solutions obtained.In summary, the proposed method is more appropriate to use natural number coding. That is, each individuali is coded as a full array of natural numbers from 1 to n, and each natural number corresponds to the number of customer nodes that need to be served in the vehicle path. Its expression sequence means the actual driving sequence of each logistics vehicle service customer point. Because the starting point and return point of vehicles’ entire mission are distribution centers, that is, the logistics vehicles depart from distribution centers and eventually return to distribution centers after performing tasks at each customer point along the way. Based on this coding method, it can be ensured that each client node will be accessed in a certain order without duplication. It can make the processing of model constraints more convenient and quicker and facilitate the reading and modification of program. If there are n customer demand points, there are k vehicles to provide services. In this way, the chromosome is composed of distribution centers and customer point together, and they are arranged in a row in the chromosome. And, each node corresponds to a locus, and finally a n + k + 1 chromosome is formed. The code of chromosome can be expressed as(16)0,i11,i12,i1s;0,i2s,…,i2r;0,ik1,…,ikn,0,where ikm represents the n customer served by vehicle k. 0 represents a distribution center, and each 0 represents a vehicle distribution path. That is to say, the first vehicle departs from distribution center 0, goes through customers i11, i12, and i1s, and then returns to distribution center 0. The same goes for vehicle k, and the end of coding means that delivery is complete. For example, chromosome 0137029805670, which means that 3 vehicles serve 9 customers, and the corresponding path is as follows:Subpath 1: distribution center—customer 1—customer 3—customer 7—distribution centerSubpath 2: distribution center—customer 2—customer 9—customer 8—distribution centerSubpath 3: distribution center—customer 5—customer 6—customer 7—distribution center(2) Generate Initial Population. GA is different from other algorithms in that it searches for populations. Therefore, the generation of initial population is very important in the entire solution process. The pros and cons of initial population have a greater impact on the operational efficiency and computational complexity of algorithm during the calculation process. The key to determining the quality of population is the generation method of initial population and population size. If the population size is small, the algorithm’s solution process will be affected to a certain extent, and generally speaking, it will not be too good. However, choosing a larger group can avoid the situation where solution obtained is a local optimum to a certain extent. But this time will increase the space complexity, and the calculation amount and calculation time will also increase. Generally, for the case where the chromosome length is not very large, the appropriate value range for the population size is 20–150. But for large-scale solving problems, the size of population should be increased to a certain extent based on actual situation.The chromosome chooses real-valued natural number coding, and the expression of solution is related to the order of arrangement. Thus, when initializing the population, a full array ofPopsize client nodes can be generated to form a two-dimensional matrix of n×Popsize, where n represents the number of client nodes and Popsize represents the number of populations.(3) Determine the Fitness Evaluation Function. Fitness is not only a criterion for judging the quality of an individual but also a key indicator for subsequent operations. The greater the adaptability of solution is, the better it is. Therefore, the possibility of being selected will be greater. Because the proposed model needs to meet the goal of minimizing total distribution cost on the one hand, on the other hand, it is also necessary to consider time window, load, and other constraints of vehicles at any customer node. So, the feasibility of solution in the population is very important. The number of vehicles in route optimization problem is undetermined, and the number of delivery vehicles can be efficiently configured in combination with realistic conditions and related restrictions. Thus, the feasibility of solution can be determined by setting an appropriate fitness function value. At the same time, the feasible solution of a single individual can be calculated, which also represents the driving route and the number of vehicles arranged for logistics vehicle. The specific process is as follows:Define a setV to store the set of vehicles in path and use “0” element to separate different vehicles. Its purpose is to show the distribution of vehicle paths by collection containing “0” elements. This means that there is only one “0” element when the collection is initialized. Then, for single individual i, starting from the first gene that has never participated in accumulation, the delivery volume, pickup volume, and time of customer point corresponding to each gene are accumulated one by one. When the accumulation to a certain customer point does not meet constraint conditions, record the genes of cycle in set V, and add a “0” element at the end to represent the end of the cycle. If there are genes in the individual that are not involved in accumulation, they will be assigned another vehicle for the next cycle operation. Otherwise, the loop calculation ends.Therefore, the fitness function used is expressed as follows:(17)fitnessi=1∑Ci,i=1,2,…,Popsize,where ∑Ci is the comprehensive total cost of individual i, that is, the sum of fixed cost, transportation cost, carbon emission cost, and time penalty cost. The total cost changes with the distance traveled by vehicles, so it is necessary to consider whether a certain customer point reached is the starting point of the vehicles’ route. If so, you need to calculate the distance from distribution centers and the distance between this point and next point. If it is not, you only need to calculate the distance between the node and node behind it. The vehicle departs from the starting point (distribution center) and serves customers on complete path along the way. After serving the last customer, it needs to return to distribution centers. Thus, the total distance of the complete path is calculated by adding the distance from the last customer point to distribution centers to the sum of distances during cycle before the end of each cycle.(4) Selection Process. The selection operation is evolutionary process of GA. It is to select individuals with large fitness values from population to perform crossover, mutation operations, or directly as the next generation of new individuals to form a new population. This process embodies the population evolution law of “natural selection, survival of the fittest.” The selection mechanism adopted by the proposed method is a mixed selection mechanism. That is, first sort the fitness of individuals in population by size. And, according to a certain proportion, the population with higher adaptability is selected to enter the next stage to prevent chromosomes with better adaptability from being mistakenly eliminated. Then, the remaining individuals are selected in a form similar to that of betting market selection. Throughout the process, the probability of an individual being selected will increase with the increase in fitness.The traditional roulette method is simple and easy to operate. However, the result of random operation is that the selection error is relatively large, and sometimes even degradation occurs. This means that there may be individuals with higher fitness levels that are not selected but eliminated, which results in the algorithm’s difficulty in obtaining high-quality solutions and large fluctuations in the subsequent process. Therefore, the proposed method operates by sorted multiround roulette method without replacement. The specific process is shown in Algorithm1.Algorithm 1: Pseudo code of multiround roulette wheel selection method with replacement based on sorting.
Assuming the population size isH, the parent population H′ is generated, which is recorded as H′=a1,a2,…,ai,…,aH, where the value of fai represents the fitness of each individual. The initial state of offspring population is X=.Begin(1)
Calculate the fitness value of all individuals and arrange them in order of size. The ranked populationH^′ is denoted as H^′=b1,b2,…,bi,…,bH, where fbi−1>fbi>fbi+1. At the same time, calculate the individual h with the greatest fitness in the parent population, namely, fah=maxfa1,fa2,…,faH(2)
Calculate the total fitness of all chromosomes inH^′, ∑Hfbii=1,2,…,H(3)
Calculate the probability of single chromosome being selected,Pbi=fbi/∑i=1Hfbi(4)
Calculate the cumulative probability of each chromosome being selectedP˜K, P˜K=∑i=1hPi, respectively(5)
Rotate the wheel to generateH uniformly distributed random numbers r between [0,1](6)
Ifr≤P˜1then(7)
Select chromosomeb1;(8)
Otherwise select the h chromosome bh2≤h≤H so that P˜h−1≤r≤P˜h;(9)
Count theζ value of each interval, ζ1,ζ2,…,ζH, ζi is the number of random numbers representing the i area, i=1,2,…,H;(10)
Take the largestζ value and individual ζj=maxζ1,ζ2,…,ζH corresponding to the interval j where bj is located is the individual Hj selected for this round, there is; Hj=bjotherwiseHj=bmini1,i2,…,ijthere are the sameζvalue interval(11)
IncorporateHj into X, that is, X0=∅Xt=Xt−1∪Hj, to get X=c1,c2,…,cH;(12)
Calculate the individual with the smallest fitness value inX, fci=minfc1,fc2,…,fcH and substitute parent individual h for offspring individual j;(13)
Record all newly selected individuals.EndThe proposed method uses a sort-based multiround roulette method without replacement, which greatly increases the number of random numbers generated. It can prevent better chromosomes from being eliminated and effectively reduce selection errors. And, the higher the fitness of individual, the higher the order of ranking, which increases the probability of fitness individual being selected. Although the traditional selection operation can reduce the difference in individual selection to a certain extent, this does not ensure that the best individual in the current generation of population can enter the next step. Therefore, in the improved operator selection method, after each operation, the solution in child with the worst fitness will be replaced. This allows high-quality solution to enter the next step, so that the convergence efficiency is improved.(5) Crossover Operation. Compared with the traditional crossover method, the partial matching crossover method does not directly exchange the crossover segments of chromosomes, but the first adds the gene segments that need to be crossed before the first gene of the opposite chromosome. Then, one by one and the original individual based on comparison, if it is the same gene, it will be retained, if it is not the same, it will be removed, so a new individual is generated. This method can not only avoid the phenomenon of premature maturity but also does not require the diversity of population. The schematic diagram of partial matching crossover operation is shown in Figure 3.Figure 3
Partial matching crossover operation.Assuming parent A and B, the chromosome codes are “235|4891|67” and “584|3192|67,” respectively, and “|” represents the cross position of chromosomes. When performing the crossover operation, first put the gene segment “4891” of the selected crossover of parent A in front. Then, arrange the chromosome of other parent B in order from the position of the fifth locus. Then, compare the elements in parent B with gene segment “4891” one by one. In this process, the same gene is removed; otherwise, the gene remains unchanged. From this, we can get offspring A as “489143267.” Similarly, child B can be obtained as “319254867.”(6) Mutation Operation. In order to maintain the diversity of individuals in population and avoid the convergence of individuals in later stages of evolution, mutation operations can repair this. The mutation operation can expand the range of local optimization comprehensive target value, thereby enhancing the global search ability of algorithm. The proposed method adopts the mutation method of reversal operation to operate on an individual. First, select individuals in the population according to a certain probability and randomly generate two mutation points in the individual’s genes, resulting in mutation. Then, perform reversal operation on variant segment. The mutation operation process is shown in Figure 4.Figure 4
Mutation operation.Randomly select a certain body in populationA = “156827394” and then randomly generate two mutations “3” and “8,” namely, A = 156|8273|94. Finally, the variant segment was reversed to obtain a new individual A “156372894.”(7) Catastrophe Operation. When there are no better individuals for N consecutive generations or the population diversity change is very low, the algorithm is very likely to fall into a local optimal solution. By using the catastrophic method to kill the optimal individual of population in the next iteration, it will jump out of the local extreme value.The catastrophe judgment process is as follows: letϕ0 be the initial setting value, and the catastrophe counter starts counting from 0. When the catastrophe counter value exceeds ϕ0, it is considered that local search has produced the optimal solution, and the catastrophe is carried out. If the new optimal solution continues to be produced, the catastrophe count value is cleared. In order to prevent insufficient local search during catastrophe, the set value is enlarged and catastrophe coefficient is introduced, and then the catastrophe judgment value ϕ is as follows.Whenϕ is less than the initial setting value, count again. When ϕ is greater than the initial set value, it indicates that the population diversity change is very low. Use this value to update the catastrophe judgment value according to equation (18) and then start the countdown.(8) Termination Conditions. The termination criterion is used to judge whether the genetic operation is over. The termination conditions of algorithm generally have the following forms: the algorithm is executed to the maximum number of iterations, when the adjacent generations are unable to evolve, and the absolute value of the difference between the two target values is less than a precision. The proposed method adopts setting the maximum number of iterations N, when the evolutionary algebra reaches the maximum number of iterations, the algorithm ends. Otherwise, continue to repeat the above process.


## 4. Experiment and Analysis

### 4.1. Simulation Data

In the experiment, the data comes from actual operating data of a company (company U). U company’s suppliers’ hand over order information to U company, and then, by U company’s internal operations, the transportation department is responsible for supermarket distribution. In the entire process after receiving the goods and before delivery, U company has control over the goods and can guarantee the full temperature control. Among them, the cold chain logistics division of U company is responsible for fresh food distribution business of 20 supermarkets of different sizes. Its current distribution status is shown in Table1.Table 1
Distribution situation of U company.
ProjectDetails of the companyTypes of foodFresh agricultural productsNumber of service stores20 large, medium and small supermarketsNumber of vehicles12Vehicle typeH-type refrigerator carIn order to facilitate business, company U plans to select several suitable cold storages as distribution centers among the existing cold storages of the group. Through the investigation of U company's cold chain logistics business, and according to its business data for the first half of 2020, the geographic location and average delivery volume of its 20 stores in charge of delivery were obtained. For case analysis, it is necessary to preprocess the collected data such as the location of alternative distribution centers and stores and the demand of each store. The main preprocessing step is to first convert latitude and longitude information of the 3 alternative distribution centers and 20 stores into the format of degrees, minutes, and seconds. Then, we use the MAPGIS software to convert latitude and longitude into abscissa and ordinate in the rectangular coordinate system and then calculate the distance between points.First, sort out the latitude and longitude information of 20 stores and 3 alternative distribution centers and number them. Among them are 1–3 four alternative distribution centers. 4–23 represents 20 stores, and the specific latitude and longitude information is shown in Table2.Table 2
Longitude and latitude coordinates of alternative distribution centers and stores.
Serial numberLongitudeLatitude1118.27938.1252118.35037.9683118.31938.0854118.41737.8425118.43737.9466118.34438.0787118.43837.8728118.46537.9349118.45937.92210118.46737.90011118.16737.76312118.52937.94713118.60937.90614118.33038.08215118.60037.91716118.58137.92917118.26737.92418118.26537.89619118.29237.87820118.29137.93021118.43337.85222118.28337.93523118.35138.045In order to more intuitively observe the geographical distribution of each store and alternative distribution center, a 23-point geographical distribution scatter plot is made, as shown in Figure5.Figure 5
Scatter plot of geographical distribution.In addition, green logistics and distribution need to meet the customer’s required time window, maximum tolerable time window, and service time. The data of 15 users are shown in Table3.Table 3
Time window, maximum tolerance time window and service time of demand points.
Distribution pointDemand time windowMaximum tolerance time windowService time (h)Requirement (t)15 : 35–17 : 105 : 00–17 : 300—26 : 00–8 : 455 : 30–9 : 100.251.637 : 20–9 : 007 : 00–9 : 300.342.546 : 00–7 : 305 : 40–8 : 000.121.356 : 20–8 : 306 : 00–9 : 000.201.666 : 50–9 : 106 : 30–9 : 300.151.977 : 10–9 : 006 : 40–9 : 200.321.187 : 30–9 : 107 : 00–9 : 400.290.597 : 30–9 : 307 : 20–10 : 000.062.4107 : 35–9 : 257 : 00–9 : 500.062.2117 : 10–9 : 007 : 00–9 : 250.111.5127 : 15–9 : 306 : 55–10 : 000.161.7137 : 20–9 : 207 : 00–9 : 450.250.8147 : 30–9 : 307 : 00–10 : 000.361.4156 : 45–9 : 006 : 30–9 : 300.082.1


### 4.2. Parameter Setting

Except for the distribution point information and alternative distribution center information, other parameter settings are shown in Tables4 and 5 . As of 2020, the historically highest carbon trading price is 83 yuan/ton, the lowest is 31.5 yuan/ton, and the historical average price is 47 yuan/ton. So, take the average price here.Table 4
Parameter setting.
ParameterValueAlternative distribution center3Supermarket20Deterioration rate during transportation0.0001Deterioration rate of door opening0.0003Service time, h0.10Average speed, km/h45Average price, yuan/kg9Initial carbon quota, kg/day560Carbon trading price, yuan/kg0.05Vehicle load, t2Number of vehicles12Table 5
Vehicle parameter setting.
ParameterValueOutside vehicle size, m5.9 × 2.2 × 3.1Deterioration coefficient0.07Heat transfer area, m240Carriage size, m4 × 2.1 × 1.8Heat transfer coefficient,kcal/h∗m2∗C°0.12Temperature difference, °C8In order to calculate the total cost of whole process, the following parameters are also needed: the electricity consumed in storage process is regarded as industrial electricity, and the electricity price is 1 yuan/kWh; the fuel consumption during idling process is 8 (liters/100-ton kilometers). The fuel consumption is 12 (liters/hundred-ton kilometers); the emission factor is the amount of carbon dioxide emissions per unit of fuel consumed. The refrigerated vehicle used is a diesel vehicle. According to data, 1 liter of diesel produces 2.65 kg of carbon dioxide, which means an emission factor of 2.65. The fixed usage cost of vehicles is 70 yuan/day; the variable cost of vehicles is 25 yuan/km.


### 4.3. Convergence Curve

The proposed method uses improved GA to solve the optimization problem. In order to demonstrate its performance, compared with standard GA. Figure6 shows the optimization iteration curves of the two GA total objective functions considering the carbon emission model.Figure 6
Optimization iteration curve of total objective function.It can be seen from Figure6 that improved GA achieves convergence after about 500 iterations, and the fitness value is 0.11. The standard GA iteration takes about 700 times to achieve convergence, and fitness value is 0.125 when it is stable. It can be seen that improved GA converges faster, and the average fitness value and best fitness value are both smaller. Thus, its optimization performance is more ideal, and it is more suitable for handling the selection and route optimization problems of distribution centers in green logistics.


### 4.4. Parameter Analysis

In terms of population size, the initial population size is usually set to 30–150. If the population is small, it will limit the search space of algorithm and easily fall into the local optimum. If population size is too high, the search will be more complicated and search efficiency will be reduced. This results in an increase in the number of iterations and a long running time, which reduces the algorithm efficiency. Therefore, in order to select a suitable population size, the crossover rate is 0.8, the mutation rate is 0.1, and the number of iterations is 1000 generations to simulate the situation of different population sizes. Observe its influence on the value of objective function, the number of iterations for the first optimal solution, and the running time of program. The change of objective function value under different population sizes is shown in Figure7.Figure 7
Objective function values under different population sizes.It can be seen from Figure7 that when population size N increases from 30 to 90, the overall change of objective function value is relatively large and fluctuates significantly due to the relatively small search space range. The overall objective function value is too large. When population size N is 120 and 150, there is little difference in the value of objective value function. The value is smaller than population size below 90, and fluctuation is relatively stable. Therefore, if only judging from the operating results of target value, it can be inferred that the population size selection 120 and 150 can make the objective function optimal.Similarly, the number of iterations to obtain the optimal solution for the first time under different population sizes is shown in Figure8.Figure 8
The number of iterations to obtain the optimal solution for the first time under different population sizes.It can be seen from Figure8 that when population size is 30 and 60, the algebra that produces the optimal solution is very unstable, and fluctuation range is also very large. And, the algebra is relatively late, basically above 1000. When the population size is 120 and 150, the number of iterations to produce the optimal solution is relatively stable, basically maintaining a slight fluctuation in the 800 to 1000 generations. Besides, its convergence is also acceptable. When the population size is 90 generations, although the overall number of iterations is relatively concentrated between 900 and 1000, the uncertainty is also very large, and the maximum number of iterations and minimum number of iterations are quite different. Compared with the population size of 120 and 150, the instability of its existence is more prominent. Therefore, it is a good choice to choose a population size of 120 or 150. The specific selection needs to be further analyzed in combination with other parameters.The change of running time under different population sizes is shown in Figure9.Figure 9
Running time under different population sizes.It can be seen from Figure9 that no matter what the value of population size N is, the running time is quite stable. However, as the population size increases, the number of selections, crossover, and mutation operations will increase, and the amount of operations on memory data will increase proportionally. When the population size increases from 30 to 150, the running time also gradually increases, and the increase is getting bigger and bigger. When the population is 150, it greatly exceeds 2 times the running time for population size of 60. Thus, when population size is within 120, the increase in operating time is relatively stable, which is acceptable.Based on the experimental results in the above three cases, when the population size is 60 and 90, the performance of target value and the first optimal solution is poor, and the stability is not good. When population size is 150, the algebraic fluctuation of optimal solution is slightly larger due to the long running time. Therefore, the effect is not ideal. When population size is 120, objective function value in solving objective function value and the algebra that produces optimal solution is better, and the running speed is appropriate. Therefore, it is more appropriate to choose 120 as population size to solve the problem.Since different crossover rates and mutation rates will have a certain impact on the objective function value, when selecting the optimal crossover rate and optimal mutation rate, more attention is paid to the change of objective function value under different values. Figure10 shows the result of objective function value under different crossover rates and mutation rates. The population size selected during the experiment is 120, and the number of termination iterations is 1000.Figure 10
Objective function values under different crossover rates and mutation rates.It can be seen from Figure10 that, for crossover rate, when crossover rate is in the interval of [0.35, 0.55], the objective function value shows a downward trend. In the interval of [0.55, 0.75], the function value showed an upward trend, and in the interval of [0.75, 0.85], the function value showed a downward trend. Therefore, the minimum value is reached when the crossover ratio is 0.55 and 0.85. By comparing the function values corresponding to the two extreme points, it is found that the objective function value is smaller when the crossover rate is 0.85. Therefore, crossover probability is 0.85, which is the easiest to find the optimal solution.For the case of mutation rate, set the mutation rate in the interval of [0.01,0.1] for calculation analysis. In the interval [0.01,0.07], the objective function value generally showed a downward trend, while in the interval [0.07,0.10], the objective function value showed an upward trend, and the average value of objective function reached the minimum when variation rate was 0.07. Thus, it is determined that the mutation probability is set at 0.07, and the objective function value is the smallest at this time. Therefore, after a comprehensive analysis of multiple running data, crossover rate and mutation rate are finally set to 0.85 and 0.07.


### 4.5. Path Planning Results

From the perspective of low-carbon and environmental protection, the path planning results obtained by proposed method are shown in Figure11. Among them, one of the three candidate locations is selected as a distribution center to distribute cold chain agricultural and sideline products to 15 demand points.Figure 11
Generation diagram of distribution path.From Figure11, the path of this delivery task can be given:(1)
Path 1: 0⟶14⟶15⟶9⟶1⟶11⟶0(2)
Path 2: 0⟶2⟶3⟶10⟶0(3)
Path 3: 0⟶13⟶5⟶7⟶6⟶4⟶12⟶0(4)
Path 4: 0⟶8⟶6⟶0According to improved GA iterative optimization, the optimal distribution route division and service order of customer demand points within the route are generated. The delivery is carried out by four vehicles, respectively. The driving distance of route 1 is 289.25 km, the driving distance of route 2 is 162.37 km, the driving distance of route 3 is 301.81, and the driving distance of Route 4 is 98.56 km. The total distance traveled for this distribution task is 851.99 km. From this, it can be estimated that the cost of distribution includes the transportation cost of 5642.75 yuan and the departure cost of 800 yuan. The penalty cost of customer satisfaction caused by the deviation from time window is 328.55 yuan, and the total distribution cost is 6771.3 yuan. It can be proved that the logistics location-routing optimization problem model using improved GA algorithm is effective.


## 4.1. Simulation Data

In the experiment, the data comes from actual operating data of a company (company U). U company’s suppliers’ hand over order information to U company, and then, by U company’s internal operations, the transportation department is responsible for supermarket distribution. In the entire process after receiving the goods and before delivery, U company has control over the goods and can guarantee the full temperature control. Among them, the cold chain logistics division of U company is responsible for fresh food distribution business of 20 supermarkets of different sizes. Its current distribution status is shown in Table1.Table 1
Distribution situation of U company.
ProjectDetails of the companyTypes of foodFresh agricultural productsNumber of service stores20 large, medium and small supermarketsNumber of vehicles12Vehicle typeH-type refrigerator carIn order to facilitate business, company U plans to select several suitable cold storages as distribution centers among the existing cold storages of the group. Through the investigation of U company's cold chain logistics business, and according to its business data for the first half of 2020, the geographic location and average delivery volume of its 20 stores in charge of delivery were obtained. For case analysis, it is necessary to preprocess the collected data such as the location of alternative distribution centers and stores and the demand of each store. The main preprocessing step is to first convert latitude and longitude information of the 3 alternative distribution centers and 20 stores into the format of degrees, minutes, and seconds. Then, we use the MAPGIS software to convert latitude and longitude into abscissa and ordinate in the rectangular coordinate system and then calculate the distance between points.First, sort out the latitude and longitude information of 20 stores and 3 alternative distribution centers and number them. Among them are 1–3 four alternative distribution centers. 4–23 represents 20 stores, and the specific latitude and longitude information is shown in Table2.Table 2
Longitude and latitude coordinates of alternative distribution centers and stores.
Serial numberLongitudeLatitude1118.27938.1252118.35037.9683118.31938.0854118.41737.8425118.43737.9466118.34438.0787118.43837.8728118.46537.9349118.45937.92210118.46737.90011118.16737.76312118.52937.94713118.60937.90614118.33038.08215118.60037.91716118.58137.92917118.26737.92418118.26537.89619118.29237.87820118.29137.93021118.43337.85222118.28337.93523118.35138.045In order to more intuitively observe the geographical distribution of each store and alternative distribution center, a 23-point geographical distribution scatter plot is made, as shown in Figure5.Figure 5
Scatter plot of geographical distribution.In addition, green logistics and distribution need to meet the customer’s required time window, maximum tolerable time window, and service time. The data of 15 users are shown in Table3.Table 3
Time window, maximum tolerance time window and service time of demand points.
Distribution pointDemand time windowMaximum tolerance time windowService time (h)Requirement (t)15 : 35–17 : 105 : 00–17 : 300—26 : 00–8 : 455 : 30–9 : 100.251.637 : 20–9 : 007 : 00–9 : 300.342.546 : 00–7 : 305 : 40–8 : 000.121.356 : 20–8 : 306 : 00–9 : 000.201.666 : 50–9 : 106 : 30–9 : 300.151.977 : 10–9 : 006 : 40–9 : 200.321.187 : 30–9 : 107 : 00–9 : 400.290.597 : 30–9 : 307 : 20–10 : 000.062.4107 : 35–9 : 257 : 00–9 : 500.062.2117 : 10–9 : 007 : 00–9 : 250.111.5127 : 15–9 : 306 : 55–10 : 000.161.7137 : 20–9 : 207 : 00–9 : 450.250.8147 : 30–9 : 307 : 00–10 : 000.361.4156 : 45–9 : 006 : 30–9 : 300.082.1


## 4.2. Parameter Setting

Except for the distribution point information and alternative distribution center information, other parameter settings are shown in Tables4 and 5 . As of 2020, the historically highest carbon trading price is 83 yuan/ton, the lowest is 31.5 yuan/ton, and the historical average price is 47 yuan/ton. So, take the average price here.Table 4
Parameter setting.
ParameterValueAlternative distribution center3Supermarket20Deterioration rate during transportation0.0001Deterioration rate of door opening0.0003Service time, h0.10Average speed, km/h45Average price, yuan/kg9Initial carbon quota, kg/day560Carbon trading price, yuan/kg0.05Vehicle load, t2Number of vehicles12Table 5
Vehicle parameter setting.
ParameterValueOutside vehicle size, m5.9 × 2.2 × 3.1Deterioration coefficient0.07Heat transfer area, m240Carriage size, m4 × 2.1 × 1.8Heat transfer coefficient,kcal/h∗m2∗C°0.12Temperature difference, °C8In order to calculate the total cost of whole process, the following parameters are also needed: the electricity consumed in storage process is regarded as industrial electricity, and the electricity price is 1 yuan/kWh; the fuel consumption during idling process is 8 (liters/100-ton kilometers). The fuel consumption is 12 (liters/hundred-ton kilometers); the emission factor is the amount of carbon dioxide emissions per unit of fuel consumed. The refrigerated vehicle used is a diesel vehicle. According to data, 1 liter of diesel produces 2.65 kg of carbon dioxide, which means an emission factor of 2.65. The fixed usage cost of vehicles is 70 yuan/day; the variable cost of vehicles is 25 yuan/km.


## 4.3. Convergence Curve

The proposed method uses improved GA to solve the optimization problem. In order to demonstrate its performance, compared with standard GA. Figure6 shows the optimization iteration curves of the two GA total objective functions considering the carbon emission model.Figure 6
Optimization iteration curve of total objective function.It can be seen from Figure6 that improved GA achieves convergence after about 500 iterations, and the fitness value is 0.11. The standard GA iteration takes about 700 times to achieve convergence, and fitness value is 0.125 when it is stable. It can be seen that improved GA converges faster, and the average fitness value and best fitness value are both smaller. Thus, its optimization performance is more ideal, and it is more suitable for handling the selection and route optimization problems of distribution centers in green logistics.


## 4.4. Parameter Analysis

In terms of population size, the initial population size is usually set to 30–150. If the population is small, it will limit the search space of algorithm and easily fall into the local optimum. If population size is too high, the search will be more complicated and search efficiency will be reduced. This results in an increase in the number of iterations and a long running time, which reduces the algorithm efficiency. Therefore, in order to select a suitable population size, the crossover rate is 0.8, the mutation rate is 0.1, and the number of iterations is 1000 generations to simulate the situation of different population sizes. Observe its influence on the value of objective function, the number of iterations for the first optimal solution, and the running time of program. The change of objective function value under different population sizes is shown in Figure7.Figure 7
Objective function values under different population sizes.It can be seen from Figure7 that when population size N increases from 30 to 90, the overall change of objective function value is relatively large and fluctuates significantly due to the relatively small search space range. The overall objective function value is too large. When population size N is 120 and 150, there is little difference in the value of objective value function. The value is smaller than population size below 90, and fluctuation is relatively stable. Therefore, if only judging from the operating results of target value, it can be inferred that the population size selection 120 and 150 can make the objective function optimal.Similarly, the number of iterations to obtain the optimal solution for the first time under different population sizes is shown in Figure8.Figure 8
The number of iterations to obtain the optimal solution for the first time under different population sizes.It can be seen from Figure8 that when population size is 30 and 60, the algebra that produces the optimal solution is very unstable, and fluctuation range is also very large. And, the algebra is relatively late, basically above 1000. When the population size is 120 and 150, the number of iterations to produce the optimal solution is relatively stable, basically maintaining a slight fluctuation in the 800 to 1000 generations. Besides, its convergence is also acceptable. When the population size is 90 generations, although the overall number of iterations is relatively concentrated between 900 and 1000, the uncertainty is also very large, and the maximum number of iterations and minimum number of iterations are quite different. Compared with the population size of 120 and 150, the instability of its existence is more prominent. Therefore, it is a good choice to choose a population size of 120 or 150. The specific selection needs to be further analyzed in combination with other parameters.The change of running time under different population sizes is shown in Figure9.Figure 9
Running time under different population sizes.It can be seen from Figure9 that no matter what the value of population size N is, the running time is quite stable. However, as the population size increases, the number of selections, crossover, and mutation operations will increase, and the amount of operations on memory data will increase proportionally. When the population size increases from 30 to 150, the running time also gradually increases, and the increase is getting bigger and bigger. When the population is 150, it greatly exceeds 2 times the running time for population size of 60. Thus, when population size is within 120, the increase in operating time is relatively stable, which is acceptable.Based on the experimental results in the above three cases, when the population size is 60 and 90, the performance of target value and the first optimal solution is poor, and the stability is not good. When population size is 150, the algebraic fluctuation of optimal solution is slightly larger due to the long running time. Therefore, the effect is not ideal. When population size is 120, objective function value in solving objective function value and the algebra that produces optimal solution is better, and the running speed is appropriate. Therefore, it is more appropriate to choose 120 as population size to solve the problem.Since different crossover rates and mutation rates will have a certain impact on the objective function value, when selecting the optimal crossover rate and optimal mutation rate, more attention is paid to the change of objective function value under different values. Figure10 shows the result of objective function value under different crossover rates and mutation rates. The population size selected during the experiment is 120, and the number of termination iterations is 1000.Figure 10
Objective function values under different crossover rates and mutation rates.It can be seen from Figure10 that, for crossover rate, when crossover rate is in the interval of [0.35, 0.55], the objective function value shows a downward trend. In the interval of [0.55, 0.75], the function value showed an upward trend, and in the interval of [0.75, 0.85], the function value showed a downward trend. Therefore, the minimum value is reached when the crossover ratio is 0.55 and 0.85. By comparing the function values corresponding to the two extreme points, it is found that the objective function value is smaller when the crossover rate is 0.85. Therefore, crossover probability is 0.85, which is the easiest to find the optimal solution.For the case of mutation rate, set the mutation rate in the interval of [0.01,0.1] for calculation analysis. In the interval [0.01,0.07], the objective function value generally showed a downward trend, while in the interval [0.07,0.10], the objective function value showed an upward trend, and the average value of objective function reached the minimum when variation rate was 0.07. Thus, it is determined that the mutation probability is set at 0.07, and the objective function value is the smallest at this time. Therefore, after a comprehensive analysis of multiple running data, crossover rate and mutation rate are finally set to 0.85 and 0.07.


## 4.5. Path Planning Results

From the perspective of low-carbon and environmental protection, the path planning results obtained by proposed method are shown in Figure11. Among them, one of the three candidate locations is selected as a distribution center to distribute cold chain agricultural and sideline products to 15 demand points.Figure 11
Generation diagram of distribution path.From Figure11, the path of this delivery task can be given:(1)
Path 1: 0⟶14⟶15⟶9⟶1⟶11⟶0(2)
Path 2: 0⟶2⟶3⟶10⟶0(3)
Path 3: 0⟶13⟶5⟶7⟶6⟶4⟶12⟶0(4)
Path 4: 0⟶8⟶6⟶0According to improved GA iterative optimization, the optimal distribution route division and service order of customer demand points within the route are generated. The delivery is carried out by four vehicles, respectively. The driving distance of route 1 is 289.25 km, the driving distance of route 2 is 162.37 km, the driving distance of route 3 is 301.81, and the driving distance of Route 4 is 98.56 km. The total distance traveled for this distribution task is 851.99 km. From this, it can be estimated that the cost of distribution includes the transportation cost of 5642.75 yuan and the departure cost of 800 yuan. The penalty cost of customer satisfaction caused by the deviation from time window is 328.55 yuan, and the total distribution cost is 6771.3 yuan. It can be proved that the logistics location-routing optimization problem model using improved GA algorithm is effective.


## 5. Conclusion

This paper proposes a location-routing optimization solution of green logistics using improved GA algorithm from the perspective of low-carbon and environmental protection. Considering the cost factors, time window, deterioration rate, carbon trading mechanism, and other influencing factors, a green logistics location-routing optimization model is constructed, which can minimize the total cost as optimization goal. Besides, the catastrophic adaptive GA algorithm is used to solve the optimization problem, so as to obtain a low-carbon, economical, and efficient distribution route. The distribution cost after route optimization is 6771.3 yuan, which meets the economic and environmental design requirements. Such a reasonable design of logistics location-routing plan can not only enhance customer satisfaction with logistics and distribution links but also reduce the distribution costs of enterprises and save ecological environment energy.In fact, the location-routing optimization problem is more complicated. In order to simplify the model, the proposed problem model has been idealized. In the next study, multiple models, multiple products, multiple cycles, and uncertain demand will be considered, as well as real-time road conditions and other influencing factors to further meet the needs of actual scenarios.



---

*Source: 6101194-2021-11-30.xml*