Location Optimization of Fresh Agricultural Products Cold Chain Distribution Center under Carbon Emission Constraints

: This paper seeks to effectively realize energy saving and emission reduction in the process of location for fresh agricultural cold chain logistics. It is based on the traditional location for distribution center optimized for both freshness and carbon emissions. A bi-objective function location model was constructed to minimize the total cost and carbon emission and a two-stage heuristic algorithm was designed to solve the model. According to analyzing the location case of Y enterprise in Zhejiang Province, the THA had an average total cost optimization rate of 7.70% and an average carbon emission optimization rate of 10.23% compared to the WOA, while it had a rate of 12.77% and 14.12%, respectively, compared to the PSO. When the unit carbon emission cost increased from 0.02 dollar/kg to 0.06 dollar/kg, the comprehensive carbon emission cost increased by 42.07% and the total cost increased by 3.05%. Therefore, logistics enterprises can achieve the reduction of logistics costs and sustainable development through reasonable location for fresh agricultural products cold chain distribution centers.


Introduction
As people's living standards improve, the demand for the cold chain for fresh agricultural products is on the rise, and the quality requirements of such products are getting more stringent [1][2][3].Since fresh agricultural products are vulnerable and need to be delivered on time, different procedures of logistics are independent of each other in the traditional cold chain system, so it is difficult to transport these goods under a stable temperature during the whole transit [4,5].Meanwhile, the unreasonable use of devices such as cold storage, refrigerated vehicles, and automated temperature control consumes a lot of energy, causes large carbon emissions, and harms the environment [6,7].As the concept of sustainable development prevails, the domestic carbon trading market has been initially launched, and the government-led policies of energy conservation and emission reduction have been gradually transformed into a market-oriented state [8], in which market regulation is made full use of to promote the green emission reduction of enterprises.
In this context, the optimal location decision of the cold chain distribution center of fresh agricultural products can promote the reasonable planning and unified scheduling of the logistics system of the cold chain distribution center.While ensuring the quality of fresh agricultural products, which can effectively improve the utilization rate of cold chain equipment, reduce the cost of carbon emissions, and promote the sustainable development of the whole society [9][10][11].

Literature Review 2.1. Carbon Emissions from Cold Chain Logistics
With the deepening of the concept of sustainable development, scholars have conducted active research on how to reduce carbon emissions while ensuring the development of cold chain logistics.Early literature studies mainly used carbon footprint to optimize the configuration of the cold chain supply chain.For example, Tian [12] first proposed the carbon footprint of Xinjiang grape cold chain logistics, and used the single-stage supply chain profit model to conduct research and analysis; Benjaafar et al. [13] applied carbon footprint to the cold chain.In the structural analysis of the supply chain, optimization measures are proposed based on the carbon footprint of the milk supply chain.Although the research on the optimal configuration of the supply chain based on carbon footprint expands the theoretical connotation of carbon emissions in cold chain logistics, it is difficult for logistics companies to evaluate the specific carbon emissions in cold chain operations.In recent years, more and more scholars have analyzed and optimized cold chain logistics activities through carbon emission constraints.For example, Pan et al. [14] introduced carbon emission constraints into the cold chain logistics distribution process and transformed them into corresponding economic indicators.The structure considers carbon emissions and a cost-based quantitative model of cold chain logistics distribution path optimization.Li [15] considered the comprehensive cost of carbon, optimized the inventory and distribution process of the fresh supply chain warehouse, and established a two-level inventory-path optimization model with a time window.However, at present, there are relatively few studies that consider carbon emission factors in the location selection of cold chain logistics distribution centers.The location of distribution centers will directly affect vehicle routing and inventory allocation, and further, affect the dynamic carbon emissions in transportation and static carbon in warehousing.emission.

Location for Logistics Distribution Centers
The theoretical study of the location problem of logistics distribution centers (LDC) originated from Weber's theory [16] in the early 20th century and was further extended by the P-median model proposed by Hakimi [17].The location selection for logistics distribution centers is a typical NP-hard problem with a lot of complex constraints and non-smooth characteristics [18].After decades of research and development, the solution to this type of problem can be divided into qualitative and quantitative location models.Qualitative models are mainly based on the construction of target evaluation systems for performance evaluation and ranking, such as the analytic hierarchy process [19], fuzzy evaluation method [20], grey relational degree method [21], etc. Qualitative models are beneficial to deal with uncertain influencing factors, and the evaluation results are rich in information, but the judgment level is more complicated, and it is difficult to meet the requirements of standardization and consistency.Quantitative models mainly calculate the location of the logistics distribution center through quantifiable indicators and are mainly divided into accurate solution methods and heuristic algorithms.The accurate algorithms for solving LDC are mainly based on the centroid method [22], linear programming method [23], Baumol-Wolfe model [24], etc.These algorithms can solve the optimal solution of the problem, but at the cost of time and space complexity in exchange for the quality of the understanding, they cannot be used to meet the needs of actually solving large-scale problems.The application of heuristic algorithms to solve LDC originated from the first use of a heuristic algorithm by Steven et al. to solve the logistics distribution center location problem considering inventory factors [25].Heuristic algorithms mainly based on genetic algorithm [26], ant colony algorithm [27], tabu search algorithm [28], etc.A balance is achieved in the quality of time, space, and reconciliation.Considering the scale and complexity of practical problems, most scholars currently use heuristic algorithms to solve the LDC problem.

