A novel bi-objective model of cold chain logistics considering location-routing decision and environmental effects

Economic, environmental, and social effects are the most dominating issues in cold chain logistics. The goal of this paper is to propose a cost-saving, energy-saving, and emission-reducing bi-objective model for the cold chain-based low-carbon location-routing problem. In the proposed model, the first objective (economic and environmental effects) is to minimize the total logistics costs consisting of costs of depots to open, renting vehicles, fuel consumption, and carbon emission, and the second one (social effect) is to reduce the damage of cargos, which could improve the client satisfaction. In the proposed model, a strategy is developed to meet the requirements of clients as to the demands on the types of cargos, that is, general cargos, refrigerated cargos, and frozen cargos. Since the proposed problem is NP-hard, we proposed a simple and efficient framework combining seven well-known multiobjective evolutionary algorithms (MOEAs). Furthermore, in the experiments, we first examined the effectiveness of the proposed framework by assessing the performance of seven MOEAs, and also verified the efficiency of the proposed model. Extensive experiments were carried out to investigate the effects of the proposed strategy and variants on depot capacity, hard time windows, and fleet composition on the performance indicators of Pareto fronts and cold chain logistics networks, such as fuel consumption, carbon emission, travel distance, travel time, and the total waiting time of vehicles.


Introduction
Logistics, which is a major contributor to carbon emissions (CE), pose challenges to global warming and climate change [1], especially in the context of road transportation. Wang et al. (2017) [2] stated that the CE produced by road transportation accounted for the entire transport sector CE by up to 70% which accounted for 14% in the global CE statistics. In response, according to the report of the United Nation Climate Conference in Copenhagen in 2009, China has promised that the CE would be reduced about 40%-45% of the amount of CE exhausted in 2005 by 2020 [3]. Due to the nature of low-temperature transportation, cold PLOS ONE PLOS ONE | https://doi.org/10.1371/journal.pone.0230867 April 9, 2020 1 / 29 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 1. Problem. In the proposed problem, several practical constraints were also considered: simultaneous pickup and delivery, hard time windows, and heterogeneous fleet, where the first two could which could increase client satisfaction and loyalty. Meanwhile, the types of cargos are classified into three modules: GC, RC, and FC. The former two can posit in one vehicle, but the latter must be separately stored in one vehicle. Moreover, the cold chain logistics were based on LRP.
2. Model. A novel bi-objective model was developed for the proposed problem considering three effects. The first objective can simultaneously concern two effects, which is more fit for practical logistics compared to the logistics using travel distance and time as travel costs. As the above mentioned, the second objective is modeled by the amount of damage in RC and FC (except for GC), aiming at improving the client satisfaction which could be further improved by another level: hard time windows.
To the best of our knowledge, the proposed model for the cold chain has not been studied thus far. The rest of this paper is structured as follows: Section 2 provides a review of related literature; Section 3 defines the formulation for the LRPLCCC considering simultaneous pickup and delivery, heterogeneous fleet, and hard time windows; Section 4 provides the solution method, solution representation, and search operators; Section 5 describes the computational experiments and simulated results; and, Section 6 outlines the study conclusions.

Literature review
Much attention and efforts have been drawn to the development of effective tools in supply chains and logistics systems [12], such as traveling salesman problem [13], vehicle-routing problem (VRP) [14,15], and LRP [16], where the latter is integrated logistics. Since the LRP was proposed by Jacobsen and Madsen (1980) [17] and Madsen (1983) [18], which deals with the combination of two types of decisions that often arise in logistics: the location of facilities and VRP [16], several LRP variants have appeared in the literature, such as capacitated LRP [16], LRP with simultaneous pickup and delivery [12], LRP with time windows [19], LRP with a heterogeneous fleet [19], etc. The works about LRP are plentiful with many exact and metaheuristic methods for tackling the variants of LRP.
The reason why LRP has been a research hotspot so far is its wide-range practicability and considerable difficulty. This paper investigates the application of LRP in the cold chain considering environmental effects and client satisfaction, hence, the following sections focus on the models of FCCE, low-carbon LRP (LCLRP), and cold chain with its variants based on VRP and LRP.

Models estimating fuel consumption
In the traditional models of logistics, travel distance and time were always used as the routing cost. However, in reality, the routing cost should be the costs of FCCE, which is much more practical. Hence, the models estimating the amount of FCCE should be developed for the logistics transportations. As stated by Demir et al. (2011) [20], measuring and reducing emissions requires good estimations to be fed into planning activities, which in turn require estimation models to be incorporated into the planning activities. Hence, it is necessary to analyze the effects of factors on the amount of FCCE. Demir et al. (2014) [21] provided the classification schemes for the factors: vehicle, environment, traffic, driver, and operations. However, factors of environment, traffic, and driver are quite difficult to obtain and evaluate. Hence, several researchers proposed some models using those factors that can be easily obtained, such as the models proposed by Xiao et al. (2012) [22] and Poonthalir and Nadarajan (2018) [23]. However, the accuracy estimating FCCE of the above models is unsatisfactory.
The above models estimating FCCE could be called as factor models, which are simple FCCE methods by using factors that can be easily obtained, such as load, distance, fuel efficiency, etc. Demir et al. (2014) [21] also summarized much more complicated models estimating FCCE: macro and micro models. The former uses average aggregate network parameters to estimate network-wide emission rates, which means that the factors affecting the FCCE keep constant, while the micro ones estimate the instantaneous vehicle FCCE rates at a more detailed level, which means that some factors could change in seconds like speed. The examples of macro ones are the methodology for calculating transportation emissions and energy consumption [24], the computer programme to calculate emissions from road transportation [25], etc. And the classical models in micro models are the comprehensive modal emission model (CMEM) [26], the models defined by Bowyer et al. (1985) [27], etc.
From the level of complexity, micro/macro views are complicated versions of macro/factor, and the above three are convertible, since some parameters in factor models and macro ones can be viewed as the specific parameters that combine the vehicle-related and average speed, and the dynamic parameters in micro models could be used as average values for simplifying the calculation of FCCE. Hence, the above three models could be classified into two modules: static and dynamic models. Moreover, in terms of estimating mechanism, they could be grouped into the load (power)-based and the regression-models. From the perspectives of estimated accuracy, the micro ones are the best, followed by macro ones, and factor model is the worst, the reason is that the estimation of FCCE depends on plentiful parameters which may change over traveling time. For a comprehensive view of the models and factors, the reader is referred to the surveys [21] and papers [20,[28][29][30][31]].

