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

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.


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 smallscale 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 decisionmaking on both strategic and tactical planning, 2E-LRP problem is put forward, but few scholars have carried out

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 Figure 1.The overall system is expressed as graph  = (, ), in which node set  =  ∪  ∪  means set of distribution centers, transfer stations, and customers and  means set of edges.Let  be a set of customer demands.The demand of each customer  (1 ≤  ≤ ) is expressed as twotuple combination   = (  ,   ), in which   ∈ + means this customer's demands quantity;   ∈ + means this customer's required date for receiving the goods before   ; that is, the time window of customer demand is [0,   ].
Distribution center set  consists of  distribution centers; any distribution center   (1 ≤  ≤ ) is expressed as two-tuple combination   = (  ,   ), in which   ∈ + means fixed cost of the distribution center; once the distribution center is used, the fixed cost occurs;   ∈ + means maximum load of the distribution center; the distribution center cannot allow overload; that is, the loading of distribution center cannot exceed   .
Transfer station set  consists of  transfer stations, any transfer station   (1 ≤  ≤ ) is expressed as two-tuple combination   = (  ,   ), in which   ∈ + means fixed cost of transfer station; once the transfer station is used, the fixed cost occurs;   ∈ + means maximum load of the transfer station; that is, the transferred goods quantity of transfer station cannot exceed   .
The transportation between distribution centers and transfer stations is the first-echelon transportation, and the limited capacity of each first-echelon vehicle is represented as  1 .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  2 ; obviously, considering the practical conditions, there is The edge between any two nodes  1 ,   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   may be from different first-echelon transporting vehicles, and the arrivals of these vehicles are not at the same time.Although the goods of   arrive at transfer station    , 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   -used same vehicle to the transfer station during the second-echelon transportation; the waiting time is expressed as . 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 and   -related vehicles and paths, which is expressed as follows.
Set  = { 1 ,  2 , ...,   } is expressed as the use condition of distribution center, in which   (1 ≤  ≤ ) is a 0-1variable.  = 1 means that the distribution center is used; otherwise, the distribution center is not used.
Set  = { 1 ,  2 , ...,   } is expressed as the use condition of transfer station, in which   (1 ≤  ≤ ) is a 0-1variable.  = 1 means that the transfer station is used; otherwise, the transfer station is not used.
Set  = {  2 1 } means the set of all paths in the result, in which 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, A correct plan must satisfy the following constraints: Constraints (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  depends on the echelon on which the vehicle lies.As for the first-echelon vehicle,  =  1 ; otherwise,  =  2 .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.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 Section 3.3.

Algorithm
(d) Algorithm termination: when the number of iterations reaches the designated value, stop the search.
The pseudocodes of HABP algorithm framework are shown in Algorithm 1.

. . Composition of Initial Solution.
A given two-echelon routing graph  with time constraints consists of  distribution centers,  transfer stations, and  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 route    →    →   satisfying any customer demand   , this route is interconnected; that is, (b) This route transportation can meet the time requirement of customer demand; that is, The algorithm of initial solution is shown in Algorithm 2.
. .Optimizing Strategy.For a given two-echelon routing graph  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  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 Input: two-echelon routing graph  Output: An initial solution  1. for every customer demand do 2. for i = 1 to m do 3.
for j = 1 to n do 4. if(

5.
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 .Allocate initial value for all nodes.
Step .Distribute probability for all demands and select optimizing node according to probability, based on the initial value and current solution of nodes.
Step .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 .Return the optimal solution obtained in the process as a result.
. . .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 nodes  1 and  2 for choice, and the two routes have the same cost, but  1 has better connectivity than  2 ; that is,  1 is connected with more other nodes and has better cost; it is obvious that  1 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 node   , its initial value is the mean value of the route of all transfer stations connected with it; that is, In the formula,    means the number of transfer station nodes with route to connect node   ; V   means the initial cost of distribution center node   , whose value is the mean cost value of the route connecting transfer stations; similarly V   means the initial time of distribution center node   , whose value is the mean time value of the route connecting transfer stations.
Correspondingly, for distribution center node   , its initial value is the mean value of the route of all distribution centers connected with it; that is, For distribution center node   , its initial values and mean values of all potential direct routes include the values from   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 V   and V   , which are formulated as follows: It 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.
. . .Selection of Optimizing Nodes.For a given graph  and a current solution ,  consists of distribution center utility state set , transfer center utility state set , and route set .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)  → ; (2)  → ; (3)  → ; and (4)  → .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 ; proportion of demand   is formulated as   : In the above formula,   means load of vehicle  on this route.For demand   of ℎ   , its   can be formulated as follows: Considering the time cost of   , in accordance with the definition, the transfer of all vehicles only occurs in the transfer from the first-echelon transportation to the secondechelon transportation; therefore,   transportation shall and shall only be realized through the two vehicles  1  and  1  , with the total time as follows: According 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 cost   /V   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   as   , which is formulated as follows: In formula (13), when the time consumption of demand   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 probability   that a demand   is selected for optimization is formulated as follows: . . .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 Mathematical Problems in Engineering 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 demand   , 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 demand  1 on the second-echelon route has a route in the current solution, expressed as { 0 ,  1 ,  2 }.When  1 is selected as the optimization node, route   1  0 becomes unnecessary; at this time, if   2  0 ̸ = 0, we should make a comparison of the cost and time results between original route { 0 ,  1 ,  2 } and new route { 0 ,  2 }; 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 {  |      ̸ = 0} of the given to-be-optimized demand   and all nodes connected with it, the probability that   is selected as the access point is formulated as follows: (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 node   is selected is formulated as follows: In Formula ( 16),    means the probability that the node   is chosen for optimization.(1+   )/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   can be chosen.
Formula (16) shows that the probability that transfer station node   is inversely proportional to the distance between   and   .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.
. . .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: In the above formula,  means acceptance probability,  means current iterations,  max 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.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.

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.

Figure 1 :
Figure 1: Example of the two-echelon logistics system.

𝑇 1 𝑖;
(2) time for the second-echelon transportation, expressed as     2  ; and (3) waiting time at transfer station    .The waiting time of vehicle means that, in the process of the first-echelon transportation, the goods to be transported by vehicle  2 . .Algorithm Framework.Since the above problems are NPhard 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.Initial-Solution(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 graph  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 Input: two-echelon routing graph  Output: New solution  1.  = InitialSolution(); 2. for every node in G : value(node) = STDistance(node); 3. for 1 to maxIterationNum do 4.  = RandomRePath(); 5. endfor; 6. return ; Algorithm 1: HABP algorithm framework.
1 ∈ + means the time for finishing the transportation between the two nodes.When  1 ,  2 ∈  ∪ , the first-echelon transportation occurs; when  2 ∈  ∧  1 ∈  ∪ , the second-echelon transportation occurs.If there is no direct transportation route between two nodes, it is expressed as   2  1 = 0. Since direct transportation cannot be realized between distribution centers and customers, it is obvious that when  1 ∈  ∧  2 ∈ ,   2 = 0 occurs.Vehicle   can be expressed as   = {ℎ  , {   }}, in which ℎ  = { ,  2 , ...,  +1 to finish the transportation task; {   } = { 1   ,  2   , ...,     } means that the goods of vehicle   contain the transportation demand of  customers.The time of vehicle   passing node   is expressed as      , which represents the total time of vehicle   that starts from  1 , passes  2 ,  3 , ...,  −1 , and reaches   ; the total time of vehicle   is expressed as    ; obviously,    =   +1   .The total goods weight of vehicle   is expressed as   = ∑   ∈{   }   .To meet any customer's demand   , 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   can be expressed as ℎ   = (   ,    ,  1  ,  2  ), which means that the goods path of   starts from distribution center    , transported by vehicle  1  to transfer station    , after being transferred by vehicle  2  to the customer.In the process,   must pass all edges of  1  before arriving at    and must pass all edges of  2

Table 1 :
Parameter values for generating instances.

Table 2 :
Results of the algorithms.