# Models and Methods for Two-Echelon Location Routing Problem with Time Constraints in City Logistics

Authors: Peng Yang; Lining Zeng
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This paper focuses on the two-echelon location routing problem with time constraints in city logistics system. The aim is to define the structure of a system which can optimize the location and the number of two different kinds of logistics facilities as well as the related routes on each echelon. A mathematic model considering the problem characteristics has been set up. Based on probability selection principle, this paper first puts forward a metaheuristic algorithm with comprehensive consideration of time and space accessibility to solve it. Then, random instances of different sizes are generated to verify the effectiveness of our method.


---

## Body

## 1. Introduction

Many cities find it challenging to set up a high-efficiency city logistics system to increase the freight efficiency of city freight vehicles and decrease the impact of city distribution on city living conditions [1]. In China, more and more cities have been making great effort to reduce traffic in the city center, by creating peripheral logistic platforms or transfer stations, from which only smaller vehicles are allowed to go downtown. Therefore, the bigger-capacity trucks are usually prohibited from entering the city center, and the smaller trucks provide most of the distribution activities within the city. Multiechelon, especially two-echelon logistics system, has been becoming the most common version in practice [2]. The service providers in this system usually need to deliver on or before the time requested by the customer. In consequence, how to select locations for the node facilities of the system, how to dispatch vehicles, and how to plan appropriate vehicle routes have become the key problems in the city logistics system optimization [3].Transportation route optimization is a long-term research hotspot in the logistics system optimization field. Early in 1959, Dantzig [4] put forward Vehicle Routing Problem (VRP). Afterwards, various problems extended from VRP attracted extensive attention, such as time constraints, abnormal shape, open type, and backhauls. Two-Echelon Vehicle Routing Problem (2E-VRP) is an extension, aiming at the transportation routing problem of two-echelon logistics system. It was raised by Jacobsen early in 1980 [5]. Since 2E-VRP is an NP problem, which is more complicated than classical VRP, its solution algorithm has got wide attention. It includes deterministic algorithm and heuristic algorithm. In 2008, Crainic et al. [6] and Taniguchi et al. [7] put forward a heuristic algorithm based on cluster, by which they applied two single echelon vehicle routing problems to replace 2E-VRP problem and then got iterative solution of the two subproblems. In 2009, Perboli et al. [8] put forward some effective inequality constraints and applied branch and bound methods to solve 2E-VRP problem; however, such deterministic algorithm is only applicable to small-scale problem. In 2011, Perboli et al. [9] have additionally built two heuristic algorithms to solve 2E-VRP problem. In brief, seeking high-quality and high-efficiency algorithm has become the striving direction in this field.The above problems only refer to tactical planning decision-making. When a problem extends to the decision-making on both strategic and tactical planning, 2E-LRP problem is put forward, but few scholars have carried out research on this issue. Alumur and Kara [10] once built a multiple-target MIP model for the two-echelon LRP problem of a hazardous substance and got a solution through CPLEX. Sterle [11] and Boccia [12] applied tabu search algorithm to solve 2E-LRP problem. Nguyen et al. [13, 14] applied four present constructive heuristics and a hybrid metaheuristic to solve 2E-LRP problem. Considering a 2E-LRP problem with capacity limitation, Contardo et al. [15] built two-subscript vehicle follow model, applied branch and cut algorithm for solution, and got lower bound of the model. They put forward self-adaption large-scale neighborhood search algorithm and got upper bound of the model. Zhao [16] built a heterogeneous fleet two-echelon capacitated location-routing model for joint delivery arising in city logistics.The existing algorithms aiming to solve the 2E-LRP problem mentioned above mainly focused on tabu search algorithm and large-scale neighborhood search algorithm. There are few heuristic algorithms taking the problem structure characteristics into account. Considering two-echelon location routing problem with time constraints in city logistics system and its characteristics, this paper has set up relevant data structure and mathematic model and put forward a metaheuristic algorithm with comprehensive consideration of time and space accessibility and application of probability selection principle. To verify effectiveness of the algorithm, this paper has built relevant calculation instances for solution. The remainder of this paper is organized as follows: the description of the problem and a mathematical model are presented in Section2; the proposed algorithm is introduced in Section 3. Experimentation and concluding remarks are discussed in Sections 4 and 5, respectively.


## 2. Problem Description

The problem can be described as follows. The two-echelon city logistics network normally consists of a series of known distribution centers, potential transfer stations, and known customers. Due to practical problems, such as city transportation environment and transportation efficiency, goods are distributed first from distribution centers, via transfer stations, and finally sent to customers within given time. In the system, distribution centers, transfer stations, and relevant routes constitute first-echelon city logistics network; transfer stations, customers, and relevant routes constitute second-echelon city logistics network. Under the conditions such as known customers’ geographical position and quantity demand, distribution centers and potential transfer stations’ geographical position, and capacity limitation and vehicles’ capacity limitation in the logistics network at all echelons, we try to solve problems, including facilities locating, rationing, and vehicle routing in the two-echelon logistics network, to minimize the total cost.An illustrative example of the two-echelon logistics system is given in Figure1. The overall system is expressed as graph G=(N,E), in which node set N=P∪S∪C means set of distribution centers, transfer stations, and customers and E means set of edges. Let C be a set of customer demands. The demand of each customer i(1≤i≤k) is expressed as two-tuple combination ci=(wi,di), in which wi∈R+ means this customer’s demands quantity; di∈R+ means this customer’s required date for receiving the goods before di; that is, the time window of customer demand is [0,di].Figure 1
Example of the two-echelon logistics system.Distribution center setP consists of m distribution centers; any distribution center pi(1≤i≤m) is expressed as two-tuple combination pi=(pci,pli), in which pci∈R+ means fixed cost of the distribution center; once the distribution center is used, the fixed cost occurs; pli∈R+ means maximum load of the distribution center; the distribution center cannot allow overload; that is, the loading of distribution center cannot exceed pli.Transfer station setP consists of n transfer stations, any transfer station si(1≤i≤n) is expressed as two-tuple combination si=(sci,sli), in which sci∈R+ means fixed cost of transfer station; once the transfer station is used, the fixed cost occurs; sli∈R+ means maximum load of the transfer station; that is, the transferred goods quantity of transfer station cannot exceed sli.The transportation between distribution centers and transfer stations is the first-echelon transportation, and the limited capacity of each first-echelon vehicle is represented asL1. Correspondingly, the transportation between transfer stations and customers is the second-echelon transportation, and the limited capacity of each second-echelon vehicle is expressed as L2; obviously, considering the practical conditions, there is L1>L2.The edge between any two nodesN1,N2∈P,∪S∪C is expressed as EN1N2=(cN1N2,tN1N2), in which cN1N2∈R+ means the transportation cost of single vehicle at the transportation echelon; tN1N2∈R+ means the time for finishing the transportation between the two nodes. When N1,N2∈P∪S, the first-echelon transportation occurs; when N2∈C∧N1∈S∪C, the second-echelon transportation occurs. If there is no direct transportation route between two nodes, it is expressed as EN1N2=0. Since direct transportation cannot be realized between distribution centers and customers, it is obvious that when N1∈P∧N2∈C,  EN1N2=EN2N1=0 occurs.VehicleTi can be expressed as Ti={pathi,{CTi}}, in which pathi={EN1N2,EN2N3,...,ENkNk+1} means that vehicle Ti must pass k edges of EN1N2,EN2N3,...,ENkNk+1 and nodes including N1,N2,...,Nk+1 to finish the transportation task; {CTi}={CTi1,CTi2,...,CTik} means that the goods of vehicle Ti contain the transportation demand of k customers. The time of vehicle Ti passing node Nj is expressed as tTiNj, which represents the total time of vehicle Ti that starts from N1, passes N2,N3,...,Nj-1, and reaches Nj; the total time of vehicle Ti is expressed as tTi; obviously, tTi=tTiNk+1. The total goods weight of vehicle Ti is expressed as WTi=∑cj∈{CTi}wj.To meet any customer’s demandci, the goods must be transported by two different vehicles, which are responsible for the first-echelon and second-echelon transportation tasks, respectively. Therefore, the path of demand ci can be expressed as pathci=(pci,sci,Ti1,Ti2), which means that the goods path of ci starts from distribution center pci, transported by vehicle Ti1 to transfer station sci, after being transferred by vehicle Ti2 to the customer. In the process, ci must pass all edges of Ti1 before arriving at sci and must pass all edges of Ti2 before arriving at the customer. In this paper, we set that the customer demands cannot be split.Since customer requires finishing delivering task within a certain time constraint, it is necessary to consider the time consumed in the whole transportation process. The total consumed time includes three parts:(1) time for the first-echelon transportation, expressed as tTi1sci; (2) time for the second-echelon transportation, expressed as tTi2ci; and (3) waiting time at transfer station sci. The waiting time of vehicle means that, in the process of the first-echelon transportation, the goods to be transported by vehicle Ti2 may be from different first-echelon transporting vehicles, and the arrivals of these vehicles are not at the same time. Although the goods of ci arrive at transfer station sci, it is required to wait for arrival of other goods and then to realize the dispatch; this period constitutes the waiting time. The length of waiting time depends on the time taken for the arrival of the goods transported by the ci-used same vehicle to the transfer station during the second-echelon transportation; the waiting time is expressed as wtci=max⁡{tTj1sci∣Tj2=Ti2}-tTi1sci. It needs to be noticed that all goods demands are known, and the goods of distribution center can be prepared in advance. Therefore, no waiting time is consumed in the distribution center.The plan of the given system is a set of all customer demands andci-related vehicles and paths, which is expressed as follows.SetUP={up1,up2,...,upm} is expressed as the use condition of distribution center, in which upi(1≤i≤m) is a 0-1variable. upi=1 means that the distribution center is used; otherwise, the distribution center is not used.SetUS={us1,us2,...,usn} is expressed as the use condition of transfer station, in which usi(1≤i≤n) is a 0-1variable. usi=1 means that the transfer station is used; otherwise, the transfer station is not used.SetR={rN1N2} means the set of all paths in the result, in which rN1N2∈N+ means the number of vehicles from node N1 to node N2.Finally, the cost of the whole plan can be expressed as the sum of fixed cost of distribution centers, fixed cost of transfer centers, and path-to-path transportation cost; that is,(1)Z=∑pi∈Pupi∗pci+∑si∈Susi∗sci+∑N1,N2∈P,S,CrN1N2∗cN1N2A correct plan must satisfy the following constraints:(2)tTi1sci+tTi2ci+wtci≤di(3)WTi≤L(4)∑Ti⊂level1WTi=∑Ti⊂level2WTi(5)∑pj⊂pathTiWTi≤plj(6)∑sj⊂pathTiWTi≤sljConstraints (2) imply that all customer demands must be completed within the specified time constraints. Constraints (3) ensure that any vehicle cannot be overloaded, in which L depends on the echelon on which the vehicle lies. As for the first-echelon vehicle, L=L1; otherwise, L=L2. Constraints (4) indicate that the total amount brought to all transfer stations by the bigger-capacity trucks must be equal to the total demand of customers assigned to these transfer stations. Constraints (5) and Constraints (6) state that any goods through distribution centers and transfer stations cannot exceed the load limit.


