Abstract
This study proposes a three-index flow-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 efficiency. 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 diversified initial solutions, and path relinking is introduced for intensification and postoptimisation. Results indicate that MHH-PR can reduce the gap between the near optimal and global optimal solutions by 1%-2%. The proposed algorithm significantly improves computational efficiency by reducing the computational time of more than 10 min for existing cases involving 20 nodes to less than 10 s.
1. 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. The 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 efficient 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 effectively increased logistics efficiency 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–4]. In China, many express enterprises have built some new distribution centres to improve the efficiency 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. The 2E-LRP network consists of a main depot, many satellites, customers, primary vehicles, and secondary vehicles. The main depot is usually assigned with a large capacity and is located away from customers. The primary vehicles depart from the main depot, head to the satellites to meet the satellites’ demand and return to the main depot. The satellites are typically equipped with a small capacity and are located close to customers. The secondary vehicles depart from the satellites, head to the customers to meet the customers’ demand and return to the satellites (Figure 1). The 2E-LRP is solved by locating satellites, assigning customers and planning routes for primary and secondary vehicles.
The 2E-LRP investigated in the present study is different from the more general version of the 2E-LRP involving location decisions at depot and satellite levels. The 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.
The 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 efficiency, its corresponding logistics costs are large. On October 18, 2018, JD shifted its logistics mode to provide delivery and pickup services simultaneously and thereby reduce the logistics costs while maximising existing logistics facilities. This 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 fulfilled 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 efficiently 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.
The rest of this paper is organised as follows. Section 2 reviews previous relevant studies. Section 3 presents the definition 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.
2. Literature Review
This 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. The literature related to the 2E-LRPSPD is also discussed.
2.1. Studies on LRP
The LRP combines location and routing decisions, and the dependency between these two decisions was confirmed 50 years ago [5–8]. Watson–Gandy and Dohrn [9] were the first 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–12].
In the last two decades, several LRP variants have been investigated, and they include the following: (1) stochastic and fuzzy LRPs; (2) dynamic, periodic, and multiperiod LRPs; (3) multiechelon LRPs (mainly the 2E-LRP); (4) location arc-routing problems; (5) generalized LRPs; (6) prize-collecting LRPs; (7) split delivery LRPs; (8) LRP with simultaneous delivery and pickup (LRPSPD); and (9) inventory LRPs. The variants of the LRP can be found in the review of Drexl and Schneider [13].
Previous studies proposed several exact algorithms for small- and medium-sized LRPs and their variants [14–18]. As the LRP and its variants are NP-hard problems, the heuristics for solving them have been established, and they include TS [19, 20], evolutionary local search (LS) [21], variable neighbourhood descent (VND) [22], simulated annealing heuristic (SA) [23–25], particle swarm optimisation (PSO) [26], variable neighbourhood search algorithm (VNS) [27], very large-scale neighbourhood search [28], Lagrangian relaxation granular tabu search (TS) heuristic [29], memetic algorithm (MA) with population management [30, 31], greedy randomized adaptive search procedure [32], and multiple ant colony optimisation algorithm [33].
2.2. Studies on LRPSPD, 2E-LRP, and 2E-LRPSPD
The 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. The delivery starts from the depot and ends at the depot after 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.
The 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. The exact method for the 2E-LRP was firstly proposed by Contardo et al. [45]. They 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. The 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–55].
In the following, we classify the above literature on the 2E-LRP on the basis of the notation introduced in [11] and [44] to consider location decisions. The location decisions of [45, 49–55] only involved one level, namely, the satellite level. In [3, 4, 46–48], the location decisions at two stages, namely, satellite and depot levels, were considered.
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–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. Three 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. The secondary vehicle must visit more than one processing centre to collect the different products and meet the customer’s demand. The authors presented two LS methods to solve this problem. Demircan-Yildiz et al. [63] proposed two mixed integer programming formulations, namely, two-index node-based 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. The 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.
2.3. Studies Related to 2E-LRPSPD
The 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–72].
The pickup and delivery problem is another research direction related to the 2E-LRPSPD. The classification, models and methods of the PDP can be found in [73, 74].
The 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 Laguna [75] to investigate the trajectory between pairs of elite solutions by gradually transforming one elite solution into another to obtain a satisfactory solution. The developed algorithm can enrich related research and accelerate the convergence process.
3. Problem Statement
This section describes the problem definition and the three-index flow-based mathematical formulation of the 2E-LEPSPD.
3.1. Problem Definition
Logistics and distribution enterprises construct their regional two-echelon distribution networks by developing satellites between depot and customers. A two-echelon distribution network is composed of a main depot, a number of candidate satellites, and a number of customers. The candidates of geographic locations of main depot, candidate satellites, and customers are given. The candidate satellites have different capacities and opening costs. The main depot and satellites form the first-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 after completing the pickup and delivery tasks for all customers during the second-level trip. Similarly, in the primary network, the primary vehicle with a large capacity departs from the main depot and returns to the main depot after completing the pickup and delivery of all opened satellites during the first-level trip.
The 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) The total load on primary and secondary vehicles at any point of a trip does not exceed the vehicle capacity.(ii) The 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), five candidate satellites (indexed from 1 to 5), and 15 customers (indexed from 6 to 20). The capacity of each satellite in the network is 25, and the capacities of the primary and secondary vehicles are 40 and 15, respectively. The delivery and pickup demands of customers are provided in brackets. The first item in the bracket is the delivery demand, and the second item is the pickup demand.
3.2. Mathematical Formulation
The 2E-LRPSPD can be defined using a complete directed graph . The parameters and variables are introduced before presenting 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 flow-based mixed integer formulation is described in Appendix 1. Figure 3 shows a flowchart with input, output and constraints.
4. Multistart Hybrid Heuristic with Path Relinking
The proposed optimisation model is NP-hard. Exact methods could be applied to small-scale networks to find the optimal solution. However, such methods for large-scale networks would either fail to find a global optimal solution or require a long computation time, which is unsuitable for real-time fleet dispatching tools. Thus, 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. The structure of the solution is complicated because the capacity constraints of satellites and vehicles are considered. This problem cannot be solved efficiently using population-based approaches. The iterative search of trajectory-based heuristics (such as SA and TS) is serial with a single state movement rather than parallel search. Diversity is insufficient when concentration and diversity are equally important. A local optimal solution can be easily obtained because the effectiveness of trajectory-based heuristics depends on the quality of the initial solution. Thus, this study proposes the MHH-PR using three constructive heuristics (, , and ), VND algorithm, two LS ( and ) mechanisms, and a PR strategy to solve the 2E-LRPSPD.
The result of heuristics based on search procedure depends on the quality of the initial solution. Three different constructive heuristics based on the greedy randomisation principle are used to generate the initial solutions to increase the diversity of the solutions. The three different heuristics can generate satisfactory initial solutions for different 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 different 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.
The main procedure of the MHH-PR executes maxiter iterations. In each iteration, the MHH-PR starts with an initial solution generated by , , and , and the best solution obtained by with three neighbourhood structures is improved by VND. In the VND phase, Relocate (N1), Exchange (N2), Or-Opt (N3), and 2-Opt (N4) are used to optimise the second-level trips of the solution. Then, LS1 is used to improve the first-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 VND, , and PR are continuously used until the stopping criterion is met. After conducting all the iterations, we call PR again to explore the trajectory between the two elite solutions in the pool. Figure 4 depicts the flowchart of the MHH-PR in detail.
4.1. Initial Solutions
A feasible solution consists of a set of opened satellites, a set of second-level trips and a set of first-level trips. This study uses three greedy randomised constructive heuristics called , , and to obtain the initial solution for each iteration and achieve enhanced diversity. The constructive heuristics are described in detail as follows:(1).
firstly calculates the absolute value of the difference between delivery demand and pickup demand for each customer , . Each customer is allocated to a satellite that is randomly selected from the first two nearest satellites with sufficient residual capacity in decreasing order of . Let denote the set of customers allocated to satellite . After arranging all the customers, a nearest neighbour heuristic (NNH) is used to construct the trips of the arranged customers of each satellite. For an opened satellite , NNH starts by selecting the nearest customer to the satellite in . Then, the nearest customer to the current customer is determined. The 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 . After 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 first-level network to build the first-level trips.
(2).
Heuristic by Nguyen et al. [3, 4] is used as in this study. This heuristic is an extended NNH. is different from . firstly assigns all customers to the nearest satellite and uses NNH to construct the trips for each satellite. 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. The heuristic terminates when all customers are inserted into the trips. The NHH described in is applied to the first-level network to build the first-level trips.
(3).
is similar to heuristic of Nguyen et al. [3, 4]. This heuristic is an insertion method. The heuristic by Mole and Jameson [76] is adopted to construct the second-level trips. Open satellites are linked using the same method as . A detailed description of this heuristic can be found in [3, 4].
4.2. Variable Neighbourhood Descent
VND is used to inspect neighbourhood structures . 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 continues to be searched. VND stops when all neighbourhood structures are browsed without improvement.
The VND used in this study contains four neighbourhood structures (): N1, N2, N3, and N4. It is used to improve the second-level trips of the solution.(i) N1—Relocate: (Figure 5) Remove the customer from its trip and reinsert it into a new position in the same trip or in another trip that is connected to the same satellite or another satellite. N1 can be evaluated in .(ii) N2—Exchange: (Figure 6) Exchange two customers from the same trip or from different trips which pertain to the same or different satellites. All Exchange moves can be explored in .(iii) N3—Or-Opt: (Figure 7) This move mainly reinserts a string with a length of 2-3 to the same trip or a different 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 .(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. Three cases are used as the two edges may be from the same trip or different trips with the same or distinct satellites.(a) Two edges are from the same trip. This situation is known as the 2-Opt move for the TSP.(b) Two edges are from different trips with the same satellite. Two different cases depend on whether the string is inverted.(c) Two edges are from two trips with two different satellites. Two different cases depend on whether the string is inverted.
All N4 moves can be tested in .
4.3. Local Search
This study uses two LS procedures, namely, and , where improves the first-level trips of the solution that mainly contains two neighbourhood structures, namely, N1 and N2. includes three neighbourhood structures, namely, closing an opened satellite (N5), changing all customers of two opened satellites (N6), and closing an opened satellite and opening a closed one (N7). The steps of are described in Algorithm 1. The phase randomly generates feasible solutions by using each neighbourhood structure of N5, N6, and N7. The best solution among the solutions is then selected. The 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. is stopped when the current solution does not improve in 10 successive steps.
| ||||||||||||||||||||||||
(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 sufficient 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 .(iii) N7—Closing an opened satellite and opening a closed one: The 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. The opened satellite customers are transferred to a randomly selected closed satellite after checking the satellite capacity constraint. Then, the selected closed satellite is opened, and the original opened satellite is turned off. The neighbourhood structure is to search for combinations of different opened satellites. All N7 moves can be explored in .
4.4. Path Relinking
PR was proposed by Glover and Laguna [75] to improve the solutions obtained via TS. This method can be added to any metaheuristic as an intensification strategy. The 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. They 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 first-level trips. This 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. Therefore, 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 .
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
In Algorithm 2, the balance procedure (lines 1–10) is called when and differ in length. Let denote the frequency of satellite in big tour . The balance procedure is to add copies of satellites at the end of a big tour to ensure for each satellite . Line 11 generates an intermediate big tour when and have the same length. Let denote the Hamming distance of and to guide the path from to .
Each iteration of the while loop (lines 12) generates one intermediate big tour to reduce . Let denote the rank of a customer in a big tour . PR firstly looks for customers with different 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 identified, and points are exchanged (lines 16–29). The obtained big tour may not satisfy the vehicle or satellite capacity constraints. Then, 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. Three candidate satellites (indexed from 1 to 3) and ten customers (indexed from 4 to 13) are identified in this example. . . Big tour (1-8-6-10-7-13-3-9-11-4-1-12-5) is converted into the corresponding second-level trips. On the basis of the dynamic load of the vehicle marked on the horizontal line, the dynamic load of the vehicle departing from satellite 1 exceeds its capacity constraint (Equations (A.18), (A.19), (A.21), and (A.23)). The pickup demand of satellite 1 exceeds the satellite capacity (Equations (A.13–A.15)). We remove the customer from the end of the trip until the dynamic load of the vehicle meets the capacity constraint. The set SC contains all the removed customers, . Each customer in the SC is inserted into an existing trip that meets the vehicle and satellite capacity constraints via insertion [77]. A new trip is constructed for the customer when such a trip is not found. The starting point of the new trip may be an opened satellite with sufficient residual capacity or a closed satellite. In this example, all trips have insufficient residual capacities for customer 7 or 13, and a new trip is constructed for customers 7 and 13. The feasible big tour and its second-level trips are obtained. A feasible big tour is converted into a 2E-LRP solution by adding first-level trips (line 35). is improved using and updated when (lines 36-37). is improved using VND and updated when (lines 38-39). The best solution (lines 42–44) is updated.
5. 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. Thus, the proposed heuristic is tested on the 2E-LRP and LRPSPD instances by changing the heuristic. This section presents the results of the computational experiments in three stages. The first stage evaluates the performance of the MHH-PR in the 2E-LRP instances. The second stage investigates the performance of the MHH-PR in the LRPSPD instances. The last stage 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. After 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.
5.1. Performance of MHH-PR on 2E-LRP Instances
The 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. Thus, the validity of the heuristic can be verified by solving the 2E-LRP instances. We adopt Prodhon’s 2E-LRP instances from [3, 4] to investigate the performance of the MHH-PR in solving the 2E-LRP. These instances are derived from Prodhon’s LRP instances by adding a main depot with coordinates (0, 0) using original depots as satellites and assuming that the distances from the main depot to the satellites are twice the Euclidean distance. This dataset contains 30 2E-LRP instances with capable satellites and vehicles. These instances have the following characteristics: number of customers , number of satellites , capacity of secondary vehicle and number of clusters . The capacity of the primary vehicle is 1.5 times the maximum satellite capacity. The cost of the second-level edge is the Euclidean distances of the edge, multiplied by 100 and rounded to the nearest integer. The cost of the first-level edge is calculated by multiplying the distance of the edge by 200 and rounding it to the nearest integer. The file names have the format . A suffix “” signals the instances with .
The performance of our heuristic is evaluated by running it five times on each Prodhon’s 2E-LRP instance and 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 first column is the name of the 2E-LRP instances. The 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. The text in bold characters indicates the found in the four papers for the corresponding instance. The last column, “Gap (%)”, shows the difference between the best solution of our heuristic and . . The 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 find 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%. The 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 is longer than the instances with . This condition is because the difficulty 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. The 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. The locations in the few years cannot be changed when the locations of opened satellites are determined because of the large cost of change. Thus, few minutes of computer time should be consumed to find the best locations of satellites.
5.2. Performance of MHH-PR in LRPSPD
The 2E-LRPSPD can be regarded as LRPSPD when the costs of first-level edges are zero. Thus, 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). These instances contain the following characteristics: the number of customers is 8–100, and the number of satellites varies between 2 and 15. The 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 number of depots, (i.e., Author–Year-). We utilise the SN demand separation approach proposed by Salhi and Nagy [78] to generate the delivery and pickup demands of each customer for each instance. In the SN approach, ratio , where and are the coordinates of customer , and the delivery and pickup demands are generated as and , where is the original demand of customer , which is called type . Each customer’s delivery demand and pickup demand are exchanged to generate type .
The results of MHH-PR on LRPSPD instances are summarised in Table 4. The results are compared with the optimal upper bound obtained by Karagoglan et al. [34]. In this table, the first column represents the names of the LRPSPD instances. Column “DSS” shows the demand separation strategy. Column “” shows the optimal upper bound obtained in [34]. The next two columns represent the Gap value and CPU of B&C in [34]. The results of the proposed heuristic are shown in the last three columns. Column “” shows the best objective value found in five runs of the proposed heuristic. Column “Gap” shows the gap between the results of the proposed heuristic and using the following formula: . Column “CPU” shows the average CPU runtime for five runs. The 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 find 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 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 efficient. MHH-PR can find optimal/near optimal solutions in a short time period. These results show that MHH-PR is a preferable approach for solving LRPSPD instances.
5.3. 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 different 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 first column denotes the names of instances and is consistent with the format in Table 3. The MHH and MHH-PR are run for five times in each instance. The 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 five runs. Column “” shows the robustness index of the heuristic which is calculated as , where “” is the worst solution of five runs. Column “CPU” shows the average CPU time used in five runs of the heuristic. The next three columns give the results of the instances solved by the MHH-PR. The last column indicates the gap of the best solutions obtained by the MHH and MHH-PR. .
As shown in 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%). This condition shows that the MHH is relatively stable in solving -type instances. The same conclusion can be derived from the MHH-PR. The average obtained by the MHH-PR in solving the -type instances is 2.08%, which is larger than that obtained in solving -type instances (1.67%). The combined PR strategy in the MHH can find a competitive solution on the basis of the value of . The 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. The 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.
6. Case Study
In this section, we use the proposed heuristic to solve a case in the work of Belgin et al. [79]. This case is the distribution operations on 25 supermarkets (customers) of a supermarket chain in the city centre of Kayseri in Turkey. The distribution operations are conducted twice a week and are based on a single-echelon distribution system in which customers are directly supplied by the main depot. Small vehicles with 11 ton capacity are used in the current distribution system because large vehicles are not allowed to traverse urban areas. The managers of the supermarket chain decide to locate one or two satellites at the city boundary in establishing a new two-echelon distribution system. The 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 using the proposed heuristic to obtain the optimal opened satellites, first-level trips and second-level trips, we assume that the opening costs of the two candidate satellites are the same. Table 7 represents the optimal result obtained by the proposed heuristic run for five times and its comparison with the results of the single-echelon distribution and two-echelon distribution plans in [79].
The conclusions derived by analysing and comparing the results in Table 7 are summarized as follows:(1)Comparison of two-echelon distribution and first-echelon distribution systems. The travelling distance of vehicles in the two-echelon distribution system is shorter than that in the first-echelon system. The opened satellite of our optimal results is S1, and the total travelling distance of the new two-echelon system is 1510.85 km. The total travelling distance is reduced by approximately 57.03% relative to that of the single-echelon distribution system (3516.30 km). This result shows that the new two-echelon system is superior to the single-echelon system in terms of travelling distance. The 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. The 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.(2)Comparison of results under different satellite open schemes. Belgin et al. [79] investigated a two-echelon vehicle routing problem with simultaneous pickup and delivery without considering the location decision and obtained all vehicle routing plans under different satellite open schemes. We compare our two-echelon distribution system with the two-echelon system with two satellites (S1 and S2) in [79]. Opening two satellites (S1 and S2) increases the input of satellite opening and vehicle travelling costs. The travelling distance of our system with S1 is 1510.85 km, and the value increases to 1520.50 km for the two-echelon system with S1 and S2.
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 efficiency. Although two large-capacity vehicles are used, the travelling distance and cost of vehicles are lower than the vehicle fixed cost, and (2) Companies should build few satellites with large capacities. The increase in the number of satellites increases the input of satellite opening and vehicle travelling costs and the difficulty of controlling and managing the distribution operations of a company.
7. Conclusion
In this study, we investigate the 2E-LRPSPD, which is an extension of the 2E-LRP. The problem includes delivery and pickup activities at each echelon with satellites and vehicles of limited capacities. A three-index flow-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. Thus, 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. The 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 significantly 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. The results reveal that the proposed algorithm can find the best and most competitive solution. The 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 efficiency.
Several directions are provided for future research on this topic. The 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.
Appendix
A. Mathematical Formulation
The three-index flow-based mixed integer programming formulation of the 2E-LRPSPD is expressed as follows:
The objective function expressed in Equation (A.1) aims to minimise the total costs, including satellite opening costs, fixed costs of primary and secondary vehicles, and travelling costs of primary and secondary vehicles.
Equation (A.2) is used to ensure that the primary vehicle only visits the opened satellites.
Equation (A.3) is used to guarantee that each opened satellite in the primary network is served by exactly one primary vehicle.
Similarly, Equation (A.4) is used to guarantee that each customer in the secondary network is served by exactly one secondary vehicle.
Equations (A.5) and (A.6) are the flow balance constraints.
Equations (A.7) and (A.8) are used to ensure that the primary and secondary vehicles can only serve at least one trip and leaves at only one node, respectively.
Equations (A.9) and (A.10) are the elimination constraints of subtours.
Equation (A.11) is used to ensure that any two satellites in the secondary network are not connected.
Equation (A.12) is used to ensure that the allocated customers of each satellite in the secondary network are visited by a secondary vehicle.
Equations (A.13) and (A.14) are used to determine the delivery and pickup demands of satellites.
Equation (A.15) is the capacity constraint of satellites. The total pickup and total delivery demands of all customers assigned to the satellite cannot exceed the capacity of the satellite.
Equations (A.16) and (A.17) represent the loads of the primary vehicle in the primary network when it departs from the main depot and returns to the main depot, respectively.
Equations (A.18) and (A.19) represent the loads of the secondary vehicle in the secondary network when it departs from the satellite and returns to the same satellite, respectively.
Equations (A.20) and (A.21) represent the changes in vehicle dynamic load in the primary and secondary networks, respectively.
Equations (A.22) and (A.23) are the vehicle capacity constraints in the primary and secondary networks, respectively.
Equations (A.24) and (A.25) are used to ensure that each customer must be allocated to an opened satellite.
Equations (A.26) and (A.27) are used to ensure that secondary vehicles are only assigned to opened satellites.
Equations (A.28)–(A.30) are the value constraints of the decision variables.
B. The Related Tables
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
Thanks go to Ismail Karaoglan for providing instances data for the paper. This work is supported by the National Nature Science of China (61903058), the Liaoning Social Science Planning Fund, China (L19BGL006), the Guidance Program of Liaoning Provincial Research and Development Plan, China (No. 2018401002).