A Hyperheuristic Approach for Location-Routing Problem of Cold Chain Logistics considering Fuel Consumption

In response to violent market competition and demand for low-carbon economy, cold chain logistics companies have to pay attention to customer satisfaction and carbon emission for better development. In this paper, a biobjective mathematical model is established for cold chain logistics network in consideration of economic, social, and environmental benefits; in other words, the total cost and distribution period of cold chain logistics are optimized, while the total cost consists of cargo damage cost, refrigeration cost of refrigeration equipment, transportation cost, fuel consumption cost, penalty cost of time window, and operation cost of distribution centres. One multiobjective hyperheuristic optimization framework is proposed to address this multiobjective problem. In the framework, four selection strategies and four acceptance criteria for solution set are proposed to improve the performance of the multiobjective hyperheuristic framework. As known from a comparative study, the proposed algorithm had better overall performance than NSGA-II. Furthermore, instances of cold chain logistics are modelled and solved, and the resulting Pareto solution set offers diverse options for a decision maker to select an appropriate cold chain logistics distribution network in the interest of the logistics company.


Introduction
Cold chain logistics are developing rapidly with constant improvement of living standard of the people and increasing demand for fresh food [1]. In order to lower distribution cost of cold chain logistics effectively, some researchers began to study optimization of cold chain logistics distribution network [2]. However, when optimizing cold chain logistics distribution network, one has to consider economic benefit of logistics distribution, distribution-derived carbon emission and other environmental benefits, and timeliness of the logistics network.
In a practical cold chain logistics network, more than one objective may be considered in the system, including the minimum total cost, the smallest number of vehicles, the maximum customer satisfaction, the shortest travel path, and the least distribution time. Numerous researchers have studied multiobjective logistics optimization problem. With the shortest travel path and carbon dioxide emission as optimization objectives, Jemai et al. [3] established a vehicle routing problem (VRP) based on green logistics and solved it by NSGA-II optimization. Based on optimization objective functions of the minimum cost and the maximum customer satisfaction, Liu et al. [4] developed a mathematical model for multiobjective location-routing problem (LRP). Molina et al. [5] developed a heterogeneous fleet VRP model with the minimized cost and emission as objectives. Golmohammadi et al. [6] proposed a mathematical model for multiobjective LRPs with the minimum storage distribution cost and difference in vehicle-traveled distance as optimization objective functions. e above literatures are related to multiobjective logistics distribution network, but multiple performance metrics of logistics network are not taken into account: economic benefit, environmental benefit, and timeliness. Moreover, there are a dearth of literatures on LRP-based models for cold chain logistics. Wang et al. [1] proposed developing a single-objective model in consideration of carbon emission for cold chain logistics and solving it using hybrid genetic algorithm. However, logistics network timeliness has not yet been considered.
Hence, in this paper, properties of cold chain logistics are considered from the perspective of LRP-based model; not only economic and environmental benefits but also distribution timeliness are taken into account; thereby a biobjective LRP of cold chain logistics considering fuel consumption and distribution period is proposed, and one multiobjective hyperheuristic (MOHH) is proposed to model and solve this problem. erefore, main contributions of this paper are described as follows: (i) Problem model: with logistics cost and distribution period as objective functions, a multiobjective mathematical model for cold chain logistics is developed considering environmental benefit arising from fuel consumption and social benefits like customer satisfaction arising from customer time window. (ii) Solution algorithm: for the above model, MOHH is designed. Based on the framework nature of hyperheuristics, four selection strategies for lowlevel heuristics (LLHs) and four solution acceptance criteria are proposed. (iii) Management suggestions and comments: based on experimental results, several insights and suggestions about logistics distribution are provided from management perspective.

Multiobjective Optimization Problem (MOP).
One typical MOP normally contains more than 2 conflicting or contradictory objectives and is also called multicriteria optimization problem [7]. Without loss of generality, one minimized MOP having D-dimensional decision variables and m-dimensional subobjectives can be expressed as where X � (x1, x2, . . . , x D) ∈ R D is a D-dimensional decision vector and F(X) ∈ R m is an m-dimensional objective vector, defining m mapping functions of decision space towards objective space; R D and R m are decision space and objective space, respectively; g i (X) ≤ 0(i � 1, 2, . . . , q) are q inequality constraints; and h i (X) � 0 (i � 1, 2, . . ., p) are p equality constraints. On this basis, the following related definitions are given [7].
Definition 1 (feasible solution). For one X ∈ R D , if two classes of constraints in equation (1) are satisfied, then X will be deemed as a feasible solution.
Definition 2 (feasible solution set). A set of all feasible solutions is referred to as feasible solution set, that is, X f .