## 3. Algorithm

### 3.1. Algorithm Framework

Since the above problems are NP-hard problems, this paper has designed a Heuristic Algorithm Based on Probability (HABP) for solution. The main ideas of the algorithm are the following:(a) In line 1, HABP seeks initial solution by function InitialSolution(G), which takes the Graph G as input. InitialSolution(G) computes a solution by FirstFit(FF) mechanism and return.(b) Weight number distribution: in line 2, according to the node-to-node comprehensive space-time distance consisting of single-vehicle transportation cost and transportation time, the function STDistance(node) makes initial processing of graphG in the system to distribute initial value for each node.(c) Iterative solution: in lines 3 to 5, based on the current solution and initial value, the HABP chooses customer nodes for optimization and replaces the current solution when finally better solution is found, where the function RandomRePath() is the core function of the algorithm, which seeks better solution in the iteration. RandomRePath() will be explained in detail in Section3.3.(d) Algorithm termination: when the number of iterations reaches the designated value, stop the search.The pseudocodes of HABP algorithm framework are shown in Algorithm1.Algorithm 1:HABP algorithm framework.
Input: two-echelon routing graphG Output: New solutionV 1.V = InitialSolution(G); 2. for every node inG: value(node) = STDistance(node); 3. for 1 tomaxIterationNum do 4.V=  RandomRePath(); 5. endfor; 6. returnV;


### 3.2. Composition of Initial Solution

A given two-echelon routing graphG with time constraints consists of m distribution centers, n transfer stations, and k customer demands. At the beginning of heuristic algorithm, we search for an initial solution one by one at first and then take it as the starting point for optimizing. For the algorithm in this paper, we adopt FirstFit (FF) algorithm to generate initial solution.The main ideas of FF algorithm are as follows: for the demand of each customer, we successively search for all potential routes. The first route satisfying the requirements is selected as the route for the demand. Here, we put aside the concentration of goods, assuming that all goods directly start off from the distribution center, arrive at the transfer station, and finally get to the customer. The optimizing is done in the later process of algorithm. In the process of searching for initial route, the two conditions should be met:(a) As for routepci→sci→ci satisfying any customer demand ci, this route is interconnected; that is, Epcisci≠0,Escici≠0.(b) This route transportation can meet the time requirement of customer demand; that is,tpcisci+tscici≤di.The algorithm of initial solution is shown in Algorithm2.Algorithm 2:FirstFit.
Input: two-echelon routing graphG Output: An initial solutionV 1. for every customer demand do 2.for i = 1 to m do 3.for j = 1 to n do 4.if(Epcisci≠0∧Escici≠0∧tpcisci+tscici≤di) 5.rpcisci+=1;rscici+=1; 6.break; 7. returnV;


### 3.3. Optimizing Strategy

For a given two-echelon routing graphG and initial solution, it is required to utilize optimizing strategy to seek better feasible solution and replace current solution. In the optimizing process, heuristic algorithm usually adopts random method to generate new solution and compare it with current solution. In this process, random strategy directly influences algorithm convergence rate. In HABP, problem solution means a vehicle operating plan satisfying all customer demands in G in transportation quantity and time. For two different solutions, on the condition that the time demand of all customers is satisfied, the solution with less total cost is the better one. In the process of seeking better solution, in order to improve convergence rate, it is required to choose a random value with bigger probability to obtain better solution when new solution is generated. Therefore, in the optimizing process, HABP puts forward the concept of selective probability and puts forward an optimizing mode based on probability. This mode contains the following steps:Step 1.
Allocate initial value for all nodes.Step 2.
Distribute probability for all demands and select optimizing node according to probability, based on the initial value and current solution of nodes.Step 3.
Replan for the selected node, generate new solution, and accept new solution with certain probability according to the comparison of new and old solutions.Step 4.
Return the optimal solution obtained in the process as a result.


#### 3.3.1. Initial Value of Node