Location for Cold Chain Logistics Distribution Centers
The location problem of cold chain logistics distribution center is developed on the basis of the research of distribution center location problem.Therefore, qualitative and quantitative research methods can be used for reference.At the same time, the factors that need to be considered are more complex and diversified.In recent years, it has gradually attracted the research and discussion of scholars.In terms of qualitative models: Bai et al. [29] used the grey correlation method and AHP method to solve the location problem of cold chain logistics distribution center; Pan [30] considered various risk factors of cold chain distribution center location and established a risk evaluation model of cold chain distribution center location based on the fuzzy network analysis method.From the literature analysis, it can be seen that most of the current qualitative models lack the objective level of inter-factor influence analysis and lack certain realistic significance.In terms of quantitative analysis: Montanari [31] proposed that temperature is an important factor affecting cold chain food, and began to apply Euler algorithm and Lagrange algorithm to solve the location model of cold chain distribution centers; Demirtaş et al. [32] began to use heuristic algorithms to conduct research on the cold chain logistics distribution centers of fresh vegetables, fruits and fresh agricultural products; Li et al. [33] proposed a particle swarm algorithm with mutation and dynamic self-adaptation to solve the location model of logistics distribution center with the least total cost and meeting the timeliness.From the analysis of the literature, it can be seen that in the quantitative model with the objective of minimizing the total cost, the construction cost, distribution cost and cold storage cost are mainly considered.The dynamic and static carbon emission costs are not taken into account, which is not conducive to reducing the environmental impact of carbon emission in the process of cold chain logistics.

Research Framework and Summary of Ideas
In view of the above, we put forward the following research questions in this study:

•
Question 1: How to build the location model of fresh agricultural products cold chain distribution center on the premise of considering carbon emissions?• Question 2: How to design an optimization algorithm to solve the constructed model?• Question 3: How to prove the correctness of the constructed model and algorithm through the solution of a case?
The objective of this research work was to construct the location model and design solution algorithm of fresh agricultural products cold chain distribution center.First, based on the consideration of cold chain storage end cost, transportation end cost, freshness penalty cost and comprehensive carbon emission cost, we constructed the location model of fresh agricultural products cold chain distribution center with minimizing the total location cost and minimizing the carbon emission in the distribution process as a double objective function.Second, based on the characteristics of the model, we designed a two-stage heuristic algorithm to solve it.In the first stage, we used DBSCAN spatial clustering algorithm to cluster and determine the service range of the alternative points of the distribution center to form the initial solution set; In the second stage, the improved whale optimization algorithm was used to optimize the initial solution set to obtain the overall optimal location scheme.Finally, the model was applied to Y enterprise in Zhejiang Province to verify the correctness and applicability of the method.
After the introduction (Section 1), the organizational structure of this research work is as follows.Section 2 makes a comprehensive literature review on the location problem.Section 3 introduces the location model of cold chain distribution center of fresh agricultural products.Section 4 introduces the solution algorithm designed in this research work in detail.Section 5 introduces the research cases and data collection.Section 6 discusses the location of Y enterprise in Zhejiang Province, compares the calculation results with other methods, and analyzes the sensitivity of carbon cost.Finally, Section 7 summarizes the conclusions and determines the future research direction.The overall research framework of this paper is shown in Figure 1.
the conclusions and determines the future research direction.The overall research framework of this paper is shown in Figure 1.

Descriptions and Hypotheses of Location Selection
In general, the location model of cold chain distribution centers of fresh agricultural products under carbon emission constraints can be described as follows.There are several demand points for fresh agricultural products in a certain region with fixed demands, from which a number of alternative locations are selected.Some of the alternative locations are selected for distribution centers, which can provide integrated services and meet each region's requirements of the quantity of fresh agricultural products and carbon emission, with the overall costs of location selection and carbon emissions of warehousing and distribution minimized [34][35][36][37][38].The hypotheses in the model are shown in Table 1.

Supply and Demand
(1) the inventory of a single distribution center can satisfy the demand of covering demand points, and the inventory is equal to the sum of the demand of all demand points (2) the demand for fresh agricultural products at each demand point is fixed and will not fluctuate greatly in the short term Transportation (1) alternative locations are equipped with sufficient refrigerated vehicles for transport, regardless of shipping batches and fixed costs of vehicles (2) refrigerated vehicles travel at a constant speed, regardless of service time window constraints and vehicles routing Carbon Emissions (1) the site selection process involves static carbon emissions from distribution centers and dynamic carbon emissions from refrigerated vehicle transportation (2) carbon emissions from product decay are not considered Freshness (1) the freshness of fresh agricultural products is only related to the transportation time (2) the effect of temperature change on freshness is not considered

Descriptions and Hypotheses of Location Selection
In general, the location model of cold chain distribution centers of fresh agricultural products under carbon emission constraints can be described as follows.There are several demand points for fresh agricultural products in a certain region with fixed demands, from which a number of alternative locations are selected.Some of the alternative locations are selected for distribution centers, which can provide integrated services and meet each region's requirements of the quantity of fresh agricultural products and carbon emission, with the overall costs of location selection and carbon emissions of warehousing and distribution minimized [34][35][36][37][38].The hypotheses in the model are shown in Table 1.
Table 1.Hypothetical premise in the model.

Supply and Demand
(1) the inventory of a single distribution center can satisfy the demand of covering demand points, and the inventory is equal to the sum of the demand of all demand points (2) the demand for fresh agricultural products at each demand point is fixed and will not fluctuate greatly in the short term Transportation (1) alternative locations are equipped with sufficient refrigerated vehicles for transport, regardless of shipping batches and fixed costs of vehicles (2) refrigerated vehicles travel at a constant speed, regardless of service time window constraints and vehicles routing Carbon Emissions (1) the site selection process involves static carbon emissions from distribution centers and dynamic carbon emissions from refrigerated vehicle transportation (2) carbon emissions from product decay are not considered Freshness (1) the freshness of fresh agricultural products is only related to the transportation time (2) the effect of temperature change on freshness is not considered

