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

Authors: Houming Fan; Jiaxin Wu; Xin Li; Xiaodan Jiang
Journal: Journal of Advanced Transportation
 (2020)

---

## 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.


---

## Body

## 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.Figure 1
Representation of traditional 2E-LRP.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. Section2 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.Table1 summarises some the literature described previously.Table 1
Summary of related studies.
Paper
Number of echelons in the network
Location decision
Customer demand
Whether the vehicle delivery and pickup
Other constraints considered
Model
Solution method
Karaoglan et al. [34]
One
Depot
Pickup and delivery demands
√
×
Two-index flow-based formulation
Branch and cut algorithm
Yu and Lin [35]
One
Depot
Pickup and delivery demands
√
×
×
Multistart simulated annealing
Karaoglan and Altiparmak [36]
One
Depot
Pickup and delivery demands
√
×
Flow-based formulation
Memetic algorithm
Zhang et al. [37]
One
Depot
Pickup and delivery demands
√
×
Three-index formulation
Variable neighbourhood scatter search algorithm
Karaoglan et al. [38]
One
Depot
Pickup and delivery demands
√
×
Node-based formulation and flow-based formulation
A two-phase heuristic based on simulated annealing
Wang et al. [39]
One
Depot and courier transceiver point
Pickup and delivery demands
√
Fuzzy time windows constraints
Formulation with two objective function: minimize distribution operation cost and maximize service level
A hybrid heuristic algorithm combining tabu search and efficient procedures
Crainic et al. [41]
Two
Platform (depot) and satellite
Delivery demand
×
×
Three formulations, differing by the type and the number of routing variables
×
Vidović et al. [42]
Two
Depot and satellite
Pickup demand
×
×
MILP model
A heuristic solution approach
Contardo et al. [45]
Two
Platform (depot) and satellite
Delivery demand
×
×
Two-index vehicle-flow formulation
Branch and cut algorithm
Boccia et al. [46]
Two
Satellite
Delivery demand
×
×
×
TS heuristic is based on the integration of the nested approach and the two-phase iterative approach
Li et al. [48]
Two
Satellite
Delivery demand
×
×
Mixed integer linear programming model
Genetic algorithm
Dalfard et al. [47]
Two
Satellite
Delivery demand
×
Maximum route length
Mathematical model regarding the constraints of the vehicle fleet capacity and the maximum route distance
Hybrid genetic algorithm and simulated annealing
Nguyen et al. [3, 4]
Two
Satellite
Delivery demand
×
×
Two-index integer linear program and three-index integer linear program
GRASP reinforced by a learning process and path relinking, multistart iterated local search with tabu list and path relinking
Govindan et al. [56]
Two
Satellite
Delivery demand
×
Soft-time windows constraint
Multi-objective optimisation model
A multi-objective hybrid approach called MHPV, a hybrid of two known multi-objective algorithms: multi-objective particle swarm optimisation (MOPSO) and adapted multiobjective variable neighbourhood search (AMOVNS)
Wang et al. [58]
Two
Satellite
Delivery demand
×
Time windows constraint
Biobjective model minimizing costs and maximizing customer satisfaction
A three-step customer clustering based approach
Ouhader and Elkyal [59]
Two
Satellite
Delivery demand
×
Multisourcing and multiproducts
Mixed integer linear programming formulation
×
Pichka et al. [60]
Two
Satellite
Delivery demand
×
Open route
Three mixed-integer linear programming models
A two-phase hybrid heuristic
Zhao et al.[61]
Two
Satellite
Delivery demand
×
Working hours and traffic restrictions
Mathematical formulation
The cooperative approximation heuristic algorithm
Dai et al. [2]
One, two, three, and four
Depot and satellite
Delivery demand
×
Maximum route length
One-echelon, two-echelon, three-echelon, and four echelon models
A two-phase method based on improved Clarke and Wright saving algorithm
Rahmani et al. [62]
Two
Satellite
Delivery demand
√
×
×
Local search methods
Demircan-Yildiz et al. [63]
Two
Satellite
Pickup and delivery demands
√
×
Node-based formulation and flow-based formulation
×
Ghatreh and Hosseini-Motlagh [64]
Two
Satellite
Pickup and delivery demands
√
Fuzzy demand
Mixed integer linear programming model
A combined heuristic method based on simulated annealing algorithm and genetic algorithm


### 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.


## 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.Table1 summarises some the literature described previously.Table 1
Summary of related studies.
Paper
Number of echelons in the network
Location decision
Customer demand
Whether the vehicle delivery and pickup
Other constraints considered
Model
Solution method
Karaoglan et al. [34]
One
Depot
Pickup and delivery demands
√
×
Two-index flow-based formulation
Branch and cut algorithm
Yu and Lin [35]
One
Depot
Pickup and delivery demands
√
×
×
Multistart simulated annealing
Karaoglan and Altiparmak [36]
One
Depot
Pickup and delivery demands
√
×
Flow-based formulation
Memetic algorithm
Zhang et al. [37]
One
Depot
Pickup and delivery demands
√
×
Three-index formulation
Variable neighbourhood scatter search algorithm
Karaoglan et al. [38]
One
Depot
Pickup and delivery demands
√
×
Node-based formulation and flow-based formulation
A two-phase heuristic based on simulated annealing
Wang et al. [39]
One
Depot and courier transceiver point
Pickup and delivery demands
√
Fuzzy time windows constraints
Formulation with two objective function: minimize distribution operation cost and maximize service level
A hybrid heuristic algorithm combining tabu search and efficient procedures
Crainic et al. [41]
Two
Platform (depot) and satellite
Delivery demand
×
×
Three formulations, differing by the type and the number of routing variables
×
Vidović et al. [42]
Two
Depot and satellite
Pickup demand
×
×
MILP model
A heuristic solution approach
Contardo et al. [45]
Two
Platform (depot) and satellite
Delivery demand
×
×
Two-index vehicle-flow formulation
Branch and cut algorithm
Boccia et al. [46]
Two
Satellite
Delivery demand
×
×
×
TS heuristic is based on the integration of the nested approach and the two-phase iterative approach
Li et al. [48]
Two
Satellite
Delivery demand
×
×
Mixed integer linear programming model
Genetic algorithm
Dalfard et al. [47]
Two
Satellite
Delivery demand
×
Maximum route length
Mathematical model regarding the constraints of the vehicle fleet capacity and the maximum route distance
Hybrid genetic algorithm and simulated annealing
Nguyen et al. [3, 4]
Two
Satellite
Delivery demand
×
×
Two-index integer linear program and three-index integer linear program
GRASP reinforced by a learning process and path relinking, multistart iterated local search with tabu list and path relinking
Govindan et al. [56]
Two
Satellite
Delivery demand
×
Soft-time windows constraint
Multi-objective optimisation model
A multi-objective hybrid approach called MHPV, a hybrid of two known multi-objective algorithms: multi-objective particle swarm optimisation (MOPSO) and adapted multiobjective variable neighbourhood search (AMOVNS)
Wang et al. [58]
Two
Satellite
Delivery demand
×
Time windows constraint
Biobjective model minimizing costs and maximizing customer satisfaction
A three-step customer clustering based approach
Ouhader and Elkyal [59]
Two
Satellite
Delivery demand
×
Multisourcing and multiproducts
Mixed integer linear programming formulation
×
Pichka et al. [60]
Two
Satellite
Delivery demand
×
Open route
Three mixed-integer linear programming models
A two-phase hybrid heuristic
Zhao et al.[61]
Two
Satellite
Delivery demand
×
Working hours and traffic restrictions
Mathematical formulation
The cooperative approximation heuristic algorithm
Dai et al. [2]
One, two, three, and four
Depot and satellite
Delivery demand
×
Maximum route length
One-echelon, two-echelon, three-echelon, and four echelon models
A two-phase method based on improved Clarke and Wright saving algorithm
Rahmani et al. [62]
Two
Satellite
Delivery demand
√
×
×
Local search methods
Demircan-Yildiz et al. [63]
Two
Satellite
Pickup and delivery demands
√
×
Node-based formulation and flow-based formulation
×
Ghatreh and Hosseini-Motlagh [64]
Two
Satellite
Pickup and delivery demands
√
Fuzzy demand
Mixed integer linear programming model
A combined heuristic method based on simulated annealing algorithm and genetic algorithm