In the two-echelon location routing environment, it is not completely equivalent between nodes. For example, for a customer demand, there are two transfer station nodess1 and s2 for choice, and the two routes have the same cost, but s1 has better connectivity than s2; that is, s1 is connected with more other nodes and has better cost; it is obvious that s1 is a better choice for this demand, for it has more possibility to share the transportation resource with other demands and bring optimization of total cost. Here, connectivity is measured from the two aspects of cost and time. About initial value, we consider that the node with less average cost and average time shall have higher probability for route accession.The cost and time initial values of various nodes are calculated in the formula below.About distribution center nodepi, its initial value is the mean value of the route of all transfer stations connected with it; that is, (7)vcpi=∑Episj≠0cpisjnumpi,vtpi=∑Episj≠0tpisjnumpiIn the formula,numpi means the number of transfer station nodes with route to connect node pi; vcpi means the initial cost of distribution center node pi, whose value is the mean cost value of the route connecting transfer stations; similarly vtpi means the initial time of distribution center node pi, whose value is the mean time value of the route connecting transfer stations.Correspondingly, for distribution center nodesi, its initial value is the mean value of the route of all distribution centers connected with it; that is,(8)vcsi=∑Esipj≠0csipjnumsi,vtsi=∑Esipj≠0tsipjnumsiFor distribution center nodeci, its initial values and mean values of all potential direct routes include the values from ci to transfer station node and from transfer station node to distribution center, since the values from transfer station node to distribution center can be expressed as vcsi and vtsi, which are formulated as follows: (9)vcci=∑Ecisj≠0ccisj+vcsjnumci,vtci=∑Ecisj≠0tcisj+vtsjnumciIt should be noted that, for transfer station node and customer node, the initial value does not take the path between the same nodes into account. This is because the path between the same nodes does not represent the actual capacity of the node but only reflects the convenience in combining goods.


#### 3.3.2. Selection of Optimizing Nodes

For a given graphG and a current solution V, V consists of distribution center utility state set UP, transfer center utility state set US, and route set R. The optimizing process means seeking a new solution with cost and time advantages by adjusting the content of current solution. For the algorithm in this paper, the optimizing process consists of two parts: determination of optimizing node and construction of new solution, in which selection of optimizing node is to determine the adjustment of current solution. In this algorithm, we assign probability to each customer demand, take random selection of customer node, and replan its route for optimization, based on node initial value and current solution.At first, we should define the cost of all demands in the current solution. For any demand, finishing the transportation may include four kinds of routes:(1)P→S; (2)S→S; (3)S→C; and (4)C→C. Among the four routes, route 1 and route 3 are the unique paths that the goods must pass, while route 2 and route 4 are not unique. For any of the routes, the vehicle may carry goods from multiple customers. Therefore, the transportation costs are shared by multiple customers. Those costs are accounted on one-route vehicle T; proportion of demand ci is formulated as λi:(10)λi=wiWTIn the above formula,WT means load of vehicle T on this route. For demand ci of pathci, its costi can be formulated as follows: (11)costi=∑EN1N2∈pathciλi∗cN1N2Considering the time cost ofci, in accordance with the definition, the transfer of all vehicles only occurs in the transfer from the first-echelon transportation to the second-echelon transportation; therefore, ci transportation shall and shall only be realized through the two vehicles Ti1 and Ti1, with the total time as follows: (12)ti=tTi1sci+tTi2ci+wtciAccording to the above parameters, we can assign the probabilities that are selected for optimization to different demands. We should make an analysis on the aspects of cost and time. As for the cost, a demand that requires high cost in the current solution, it must get distribution of higher probability, for such demand may have higher possibility for optimization. As a reference, the ratio of initial value to current costcosti/vcci can be considered as one of basic parameters for probability assignment. As for the time, any demand must meet the requirement of deadline; any delivery within the deadline shall be consistent in effect. However, generally, it is hoped to realize the delivery ahead of the deadline, under the condition that the cost is not raised. In fact, it can also improve the customer satisfaction. In this paper, we define the time impact factor of ci as χi, which is formulated as follows: (13)χi=1ti≤vtcitivtciti>vtciIn formula (13), when the time consumption of demand ci is higher than the mean value in the current solution, a time impact factor (>1) is assigned to it, so that it has higher probability to be selected for optimization.To sum up, in the HABP algorithm, the probabilitypi that a demand ci is selected for optimization is formulated as follows: (14)pi=costi∗χi∑cj∈Ccostj∗χj


#### 3.3.3. Construction and Acceptance Strategy of New Solution

After all demands are assigned with probabilities, we make a random selection of a demand for optimization on the probability base. The optimizing method is to cancel existing route and reselect a new route for this demand. In the selection process, the following factors should be considered:(a) Access point: for the route of a demandci, in the final section of the route, two possibilities may occur: starting from the transfer station node or starting from the customer node. It is clear that, in view of cost, “vehicle pooling” is useful for control of comprehensive cost. Therefore, during the planning of new route, “vehicle pooling” should be realized as much as possible under the premise that time and load constraints are not exceeded.(b) Treatment of original route: for the generation of new route, it is required to cancel the original route. Considering the following conditions, the vehicle transporting for a certain demandc1 on the second-echelon route has a route in the current solution, expressed as {s0,c1,c2}. When c1 is selected as the optimization node, route Es0c1 becomes unnecessary; at this time, if Es0c2≠0, we should make a comparison of the cost and time results between original route {s0,c1,c2} and new route {s0,c2}; the better one shall be selected.(c) Load: during the construction of new route, load is an important factor for consideration. At first, it should be ensured that no overload occurs. Here, “overload” includes overload of all nodes and overload of all vehicles in the whole transportation. For the overload route, probability should not be assigned. With no overload, we hope to select the route loaded with bigger loading capacity, because this means a reduction of unit transportation cost.(d) Time: the time consumption must be lower than the time constraint of the demand.(e) Selection of transfer station and distribution center nodes. During the construction of new route, we should select the transfer station and distribution center route which can bring better solution.Considering the above constraints, the construction of new route can be implemented in the following steps:(a) Selecting access point: at first, we should select access point for the demand. In this algorithm, we also apply the mode based on probability; that is, the node with shorter distance shall get higher probability to become the access point for node optimization. For the set{Nj∣ENjci≠0} of the given to-be-optimized demand ci and all nodes connected with it, the probability that Nj is selected as the access point is formulated as follows: (15)pNj=1/cNjci∑ENci≠01/cNci(b) Constructing route: after the access point is selected, we shall reconstruct the route. During the route construction, we should consider the influence of different access points. When the selected access point is the customer node, in fact, it is unnecessary to reconstruct a new route; it is only required to add an edge based on the route of the access point. Three conditions need to be met: the load is sufficient; the time constraint is satisfied; the access point is the end point of the vehicle path. If and only if the three conditions are satisfied, we add the edge between the access point and demand node on the vehicle route of the access point, which is regarded as the new route of the demand. When the selected access point is the transfer station node, the route shall be constructed in the following steps.At first, we assign probability to the transfer station node. Similar to the preceding method, the route construction also applies probability-based method. The probability that transfer station nodesj is selected is formulated as follows: (16)p=1/ccisj∗1+pSj/2∑Ecisk≠01/ccisk∗1+psk/2=1+pSj∗∑Ecisk≠0cciskccisj∗∑Ecisk≠01+pskIn Formula (16), psj means the probability that the node Sj is chosen for optimization. (1+psj)/2 is used to standardize the probability and balance the probability of all nodes. It can avoid the local optimization where the top of the formula means the result of standard probability that the node can be chosen. The bottom of the formula means the total result of all nodes, so the formula presents the final probability that the note Sj can be chosen.Formula (16) shows that the probability that transfer station node sj is inversely proportional to the distance between ci and sj. In addition, since transfer station’s fixed cost occupies a fairly large proportion in the comprehensive cost, we reduce by half the probability that unused transfer station node is selected. In random method, we are able to select transfer station node. Similarly, after selection of transfer station node, we can also search for distribution center node and finally construct a complete route.


#### 3.3.4. Acceptance of Route