Definition 3 (Pareto dominant).
Suppose that X A and X B are two feasible solutions in the feasible solution set. When and only when the following constraints are satisfied, where X A dominate X B ; in other words, compared with X B , X A is a Pareto dominant, which can be expressed as X B ≻ X A .
Definition 4 (Pareto optimal solution). One solution X * ∈ X f is referred to as Pareto optimal solution (or nondominated solution), when and only when the following condition is satisfied: Definition 5 (Pareto optimal solution set). As a set of all Pareto optimal solutions, Pareto optimal solution set is defined as follows: Definition 6 (Pareto front). A curved surface, composed of objective function values corresponding to all Pareto optimal solutions in Pareto optimal solution set P * , is called Pareto front and denoted as PF * : 2.2. Hyperheuristics. Metaheuristics (e.g., grey wolf optimization algorithm [8], nature-inspired optimization algorithm [9,10], and quantum-behaved particle swarm optimization [11]) have been widely applied in optimization problems. However, to find solutions of bespoke nature in a new problem-based domain, it takes huge efforts and lengthy time to have a good command of effectiveness of domain operator combination [12]. Moreover, it is difficult to select the promising optimal parameters [13]. Although the machine learning methods have also been used to solve the above dilemmas, such as obstacle recognition based on machine learning [14], this paper focuses on the another strategy, namely, hyperheuristics, which are capable of solving the above problems effectively. Cowling et al. [15] defined "hyperheuristic" as "heuristics to choose heuristics" for the first time. Later, Burke et al. [16] expanded its definition: (i) heuristic selection and (ii) heuristic generation. erefore, selection of optional LLHs is generally decided from indicators of spread or convergence for a candidate solution versus parent solution. As known from the above literature, an LLH can be defined as an operator (crossover or mutation operator) or a metaheuristic (e.g., NSGA-II and SPEA2). For a selective heuristic, there exists one universal simple framework, as shown in Figure 1.
e black zone is a domain barrier for isolating highlevel heuristics (HLHs) from LLHs. e low-level problem domain contains data information of the real problem, LLHs used to search directly for problem space, and population information of the problem domain (chromosome, fitness, etc.). In high-level control domain, there are two strategies with different objectives: operator selection strategy and solution acceptance mechanism. e selection strategy is used to search for LLH-based space and monitor history information of LLH performance so as to select high-quality operators (for isolating any information related to the real problem), while the acceptance criteria judge whether parent solution S p should be replaced depending on quality of children solution (or solution set) S c or not and control search direction and convergence speed of the algorithm, where S cu is the current solution. In addition, between HLHs and LLHs, there are information transmitters used to exchange and transmit information irrelevant to the problem domain, including choice information, acceptance strategy judgment information, improvement rate contributed by LLH, operator runtime and frequency, and number of consecutive improvement failures of current solution. Moreover, it is to highlight that any decision for high-level domain (selection or acceptance strategy) has to be isolated from real problem domain (such as encoding, crossover, and mutation methods); otherwise, hyperheuristic principles and objectives will be disrupting; in other words, generality and transplantability of the highlevel control strategy are improved in order to apply to any problem domain.
Although there are plenty of multiobjective algorithms (such as NSGA-II, SPEA2, and multiobjective optimization based on an improved cross-entropy method [17]) for solving the practical problems in the real-world applications, this paper focuses on the MOHH, which has critical problems in the following aspects in its design: (1) Design of high-level selection strategy: performance of LLH may exhibit phased feature, that is, varying with search progress; for example, operator performance has inconsistent manifestations in the early and late phases of search. On this basis, researchers have designed different selection strategies depending on demand, such as choice function [18,19] and sliding window-based FRR-MAB [20] (2) Approach to children solution evaluation: for a single-objective problem, it suffices to acquire fitness improvement rate (FIR) of children relative to parent. But, for an MOP, objective space is expanded to two or more dimensions, so FIR cannot be used directly to characterize an approach of replacing parent with children. Dominance relationship of children versus parent is utilized [20]. Two-stage sorting method is used to evaluate children solutions and reward corresponding LLHs [18,19] (3) Solution acceptance mechanism: quality of solution acceptance mechanism has immediate impact on population diversification and intensification, thus influencing whether an algorithm is able to obtain satisfactory solutions. For single-objective problems, there are deterministic acceptance mechanisms such as All Moves (AM), Only Improving (OI), and noninferior solution acceptance and probabilistic (nondeterministic) acceptance mechanisms such as great deluge (GD), delay acceptance, simulated annealing (SA), and Monte Carlo (4) Domain barrier and information transmitter control: hyperheuristics are aimed at improving algorithm universality and transplantability so that the algorithm can be applied to real problems in different domains. So, there is crucial information in multiobjective problem domain: chromosome fitness; in other words, chromosome fitness can be used by high-level control strategy during multiobjective optimization likewise erefore, design methods are given in this paper for the above problems in MOHH design, respectively, so as to satisfy effectiveness and high performance of solving a biobjective LRP of cold chain logistics considering fuel consumption and distribution period.