## 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.Figure2 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.Figure 2
Example of 2E-LRPSPD.


### 3.2. Mathematical Formulation

The 2E-LRPSPD can be defined using a complete directed graphG
=
N
,
E. The parameters and variables are introduced before presenting the formulation. Table 2 provides the parameters and variables in our mathematical formulation.Table 2
Parameters and variables used in the mathematical formulation.
Sets
N
Set of nodes,N
=
N
0
∪
N
S
∪
N
C
N
0
Subset of one main depot (node 0)
N
S
Subset ofs candidate satellites, N
S
=
1
,
2
,
⋯
,
s
N
C
Subset ofn customers, N
C
=
s
+
1
,
s
+
2
,
⋯
,
s
+
n
K
Set of vehicles,K
=
K
F
∪
K
S
K
F
Subset of primary vehicles
K
S
Subset of secondary vehicles
E
Set of valid arcs,E
=
E
F
∪
E
S
E
F
Subsets of valid arcs for first-level trips,E
1
=
i
,
j
,
i
,
j
∈
N
0
∪
N
S
,
i
≠
j
E
S
Subsets of valid arcs for second-level trips,E
2
=
l
,
m
l
,
m
∈
N
S
∪
N
C
,
l
≠
m
Indicators
i
,
j
Node in primary network,i
,
j
∈
N
0
∪
N
S
,
i
≠
j
l
,
m
Node in secondary network,l
,
m
∈
N
S
∪
N
C
,
l
≠
m
k
Vehicle,k
∈
K
Parameters
c
i
j
Nonnegative travelling cost of arci
,
j
c
l
m
Nonnegative travelling cost of arcl
,
m
C
K
F
Capacity of primary vehicle
C
K
S
Capacity of secondary vehicle
F
K
F
Fixed cost of primary vehicle
F
K
S
Fixed cost of secondary vehicle
C
S
l
Capacity of satellitel
l
∈
N
S
F
S
l
Opening cost of satellitel
l
∈
N
S
P
m
Pickup demand of customerm
m
∈
N
C
D
m
Delivery demand of customerm
m
∈
N
C
Additional variables
U
i
j
k
Dynamic load of primary vehiclek
k
∈
K
F in arc i
,
j
U
l
m
k
2
Dynamic load of secondary vehiclek
k
∈
K
S in arc l
,
m
P
l
Pickup demand of opened satellitel
l
∈
N
S
D
l
Delivery demand of opened satellitel
l
∈
N
S
Decision variables
x
i
j
k
Binary variable with a value of 1 when primary vehiclek
k
∈
K
F travels directly from node i to node j, i
,
j
∈
N
0
∪
N
S
y
l
m
k
Binary variable with a value of 1 when secondary vehiclek
k
∈
K
S travels directly from node l to node m, l
,
m
∈
N
S
∪
N
C
z
l
m
Binary variable with a value of 1 when customerm
m
∈
N
C is allocated to satellite l
l
∈
N
S
y
l
Status of satellite with a value of 1 when satellitel
l
∈
N
S is openTo better highlight the contribution of the proposed heuristic, the three-index flow-based mixed integer formulation is described in Appendix1. Figure 3 shows a flowchart with input, output and constraints.Figure 3
Flowchart with input, output, and constraints.


## 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.Figure2 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.Figure 2
Example of 2E-LRPSPD.


## 3.2. Mathematical Formulation

The 2E-LRPSPD can be defined using a complete directed graphG
=
N
,
E. The parameters and variables are introduced before presenting the formulation. Table 2 provides the parameters and variables in our mathematical formulation.Table 2
Parameters and variables used in the mathematical formulation.
Sets
N
Set of nodes,N
=
N
0
∪
N
S
∪
N
C
N
0
Subset of one main depot (node 0)
N
S
Subset ofs candidate satellites, N
S
=
1
,
2
,
⋯
,
s
N
C
Subset ofn customers, N
C
=
s
+
1
,
s
+
2
,
⋯
,
s
+
n
K
Set of vehicles,K
=
K
F
∪
K
S
K
F
Subset of primary vehicles
K
S
Subset of secondary vehicles
E
Set of valid arcs,E
=
E
F
∪
E
S
E
F
Subsets of valid arcs for first-level trips,E
1
=
i
,
j
,
i
,
j
∈
N
0
∪
N
S
,
i
≠
j
E
S
Subsets of valid arcs for second-level trips,E
2
=
l
,
m
l
,
m
∈
N
S
∪
N
C
,
l
≠
m
Indicators
i
,
j
Node in primary network,i
,
j
∈
N
0
∪
N
S
,
i
≠
j
l
,
m
Node in secondary network,l
,
m
∈
N
S
∪
N
C
,
l
≠
m
k
Vehicle,k
∈
K
Parameters
c
i
j
Nonnegative travelling cost of arci
,
j
c
l
m
Nonnegative travelling cost of arcl
,
m
C
K
F
Capacity of primary vehicle
C
K
S
Capacity of secondary vehicle
F
K
F
Fixed cost of primary vehicle
F
K
S
Fixed cost of secondary vehicle
C
S
l
Capacity of satellitel
l
∈
N
S
F
S
l
Opening cost of satellitel
l
∈
N
S
P
m
Pickup demand of customerm
m
∈
N
C
D
m
Delivery demand of customerm
m
∈
N
C
Additional variables
U
i
j
k
Dynamic load of primary vehiclek
k
∈
K
F in arc i
,
j
U
l
m
k
2
Dynamic load of secondary vehiclek
k
∈
K
S in arc l
,
m
P
l
Pickup demand of opened satellitel
l
∈
N
S
D
l
Delivery demand of opened satellitel
l
∈
N
S
Decision variables
x
i
j
k
Binary variable with a value of 1 when primary vehiclek
k
∈
K
F travels directly from node i to node j, i
,
j
∈
N
0
∪
N
S
y
l
m
k
Binary variable with a value of 1 when secondary vehiclek
k
∈
K
S travels directly from node l to node m, l
,
m
∈
N
S
∪
N
C
z
l
m
Binary variable with a value of 1 when customerm
m
∈
N
C is allocated to satellite l
l
∈
N
S
y
l
Status of satellite with a value of 1 when satellitel
l
∈
N
S is openTo better highlight the contribution of the proposed heuristic, the three-index flow-based mixed integer formulation is described in Appendix1. Figure 3 shows a flowchart with input, output and constraints.Figure 3
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 (C
H
1, C
H
2, and C
H
3), VND algorithm, two LS (L
S
1 and L
S
2) 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 executesmaxiter iterations. In each iteration, the MHH-PR starts with an initial solution generated by C
H
1, C
H
2, and C
H
3, and the best solution obtained by L
S
2 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 L
S
2
, VND, L
S
1, 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.Figure 4
Flowchart of MHH-PR.