After a route is completed, we should judge whether the route is acceptable. Firstly, the route must satisfy the load and time constraints; if not, such route should be given up. Secondly, after the route is constructed, if the total cost is better than the current solution, such current solution shall be replaced. However, we should not simply give up the new solution whose total cost fails to be better than the current solution. Otherwise, the algorithm may fall into local optimum. Here we assign accepted probability to the new route. In the algorithm operating process, with the increase of iterations, the current solution is continuously approaching the optimal solution. Thus, in accordance with the relations between current iterations and maximum iterations, we define acceptance probability as follows:(17)λ=IImax∗αIn the above formula,λ means acceptance probability, I means current iterations, Imax means maximum iterations of the system, and α is a constant, which means the datum acceptance probability of the system. On the basis of probability, if the system accepts the route mentioned above, the current solution shall be replaced; if not, this route shall be given up.


#### 3.3.5. Improvement Process

Algorithm3 shows the pseudocodes in the improvement process.Algorithm 3:RandomRePath.
Input: Two-echelon routing graphG, Current solution V Output: New solution V′ 1. Compute the initial value for every node by Formula (15) 2. for 1 tomaxIterationNum do 3.c = ChooseCustomer(); 4.n = ChooseLinkNode(c); 5.if (n ∈ C) 6.AddPath(c,n); 7.if(acceptPath() == true) 8.uptate(V); 9.Endif; 10.else 11.RePath(c,n); 12.if(acceptPath() == true) 13.uptate(V); 14.Endif; 15.Endif; 16.Endfor; 14. returnV;In Algorithm3, line 1 computes the initial value of all nodes. The optimization iteration is executed in lines 2-16. In line 3, the function ChooseCustomer() chooses a customer to be optimized randomly as described in Section 3.3.2. The function ChooseLinkNode() chooses the previous link node based on the connection graph in line 4, as described in Section 3.3.3. In line 5, if the previous link node is a customer node, the algorithm constructs a new path and compares the result with the current result by function acceptPath(). If acceptPath() returns “true,” the algorithm takes the new result as the current result. In lines 10-14, if the previous link node is a satellite node, the algorithm uses function RePath() to construct a new result and compares it to decide whether to accept it as described in Section 3.3.4.


## 3.1. Algorithm Framework

Since the above problems are NP-hard problems, this paper has designed a Heuristic Algorithm Based on Probability (HABP) for solution. The main ideas of the algorithm are the following:(a) In line 1, HABP seeks initial solution by function InitialSolution(G), which takes the Graph G as input. InitialSolution(G) computes a solution by FirstFit(FF) mechanism and return.(b) Weight number distribution: in line 2, according to the node-to-node comprehensive space-time distance consisting of single-vehicle transportation cost and transportation time, the function STDistance(node) makes initial processing of graphG in the system to distribute initial value for each node.(c) Iterative solution: in lines 3 to 5, based on the current solution and initial value, the HABP chooses customer nodes for optimization and replaces the current solution when finally better solution is found, where the function RandomRePath() is the core function of the algorithm, which seeks better solution in the iteration. RandomRePath() will be explained in detail in Section3.3.(d) Algorithm termination: when the number of iterations reaches the designated value, stop the search.The pseudocodes of HABP algorithm framework are shown in Algorithm1.Algorithm 1:HABP algorithm framework.
Input: two-echelon routing graphG Output: New solutionV 1.V = InitialSolution(G); 2. for every node inG: value(node) = STDistance(node); 3. for 1 tomaxIterationNum do 4.V=  RandomRePath(); 5. endfor; 6. returnV;


## 3.2. Composition of Initial Solution

A given two-echelon routing graphG with time constraints consists of m distribution centers, n transfer stations, and k customer demands. At the beginning of heuristic algorithm, we search for an initial solution one by one at first and then take it as the starting point for optimizing. For the algorithm in this paper, we adopt FirstFit (FF) algorithm to generate initial solution.The main ideas of FF algorithm are as follows: for the demand of each customer, we successively search for all potential routes. The first route satisfying the requirements is selected as the route for the demand. Here, we put aside the concentration of goods, assuming that all goods directly start off from the distribution center, arrive at the transfer station, and finally get to the customer. The optimizing is done in the later process of algorithm. In the process of searching for initial route, the two conditions should be met:(a) As for routepci→sci→ci satisfying any customer demand ci, this route is interconnected; that is, Epcisci≠0,Escici≠0.(b) This route transportation can meet the time requirement of customer demand; that is,tpcisci+tscici≤di.The algorithm of initial solution is shown in Algorithm2.Algorithm 2:FirstFit.
Input: two-echelon routing graphG Output: An initial solutionV 1. for every customer demand do 2.for i = 1 to m do 3.for j = 1 to n do 4.if(Epcisci≠0∧Escici≠0∧tpcisci+tscici≤di) 5.rpcisci+=1;rscici+=1; 6.break; 7. returnV;


## 3.3. Optimizing Strategy

For a given two-echelon routing graphG and initial solution, it is required to utilize optimizing strategy to seek better feasible solution and replace current solution. In the optimizing process, heuristic algorithm usually adopts random method to generate new solution and compare it with current solution. In this process, random strategy directly influences algorithm convergence rate. In HABP, problem solution means a vehicle operating plan satisfying all customer demands in G in transportation quantity and time. For two different solutions, on the condition that the time demand of all customers is satisfied, the solution with less total cost is the better one. In the process of seeking better solution, in order to improve convergence rate, it is required to choose a random value with bigger probability to obtain better solution when new solution is generated. Therefore, in the optimizing process, HABP puts forward the concept of selective probability and puts forward an optimizing mode based on probability. This mode contains the following steps:Step 1.
Allocate initial value for all nodes.Step 2.
Distribute probability for all demands and select optimizing node according to probability, based on the initial value and current solution of nodes.Step 3.
Replan for the selected node, generate new solution, and accept new solution with certain probability according to the comparison of new and old solutions.Step 4.
Return the optimal solution obtained in the process as a result.


### 3.3.1. Initial Value of Node

In the two-echelon location routing environment, it is not completely equivalent between nodes. For example, for a customer demand, there are two transfer station nodess1 and s2 for choice, and the two routes have the same cost, but s1 has better connectivity than s2; that is, s1 is connected with more other nodes and has better cost; it is obvious that s1 is a better choice for this demand, for it has more possibility to share the transportation resource with other demands and bring optimization of total cost. Here, connectivity is measured from the two aspects of cost and time. About initial value, we consider that the node with less average cost and average time shall have higher probability for route accession.The cost and time initial values of various nodes are calculated in the formula below.About distribution center nodepi, its initial value is the mean value of the route of all transfer stations connected with it; that is, (7)vcpi=∑Episj≠0cpisjnumpi,vtpi=∑Episj≠0tpisjnumpiIn the formula,numpi means the number of transfer station nodes with route to connect node pi; vcpi means the initial cost of distribution center node pi, whose value is the mean cost value of the route connecting transfer stations; similarly vtpi means the initial time of distribution center node pi, whose value is the mean time value of the route connecting transfer stations.Correspondingly, for distribution center nodesi, its initial value is the mean value of the route of all distribution centers connected with it; that is,(8)vcsi=∑Esipj≠0csipjnumsi,vtsi=∑Esipj≠0tsipjnumsiFor distribution center nodeci, its initial values and mean values of all potential direct routes include the values from ci to transfer station node and from transfer station node to distribution center, since the values from transfer station node to distribution center can be expressed as vcsi and vtsi, which are formulated as follows: (9)vcci=∑Ecisj≠0ccisj+vcsjnumci,vtci=∑Ecisj≠0tcisj+vtsjnumciIt should be noted that, for transfer station node and customer node, the initial value does not take the path between the same nodes into account. This is because the path between the same nodes does not represent the actual capacity of the node but only reflects the convenience in combining goods.


### 3.3.2. Selection of Optimizing Nodes

