Model and Algorithm of Two-Stage Distribution Location Routing with Hard Time Window for City Cold-Chain Logistics

: Taking cold-chain logistics as the research background and combining with the overall optimisation of logistics distribution networks, we develop two-stage distribution location-routing model with the minimum total cost as the objective function and varying vehicle capacity in di ﬀ erent delivery stages. A hybrid genetic algorithm is designed based on coupling and collaboration of the two-stage routing and transfer stations. The validity and feasibility of the model and algorithm are veriﬁed by conducting a randomly generated test. The optimal solutions for di ﬀ erent objective functions of two-stage distribution location-routing are compared and analysed. Results turn out that for di ﬀ erent distribution objectives, di ﬀ erent distribution schemes should be employed. Finally, we compare the two-stage distribution location-routing to single-stage vehicle routing problems. It is found that a two-stage distribution location-routing system is feasible and e ﬀ ective for the cold-chain logistics network, and can decrease distribution cost for cold-chain logistics enterprises. been


Introduction
Many cities find it challenging to set up a high-efficiency city logistics system to increase freight efficiency and decrease the impacts of city distribution on city living conditions [1]. Macharis et al. [2] pointed out that urban goods distribution (UGD) has an important impact on the sustainable development of cities. It is necessary to find new solutions for the management of freight distribution in order to reach a higher level of efficiency. So a city logistics model and some urban freight traffic policies have been studied by many researchers such as Taniguchi et al. [3], Allen et al. [4], Taniguchi et al. [5], Anand et al. [6], Danielis et al. [7]. Allen et al. [8] showed that the use of urban consolidation centers (UCC) is presumed to provide more efficient distribution in an urban area, and it can decrease energy use and environmental impact. A UCC can be described as a logistics facility located in relatively close proximity to the geographic area that it serves, and with the range of terms used to refer to the UCC concept main including a public distribution depot, urban trans-shipment center, freight platforms, cooperative delivery system, urban distribution center, consolidation center (sometimes specific, e.g., retail, construction), pick-up drop-off location, offsite logistics support concept and so on [9]. UCCs are one of the most frequently implemented and studied city logistics initiatives, and according to Lagorio