Symbols and Variables
The symbols and variables used in the modeling process of this paper are shown in Table 2.In addition, the intermediate variables that appear for the first time in this paper are also defined.Fresh agricultural products are vulnerable during the process of warehousing and distribution and need to be delivered on time.In addition, the rotting of such products accelerates as time goes by.Φ is the sensitivity coefficient of freshness to time.The greater Φ is, the faster the products rot in an exponential manner [39].In line with the hypothesis that the freshness of fresh agricultural products is only correlated with the time of transportation after the products have left the refrigerators in the alternative locations of the distribution centers, the function of freshness is R = 1 − e −φ∆t and Figure 2 shows the function relation of fresh agricultural products' freshness with the time change.Specifically, R0, Rbest, Rworst represent the freshness of agricultural p transported to a demand point, the freshness expected by the demand poin est level of freshness that the demand point can accept respectively.R0 > refrigerated vehicle arrives at the demand point at t1, the demand point is 1 with a corresponding potential incentive cost β2 generated.Rworst < R0 < Rbe frigerated vehicle arrives at the demand point at t2, the demand point is not 1 with a corresponding penalty cost β1 generated.R0 < Rworst when the refrig arrives at the demand point at t3, the demand point is not satisfied at all, w penalty cost generated.The penalty cost of freshness (CR) is calculated as sh tion (1).

The Overall Carbon Emission Costs
The cold chain is needed in various parts of the storage and transpor agricultural products, and devices such as cold storage, refrigerated truc mated temperature control lead to large quantities of carbon emissions.T vided the carbon emissions during the whole process of cold chain wareho tribution into static carbon emissions of cold chain distribution centers and bon emissions during refrigerated vehicle transportation [40,41].
(1) The static carbon emissions of cold chain distribution centers This type of emissions mainly comes from daily storage and tempera the cold chain.The former includes diesel and electricity energy consumed and equipment, and the latter part is largely electricity energy consumptio tion of the carbon emissions generated by the two parts of electricity energy CE1 is shown in Equation ( 2).