### 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 calledC
H
1, C
H
2, and C
H
3 to obtain the initial solution for each iteration and achieve enhanced diversity. The constructive heuristics are described in detail as follows:(1)
C
H
1.C
H
1 firstly calculates the absolute value diff
m of the difference between delivery demand D
m and pickup demand P
m for each customer m
,
m
∈
N
C, diff
m
=
D
m
−
P
m. 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 diff. Let C
s denote the set of customers allocated to satellite s. 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 s, NNH starts by selecting the nearest customer to the satellite in C
s. 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 s, and another trip is constructed for s considering that the nearest customer is not visited in C
s. After all customers are visited in C
s, 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)
C
H
2.Heuristic
H
2 by Nguyen et al. [3, 4] is used as C
H
2 in this study. This heuristic is an extended NNH. C
H
2 is different from C
H
1. C
H
1 firstly assigns all customers to the nearest satellite and uses NNH to construct the trips for each satellite. C
H
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. The heuristic terminates when all customers are inserted into the trips. The NHH described in C
H
1 is applied to the first-level network to build the first-level trips.(3)
C
H
3.C
H
3 is similar to heuristic
H
3 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 C
H
1. A detailed description of this heuristic can be found in [3, 4].


### 4.2. Variable Neighbourhood Descent

VND is used to inspect neighbourhood structuresN
k
,
k
=
1
,
2
,
⋯
,
k
max. If an improving move can be found in neighbourhood structure N
k, then this move is executed, and k is reset to 1; otherwise, the next neighbourhood structure N
k
+
1 continues to be searched. VND stops when all neighbourhood structures are browsed without improvement.The VND used in this study contains four neighbourhood structures (k
max
=
4): 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 O
n
2.
(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 Exchangemoves can be explored in O
n
2.
(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 O
n
2.
(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.Figure 5
Examples of N1.Figure 6
Examples of N2.Figure 7
Examples of N3.Figure 8
Examples of N4.All N4 moves can be tested inO
n
2.


### 4.3. Local Search

This study uses two LS procedures, namely,L
S
1 and L
S
2, where L
S
1 improves the first-level trips of the solution that mainly contains two neighbourhood structures, namely, N1 and N2. L
S
2 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 L
S
2 are described in Algorithm 1. The L
S
2 phase randomly generates 2
N
S feasible solutions by using each neighbourhood structure of N5, N6, and N7. The best solution among the 6
N
S 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. L
S
2 is stopped when the current solution does not improve in 10 successive steps.Algorithm 1: LS2.
1 n
u
i
_
l
s
2
←
0
2  repeat
3  Generate 2
N
S feasible solutions using N5;
4  Generate 2
N
S feasible solutions using N6;
5  Generate 2
N
S feasible solutions using N7;
6  Select best solution S
′ among 6
N
S feasible solutions
7    if
f
S
≥
f
S
′
then
8     begin
S
←
S
′; f
S
←
f
S
′; n
u
i
_
l
s
2
←
0
end
9    else
10    n
u
i
_
l
s
2
←
n
u
i
_
l
s
2
+
1;
11    end if
12 until
n
u
i
_
l
s
2
>
10(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 O
s.
(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 O
s
2.
(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 O
s
2.


### 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 A and B.Algorithm 2: PR procedure PR.A
,
B.
1 if
A
≠
B
then
2  for
i←1 to s
do
3d
←
f
A
,
i
−
f
B
,
i;
4   if
d
>
0
then
5     Add d copies i at the end of B;
6   else if
d
<
0
then
7    Add d copies i at the end of A;
8   end if
9   end for
10 end if
11
U
←
A;
12 while
Δ
U
,
B
≠
0
do
13  seek a customer m such that p
o
s
U
,
m
≠
p
o
s
B
,
m;
14   if found then
15   swap m and U
p
o
s
B
,
m in U;
16   else
17    /Here ∀
m
∈
N
C, p
o
s
U
,
m
=
p
o
s
B
,
m;
18   seek an index j such that U
j
≠
B
j;
19    if satellite l
=
U
j occurs in B then
20   seek index z such that B
z
=
l
a
n
d
U
z
≠
l;
21     if found then
22    Swap U
j and U
z;
23     else
24     U
j
←
B
j;
25     end if
26    else
27     U
j
←
B
j;
28    end if
29   end if
30   Select p
r
o
b at random in 0
,
1;
31   if (U is unfeasible) and p
r
o
b
<
δ
then
32    repair (U);
33   end if
34   if
U feasible then
35   deduce an 2E-LRP solution S from U;
36    if
c
o
s
t
S
≤
1
+
μ
c
o
s
t
S
∗
then
37   improve S using LS2;
38    if
c
o
s
t
S
≤
1
+
η
c
o
s
t
S
∗
then
39     improve S using VND;
40      end if
41     end if
42    if
c
o
s
t
S
<
c
o
s
t
S
∗
then
43   S
∗
←
S
44    end if
45   end if
46 end whileInAlgorithm 2, the balance procedure (lines 1–10) is called when A and B differ in length. Let f
X
,
i denote the frequency of satellite i in big tour A. The balance procedure is to add copies of satellites at the end of a big tour to ensure f
A
,
i
=
f
B
,
i for each satellite i. Line 11 generates an intermediate big tour U when A and B have the same length. Let Δ
A
,
B denote the Hamming distance of A and B to guide the path from A to B
.Each iteration of thewhile loop (lines 12) generates one intermediate big tour U to reduce Δ
A
,
B. Let p
o
s
U
,
m denote the rank of a customer m in a big tour U. PR firstly looks for customers with different locations in U and B (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. C
K
F
=
30. C
S
1
=
C
S
2
=
C
S
3
=
60. Big tour U (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, S
C
=
7
,
13. 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 U is converted into a 2E-LRP solution S by adding first-level trips (line 35). S is improved using L
S
2 and updated when c
o
s
t
S
≤
1
+
μ
c
o
s
t
S
∗ (lines 36-37). S is improved using VND and updated when c
o
s
t
S
≤
1
+
η
c
o
s
t
S
∗ (lines 38-39). The best solution S
∗(lines 42–44) is updated.Figure 9
Examples ofrepair procedure.


## 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 calledC
H
1, C
H
2, and C
H
3 to obtain the initial solution for each iteration and achieve enhanced diversity. The constructive heuristics are described in detail as follows:(1)
C
H
1.C
H
1 firstly calculates the absolute value diff
m of the difference between delivery demand D
m and pickup demand P
m for each customer m
,
m
∈
N
C, diff
m
=
D
m
−
P
m. 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 diff. Let C
s denote the set of customers allocated to satellite s. 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 s, NNH starts by selecting the nearest customer to the satellite in C
s. 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 s, and another trip is constructed for s considering that the nearest customer is not visited in C
s. After all customers are visited in C
s, 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)
C
H
2.Heuristic
H
2 by Nguyen et al. [3, 4] is used as C
H
2 in this study. This heuristic is an extended NNH. C
H
2 is different from C
H
1. C
H
1 firstly assigns all customers to the nearest satellite and uses NNH to construct the trips for each satellite. C
H
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. The heuristic terminates when all customers are inserted into the trips. The NHH described in C
H
1 is applied to the first-level network to build the first-level trips.(3)
C
H
3.C
H
3 is similar to heuristic
H
3 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 C
H
1. A detailed description of this heuristic can be found in [3, 4].


## 4.2. Variable Neighbourhood Descent

VND is used to inspect neighbourhood structuresN
k
,
k
=
1
,
2
,
⋯
,
k
max. If an improving move can be found in neighbourhood structure N
k, then this move is executed, and k is reset to 1; otherwise, the next neighbourhood structure N
k
+
1 continues to be searched. VND stops when all neighbourhood structures are browsed without improvement.The VND used in this study contains four neighbourhood structures (k
max
=
4): 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 O
n
2.
(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 Exchangemoves can be explored in O
n
2.
(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 O
n
2.
(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.Figure 5
Examples of N1.Figure 6
Examples of N2.Figure 7
Examples of N3.Figure 8
Examples of N4.All N4 moves can be tested inO
n
2.


## 4.3. Local Search

This study uses two LS procedures, namely,L
S
1 and L
S
2, where L
S
1 improves the first-level trips of the solution that mainly contains two neighbourhood structures, namely, N1 and N2. L
S
2 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 L
S
2 are described in Algorithm 1. The L
S
2 phase randomly generates 2
N
S feasible solutions by using each neighbourhood structure of N5, N6, and N7. The best solution among the 6
N
S 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. L
S
2 is stopped when the current solution does not improve in 10 successive steps.Algorithm 1: LS2.
1 n
u
i
_
l
s
2
←
0
2  repeat
3  Generate 2
N
S feasible solutions using N5;
4  Generate 2
N
S feasible solutions using N6;
5  Generate 2
N
S feasible solutions using N7;
6  Select best solution S
′ among 6
N
S feasible solutions
7    if
f
S
≥
f
S
′
then
8     begin
S
←
S
′; f
S
←
f
S
′; n
u
i
_
l
s
2
←
0
end
9    else
10    n
u
i
_
l
s
2
←
n
u
i
_
l
s
2
+
1;
11    end if
12 until
n
u
i
_
l
s
2
>
10(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 O
s.
(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 O
s
2.
(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 O
s
2.


## 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 A and B.Algorithm 2: PR procedure PR.A
,
B.
1 if
A
≠
B
then
2  for
i←1 to s
do
3d
←
f
A
,
i
−
f
B
,
i;
4   if
d
>
0
then
5     Add d copies i at the end of B;
6   else if
d
<
0
then
7    Add d copies i at the end of A;
8   end if
9   end for
10 end if
11
U
←
A;
12 while
Δ
U
,
B
≠
0
do
13  seek a customer m such that p
o
s
U
,
m
≠
p
o
s
B
,
m;
14   if found then
15   swap m and U
p
o
s
B
,
m in U;
16   else
17    /Here ∀
m
∈
N
C, p
o
s
U
,
m
=
p
o
s
B
,
m;
18   seek an index j such that U
j
≠
B
j;
19    if satellite l
=
U
j occurs in B then
20   seek index z such that B
z
=
l
a
n
d
U
z
≠
l;
21     if found then
22    Swap U
j and U
z;
23     else
24     U
j
←
B
j;
25     end if
26    else
27     U
j
←
B
j;
28    end if
29   end if
30   Select p
r
o
b at random in 0
,
1;
31   if (U is unfeasible) and p
r
o
b
<
δ
then
32    repair (U);
33   end if
34   if
U feasible then
35   deduce an 2E-LRP solution S from U;
36    if
c
o
s
t
S
≤
1
+
μ
c
o
s
t
S
∗
then
37   improve S using LS2;
38    if
c
o
s
t
S
≤
1
+
η
c
o
s
t
S
∗
then
39     improve S using VND;
40      end if
41     end if
42    if
c
o
s
t
S
<
c
o
s
t
S
∗
then
43   S
∗
←
S
44    end if
45   end if
46 end whileInAlgorithm 2, the balance procedure (lines 1–10) is called when A and B differ in length. Let f
X
,
i denote the frequency of satellite i in big tour A. The balance procedure is to add copies of satellites at the end of a big tour to ensure f
A
,
i
=
f
B
,
i for each satellite i. Line 11 generates an intermediate big tour U when A and B have the same length. Let Δ
A
,
B denote the Hamming distance of A and B to guide the path from A to B
.Each iteration of thewhile loop (lines 12) generates one intermediate big tour U to reduce Δ
A
,
B. Let p
o
s
U
,
m denote the rank of a customer m in a big tour U. PR firstly looks for customers with different locations in U and B (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. C
K
F
=
30. C
S
1
=
C
S
2
=
C
S
3
=
60. Big tour U (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, S
C
=
7
,
13. 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 U is converted into a 2E-LRP solution S by adding first-level trips (line 35). S is improved using L
S
2 and updated when c
o
s
t
S
≤
1
+
μ
c
o
s
t
S
∗ (lines 36-37). S is improved using VND and updated when c
o
s
t
S
≤
1
+
η
c
o
s
t
S
∗ (lines 38-39). The best solution S
∗(lines 42–44) is updated.Figure 9
Examples ofrepair procedure.


## 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 adoptProdhon’s 2E-LRP instancesfrom [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 n
∈
20
,
50
,
100
,
200, number of satellites m
∈
5
,
10, capacity of secondary vehicle Q
=
70
,
150 and number of clusters β
=
1
,
2
,
3. 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 n
−
m
−
β. A suffix “b” signals the instances with Q
=
150.The performance of our heuristic is evaluated by running it five times on eachProdhon’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 “C
o
s
t” 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 B
e
s
t
R
e
s
u
l
t found in the four papers for the corresponding instance. The last column, “Gap (%)”, shows the difference between the best solution of our heuristic and B
e
s
t
R
e
s
u
l
t. G
a
p
=
M
H
H
_
P
R
b
e
s
t
−
B
e
s
t
R
e
s
u
l
t
/
B
e
s
t
R
e
s
u
l
t
×
100. The average of all instances for each column is given in the last row of the table.Table 3
Comparison of MHH-PR and other heuristics inProdhon’s 2E-LRP instances.
Instance
Nguyen et al. [3, 4]
Contardo et al. [45]
Pichka et al. [60]
This paper
C
o
s
t
CPU (s)
C
o
s
t
CPU (s)
C
o
s
t
CPU (s)
C
o
s
t
CPU (s)
Gap (%)
20-5-1
89075
2.4
89075
2.8
89075
2.0
89075
4.3
0.00
20-5-1b
61863
2.6
61863
0.2
61838
2.0
61838
4.8
0.00
20-5-2
85290
1.6
84478
2.8
84478
2.0
84478
5.0
0.00
20-5-2b
60838
1.4
60838
0.0
60838
3.0
60838
4.7
0.00
50-5-1
134855
9.2
130843
8.2
135008
47.0
134481
21.7
2.78
50-5-1b
101530
15.8
101530
8.7
101879
32.0
102690
9.4
1.14
50-5-2
132159
12.9
131825
37.2
132131
48.0
133666
32.0
1.40
50-5-2b
110547
18.6
110332
21.1
110395
34.0
111242
12.0
0.82
50-5-2Bis
122654
27.6
122599
43.2
122809
48.0
123286
67.7
0.56
50-5-2bBis
105776
30.0
105696
57.4
105835
34.0
107365
45.2
1.58
50-5-3
128379
16.2
128379
17.8
128379
49.0
129408
29.5
0.80
50-5-3b
104006
17.0
104006
14.2
104006
33.0
104006
22.9
0.00
100-5-1
320130
37.3
319137
445.0
323693
158.0
330437
272.6
3.54
100-5-1b
258205
32.0
257349
607.9
260754
109.0
263363
114.8
2.34
100-5-2
234179
29.6
231305
43.7
235752
157.0
235023
141.0
1.61
100-5-2b
195426
24.9
194729
602.4
195450
107.0
196897
70.3
1.11
100-5-3
245944
68.0
244194
167.9
250992
156.0
249336
111.7
2.11
100-5-3b
195254
54.3
194110
210.9
198343
108.0
199126
105.8
2.58
100-10-1
358939
336.4
358068
65.6
356669
159.0
370088
157.3
3.76
100-10-1b
302584
335.5
297167
118.7
301509
114.0
316786
146.6
6.60
100-10-2
306303
362.0
305402
163.0
308712
157.0
307599
191.9
0.72
100-10-2b
264389
294.1
265138
225.9
266692
116.0
266535
125.4
0.81
100-10-3
313249
370.8
313517
106.0
322096
160.0
321605
287.6
2.67
100-10-3b
266383
129.1
264096
181.6
275951
114.0
272119
190.0
3.04
200-10-1
554598
700.9
552816
648.7
561835
613.0
559353
1878.8
1.18
200-10-1b
452286
723.7
448236
927.6
460476
448.0
458935
906.5
2.39
200-10-2
502173
220.2
498199
339.3
506124
623.0
509329
1766.7
2.23
200-10-2b
425311
267.3
423048
1161.5
429033
424.0
431368
653.4
1.97
200-10-3
533732
632.9
534569
285.2
543312
605.0
545558
1368.0
2.22
200-10-3b
418800
323.6
404284
439.1
409961
413.0
422903
622.4
4.61
Average
246161.9
169.9
244560.9
231.8
248134.2
169.2
249957.8
312.3
1.82
The bold values in the first column is the name of the instance.n denotes number of customers. m denotes number of satellites. Q denotes capacity of secondary vehicle.  β denotes number of clusters. The names have the format n-m-β. A suffix “b” signals the instances with Q
=
150. The name of instance can not be set in bold.
The bold in other columns are the best objective values, and they need to be saved in bold.As shown in Table3, 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 Q
=
70 is longer than the instances with Q
=
150. 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 inhttp://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, N
C; and the number of depots, N
S (i.e., Author–Year-N
C
N
S). 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 r
m
=
min
x
m
/
y
m
,
y
m
/
x
m, where x
m and y
m are the coordinates of customer m, and the delivery and pickup demands are generated as D
m
=
r
m
d
m and P
m
=
d
m
−
D
m, where d
m is the original demand of customer m, which is called type X. Each customer’s delivery demand and pickup demand are exchanged to generate type Y.The results of MHH-PR on LRPSPD instances are summarised in Table4. 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 “B
e
s
t
R
e
s
u
l
t” 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 “C
o
s
t” 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 B
e
s
t
R
e
s
u
l
t using the following formula: G
a
p
=
C
o
s
t
−
B
e
s
t
R
e
s
u
l
t
/
B
e
s
t
R
e
s
u
l
t
×
100. Column “CPU” shows the average CPU runtime for five runs. The last row of the table gives the average of the columns.Table 4
Comparison  of proposed heuristic and B&C on Barreto’s LRPSPD instances.
Instance
DSS
B
e
s
t
R
e
s
u
l
t
B&C [34]
MHH-PR
Gap (%)
CPU
C
o
s
t
Gap (%)
CPU
Srivastava86_8∗2
X
625.43
0.00
0.38
625.43
0.00
1.09
Y
625.43
0.00
0.26
625.43
0.00
1.03
Perl83_12∗2
X
242.41
0.00
0.60
242.41
0.00
3.05
Y
242.41
0.00
0.55
242.41
0.00
2.21
Gaskell67_21∗5
X
454.48
0.00
13.34
454.48
0.00
7.51
Y
454.48
0.00
8.98
454.48
0.00
6.08
Gaskell67_22∗5
X
629.51
0.00
109.62
629.51
0.00
5.17
Y
629.51
0.00
64.99
629.51
0.00
5.04
Min92_27∗5
X
2998.80
0.00
91.87
2998.80
0.00
4.95
Y
2998.80
0.00
190.62
2998.80
0.00
6.61
Gaskell67_29∗5
X
490.34
0.00
6.63
490.34
0.00
4.65
Y
490.34
0.00
6.57
490.34
0.00
5.99
Gaskell67_32∗5_1
X
563.47
0.00
18.77
563.47
0.00
6.83
Y
563.47
0.00
27.31
563.47
0.00
6.98
Gaskell67_32∗5_2
X
507.03
0.00
13.72
507.03
0.00
5.88
Y
507.03
0.00
7.87
507.03
0.00
5.30
Gaskell67_36∗5
X
494.86
0.00
18.36
503.16
1.68
7.20
Y
494.86
0.00
18.13
501.05
1.25
13.62
Ghristofides69_50∗5
X
576.69
2.40
14400.00
612.58
6.22
15.38
Y
578.97
2.69
14400.00
600.25
3.67
16.42
Perl83_55∗15
X
985.56
2.94
14400.00
999.06
1.37
57.50
Y
988.06
3.14
14400.00
1002.17
1.43
55.07
Christofides69_75∗10
X
888.16
10.33
14400.00
896.51
0.94
47.96
Y
884.00
8.90
14400.00
890.80
0.77
62.06
Perl83_85∗7
X
1381.57
5.88
14400.00
1386.52
0.36
131.55
Y
1368.64
4.96
14400.00
1412.78
3.22
114.98
Daskin95_88∗8
X
399.17
11.20
14400.00
406.68
1.88
109.69
Y
401.98
0.00
14068.40
419.60
4.38
134.31
Christofides69_100∗10
X
880.18
7.45
14400.00
884.16
0.45
151.00
Y
880.89
6.17
14400.00
902.04
2.40
115.67
Average
2.20
5768.90
814.68
1.00
37.03
In this table, the first column represents the names of the LRPSPD instances. 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,N
C; and the number of depots, N
S  (i.e., Author–Year-N
C
∗
N
S). The names of instances can not be set in bold.As shown in Table4, 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 20Prodhon’s 2E-LRP instances and the SN demand separation approach to construct the instances for the 2E-LRPSPD.X and Y represent two different demand separation strategies. Tables 5 and 6 present the results of the MHH and MHH-PR in solving the X-type and Y-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 X-type (or Y-type) 2E-LRPSPD instances solved by the MHH. Column “B
e
s
t” shows the objective value of the best feasible solution in five runs. Column “R
o” shows the robustness index of the heuristic which is calculated as R
o
=
W
o
r
s
t
−
B
e
s
t
/
W
o
r
s
t
×
100, where “W
o
r
s
t” 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. G
a
p
=
B
e
s
t
M
H
H
−
B
e
s
t
M
H
H
_
P
R
/
B
e
s
t
M
H
H
×
100.Table 5
Comparison of MHH and MHH-PR inX-type 2E-LRPSPD instances.
Instance
MHH
MHH_PR
Gap
B
e
s
t
R
o
CPU
B
e
s
t
R
o
CPU
20-5-1
57728
1.87
3.11
57602
2.09
4.58
0.22
20-5-1b
42145
0.00
1.76
42145
0.00
2.79
0.00
20-5-2
59845
2.30
2.75
59845
1.47
5.66
0.00
20-5-2b
44785
0.45
2.29
44785
0.00
3.37
0.00
50-5-1
100173
1.15
16.31
98492
2.37
24.61
1.68
50-5-1b
72576
4.47
13.53
72576
4.43
19.42
0.00
50-5-2
96072
3.03
17.40
96398
2.22
21.82
-0.34
50-2-2b
76387
4.71
19.70
76272
3.01
24.56
0.15
100-5-1
242389
0.67
111.95
238508
2.02
148.22
1.60
100-5-1b
187373
1.97
100.54
183865
2.22
164.95
1.87
100-5-2
208932
1.39
130.36
206287
1.60
222.44
1.27
100-5-2b
175166
0.54
112.28
173133
1.05
247.50
1.16
100-10-1
227042
2.32
143.96
224976
2.65
238.13
0.91
100-10-1b
183507
2.39
130.15
181365
2.37
178.42
1.17
100-10-2
216237
2.15
128.35
213988
1.85
185.74
1.04
100-10-2b
188149
3.19
107.34
187568
1.29
131.90
0.31
200-10-1
404416
1.02
1584.96
398134
2.54
2366.07
1.55
200-10-1b
308476
1.63
1170.04
308543
2.58
1420.81
-0.02
200-10-2
390336
2.51
1566.20
389021
2.77
2591.71
0.34
200-10-2b
317749
4.12
955.96
314760
3.01
1270.40
0.94
Average
179974.15
2.09
315.95
178413.15
2.08
463.65
0.69
The bold values in the first column is the name of the instance.n denotes number of customers. m denotes number of satellites. Q denotes capacity of secondary vehicle.  β denotes number of clusters. The names have the format n-m-β. A suffix “b” signals the instances with Q
=
150. The name of instance may not be set in bold.Table 6
Comparison of MHH and MHH-PR inY-type 2E-LRPSPD instances.
Instance
MHH
MHH-PR
Gap
B
e
s
t
R
o
center
B
e
s
t
R
o
CPU
20-5-1
57744
1.83
3.57
57602
1.76
4.93
0.25
20-5-1b
42145
0.00
1.67
42145
0.00
2.63
0.00
20-5-2
59845
1.95
2.68
59845
1.74
6.31
0.00
20-5-2b
44785
1.11
2.26
44785
0.00
3.46
0.00
50-5-1
97701
1.93
16.56
96592
3.04
19.45
1.14
50-5-1b
74887
2.18
13.61
73586
2.97
16.38
1.74
50-5-2
97792
0.88
18.53
95291
3.25
27.61
2.56
50-2-2b
75291
3.85
19.85
76812
1.20
29.42
−2.02
100-5-1
239130
2.09
108.57
238406
1.47
159.51
0.30
100-5-1b
187729
1.39
99.26
186279
0.77
125.31
0.77
100-5-2
208256
0.86
130.06
205307
1.42
206.52
1.42
100-5-2b
174316
1.44
113.76
172837
1.53
188.09
0.85
100-10-1
227021
2.35
149.30
224100
1.73
246.36
1.29
100-10-1b
183236
2.77
127.31
183010
1.28
178.47
0.12
100-10-2
215486
0.90
146.94
213777
1.63
232.75
0.79
100-10-2b
187990
2.02
112.06
188729
0.54
135.64
−0.39
200-10-1
405217
2.27
1590.35
397631
1.89
2319.23
1.87
200-10-1b
310766
1.90
1116.04
306141
1.73
1995.94
1.49
200-10-2
396124
1.59
1516.11
386545
2.53
2123.58
2.42
200-10-2b
314432
4.34
979.19
314427
2.89
1342.10
0.00
Average
179994.65
1.88
313.39
178192.35
1.67
468.18
0.73
The bold values in the first column is the name of the instance.n denotes number of customers. m denotes number of satellites. Q denotes capacity of secondary vehicle. β denotes number of clusters. The names have the format n-m-β. A suffix “b” signals the instances with Q
=
150. The name of instance can not be set in bold.As shown in Tables5 and 6, the average R
o obtained by the MHH in solving all Y-type instances is 1.88%, which is less than that obtained in solving all X-type instances (2.09%). This condition shows that the MHH is relatively stable in solving Y-type instances. The same conclusion can be derived from the MHH-PR. The average R
o obtained by the MHH-PR in solving the X-type instances is 2.08%, which is larger than that obtained in solving Y-type instances (1.67%). The combined PR strategy in the MHH can find a competitive solution on the basis of the value of G
a
p. The improvement rates of PR for the X-type and Y-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 X-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 Y-type instances. Although the runtime increases, a good solution is obtained by the cost of consumption at a certain time.


## 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 adoptProdhon’s 2E-LRP instancesfrom [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 n
∈
20
,
50
,
100
,
200, number of satellites m
∈
5
,
10, capacity of secondary vehicle Q
=
70
,
150 and number of clusters β
=
1
,
2
,
3. 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 n
−
m
−
β. A suffix “b” signals the instances with Q
=
150.The performance of our heuristic is evaluated by running it five times on eachProdhon’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 “C
o
s
t” 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 B
e
s
t
R
e
s
u
l
t found in the four papers for the corresponding instance. The last column, “Gap (%)”, shows the difference between the best solution of our heuristic and B
e
s
t
R
e
s
u
l
t. G
a
p
=
M
H
H
_
P
R
b
e
s
t
−
B
e
s
t
R
e
s
u
l
t
/
B
e
s
t
R
e
s
u
l
t
×
100. The average of all instances for each column is given in the last row of the table.Table 3
Comparison of MHH-PR and other heuristics inProdhon’s 2E-LRP instances.
Instance
Nguyen et al. [3, 4]
Contardo et al. [45]
Pichka et al. [60]
This paper
C
o
s
t
CPU (s)
C
o
s
t
CPU (s)
C
o
s
t
CPU (s)
C
o
s
t
CPU (s)
Gap (%)
20-5-1
89075
2.4
89075
2.8
89075
2.0
89075
4.3
0.00
20-5-1b
61863
2.6
61863
0.2
61838
2.0
61838
4.8
0.00
20-5-2
85290
1.6
84478
2.8
84478
2.0
84478
5.0
0.00
20-5-2b
60838
1.4
60838
0.0
60838
3.0
60838
4.7
0.00
50-5-1
134855
9.2
130843
8.2
135008
47.0
134481
21.7
2.78
50-5-1b
101530
15.8
101530
8.7
101879
32.0
102690
9.4
1.14
50-5-2
132159
12.9
131825
37.2
132131
48.0
133666
32.0
1.40
50-5-2b
110547
18.6
110332
21.1
110395
34.0
111242
12.0
0.82
50-5-2Bis
122654
27.6
122599
43.2
122809
48.0
123286
67.7
0.56
50-5-2bBis
105776
30.0
105696
57.4
105835
34.0
107365
45.2
1.58
50-5-3
128379
16.2
128379
17.8
128379
49.0
129408
29.5
0.80
50-5-3b
104006
17.0
104006
14.2
104006
33.0
104006
22.9
0.00
100-5-1
320130
37.3
319137
445.0
323693
158.0
330437
272.6
3.54
100-5-1b
258205
32.0
257349
607.9
260754
109.0
263363
114.8
2.34
100-5-2
234179
29.6
231305
43.7
235752
157.0
235023
141.0
1.61
100-5-2b
195426
24.9
194729
602.4
195450
107.0
196897
70.3
1.11
100-5-3
245944
68.0
244194
167.9
250992
156.0
249336
111.7
2.11
100-5-3b
195254
54.3
194110
210.9
198343
108.0
199126
105.8
2.58
100-10-1
358939
336.4
358068
65.6
356669
159.0
370088
157.3
3.76
100-10-1b
302584
335.5
297167
118.7
301509
114.0
316786
146.6
6.60
100-10-2
306303
362.0
305402
163.0
308712
157.0
307599
191.9
0.72
100-10-2b
264389
294.1
265138
225.9
266692
116.0
266535
125.4
0.81
100-10-3
313249
370.8
313517
106.0
322096
160.0
321605
287.6
2.67
100-10-3b
266383
129.1
264096
181.6
275951
114.0
272119
190.0
3.04
200-10-1
554598
700.9
552816
648.7
561835
613.0
559353
1878.8
1.18
200-10-1b
452286
723.7
448236
927.6
460476
448.0
458935
906.5
2.39
200-10-2
502173
220.2
498199
339.3
506124
623.0
509329
1766.7
2.23
200-10-2b
425311
267.3
423048
1161.5
429033
424.0
431368
653.4
1.97
200-10-3
533732
632.9
534569
285.2
543312
605.0
545558
1368.0
2.22
200-10-3b
418800
323.6
404284
439.1
409961
413.0
422903
622.4
4.61
Average
246161.9
169.9
244560.9
231.8
248134.2
169.2
249957.8
312.3
1.82
The bold values in the first column is the name of the instance.n denotes number of customers. m denotes number of satellites. Q denotes capacity of secondary vehicle.  β denotes number of clusters. The names have the format n-m-β. A suffix “b” signals the instances with Q
=
150. The name of instance can not be set in bold.
The bold in other columns are the best objective values, and they need to be saved in bold.As shown in Table3, 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 Q
=
70 is longer than the instances with Q
=
150. 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 inhttp://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, N
C; and the number of depots, N
S (i.e., Author–Year-N
C
N
S). 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 r
m
=
min
x
m
/
y
m
,
y
m
/
x
m, where x
m and y
m are the coordinates of customer m, and the delivery and pickup demands are generated as D
m
=
r
m
d
m and P
m
=
d
m
−
D
m, where d
m is the original demand of customer m, which is called type X. Each customer’s delivery demand and pickup demand are exchanged to generate type Y.The results of MHH-PR on LRPSPD instances are summarised in Table4. 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 “B
e
s
t
R
e
s
u
l
t” 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 “C
o
s
t” 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 B
e
s
t
R
e
s
u
l
t using the following formula: G
a
p
=
C
o
s
t
−
B
e
s
t
R
e
s
u
l
t
/
B
e
s
t
R
e
s
u
l
t
×
100. Column “CPU” shows the average CPU runtime for five runs. The last row of the table gives the average of the columns.Table 4
Comparison  of proposed heuristic and B&C on Barreto’s LRPSPD instances.
Instance
DSS
B
e
s
t
R
e
s
u
l
t
B&C [34]
MHH-PR
Gap (%)
CPU
C
o
s
t
Gap (%)
CPU
Srivastava86_8∗2
X
625.43
0.00
0.38
625.43
0.00
1.09
Y
625.43
0.00
0.26
625.43
0.00
1.03
Perl83_12∗2
X
242.41
0.00
0.60
242.41
0.00
3.05
Y
242.41
0.00
0.55
242.41
0.00
2.21
Gaskell67_21∗5
X
454.48
0.00
13.34
454.48
0.00
7.51
Y
454.48
0.00
8.98
454.48
0.00
6.08
Gaskell67_22∗5
X
629.51
0.00
109.62
629.51
0.00
5.17
Y
629.51
0.00
64.99
629.51
0.00
5.04
Min92_27∗5
X
2998.80
0.00
91.87
2998.80
0.00
4.95
Y
2998.80
0.00
190.62
2998.80
0.00
6.61
Gaskell67_29∗5
X
490.34
0.00
6.63
490.34
0.00
4.65
Y
490.34
0.00
6.57
490.34
0.00
5.99
Gaskell67_32∗5_1
X
563.47
0.00
18.77
563.47
0.00
6.83
Y
563.47
0.00
27.31
563.47
0.00
6.98
Gaskell67_32∗5_2
X
507.03
0.00
13.72
507.03
0.00
5.88
Y
507.03
0.00
7.87
507.03
0.00
5.30
Gaskell67_36∗5
X
494.86
0.00
18.36
503.16
1.68
7.20
Y
494.86
0.00
18.13
501.05
1.25
13.62
Ghristofides69_50∗5
X
576.69
2.40
14400.00
612.58
6.22
15.38
Y
578.97
2.69
14400.00
600.25
3.67
16.42
Perl83_55∗15
X
985.56
2.94
14400.00
999.06
1.37
57.50
Y
988.06
3.14
14400.00
1002.17
1.43
55.07
Christofides69_75∗10
X
888.16
10.33
14400.00
896.51
0.94
47.96
Y
884.00
8.90
14400.00
890.80
0.77
62.06
Perl83_85∗7
X
1381.57
5.88
14400.00
1386.52
0.36
131.55
Y
1368.64
4.96
14400.00
1412.78
3.22
114.98
Daskin95_88∗8
X
399.17
11.20
14400.00
406.68
1.88
109.69
Y
401.98
0.00
14068.40
419.60
4.38
134.31
Christofides69_100∗10
X
880.18
7.45
14400.00
884.16
0.45
151.00
Y
880.89
6.17
14400.00
902.04
2.40
115.67
Average
2.20
5768.90
814.68
1.00
37.03
In this table, the first column represents the names of the LRPSPD instances. 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,N
C; and the number of depots, N
S  (i.e., Author–Year-N
C
∗
N
S). The names of instances can not be set in bold.As shown in Table4, 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 20Prodhon’s 2E-LRP instances and the SN demand separation approach to construct the instances for the 2E-LRPSPD.X and Y represent two different demand separation strategies. Tables 5 and 6 present the results of the MHH and MHH-PR in solving the X-type and Y-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 X-type (or Y-type) 2E-LRPSPD instances solved by the MHH. Column “B
e
s
t” shows the objective value of the best feasible solution in five runs. Column “R
o” shows the robustness index of the heuristic which is calculated as R
o
=
W
o
r
s
t
−
B
e
s
t
/
W
o
r
s
t
×
100, where “W
o
r
s
t” 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. G
a
p
=
B
e
s
t
M
H
H
−
B
e
s
t
M
H
H
_
P
R
/
B
e
s
t
M
H
H
×
100.Table 5
Comparison of MHH and MHH-PR inX-type 2E-LRPSPD instances.
Instance
MHH
MHH_PR
Gap
B
e
s
t
R
o
CPU
B
e
s
t
R
o
CPU
20-5-1
57728
1.87
3.11
57602
2.09
4.58
0.22
20-5-1b
42145
0.00
1.76
42145
0.00
2.79
0.00
20-5-2
59845
2.30
2.75
59845
1.47
5.66
0.00
20-5-2b
44785
0.45
2.29
44785
0.00
3.37
0.00
50-5-1
100173
1.15
16.31
98492
2.37
24.61
1.68
50-5-1b
72576
4.47
13.53
72576
4.43
19.42
0.00
50-5-2
96072
3.03
17.40
96398
2.22
21.82
-0.34
50-2-2b
76387
4.71
19.70
76272
3.01
24.56
0.15
100-5-1
242389
0.67
111.95
238508
2.02
148.22
1.60
100-5-1b
187373
1.97
100.54
183865
2.22
164.95
1.87
100-5-2
208932
1.39
130.36
206287
1.60
222.44
1.27
100-5-2b
175166
0.54
112.28
173133
1.05
247.50
1.16
100-10-1
227042
2.32
143.96
224976
2.65
238.13
0.91
100-10-1b
183507
2.39
130.15
181365
2.37
178.42
1.17
100-10-2
216237
2.15
128.35
213988
1.85
185.74
1.04
100-10-2b
188149
3.19
107.34
187568
1.29
131.90
0.31
200-10-1
404416
1.02
1584.96
398134
2.54
2366.07
1.55
200-10-1b
308476
1.63
1170.04
308543
2.58
1420.81
-0.02
200-10-2
390336
2.51
1566.20
389021
2.77
2591.71
0.34
200-10-2b
317749
4.12
955.96
314760
3.01
1270.40
0.94
Average
179974.15
2.09
315.95
178413.15
2.08
463.65
0.69
The bold values in the first column is the name of the instance.n denotes number of customers. m denotes number of satellites. Q denotes capacity of secondary vehicle.  β denotes number of clusters. The names have the format n-m-β. A suffix “b” signals the instances with Q
=
150. The name of instance may not be set in bold.Table 6
Comparison of MHH and MHH-PR inY-type 2E-LRPSPD instances.
Instance
MHH
MHH-PR
Gap
B
e
s
t
R
o
center
B
e
s
t
R
o
CPU
20-5-1
57744
1.83
3.57
57602
1.76
4.93
0.25
20-5-1b
42145
0.00
1.67
42145
0.00
2.63
0.00
20-5-2
59845
1.95
2.68
59845
1.74
6.31
0.00
20-5-2b
44785
1.11
2.26
44785
0.00
3.46
0.00
50-5-1
97701
1.93
16.56
96592
3.04
19.45
1.14
50-5-1b
74887
2.18
13.61
73586
2.97
16.38
1.74
50-5-2
97792
0.88
18.53
95291
3.25
27.61
2.56
50-2-2b
75291
3.85
19.85
76812
1.20
29.42
−2.02
100-5-1
239130
2.09
108.57
238406
1.47
159.51
0.30
100-5-1b
187729
1.39
99.26
186279
0.77
125.31
0.77
100-5-2
208256
0.86
130.06
205307
1.42
206.52
1.42
100-5-2b
174316
1.44
113.76
172837
1.53
188.09
0.85
100-10-1
227021
2.35
149.30
224100
1.73
246.36
1.29
100-10-1b
183236
2.77
127.31
183010
1.28
178.47
0.12
100-10-2
215486
0.90
146.94
213777
1.63
232.75
0.79
100-10-2b
187990
2.02
112.06
188729
0.54
135.64
−0.39
200-10-1
405217
2.27
1590.35
397631
1.89
2319.23
1.87
200-10-1b
310766
1.90
1116.04
306141
1.73
1995.94
1.49
200-10-2
396124
1.59
1516.11
386545
2.53
2123.58
2.42
200-10-2b
314432
4.34
979.19
314427
2.89
1342.10
0.00
Average
179994.65
1.88
313.39
178192.35
1.67
468.18
0.73
The bold values in the first column is the name of the instance.n denotes number of customers. m denotes number of satellites. Q denotes capacity of secondary vehicle. β denotes number of clusters. The names have the format n-m-β. A suffix “b” signals the instances with Q
=
150. The name of instance can not be set in bold.As shown in Tables5 and 6, the average R
o obtained by the MHH in solving all Y-type instances is 1.88%, which is less than that obtained in solving all X-type instances (2.09%). This condition shows that the MHH is relatively stable in solving Y-type instances. The same conclusion can be derived from the MHH-PR. The average R
o obtained by the MHH-PR in solving the X-type instances is 2.08%, which is larger than that obtained in solving Y-type instances (1.67%). The combined PR strategy in the MHH can find a competitive solution on the basis of the value of G
a
p. The improvement rates of PR for the X-type and Y-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 X-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 Y-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. Table7 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].Table 7
Optimal distribution plan.
Distribution system
Opened satellites
Travelling distance (km)
Distribution plan
Single-echelon [79]
—
3516.30
D-C4-C8-C25-C2-C1-C18-D
D-C5-C17-C16-D
D-C7-C15-C14-C19-C22-D
D-C9-C13-C21-C23-C12-D
D-C20-C11-C3-C6-C10-C24-D
Two-echelon [79]
S1
1512.90
D-S1-D
D-S1-D
S1-C17-C16-C5-S1
S1-C18-C1-C2-C8-C25-C6-S1
S1-C20-C11-C3-C10-C24-C4-S1
S1-C21-C13-C9-C7-C12-S1
S1-C23-C22-C19-C14-C15-S1
S1, S2
1520.50
D-S1-D
D-S2-D
S1-C5-C16-C17- S1
S1-C7-C9-C18-C15-S1
S1-C23-C21-C22-C19-C12- S1
S2-C6-C11-C19-C3-C24-C10-S2
S2-C25-C8-C2-C13-C1-C4-S2
Two-echelon
S1
1510.85
D-S1-D
D-S1-D
S1-C7-C9-C13-C18-C14-C19-S1
S1-C21-C22-C23-C12-C5-S1
S1-C20-C11-C3-C10-C24-C4-S1
S1-C1-C2-C6-C25-C8-C15-S1
S1-C17-C16-S1The conclusions derived by analysing and comparing the results in Table7 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.



---

*Source: 9743841-2020-01-20.xml*