Literature Review
The concept of the location-routing problem (LRP) can be traced back to 1961. Von Boventer first discussed the relationship between location selection and transportation cost in transportation problems [15]. LRP and its variants have been studied extensively in the past. Perl and Daskin proposed a multi-vehicle and multi-facility warehouse location-routing problem model with vehicle capacity constraints [16]. Li et al. [17] studied a three-tier distribution system with suppliers, warehouses, and multiple geographically dispersed retailers, where retailers could replenish goods from warehouses and suppliers. A model was established to minimize the long-term average cost within the system while meeting the demand of each retailer. Prodhon [18] put forward a linear programming model for multi-plan periodic location paths by reasonably defining variables and designed a hybrid evolutionary algorithm to solve it. Xu [19] introduced the damage cost of goods incurred in transit and established a location-route optimisation model considering two-stage transportation cost, damage cost, and penalty cost, and service level. A Genetic Algorithm-Particle Swarm Optimization algorithm was designed to solve the problem. However, the energy consumption cost was not considered in the model. Song et al. [20] studied the LRP of a multi-journey vehicle path, i.e., considering that the vehicle runs multiple routes within the travel constraint time, and a three-stage heuristic algorithm was designed to solve the problem. Zhao et al. [21] build a heterogeneous fleet two-echelon capacitated location-routing model for joint delivery in city logistics. Wang et al. [22] proposed a bi-objective model, and designed an improved algorithm. The effectiveness of the improved algorithm was demonstrated through a comparison. Koç et al. [23] established a location-routing problem model considering a heterogeneous fleet and time windows. A hybrid evolutionary algorithm combining multiple heuristics was designed to solve the problem. Leng et al. [24] considered multiple conditions in the regional low-carbon location-routing problem, such as simultaneous pickup and delivery, time windows. Yu et al. [25] and Zhao et al. [26] studied location-routing problem with simultaneous pick-up and delivery. To solve the problem, a simulated annealing (SA) heuristic algorithm was designed and a hyper-heuristic approach based on iterated local search was proposed, respectively.
For the cold-chain logistics problems, Yang et al. [27] studied a location model for perishable products in two-tier distribution centers, introduced the cost of goods damage into the objective function, designed an improved genetic algorithm to solve the model, and demonstrated the effectiveness of the model and algorithm through an example. However, the distance studied was calculated using the radial shape of the location and demand points. Zhao et al. [28] designed a satisfaction degree function according to service time windows and introduced a minimum envelope clustering analysis method and tabu search algorithm to solve the problem. Zheng et al. [29] constructed location inventory routing problem under a demand environment in a cold-chain logistics network, and non-dominated sorting introduced a multi-objective genetic algorithm (GA). Wang et al. [30] studied a cold-chain logistics distribution network considering carbon footprint, the model of minimum total cost including carbon emission cost was constructed, and a hybrid algorithm was designed to solve the model.
Most of the existing LRP research focused on a general logistics distribution network and single-stage location-routing problems. When it comes to two-stage LRP, most of the models assume that the transportation path between the distribution center and the transfer stations is radial, i.e, the vehicle only provides services for one transfer station at a time and then returns to the distribution center, and most of the algorithms are involved in single-stage optimisation of transfer station selection and paths without considering the coupling and collaborative optimisation of the two stages. In view of this, we studied a two-stage cold-chain logistics distribution network, including two-stage optimisation of transfer station locations and routing. A two-stage location-routing model with the minimum total cost associated with the hard time window is constructed, and different types of vehicles are considered for distribution tasks. An integrated approach is employed to design the algorithm to ensure the quality of the solution. Finally, through a case study, we verify the feasibility and effectiveness of the model and algorithm.

Problem Description
The cold-chain logistics problem considered in this study comprises a cold-chain logistics distribution center, multiple potential cold-chain logistics transfer stations and customer points. The distribution center need to deliver goods to the customers through transfer stations within a specified time window. In the system, the distribution center, transfer stations, and routing between them constitute the first-stage city logistics network. The transfer stations, customers, and routing between them constitute the second-stage city logistics network. As shown in Figure 1, we optimise locations of transfer stations and the two-stage distribution network vehicle routing problem under the condition of minimum total cost. effectiveness of the model and algorithm through an example. However, the distance studied was 91 calculated using the radial shape of the location and demand points. Zhao et al. [28] designed a 92 satisfaction degree function according to service time windows and introduced a minimum 93 envelope clustering analysis method and tabu search algorithm to solve the problem. Zheng et al. [29] constructed location inventory routing problem under a demand environment in a cold-chain 95 logistics network, and non-dominated sorting introduced a multi-objective genetic algorithm (GA).

96
Wang et al. [30] studied a cold-chain logistics distribution network considering carbon footprint, 97 the model of minimum total cost including carbon emission cost was constructed, and a hybrid 98 algorithm was designed to solve the model.

99
Most of the existing LRP research focused on a general logistics distribution network and 100 single-stage location-routing problems. When it comes to two-stage LRP, most of the models 101 assume that the transportation path between the distribution center and the transfer stations is

Problem Assumptions
To facilitate the study, the following assumptions are made: (1) the geographical location of distribution center, potential transfer stations, time window, demand of customers are known; (2) each customer can only be served by one delivery vehicle; (3) the demand of a single customer is less than the vehicle capacity. (4) traffic congestion is not considered; (5) the cargo load of a vehicle should not exceed its rated load; (6) the unit transfer cost for each transfer station is known and is a constant; temperature changes and the first stage-cargo losses are not considered; (7) transfer stations compete with each other, and enterprises can choose different transfer stations to provide services according to their own cost minimization; (8) the capacity of refrigerated transport vehicles for distribution center (first-stage route) is known, and the demand of a transfer station can be greater than the capacity of a transport vehicle, i.e. the demand of the first-stage distribution path can be split; (9) the capacity of refrigerated transport vehicles (second-stage route) for transfer stations is known, and can not be exceeded. Each vehicle will return to the starting point after completing tasks; (10) types of vehicles are different for deliveries from distribution center and transfer stations, and capacities of same type of vehicles are the same, and there are enough transport vehicles for the system.