For a given graphG and a current solution V, V consists of distribution center utility state set UP, transfer center utility state set US, and route set R. The optimizing process means seeking a new solution with cost and time advantages by adjusting the content of current solution. For the algorithm in this paper, the optimizing process consists of two parts: determination of optimizing node and construction of new solution, in which selection of optimizing node is to determine the adjustment of current solution. In this algorithm, we assign probability to each customer demand, take random selection of customer node, and replan its route for optimization, based on node initial value and current solution.At first, we should define the cost of all demands in the current solution. For any demand, finishing the transportation may include four kinds of routes:(1)P→S; (2)S→S; (3)S→C; and (4)C→C. Among the four routes, route 1 and route 3 are the unique paths that the goods must pass, while route 2 and route 4 are not unique. For any of the routes, the vehicle may carry goods from multiple customers. Therefore, the transportation costs are shared by multiple customers. Those costs are accounted on one-route vehicle T; proportion of demand ci is formulated as λi:(10)λi=wiWTIn the above formula,WT means load of vehicle T on this route. For demand ci of pathci, its costi can be formulated as follows: (11)costi=∑EN1N2∈pathciλi∗cN1N2Considering the time cost ofci, in accordance with the definition, the transfer of all vehicles only occurs in the transfer from the first-echelon transportation to the second-echelon transportation; therefore, ci transportation shall and shall only be realized through the two vehicles Ti1 and Ti1, with the total time as follows: (12)ti=tTi1sci+tTi2ci+wtciAccording to the above parameters, we can assign the probabilities that are selected for optimization to different demands. We should make an analysis on the aspects of cost and time. As for the cost, a demand that requires high cost in the current solution, it must get distribution of higher probability, for such demand may have higher possibility for optimization. As a reference, the ratio of initial value to current costcosti/vcci can be considered as one of basic parameters for probability assignment. As for the time, any demand must meet the requirement of deadline; any delivery within the deadline shall be consistent in effect. However, generally, it is hoped to realize the delivery ahead of the deadline, under the condition that the cost is not raised. In fact, it can also improve the customer satisfaction. In this paper, we define the time impact factor of ci as χi, which is formulated as follows: (13)χi=1ti≤vtcitivtciti>vtciIn formula (13), when the time consumption of demand ci is higher than the mean value in the current solution, a time impact factor (>1) is assigned to it, so that it has higher probability to be selected for optimization.To sum up, in the HABP algorithm, the probabilitypi that a demand ci is selected for optimization is formulated as follows: (14)pi=costi∗χi∑cj∈Ccostj∗χj


### 3.3.3. Construction and Acceptance Strategy of New Solution

After all demands are assigned with probabilities, we make a random selection of a demand for optimization on the probability base. The optimizing method is to cancel existing route and reselect a new route for this demand. In the selection process, the following factors should be considered:(a) Access point: for the route of a demandci, in the final section of the route, two possibilities may occur: starting from the transfer station node or starting from the customer node. It is clear that, in view of cost, “vehicle pooling” is useful for control of comprehensive cost. Therefore, during the planning of new route, “vehicle pooling” should be realized as much as possible under the premise that time and load constraints are not exceeded.(b) Treatment of original route: for the generation of new route, it is required to cancel the original route. Considering the following conditions, the vehicle transporting for a certain demandc1 on the second-echelon route has a route in the current solution, expressed as {s0,c1,c2}. When c1 is selected as the optimization node, route Es0c1 becomes unnecessary; at this time, if Es0c2≠0, we should make a comparison of the cost and time results between original route {s0,c1,c2} and new route {s0,c2}; the better one shall be selected.(c) Load: during the construction of new route, load is an important factor for consideration. At first, it should be ensured that no overload occurs. Here, “overload” includes overload of all nodes and overload of all vehicles in the whole transportation. For the overload route, probability should not be assigned. With no overload, we hope to select the route loaded with bigger loading capacity, because this means a reduction of unit transportation cost.(d) Time: the time consumption must be lower than the time constraint of the demand.(e) Selection of transfer station and distribution center nodes. During the construction of new route, we should select the transfer station and distribution center route which can bring better solution.Considering the above constraints, the construction of new route can be implemented in the following steps:(a) Selecting access point: at first, we should select access point for the demand. In this algorithm, we also apply the mode based on probability; that is, the node with shorter distance shall get higher probability to become the access point for node optimization. For the set{Nj∣ENjci≠0} of the given to-be-optimized demand ci and all nodes connected with it, the probability that Nj is selected as the access point is formulated as follows: (15)pNj=1/cNjci∑ENci≠01/cNci(b) Constructing route: after the access point is selected, we shall reconstruct the route. During the route construction, we should consider the influence of different access points. When the selected access point is the customer node, in fact, it is unnecessary to reconstruct a new route; it is only required to add an edge based on the route of the access point. Three conditions need to be met: the load is sufficient; the time constraint is satisfied; the access point is the end point of the vehicle path. If and only if the three conditions are satisfied, we add the edge between the access point and demand node on the vehicle route of the access point, which is regarded as the new route of the demand. When the selected access point is the transfer station node, the route shall be constructed in the following steps.At first, we assign probability to the transfer station node. Similar to the preceding method, the route construction also applies probability-based method. The probability that transfer station nodesj is selected is formulated as follows: (16)p=1/ccisj∗1+pSj/2∑Ecisk≠01/ccisk∗1+psk/2=1+pSj∗∑Ecisk≠0cciskccisj∗∑Ecisk≠01+pskIn Formula (16), psj means the probability that the node Sj is chosen for optimization. (1+psj)/2 is used to standardize the probability and balance the probability of all nodes. It can avoid the local optimization where the top of the formula means the result of standard probability that the node can be chosen. The bottom of the formula means the total result of all nodes, so the formula presents the final probability that the note Sj can be chosen.Formula (16) shows that the probability that transfer station node sj is inversely proportional to the distance between ci and sj. In addition, since transfer station’s fixed cost occupies a fairly large proportion in the comprehensive cost, we reduce by half the probability that unused transfer station node is selected. In random method, we are able to select transfer station node. Similarly, after selection of transfer station node, we can also search for distribution center node and finally construct a complete route.


### 3.3.4. Acceptance of Route

After a route is completed, we should judge whether the route is acceptable. Firstly, the route must satisfy the load and time constraints; if not, such route should be given up. Secondly, after the route is constructed, if the total cost is better than the current solution, such current solution shall be replaced. However, we should not simply give up the new solution whose total cost fails to be better than the current solution. Otherwise, the algorithm may fall into local optimum. Here we assign accepted probability to the new route. In the algorithm operating process, with the increase of iterations, the current solution is continuously approaching the optimal solution. Thus, in accordance with the relations between current iterations and maximum iterations, we define acceptance probability as follows:(17)λ=IImax∗αIn the above formula,λ means acceptance probability, I means current iterations, Imax means maximum iterations of the system, and α is a constant, which means the datum acceptance probability of the system. On the basis of probability, if the system accepts the route mentioned above, the current solution shall be replaced; if not, this route shall be given up.


### 3.3.5. Improvement Process

Algorithm3 shows the pseudocodes in the improvement process.Algorithm 3:RandomRePath.
Input: Two-echelon routing graphG, Current solution V Output: New solution V′ 1. Compute the initial value for every node by Formula (15) 2. for 1 tomaxIterationNum do 3.c = ChooseCustomer(); 4.n = ChooseLinkNode(c); 5.if (n ∈ C) 6.AddPath(c,n); 7.if(acceptPath() == true) 8.uptate(V); 9.Endif; 10.else 11.RePath(c,n); 12.if(acceptPath() == true) 13.uptate(V); 14.Endif; 15.Endif; 16.Endfor; 14. returnV;In Algorithm3, line 1 computes the initial value of all nodes. The optimization iteration is executed in lines 2-16. In line 3, the function ChooseCustomer() chooses a customer to be optimized randomly as described in Section 3.3.2. The function ChooseLinkNode() chooses the previous link node based on the connection graph in line 4, as described in Section 3.3.3. In line 5, if the previous link node is a customer node, the algorithm constructs a new path and compares the result with the current result by function acceptPath(). If acceptPath() returns “true,” the algorithm takes the new result as the current result. In lines 10-14, if the previous link node is a satellite node, the algorithm uses function RePath() to construct a new result and compares it to decide whether to accept it as described in Section 3.3.4.