Mathematical Model of Biobjective LRP of Cold
Chain Logistics considering Fuel Consumption and Distribution Period

Problem Formulation, Assumptions, and Definitions of
Variables. LRP of cold chain logistics studied in this paper is a distribution centre-to-customer supply chain. Suppose that the cold chain logistics network is a directed network composed of distribution centres and customer points; that is, N � {N C , N D } represents a set of all points on the logistics network, where N C represents a set of customer points; it is a weighted directed network where point-to-point distances are the weights, including the distance between a distribution centre and a customer and the distance between customers. Optimization of a logistics network consists of two major aspects: reasonable distribution centre selection and vehicle routing. To study instances of cold chain logistics network reasonably, the following assumptions are made: (3) Other Assumptive Conditions. Each distribution centre is equipped with vehicles; a distribution vehicle transports cold chain logistics products from each distribution centre to various customers and returns to the original distribution centre. e maximum load weight of each vehicle is Q, no demand of one customer exceeds the maximum load weight of each vehicle, demand of one customer can be satisfied by one vehicle only through distribution service, and one vehicle is able to serve many times. If drivers have the same driving habit and driving skill and each distribution centre has abundant inventory, then cargo shortage will not occur. Notations and variables of the model are defined as follows: candidate distribution centre set N D � {1, 2, . . ., m}, transportation distance corresponding to each edge is d ij and travel speed is v ij ; customer i ∈ N C has nonnegative cargo allocation demand D i and time window requirement (ET i , LT i ), while w i is customer service time; vehicle type is the same and the maximum capacity is CV; distribution centre i ∈ N D has an enabling cost of O i and the maximum capacity of CD i . e following are decision variables: x ijk is 1 if edge (i, j) serves vehicle k; otherwise x ijk is 0; y j is 1 when distribution centre j is enabled; otherwise y j is 0; z ij is 1 when customer i is served by distribution centre j; otherwise z ij is 0. e following are additional variables: L ijk represents dynamic load of vehicle k on edge (i, j); AT ik is the time point at which vehicle k arrives at node i. erefore, based on the above model assumptions and definitions of notations and variables, a multiobjective mathematical model can be given in this paper.