Parameter and Variables
To build the model, the following parameters and variables are defined: Q 2 : Vehicle capacity in the second-stage route; X k ij is a 0-1 variable: when X k ij = 1, the vehicle k passes the road between transfer station(or distribution center) i and transfer station(or distribution center) j; otherwise, X k ij = 0; x k j is a 0-1 variable: when x k j = 1, the vehicle k services for customer i; otherwise, x k j = 0; x l ijg is a 0-1 variable: when x l ijg = 1, the vehicle l of transfer station g passes the road between customer (or transfer station ) i and customer (or transfer station) j; otherwise, X l ijg = 0; Z g is a 0-1 variable: when Z g = 1, the transfer station g is used; otherwise, Z g = 0; y l ig is a 0-1 variable: when y l ig = 1, the vehicle l of transfer station g provides service for customer i; otherwise, y l ig = 0.

Model Development
The two-stage distribution location-routing with hard time window for city cold-chain logistics model constructed in this paper takes the minimum total cost as the objective function. In consequence, the sub-cost should be analyzed firstly. Then the total cost of the two-stage distribution location-routing is obtained by the various sub-cost.

Objective Function Analysis of Model (1) Transportation Cost
The transportation cost mainly includes vehicle use cost, fuel used in the process of transportation, driver's cost, and other factors. For the convenience of research, the transportation cost is considered in two parts, referring to the use and consumption cost of vehicles per kilometre and driver's cost. The transportation cost of a refrigerated vehicle in the first-stage and second-stage routes can be expressed as follows: (2) Damage Cost Commodities in the cold-chain distribution are easily spoiled and therefore need to be kept in an appropriate low-temperature environment. The quality of perishable goods gradually declines or can lose value with time. When the quality of the product declines to a certain extent, spoilage cost will be incurred. The cost of damage is divided into two parts: the damage cost of goods accumulated over time in the process of transportation and the damage cost incurred when opening the door in the process of unloading. Because the first-stage delivery is usually carried out in an enclosed environment, spoilage cost will be minimal. We only considered the damage cost in the second-stage distribution. Hence, the total damage cost can be expressed as: The first part represents the cost of cargo damage in the transportation process, and the second part represents the cost of cargo damage in the unloading process. (

3) Refrigeration Cost
The cost of storage associated with maintaining the temperature and humidity inside the carriage is called the refrigeration cost. The refrigeration methods employed in refrigerated vehicles on the market today mainly include liquid nitrogen refrigeration, mechanical refrigeration, dry ice refrigeration, and cold plate refrigeration. The refrigeration method considered in this study is mechanical refrigeration, and the cost generated in the refrigeration process is related to fuel consumption, time, and weight of goods. Literature [31,32] proposed a formula for calculating the fuel consumption: W = ω εP e ξ * 10 −3 . The refrigeration cost includes the cost associated with the energy consumption of the vehicle in maintaining a low-temperature environment during delivery and the cost of additional energy supplied to the refrigeration system during the unloading process.
The refrigeration cost in the first-stage and second-stage route can be expressed as follows: where, W is the fuel consumption (g/h); ω is the power utilization coefficient of the refrigerator; ε is the fuel consumption rate (g/kW·h); P e is the effective power of the refrigerator (kW); ξ is the specific gravity of the fuel; W 1 is the fuel oil consumption during transportation; and W 2 is the fuel consumption during unloading. The first part represents the refrigeration cost of refrigerated vehicles during transportation, and the second part represents the refrigeration cost of the refrigerated vehicles during unloading.