## 3.3.1. Initial Value of Node

In the two-echelon location routing environment, it is not completely equivalent between nodes. For example, for a customer demand, there are two transfer station nodess1 and s2 for choice, and the two routes have the same cost, but s1 has better connectivity than s2; that is, s1 is connected with more other nodes and has better cost; it is obvious that s1 is a better choice for this demand, for it has more possibility to share the transportation resource with other demands and bring optimization of total cost. Here, connectivity is measured from the two aspects of cost and time. About initial value, we consider that the node with less average cost and average time shall have higher probability for route accession.The cost and time initial values of various nodes are calculated in the formula below.About distribution center nodepi, its initial value is the mean value of the route of all transfer stations connected with it; that is, (7)vcpi=∑Episj≠0cpisjnumpi,vtpi=∑Episj≠0tpisjnumpiIn the formula,numpi means the number of transfer station nodes with route to connect node pi; vcpi means the initial cost of distribution center node pi, whose value is the mean cost value of the route connecting transfer stations; similarly vtpi means the initial time of distribution center node pi, whose value is the mean time value of the route connecting transfer stations.Correspondingly, for distribution center nodesi, its initial value is the mean value of the route of all distribution centers connected with it; that is,(8)vcsi=∑Esipj≠0csipjnumsi,vtsi=∑Esipj≠0tsipjnumsiFor distribution center nodeci, its initial values and mean values of all potential direct routes include the values from ci to transfer station node and from transfer station node to distribution center, since the values from transfer station node to distribution center can be expressed as vcsi and vtsi, which are formulated as follows: (9)vcci=∑Ecisj≠0ccisj+vcsjnumci,vtci=∑Ecisj≠0tcisj+vtsjnumciIt should be noted that, for transfer station node and customer node, the initial value does not take the path between the same nodes into account. This is because the path between the same nodes does not represent the actual capacity of the node but only reflects the convenience in combining goods.


## 3.3.2. Selection of Optimizing Nodes

For a given graphG and a current solution V, V consists of distribution center utility state set UP, transfer center utility state set US, and route set R. The optimizing process means seeking a new solution with cost and time advantages by adjusting the content of current solution. For the algorithm in this paper, the optimizing process consists of two parts: determination of optimizing node and construction of new solution, in which selection of optimizing node is to determine the adjustment of current solution. In this algorithm, we assign probability to each customer demand, take random selection of customer node, and replan its route for optimization, based on node initial value and current solution.At first, we should define the cost of all demands in the current solution. For any demand, finishing the transportation may include four kinds of routes:(1)P→S; (2)S→S; (3)S→C; and (4)C→C. Among the four routes, route 1 and route 3 are the unique paths that the goods must pass, while route 2 and route 4 are not unique. For any of the routes, the vehicle may carry goods from multiple customers. Therefore, the transportation costs are shared by multiple customers. Those costs are accounted on one-route vehicle T; proportion of demand ci is formulated as λi:(10)λi=wiWTIn the above formula,WT means load of vehicle T on this route. For demand ci of pathci, its costi can be formulated as follows: (11)costi=∑EN1N2∈pathciλi∗cN1N2Considering the time cost ofci, in accordance with the definition, the transfer of all vehicles only occurs in the transfer from the first-echelon transportation to the second-echelon transportation; therefore, ci transportation shall and shall only be realized through the two vehicles Ti1 and Ti1, with the total time as follows: (12)ti=tTi1sci+tTi2ci+wtciAccording to the above parameters, we can assign the probabilities that are selected for optimization to different demands. We should make an analysis on the aspects of cost and time. As for the cost, a demand that requires high cost in the current solution, it must get distribution of higher probability, for such demand may have higher possibility for optimization. As a reference, the ratio of initial value to current costcosti/vcci can be considered as one of basic parameters for probability assignment. As for the time, any demand must meet the requirement of deadline; any delivery within the deadline shall be consistent in effect. However, generally, it is hoped to realize the delivery ahead of the deadline, under the condition that the cost is not raised. In fact, it can also improve the customer satisfaction. In this paper, we define the time impact factor of ci as χi, which is formulated as follows: (13)χi=1ti≤vtcitivtciti>vtciIn formula (13), when the time consumption of demand ci is higher than the mean value in the current solution, a time impact factor (>1) is assigned to it, so that it has higher probability to be selected for optimization.To sum up, in the HABP algorithm, the probabilitypi that a demand ci is selected for optimization is formulated as follows: (14)pi=costi∗χi∑cj∈Ccostj∗χj


## 3.3.3. Construction and Acceptance Strategy of New Solution

After all demands are assigned with probabilities, we make a random selection of a demand for optimization on the probability base. The optimizing method is to cancel existing route and reselect a new route for this demand. In the selection process, the following factors should be considered:(a) Access point: for the route of a demandci, in the final section of the route, two possibilities may occur: starting from the transfer station node or starting from the customer node. It is clear that, in view of cost, “vehicle pooling” is useful for control of comprehensive cost. Therefore, during the planning of new route, “vehicle pooling” should be realized as much as possible under the premise that time and load constraints are not exceeded.(b) Treatment of original route: for the generation of new route, it is required to cancel the original route. Considering the following conditions, the vehicle transporting for a certain demandc1 on the second-echelon route has a route in the current solution, expressed as {s0,c1,c2}. When c1 is selected as the optimization node, route Es0c1 becomes unnecessary; at this time, if Es0c2≠0, we should make a comparison of the cost and time results between original route {s0,c1,c2} and new route {s0,c2}; the better one shall be selected.(c) Load: during the construction of new route, load is an important factor for consideration. At first, it should be ensured that no overload occurs. Here, “overload” includes overload of all nodes and overload of all vehicles in the whole transportation. For the overload route, probability should not be assigned. With no overload, we hope to select the route loaded with bigger loading capacity, because this means a reduction of unit transportation cost.(d) Time: the time consumption must be lower than the time constraint of the demand.(e) Selection of transfer station and distribution center nodes. During the construction of new route, we should select the transfer station and distribution center route which can bring better solution.Considering the above constraints, the construction of new route can be implemented in the following steps:(a) Selecting access point: at first, we should select access point for the demand. In this algorithm, we also apply the mode based on probability; that is, the node with shorter distance shall get higher probability to become the access point for node optimization. For the set{Nj∣ENjci≠0} of the given to-be-optimized demand ci and all nodes connected with it, the probability that Nj is selected as the access point is formulated as follows: (15)pNj=1/cNjci∑ENci≠01/cNci(b) Constructing route: after the access point is selected, we shall reconstruct the route. During the route construction, we should consider the influence of different access points. When the selected access point is the customer node, in fact, it is unnecessary to reconstruct a new route; it is only required to add an edge based on the route of the access point. Three conditions need to be met: the load is sufficient; the time constraint is satisfied; the access point is the end point of the vehicle path. If and only if the three conditions are satisfied, we add the edge between the access point and demand node on the vehicle route of the access point, which is regarded as the new route of the demand. When the selected access point is the transfer station node, the route shall be constructed in the following steps.At first, we assign probability to the transfer station node. Similar to the preceding method, the route construction also applies probability-based method. The probability that transfer station nodesj is selected is formulated as follows: (16)p=1/ccisj∗1+pSj/2∑Ecisk≠01/ccisk∗1+psk/2=1+pSj∗∑Ecisk≠0cciskccisj∗∑Ecisk≠01+pskIn Formula (16), psj means the probability that the node Sj is chosen for optimization. (1+psj)/2 is used to standardize the probability and balance the probability of all nodes. It can avoid the local optimization where the top of the formula means the result of standard probability that the node can be chosen. The bottom of the formula means the total result of all nodes, so the formula presents the final probability that the note Sj can be chosen.Formula (16) shows that the probability that transfer station node sj is inversely proportional to the distance between ci and sj. In addition, since transfer station’s fixed cost occupies a fairly large proportion in the comprehensive cost, we reduce by half the probability that unused transfer station node is selected. In random method, we are able to select transfer station node. Similarly, after selection of transfer station node, we can also search for distribution center node and finally construct a complete route.


