Presenting a Multi-Start Hybrid Heuristic for Solving the Problem of Two-Echelon Location-Routing Problem with Simultaneous Pickup and Delivery (2E-LRPSPD)

is study proposes a three-index ow-based mixed integer formulation to solve a two-echelon location routing problem with simultaneous pickup and delivery. In this formulation, pickup and delivery demands can be addressed using the same vehicle in each echelon of the network to reduce costs and increase logistics eciency. We solve such NP-hard problem by developing a multistart hybrid heuristic with path relinking (MHH-PR) which is composed of local search and a variable neighbourhood descent algorithm. In the algorithm, three constructive heuristics are applied to generate diversied initial solutions, and path relinking is introduced for intensication and postoptimisation. Results indicate that MHH-PR can reduce the gap between the near optimal and global optimal solutions by 1%-2%. e proposed algorithm signicantly improves computational eciency by reducing the computational time of more than 10 min for existing cases involving 20 nodes to less than 10 s.


Introduction
With the rapid growth of many e-retailers, including Amazon (US), Taobao (China), and Rakuten (Japan), many consumers have embraced online shopping in the past few years. e issues of picking, transporting, and delivering hundreds of millions of parcels consequently turn to be critical in the city logistics system. In 2018, the volume of express delivery of China reached 50.71 billion packages, according to China's State Post Bureau. Almost 140 million packages were handled in a single day. Express enterprises only have an e cient logistics network to deal with such a large number of express deliveries every day [1].
At present, the development of satellites around cities to establish a peripheral logistics mode has e ectively increased logistics e ciency and reduced costs because such logistics mode transports the goods dropped at these satellites instead of assigning vehicles from remote depots to directly conduct deliveries to end customers [2][3][4]. In China, many express enterprises have built some new distribution centres to improve the e ciency of the logistics network. For example, S.F. Express, one of the largest express delivery companies in China, built distribution centres in Inner Mongolia, Harbin, Yichun, Heihe, and Wuchang in 2018. A two-echelon location routing problem (2E-LRP) has been carefully investigated in designing a peripheral logistics mode. e 2E-LRP network consists of a main depot, many satellites, customers, primary vehicles, and secondary vehicles. e main depot is usually assigned with a large capacity and is located away from customers. e primary vehicles depart from the main depot, head to the satellites to meet the satellites' demand and return to the main depot. e satellites are typically equipped with a small capacity and are located close to customers. e secondary vehicles depart from the satellites, head to the customers to meet the customers' demand and return to the satellites (Figure 1). e 2E-LRP is solved by locating satellites, assigning customers and planning routes for primary and secondary vehicles.
e 2E-LRP investigated in the present study is di erent from the more general version of the 2E-LRP involving location decisions at depot and satellite levels. e location decision in this study only involves the satellite level because of two reasons. Firstly, most logistics enterprises in China have built their primary depots in major areas, and building new depots is unnecessary. Secondly, the cost of opening a depot is large, and locating a new depot is expensive for enterprises. e traditional 2E-LRP only covers the delivery aspect and cannot satisfy the increasing demand for simultaneous delivery and pickup of parcels. Consequently, the traditional 2E-LRP does not correspond with the current trends in logistics system development. JD mall, one of the top three largest 3C online retailers in China, has previously implemented a self-run logistics mode to provide delivery services to customers. Although such simple operational mode has achieved high e ciency, its corresponding logistics costs are large. On October 18, 2018, JD shi ed its logistics mode to provide delivery and pickup services simultaneously and thereby reduce the logistics costs while maximising existing logistics facilities. is change has provided convincing evidence to justify the importance of extending the traditional 2E-LRP to the 2E-LRP with simultaneous pickup and delivery (2E-LRPSPD). Under the 2E-LRPSPD, the need to simultaneously deliver and pickup parcels for customers is ful lled by the same vehicle. We present a three-index mixed integer linear programming for the 2E-LRPSPD. A multistart hybrid heuristic with path relinking (MHH-PR) is proposed to eciently solve the 2E-LRPSPD. Our model and heuristic can help enterprises, such as JD, to optimise their logistics network and can be applied to supermarkets, such as CR Vanguard, to plan their distribution network. Vanguard is one of the largest retail chain groups in China. It has many stores and suppliers in China, and it implements the 'backhaul' program by utilising these resources to reduce logistics costs. In this program, vehicles leave the distribution centre to complete delivery tasks and pick up goods from suppliers. Our model and heuristic can be used by Vanguard to appropriately select the location of its distribution centre and plan its distribution routes. e rest of this paper is organised as follows. Section 2 reviews previous relevant studies. Section 3 presents the de nition and mathematical formulation of the 2E-LEPSPD. Section 4 describes the proposed MHH-PR in detail. Section 5 reports the computational results. Section 6 presents a case study. Section 7 provides the conclusion and future directions.