(4) Penalty Cost
In an actual distribution process, the distribution vehicles may not arrive on time for various reasons, this will incur a penalty cost. The concept of a time window is introduced. As the timing for first-stage distribution is flexible, we only consider the time window requirements of the customers in the second-stage distribution. The time window is divided into a hard time window and a soft time window. Considering the characteristics of cold-chain distribution, we calculated the penalty cost associated with the hard time window. Assuming that the earliest service time allowed by customer i is ET i , the latest service time allowed by customer i is LT i , the service time window required by the customer i is [ET i , LT i ]. The penalty function equation can be expressed as follows [33]: where M is infinite. t is the time when the vehicle arrives at the customer. [ET, LT] is the service time window required by the customer.

(5) Transfer Cost
The cost at the transfer station mainly consists of the time cost incurred in the transfer of goods from large refrigerated vehicles to small refrigerated vehicles. Therefore, the transfer cost can be expressed as follows:

Model Setting
Based on the analysis of sub-cost in Section 3.4.1, a two-stage distribution location-routing model with the hard time window constrains for city cold-chain logistics is established as follows: Subject to: i∈N g g∈M 1 x l ijg = y l jg j ∈ N g , l ∈ L g (12) j∈N g g∈M 1 x l ijg = y l ig i ∈ N g , l ∈ L g (13) j∈N g j∈M 1 i∈N g q i y ig ≤ Z g q g g ∈ M 1 (17) The objective function of the model is shown in (7). Constraint (8) shows that a vehicle cannot exceed its maximum load in the first level distribution. Constraint (9) shows that vehicle can not exceed its maximum load in the second-level distribution. The first stage of distribution flow balance is shown in (10) and (11). The second-stage of distribution flow balance is shown in (12) and (13). Constraint (14) shows that there is only one refrigerated vehicle providing a delivery service for one customer. As per constraints (15) and (16), the refrigerated vehicles starting from a transfer station must return to the same after serving the customer, and the refrigerated vehicles starting from the distribution center must return to the distribution center after serving transfer stations. The total customer requirements assigned to a transfer station must be lower than or equal to the storage capacity of the transfer station, which is imposed by (17). Constraint (18) represents the time window for customer i.

Algorithm Design
Two-stage distribution location-routing belong to the NP-Hard problem, so we use heuristic algorithms to solve this problem. A hybrid genetic algorithm is proposed in the paper combining heuristic rules and distance clustering. The flow chart of the algorithm is shown in Figure 2.

Chromosome Coding
Each chromosome consists of three sub-strings: Sub-string 1 involves encoding length J. The integer of gene 1 − J is a non-repetitive sequence, representing the priority selection order and delivery order of the transfer station. Sub-string 2 is encoded as an integer with length N and gene 1 − N, representing the distribution order of the demand points. Sub-string 3 is encoded as an integer code of length 1, indicating the number of transfer station selected. For example, when J = 5 and N = 6, namely, there are 5 transfer station and 6 customer point, for the chromosome shown in Figure 3, the priority of the transition to be selected is 5, 4, 2, 1, 3. The distribution path of the demand points is 2-3-4-1-5-6; 2 indicates that the first two bits of the first layer of coding are selected to set up a transfer station. Moreover, 5, 4 coded in the first layer is the order of distribution of the transfer station.

Generate the Initial Population at Random
Genetic algorithm starts from the population of the problem solutions, so it is necessary to generate an initial population as the starting point of evolution. According to the encoding method, N chromosomes are randomly generated.
(3) Crossover operation of sub-string 3 Two-point crossover: random selection of two chromosomes as parents and two offspring chromosomes are obtained by direct exchange of two parent chromosomes. For example, select two parent chromosomes 1 and 2, the chromosomes of the crossed offspring are 2 and 1.