## 3.3.4. Acceptance of Route

After a route is completed, we should judge whether the route is acceptable. Firstly, the route must satisfy the load and time constraints; if not, such route should be given up. Secondly, after the route is constructed, if the total cost is better than the current solution, such current solution shall be replaced. However, we should not simply give up the new solution whose total cost fails to be better than the current solution. Otherwise, the algorithm may fall into local optimum. Here we assign accepted probability to the new route. In the algorithm operating process, with the increase of iterations, the current solution is continuously approaching the optimal solution. Thus, in accordance with the relations between current iterations and maximum iterations, we define acceptance probability as follows:(17)λ=IImax∗αIn the above formula,λ means acceptance probability, I means current iterations, Imax means maximum iterations of the system, and α is a constant, which means the datum acceptance probability of the system. On the basis of probability, if the system accepts the route mentioned above, the current solution shall be replaced; if not, this route shall be given up.


## 3.3.5. Improvement Process

Algorithm3 shows the pseudocodes in the improvement process.Algorithm 3:RandomRePath.
Input: Two-echelon routing graphG, Current solution V Output: New solution V′ 1. Compute the initial value for every node by Formula (15) 2. for 1 tomaxIterationNum do 3.c = ChooseCustomer(); 4.n = ChooseLinkNode(c); 5.if (n ∈ C) 6.AddPath(c,n); 7.if(acceptPath() == true) 8.uptate(V); 9.Endif; 10.else 11.RePath(c,n); 12.if(acceptPath() == true) 13.uptate(V); 14.Endif; 15.Endif; 16.Endfor; 14. returnV;In Algorithm3, line 1 computes the initial value of all nodes. The optimization iteration is executed in lines 2-16. In line 3, the function ChooseCustomer() chooses a customer to be optimized randomly as described in Section 3.3.2. The function ChooseLinkNode() chooses the previous link node based on the connection graph in line 4, as described in Section 3.3.3. In line 5, if the previous link node is a customer node, the algorithm constructs a new path and compares the result with the current result by function acceptPath(). If acceptPath() returns “true,” the algorithm takes the new result as the current result. In lines 10-14, if the previous link node is a satellite node, the algorithm uses function RePath() to construct a new result and compares it to decide whether to accept it as described in Section 3.3.4.


## 4. Computational Experiments

In this section, we present the results of the computational testing. To the best of our knowledge, there is no suitable standard benchmark for the problem considered here. Therefore, the main experiments were carried out on a number of randomly generated test instances, specially designed to reflect the sizes and peculiarities of real-world distribution problems.


### 4.1. The Generation of the Instances

The main data generated by the instances include (a) the load and fixed cost of the distribution center, (b) the load and fixed cost of the satellite, (c) the weight and time constraints of the customer, and (d) the connected graph G, where each element contains the distance between two nodes and the transportation time.Other major parameters that need to be set autonomously in the experiment include (a) the number of parking lots P, (b) the number of satellites S, (c) the number of customers C, (d) the capacity of the first-floor vehicles which is set as 1500, (e) the capacity of the second layer vehicle which is set as 900, and (f) the number of algorithm iterations which is set as 1000 times.It should be noted that the scale of the experiment in this paper is set from 0 to 500. Generally speaking, for any optimization algorithm, one must supply some termination conditions specifying when the algorithm halts. In the experiment of this paper, we have tested that, after 1000 iterations of the system, the optimization effect becomes weaker. That is, it can obtain a better solution. In addition, the system sets a uniform number of iterations, which can eliminate noises in different iterations and pay more attention to the optimization effect of the algorithm itself. Therefore 1000 is set as the number of iterations. When the system needs to deal with large-scale problems, it only needs to add the termination condition to the algorithm. The system will automatically terminate when the convergence speed is reduced. All randomly generated values follow the following methods:(18)value=a+b∗rand∗swherea is the lowest value and rand() is the random value, whose values are real numbers between 0 and 1. b represents the coefficients of two random values; s denotes the scaling factor which is used to expand the random range and calculated by s=c±d∗rand(), where c and d represent the minimum size coefficient and the scale factor parameter, respectively. Therefore, the value of s is a real number whose value is between c-d and c+d. When c=1 and d=0, it means that the value does not need to be scaled. The values of the parameters are shown in Table 1.Table 1
Parameter values for generating instances.
Value a b c d p l i 2000 2000 k/m 0 p c i 1000 1000 1 0.15 s l i 1000 1000 k/n 0 s c i 500 1000 1 0.15 w i 50 80 1 0 d i 150 170 1 0 c N 1 N 2 30 30 1 0 t N 1 N 2 20 30 1 0.2


### 4.2. Experiments and Results

In this paper, the purpose of proposing the HABP algorithm is to solve problems of insufficient cost optimization and long calculation time. Therefore, this paper mainly compares the change of cost optimization and algorithm running time after using the HABP algorithm. In the simulation experiment, we analyzed the total cost and running time of HABP when using the FirstFit algorithm and the HABP algorithm in instances of different scales. In the algorithm, the following conditions need to be met:(a)
The capacity constraints of the vehicle(b)
Load constraints of the distribution center(c)
The load constraints of the transfer station(d)
Customer's time constraints(e)
Path access constraintsAll tests were performed on an Intel i5-3230M 2.6GHz processor. The algorithm is implemented in C++. The results of the algorithms are shown in Table2.Table 2
Results of the algorithms.
Computational instances Number of distribution centers Number of transfer stations Number of customers Name of instances Finalcosts Runtime  (s) Initialsolution  cost Optimization(%) 1 3 8 20 P3-S8-C10 3051.93 0.09 4473.99 31.785 2 3 8 50 P3-S8-C50 6078.23 0.12 7866.92 22.737 3 3 8 100 P3-S8-C100 9274.95 0.16 12489.2 25.737 4 3 8 150 P3-S8-C150 12799.8 0.19 17820.7 28.174 5 3 8 200 P3-S8-C200 14959.7 0.22 19332.7 22.580 6 3 8 250 P3-S8-C250 19391.5 0.23 26043.9 25.543 7 3 8 300 P3-S8-C300 23544.5 0.23 31092 26.198 8 3 8 350 P3-S8-C350 25819.1 0.27 32973.6 21.698 9 3 8 400 P3-S8-C400 31814.9 0.26 41883.6 24.04 10 3 8 450 P3-S8-C450 34776.4 0.29 46474.6 25.171 11 3 8 500 P3-S8-C500 35206.1 0.31 44403.3 20.713 12 3 8 550 P3-S8-C550 43403.2 0.33 56309.2 22.92 13 3 8 600 P3-S8-C600 44031.1 0.35 56052 21.446 14 3 8 650 P3-S8-C650 51651.6 0.34 67129.7 23.057 15 3 8 700 P3-S8-C700 56875.8 0.35 72771.3 21.843 16 3 8 750 P3-S8-C750 61195.3 0.36 77920.5 21.464 17 3 8 800 P3-S8-C800 62861 0.39 79422.5 20.852 18 3 8 850 P3-S8-C850 69820.9 0.41 89022.5 21.569 19 3 8 900 P3-S8-C900 68511.9 0.43 85313.8 19.694 20 3 8 950 P3-S8-C950 78655.7 0.44 97806.6 19.581 21 3 8 1000 P3-S8-C1000 83872.3 0.46 104377 19.645In the actual environment, the number of distribution centers or transfer stations is normally fixed, while the number of customers usually varies. Therefore, in the simulation experiment, the value of C is variable, ranging from 50 to 1000 with a span of 50. From the experimental results, the following conclusions can be drawn.(a) Compared with the FirstFit algorithm, HABP can effectively optimize the cost, and the optimized ratio is over 20%. What is worth mentioning is that the FirstFit (FF) algorithm is a typical optimization algorithm in computer operating system, which is widely used in industry and is thought to be efficient. The simulation that used FF as the comparison object is meaningful. The experiment was based on the statistical result. We generated multiple graphs of given scale, optimized the result, and took the average result as the final result.(b) The optimization result of HABP is slightly reduced when the size of customer increases, but the decrease is not obvious. In the experiment, the effect of the number of iterations on the experimental results is more reflected. This is because the setting is 1000 times. Under the premise of the number of iterations, the larger the scale, the lower the coverage of the randomly selected optimization node.(c) The convergence speed of HABP is fast. When the number of iterations is determined, similar optimization values can be obtained for different scale instances, indicating that the algorithm can obtain a relatively better solution within a small number of iterations.(d) HABP has extremely fast operation speed. This speed is not affected by scale. Therefore, a significant advantage of HABP is that it can be applied to large-scale calculations, which cannot be realized by most current researches and algorithms.


