Optimization of a Low-Carbon Two-Echelon Heterogeneous-Fleet Vehicle Routing for Cold Chain Logistics under Mixed Time Window

: Due to the rise of social and environmental concerns on global climate change, developing the low-carbon economy is a necessary strategic step to respond to greenhouse e ﬀ ect and incorporate sustainability. As such, there is a new trend for the cold chain industry to establish the low-carbon vehicle routing optimization model which takes costs and carbon emissions as the measurements of performance. This paper studies a low-carbon vehicle routing problem (LC-VRP) derived from a real cold chain logistics network with several practical constraints, which also takes customer satisfaction into account. A low-carbon two-echelon heterogeneous-ﬂeet vehicle routing problem (LC-2EHVRP) model for cold chain third-party logistics servers (3PL) with mixed time window under a carbon trading policy is constructed in this paper and aims at minimizing costs, carbon emissions and maximizing total customer satisfaction simultaneously. To ﬁnd the optimal solution of such a nondeterministic polynomial (NP) hard problem, we proposed an adaptive genetic algorithm (AGA) approach validated by a numerical benchmark test. Furthermore, a real cold chain case study is presented to demonstrate the inﬂuence of the mixed time window’s changing which a ﬀ ect customers’ ﬁnal satisfaction and the carbon trading settings on LC-2EHVRP model. Experiment of LC-2EHVRP model without customer satisfaction consideration is also designed as a control group. Results show that customer satisfaction is a critical inﬂuencer for companies to plan multi-echelon vehicle routing strategy, and current modest carbon price and trading quota settings in China have only a minimal e ﬀ ect on emissions’ control. Several managerial suggestions are given to cold chain logistics enterprises, governments, and even consumers to help improve the development of cold chain logistics.


Introduction
Along with the increasing development of people's living standards and the adjustment of health concept, the pursuit of fresh and perishable food has become a new normal and created a huge market for the cold chain industry. Cold chains are primarily employed for the preservation and transportation of temperature sensitive products such as foodstuffs and medical goods in the proper temperature range to slow biological decay processes and deliver safe and high-quality foods from the supplier to final consumption stage [1]. Due to the perishability of shipment, it is necessary to improve the efficiency of cold chain logistics to keep supplies fresh and cut waste. Meanwhile, the negative external impacts of freight transportation such as carbon emissions are equally deserving of attention [2]. Since both transporting and keeping the products at given temperatures are energy-intensive processes and generate lots of greenhouse gases (GHG) which are responsible for global climate change, there is a new trend and focus to establish a LC-VRP model of cold chain logistics to reduce CO2 and other GHG emissions through optimizing vehicle routes [3].
So far, although many scholars have spent a great deal of time on studying the vehicle routing problem (VRP) to optimize vehicle driving routes in cold chain logistics, the majority of them still took one echelon VRP as a primary concern, and complex multi-distribution layers in a real-life scenario did not get much attention. For instance, a flow of milk supply chain could even contain nine stages from cows to consumers associated with a single milk-processing facility [4]. Under such circumstances, investigating multi-echelon VRP for cold chain is necessary.
In addition, to help meet the goals of low-carbon economy, government regulators have decided to implement carbon trading policy which is one of the most effective carbon emissions' controlling mechanisms that creates both pressures and incentives to encourage companies to take action to reduce carbon emissions through setting a limited number of tradable emission allowances [5]. Companies also have implemented various internal and external carbon management practices in response to growing pressure from institutions and stakeholders [6]. To develop sustainable cold chain logistics, we shall take economic and carbon emission considerations simultaneously into an optimization model based on the appropriate carbon trading policy.
Confronting the flourishing, prosperous and immerse market demands, many 3PL are rushing to launch and conduct cooperation with cold chain suppliers. In China, ST Express established Shen-Xue Cold Chain Management Co., Ltd., last year, to provide same-day delivery and next-day delivery services for refrigerated and frozen items; SF Express recently announced the establishment of a joint venture with the American cold chain logistics company, Xia Hui, which is the main logistics service provider of McDonald's [7]. Furthermore, in today's fierce market, how to sustain customer loyalty is crucial for an enterprise's development. Surveys show that the customer loyalty is closely bound up with satisfaction, that is, the higher customer satisfaction of service or product, the higher customer loyalty a company will achieve [8,9]. For this point, corporate leaders set sights on satisfying the consumer need and improve the service quantity to raise customer satisfaction [10]. It is obvious that delivery time of goods is a non-negligible factor affected by customer satisfaction in cold chain logistics and the vehicle routing planning takes great effect on efficiency and promptness of delivery service. As such, it is wise for 3PL to pay attention to customer satisfaction when planning cold chain distribution routes.
In summary, an optimal low-carbon multi-echelon VRP model in cold chain logistics would help 3PL, not only in saving economic costs but also to reduce carbon emissions and raise customer satisfaction. In this work, a two-layer logistic network is taken into account. As the optimization process is NP-hard, AGA is employed to solve the proposed LC-2EHVRP model. Detailed research contributions are given as follows: • This paper proposes an adaptive genetic algorithm for optimizing LC-2EHVRP with mixed time window while considering economic costs, carbon emissions and customer satisfaction simultaneously.

•
Compared with previous VRP studies in the cold chain field, the factors that could influence carbon emissions such as multistage distribution, alternative location and path selection, different vehicle type, and time window are considered, in order to conform to real-world scenes. • A mixed time window which adds a penalty mechanism of soft time window to the strict rejection of hard time window is adopted to accurately evaluate customer satisfaction.
The rest of the paper is arranged as follows. The literature review associated with the problem addressed is presented in Section 2. Then, the construction of LC-2EHVRP model is established in Section 3. The fourth section proposed an AGA and its solution procedure. Computational results of Sustainability 2020, 12,1967 3 of 22 benchmark test and a real cold chain case study are shown in Section 5 followed by discussion and managerial implication. Lastly, conclusion and future suggestions are drawn in Section 6.