Model Composition.
e biobjective model in this paper is designed to optimize total logistics cost and distribution period; the total logistics cost consists of fuel consumption cost, cargo damage cost, transportation cost, refrigeration cost of refrigerated trucks, penalty cost of time window, and operation cost of distribution centres, as detailed below.
(1) Fuel Consumption Cost. According to a report on the website of the government of Japan [21], fuel consumption of one type of vehicles is directly proportional to load weight. In light of this, Xiao et al. [22] assumed that a vehicle runs at a constant speed, only considered impacts of load weight and traveled distance on vehicle fuel emission, and split vehicle weight into dead weight M 0 and load weight M 1 ; then Fuel Consumption Rate (FCR) can be calculated using the following formula: e maximum load weight that can be borne by a vehicle is defined as CV; then, FCRs at zero load and full load can be expressed as ρ 0 � αM 0 + b and ρ * � α(M 0 + CV) + b, respectively. us, FCR per unit distance is erefore, the fuel consumption generated by vehicle k under load L ijk is   Computational Intelligence and Neuroscience and, for this problem mentioned, fuel consumption cost is where c 0 represents unit fuel cost (Yuan/L).
(2) Cargo Damage Cost of Cold Chain Logistics. Based on the rationale of previous studies on cold chain logistics products, researchers incorporated cargo damage cost analysis into distribution routing optimization of cold chain logistics network. To facilitate the study, suppose that cold chain logistics products in the whole cold chain logistics distribution network are stored at a constant temperature in stable environment. e top consideration is quality degradation of cold chain logistics products over time, and according to food science principles combined with deterioration rate function of cold chain logistics products, Qi [23] induced a mathematical relationship between product quality change ∆Q and original quality Q 0 as follows: where t ij � d ij /v ij . Cargo damage of cold chain logistics is related to logistics distribution time only, and cargo damage cost is expressed as follows: Let φ be time sensitivity coefficient; its expression is (3) Refrigeration Cost of Refrigeration Equipment. In the course of distribution in cold chain logistics, a refrigeration equipment is not fully closed; its external factors such as solar radiation will result in a temperature gradient between inner and outer surfaces so as to increase thermal load inside the refrigeration equipment. e evaporator in the refrigeration equipment will refrigerate the cargo until overall thermal load is balanced, and eventually, constant temperature is achieved. Fu [24] performed thermal equilibrium analysis by thermal equilibrium method frequently used in energy consumption analysis of building system, and in cold chain logistics distribution, additional heat of refrigeration equipment arises from such factors as heat exchange between internal and external walls of refrigeration equipment, heat diffusion of products and equipment in cold chain logistics, and hot air soaking. For the brevity of study, only heat load H 1 due to high-ratio heat exchange between internal and external walls of refrigeration equipment and heat load H 2 due to high-ratio heat exchange upon door opening are considered. H 1 is expressed as where R 1 is the thermal conductivity of refrigeration equipment in W/(m 2 ·K); S is the heat transfer area in m 2 of refrigeration equipment on a refrigerated truck; T is the temperature outside refrigeration equipment in K; and T 0 is the product storage temperature of cold chain logistics in K. H 2 is expressed as follows: where R 2 is the thermal conductivity of air in W/(m 2 ·K), S d represents the door area of refrigeration equipment in m 2 , and w 2 is a parameter of unit refrigeration cost; the refrigeration cost during cold chain logistics network distribution can be expressed as (4) Transportation Cost. When cold chain logistics products are distributed from a distribution centre to a customer, transportation cost will be generated, including management cost of distribution vehicle drivers and transportation staff, toll, and vehicle loss cost, while unit transportation cost per vehicle is assumed to be known and fixed in the logistics system. Transportation cost of all vehicles is related to transportation distance and calculated using the following formula: where φ is the cost per traveling distance.
(5) Penalty Cost of Time Window. As products in cold chain logistics are perishable, the customer demands that the logistics company distribute cold chain logistics products in a specified time frame. To improve logistics service level and customer satisfaction, penalty cost of time window is incorporated into cold chain logistics distribution routing objective function as follows: where P k (i) represents penalty cost of time window that may be paid when vehicle k is serving customer i; AT ik represents the time point at which vehicle k arrives at customer point i; α is a penalty coefficient for the case that a distribution vehicle arrives at a customer point at a time earlier than the start time of time window; β is a penalty coefficient for the case that a distribution vehicle arrives at a customer point at a time later than the start time of time window; ET i and LT i represent the start time and end time of time window, Computational Intelligence and Neuroscience 5 respectively. erefore, penalty cost of time window can be expressed as (6) Operation Cost of Distribution Centres. To ensure normal operation of cold chain logistics network, labour cost, equipment maintenance cost, and utilities have to be invested; the above costs altogether constitute operation cost of distribution centres in cold chain logistics. It is assumed that daily investment is basically unchanged in practical operation; thus distribution centres and operation costs are assumed to be fixed known conditions. Hence operation cost of distribution centres can be expressed as

Location-Routing Optimization Modelling for Cold
Chain Logistics. Based on the above model assumptions and definitions of notations and variables, a mathematical model was developed as follows: Equations (20) and (21) are objective functions. Objective (20) represents minimization of total cost during optimization of cold chain logistics LRP, including fuel consumption cost, cargo damage cost, transportation cost, refrigeration cost of refrigerated truck, penalty cost of time window, and operation cost of distribution centre; objective (21) represents the shortest vehicle distribution time. Constraint (22) ensures that each customer is served once; constraint (23) ensures that the vehicle has to leave after serving a customer; constraint (24) indicates that each customer can be assigned to one distribution centre only; constraints (25)- (27) indicate that the enabled distribution centre has to serve a customer; constraint (28) ensures that total demand of customers served by each distribution centre is not more than capacity of the distribution centre; constraint (29) ensures that total volume of cargoes carried by each vehicle is not more than its loading capacity; constraint (30) represents temporal connection between two consecutive customer points in a vehicle distribution path; and constraint (31) indicates that the time when a vehicle arrives at a customer point may not be out of the permissible time interval.