CE EC e
According to the mass balance model revised by the Intergovernmenta mate Change (IPCC), the carbon emissions generated by the diesel consu measured [42]: carbon emissions = (diesel input × carbon content per unit Specifically, R 0 , R best , R worst represent the freshness of agricultural products when transported to a demand point, the freshness expected by the demand point and the lowest level of freshness that the demand point can accept respectively.R 0 > R best when the refrigerated vehicle arrives at the demand point at t 1 , the demand point is 100% satisfied, with a corresponding potential incentive cost β 2 generated.R worst < R 0 < R best when the refrigerated vehicle arrives at the demand point at t 2 , the demand point is not 100% satisfied, with a corresponding penalty cost β 1 generated.R 0 < R worst when the refrigerated vehicle arrives at the demand point at t 3 , the demand point is not satisfied at all, with an infinite penalty cost generated.The penalty cost of freshness (CR) is calculated as shown in Equation (1).

The Overall Carbon Emission Costs
The cold chain is needed in various parts of the storage and transportation of fresh agricultural products, and devices such as cold storage, refrigerated trucks, and automated temperature control lead to large quantities of carbon emissions.This paper divided the carbon emissions during the whole process of cold chain warehousing and distribution into static carbon emissions of cold chain distribution centers and dynamic carbon emissions during refrigerated vehicle transportation [40,41].
(1) The static carbon emissions of cold chain distribution centers This type of emissions mainly comes from daily storage and temperature control of the cold chain.The former includes diesel and electricity energy consumed by machinery and equipment, and the latter part is largely electricity energy consumption.The calculation of the carbon emissions generated by the two parts of electricity energy consumption CE 1 is shown in Equation ( 2).
According to the mass balance model revised by the Intergovernmental Panel on Climate Change (IPCC), the carbon emissions generated by the diesel consumption can be measured [42]: carbon emissions = (diesel input × carbon content per unit calorific value of diesel × the transforming factors per unit calorific value of diesel-fixed carbon content of the product-fixed carbon content of waste) × carbon oxidation factor × 44/12.This paper does not include the calculation of fixed carbon content of products and waste during the warehousing and distribution process of fresh agricultural products because they are not easy to measure.The amount of diesel consumed in daily storage operations CE 2 is calculated as shown in Equation (3).
In summary, the amount of static carbon emissions of cold chain distribution centers CE s is calculated as shown in Equation (4).
(2) Dynamic carbon emissions during refrigerated vehicle transportation This type of emissions mainly comes from the diesel consumption and refrigeration equipment of the vehicles.In this paper, both parts are regarded as diesel consumption that is directly related to factors such as the load, weight and transport distance of the vehicles.
According to the linear relationship between the diesel consumption of refrigerated vehicles driving per unit distance (VP) and the load of the vehicles (Q), the diesel consumption of empty vehicles is VP 0 = a × Q 0 + b, and that of fully-loaded vehicles is Based on the diesel consumption of refrigerated vehicles during transportation, the amount of carbon emissions can be calculated as follows: carbon emissions = diesel consumption × the emission per unit of diesel consumption.According to the real-time carbon emission model [43], the diesel consumption of refrigerated vehicles driving per unit distance (VP) is calculated as shown in Equation ( 5).

VP = (
On the premise that the longitude and latitude plane coordinates of demand points are known, and in order to facilitate the study and avoid the influence of the non-ideal state of the roads for transportation, this paper adopted the Euclidean distance measurement method, and the straight-line distance from i to j is used to represent the transportation distance, as shown in Equation ( 6).
In summary, the dynamic carbon emissions during refrigerated vehicle transportation are calculated as shown in Equation (7).

Model Construction
The bi-objective function includes the minimum total costs MinC 1 and minimum carbon emissions MinC 2 , as shown in Equations ( 8) and ( 9), respectively.
The objective function (8) minimises the total costs of location selection, which consists of four parts: the storage costs (the reconstruction costs of cold chain distribution centers and the storage costs of fresh agricultural products), transport costs (the freight cost of transporting fresh agricultural products by refrigerated vehicles), the penalty costs of the freshness of products when delivered to the demand points and the overall carbon emission cost (static carbon emissions of cold chain distribution centers and dynamic carbon emissions during refrigerated vehicle transportation).The objective function ( 9) minimises the overall carbon emissions of the warehousing and distribution process, including static carbon emissions of cold chain distribution centers and dynamic carbon emissions during refrigerated vehicle transportation.Constraint (10) ensures that the number of cold chain distribution centers does not exceed the number of alternative locations.Constraint (11) ensures that the demands of all demand points for fresh agricultural products should be met.Constraint (12) ensures that if an alternative location provides goods for some demand points, a distribution center needs to be established at that location, where M is a large positive number.Constraint (13) ensures that at least one alternative location i ensures warehousing and distribution services for the demand point j.Constraint (14) ensures that the arrival time of the vehicles in the penalty costs of fresh agricultural costs lies between the time of the most satisfactory freshness of the demand point (t best ) and the time of the lowest level of freshness that the demand point can accept (t worst ).Constraint (15) ensures the restrictions on the values of the variables.

Model Predictions
The location model of cold chain distribution centers of fresh agricultural products under carbon emission constraints is an NP-hard problem with multiple constraints and multiple objective functions, which is difficult to solve by a single exact algorithm.Therefore, this paper designed a two-phase heuristic algorithm to solve the problem.In the first phase of the heuristic algorithm, the DBSCAN spatial clustering algorithm was used to cluster and divide the demand points within the service coverage of the alternative locations for distribution centers to obtain the initial solution set.In the second phase of the heuristic algorithm operation, an improved WOA was used to carry out location iterations for the initial solution set, reduce the local optimum and convergence error, and obtain the overall optimal location selection scheme.

DBSCAN Spatial Clustering Algorithm for Initial Solution Clustering
This paper introduced the DBSACAN spatial clustering algorithm for the clustering process of demand points of fresh agricultural products.On the one hand, the algorithm can accurately determine the core points, boundary points and noise points, so as to further select the standard points and improve the iterative fault tolerance and solution accuracy of the algorithm [44,45], as shown in Figure 3; on the other hand, the algorithm has no centroid, and only forms a solution set by continuously connecting high-density points in the neighborhood, which can form adaptive clustering results of different shapes and sizes.
process of demand points of fresh agricultural products.On the one hand, the algorithm can accurately determine the core points, boundary points and noise points, so as to further select the standard points and improve the iterative fault tolerance and solution accuracy of the algorithm [44,45], as shown in Figure 3; on the other hand, the algorithm has no centroid, and only forms a solution set by continuously connecting high-density points in the neighborhood, which can form adaptive clustering results of different shapes and sizes.{ } Equation ( 16) represents the neighbourhood of the initial solution x within the radius Epsilon and y is the reference point in the set of initial solution.d(x, y) stands for the distance from the initial solutions to the reference points in the set.Equation ( 17) represents the neighbourhood density of x and Equation ( 18) is the condition for the clustering of core point xc.All core points in set X were set Xc and the set composed of non-core points was Xnc.

Border, Standard and Noise Point Clustering
Based on the core points, when the number of points contained in a demand point within the neighbourhood of radius Epsilon is lower than the neighbourhood density threshold minPts, there are two cases: if the point is equal to the neighbourhood range of a certain core point xc, then the point is a boundary point xbd and the constituted set is set to Xbd, as shown in Equation (19); if the point is smaller than the neighbourhood range of a certain core point xc then the point is a standard point xst and the set formed is set to Xst, as shown in Equation (20).If a demand point does not belong to the set of core points Xc, does not belong to the set of boundary points Xbd, and also does not belong to the set of standard points Xst, then the point is a noise point xn, as shown in Equation (21).N epsilon (x) = { y ∈ X|d(x, y) ≤ epsilon} ( 16) Equation ( 16) represents the neighbourhood of the initial solution x within the radius Epsilon and y is the reference point in the set of initial solution.d(x, y) stands for the distance from the initial solutions to the reference points in the set.Equation ( 17) represents the neighbourhood density of x and Equation ( 18) is the condition for the clustering of core point x c .All core points in set X were set X c and the set composed of non-core points was X nc .

Border, Standard and Noise Point Clustering
Based on the core points, when the number of points contained in a demand point within the neighbourhood of radius Epsilon is lower than the neighbourhood density threshold minPts, there are two cases: if the point is equal to the neighbourhood range of a certain core point x c , then the point is a boundary point x bd and the constituted set is set to X bd , as shown in Equation (19); if the point is smaller than the neighbourhood range of a certain core point x c then the point is a standard point x st and the set formed is set to X st , as shown in Equation (20).If a demand point does not belong to the set of core points X c , does not belong to the set of boundary points X bd , and also does not belong to the set of standard points X st , then the point is a noise point x n , as shown in Equation ( 21).
x The whale optimization algorithm (WOA) is a novel heuristic algorithm of swarm intelligence proposed by Australian scholar Mirjalili in 2016 according to the behaviour of humpback whales rounding up prey [46].The search process of the algorithm is largely divided into three stages: encircling, bubble-net attacking and spiral-shaped hunting.A solution can be represented by the position coordinates of a whale, and the process of searching for a solution can be seen as the process of several whales continuously changing their positions until the optimal solution is obtained.
(1) Encircling The algorithm simulates a group of humpback whales working together to locate their prey and encircle the initial solution of the target.The algorithmic model for this stage is as follows.
Equation ( 22) represents the updated position vector of other humpback whales x j after being influenced by the best humpback whale X * , and k represents the number of iterations.Equation (23) represents the distance between other whales and the best whale X * before the encirclement.Equation (24) represents that A and C are the iteration matrix coefficients.r 1 and r 2 are random numbers from 0 to 1. Equation (25) represents that a is an iteration parameter that decreases linearly from 2 to 0 as the number of iterations t increases.Humpback whales obtain food by controlling the |A 1 | vector.When |A 1 | ≤ 1, the whales in other positions approach the current best one.
(2) Bubble-net attacking As the other humpback whales move towards the best one X * to encircle the prey, they not only contract the encirclement but also follow a logarithmic spiral towards the prey.To synchronize the algorithmic operations of encircling the prey and swimming up the spiral, it was assumed that encircling is selected with 50% probability p and swimming in a spiral with 50% probability 1 − p. Equation (26) shows the mathematical model.
Parameter b is a logarithmic spiral constant and l is a random number between −1 and 1.When l is close to −1, humpback whales approach the prey.Otherwise, the whales swim away.
(3) Spiral-shaped hunting The swimming range of a humpback whale may be |A 1 | > 1.At this time, other humpback whales randomly select individuals in the group to approach, which facilitates the algorithm to obtain the overall optimal solution.From this, the position iteration mechanism of the WOA algorithm under the influence of |A 1 | can be obtained, as shown in Figure 4.
The random individual positions of the current humpback whales were assumed as X rand and the location of the whale to swim as X j .The mathematical model of the behaviour of spiral-shaped hunting is shown in Equations ( 27) and (28).
mechanism of the WOA algorithm under the influence of |A1| can be obtained, as shown in Figure 4.The random individual positions of the current humpback whales were assumed as X rand and the location of the whale to swim as X j .The mathematical model of the behaviour of spiral-shaped hunting is shown in Equations ( 27) and ( 28).

Using Sort Matrices to Calculate the Overall Fitness of the Bi-objective Function
During the process of the location selection for cold chain distribution centers of fresh agricultural products, multiple constraints and multiple objective functions are involved, so it is difficult to measure the comparative priority of minimum total costs and minimum carbon emissions.Therefore, this paper used the comparison sort performance matrix to calculate the overall fitness of the alternative locations, in order to measure the relationship between the total costs and carbon emissions [47].First of all, the initial value of the bi-objective function of the initial solution set Xc of the location selection for cold chain distribution centers was calculated, and the performance comparison sort of the sort matrices was carried out, which was further transformed into the overall fitness of the alternative locations.Assuming that the total costs of the location selection and the total carbon emissions during the warehousing and distribution process are Obj1 and Obj2.Their ranking numbers are OR1 and OR2 in ascending order, and the values of fitness are f1 and f2, as shown in Table 3.

Object Alternative
The results of the performance comparison sort of the sort matrices were transformed into the individual fitness of Equations ( 8) and (9), in order to calculate the overall fitness of the alternative locations.The algorithm model at this stage is as follows:

Using Sort Matrices to Calculate the Overall Fitness of the Bi-objective Function
During the process of the location selection for cold chain distribution centers of fresh agricultural products, multiple constraints and multiple objective functions are involved, so it is difficult to measure the comparative priority of minimum total costs and minimum carbon emissions.Therefore, this paper used the comparison sort performance matrix to calculate the overall fitness of the alternative locations, in order to measure the relationship between the total costs and carbon emissions [47].First of all, the initial value of the bi-objective function of the initial solution set X c of the location selection for cold chain distribution centers was calculated, and the performance comparison sort of the sort matrices was carried out, which was further transformed into the overall fitness of the alternative locations.Assuming that the total costs of the location selection and the total carbon emissions during the warehousing and distribution process are Obj 1 and Obj 2 .Their ranking numbers are OR 1 and OR 2 in ascending order, and the values of fitness are f1 and f2, as shown in Table 3.
Table 3. Performance comparison ranking of ranking matrices for bi-objective function.

Object Alternative
Al(1) Al( 2) The results of the performance comparison sort of the sort matrices were transformed into the individual fitness of Equations ( 8) and ( 9), in order to calculate the overall fitness of the alternative locations.The algorithm model at this stage is as follows: Equation ( 29) represents the fitness value of the jth objective function.OR j (i) stands for the results of comparison sort.M is the total number of alternative locations.Equation (30) represents the overall fitness of alternative locations.K is a variable within the range of 1 to 2, which is used to increase the fitness of the optimal individual.

Adopting the Mutation-Induced Perturbation of the Cauchy Distribution
When the humpback whales conduct a global search in the early iterations of the algorithm, the position of the best humpback whale is needed as a guide for other whales to approach and encircle.Since the location of the best one is selected randomly by the algorithm, there is a possibility that the global optimal solution is missed.Thus, this paper adopted the mutation-induced perturbation that conforms to the Cauchy distribution [48][49][50].Because the Cauchy distribution features wide distribution at both ends, it can make the distribution at both ends of the mutation more easily jump out of the local optimum and give full play to the effect of the mutation at both ends.In addition, during the global search by WOA, the target humpback whales are changing their coordinate positions to effectively compensate for the blindness of the mutations that conforms to the Cauchy distribution.Equation (31) shows the calculation of the Cauchy inverse cumulative distribution function. )) Therefore, Equation ( 22) can therefore be improved to Equation (32).
In particular, F −1 represents the inverse cumulative distribution function of the Cauchy distribution.X mn is the nth location of the mth humpback whale between the Cauchy mutation.γ = A. p is uniformly distributed within the range of 0 to 1.