Literature Reviews
The main goal of this paper is to develop a LC-2EHVRP model with mixed time window, simultaneously considering economic cost, environmental issue (carbon emissions) and customer satisfaction for 3PL in cold chain logistics and obtain an optimal solution to deal with it. Consequently, the focus of the literature survey in this study is on the research about VRP for cold chain logistics, the research about VRP with environmental concern, and the research about VRP with customer satisfaction concern.

Vehicle Routing Problem for Cold Chain Logistics
Timely delivery of freight and reduction of transportation costs are the pressing tasks which the cold chain industry are facing, thus this has been attracting a great number of scholars devote to VRP in the cold chain field. VRP originated from a salesman's problem research in 1954 [11], and was first introduced by Golden et al. as the phrase "vehicle routing" [12]. In the past years, scholars used a variety of methods such as genetic algorithm (GA), tabu algorithm, particle swarm optimization (PSO) and other heuristics or hybrid methods to find optimal solutions of classic VRP, open VRP, capacitated VRP and heterogeneous VRP applied in city logistics, waste collection and other related fields [13][14][15][16]. Later, some intermediate distribution centers were embedded into the delivery path as an extension of the classical VRP, named multi-echelon VRP, which is quite popular among supply chain and logistics researches. Wang et al. established a two-echelon logistics distribution network model adopting two types of vehicles and paid a foundation for two-or multi-echelon logistics distribution network optimization [17]. Dondo et al. introduced a mixed-integer linear programming formulation to address the N-echelon vehicle routing and scheduling problem with cross-docking considering time constraints and multiple items [18]. Wang et al. combined a k-means clustering algorithm and an improved non-dominated sorting genetic algorithm-II to study the two-echelon collaborative multiple centers VRP [19]. Huai et al. solved a multi-type VRP in the cold chain logistics system with genetic algorithms based on two decoding methods [20].
Related to cold chain researches, Zhang et al. employed a seeker genetic algorithm to find the path with lower cost in a one echelon cold chain logistics model [21]. Chen et al. formulated the production scheduling and vehicle routing problem under stochastic demand for perishable food products as a mixed-integer nonlinear programming model to maximize the expected total profit of the supplier [22]. Lu et al. studied a location inventory routing problem for cold chain logistics network problem under an uncertain demand environment with a new discrete particle swarm optimization [23]. Govindan et al. proposed a novel hybrid approach called MHPV to solve the two-echelon location routing problem considering heterogeneous fleet of perishable food supply chain [24].
From the above studies, we can find that most scholars tend to employ various methods to solve one echelon VRP model in cold chain logistics, while cold chain consisting of multiple layers did not receive enough concerns.

Vehicle Routing Problem with Environmental Concerns
Of recent years, the issue of global warming has drawn widespread attention and companies have adopted the sustainability logic in response to reflect the environmental concern [25]. Several researchers have set their sights on reducing GHG generated from vehicle transporting and distribution process. Li et al. [26] took the full set of GHG emissions composed of CO2, CH4 and N2O into account when modeling the green cold chain VRP. Kazemian et al. developed a novel model for capacitated VRP model with time-dependent speed limits aiming to minimize two green objectives, total fuel consumption and GHG emission [27].
Since carbon emissions are mainly blamed for global warming, governments all over the world have made regulation policies such as carbon tax, and the carbon trading (cap-and-trade) scheme were adopted in the researches [28]. Liu et al. investigated a minimal-carbon-footprint time-dependent heterogeneous-fleet VRP and solved by GA [29]. Hariga et al. considered the impact of accounting for carbon emissions based on carbon tax scheme when analyzing a multistage supply chain with a single facility, a single distribution center and a retailer [3]. Lin et al. built a mode of low-carbon urban distribution path optimization and analyzed various cost components of urban distribution for perishable food and corresponding constraints [30]. Shen el al. introduced a carbon trading mechanism in multiple depots of VRP with time window to achieve carbon reduction goals [31]. Qin et al. also introduced a carbon trading mechanism to calculate carbon emissions' costs in cold chain logistics [32].
To sum up, carbon trading and carbon tax policy are two main mechanisms to evaluate carbon emissions in the cold chain.

Vehicle Routing Problem with Time Window Constraints
There is no specific requirement of delivery time for VRP without time windows; however, in the case of time window constraints, the model must not only consider optimization of vehicle routing, but consider early or delay delivery and its corresponding loss. The time window constraints can be divided into three types: hard time window, soft time window and mixed time window. Hard time window characterizes an allowable time period with the earliest starting and the latest ending point, and services that cannot be fulfilled within this time window are rejected [33]. Soft time window released hard constraints through accepting but applying penalty on early or late arrivals [34,35]. Mixed time window added a penalty mechanism of soft time window to the strict rejection of hard time window by setting both inner time window and outer time window. A service is allowed to conduct within the outer time window with a certain penalty when it cannot be fulfilled within the inner one. However, services also should be rejected when it arrived out of the outer time window [35]. Leng et al. applied a multi-objective hyper-heuristic approach (MOHH) to obtain Pareto-optimal solutions of a multi-objective regional low-carbon location routing problem (MORLCLRP) with three practical constraints: simultaneous pickup and delivery, heterogeneous fleet, and hard time windows [36]. Du et al. investigated a real-time vehicle-dispatching system with mixed time window for consolidating milk runs [35].
To the best of our knowledge, only few researchers have taken consideration of mixed time window constraints in the VRP study.