MOHH.
Conventional methods for solving MOP include weighting method, constraint method, and hybrid method. In recent years, Multiobjective Evolutionary Algorithms (MOEAs), as one kind of new methods for solving MOP, have been developed gradually. e most frequently used MOEAs are Genetic Algorithms (GAs), or multiobjective GAs, and representative GAs include VEGA, HLGA, MOGA, NPGA, and NSGA [25]. In this paper, MOHH is proposed for solving cold chain logistics-based multiobjective LRP model and distribution routing for cold chain logistics and is compared with Nondominated Sorting Genetic Algorithm II (NSGA-II) often used to solve MOP using elitist strategy [26]. NSGA-II employs simple yet efficient nondominated sorting mechanism so as to capture Pareto front and Pareto solution set with good distribution.
Problem domain and algorithm domain are designed below. e problem domain involves chromosome encoding and LLHs for optimization, while in high-level domain, four high-level selection strategies and acceptance mechanisms are designed; and finally, a multiobjective hyperheuristic framework is given.

Problem Domain
(1) Chromosome Encoding. Direct encoding is used. If there are three distribution centres and 10 customers, 0 represents the distribution centre, and 1-10 represent customer numbers; then a distribution centre can be represented by [11][12][13]. Customer sequence is generated randomly, for example, 9-5-3-2-4-1-6-8-7-10, and sequence location order represents the order of a vehicle arriving at a customer. ereafter, according to principles "customer demand served by a vehicle may not exceed load weight of the vehicle" and "customer demand served by a distribution centre may not exceed capacity of the distribution centre," distribution centres are allocated by centroid method; hence initialization path is converted to the following: 12-9-5-3-12; 11-2-4-1-6-11; 13-8-7-10-13. e first path represents that a vehicle departs from distribution centre no. 2 and distributes goods to customers nos. 9, 5, and 3 and then returns to the original distribution centre.
(2) LLHs. LLHs can be classified into mutation heuristic, ruin-recreate heuristic, local search heuristic, and crossover heuristic. LLHs are used to manipulate solution space directly as follows: (i) Mutation heuristic: LLH 1 : 2-opt: select arbitrarily one path from parent solution and swap locations of adjacent customer points so as to generate a new children solution LLH 2 : Or-opt: select arbitrarily one path from parent solution and insert two adjacent customer points into other locations so as to generate a new children solution LLH 3 : interchange: select randomly two paths from parent solution and swap locations of any two customer points so as to generate a new children solution LLH 4 : replace: select randomly one path and swap locations of any two customer points LLH 5 : shift: select randomly one path from parent solution while enabling a new distribution centre for this path so as to generate a new children solution LLH 6 : interchange: select randomly two paths from parent solution and swap their distribution centres so as to generate a new children solution (ii) Ruin-recreate heuristic: LLH 7 : location-based radial ruin: select any one customer point as "reference customer," and remove other customers at a probability of 1%-10% according to a rule of approaching "reference customer" in location, so as to generate a new children solution (iii) Local search heuristic: LLH 8 : interchange: as in the case of LLH 4 , but this operator returns the improved children solution only LLH 9 : shift: the operation is basically the same as that of LLH 3 , except that LLH 9 selects a customer point according to "the best improvement" criterion, inserts the customer point into random location of the original path or any other path, and returns an improved children solution only LLH 10 : 2-opt * : select randomly two paths from parent solution, choose one location according to "the best improvement" criterion, swap all customer points behind the location, and return an improved solution only LLH 11 : GENI: calculate the distance between any two customers on different paths in parent solution, select the shortest distance as reference distance, remove customer points nearest to the reference distance, rearrange these customers into a new path, and return an improved children solution only (iv) Crossover heuristic: LLH 12 : combine: select randomly two parent paths, copy one of the parent paths at a probability of 25%-75% to generate a subpath, add paths that originate from another parent path and are not contradictory against this subpath, and eventually arrange arbitrarily the remaining customer points LLH 13 : longest combine: select randomly two parent paths, sort all paths in descending order of number of customer points served, add all paths without repeated customers to generate a subpath, and eventually arrange arbitrarily the remaining customer points