Literature Review
is section is organised as follows. We provide an overview of the literature on the location routing problem (LRP) and then present a review on the LRPSPD, 2E-LRP, and 2E-LRPSPD. e literature related to the 2E-LRPSPD is also discussed.

Studies on LRP.
e LRP combines location and routing decisions, and the dependency between these two decisions was con rmed 50 years ago [5][6][7][8]. Watson-Gandy and Dohrn [9] were the rst to consider customer visits in locating depots, and many studies on LRP followed. A survey of the research on the standard LRP can be found in [10][11][12].

Studies on LRPSPD, 2E-LRP, and 2E-LRPSPD.
e LRPSPD is an extension of the LRP in terms of the types of customer demand. In the LRPSPD, the demand for delivery and pickup is integrated into the traditional LRP such that the pickup and delivery demands of customers are fulfilled at the same time. e delivery starts from the depot and ends at the depot a er pickup. Numerous approaches, such as branch-and-cut algorithm (B&C) [34], multistart SA [35], MA [36], variable neighbourhood scatter search algorithm [37], and hybrid heuristic approaches [38,39], have been proposed to solve the LRPSPD. e 2E-LRP is another extension of the LRP in terms of the network level. In 1980, Jacobsen and Madsen [40] officially presented the 2E-LRP, in which they assumed that each customer can be selected to serve as a satellite, and proposed three constructive heuristics to solve the model. Subsequently, various types of mixed integer programming models were developed by Crainic et al. [41], Vidović et al. [42], and Winkenbach et al. [43]. Cuda et al. [44] provided an overview of a two-echelon distribution system, including the 2E-LRP, two-echelon vehicle routing problem and truck and trailer routing problem.
Few exact solution methods have been reported because of the computational complexity of solving the 2E-LRP. e exact method for the 2E-LRP was firstly proposed by Contardo et al. [45]. ey embedded a two-index vehicle flow formulation and several families of valid inequalities into the B&C to obtain the lower bounds of small-and medium-sized instances of the 2E-LRP. As the 2E-LRP is an NP-hard problem, metaheuristic approaches and their hybrids have been well developed to solve the 2E-LRP and thereby improve computational efficiency, especially for large-scale networks. e frequently used heuristics can be well summarised as TS [46], adaptive large neighbourhood search [45], hybrid genetic algorithm (GA) [47,48], and hybrid heuristics [3,4,[49][50][51][52][53][54][55].
Various types of extensions of the conventional 2E-LRP have been investigated in the past decade to meet real-world requirements. Time constraints have been added to form the 2E-LRP with time window [56][57][58]. Ouhader and Elkyal [59] investigated a multisourcing and multiproduct 2E-LRP to design a pooled distribution supply chain, in which the route can end at different depots rather than the starting one of the classical 2E-LRP. Pichka et al. [60] addressed the two-echelon open location routing problem (2E-OLRP), in which the first-echelon and second-echelon vehicles do not return to the starting point. ree new mathematical models and a hybrid heuristic were proposed for the 2E-OLRP. Zhao et al. [61] constructed an optimisation model with consideration of practical constraints, such as vehicle capacity, working hours, and traffic restrictions, and presented a cooperative approximation heuristic algorithm for a two-echelon urban logistics network design under joint delivery alliances. Dai et al. [2] increased the dimensions of the LRP, thus obtaining one-echelon, two-echelon, three-echelon, and four-echelon LRPs.
Some scholars combined the LRPSPD and 2E-LRP to form the 2E-LRPSPD and thereby improve logistics efficiency. Few studies have focused on this problem. Rahmani et al. [62] firstly investigated the 2E-LRPPD. In their work, each processing centre (satellite) only provides one type of product, but the customer can ask for several types of products. e secondary vehicle must visit more than one processing centre to collect the different products and meet the customer's demand. e authors presented two LS methods to solve this problem. Demircan-Yildiz et al. [63] proposed two mixed integer programming formulations, namely, two-index nodebased and two-index flow-based formulations, and a family of valid inequalities to strengthen the formulations. However, they did not propose any algorithm to solve the model. Under such circumstance, Ghatreh and Hosseini-Motlagh [64] proposed a fuzzy programming approach and a combined heuristic method (GA and SA) to solve the 2E-LRPSPD under fuzzy demand. e low computational efficiency did not satisfy the basic requirement of real-time scheduling because solving a small-sized network (20 customers) takes 14 min on average, whereas a medium-sized network (50 customers) takes 20 min. Table 1 summarises some the literature described previously.

Studies Related to 2E-LRPSPD.
e 2E-LRPSPD deals with the location of satellites. For a facility location problem, Bruno et al. [65] presented a mathematical model for facility location systems in noncompetitive contexts. An and Svensson [66] reviewed the approximation algorithms for a facility location problem. Murray and Church [67] designed a SA for p-median and maximal covering location problems. Masone et al. [68] developed a three-stage solution method for the optimal diversity management problem. Additional details can be found in [69][70][71][72]. e pickup and delivery problem is another research direction related to the 2E-LRPSPD. e classification, models and methods of the PDP can be found in [73,74]. e contributions of the present study to the improvement of computational efficiency in solving the 2E-LRPSPD are summarised as follows: (1) development of a three-index mixed integer linear programming; and (2) development of a hybrid heuristic that combines VND and LS and is embedded with the path relinking (PR) strategy proposed by Glover and Journal of Advanced Journal of Advanced Transportation 6 the formulation. Table 2 provides the parameters and variables in our mathematical formulation.
To better highlight the contribution of the proposed heuristic, the three-index ow-based mixed integer formulation is described in Appendix 1. Figure 3 shows a owchart with input, output and constraints.

Multistart Hybrid Heuristic with Path
Relinking e proposed optimisation model is NP-hard. Exact methods could be applied to small-scale networks to nd the optimal solution. However, such methods for large-scale networks would either fail to nd a global optimal solution or require a long computation time, which is unsuitable for real-time eet dispatching tools. us, the development of heuristics is the key issue in yielding appropriately good solutions to the model in a short time period. Population-based approaches (such as GA, ant colony algorithm, and PSO) are not applicable for this problem because of the large runtime and memory storage required by genetic operators. e structure of the solution is complicated because the capacity constraints of satellites and vehicles are considered.
is problem cannot be solved e ciently using population-based approaches. e iterative search of trajectory-based heuristics (such as SA and TS) is serial with a single state movement rather than parallel search. Diversity is insu cient when concentration and diversity are equally important. A local optimal solution can be easily obtained because the e ectiveness of trajectory-based heuristics depends on the quality of the initial solution. us, this study proposes the MHH-PR using three constructive heuristics ( 1, 2, and 3), VND algorithm, two LS ( 1 and 2) mechanisms, and a PR strategy to solve the 2E-LRPSPD. e result of heuristics based on search procedure depends on the quality of the initial solution. ree di erent constructive heuristics based on the greedy randomisation principle are used to generate the initial solutions to increase the diversity of the solutions. e three di erent heuristics can generate satisfactory initial solutions for di erent distributions of customers. VND and LS are adopted in the proposed hybrid to achieve a balance between the development of the current optimal solution and the exploration of the unknown search space. VND is widely used to solve large discrete optimisation problems because it requires a small number of parameters and short runtime. VND expands the search space while browsing di erent neighbourhoods, thus causing the current solution to continuously approach the optimal solution. LS is a simple heuristic method that searches for the neighbours of the current solution until an improved solution is found. Using PR in the proposed heuristic can help to explore the trajectory between pairs of elite solutions and ultimately obtain a competitive solution.
e main procedure of the MHH-PR executes maxiter iterations. In each iteration, the MHH-PR starts with an initial solution generated by 1, 2, and 3, and the best solution obtained by 2 with three neighbourhood structures is improved by VND. In the VND phase, Relocate (N1), Laguna [75] to investigate the trajectory between pairs of elite solutions by gradually transforming one elite solution into another to obtain a satisfactory solution. e developed algorithm can enrich related research and accelerate the convergence process.

Problem Statement
is section describes the problem de nition and the three-index ow-based mathematical formulation of the 2E-LEPSPD.

Problem De nition.
Logistics and distribution enterprises construct their regional two-echelon distribution networks by developing satellites between depot and customers. A twoechelon distribution network is composed of a main depot, a number of candidate satellites, and a number of customers. e candidates of geographic locations of main depot, candidate satellites, and customers are given. e candidate satellites have di erent capacities and opening costs. e main depot and satellites form the rst-level or primary network, and satellites and customers form the second-level or secondary network. In the network, customers have pickup and delivery demands, and these demands are met simultaneously by one vehicle. Two types of vehicles, namely, primary and secondary vehicles, are used in this network. In the secondary network, the secondary vehicle with a small capacity departs from a satellite and returns to the same satellite a er completing the pickup and delivery tasks for all customers during the secondlevel trip. Similarly, in the primary network, the primary vehicle with a large capacity departs from the main depot and returns to the main depot a er completing the pickup and delivery of all opened satellites during the rst-level trip. e 2E-LRPSPD is investigated in this work to determine the locations of satellites by selecting the candidates and trips of vehicles to minimise the total cost in the two-echelon system under the following assumptions.
(i) e total load on primary and secondary vehicles at any point of a trip does not exceed the vehicle capacity.
(ii) e total pickup and delivery demands of customers assigned to a satellite do not exceed the capacity of satellites. (iii) Each customer (satellite) in the primary (secondary) network is served by exactly one primary (secondary) vehicle. Figure 2 depicts a two-echelon network composed of one depot (node 0), ve candidate satellites (indexed from 1 to 5), and 15 customers (indexed from 6 to 20). e capacity of each satellite in the network is 25, and the capacities of the primary and secondary vehicles are 40 and 15, respectively. e delivery and pickup demands of customers are provided in brackets. e rst item in the bracket is the delivery demand, and the second item is the pickup demand.

Mathematical Formulation.
e 2E-LRPSPD can be de ned using a complete directed graph = ( , ).
e parameters and variables are introduced before presenting starts by selecting the nearest customer to the satellite in . en, the nearest customer to the current customer is determined. e customer is appended to the trip when the vehicle capacity is not violated. Otherwise, the vehicle returns to , and another trip is constructed for considering that the nearest customer is not visited in . A er all customers are visited in , the same procedure is used to construct the trips for other opened satellites until all customers are covered. NHH is also applied to the rst-level network to build the rst-level trips. (2) 2. Heuristic 2 by Nguyen et al. [3,4] is used as 2 in this study. is heuristic is an extended NNH. 2 is di erent from 1. 1 rstly assigns all customers to the nearest satellite and uses NNH to construct the trips for each satellite. 2 randomly selects a satellite and constructs a second-level trip for the opened satellite by adding the nearest customer that is not visited. Its delivery and pickup demands meet the capacity of the vehicle until no additional customers can be added. In this way, trips are continuously constructed until the satellite has no remaining capacity to satisfy the unvisited customers, and a new satellite is randomly opened. e heuristic terminates when all customers are inserted into the trips. e NHH described in 1 is applied to the rst-level network to build the rst-level trips. 3.
3 is similar to heuristic 3 of Nguyen et al. [3,4]. is heuristic is an insertion method. e heuristic by Mole and Jameson [76] is adopted to Exchange (N2), Or-Opt (N3), and 2-Opt (N4) are used to optimise the second-level trips of the solution. en, LS1 is used to improve the rst-level trips of the improved solution. We call PR to explore the trajectory of the current and elite solutions in the solution space and attempt to obtain a satisfactory solution, which we add to the pool of elite solutions. Constructive heuristics 2, VND, 1, and PR are continuously used until the stopping criterion is met. A er conducting all the iterations, we call PR again to explore the trajectory between the two elite solutions in the pool. Figure 4 depicts the owchart of the MHH-PR in detail.

Initial Solutions.
A feasible solution consists of a set of opened satellites, a set of second-level trips and a set of rst-level trips. is study uses three greedy randomised constructive heuristics called 1, 2, and 3 to obtain the initial solution for each iteration and achieve enhanced diversity. e constructive heuristics are described in detail as follows: (  (ii) N2-Exchange: ( Figure 6) Exchange two customers from the same trip or from di erent trips which pertain to the same or di erent satellites. All Exchange moves can be explored in 2 .
(iii) N3-Or-Opt: (Figure 7) is move mainly reinserts a string with a length of 2-3 to the same trip or a di erent trip. Regardless of whether the inserted trip is the same or not, two cases depend on whether the string is inverted. All N3 moves can be tested in 2 .
construct the second-level trips. Open satellites are linked using the same method as 1. A detailed description of this heuristic can be found in [3,4].

Variable Neighbourhood Descent.
VND is used to inspect neighbourhood structures , = 1, 2, ⋅ ⋅ ⋅ , max . If an improving move can be found in neighbourhood structure , then this move is executed, and is reset to 1; otherwise, the next neighbourhood structure +1 continues to be searched. VND stops when all neighbourhood structures are browsed without improvement.
e VND used in this study contains four neighbourhood structures ( max = 4): N1, N2, N3, and N4. It is used to improve the second-level trips of the solution. Binary variable with a value of 1 when primary vehicle ∈ travels directly from node to node , , ∈ 0 ∪ Binary variable with a value of 1 when secondary vehicle ∈ travels directly from node to node , , ∈ ∪ Binary variable with a value of 1 when customer ∈ is allocated to satellite ∈ Status of satellite with a value of 1 when satellite ∈ is open best solution among the 6 ᐈ ᐈ ᐈ ᐈ ᐈ ᐈ ᐈ ᐈ solutions is then selected. e best solution is taken as the new current solution when it is better than the current solution, and the neighbours of the new current solution are continuously searched. 2 is stopped when the current solution does not improve in 10 successive steps.
(iv) N4-2-Opt: (Figure 8) Remove two noncontiguous edges from one or two trips, and add the two other edges to reconnect the resulting string. ree cases are used as the two edges may be from the same trip or di erent trips with the same or distinct satellites.
(a) Two edges are from the same trip. is situation is known as the 2-Opt move for the TSP. (b) Two edges are from di erent trips with the same satellite. Two di erent cases depend on whether the string is inverted. (c) Two edges are from two trips with two di erent satellites. Two di erent cases depend on whether the string is inverted.
All N4 moves can be tested in 2 .   combinations of di erent opened satellites. All N7 moves can be explored in 2 .

Path
Relinking. PR was proposed by Glover and Laguna [75] to improve the solutions obtained via TS. is method can be added to any metaheuristic as an intensi cation strategy. e principle of this method is to generate a pool of elite solutions and explore the trajectory between pairs of elite solutions in the pool to obtain the best solution. Nguyen et al. [3,4] applied this method to solve the 2E-LRP. ey believed that such method is easier to apply to big tours compared with a 2E-LRP solution because of the capacity constraints of satellites and vehicles. A big tour is obtained by concatenating the second-level trips in the order of the rst-level trips.
is study uses the PR method designed in [3,4]. Nguyen et al. [3,4] investigated the 2E-LRP without considering the possibility of a customer having pickup and delivery demands simultaneously. erefore, we design a repair procedure and improved procedures that are used to improve S according to the characteristics of our problem. Algorithm 2 summarises the PR for two big tours and .
(i) N5-Closing an opened satellite: A satellite is randomly selected from two satellites with the least and second least quantity of allocated customers, and its customers are allocated to other opened satellites with su cient residual capacity. All N5 moves can be evaluated in ( ). (ii) N6-Changing all customers of two opened satellites: Two opened satellites are randomly selected, and their customers are exchanged when their capacities can meet the demands of all customers of the other satellite. All N6 moves can be tested in 2 .
(iii) N7-Closing an opened satellite and opening a closed one: e more customers to each opened satellite are allocated, the more likely the opened satellite will close. An opened satellite is selected on the basis of the probability that it is proportional to the number of customers assigned to it. e opened satellite customers are transferred to a randomly selected closed satellite a er checking the satellite capacity constraint. en, the selected closed satellite is opened, and the original opened satellite is turned o . e neighbourhood structure is to search for   In Algorithm 2, the balance procedure (lines 1-10) is called when and di er in length. Let ( , ) denote the frequency of satellite in big tour . e balance procedure is to add copies of satellites at the end of a big tour to ensure denote the rank of a customer in a big tour . PR rstly looks for customers with di erent locations in and (line 13). If a customer is found, then it is swapped with another point (lines [14][15]. Otherwise, a misplaced satellite is identi ed, and points are exchanged (lines [16][17][18][19][20][21][22][23][24][25][26][27][28][29]. e obtained big tour may not satisfy the vehicle or satellite capacity constraints. en, the repair procedure is called to repair the big tour on the basis of a certain probability to limit the runtime (lines 30-33). Figure 9 shows a small example of the repair procedure. ree candidate satellites (indexed from 1 to 3) and ten customers (indexed from 4 to 13) are identi ed in this example.

Computational Results
Considering that the 2E-LRPSPD is a new problem, we cannot use existing results to compare the proposed heuristic's performance with those of other approaches. us, the proposed heuristic is tested on the 2E-LRP and LRPSPD instances by changing the heuristic.

18
seek an index such that ̸ = ;

19
if satellite = occurs in then 20 seek index such that = ̸ = ; comparing the results with the results of Nguyen et al. [3,4], Contardo et al. [45], and Pichka et al. [60]. Table 3 reports the results. In this table, the rst column is the name of the 2E-LRP instances. e next two columns report the best results in [3,4]. Column " " shows the best objective values found by the heuristic. Column "CPU(s)" shows the average CPU runtime. Columns 4-5, 6-7, and 8-9 are the same as columns 2-3 and show the results in [45,60] and this study. e text in bold characters indicates the found in the four papers for the corresponding instance. e last column, "Gap (%)", shows the di erence between the best solution of our heuristic and . = MHH_PR − / × 100. e average of all instances for each column is given in the last row of the table. As shown in Table 3, the proposed heuristic obtains 5 optimal solutions out of 30; this result is less than the eight solutions by Nguyen et al. [3,4] and the seven solutions by Pichka et al. [60]. We can nd the optimal solution for all small-sized instances with 20 customers and a medium-sized instance with 50 customers. Five instances have Gap rates larger than 3%. e overall 1.82% gap shows the good performance of the proposed heuristic. Regarding CPU runtime, we can observe that the CPU time consumed by the proposed heuristic to solve the instances with = 70 is longer than the instances with = 150. is condition is because the di culty of solving the problem increases when the secondary vehicle has a small capacity and many vehicles are needed to complete the delivery tasks of all customers. e average CPU time of our results is greater than those in other studies; however, this is a reasonable time to solve our problem involving strategic location decisions. e locations in the few years cannot be changed when the locations of opened satellites are determined because of the large cost of change. us, few minutes of computer time should be consumed to nd the best locations of satellites.
compares the performances of MHH and MHH-PR in the 2E-LRPSPD instances. MHH-PR is coded on MATLAB R2017a, and all experiments are performed on a computer with Intel Xeon 2.80 GHz and 8 GB RAM. A er some preliminary experiments, the following parameters are implemented in the proposed heuristic to solve the 2E-LRP, LRPSPD, and 2E-LRPSPD test instances: maxiter is set to 20; and are taken as 0.15 and 0.10, respectively; is taken as 0.002; is set to 1/3.

Performance of MHH-PR on 2E-LRP Instances.
e pickup demand of the customer is zero, and the 2E-LRPSPD can be considered a 2E-LRP when the customer's original demand is regarded as the delivery demand. us, the validity of the heuristic can be veri ed by solving the 2E-LRP instances. We adopt Prodhon's 2E-LRP instances from [3,4] Table 4. e results are compared with the optimal upper bound obtained by Karagoglan et al. [34]. In this table, the rst column represents the names of the LRPSPD instances.

Performance of MHH-PR in LRPSPD.
e 2E-LRPSPD can be regarded as LRPSPD when the costs of rst-level edges are zero. us, the proposed heuristic can be utilised to solve LRPSPD instances.
We provide brief information about the LRPSPD test instances proposed by Barreto. Barreto's LRP instances include 15 test instances obtained from other LRP literature or by adding depots to classic VRP instances (these instances can be found in http://sweet.us.pt/~iscf143). ese instances contain the following characteristics: the number of customers is 8-100, and the number of satellites varies between 2 and 15. e name of each instance is represented by the coding structure of the author, Author, who generated the instance; year of the instance, Year; the number of customers, ᐈ ᐈ ᐈ ᐈ ᐈ ᐈ ᐈ ᐈ ; and the problem when the number of customers is greater than 50 within 4 h. MHH-PR does not exceed 160 s in solving the instances with less than 100 customers. MHH-PR obtains these results by consuming 37.03 s on average for 30 instances , making MHH-PR more e cient. MHH-PR can nd optimal/near optimal solutions in a short time period. ese results show that MHH-PR is a preferable approach for solving LRPSPD instances.

Performance of MHH-PR in 2E-LRPSPD Instances.
Here, we compare the MHH and MHH-PR in solving the 2E-LRPSPD. We utilise 20 Prodhon's 2E-LRP instances and the SN demand separation approach to construct the instances for the 2E-LRPSPD. and represent two di erent demand separation strategies. Tables 5 and 6 present the results of the MHH and MHH-PR in solving the -type and -type 2E-LRPSPD instances. In Tables 5 and 6, the rst column denotes the names of instances and is consistent with the format in Table 3. e MHH and MHH-PR are run for ve times in Column "DSS" shows the demand separation strategy. Column " " shows the optimal upper bound obtained in [34]. e next two columns represent the Gap value and CPU of B&C in [34]. e results of the proposed heuristic are shown in the last three columns. Column " " shows the best objective value found in ve runs of the proposed heuristic. Column "Gap" shows the gap between the results of the proposed heuristic and using the following formula: Gap = (( − )/ ) × 100. Column "CPU" shows the average CPU runtime for ve runs. e last row of the table gives the average of the columns.
As shown in Table 4, the average Gap of MHH-PR in this study is 1.00%, which is less than 2.20% in [34], and the optimal solutions for all small-sized instances except Gaskell67_36 * 5 can be found via MHH-PR. Four instances have Gap larger than 3%. For CPU time with small-sized instances less than 36 customers, the proposed heuristic and B&C can nd the optimal solutions, but the proposed heuristic consumes less than CPU time than that of B&C in most instances. B&C cannot solve the in solving -type instances (1.67%). e combined PR strategy in the MHH can nd a competitive solution on the basis of the value of Gap. e improvement rates of PR for the -type and -type instances are 0.69% and 0.73%, respectively. PR considers the quality and dispersion of the solution in the search process and thus helps the MHH to obtain a good solution. As a result, the search performance of the entire heuristic is improved. Regarding the runtime for -type instances, the average CPU time of the MHH is 315.95 s, whereas that the MHH-PR is 463.65 s. e same phenomenon can be observed from the average CPU times of the MHH and MHH-PR in the -type instances. Although the runtime increases, a good solution is obtained by the cost of consumption at a certain time.

Case Study
In this section, we use the proposed heuristic to solve a case in the work of Belgin et al. [79]. is case is the distribution operations on 25 supermarkets (customers) of a supermarket each instance. e next three columns show the results of the -type (or -type) 2E-LRPSPD instances solved by the MHH. Column " " shows the objective value of the best feasible solution in ve runs. Column " " shows the robustness index of the heuristic which is calculated as where " " is the worst solution of ve runs. Column "CPU" shows the average CPU time used in ve runs of the heuristic. e next three columns give the results of the instances solved by the MHH-PR. e last column indicates the gap of the best solutions obtained by the MHH and MHH-PR. Tables 5 and 6, the average obtained by the MHH in solving all -type instances is 1.88%, which is less than that obtained in solving all -type instances (2.09%). is condition shows that the MHH is relatively stable in solving -type instances. e same conclusion can be derived from the MHH-PR. e average obtained by the MHH-PR in solving the -type instances is 2.08%, which is larger than that obtained distance of vehicles in the two-echelon distribution system is shorter than that in the first-echelon system. e opened satellite of our optimal results is S1, and the total travelling distance of the new two-echelon system is 1510.85 km. e total travelling distance is reduced by approximately 57.03% relative to that of the single-echelon distribution system (3516.30 km).

As shown in
is result shows that the new two-echelon system is superior to the single-echelon system in terms of travelling distance. e same conclusion can be obtained in the results of [79]. Assume that the travelling cost of a unit distance and the fixed cost of small vehicles are 1 and 50, respectively; and that the travelling cost of a unit distance and the fixed cost of large vehicles are 2 and 100, respectively. e total costs (including fixed and travelling costs) of the single-echelon and two-echelon systems are 3766.3 and 3296.85, respectively, and the total cost of the two-echelon system was approximately 12.46% less than that of the single-echelon system. Although the establishment of satellites requires a certain cost, the benefits of reduced vehicle travelling costs are large. capacity are used in the current distribution system because large vehicles are not allowed to traverse urban areas. e managers of the supermarket chain decide to locate one or two satellites at the city boundary in establishing a new two-echelon distribution system. e satellites have the same capacities and meet the delivery and pickup demands of all customers. In the new system, large vehicles with a 26 ton capacity are used in the primary distribution network to complete the distribution operations between the main depot and the satellites. Small vehicles with 11 ton capacity are used in the secondary network to complete the distribution operations between the satellites and the customers. Belgin et al. [79] provided the distances between two candidate satellites (S1, S2), customers (C1-C25) and average delivery and pickup demands of customers, as shown in Tables 8 and 9 in Appendix 2.
In addition, our results are better than the best results in [79] in terms of travelling distance. We summarize the following managerial insights: (1) Using large-capacity vehicles in the primary network for long distances can take advantage of the economies of scale by large vehicles and reduce the costs while increasing e ciency. Although two large-capacity vehicles are used, the travelling distance and cost of vehicles are lower than the vehicle xed cost, and (2) Companies should build few satellites with large capacities. e increase in the number of satellites increases the input of satellite opening and vehicle travelling costs and the di culty of controlling and managing the distribution operations of a company.

Conclusion
In this study, we investigate the 2E-LRPSPD, which is an extension of the 2E-LRP. e problem includes delivery and pickup activities at each echelon with satellites and vehicles of limited capacities. A three-index ow-based mixed integer programming formulation is proposed to deal with such problem. As the proposed problem is an NP-hard problem, a global optimal solution cannot be obtained at an acceptable runtime, especially for large-sized networks. us, the MHH-PR with three greedy randomised construction heuristics is presented to increase the diversity of solutions. In particular, two LS algorithms, a VND and a PR are adopted to enhance the algorithm's search ability. We evaluate the performance of the MHH-PR in 2E-LRP and LRPSPD instances. e computational results indicate that the MHH-PR can obtain good quality solutions (the average gaps are 1.82% and 1.00% for the 2E-LRP and LRPSPD instances, respectively) and that the computational time is signi cantly reduced (the average CPU times are 312.3 and 37.03 s for the 2E-LRP and LRPSPD instances, respectively). We compare the performance of the MHH and MHH-PR in the 2E-LRPSPD instances. e results reveal that the proposed algorithm can nd the best and most competitive solution. e implementation of the MHH-PR in a distribution system from the literature shows that relative to existing single-echelon distribution systems, the two-echelon distribution system obtained by the proposed heuristic can reduce the total travelling distance by 57.03% and the costs while increasing the e ciency.
Several directions are provided for future research on this topic. e problem can be studied by considering many realistic constraints, such as time windows and vehicle synchronisation. Exact algorithms can be designed for the 2E-LRPSPD to evaluate the proposed heuristic.     Tables   Tables 8-10.

B. The Related
Data Availability e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.