Vehicle Routing Problem Considering Customer Satisfaction
Besides environmental concerns in VRP, researchers also have seriously taken customer satisfaction on board. Since customer satisfaction represents a measure of the fulfillment of customer expectations and generally depends on the delivery time of goods compared to the required time window [37,38], several customer satisfaction measurements with time constraint settings have been developed. Researches reveal that customer satisfaction is inversely proportional to the vehicle waiting time, and the shorter the waiting time is, the higher satisfaction level is achieved [33,39]. Zhang et al. used an improved fuzzy due-time window to represent customer satisfaction degree [40]. Sivaramkumar et al. introduced a gap time which is difference between ready time and issuing time to improve customer satisfaction which helps managers to retain customers [41]. Guerriero et al. evaluated customer satisfaction according to the instance of a drone's arriving time at event's location under a soft time condition [38]. Qin et al. assumed that customers will be completely satisfied if the vehicle arrives within the time window requested by them, otherwise, they will be dissatisfied [32].
The above survey gave a holistic view of the research trends in VRP of cold chain logistics, and it is evident that modeling and solving the multi-echelon VRP for cold chains while considering economic, environmental and customer satisfaction issues could be an interesting and relevant research topic. In this paper, from an energy-saving and cost reduction perspective, we proposed a LC-2EHVRP model Sustainability 2020, 12, 1967 5 of 22 and followed the scholarly article by [32] to considering the concerns of environmental and customer satisfaction of the cold chain distribution network for 3PL. A mixed time window is introduced to measure customer satisfaction more precisely. At last, the extended model is solved using an AGA which is tested and verified by numerical experiments.

Problem Description
This study analyzed a LC-2EHVRP for 3PL in cold chain logistics as shown in Figure 1 that can be described as follows: There are several depots and multiple desired customers in the cold chain networks, between depots and customers there are some smaller intermediate depots, that is, satellites. Each depot has certain number of refrigerated vehicles which are identical, meanwhile, each satellite also has a certain number of another type of refrigerated vehicles which are identical, too. Freights are delivered from depots to satellites (1st-echelon) and then, transferred from satellites to final customers (2nd-echelon) by vehicles based on customers' demands along with time requirements. At satellites, when the 1st-echelon vehicles arrive at and unload their cargoes, products must be immediately loaded into the 2nd-echelon vehicles.
During the transportation process, the value of products was being decayed with time. In addition, the evaluation of customer satisfaction level is according to delivery time of the 2nd-echelon vehicle. This two-echelon low-carbon cold chain model is constructed to optimize economic and environmental objectives of cold chain logistics network for the third-party logistic service provider (3PL), moreover, maximize total customer satisfaction. Other assumptions related to this model are given as follows:

Problem Description
This study analyzed a LC-2EHVRP for 3PL in cold chain logistics as shown in Figure 1 that can be described as follows: There are several depots and multiple desired customers in the cold chain networks, between depots and customers there are some smaller intermediate depots, that is, satellites. Each depot has certain number of refrigerated vehicles which are identical, meanwhile, each satellite also has a certain number of another type of refrigerated vehicles which are identical, too. Freights are delivered from depots to satellites (1st-echelon) and then, transferred from satellites to final customers (2nd-echelon) by vehicles based on customers' demands along with time requirements. At satellites, when the 1st-echelon vehicles arrive at and unload their cargoes, products must be immediately loaded into the 2nd-echelon vehicles.
During the transportation process, the value of products was being decayed with time. In addition, the evaluation of customer satisfaction level is according to delivery time of the 2nd-echelon vehicle. This two-echelon low-carbon cold chain model is constructed to optimize economic and environmental objectives of cold chain logistics network for the third-party logistic service provider (3PL), moreover, maximize total customer satisfaction. Other assumptions related to this model are given as follows: Any vehicle that arrives early or late will incur a loss of customer satisfaction; • The freight cannot be managed by direct shipment from depots to final customers.

Notations
For easy reference, this section gives the description of notations that will be used in this paper which are listed in Table 1.

Sets
Description K 1 The set of vehicles in the 1st-echelon, The set of vehicles in the 2nd-echelon, k 2 ∈K 2 .

Notations
For easy reference, this section gives the description of notations that will be used in this paper which are listed in Table 1. Table 1. Description of notations.

Sets Description
The set of vehicles in the 1st-echelon, The set of vehicles in the 2nd-echelon, k 2 ∈ K 2 . D The set of candidate depots, d ∈ D. S The set of candidate satellites, s ∈ S. C The set of candidate customers, c ∈ C. N 1 The set of locations consist of {D ∪ S}. N 2 The set of locations consist of {S ∪ C}. Carbon emissions generated by per unit electricity consumption.

Decision variables Description
x i jk 1 0-1 variable, 1 if vehicle k 1 traverses arc (i, j), i, j ∈ N 1 and 0, otherwise. y i jk 2 0-1 variable, 1 if vehicle k 2 traverses arc (i, j), i, j ∈ N 2 and 0, otherwise. Q i jk 1 Quantity of goods delivered by vehicle k 1 traversing arc (i, j), i, j ∈ N 1 Q i jk 2 Quantity of goods delivered by vehicle k 2 traversing arc (i, j), i, j ∈ N 2 T ik 1 Arrival time of vehicle k 1 on the 1st-echelon, i, j ∈ N 1 . t ik 2 Arrival time of vehicle k 2 on the 2nd-echelon, i, j ∈ N 2 .

Fixed Costs
The fixed cost of operational vehicles traversed in two-echelon arcs should be considered, such as trucks' daily maintenance and repairing costs, driver wages and other expenses that have no concern with traveling distances. The fixed cost can be expressed as:

Deterioration Costs
Deterioration costs are losses caused by products' spoilage, decay and deterioration during the transportation and service (doors are being opened) processes in cold chain logistics. Based on Wang et al.'s model [42], percentage of freight deterioration in cold chains can be measured as 1 − e −θt × 100%. Thus, the deterioration cost in the cold chain network can be expressed as:

Energy Costs
The energy cost of entire cold chain logistics is composed of two types of resources' costs, fuel cost and electricity cost associated with transmitting, servicing and refrigerating process.
Vehicles burn a lot of fuel resources in road transportation. According to Soysal et al.'s model [43], the total amount of fuel consumption for a truck carrying load, load, with an average speed speed for traversing a distance, speed, is evaluated as follows: where k denotes vehicle engine friction index (kJ/rev/liter), N e is the speed of engine (rev/s), V is engine's displacement (L), µ is the vehicle curb weight (kg), λ = ξ κψ , ξ is fuel-to-air mass ratio, κ is the heating value of a typical diesel fuel (kJ/g), ψ is a conversion factor from grams (g/s) to liters (L/s), γ = 1 1000ε , ε is vehicle drive train efficiency and is an efficiency parameter for diesel engines, β = 0.5 C d A ρ, C d is the coefficient of aerodynamic drag, A is the frontal surface area (m 2 ), ρ is the air density (kg/m 3 ), α = gsinφ + g C r cos φ, ϕ is the road angle, Cr is the coefficient of rolling resistance.
We introduced ρ 1 , ρ 2 , ρ 3 to simplify the formulation, where ρ 1 = λ k N e V, ρ 2 = λ γ β, ρ 3 = λ γ α. Then, the fuel consumption of refrigeration in cold chain logistics is presented by: The fuel cost in cold chain transportation process is presented by: In order to prevent perishability of cool goods from transportation and unloading processes, vehicles need to install refrigeration equipment to keep products always in an optimal low-temperature environment. The electricity cost of cold chain distribution for 3PL are composed of refrigeration consumption costs related to the heat load generated during transportation and uploading periods. Here we have learned from the calculation of carriage's heat load in the model developed by Xie et al. [44].
where G(t) represents the heat load of carriage (kcal/h); R denotes heat-transfer coefficient of carriage (kcal/ (hm 2 • C), S e , S i are average external and internal surfaces of carriages (m 2 ), respectively; ∆T is the difference between the outside temperature and inside temperature of the refrigerated vehicle. Hence, the electricity consumption of refrigeration in cold chain logistics is presented by: The electricity cost generated during transportation and unloading processes is presented by: Thus, the total energy consumption costs in the two-echelon cold chain distribution logistics can be expressed as:

Carbon Emissions
In order to calculate carbon emissions' cost, we applied a carbon trading mechanism which is created by government to promote carbon reduction. Under carbon trading policy, a company is distributed tradable carbon emissions' allowances by government and could freely trade those allowances according to their total emissions [5]. Total carbon emissions are caused by fuel consumption and electricity consumption, and the carbon emissions' cost in the cold chain distribution logistics can be expressed as:

Customer Satisfaction
In cold chain logistics, delivering frozen or perishable goods on time prevents the quality decay and enables corporations to achieve a high level of customers' satisfaction and trustworthiness [45]. In this paper, we introduced a mixed time window [21] which is fuzzed to evaluate a consumer's satisfaction with a membership function, as expressed in Formula (11).
where [e i , l i ] represents the inner time window and the expected service time of customer i and e i , l i denotes the maximum tolerated hard time range. Since the closer the vehicle arrival's time to the customer required time window, the higher level of customer satisfaction will be achieved. ð is the sensitivity index. In this paper, we set it as "3", see Figure 2.

Model Establishment
In summary, our model aims to minimize the economic cost, the carbon emissions, and maximize the total customer satisfaction simultaneously. Through the comprehensive analysis above, carbon emissions can be transferred to emission cost according to carbon trading mechanism and customer satisfaction also can be calculated. Therefore, we integrate these three objectives into a

Model Establishment
In summary, our model aims to minimize the economic cost, the carbon emissions, and maximize the total customer satisfaction simultaneously. Through the comprehensive analysis above, carbon emissions can be transferred to emission cost according to carbon trading mechanism and customer satisfaction also can be calculated. Therefore, we integrate these three objectives into a mathematical model which can be computed as: Subject to: i∈s j∈s (29) Constraint (13) ensures that each customer in the 2nd-echelon layer is only visited once. Constraints (14) and (15) ensure that each vehicle, at most, departs from one starting point. Constraints (16) and (17) ensure that both types of vehicles' capacity should be respected. Constraints (18) and (19) ensure that the capacity of each depot or satellite should be respected. Constraints (20) and (21) eliminate sub-tours and ensure that each vehicle must return to the original starting location after assignment completed. Constraints (22) and (23) indicate the traveling time of each vehicle. Constraint (24) ensures that total quantity of delivered goods that can satisfy all customers' requirements in cold chain network. Constraint (25) ensures each vehicle's uploading quantity meets its visited customers' demands. Constraint (26) guarantees that traveling among depots is not allowed, respectively. Constraints (27) and (28) indicate products' balancing flow between two echelon layers. Constraint (29) shows before-mentioned decision variables and parameters are non-negative.

Algorithm
Since the LC-2EHVRP we proposed is NP-hard and it is too large and complicated to be solved by CPLEX, GA, a search heuristic, adopts the evolutionary behavior of biological systems, and is typically employed as a versatile optimization method [46]. However, the difficulty in GA's application is its premature convergence problem, in many cases, the algorithm generates a local optimum rather than global [47]. In order to combat premature convergence and solve our model efficiently, an AGA is adopted in this work.

Initial Population Generation and Chromosome Representation
Generating initial population is a critical part of GA which directly conditions the speed and convergence of the algorithm. In this paper, population was generated randomly according to the population size. In GA, a chromosome is a set of parameters which represents a proposed solution to the problem that needs to be solved [48]. In this work, each chromosome consists of three sub-strings as shown in Figure 3. The first sub-string has K 1 * (S + 1) genes which represents the vehicle routes in the 1st-echelon; the second sub-string lists total customers with random sorting and the last sub-string shows cutting points of vehicles driving in the 2nd-echelon. Each gene of vehicle starting in the first sub-string could be an assignment of depot number or '0', the same is true with the gene of distribution centers' vehicles passed and genes in last sub-string. The second and the third sub-string correspond to each other and constitute the driving routes in the 2nd-echelon. Generating initial population is a critical part of GA which directly conditions the speed and convergence of the algorithm. In this paper, population was generated randomly according to the population size. In GA, a chromosome is a set of parameters which represents a proposed solution to the problem that needs to be solved [48]. In this work, each chromosome consists of three sub-strings as shown in Figure 3. The first sub-string has K 1 * (S+1) genes which represents the vehicle routes in the 1st-echelon; the second sub-string lists total customers with random sorting and the last substring shows cutting points of vehicles driving in the 2nd-echelon. Each gene of vehicle starting in the first sub-string could be an assignment of depot number or '0', the same is true with the gene of distribution centers' vehicles passed and genes in last sub-string. The second and the third sub-string correspond to each other and constitute the driving routes in the 2nd-echelon. For instance, we consider the problem that has 2 candidate depots (numbers 1 to 2), 2 candidate satellites (numbers 3 to 4), 12 customers (numbers 5 to 16) and 3 available vehicles in the 1st-echelon, 4 available vehicles in the 2nd-echelon. The chromosome representation of the three parts is shown in Figure 4: The sub-string 1 indicates that there are two paths for two vehicles departing from depot 1 and depot 2, respectively, in the first layer: a truck starting from depot 1 directly to satellite 4 is coded as 1-0-4 and a truck starting from depot 2, then goes through satellite 3 and satellite 4 is coded as 2-3-4. Afterwards, in the second layer, the delivery routes of 12 customers are divided by cutting points in the sub-string 3. That is, 4 driving paths of 4 trucks which are 10-13, 11-6-15-7, 12-9-14, 5-16-8. For instance, we consider the problem that has 2 candidate depots (numbers 1 to 2), 2 candidate satellites (numbers 3 to 4), 12 customers (numbers 5 to 16) and 3 available vehicles in the 1st-echelon, 4 available vehicles in the 2nd-echelon. The chromosome representation of the three parts is shown in Figure 4: The sub-string 1 indicates that there are two paths for two vehicles departing from depot 1 and depot 2, respectively, in the first layer: a truck starting from depot 1 directly to satellite 4 is coded as 1-0-4 and a truck starting from depot 2, then goes through satellite 3 and satellite 4 is coded as 2-3-4. Afterwards, in the second layer, the delivery routes of 12 customers are divided by cutting points in the sub-string 3. That is, 4 driving paths of 4 trucks which are 10-13, 11-6-15-7, 12-9-14, 5-16-8.
For instance, we consider the problem that has 2 candidate depots (numbers 1 to 2), 2 candidate satellites (numbers 3 to 4), 12 customers (numbers 5 to 16) and 3 available vehicles in the 1st-echelon, 4 available vehicles in the 2nd-echelon. The chromosome representation of the three parts is shown in Figure 4: The sub-string 1 indicates that there are two paths for two vehicles departing from depot 1 and depot 2, respectively, in the first layer: a truck starting from depot 1 directly to satellite 4 is coded as 1-0-4 and a truck starting from depot 2, then goes through satellite 3 and satellite 4 is coded as 2-3-4. Afterwards, in the second layer, the delivery routes of 12 customers are divided by cutting points in the sub-string 3. That is, 4 driving paths of 4 trucks which are 10-13, 11-6-15-7, 12-9-14, 5-16-8.

Fitness Evaluation
The algorithm searched for optimal vehicle routes with the highest fitness value where the fitness function is established based on objective function: where f i represents the fitness value of chromosome i and represents the objective function value of chromosome i.

Crossover
Different crossover methods are applied to different sub-strings of the chromosome, in this paper, according to their structures. For the first sub-string, we use 2 points crossover operator which is to choose any two points parent chromosomes and exchange the nodes between them to generated children [49]. For the second sub-string, an inver-over operator is employed which is the number of inversions applied to a single chromosome and then, which segment to be inverted is determined by

Fitness Evaluation
The algorithm searched for optimal vehicle routes with the highest fitness value where the fitness function is established based on objective function: where f i represents the fitness value of chromosome i and z i represents the objective function value of chromosome i.

Crossover
Different crossover methods are applied to different sub-strings of the chromosome, in this paper, according to their structures. For the first sub-string, we use 2 points crossover operator which is to choose any two points parent chromosomes and exchange the nodes between them to generated children [49]. For the second sub-string, an inver-over operator is employed which is the number of inversions applied to a single chromosome and then, which segment to be inverted is determined by another individual [50]. For the third sub-string, a kind of partially mapped crossover (PMX) operator is introduced [51][52][53].
In addition, we also employed an adaptive crossover mechanism [51] which can be expressed as: where P t c represents crossover probability of generation t; P cmax and P cmin represent the maximum and minimum probability among all historic crossover; t max is the maximum iteration number during evolution.

Mutation
Similar to crossover process, the mutation operators applied on each sub-string are also not the same. First sub-string and last sub-string of chromosomes need to regenerate during mutation process, while second sub-string needs to mutate with an inver-over operator [50]. Then, an adaptive mutation mechanism [54] is introduced as follows: where P t m represents mutation probability of generation t; P cmax and P cmin represent the maximum and minimum probability among all historic crossover; t max is the maximum iteration number during evolution; A is a constant suggested to be valued as "10".

Selection
Selection in GA significantly affects its convergence; the basic idea is the higher fitness value, the larger probability of surviving [55]. In this work, we chose a roulette-wheel selection which assumes that the probability of selection is proportional to an individual's fitness as a selection strategy [51].

Termination Condition
The algorithm evolves the initial population iteratively until the termination condition is matched. The termination condition is set as to whether the number of evolution generation is greater than predetermined iterations. If it is matched, the operating loop will be stopped, otherwise, the loop will be continued.

Computational Experiment
In this section, the proposed AGA is tested via a classical two-echelon capacitated vehicle routing problem (2E-CVRP) dataset developed by Perboil et al. [56]. Then, a real cold chain logistics case study is introduced to investigate the influence of external mixed time window and carbon policy setting on LC-HVRP model. Finally, discussions and managerial suggestions are given according to experimental results.
The previously mentioned algorithm has been implemented in MATLAB 2012a. All the experiments in this paper have been performed on a PC with Intel ® Core TM (Santa Clara, CA, USA) i7-3610QM CPU@ 2.30 GHz 2.30 GHz and memory of 8 GB.

Algorithm Experiment
Computational experiments of 15 benchmark instances of a two-echelon capacitated vehicle routing problem (2E-VRP) database [56] are carried out to investigate the effectiveness of proposed algorithm. The test instances in database cover up to 51 nodes and are organized in 4 sets: the first three are generated from VRP instances by Christofides and Eilon, which are denoted as the form of E-n13-k4-2 or E-n22-s8-14 [57]; while the last one is built from the dataset built by Cranic et al., which comprises customers and varying number of satellite distributions [58]. Fifteen test instances are selected from set 2 and set 3 and performed 20 times to obtain the optimal solution as the final solution by proposed AGA and classic genetic algorithm (GA). Table 2 presents the test results including final solutions and execution time, and compares them with the optimal solutions that are reported in the literature. The experiment's results show that AGA can yield a better solution in short running time compared with GA and standard benchmark in most cases. Although GA also reaches an optimal or near-optimal solution, it takes more time. Thus, the algorithm proposed in this paper has an advantage for solving two-echelon vehicle routing problems.
Then, in order to verify the applicability of proposed AGA approach, average optimization rate of three groups (nodes = 24, nodes = 35 and nodes = 53) of 15 instances and their corresponding average execution time are also calculated. Figure 5 gives an insight of the performance of AGA: it shows with the number of size increases in test instances, that the AGA needs to take more and more time to reach the optimal result and the best average optimization rate of objective function of 13.99% is occurred when size of 2E-CVRP instances is 24 nodes. While for a large group, where the nodes equal to 53, it reduces to 4.96%. The figure confirms that proposed AGA performs better for instances with small sizes. shows with the number of size increases in test instances, that the AGA needs to take more and more time to reach the optimal result and the best average optimization rate of objective function of 13.99% is occurred when size of 2E-CVRP instances is 24 nodes. While for a large group, where the nodes equal to 53, it reduces to 4.96%. The figure confirms that proposed AGA performs better for instances with small sizes.

Problem Data and Experimental Setting
The experiment is designed based on a distribution process of the cold chain logistics company in China with two layers and different vehicle fleet configuration. A cold chain logistics company called KY Express is responsible for transferring dairy products from 2 depots (factories), denotes as D1, D2 to 15 customers from C1 to C15, passing through three satellites S1, S2 and S3. Table 3 shows that the geographic location, demand, required time window and service time of a total 20 location points. Tables 4 and 5 illustrate the parameter setting of each vehicle fleet in the first echelon and the second echelon of logistics network, respectively. Other parameters related to final objectives are listed in Table 6. All data and parameters mentioned above in this section are employed in the experiment to verify the LC-2EHVRP model.

Problem Data and Experimental Setting
The experiment is designed based on a distribution process of the cold chain logistics company in China with two layers and different vehicle fleet configuration. A cold chain logistics company called KY Express is responsible for transferring dairy products from 2 depots (factories), denotes as D1, D2 to 15 customers from C1 to C15, passing through three satellites S1, S2 and S3. Table 3 shows that the geographic location, demand, required time window and service time of a total 20 location points. Tables 4 and 5 illustrate the parameter setting of each vehicle fleet in the first echelon and the second echelon of logistics network, respectively. Other parameters related to final objectives are listed in Table 6. All data and parameters mentioned above in this section are employed in the experiment to verify the LC-2EHVRP model.

Result Analysis
In an established LC-2EHVRP model, objectives are three-fold. In other words, final optimal solutions are not only associated with economic costs but also with the carbon trading mechanism adopted and customer satisfaction consideration.

Results with mixed time window changing
The optimal solution of LC-2EHVRP is affected by the mixed time window changing. For in the case of Section 5.2, inner time window can be traced by Table 3; in this part, we aim to set the external time window rationally. The case study is tested with every 5-min external time window changes and before the test starting, we set the carbon price (C e ) as 30 CNY/t which is close to annual carbon trading price in China in 2018 according to China carbon trading official website [59]. Meanwhile, as carbon quota in carbon trading mechanism is determined based on historic data, we assign carbon quota (Q) as 0 initially. The results of above test are illustrated in Figure 6, meanwhile, the optimal routes without any time window (hard time window) and under 20 min external time window are shown in Figures 7 and 8 In an established LC-2EHVRP model, objectives are three-fold. In other words, final optimal solutions are not only associated with economic costs but also with the carbon trading mechanism adopted and customer satisfaction consideration.

Results with mixed time window changing
The optimal solution of LC-2EHVRP is affected by the mixed time window changing. For in the case of Section 5.2, inner time window can be traced by Table 3; in this part, we aim to set the external time window rationally. The case study is tested with every 5-min external time window changes and before the test starting, we set the carbon price (C e ) as 30 CNY/t which is close to annual carbon trading price in China in 2018 according to China carbon trading official website [59]. Meanwhile, as carbon quota in carbon trading mechanism is determined based on historic data, we assign carbon quota (Q) as 0 initially. The results of above test are illustrated in Figure 6, meanwhile, the optimal routes without any time window (hard time window) and under 20 min external time window are shown in Figures 7 and 8, respectively.   In an established LC-2EHVRP model, objectives are three-fold. In other words, final optimal solutions are not only associated with economic costs but also with the carbon trading mechanism adopted and customer satisfaction consideration.

Results with mixed time window changing
The optimal solution of LC-2EHVRP is affected by the mixed time window changing. For in the case of Section 5.2, inner time window can be traced by Table 3; in this part, we aim to set the external time window rationally. The case study is tested with every 5-min external time window changes and before the test starting, we set the carbon price (C e ) as 30 CNY/t which is close to annual carbon trading price in China in 2018 according to China carbon trading official website [59]. Meanwhile, as carbon quota in carbon trading mechanism is determined based on historic data, we assign carbon quota (Q) as 0 initially. The results of above test are illustrated in Figure 6, meanwhile, the optimal routes without any time window (hard time window) and under 20 min external time window are shown in Figures 7 and 8, respectively.  From the line chart, we can see that every five minutes extends to customers' acceptable time range (no more than 30 min); carbon emissions and customer total satisfaction obtained in the cold chain logistics will change accordingly. Among all external time window points, 20 min earlier or later Sustainability 2020, 12, 1967 16 of 22 time window extending seems to be a noticeable one as carbon emissions is the relatively lower, and total customer satisfaction is highest in the experimental model. From Figures 7 and 8, we can see vehicle routes are changed under different time window constraints. In order to obtain better results, follow-up experiments will be done under 20-min external time window extending. From the line chart, we can see that every five minutes extends to customers' acceptable time range (no more than 30 min); carbon emissions and customer total satisfaction obtained in the cold chain logistics will change accordingly. Among all external time window points, 20 min earlier or later time window extending seems to be a noticeable one as carbon emissions is the relatively lower, and total customer satisfaction is highest in the experimental model. From Figures 7 and 8, we can see vehicle routes are changed under different time window constraints. In order to obtain better results, follow-up experiments will be done under 20-min external time window extending.

Influences of carbon trading and customer satisfaction considerations
Both carbon trade price and carbon quota in carbon trading mechanism affect total costs and carbon emissions in VRP model [31]. However, in this work, additionally considering customer satisfaction is also possible to affect the final environmental influences. Two models are tested in this paper to investigate the influences of those factors, the first one is the LC-2EHVRP model which aims at minimizing the total cost including economic costs and carbon emissions' costs, and another one is the LC-2EHVR model, additionally with added customer satisfaction consideration into the final objective. Based on the carbon emissions (406 kg) from external time window test, we set four carbon quotas as none, small (100 kg), medium (300 kg) and large (500 kg) with three carbon prices, Table 7 lists final computational results of two objective models, Figures 9 and 10 represent changing trends of carbon emissions and total costs in two models under different carbon trading mechanisms.

Influences of carbon trading and customer satisfaction considerations
Both carbon trade price and carbon quota in carbon trading mechanism affect total costs and carbon emissions in VRP model [31]. However, in this work, additionally considering customer satisfaction is also possible to affect the final environmental influences. Two models are tested in this paper to investigate the influences of those factors, the first one is the LC-2EHVRP model which aims at minimizing the total cost including economic costs and carbon emissions' costs, and another one is the LC-2EHVR model, additionally with added customer satisfaction consideration into the final objective. Based on the carbon emissions (406 kg) from external time window test, we set four carbon quotas as none, small (100 kg), medium (300 kg) and large (500 kg) with three carbon prices, Table 7 lists final computational results of two objective models, Figures 9 and 10 represent changing trends of carbon emissions and total costs in two models under different carbon trading mechanisms.    Table 7, Figures 9 and 10, we can obtain findings as follows: (1) Customer satisfaction is an importation influencer of VRP model. Comparing with LC-2EHVRP model, tested LC-2EHVRP model with customer satisfaction consideration shows higher figures on both total cost and carbon emissions. It reveals that taking customer satisfaction into account will simultaneously increase the total cost and cause more carbon emissions in the cold chain logistics network. However, the optimal solution of LC-2EHVRP model without customer satisfaction consideration will lead to more customers' complaints. It is in line with the findings in cold chain VRP study of Qin et al [32]. Hence, developing a low-carbon VRP model with customer satisfaction consideration is necessary and effective.
(2) No matter how carbon price and quota change, carbon generations of LC-2EHVRP model without customer satisfaction consideration are the same all the time, which means its optimal    Table 7, Figures 9 and 10, we can obtain findings as follows: (1) Customer satisfaction is an importation influencer of VRP model. Comparing with LC-2EHVRP model, tested LC-2EHVRP model with customer satisfaction consideration shows higher figures on both total cost and carbon emissions. It reveals that taking customer satisfaction into account will simultaneously increase the total cost and cause more carbon emissions in the cold chain logistics network. However, the optimal solution of LC-2EHVRP model without customer satisfaction consideration will lead to more customers' complaints. It is in line with the findings in cold chain VRP study of Qin et al [32]. Hence, developing a low-carbon VRP model with customer satisfaction consideration is necessary and effective.
(2) No matter how carbon price and quota change, carbon generations of LC-2EHVRP model without customer satisfaction consideration are the same all the time, which means its optimal  Table 7, Figures 9 and 10, we can obtain findings as follows: (1) Customer satisfaction is an importation influencer of VRP model. Comparing with LC-2EHVRP model, tested LC-2EHVRP model with customer satisfaction consideration shows higher figures on both total cost and carbon emissions. It reveals that taking customer satisfaction into account will simultaneously increase the total cost and cause more carbon emissions in the cold chain logistics network. However, the optimal solution of LC-2EHVRP model without customer satisfaction consideration will lead to more customers' complaints. It is in line with the findings in cold chain VRP study of Qin et al [32]. Hence, developing a low-carbon VRP model with customer satisfaction consideration is necessary and effective.
(2) No matter how carbon price and quota change, carbon generations of LC-2EHVRP model without customer satisfaction consideration are the same all the time, which means its optimal routes under each condition are identical. This is owing to carbon emissions' costs based on carbon trading policy which makes up an extremely small share of final cost compared to other fixed costs and energy costs. Whereas when taking customer satisfaction into the final objectives of LC-2EHVRP, influenced by customer satisfaction consideration, carbon emissions decreased with the rising of carbon price under a large quota. When none or small carbon quota is given, tested LC-2EHVRP model always generates 371.86 kg carbon emissions and cannot be influenced by carbon price changing. When under a medium carbon quota condition, carbon emissions will reduce when carbon price is high enough. (3) Since surplus carbon quota can be traded in the market, total costs of the two models are declining in this case, with the escalation of carbon price under a large trading quota. Meanwhile, total cost of LC-2EHVRP model without customer satisfaction consideration proportionally decreases in this case, with the reducing of carbon price when carbon quota is not sufficient. However, for the model considering customer satisfaction, higher carbon price and larger trading quota will change costs nonlinearly as optimal solutions are different.

Discussion and Managerial Implication
In this paper, a proposed adaptive algorithm approach which is validated by dataset benchmarks is employed to the establishing of a LC-2EHVRP model based on a real case study of cold chain logistics. As compared with the traditional supply chain, cold chain logistics is more complicated and costly, not only at the economic level but also at the environment level. Additional refrigerating process of cold chain logistics with transportation will consume a lot of power and generate a considerable sum of carbon emissions. Thus, developing a low-carbon multi-echelon VRP cold chain model by searching optimal routes for vehicles in each echelon is the main goal of our study. In addition, the features of perishable or frozen stuff made customers more sensitive with delivery time. Therefore, customer satisfaction could be another consideration of the low-carbon VRP model. Several numerical tests are conducted on 5.1 and 5.2, and the main summaries are listed as follows: (1) Customer Satisfaction is a critical influencer of cold chain VRP model under mixed time window. Generally, the optimal routes for customers are not the best choice at economic or environment levels.
(2) Introducing mixed time window helps to measure customer satisfaction more accurately and an appropriate external setting time window contributes to reducing carbon emissions and improving satisfaction. (3) Modest carbon trading setting is difficult to prompt 3PL to control carbon emissions. Low-carbon model is only affected by carbon trading price and quota which reach certain higher figures.
Based on the above summaries, some suggestions are given to cold chain logistics enterprises, governments and even consumers and are proposed as follows: For cold chain logistics enterprises, trying to figure out the optimal routes and distribution process for cold chain logistics is an effective and cheaper way to reduce total costs, improve delivering efficiency and protect the environment. Besides, there is a trade-off between the environmental protection and improving of customer satisfaction. Customer requirements should not be neglected in the low-carbon cold chain. Furthermore, companies can launch a campaign to introduce their low-carbon logistics to get consumers' support of extending the time limit of delivery.
For governments, developing incentive or penalty mechanisms are crucial ways to control carbon emissions directly. As in our case study, according to current trading price from 30 CNY/t to 70 CNY/t in the carbon trading market, final carbon costs for cold chain logistics companies are too small to affect the distribution planning. Carbon price still has a huge expanding area in current status.
Lastly, protecting the environment is a common task requiring concerted efforts from all of mankind. In the cold chain logistics, excess dioxide is produced according to consumers' urgent requirements of delivery time. Extending limitations of delivery time properly can relieve stress of logistics enterprises and encourage them to be more environmentally friendly.

Conclusions
Facing the growing menace of global warming, adopting environmentally conscious practices which focus on carbon emissions' reductions has become a crucial social concern. This paper tackled the low-carbon VRP for a two-echelon cold chain logistics network which comprised of depots, satellites and customers. A LC-2EHVRP model derived from 3PL of cold chain logistics in China was presented which seeks to minimize costs, carbon emissions and maximize total customer satisfaction of the cold chain logistics simultaneously under a mixed time window. Moreover, an AGA approach was employed in this model, which is firstly validated by a benchmark test, then investigation of the influence of mixed time window changing and carbon price and trading quota settings in carbon trading policy on the optimal solution of proposed model.
With the numerical experiment of real cases in this article, we find that there is a trade-off among customer satisfaction, economic costs and environmental protection, and setting an appropriate external time window helps to reduce carbon emissions and satisfy customers in a multi-echelon VRP model based on mixed time window. Among all external time window points, 20 min earlier or later time window extending seems to be a noticeable one as carbon emissions are relatively lower and total customer satisfaction is highest under this condition. In addition, the sensitive analysis also provides valuable information that optimal solutions in low-carbon VRP model, without customer satisfaction consideration, are not sensitive to current modest carbon trading setting. While when taking customer satisfaction into the final objectives of LC-2EHVRP, carbon emissions decreased with the rising of the carbon price under a large quota. When none or small carbon quota is given, tested LC-2EHVRP model always generates 371.86 kg carbon emissions. When under a medium carbon quota condition, carbon emission will reduce when carbon price is high enough.
The findings of this work are expected to provide practical managerial implications for cold chain logistics enterprises, governments and customers. Cold chain logistics enterprises should try to find a balance between the customer satisfaction and carbon emissions when they are planning vehicle routing strategy. It is necessary for governments to make a proper carbon trading mechanism, as in China, current carbon price and carbon quota allowance still have a huge expanding area. By releasing time requirements properly, customers can encourage cold chain logistics enterprises to make more effective low-carbon vehicle routes when considering customer satisfaction.
Compared to traditional (one echelon) VRP model in cold chain logistics, a ME-VRP model can schedule multiple types of vehicles in each echelon, therefore, an interesting line of future research is investigating the effect of joint distribution of different vehicles. Another interesting one is to further develop a rational carbon control mechanism to provide incentive to enterprises to lower carbon emissions.