High-Level Selection Strategies
(1) Simple Random (SR). SR selects randomly LLHs in each iteration and is usually used as a reference for comparison with any other strategy.
(2) Tabu Search (TS). TS is to prohibit repeating previous operation. To address a defect that local neighbourhood search algorithm gets readily trapped in local optimal points, TS algorithm formulates a tabu list to record searched local optimal points and not to search or search selectively elements in the tabu list in next iteration; thereby trap in local optimum is eliminated and global optimization objective is achieved. TS algorithm is an extension of local neighbourhood search and a heuristic for global neighbourhood search and gradual optimization.
As an HLH strategy for MOHH, TS scores performance of each LLH and thereby selects an operator for current iteration and updates scores using Reinforcement Learning mechanism; in other words, if a children solution arising from manipulation of the LLH improves parent solution, then score of this LLH will be added; otherwise it will be deducted.
(3) Choice Function (CF). CF selects an LLH based on three different metrics; the first metric records previous performance of each LLH, denoted as f 1 , and LLH performance of each LLH can be expressed as where I n (LLH j ) and T n (LLH j ) are the improvement value of each previously called heuristic and time to call the heuristic, respectively, and η is a coefficient ranging from 0 to 1.
Computational Intelligence and Neuroscience e second metric attempts to capture connection between LLHs and can be expressed using the following formula: where I n (LLH k , LLH j ) and T n (LLH k , LLH j ) represent FIR of LLH j following LLH k and call time of LLH j following LLH k , respectively. Likewise, μ is a coefficient ranging from 0 to 1. e third metric is the time period since CF selects the last LLH j : To rank LLHs, a score is given to each LLH: where e, f, and g represent weights of three metrics, respectively.
(4) Ant-Inspired Selection (AS) [27]. Ant colony optimization algorithm is a native heuristic of ant behaviour or ant population behaviour, used to solve hard combinatorial optimization problem. Ant colony algorithm applied in hyperheuristics is able to provide an HLH strategy so as to construct good LLH sequence. Each ant determines the next fixed point based on the pheromone in each path (each vertex represents an LLH), and once the ant arrives at a new vertex, it will apply the LLH corresponding to the vertex. Unlike the ant colony algorithm in combinatorial optimization problem solving mode, an ant colony manipulates LLHs; thus each ant is permitted to visit one point with multiple times, indicating that each heuristic may operate multiple times in one iteration. Details of the design are described in literature [27].

High-Level Acceptance
Mechanisms. Acceptance criteria are used to judge whether a children solution can replace parent solution or not, and its performance has immediate impact on convergence speed and optimization accuracy of a hyperheuristic [28]. It is therefore very important to determine superiority or inferiority of a children solution to parent solution in design of acceptance criteria. e following strategy is used in this paper: in a randomly selected objective function, as long as a children solution has a smaller fitness than parent solution, this children solution will be regarded to be superior to the parent solution and can go to the next iteration in place of the parent solution. is strategy is also used in the above three selection strategies. In a LRP model, the method focuses on designing distribution centres and customer groups, and the first job is to address distribution path between customers. e process is detailed as follows: Step 1: initialize the population by encoding customers directly: suppose each chromosome has 10 customer points and 3 distribution centres, 1-10 represent customer numbers, and 0 represents a distribution centre. Distribution centres can be represented by 0(1)-0(3), respectively.
Step 2: initialize parameters of an HLH strategy.
Step 3: optimize population P: Step 3.1: select one of two objectives (f 1 , f 2 ) at an equal probability randomly.
Step 3.2: select an LLH based on its performance.
Step 3.3: calculate the value of the selected objective function f n .
Step 3.4: calculate the value of another objective function.
Step 3.5: update performance metrics of the LLH.
Step 4: update Pareto solution set by nondominated sorting.
Step 5: evaluate a children solution according to an acceptance mechanism in acceptance criteria and opt to accept the children solution or not.
Step 6: update parameters related to HLH strategy and LLH score.
Step 7: judge whether termination condition is satisfied or not; if yes, then stop iteration and output optimal solution; otherwise, return to step 3.

Parameter Configurations.
e MOHH was coded in parallel in MATLAB 2015b using a 2.60 GHz Intel Core i5-3230K with 4 GB of RAM and running Windows 10.
e parameter configurations were set as follows. e size of population was set to 100; the maximum number of iterations was set to 5000. In the choice function, three 8 Computational Intelligence and Neuroscience weights were 0.5, 0.2, and 0.3. In the ant-inspired selection, we followed the default values proposed by Wang et al. [27]: length of route LP � 11; number of ants m � 15; factors ε � 0.001 and σ � 1.001.

Performance Metrics.
To validate performance of the proposed MOHH algorithm, the following three metrics were used: (1) Number of Pareto solutions (NPS): this metric is used to determine the number of Pareto solutions obtained by the algorithm. (2) Spacing metric (SM): this metric is used to determine distribution uniformity of Pareto solution set: where ‖x i t − y i t ‖ represents the Euclidean distance between nondominated solutions x i t and y i t .