Introducing Sine and Cosine Inertia Weights
Humpback whales control the |A 1 | vector to move closer or farther away from the best humpback whales so as to realize the three stages of encircling, bubble-net attacking and spiral-shaped hunting.In other words, the WOA controls the search capability through the |A 1 | vector.The key iteration parameter a in the |A 1 | vector decreases linearly as the number of iterations increases.The law of the linear decreasing leads to a reduction in the convergence range in the early iterations that require a larger search range, while causing a slowdown in the rate of convergence in the later iterations that require a smaller search range, thus affecting the overall quality of convergence of the algorithm.This paper, based on the characteristics of sine and cosine curves, introduced non-linear sine and cosine inertia weights to the iteration parameter a.The search range was expanded by the sine inertia weights and the rate of convergence was enhanced by the cosine inertia weights, in order to improve the quality of convergence and population diversity of the algorithm [51][52][53].The formulas of sine and cosine inertia weights are shown in Equations ( 33) and (34).
The control parameter r 3 represents a random number in the range of 0 to 2, which is used to control the influence of the relative positions of the random humpback whales and the optimal humpback whale.The random parameter rand represents a random number in the range of 0 to 1, which is used to improve the population diversity of the algorithm.t is the current number of iterations and T is the total number of iterations, and the division of the two makes the iteration parameter a somewhat adaptive.Specifically, in the early iterations of the algorithm, the value of t/T is small and the inertia weight w is large; in the later iterations of the algorithm, t gradually converges to T, so the value of t/T is large and the inertia weight w is small.Equation ( 32) can be therefore improved to Equation (35).