LRP considering environmental effects
With the implementation of energy conservation and emission reduction, more and more experts have indicated that logistics should not only consider economic indicators but also can reduce energy and emissions to provide a better ecological environment. Lin et al. (2014) [32] provides a comprehensive survey of VRPs considering FCCE. In their paper, a classification was proposed for various variants of VRP and solution methods. However, the models estimating FCCE were not reviewed, and they only listed and analyzed the papers about VRP considering FCCE. Moreover, as an extensive variant of VRP, the investigation of the green/lowcarbon LRP was not provided. Our previous work [29,30] provided a detailed assessment of the FCCE model for green/low carbon LRP, the solution approach, and the number of optimization objectives. This section provides additional papers on the recently published LCLRP, and other papers are also available in our previous paper [29,30].
In the last few months, several researchers have provided new sights in LCLRP. Dukkanci et al. (2019) [33] provided a single-objective model defined by logistics costs consisting fixed cost of depots and fuel consumption costs, while no fixed costs of vehicles were optimized. In their model of estimating fuel consumption, the micro view CMEM were used. However, no extra constraints were analyzed like the proposed ones in this paper. Koc (2019) [34] provided a single-objective model consisting of costs of depots, vehicles, and fuel consumptions, like their work published in 2016 [35]. Li et al. (2018) [36] proposed a many-objective model for a variant of LRP, that is, multi-depot green vehicle routing problem. In their paper, the objectives were revenue, logistics cost, traveling time, and carbon emission, where the latter used a simple estimation of carbon emission like the total traveling distance, hence, the reduction in carbon emission was not enough better than those using travel distance and time as traveling cost.

Cold chain logistics
Cold chain logistics is a low-temperature transport that helps improve the preservation of goods such as table grapes [37,38], salmon delivery [39], vaccine [40] and blood [41]. Therefore, the cold chain stream consumes more fuel and emits more emissions to maintain the quality and safety of perishable goods. Therefore, it is necessary to save energy and reduce carbon emissions in cold chain logistics to seek a win-win situation of economic and environmental sustainability [2].
Recently, the cold chain logistics considering FCCE has received increasing attention from research. Details are as follows: 1. Wang et al. (2017) [2] investigates the optimization of VRP with time windows for the lowcarbon cold-chain logistics based on carbon tax in China. In their paper, the objective consists of six parts: fixed costs of vehicles, transportation cost, refrigeration costs, penalty costs (soft time window), damage cost, and carbon emission. [44] proposed a low-carbon VRP-based model for the cold chain, which was very similar to the model proposed by Wang et al. (2017) [2]. The solution method in their paper is an acid-ant colony optimization algorithm.

Zhang et al. (2019)
The above four paper investigated the cold chain considering environmental effects. However, they have several common characteristics: (1) same FCCE model, that is, the factor model proposed by Xiao et al. (2012) [22]; (2) Refs. [2,42,44] applied the travel distance as the routing cost and used penalty costs (i.e., soft time windows), only the Ref. [43] used the FCCE costs as routing cost and did not use the penalty costs; (3) Swam heuristic methods. Refs. [2,42,43] designed the genetic algorithm and Refs. [44] applied an acid-ant colony optimization algorithm; (4) Refs. [2,43,44] defined the model for the cold chain logistics based on VRP, only the Ref. [42] is based on LRP.
However, in the reality, the routing cost is the FCCE cost like the models [34,35], and the factor model estimating FCCE might result in selecting the unpleasant routes with high pollution and high emissions since this model neglect the impacts of speed traveled over each arc and other parameters affecting FCCE. Furthermore, it is known through optimization experience that genetic algorithms and ant colony algorithms are not satisfactory for LRP and VRP. Since VRP and LRP are discrete combinatorial optimization problems, it is difficult to obtain the performance of the solution method by the ant colony algorithm and the discrete steps in the crossover and mutation processes.
The difference between this paper and the above papers are as follows: (1) model. This paper defined a novel bi-objective model for the LRPLCCC with minimizing the total logistics costs and the total amount of damage. The first objective consists of three parts: fixed costs of depots to open, fixed costs of vehicles to rent, and costs of FCCE used as routing cost, and the second objective combined hard time windows are used to improve the quality and safety of food product to improve the client satisfaction; (2) Constraints. The proposed model is limited by three practical constraints, which are first used in the cold chain. (3) Complexity. The above four papers only focus on a single type of cargos, but this paper studies the mixed cargo combined with integrated logistics and cold chain, which is the first study in the mixed logistics. (4) solution method. This paper applies local search and mutation procedures to solve the proposed model and develops an effective framework.
As for the cold chain only considering economic effects were studied by several researchers, and the reader is referred to the Refs. [45][46][47][48][49] based on VRP and Refs. [50,51] based on LRP, which are the latest papers after 2013, and the earlier works can be seen in the paper of Wang et al. (2018) [42].