Experimental Comparison of HLH Strategies.
To obtain better results in solving a biobjective LRP of cold chain logistics, MOHH algorithm was optimized by HLH strategy first; in other words, the best combination of selection strategy and acceptance criteria was found by comparing experimental results. One of Solomon reference instances was chosen for model solving, and as a standard LRP problem has neither refrigeration cost nor cargo damage cost, the total cost calculated in this section excludes refrigeration cost and cargo damage cost. Model solving results of HLH strategies are shown in Table 1.
Statistics NPS, SM, and DM are compared in Table 2. As can be found in Tables 1 and 2, AS-GD and TS-SA gave rise to more Pareto solutions than other combinations of HLH strategies; HLH strategy combination TS-SA resulted in smaller SM metric of nondominated solution set and smaller  Computational Intelligence and Neuroscience spacing, indicating that HLH strategy TS-SA results in more uniform distribution of nondominated solution set, as compared with other HLH strategies; in comparison, TS-GD, TS-OI, TS-SA, and AS-GD metrics had better DM values than other HLH strategy combinations. Figure 2 shows Pareto fronts obtained when selection strategies AS, CF, SR, and TS are combined with four different acceptance criteria, respectively. As can be found in the figure, when TS acts as a selection strategy, more Pareto solutions were obtained in more uniform distribution, and when SA acts as acceptance criteria, solution metrics NPS, SM, and DM were superior to other acceptance criteria.
In summary, when MOHH was applied to biobjective LRP optimization, the Pareto solution set obtained subject to selection strategy TS and acceptance criteria SA gave rise to good values of various evaluation metrics; thus TS-SA as a HLH strategy has better effect of solving biobjective LRPs.

Statistical Analysis of LLH Utilization.
Two new LLHs were added to improve solving biobjective problems by optimization, utilization data and acceptance rates of LLHs were statistically analysed during algorithm optimization, and performances of HLH strategies and LLHs in optimization were analysed in order to study the impact of MOHH algorithm framework on multiobjective optimization. Figure 3 shows the statistical results of LLH utilization. Figure 3 shows that, among all LLHs, LLH 6 had the highest utilization with a mean utilization of 11.3%; utilization of LLH 8 reached high values, too, with a mean utilization of 10.2%, whereas mean utilization of LLH 7 was merely 5.8%; except that LLHs were randomly selected in SR selection strategy, LLHs were selected by probability in all the remaining selection strategies, and the probability was determined based on performance of an LLH during algorithm iteration; it can therefore be inferred that LLH 6 and LLH 8 have better performance than other LLHs. In utilization statistics of LLH 6 , selection strategy TS showed the highest utilization of LLH 6 ; thus it can be inferred that TS enables better performance of LLH 6 , a LLH 7 utilization of just 4.6%, and lower utilization of worse-performance LLHs; therefore, in solving biobjective LRP, selection strategy TS is superior to three other selection strategies, which also validates experimental conclusions in this paper.

Statistical Analysis of Acceptance Rates in LLHs.
To study effects of acceptance criteria on solution results, acceptance rates of LLHs as per four acceptance criteria     Computational Intelligence and Neuroscience using selection strategy TS were statistically analysed (shown in Figure 4).
As can be found in Figure 4, according to acceptance criterion AM, all solutions are accepted, so all acceptance rates were 1; as per acceptance criteria GD and SA, goodperformance LLHs (LLH 6 and LLH 8 ) had mean acceptance rates above 0.9; according to acceptance criterion OI, only improved solutions are accepted and higher acceptance rates of high-quality LLHs were achieved, but acceptance rates of ordinary-performance LLHs were lower and readily trapped in local optimum. For poor-performance LLH (LLH 7 ), acceptance rates as per acceptance criteria SA and GD were 0.50 and 0.53, respectively. According to acceptance criterion SA, it is easier to reject poor-performance LLHs and accept good-performance LLHs; therefore, acceptance criterion SA has more positive impact on improvement of solutions to a biobjective LRP, which validates experimental conclusions in this paper.