−−→
)) (35) Figure 5 shows the influence of the |A 1 | vector under sine and cosine inertia weights on the global search and local convergence of the hunting of the humpback whales in different iteration cycles of the one-dimensional hunting space.
the early iterations of the algorithm, the value of t/T is small and the large; in the later iterations of the algorithm, t gradually converges to T is large and the inertia weight w is small.Equation ( 32) can be therefore tion (35).
Figure 5 shows the influence of the |A1| vector under sine and co on the global search and local convergence of the hunting of the hu different iteration cycles of the one-dimensional hunting space.The solid line represents the situation of shunting under the influ weights.The target humpback whale moves in a sinusoidal logarithm guided by the random individual X rand of the population.A global se before the stage of spiral-shaped hunting, reducing the blind spot of vergence range of the cosine inertia weights and avoiding the situatio tial optimal solution is missed.The dotted line represents the situation the influence of cosine inertia weights.The target humpback whale mo arithmic spiral trajectory guided by the best humpback whale X*.A takes place at the late stage of spiral-shaped hunting, which compe convergence of the global search of sine inertia weights and improves ciency of the algorithm.

Solving Steps of Two-Stage Heuristic Algorithm
The specific algorithm flow is shown in Figure 6.The solid line represents the situation of shunting under the influence of sine inertia weights.The target humpback whale moves in a sinusoidal logarithmic spiral trajectory guided by the random individual X rand of the population.A global search is carried out before the stage of spiral-shaped hunting, reducing the blind spot of the small local convergence range of the cosine inertia weights and avoiding the situation where the potential optimal solution is missed.The dotted line represents the situation of shunting under the influence of cosine inertia weights.The target humpback whale moves in a cosine logarithmic spiral trajectory guided by the best humpback whale X*.A local convergence takes place at the late stage of spiral-shaped hunting, which compensates for the slow convergence of the global search of sine inertia weights and improves the operation efficiency of the algorithm.

Solving Steps of Two-Stage Heuristic Algorithm
The specific algorithm flow is shown in Figure 6.

Case Description and Data Acquisition
To test the effectiveness and feasibility of the two-phase heuristic algorithm in the actual location selection process, a case study was selected as follows.Y enterprise in Zhejiang Province, in order to meet the daily demands of fresh agricultural products of 36 regional demand points [54,55], decided to establish a number of cold chain distribution centers in order to provide integrated services, taking carbon emission and sustainable development into account.Among the 36 demand points, 10 were selected as alternative locations, including 1 Yuhang District, 5 Qiantang District, 7 Tonglu County, 10 Fenghua District, 11 Yuyao City, 22 Tongxiang City, 26 Pujiang County, 29 Linhai City, 32 Tiantai County and 34 Jinyun County, among which six locations had to be selected to build a cold chain distribution center and optimize the address.Table 4 shows the Zhejiang Y enterprise's demand points for fresh agricultural products and the corresponding demands.The demand points of the Zhejiang Y enterprise distributed in 36 regions were divided, and the demands in each region were summed up as the transportation volume of the demand point.To facilitate the solution of the straight-line distance dij from the alternative location i of a cold chain distribution center to the target demand point j, this paper used the ArcGIS software to convert the latitude and longitude coordinate data into plane coordinate data, as shown in Figure 7.In the actual location selection, high storage cost per unit of goods and the fixed reconstruction costs of the alternative locations are important factors affecting the total costs.The data related to the alternative locations for cold chain distribution centers of fresh agricultural products of the Zhejiang Y enterprise are shown in Table 5.The demand points of the Zhejiang Y enterprise distributed in 36 regions were divided, and the demands in each region were summed up as the transportation volume of the demand point.To facilitate the solution of the straight-line distance d ij from the alternative location i of a cold chain distribution center to the target demand point j, this paper used the ArcGIS software to convert the latitude and longitude coordinate data into plane coordinate data, as shown in Figure 7.In the actual location selection, high storage cost per unit of goods and the fixed reconstruction costs of the alternative locations are important factors affecting the total costs.The data related to the alternative locations for cold chain distribution centers of fresh agricultural products of the Zhejiang Y enterprise are shown in Table 5.According to the relevant data, the relevant parameters assumed of refrigerated vehicles in the process of cold chain warehousing and distribution of fresh agricultural products are shown in Table 6.According to IPCC's proposal, the carbon price was USD 30 per ton in 2020 [56] (USD 0.03/kg•CO 2 ).Therefore, this paper assumes that the current cost per unit of carbon emissions in Zhejiang Province is 0.03 dollar/kg•CO 2 .Based on relevant data, the values of other constants are assumed as shown in Table 7.In order to verify the efficiency of the two-phase heuristic algorithm (THA) proposed in this paper for solving the location selection for cold chain distribution centers of fresh agricultural products, a simulation comparison was conducted on THA, the basic whale optimization algorithm [57,58] (WOA) and the basic particle swarm optimization algorithm [59,60] (PSO) through the case of Y enterprise in Zhejiang Province.The experimental platform for the simulation was Win11, and MATLAB R2021b software was used for the operation.
The number of population N was 36, and the maximum number of iterations T was 200.The neighbourhood density radius Epsilon was 100 km, and the neighbourhood density threshold minPts was a natural number 6 in the DBSCAN algorithm.The probability of encirclement p was 0.5 and the probability of mutation-induced perturbation of the Cauchy Distribution pγ was 0.3 in the IWOA algorithm.The locations of cold chain distribution centers obtained by the three algorithms and the final optimization results of the experiment are shown in Table 8.The effect pictures of the location selection and the iterative convergence curves of the three algorithms are shown in Figures 8-10. Figure 11 presents the iterative convergence contrast curves of the three algorithms, and Figure 12 shows the changing curves of the overall carbon emissions during the location selection process of the three algorithms.According to the comparison of the performance of the three algorithms in the location selection shown in Table 8, in terms of the total costs of the location selection, the THA algorithm's average optimization is 7.70% compared with the WOA algorithm and 12.77% compared with the PSO algorithm.In terms of the overall carbon emissions, the THA algorithm's average optimization rate is 10.23% compared with the WOA algorithm and 14.12% compared with the PSO algorithm.It is suggested that the THA algorithm matches a better location selection scheme with fewer iterations under the constraint of the same cost per unit of carbon emissions.
According to the location selection scheme and the comparison of iterative convergence shown in Figures 8-12, the search process of the WOA and PSO algorithms is significantly broken and stagnant in the early iterations, indicating that the local optimum has occurred, while the iterative performance of the THA algorithm is more continuous and smooth, and the fitness decreases more quickly, which means it has a more efficient iteration capability.At the late stage of algorithm iteration, both the WOA and PSO algorithms enter a bottleneck to some extent and differ significantly from the THA algorithm in terms of effectiveness and continuity of iteration.