## 4.1. The Generation of the Instances

The main data generated by the instances include (a) the load and fixed cost of the distribution center, (b) the load and fixed cost of the satellite, (c) the weight and time constraints of the customer, and (d) the connected graph G, where each element contains the distance between two nodes and the transportation time.Other major parameters that need to be set autonomously in the experiment include (a) the number of parking lots P, (b) the number of satellites S, (c) the number of customers C, (d) the capacity of the first-floor vehicles which is set as 1500, (e) the capacity of the second layer vehicle which is set as 900, and (f) the number of algorithm iterations which is set as 1000 times.It should be noted that the scale of the experiment in this paper is set from 0 to 500. Generally speaking, for any optimization algorithm, one must supply some termination conditions specifying when the algorithm halts. In the experiment of this paper, we have tested that, after 1000 iterations of the system, the optimization effect becomes weaker. That is, it can obtain a better solution. In addition, the system sets a uniform number of iterations, which can eliminate noises in different iterations and pay more attention to the optimization effect of the algorithm itself. Therefore 1000 is set as the number of iterations. When the system needs to deal with large-scale problems, it only needs to add the termination condition to the algorithm. The system will automatically terminate when the convergence speed is reduced. All randomly generated values follow the following methods:(18)value=a+b∗rand∗swherea is the lowest value and rand() is the random value, whose values are real numbers between 0 and 1. b represents the coefficients of two random values; s denotes the scaling factor which is used to expand the random range and calculated by s=c±d∗rand(), where c and d represent the minimum size coefficient and the scale factor parameter, respectively. Therefore, the value of s is a real number whose value is between c-d and c+d. When c=1 and d=0, it means that the value does not need to be scaled. The values of the parameters are shown in Table 1.Table 1
Parameter values for generating instances.
Value a b c d p l i 2000 2000 k/m 0 p c i 1000 1000 1 0.15 s l i 1000 1000 k/n 0 s c i 500 1000 1 0.15 w i 50 80 1 0 d i 150 170 1 0 c N 1 N 2 30 30 1 0 t N 1 N 2 20 30 1 0.2


## 4.2. Experiments and Results

In this paper, the purpose of proposing the HABP algorithm is to solve problems of insufficient cost optimization and long calculation time. Therefore, this paper mainly compares the change of cost optimization and algorithm running time after using the HABP algorithm. In the simulation experiment, we analyzed the total cost and running time of HABP when using the FirstFit algorithm and the HABP algorithm in instances of different scales. In the algorithm, the following conditions need to be met:(a)
The capacity constraints of the vehicle(b)
Load constraints of the distribution center(c)
The load constraints of the transfer station(d)
Customer's time constraints(e)
Path access constraintsAll tests were performed on an Intel i5-3230M 2.6GHz processor. The algorithm is implemented in C++. The results of the algorithms are shown in Table2.Table 2
Results of the algorithms.
Computational instances Number of distribution centers Number of transfer stations Number of customers Name of instances Finalcosts Runtime  (s) Initialsolution  cost Optimization(%) 1 3 8 20 P3-S8-C10 3051.93 0.09 4473.99 31.785 2 3 8 50 P3-S8-C50 6078.23 0.12 7866.92 22.737 3 3 8 100 P3-S8-C100 9274.95 0.16 12489.2 25.737 4 3 8 150 P3-S8-C150 12799.8 0.19 17820.7 28.174 5 3 8 200 P3-S8-C200 14959.7 0.22 19332.7 22.580 6 3 8 250 P3-S8-C250 19391.5 0.23 26043.9 25.543 7 3 8 300 P3-S8-C300 23544.5 0.23 31092 26.198 8 3 8 350 P3-S8-C350 25819.1 0.27 32973.6 21.698 9 3 8 400 P3-S8-C400 31814.9 0.26 41883.6 24.04 10 3 8 450 P3-S8-C450 34776.4 0.29 46474.6 25.171 11 3 8 500 P3-S8-C500 35206.1 0.31 44403.3 20.713 12 3 8 550 P3-S8-C550 43403.2 0.33 56309.2 22.92 13 3 8 600 P3-S8-C600 44031.1 0.35 56052 21.446 14 3 8 650 P3-S8-C650 51651.6 0.34 67129.7 23.057 15 3 8 700 P3-S8-C700 56875.8 0.35 72771.3 21.843 16 3 8 750 P3-S8-C750 61195.3 0.36 77920.5 21.464 17 3 8 800 P3-S8-C800 62861 0.39 79422.5 20.852 18 3 8 850 P3-S8-C850 69820.9 0.41 89022.5 21.569 19 3 8 900 P3-S8-C900 68511.9 0.43 85313.8 19.694 20 3 8 950 P3-S8-C950 78655.7 0.44 97806.6 19.581 21 3 8 1000 P3-S8-C1000 83872.3 0.46 104377 19.645In the actual environment, the number of distribution centers or transfer stations is normally fixed, while the number of customers usually varies. Therefore, in the simulation experiment, the value of C is variable, ranging from 50 to 1000 with a span of 50. From the experimental results, the following conclusions can be drawn.(a) Compared with the FirstFit algorithm, HABP can effectively optimize the cost, and the optimized ratio is over 20%. What is worth mentioning is that the FirstFit (FF) algorithm is a typical optimization algorithm in computer operating system, which is widely used in industry and is thought to be efficient. The simulation that used FF as the comparison object is meaningful. The experiment was based on the statistical result. We generated multiple graphs of given scale, optimized the result, and took the average result as the final result.(b) The optimization result of HABP is slightly reduced when the size of customer increases, but the decrease is not obvious. In the experiment, the effect of the number of iterations on the experimental results is more reflected. This is because the setting is 1000 times. Under the premise of the number of iterations, the larger the scale, the lower the coverage of the randomly selected optimization node.(c) The convergence speed of HABP is fast. When the number of iterations is determined, similar optimization values can be obtained for different scale instances, indicating that the algorithm can obtain a relatively better solution within a small number of iterations.(d) HABP has extremely fast operation speed. This speed is not affected by scale. Therefore, a significant advantage of HABP is that it can be applied to large-scale calculations, which cannot be realized by most current researches and algorithms.


## 5. Conclusion

In this paper, a mixed-integer linear programming model is put forward for a two-echelon location routing problem with time constraints in city logistics. To solve this problem, we propose a metaheuristic algorithm with comprehensive consideration of time and space accessibility and application of probability selection principle. The methods are tested using randomly generated instances with up to 1000 customers. Experimental results show that the proposed algorithm is effective in terms of quality of solutions and computation times in most of the solved instances. In future research, we aim to consider more real-life constraints to deal with some other realistic problems, such as uncertain demand of customers and drop shipping measure. We are also interested in improving the algorithm with other metaheuristic techniques and an iterative process.



---

*Source: 2549713-2018-12-13.xml*