Problem description
The problem is defined on a complete and directed graph O = (V, E), where V consists of N clients and M candidate depots, and the E is an edge set. Each client i2N have pickup demand p i and delivery demand d i , and hard time window and service time are (e i , l i ) and st i , respectively. Besides, the cargo types are ξ = {GC, RC, FC}. Each candidate depot j2M has a capacity CD j , a renting fee FD j , and a time windows DTW j . The clients are served by a heterogeneous fleet H = {h 1 , h 2 , h 3 }, and each vehicle type have corresponding parameters (see Section 5.2) which are used to estimate the amount of FCCE, including fixed cost FV h (h2H) and a capacity CV h . In terms of edge set, there are a travel distance D ij and speed limitation SP ij for each edge (i, j)2E = {(i, j): i, j2V, i6 ¼j}\{(i, j): i, j2M}. The goal is to determine the set of depots to open and the tracing of the routes to minimize the total logistics cost and the amount of damage on the quality of the cargos. Before defining our model, several assumptions should be described: 1. Each client can only be served once by one vehicle and depot; 2. Each client has the same type of delivery and delivery demands; 3. The GC and RC can be delivered or picked up by the same vehicle for the clients, but the vehicle state must adaptively change the state of the refrigerator according to the type of cargo on each edge, and the FC must be supplied in the vehicle in the frozen state; 4. Each vehicle must return to the original depot before the depot closes; 5. The load of each vehicle on each edge must be less than its capacity; 6. The depot must serve the clients assigned to it; 7. The load of each depot must be less than its capacity; 8. The vehicles travel on each edge with the known speed limitation without considering the road and traffic conditions; 9. The vehicle must wait until the moment reaches the opening time window of each client if it arrives early.

Model development
This section provides the model estimating FCCE and a variable function of refrigerated goods quality.

Model estimating FCCE
This section presents a microscopic model to estimate the amount of FCCE, namely CMEM, which is easily applicable to calculate FCCE, and it has been extensively used in LCLRP [28-30, 34, 35]. This paper applies the version of CMEM which is used in the cases with a heterogeneous fleet. The fuel consumption rate FCR h (g/s) of a vehicle type h2H is given by where φ is the fuel-to-air mass ratio; k h is engine fiction parameter (kJ/rev/L); N h is engine speed (rev/s); T h is engine displacement (L); P h is the second-by-second engine power putout (in kW); η is an efficiency parameter for diesel engines; and κ is the heating value of a typical diesel fuel (kJ/g).
where P h tract represents the total tractive power putout (kW); η tf is the vehicle drive train efficiency; P h rc is the engine power demand associated with running losses of the engine and the operation of vehicle accessories such as air conditions and refrigerating compressor in the cold chain logistics. The estimation of P h tract is given by: where G is the sum of the total weight of the vehicle weight (w h ) and the vehicle load (kg); a and g are acceleration of vehicle and gravitation (m/s 2 ), respectively; θ is the road angle; C h d and C r are, respectively, coefficient of aerodynamic drag of a type h2H and coefficient of rolling resistance; ρ/A h is the air density (kg/m 3 )/the front surface area (m 2 ); s is the instantaneous traveling speed (m/s). Then the fuel consumption F h,g (g) of vehicle type h over traveling a time t from t 1 to t 2 at instantaneous traveling speed s(t) is calculated as where λ = φ/κψ, γ = 1/(1000×n tf η), w = a+gsinθ+gC r cosθ and β v = 0.5C h d ρA h , and ψ is conversion factor which converts F h,g in gram to F h,L in liter. The above microscopic equation is a basic estimate of FCCE and consists of three modules: (1) engine module, linearly proportional to the travel time; (2) moment module, linearly proportional to the product of vehicle weight and travel distance; (3) speed module, linearly proportional to the integral of s(t) 3 . If a vehicle of type h travels a distance d at a constant speed s, then Eq (6) can be rewritten as: Eq (7) is the general model estimating the fuel consumption during s6 ¼0. However, for RC and FC, the vehicles must wait for the clients until the opening time window is open. Hence, the refrigerating compressor must run for keeping the quality and safety of food production during the waiting time of vehicles and service time for RC and FC. Therefore, for the cases in this paper, the F ijh over arc (i, j)2E can be written as where dv1 ij is the state variant of the refrigerating compressor depending on type of cargos (see Eq 9), indicating that each type of vehicle has three states.
AT ih is the arriving moment on the node i2V of a vehicle type h; L ijh is the load of vehicle type h over arc (i, j)2E; dv2 ij is a variant equaling to 1 if the RC/FC is a part of L ijh and to 0 otherwise.
x ijh is a decision variant which equals to 1 if a vehicle of type h2H travels on edge (i, j)2E and to 0 otherwise.
The assumption H3 is imposed by dv1 ij ×dv1 jk 2{1,2,4,9} if the clients i, j, and k are served by the same vehicle. Moreover, as reported by Refs. [28][29][30]34,35], the exhausted CE (in kilogram) is estimated through the amount of fuel consumption, in other words, the amount of carbon emission is directly proportional to the amount of fuel consumption, namely, 1 liter of fuel can produce 2.32 kilogram of CO 2 . Hence, the cost of FCCE is as follows: where c fc and c cc are, respectively, the price of 1-L fuel and 1-kg CE.