Sensitivity Analysis Based on the Cost per Unit of Carbon Emissions
In order to verify the impact of the cost per unit of carbon emissions on the total costs of location selection during the process of warehousing and distribution, this paper used an analysis method with control variables and a two-phase heuristic algorithm to change the cost per unit of carbon emissions to carry out a multi-group controlled experiment while other constant parameters remained unchanged.Assuming a cost per unit of carbon emissions of 0.02 dollar/kg as the starting point and an incremental rule of 0.01 dollar/kg, a According to the comparison of the performance of the three algorithms in the location selection shown in Table 8, in terms of the total costs of the location selection, the THA algorithm's average optimization is 7.70% compared with the WOA algorithm and 12.77% compared with the PSO algorithm.In terms of the overall carbon emissions, the THA algorithm's average optimization rate is 10.23% compared with the WOA algorithm and 14.12% compared with the PSO algorithm.It is suggested that the THA algorithm matches a better location selection scheme with fewer iterations under the constraint of the same cost per unit of carbon emissions.
According to the location selection scheme and the comparison of iterative convergence shown in Figures 8-12, the search process of the WOA and PSO algorithms is significantly broken and stagnant in the early iterations, indicating that the local optimum has occurred, while the iterative performance of the THA algorithm is more continuous and smooth, and the fitness decreases more quickly, which means it has a more efficient iteration capability.At the late stage of algorithm iteration, both the WOA and PSO algorithms enter a bottleneck to some extent and differ significantly from the THA algorithm in terms of effectiveness and continuity of iteration.

Sensitivity Analysis Based on the Cost per Unit of Carbon Emissions
In order to verify the impact of the cost per unit of carbon emissions on the total costs of location selection during the process of warehousing and distribution, this paper used an analysis method with control variables and a two-phase heuristic algorithm to change the cost per unit of carbon emissions to carry out a multi-group controlled experiment while other constant parameters remained unchanged.Assuming a cost per unit of carbon emissions of 0.02 dollar/kg as the starting point and an incremental rule of 0.01 dollar/kg, a sensitivity analysis based on carbon emissions was conducted and the changes in the overall carbon emission costs and the total costs of location selection are shown in Table 9. Note: The difference in the total costs is calculated based on the total costs of the location selection at the price per unit of carbon emissions of 0.03 dollar/kg.
According to Table 9, as the cost per unit of carbon emissions increases, both the overall carbon emission costs and the total costs of location selection show an upward trend, with the former accounting for a bigger proportion of the total costs.It is further verified that the cost per unit of carbon emissions and the total costs are significantly positively correlated, and the two-phase heuristic algorithm can effectively obtain locations for cold chain distribution centers of fresh agricultural products under carbon emission constraints.As people pay much attention to environmental protection, the cost per unit of carbon emissions will continue to rise [61].It will be an inevitable choice for enterprises to take efficient energy-saving and emission-reduction measures in order to reduce costs and increase efficiency in the future.

Conclusions
Based on the traditional location for distribution centers, this paper optimized for both freshness and carbon emissions.A bi-objective function model was constructed to minimize the total cost and the carbon emission.The total cost consists of four com-ponents: the cost of storage, the cost of transportation, the cost of freshness penalty and the comprehensive carbon emission cost.The carbon emission consists of two compo-nents: the static carbon emission of the cold chain distribution centers and the dynamic carbon emission of the refrigerated vehicles transportation.A two-stage heuristic algo-rithm was designed to solve the model: the first stage used the DBSCAN algorithm to divide the initial solution set and the second stage used the IWOA to obtain the overall optimal solution.
The results of the case of enterprise Y in Zhejiang Province verified the feasibility of the model and the main conclusions obtained include: first, the THA had an average total cost optimization rate of 7.70% and an average carbon emission optimization rate of 10.23% compared to the WOA, while it had rates of 12.77% and 14.12%, respectively, compared to the PSO.This shows that the THA can significantly optimize the solution results.Second, the unit carbon emission cost affected the location for cold chain distribution centers.When the unit carbon emission cost increased from 0.02 dollar/kg to 0.06 dollar/kg, the overall carbon emission cost increased by 42.07%, the total cost increased by 3.05% and the proportion of overall carbon emission cost in the total cost increased from 2.21% to 3.05%.Therefore, logistics enterprises can reduce the total cost through reasonable location for fresh agricultural products cold chain distribution centers [45,46].
Further research can apply the two-stage heuristic algorithm to the two-level planning for location of fresh agricultural products cold chain distribution centers and the route optimization of refrigerated vehicles [62][63][64][65].The impact of time window constraints at demand points [66] and traffic condition constraints [67] on the planning of cold chain logistics systems can also be considered.In addition, the author will focus on this type of problem in the future.