Mutation Operation
To maintain the diversity of the population, we conducted mutation operation for each sub-string in the chromosome. Interchange mutation operation is carried out in sub-string 1. Inversion mutation operation is chosen to perform on sub-string 2. Single-point mutation operation is used in sub-string 3. The detailed description of mutation operation is as follows: (1) Mutation operation of sub-string 1 Interchange mutation: Two random points are selected in the encoding string, and their positions are swapped. For example, select exchange points 3 and 7 from chromosomes 1 4 3| 9 2 6 7| 5 8, the chromosome obtained after exchange mutation is 1 4 7| 9 2 6 3| 5 8.
(2) Mutation operation of sub-string 2 Insertion mutation: randomly select two position points and insert the number of the first position point of the code string into the front of the second position in the code string. For example, we randomly select gene 4 from parent chromosome 1 4 |3 5 7 6| 2 9 8, and randomly select an insertion point 6, then the chromosome obtained after inserting the mutation is 1 3 5 7 4 6 2 9 8.
(3) Mutation operation of sub-string 3 Single point mutation: the gene is mutated by random variation. For example, the parent chromosome is 3, the chromosome of offspring is 2 after single point mutation.

Decoding Chromosome
The chromosomes are decoded to determine the priority selection order and distribution order of the transit point, distribution order of each demand point, and number of transit points. According to the principle of distance clustering and capacity of the transit points, each demand point is allocated to the selected transition point successively. Moreover, the route is divided according to the limit of the vehicle capacity, and vehicles are assigned to each transit point to cover all the demand points allocated to the transit point.

Calculating Individual Fitness Value
The calculation formula for the individual fitness value is F i = 1/Z i , where F i represents the fitness value of the i-th individual, and Z i represents the target function value of the i-th individual.

Selection Strategy Operation
Selection is the process of selecting the best individuals from the current population to produce a new population. We use the roulette selection strategy to select in this paper. The roulette selection is also known as proportional selection operator, and is to calculate the probability of each individual based on the fitness value of the individual. According to this probability, individuals are randomly selected to form the offspring population. This operation is repeated until the size of new population is identical to the original population size.

Algorithm Termination Condition
If gen > maxgen, output the best solution; this marks the end of the algorithm. Otherwise, let gen = gen + 1 and proceed to Step 4.

Description of Experiment and Analysis of Experimental Result
In the section, we use a randomly generated example to verify the validity of the model and the algorithm. Section 5.1 describes the experimental data generation and the parameters selection of the model. The experimental results are analyzed in Section 5.2. Comparative analysis of the results using different objective function of two-stage location-routing is shown in Section 5.3. In Section 5.4, we analyze the difference of distribution cost between two-stage vehicle routing and single-stage vehicle routing. The proposed algorithm and computations are carried out using MATLABR2017a package; all numerical experiments referenced in this paper are evaluated by a computer with windows7-x32 operating system, intel core i7 CPU @ 3.4 GHZ and 4 GB memory.

Experimental Data and Parameter Values
The test data for the calculation example are obtained using random generation. The customer demand is generated randomly within the interval [0, 3], and the x and y coordinates of the customer's position are generated randomly in [0, 100]. The service time is set to 10 min. The data pertaining to the distribution center and transit point are set according to practical significance. Tables 1 and 2 list the detailed data. The distribution center is represented by DC, the transit point by Pi (i = 1,2,3,4), and the customer point by natural numbers 1, 2, 3, . . . , 30. The values of the parameters in the model are: Q 1 = 25 t; Q 2 = 8 t; c 1 = 2 yuan/km; c 2 = 1.5 yuan/km; c 1 = 100 Yuan/h; c 2 = 60 yuan/h; when Q 1 = 25 t, W 1 = 1.2 kg/h, W 2 = 2.2 kg/h, v 1 = 60 km/h; when Q 2 = 8 t, W 1 = 0.8 kg/h, W 2 = 1.8 kg/h; v 2 = 40 km/h; P 1 = 4000 yuan/t; θ 1 = 0.003; θ 2 = 0.002; θ 3 = 9.3 yuan/kg; λ 1 = 50 yuan; and v 3 = 40 t/h. The values of the parameters in the algorithm are: the initial population size N = 100, gen = 1, max iteration number maxgen = 500, crossover probability p c = 0.8, and mutation probability p m = 0.1.