PLOS ONE
This paper applies the variable function of the quality of refrigerated goods used in [2,[42][43][44]: D(t) = D 0 e -at , where D(t) is the quality of the cargo at time t; t is the transportation time; D 0 is the quality of the product from the depot; and a is, related to the characteristics of cargos and temperature, the spoilage rate of the product. Combined with the characteristics of this paper, we used the following equations to calculate the amount of damage on the quality: where D1 ijh is the damage when the vehicle h2H leaves node i2V and arrives at node j2V and the door is not opened in the process of transportation; D2 ijh is the damage when the vehicle h2H serves client j2V (the door is open); L � ijh is the RC/FC weight of the vehicle h2H over arc (i, j)2E, if the total load belongs to GC, its value equals to 0; a 1,lx and a 2,lx are the spoilage rates for the cargos type lx2{GC, RC, FC}, respectively; Let r equal to 0 if the cargos type of client j2V is GC and to 1 otherwise. Hence, the total damage on the quality of cargos are as

Establishment of the formulation
Section 3.2 provides a model for estimating FCCE and a method for calculating quality damage. Therefore, the formal model of LRPLCCC can be defined as follows: where y i be equal to 1 if a depot i2M is selected and to 0 otherwise. Objective (14) is the total cost of three parts: fixed costs of the open depots, the fixed costs of the leased vehicles, and the FCCE costs. Objective (15) is the total damage of RC and FC quality.
The following constraints are primarily used to satisfy and guarantee the above hypotheses in Section 3.1. The following two are degree restrictions. In particular, constraint (16) states that each customer must serve only once; constraint (17) provides a balance between entering arcs and leaving arcs.
The following three ensure that a load on the depot must not exceed its capacity and the additional restrictions on the depot load at the beginning and end of the service.
where z ij indicates the assignment of the client i2N, if it is assigned to the depot j2M, then z ij = 1, otherwise z ij = 0. Corresponding constraints on a load of vehicles on each edge are as follows: In detail, constraint (21) imposes that the load on each edge must be less than the capacity of the assigned vehicle; constraint (22) is the bound on the load variables; constraint (23) and (24) are the restrictions on the load variables at the starting and finishing stage; constraint (25) implies that the pickup and delivery demand of each client is met.
The constraints on hard time window for each client and each depot are as follows: Eq (26) is the continuity of the vehicle travel time. Constraint (27) stipulates that the arrival time of each vehicle must not exceed the closing time window of each client. Constraint (28) limits each vehicle to return before the closing time of each depot. The following three prohibits the formation of routes that do not start and end in the same depot: Constraints (32) and (33) guarantee that each client must be served by only one vehicle and depot to open: X Constraints (34) and (35) make sure that the depot to open must serve at least one client: In particular, constraint (34) indicates that the unselected depots must not serve any of the clients; constraint (35) states that the depot to be opened must serve at least one client.

Extra valid restriction
Constraints (16)-(35) are the key and necessary for the proposed problem, and this paper also provides several polynomial-size, valid, not necessary, and alternative inequalities described below. First, the subtour must not exist, that is, subtour elimination, which could be seen as a complementary one for constraints (16) and (17): Constraint (36) is special case of classical subtour elimination constraints used by   [19]. The following inequalities are used to restrict the depots to open and vehicles: X In particular, constraint (37) indicates that the unselected depot must not assign a vehicle; constraint (38) states that the depots to be opened must allocate vehicles to serve the clients. The above can be used as a supplementary restriction for constraints (16), (17), (34) and (35), since constraints (16), (17), (34) and (35) mandate that each client must serve only once by a vehicle and opened depot, indicating that the open depot must allocate the vehicle.
Constraint (39) is the bound on the number of the assigned vehicles, including upper bound and lower bound.
where d•e is the smallest integer larger than •. The next inequality imposes the number of depots to open: However, constraint (40) is not always feasible, since it also depends on the number of clients, the example is the instance with 55 clients and 15 depots introduced by Barreto et al. (2007) [52]. The final inequality is a complementary restriction for constraint (33) to forbid the different depots in a single route. X The above constraints (16)-(41) could also be used in our previous works [29][30] for the another variant of LCLRP except for the constraint (28).

Proposed methods
Since this paper tackles a bi-objective model for the LRPLCCC, the corresponding MOEAs were used to obtain the Pareto solutions. In seven MOEAs, a practical framework was developed to generate the children solutions which is fit to the local search and mutational procedures (Section 4.1). In the proposed framework, a simple and effective chromosome representation for the proposed LRPLCCC were developed (Section 4.2). Moreover, we designed 14 neighborhood operators (Section 4.3). The details are provided as follows.

Framework used in MOEAs for the LRPLCCC
An overview of the pseudocode for the proposed framework used in MOEAs is given in Algorithm 1 (Fig 1). Input data includes the maximum number of iterations (T max ), mutation probability (p m ), and initial population (Pop) made of feasible random chromosomes (Section 4.2). The output data is an elite population (EP) that consists of non-dominated solutions by removing duplicates and dominating solutions.
Step 2 is the parameter settings in the seven MOEAs, such as the grid density in GrEA, the scaling parameter in the IBEA, and the reference points (or uniform points) in the NSGA-III. Then the main loop is performed, stopping when a maximum number T max of iterations.
In each iteration of the main loop, the evolutionary phase is first performed in Steps 6-12. In other words, for each individual in the current population, a mutational operator is randomly selected from the pool of perturbation to provide slight randomness, depending on the p m value, then a local search procedure is also chosen from the pool of local search procedures to improve the quality of the obtained solution.
After improving all individuals in the current population, the children population (CP) is merged with the parent population (i.e., the current population) to obtain the next population with N p individuals by utilizing the update mechanism used in seven MOEAs (Steps 13-14). Moreover, this framework also provides the external archive to save the best non-dominated solutions with 5×N p individuals which is the data outputted (Steps [16][17].
When the main loop stops, the algorithm ends and returns to EP.

Chromosome representation
Since LRPLCCC is one of the discrete combinatorial optimization problems, a simple and efficient representation is used to represent the solution, which was also used in our previous paper [12,[28][29][30]. In the construction of solutions, each route (corresponding to the vehicle route) consists of a sequence of clients and depots inserted at both ends. Hence, a complete solution is represented by R = {r 1 , r 2 , . . ., r k }, where r i is a complete route. Besides, we also save several attributes, such as the load, type, and state of vehicles and two objective values, aiming at allowing a fast calculation in the implementation process in the operators and objective space. Besides, the adopted representation, together with the operators (detailed in Section 4.2), could provide feasible children solutions by meting constraints (16)- (41).

PLOS ONE
listed in the right portion. Besides, the initial population is randomly generated, which satisfies all constraints (16)-(41).

Neighborhood search mechanisms
In the proposed framework, 14 operators are designed to obtain the Pareto solutions, which can be classified into three modules: dominated operators (DO), Non-dominated operators (NDO), and mutational operators (MO). The execution mechanism of DOs is to obtain the children solutions dominated the parent individuals, which can accelerate convergence in the early stages of the algorithm. Since EP consists of non-dominated solutions, NDOs are used to obtain a large number of non-dominated solutions. While often insufficient to achieve competitive results, MOs can be used to provide randomization when searching globally for performing simple random moves. In this paper, the DO and NDO pools consist of five operators: 2-opt, swap, insert, segment-based swap, and segment-based insert. The MO pool consists of four operators: add, decompose, insert +, and swap+ moves. Details are as follows.
The 2-opt move is executed by removing two edges from two routes and reconnecting the new four paths created (see Refs. 52). The swap move swaps the location of two clients from different routes. The insert move is implemented by moving a client into a position of another route. The segment-based swap and segment-based insert moves like swap and insert moves, but the objects of them are two or three clients instead of one client.
The schematic diagrams of the above three moves are provided in Fig 3. As to the schematic diagrams of insert+ and swap+ moves also like swap and insert moves, but the output of these moves usually are dominated by the parent solutions.
The above operators are implemented by a complete process [53] rather than random selection, avoiding costing CPU time. Moreover, considering the hard time windows, the schematic diagrams of the above 12 operators implement the moves between different routes rather than those inside the routes.
As to the add move, it seeks to avoid a fast convergence to solutions with few depots (prone to happen due to 2-opt, insert, and segment-based insert moves). Meanwhile, it diversifies the opened depots by opening a new one and randomly reassigning between 1 and 2/3 of the routes to it [16] or randomly choosing a depot from the set of all depots to open after closing one opened depot.
The decompose move is executed by decomposing one route into two routes, and the breakpoint is randomly selected only if it has more than one clients to serve. The benefit is to avoid a fast convergence since the 2-opt, insert, and segment-based insert moves could easily result in the long tracing of the route using vehicles with higher capacity.
Since a forthcoming article (2020) has tested and analyzed the performance of the first found and the best improvement, the results showed that the performance of the operators using the first found improvement is superior to those using the best improvement. Hence, this paper used first found improvements as the mechanism in operators. In other words, if a better solution is found, the operator will stop and return to the child solution. In the DO pool, the definition of "better solution" is that the child solution dominates the parent solution, but in the NDO pool, it is defined that the parent solution cannot be dominated by the child.

Optimization simulation and analysis of results
Implementation aspects and evaluations of the proposed problem and algorithms are provided and discussed in the following sections.

Implementation aspects and parameters configurations
The proposed algorithms and problem were coded in MATLAB and results were obtained using 4.0 GHz Intel Core i7-6700K CPU with 12 GB of RAM and running Windows 10.
Parameters configuration plays an important role in affecting the performance of the proposed problem and algorithms. However, this paper favors the default values which has tuned by the published papers, such as the scaling factor in GrEA (0.05) [9] and the number of the divisions of the objective space in each dimension (45) [10].
The size of the initial population N p is 100, just like the traditional setting of the bi-objective problem. And the maximum number of iterations, T max , is directly dependent on the number of nodes (i.e., clients and depots): For all instances, the multiplier α is set to 10. As for the size of the external archive, we return 5×N p elite individuals to prevent the loss of non-dominated solutions during the search process. For the mutational probability p m , we conducted an initial experiment with different values in the interval [0, 1] combined with seven MOEAs. The reason can be drawn that each MOEA may favor different p m . And the results for the initial experiment were presented in Section 5.4.

Test instances
Since the LRPLCCC is first studied in this paper, the instances proposed by the existing papers are unable to be reused. Therefore, we randomly generated 28 instances with different number of clients and depots, called as C|N|-|M|-No., where |N| and |M| are, respectively, the number of clients and depots, and No. is the serial number of instances. The number of clients is |N|2 {20,30,40,50,60} and the number of candidate depots is |M|2{4,5,6,7,8}. All nodes are randomly located in the square interval [0, 50] 2 km. The delivery and pickup demand of each client is generated using a uniform distribution in the range [100, 1600] kg, and the opening time window e i (i2N) of each client is obtained from the instance C101 [54], which decreased by 50%, and the closing time window l i (i2N) for each client depends on the service time (i.e., l i = e i +st i ): The capacity and fixed cost of each candidate depot are, respectively, generated using a uniform distribution in the range [10,15] tons and [500, 1000] yuan. And the closing time window for all depots is set to 12 hours.
Regarding the parameters of the heterogonous fleet, they are given in Tables 1 and 2, which provide vehicle-specific parameters and general parameters, including fixed costs, freezer power for each type of vehicle, rate of corruption, etc.

Performance metrics
To validate the reliability of the proposed algorithms, the four well-known performance metrics, namely generational distance (GD), inverted generational distance (IGD), hypervolume (HV), and the ratio of non-dominated solutions (RNI), are applied.
GD describes the quality of an approximate Pareto front (APF) by measuring the distance between APF and real Pareto front (RPF), the smaller the GD value, the better the quality. The GD equals to 0 indicates that APF is a part of RPF.
The IGD describes the quality and uniformity of APF by measuring the distance between RPF and APF. The smaller the IGD value, the better the distribution and convergence. The IGD equals to 0 indicates that the APF equals to RPF.
The HV measures the size of the coverage space between the APF and the reference point. The larger the HV value, the better the diversity and distribution. HV is in the range [0, 1].
RNI describes the contribution rate of the APF for the RPF. The larger the RNI value is, the better the APF. The RNI equals to 1 indicates that the APF equals to RPF.
However, in our case RPF is unknown. For this reason, and following the paper [55], this set RPF was made of all sets APF obtained by all MOEAs by removing dominated and repeated individuals. The set RPF is in fact an approximation of the real Pareto front.

Parameter turning (p m ) and analyses on the pairs of MOEAs
The mutational probability p m has significant impacts on the performance of MOEAs. For the proposed seven MOEAs (BiGE, GrEA, IBEA, NSGAIII, NSGAII, NSLS, and SPEA2), a situation may exist: a special pair of each algorithm and p m value could maximize performance in obtaining the Pareto solutions. Hence, we evaluated the effects of the p m values on the performance of seven MOEAs. The instances used in this section were the C20-4-no., C30-5-no., and C40-6-no., where no.2 {1,2,3,4}, i.e., the total 12 instances. Each MOEA performed four runs on each instance. Moreover, we tested the impacts of the p m values in the range [0,1], i.e., p m 2 {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1}. Therefore, a total of 3696 runs were performed. To rank MOEAs, we ranked 77 pairs using the scoring system of the CHESC-crossdomain heuristic search challenge (http://www.asap.cs.nott.ac.uk/external/chesc2011/). In this system, the top eight pairs score 10, 8, 6, 5, 4, 3, 2 and 1 respectively. The media used are performance metrics (HV, IGD, and GD) obtained from the average of four runs. Therefore, the highest score for each pair is 120 (12 instances). The results are shown in Table 3. Looking at Table 3, the first two pairs with the best HV, IGD, and GD are NSGA-II (0.7), NSGA-III (0.9), NSGA-II (0.6), SPEA2 (0.9), BiGE (0.1), and GrEA (0.1). They were used in the following experiments to analyze the effects of problem parameters on performance parameters of the Pareto frontier and the cold chain network. We could find that: (1)

PLOS ONE
Bi-objective model of cold chain based location-routing decision and environmental effects best performance in terms of GD values. The reason is that the fewer mutation operators are used, the higher the intensification of the algorithm, but the performance of MOEAs using 0 as the p m value will seriously deteriorate the performance of the algorithm. Therefore, it is very important to use the perturbation operator in the search process. Aiming at analyzing the effects of problem parameters on the performance indicators of Pareto fronts and cold chain logistics, we used the top two MOEAs in terms of average HV, IGD, and GD values, that is, NSGA-II (0.7), NSGA-III (0.9), NSGA-II (0.6), SPEA2 (0.9), BiGE (0.1), and GrEA (0.1).

The effect of the cargo type
In this paper, the proposed strategy (S1) is assumed: the load of several vehicles may be mixed cargos (GC and RC), and the FC must be served separately, i.e., dv1 ij ×dv1 jk 2 {1,2,4,9}. At the same time, there is another version (S2): the type of cargos loaded in each vehicle must keep the same, in other words, each vehicle only loads one type of cargos, i.e., dv1 ij ×dv1 jk 2 {1,4,9}. To analyze these two strategies, this section uses the six MOEAs above to optimize instances used S2, which are also used in Section 5.4. The results are shown in Table 4.
As shown in Table 4, the performance of S1 is much better than S2, especially C20-4-2, C20-4-3, and C20-4-4. S1 could help increase the HV value by 1.75%, reduce the GD value by 91.62% and the IGD value of 99.46%. Based on the RNI values, S1 and S2 will provide 98% and 0.09% non-dominated solutions for RPF, respectively. The proposed strategy could obtain better Pareto fronts than S2, as shown in Fig 5, and it can be concluded that most of the solutions obtained by S1 dominate the solutions obtained by S2. Therefore, S1 will help provide decision makers with better solutions to choose the right solutions. Although the positions of Pareto fronts obtained by the two strategies are very close, the RNI indicator can provide the proportion of non-dominant individuals in the construction of RPF. Considering the length of this article, we also provided a partial enlargement of the Pareto frontier in the S1 File and analyzed the impacts of two strategies on CE/fuel consumption, vehicle travel distance, travel time, and total waiting time (VWT). And we could draw the following conclusions for most instances: (1) S1 seems to produce more but not large CE than S2; (2) The travel distance of S1 is less than the travel distance of S2; (3) The travel time of S1 is much lower than the travel time of S2; (4) The VWT of S1 is also much lower than the VWT of S2. Therefore, we conclude that for LRPLCCC, S1 is better than S2.

Efficiency of the proposed model
Since this paper first studied the LRPLCCC, the proposed model utilizing FCCE costs as routing cost should be analyzed and compared with the traditional models using travel distance (TD.) and travel time (TT.) as routing cost. In the models with TD. and TT., we assume that prices per kilometer and minute of both models are 5 Yuan/km and 2.5 Yuan/min, respectively. After the Pareto front of each instance was obtained by the models with TD. and TT., we used the proposed model to recalculate two objective values for comparing the Pareto fronts. Table 5 is the performance indicators for the three models, where all individuals (including the dominated solution) obtained by models using TD and TT were used to recalculate the Pareto fronts under FCCE. In addition, the non-dominated solution obtained by the models with TD and TT are also used to recalculate the Pareto fronts under FCCE, the corresponding data and figures are shown in the S1 File. Looking at Table 5, for the HV values, the increment of the Pareto front obtained by the model using the FCCE averaged 14.05% and 5.57% compared to the model with TD. and TT. For the GD and IGD values, the proposed model can reduce the GD values of 96.61% and 96.29%, and the IGD values of 99.43% and 98.86%, respectively, compared with the model with TD and TT. For the RNI values, the models with TD and TT could reduce by 98.41% and 97.63% when compared to the RNI obtained by the proposed model. Hence, the proposed model in this paper is efficient in terms of Pareto solutions. Fig 6 is the Pareto fronts obtained by three models.
From the Pareto fronts shown in Fig 6, we could conclude that the Pareto fronts of the proposed model in this paper could dominate most solutions obtained by the models with TD and TT, especially when the logistics costs are less than 0.5 the maximum costs. Moreover, the Pareto fronts obtained by the model using TT are better than those obtained by the model using TD. Besides, we also analyzed the effects of three models on the CE/fuel consumption, travel distance, and travel time in the S1 File. As the results showed, the following can be obtained for the most instances: (1) The model using FCCE costs as routing cost could produce the minimum CEs among three models, followed by the model with TT; (2) The model using TD could serve clients within the minimum travel distance, which matches its nature, followed by our proposed model; (3) The model using TT could serve the clients within the

PLOS ONE
minimum travel time and VWT, which matches its nature, followed by our proposed model. Hence, the proposed model using the FCCE cost as routing cost could improve economic, environmental, and social effects.

The effect of depot capacity
In the literature, some studies on LRP solved the instances of an unlimited capacitated depots [56]. We now analyze the impact of depot capacity variants on the Pareto frontier. To this end, we have increased the capacity of candidate depots by five times (+ 10, + 20, + 30, + 40 and + 50%) and a variant with unlimited depot capacity. This section applied eight newly generated instances with 50 and 60 clients. Table 6 lists the performance indicators for the eight cases. Looking at Table 6, in terms of HV values, the HV values increase as the depot capacity increases. As shown in Fig 7(A), the increments reach 1.23, 1.70, 2.71, 4.05, 4.51, and 7.61% on average when compared to the HV values of the base instances. As to the average IGD values, the IGD values decrease as the depot capacity increases. As shown in Fig 7(B), the decrements reach 23.74%, 34.76%, 51.39%, 70.54%, 77.60%, and 100% on average when compared to the IGD values of the base instances. From the perspectives of average RNI values, the RNI values increase as the depot capacity increases, as shown in Fig 7(C). Fig 8 is the Pareto fronts of seven variants of the eight instances. We could find that the Pareto fronts of instances with unlimited depot capacity will dominate others. However, the differences between the seven variants are not particularly large. The reason could be hard time windows. The hard time windows for clients and depots limit the LRPLCCC logistics network, so while the depot capacity increases, there are some clients that cannot be served by some depots.
In the S1 File, we also analyzed and compared the effects of the depot capacity on the FCCE, travel distance, travel time, VWT, and total costs of depots to open. As the results of most instances showed that: (1) The exhausted CEs, travel distances, and travel time increase as the depot capacity increases; (2) The total waiting time of vehicles and costs of depots to open decrease as the depot capacity increases. The reason is the higher the depot capacity, the more clients a depot serves, resulting in that fewer depots can be opened and vehicle travel distance, CE, and travel time are increased. While the vehicle waiting time is related to the client time windows, the results showed that the time windows match the cases with unlimited depot capacity.
As stated as Koc et al. (2019) [34], in practice, for some logistics enterprise and their logistics items, the depot capacity is not very important, and our results of the instances with unlimited depot capacity showed that Pareto fronts have better HV, IGD, and RNI values.

The effect of hard time windows
Section 5.7 has analyzed the effects of the depot capacity on the Pareto fronts and several performance indicators of the LRPLCCC, and we could obtain that the hard time windows of clients impact the selection of the sets of the depots to open and the tracings of the routes. Hence, this section is provided to analyze the effects of the time windows of clients on the Pareto fronts and other characteristics of the LRPLCCC. Again, we also generated six variants of the base cases by changing the time window margins like our previous work [28], i.e., l i = e i +(1+δ) ×st i (i2N) where δ2 {10%, 20%, 30%, 40%, 50%}, and a special case l i = 1 (i.e., δ = 1). Table 7 presents the results on performance indicators of seven cases, i.e., HV, IGD, and RNI. Looking at Table 7, we could obtain that: (1) The HV and RNI values of eight cases increase as the δ values increase; (2) The IGD values of eight instances reduce as the δ values increase.  28.40, 35.47, 41.68, 49.33, and 100% of IGD values. In terms of RNI values, six δ values could provide 0, 2.69, 3.62, 4.85, 8.24, 11.56, 14.29, and 100% non-dominated solutions for the RPF. Hence, we could conclude that the hard time windows of clients have significant impacts on the performance indicators of the Pareto fronts. Fig 10 provides the Pareto fronts of seven variants of each instances.
As shown in Fig 9, the Pareto fronts of the instances using δ = 1 values dominate most of the solutions of the other cases. As to the others, we could find that there are small (but not large) gaps among the Pareto fronts of the cases using the other six δ values. However, from the results in Table 7, larger δ values could obtain a better Pareto front, but the range of variation gradually approaches a constant like the results in Section 5.7. Moreover, we could find that the two poles of the Pareto front are extended with the increase of δ values, especially when δ = 1. The reason is that the larger closing time window of clients allows much more selection of routings since the closing time window is a strong constraint.
We also analyzed the impacts of the hard time windows on CE/fuel consumption, travel distance, and travel time in the S1 File. Let the second objective be the abscissa, we found that when the damage quality is small (that is, the total cost is large), the above three curves almost coincide, but when the damage quality is large, the curve will change greatly (i.e., the total cost is small), especially CE/fuel consumption, travel distance, and travel time is the smallest in Moreover, a novel conclusion could be reached by analyzed the Pareto fronts of this section and Section 5.7: The Pareto front could rotate around a fixed point. This behavior is not shared by most other problems and methods in the literature.

The effect of fleet composition
This section analyzes the benefits of using a heterogeneous fleet of vehicles over a homogenous one. To this end, we have conducted four sets of experiments on the newly generated 50-and   Table 8.
Looking at Table 8, in terms of the HV values, the instances using HF could obtain the highest values, which provide more 7.66, 3.56, and 18.63% of the HV values compared to the others. From the IGD and GD values, the instances with HF could obtain the best values with all indicators of 0. Fig 11 is the Pareto fronts of eight instances, and the Pareto fronts of the instances using HF could dominate most solutions of other variants. Moreover, the Pareto fronts of the instances using H1 are close to those of instances with HF when the total logistics costs are high, this could obtain from the GD values which are much lower than others. Hence, the HF could provide much better solutions for the decision-makers. Table 9 is the composition of the vehicle types used in the HF case. We could find that the most commonly used fleet consists of H1 and H2 in the non-dominated solutions of the HF cases, but the combination of H2 and H3 does not provide any non-dominated solutions for the Pareto fronts. Therefore, the instances in this section favor H1 and H2 except for the H2 and H3, other combinations could provide more or less non-dominated solutions. The best fleet composition depends on many factors, such as delivery and pick-up demands and hard time for customers.  In the S1 File, we also analyzed the impacts of the fleet composition on CE/fuel consumption, driving distance, travel time, and total vehicle waiting time. In most instances, the results showed that: (1) The HF in will achieve the lowest CE/fuel consumption, while the use of H3 will bring the highest CE/fuel consumption in the same damage. The reason is that HF's strategy could choose the best types of vehicles for LLRLCCC based on FCCE and vehicle load and total cost. (2) The use of H3 can complete the task at the minimum travel distance, travel time and VWT, which can be explained by the fact that the larger the capacity of the vehicle, the more the vehicle can serve more clients.
Moreover, the limitation of this section is that the fleet composition depends on the nature of the instance, such as vehicle load, hard time window, etc. In other words, if the instance meets the specific requirements, any fleet composition can get the most benefit. However, the generated instances using HF could provide the best Pareto fronts in this paper.
In short, the above sections analyzed efficiency of the proposed strategy and model, and analyzed the effects of the problem parameters such as depot capacity, hard time windows, and fleet composition on the performance indicators of Pareto fronts and LRPLCCC logistics network, such as CE/fuel consumption, travel distance, travel time, and the total waiting time of vehicles. Hence, the logistics enterprises should also analyze the problem parameters on key performance indicators of the logistics network, aiming at improving the sustainable development in economy, environment, and society.

Concluding remarks
In this work, a novel bi-objective mathematical model for cold chain-based low-carbon location-routing problem was developed. In the proposed model, the first objective consisted of three parts: fixed costs of depots to open, fixed costs of renting vehicles, and the total routing costs, where the latter can be defined with respect to the cost of fuel consumption and carbon emission, and the second objective was used to minimize the total damage on the quality of cargos. Therefore, the bi-objective model can improve economic, environmental and social benefits. This paper also proposes a strategy for mixing cargos, which are the clients' needs. Besides, several practical constraints were considered in the proposed model: simultaneous pickup and delivery, hard time windows, and a heterogeneous fleet. In the proposed method, we proposed a simple and effective framework inserted in seven well-known multi-objective evolutionary algorithms to solve the proposed model.
In the experiments, we first evaluated the effects of mutational probability in the seven MOEAs, and we determined the top two MOEAs in terms of average HV, IGD, and GD values to optimize the rest experiments. Then, we examined the efficiency of the proposed strategy and model. Extensive analyses are performed to empirically assess the effect of various problem parameters, such as depot capacity, hard time windows, and fleet composition on key performance indicators of Pareto fronts and LRPLCCC network, including fuel consumption, carbon emission, travel distance, travel time, and the total waiting time of vehicles.
This study has some limitations. The model presented in this paper was based on the assumption that the type of cargos for each client remains the same. However, in real-world applications, the client's cargo types could be mixed, i.e., general, refrigerated, and frozen cargos. Hence, the future works may focus on the cold chain logistics considering environmental effects that the types of cargos of each client are mixed. Moreover, the proposed framework randomly selects operators to guide the search rather than adaptively selection based on the performance of each operator. Hence, the future works also may develop a strategy that can monitor the efficiency of operators and adaptively select the promising operators to optimize the problems.