Figure 1 .
Figure 1.The overall research framework diagram.

Figure 1 .
Figure 1.The overall research framework diagram.

Figure 2 .
Figure 2. The functional relationship between freshness of fresh agricultural produ

Figure 2 .
Figure 2. The functional relationship between freshness of fresh agricultural products and time.

Figure 3 .
Figure 3.The relative distribution position of the core parameters of the DBSCAN algorithm.4.1.1.Core Point Clustering With each initial solution x as the center of a circle and Epsilon as the radius to form a neighbourhood, the solutions contained in the neighbourhood were counted.The core point xc was generated by comparison with min Pts.The equations are as follows.

Figure 3 .
Figure 3.The relative distribution position of the core parameters of the DBSCAN algorithm.4.1.1.Core Point ClusteringWith each initial solution x as the center of a circle and Epsilon as the radius to form a neighbourhood, the solutions contained in the neighbourhood were counted.The core point x c was generated by comparison with min Pts.The equations are as follows.

Figure 4 .
Figure 4.The position iteration mechanism controlled by |A 1 | in WOA algorithm.

Figure 5 .
Figure 5.The solution mechanism of WOA algorithm under sin-cos inertial w sional space.

Figure 5 .
Figure 5.The solution mechanism of WOA algorithm under sin-cos inertial weight in onedimensional space.

Figure 6 .
Figure 6.The two-stage heuristic algorithm flow.Figure 6.The two-stage heuristic algorithm flow.

Figure 6 .
Figure 6.The two-stage heuristic algorithm flow.Figure 6.The two-stage heuristic algorithm flow.

Figure 7 .
Figure 7.The Plane cartesian coordinate system based on latitude and longitude of Y enterprise's demand points.

Figure 7 .
Figure 7.The Plane cartesian coordinate system based on latitude and longitude of Y enterprise's demand points.

Figure 8 .
Figure 8.The effect pictures of location selection and iterative convergence curves of the THA algorithm.

Figure 9 .
Figure 9.The effect pictures of location selection and iterative convergence curves of the WOA algorithm.

Figure 10 .
Figure 10.The effect pictures of location selection and iterative convergence curves of the PSO algorithm.

Figure 8 .
Figure 8.The effect pictures of location selection and iterative convergence curves of the THA algorithm.

Figure 8 .
Figure 8.The effect pictures of location selection and iterative convergence curves of the THA algorithm.

Figure 9 .
Figure 9.The effect pictures of location selection and iterative convergence curves of the WOA algorithm.

Figure 10 .
Figure 10.The effect pictures of location selection and iterative convergence curves of the PSO algorithm.

Figure 9 .
Figure 9.The effect pictures of location selection and iterative convergence curves of the WOA algorithm.

Figure 8 .
Figure 8.The effect pictures of location selection and iterative convergence curves of the THA algorithm.

Figure 9 .
Figure 9.The effect pictures of location selection and iterative convergence curves of the WOA algorithm.

Figure 10 .
Figure 10.The effect pictures of location selection and iterative convergence curves of the PSO algorithm.Figure 10.The effect pictures of location selection and iterative convergence curves of the PSO algorithm.

Figure 10 . 24 Figure 11 .
Figure 10.The effect pictures of location selection and iterative convergence curves of the PSO algorithm.Figure 10.The effect pictures of location selection and iterative convergence curves of the PSO algorithm.Sustainability 2022, 14, x FOR PEER REVIEW 20 of 24

Figure 11 .
Figure 11.The iterative convergence contrast curves of the three algorithms.Figure 11.The iterative convergence contrast curves of the three algorithms.

Figure 11 .
Figure 11.The iterative convergence contrast curves of the three algorithms.

Figure 12 .
Figure 12.The overall carbon emissions contrast curves of the three algorithms.

Figure 12 .
Figure 12.The overall carbon emissions contrast curves of the three algorithms.

Author Contributions:
All of the authors contributed to the work in the paper.Specifically: Conceptualization, H.W. and H.R.; methodology, H.W. and H.R.; software, H.W., H.R. and X.D.; validation, H.W., H.R. and X.D.; data curation, H.W. and X.D.; writing-original draft preparation, H.W. and H.R.; writing-review and editing, H.W. and H.R.; funding acquisition, H.W. and X.D.All authors have read and agreed to the published version of the manuscript.Funding: The research was supported by Shandong Province Social Science Project (No. 18CGLJ35), 2019 Shandong Humanities and Social Sciences Project (No. 19-ZZ-GL-09), 2019 Qingdao Social Science Planning Research Project (QDSKL1901168).

Table 1 .
Hypothetical premise in the model.

Table 2 .
Symbols and variables in the model.

Table 5 .
Data related to Y enterprise's alternative locations for cold chain distribution centers.

Table 4 .
Y enterprise's demand points for fresh agricultural products and the corresponding demands.

Table 5 .
Data related to Y enterprise's alternative locations for cold chain distribution centers.

Table 6 .
Relevant parameters of Y enterprise's refrigerated vehicles.

Table 7 .
Other constant values related to the model.

Table 8 .
Comparison of the performance of the three algorithms in the location selection.

Table 9 .
The impact of the changes in cost per unit of carbon emissions on the total costs.