Analysis of Experimental Result
MATLAB programming was used to implement the above algorithm, and the test example was solved for 30 iterations. The calculation time was within 30 s, and the calculation efficiency was high.

Comparative Analysis of the Results under Different Objective Function of Two-Stage Location Routing
The proposed algorithm is employed to optimise the objective functions of the minimum total cost and minimum travel distance. MATLAB software was run 30 times to obtain the optimal distribution path number and cost of all levels, as listed in Table 3   Different objective functions have no impact on the transfer cost, total number of paths, and cost of the first-level paths. The transit cost is 2212.5 yuan, the total number of routes is 9, and the cost value, delivery distance, and travel time in the first-stage routes the same. The second-stage cost and routes are different for different objective functions. The objective function with the shortest distance saves 17.95 km in distance and 26.92 min in time compared to that with the minimum total cost, incurring 131.25 yuan more in total cost. Therefore, cold-chain participants can choose different route schemes to maximise their interests. If the carrier is the main participant, the route scheme with the lowest total cost may be the best choice. If logistics outsourcing is involved, the shortest travel distance may be the best choice for logistics providers. Tables 3 and 4 shows that, for single-stage vehicle paths, when the distribution vehicle capacity is 25 t, the vehicle path cost is 8001.18, when the distribution vehicle capacity is 8 t, the vehicle path cost is 8685.68, and for two-stage distribution location-routing, vehicle routing cost is 5757.35, transfer cost is 2212.5, total cost is 7969.85. Thus, the single-stage vehicle path increases the vehicle path cost. Despite the cost savings representing a small percentage of the total cost, the large activity of two-stage distribution in a big city may represent a substantial benefit per year. Beside the significant economic advantages yielded by this kind of distribution approach, two-stage distribution systems may result in being more efficient with respect to standard approaches and also from an environmental point of view. In addition, with the gradual implementation of vehicle restriction policy in China, the two-stage distribution scheme established can avoid large vehicles in the city center in this paper, thus saving distribution cost and improving operating profits.

Conclusions
Based on the characteristics of city logistics distribution, we have established two-stage distribution location-routing model with the hard time window for city cold-chain logistics. Different vehicle types were considered for the different distribution tasks. In order to solve the model, we designed a hybrid genetic algorithm, and the model and algorithm were verified through a random generated example. The effects of different objective functions on the location routing and total distribution cost were analysed, and the objective functions had no effect on the transfer location selection, but did have on the total cost. Therefore, different distribution optimisation schemes were proposed for different distribution subjects. A two-stage distribution location-routing problem and a single-stage vehicle-routing problem were compared and analysed. The results showed two-stage distribution scheme can help reduce overall distribution cost, indicating that the two-stage location-routing system is more feasible and effective. In conclusion, the two-stage location-routing model can provide a theoretical reference for the planning of distribution networks in the field of cold-chain logistics.
The operating conditions of city logistics distribution networks are very complex. Therefore, in future research, real traffic conditions, including potential congestion, and feasible geographical locations for transfer stations should be considered in the distribution planning. For example, the introduction of a road congestion index into the model can help to better reflect the actual feasible solution. Also, additional problem formulations in two-stage distribution, such as bi-level optimization where lower and upper optimization tasks are considered, will be faced in the next steps of future works jointly with comparison exercises with exact solutions and solvers (e.g., CPLEX).

Data Availability:
The data used to support the finding of this study are available from the corresponding author upon request.