Experimental Comparison between Algorithms MOHH
and NSGA-II. To further illustrate effectiveness of MOHH algorithm in solving biobjective LRPs, TS-SA serves as an HLH strategy combination of MOHH for instance solving, as compared with algorithm NSGA-II. For the sake of fairness, NSGA-II adopts the same heuristic that MOHH adopts (Table 3), while multiobjective performance metrics were used to assess the resulting Pareto fronts (Table 4).
In nine instances, MOHH had 5.2 Pareto solutions on average, while NSGA-II obtained 4.6 Pareto solutions. In terms of performance metrics, MOHH was slightly superior to NSGA-II in NPS and mean SM of MOHH-derived Pareto  SA. In cold chain logistics network programming, there are usually more than one decision objective; optimization of cost alone would fail to meet benefits of the logistics company and customers; therefore, it is essential to perform multiobjective optimization of cold chain logistics. is section focuses on the cold chain logistics of seafood of a supermarket in Hangzhou City, and the subbranches of the supermarket have a wide distribution and range of services. e distribution of subbranches (i.e., customers) is shown in Figure 5.
Aiming at improving competitiveness of cold chain logistics, the supermarket selected five distribution centres in Hangzhou. Five distribution centres met the requirements of the cold chain. In this paper, 30 subbranches are selected:    Tables 5 and 6.
Experimentally obtained Pareto solutions are shown in Figure 6. Solving biobjective cold chain logistics LRP by MOHH resulted in a total of 12 Pareto optimal solutions. Distribution time and costs corresponding to Pareto solutions were analysed in Table 7, where Cost 1 represents path cost, Cost 2 represents fuel cost, Cost 3 represents refrigeration cost of the refrigerated truck, Cost 4 represents cargo damage cost of cold chain logistics, Cost 5 represents penalty cost of time window, and Cost 6 represents operation cost of each distribution centre.
As described in Table 7, all the above 12 distribution schemes are Pareto optimal solutions. Pareto No. 1 corresponds to a solution with the shortest distribution time (4.2 h), but corresponding total cost reached 22126.0 Yuan; Pareto No. 12 corresponds to a solution with the minimum total cost (5522.0 Yuan), but its distribution time reached 7.5 h.
To validate the effect of MOHH in solving instances, LLH utilization during algorithm iteration was statistically analysed, as shown in Figure 7. Based on Figure 7, in selection strategy TS, utilization rates of LLH 6 and LLH 8 were      Computational Intelligence and Neuroscience the highest (13.1% and 12.1%, respectively), while utilization rates of LLH 3 and LLH 7 were only 5.1% and 5.5%, respectively. Moreover, according to acceptance criterion SA, acceptance rate of LLH 6 was more than 0.9, while acceptance rates of LLH 3 and LLH 7 were 0.53 and 0.52, respectively. According to statistical analysis in Section 3.3, LLH 6 and LLH 8 had the best optimization effect, while LLH 7 had worse performance. In process of solving cold chain logistics instances, in MOHH, LLH 6 and LLH 8 were still utilized and accepted at higher rates, while LLH 7 was utilized and accepted at lower rates; it can therefore be concluded that MOHH has better utilization of good-performance LLHs and better improvement of solution quality.
Practicality of MOHH was validated by instance analysis. For multiobjective LRPs, difference in objective preference will result in difference in distribution scheme and great differences in distribution centre selection and routing. Nonetheless, multiobjective optimization is capable of balancing between the logistics company and customer demand, and to meet different demands, different distribution schemes can be used to enable balance of interests between the logistics company and customers where possible. Furthermore, with total cost and distribution time as study objectives, it is in favour of decision maker of the logistics company in offering diverse path options.

Conclusions
In this paper, one biobjective mathematical model is established for cold chain logistics distribution system and is used to improve economic benefit, environmental benefit, timeliness, and customer satisfaction of the logistics distribution system. e first primary objective of the established multiobjective model is to optimize logistics distribution cost, consisting of path cost, fuel cost, refrigeration cost of the refrigerated truck, cargo damage cost of cold chain logistics, penalty cost of time window, and operation cost of distribution centres, while the second primary objective is to optimize timeliness of the distribution system or distribution time. To address this problem, MOHH algorithm was designed for model solving. As to HLH strategy for MOHH, the best HLH strategy combination was obtained from combinations of four selection strategies and four acceptance criteria.
Based on experimental results, TS-SA was the best HLH strategy. Further comparison with solution results of classic multiobjective solving algorithm NSGA-II validated effectiveness of MOHH, and, finally, solving instances of cold chain logistics by MOHH validated practicality of the algorithm, while the resulting Pareto solution set offers diverse options for a decision maker to select an appropriate distribution scheme depending on actual demand. e next study will focus on considering the use of more actual constraints in cold chain logistics such as heterogeneous vehicle [29] and simultaneous pickup and delivery [30], so that the problem can be more realistic. Higherperformance HLH strategy will be designed to improve solution quality and stability of the currently proposed MOHH. Moreover, we will focus on the alternative frameworks of MOHH like [31,32] and the framework that has a digital twin to repatriate at runtime [33]. In addition, consideration of other carbon emission models will become one of the focuses in the next work.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.