# The Two-Echelon Vehicle Routing Problem with Transshipment Nodes and Occasional Drivers: Formulation and Adaptive Large Neighborhood Search Heuristic

Authors: Vincent F. Yu; Minh P. K. Nguyen; Kuza Putra; Aldy Gunawan; I Gusti Bagus Budi Dharma
Journal: Journal of Advanced Transportation
 (2022)

---

## Abstract

This research introduces a new variant of the two-echelon vehicle routing problem (2EVRP) called the two-echelon vehicle routing problem with transshipment nodes and occasional drivers (2EVRP-TN-OD). In addition to city freighters in the second-echelon network, a set of occasional drivers (ODs) is available to serve customers. ODs are the basis of a crowd-shipping system in which crowds with planned trips are willing to take detours to deliver packages in exchange for some compensation. To serve customers, ODs collect the assigned packages at either satellite served by first-echelon trucks or transshipment nodes served by city freighters. We formulate this problem as a mixed-integer nonlinear programming model and develop an adaptive large neighborhood search (ALNS) to solve it. New problem-specific destroy and repair operators and a tailored local search procedure are embedded into ALNS to deal with the problem’s unique characteristics. The experiments show that the proposed ALNS effectively solves 2EVRP-TN-OD by outperforming Gurobi in terms of both solution quality and computational time. Moreover, the experiments confirm that employing occasional drivers leads to lower operational costs. Sensitivity analyses on the characteristics of occasional drivers and the impact of transshipment nodes are presented as interesting managerial insights from 2EVRP-TN-OD.


---

## Body

## 1. Introduction

City logistics has recently become an emerging branch of supply chain management due to rapid population growth in urban areas around the world. Consequently, rising consumer demand to fulfill the needs of this population has become unavoidable, leading to an increasing number of e-commerce industries. In 2019, e-commerce sales globally reached a value of US$3.53 trillion and were projected to grow to US$6.54 trillion in 2022 [1]. This market’s sales have led to new challenges in delivery operations. While raising the number of operational vehicles is one direct solution, it results in many freight transportation issues, such as traffic congestion, higher emissions and air pollution, and noise, to mention a few.Recent technological advancements called crowd-shipping have resulted in a sharing economy-based delivery concept for increasing delivery effectiveness and reducing operational costs [2]. Crowd-shipping makes use of idle resources to perform the delivery task that would otherwise be performed by delivery companies. Several large retailers, such as Walmart and Amazon, have started to develop and implement crowd-shipping platforms [3]. These recent business practices performed by well-known enterprises show that crowd-shipping is growing in popularity as a new and promising delivery system.In order to improve the living quality in a city, a two-echelon distribution strategy is commonly implemented to reduce the volume of vehicles travelling around the urban area by setting up intermediate facilities called satellites [4]. A satellite represents a secondary facility located nearby the city center, where the freights delivered by municipal trucks from the central depots are transferred and consolidated to smaller vehicles, commonly called city freighters [5, 6]. This strategy leads to an interesting topic in vehicle routing problem variants called the two-echelon vehicle routing problem (2EVRP).Due to the rise of 2EVRP implementations and the crowd-shipping concept, one of the main aims of this research is to introduce a new variant of 2EVRP called 2EVRP with transshipment nodes and occasional drivers (2EVRP-TN-OD). The first-echelon network connects the central depot with satellites, and the second-echelon network connects satellites with customers and transshipment nodes (TNs). In this research, every occasional driver (OD) has his/her origin and destination nodes, available capacity, and flexibility. ODs are able to utilize their available capacity and collect the assigned packages from either a satellite or a TN. Moreover, every employed OD receives a particular amount of compensation, depending on the total distances travelled. The contributions of this research work are summarized as follows:(1)
Introduce a new variant of 2EVRP called the two-echelon vehicle routing problem with transshipment nodes and occasional drivers (2EVRP-TN-OD).(2)
Formulate a mixed-integer nonlinear programming (MINLP) for 2EVRP-TN-OD.(3)
Develop an adaptive large neighborhood search algorithm and show the effectiveness of the proposed algorithm.(4)
Perform sensitivity analyses on important characteristics related to the crowd-shipping concept to better understand the impact of the integration into the two-echelon distribution network.The remaining parts of this article are described as follows: Section2 provides a description of previous related works. Section 3 explains the problem definition and formulation. Section 4 describes the detailed procedure of our developed ALNS. Section 5 shows the generation of benchmark instances, the computational results of ALNS, and the sensitivity analyses related to 2EVRP-TN-OD. Lastly, Section 6 concludes the findings of this work and possible future direction.


## 2. Literature Review

### 2.1. Two Echelon Vehicle Routing Problem

Crainic et al. [7] first introduced 2EVRP by considering the setting up of a two-tier distribution network. The model aims to mainly address the day-before-planning problem with several key decisions related to freight assignment to the satellites as well as the routing and scheduling of vehicles that operate in both first- and second-echelon networks. Perboli et al. [8] pioneered the method’s development by addressing the basic problem of 2EVRP, where time windows and satellite synchronizations are not considered. This problem is called a two-echelon capacitated vehicle routing problem (2E-CVRP). Recent work has focused on integrating trends and realistic scenarios into the 2EVRP structure. Wang et al. [9] considered the environmental impact of 2EVRP. In addition, Belgin et al. [10] developed an extension of 2EVRP by incorporating not only delivery but also pickup demand.Due to the rise in the utilization of electric vehicles, several works tackled similar networks by considering electric vehicles [11, 12]. Anderluh et al. [13] explicitly integrated emission minimization by formulating 2EVRP into a multiobjective problem in which both operational and emission costs are considered. In light of rapid technological development, several helpful innovative alternatives have become available, such as parcel lockers, crowd-shipping, delivery robots, and drones. These technology-enabled options led to the fruitful developments of 2EVRP [14–18]. Mühlbauer and Fontaine [19] adopted swap containers for solving 2E-CVRP, where the first echelon’s requests are fulfilled by a set of vans and the second echelon’s requests are serviced by a set of cargo bicycles or city freighters provided by the logistic companies. Various studies were collected and summarized in the review of 2EVRP variants by Sluijk et al. [20]. It can be seen that the existing literature has rarely considered heterogeneous vehicles, and so far, only Yu et al. [17] have introduced occasional drivers to 2EVRP with covering locations.


### 2.2. Location-Routing Problem

The location-routing problem (LRP) is an extension of VRP that involves two key decisions: the selection of depots and the routing of vehicles originating from a particular depot [21]. The network setting of LRP is similar to the second-echelon network of 2EVRP. The main differences between LRP and 2EVRP can be found in their objective functions and network structures. Various works have proposed solution methods such as heuristics and exact algorithms [22–26]. Schneider and Löffler [27] constructed a three-phase tree-based search algorithm for solving CLRP. Their numerical analysis showed the effectiveness of the algorithm in solving LRPs. New variants of LRP were also proposed to cope with emerging challenges and trends. Karaoglan et al. [28] set up a model and a heuristic to deal with LRP under simultaneous pickup and delivery to handle the existence of returned products from customers. Yu and Lin [29] tackled LRP under the condition that all operational fleets are supplied by a third-party logistics company. Nowadays, due to the efforts to reduce the number of operating vehicles within cities, LRP has been extended to two-tier distribution networks, namely, 2ELRP [4]. Several recent works have addressed the case of 2ELRP under specific real-world characteristics [30, 31]. Arnold and Sörensen [32] proposed a progressive filtering heuristic adopting the regret Clarke-Wright and knowledge-guided local search heuristic for solving 2ELRP and its special case, the single truck and trailer routing problem with satellite depots (STTRPSD), which have a configuration similar to the settings of 2EVRP.


### 2.3. Crowdshipping

The trend of providing crowd-shipping models and analyzing the positive impacts of considering this concept gives rise to a new VRP variant. Archetti et al. [33] pioneered this strand of research by introducing the vehicle routing problem with occasional drivers (VRPOD), where deliveries are conducted by two types of fleets: company vehicles and ODs. Macrina et al. [34] introduced three different characteristics to the VRPOD problem: time window constraints, multiple deliveries, and split delivery that can be performed by ODs. Chen et al. [35] proposed a multihop problem where parcels can be delivered using multiple consecutive tours to reach their destination within a designed time window and a maximum rate of detour compared to the drivers’ shortest paths.Based on the definition of crowdsourcing, Sampaio et al. [36] presented a system with a large number of drivers, a uniquely owned depot, and short working shifts when performing the pickup and delivery activities. These drivers do not have a desired destination, and they will return to the depot after delivering the requests. Moreover, they limit the number of transfers to once per request only, which distinguishes their work from the proposed study of Chen et al. [35]. These configurations have turned their problem into the classic pickup and delivery problem with time windows and transfer (PDPTW-T). Voigt and Kuhn [37] constructed their model based on the work of Sampaio et al. [36] with a modification in the shipping fleet team by categorizing them into occasional drivers and regular drivers. In their settings, the latter team will be the main force handling the delivery operation, and the former ones will accept the tasks based on their convenience. To sum up, recent works on the crowd-shipping system consider the existence of transshipment nodes so that ODs have more options for picking up the assigned demand [38–40]. The main advantage of a crowd-shipping system, as mentioned in all the aforementioned works, is savings on operational costs.


### 2.4. Transshipment Nodes

Transshipment nodes (TNs) are considered important factors when handling the vehicle routing problem with crowd shipping and are usually referred to as intermediate depots [37, 40]. The benefits TNs bring to the context of ODs have been proposed by various scholars. Macrina et al. [40] believed the adoption of these intermediate transfer points would gain interest and relax the barriers for more ODs to performing delivery tasks due to the advantages in the geography context, with a total savings cost of up to 31% compared to the nontransshipment nodes model in VRPOD. This rate will significantly increase in proportion to the number of available TNs and their capacities. In addition, by allowing the ODs to visit the depot, a saving of 6% in the consumed operation costs has been found in their study. Voigt and Kuhn [37] proposed a similar model to our work with the involvement of both regular drivers and ODs in their study and performed some experiments with several scenarios: VRPOD without TNs (T0), VRPOD with random TNs (TR), and VRPOD with fixed assigned TNs (T4). Though there is no significant difference in the outcomes when comparing T0 and TR in scenarios with 100 requests and 100 occasional drivers, other computational results show that the adoption of TNs usually has huge advantages in the consumed cost when compared with the reverse scenarios in their study. Yu et al. [41] also confirmed the benefit to total operation cost with the setup of TNs in the crowd shipping delivery system with four designed scenarios: alternative delivery options are applied/not applied, and TNs are applied/not applied. The mentioned works have proven the effectiveness of TNs in advancing savings in VRP in the crowd-shipping context and are worth considering when working with the crowd-shipping system.


### 2.5. Main Contributions

Our work contributes to the 2EVRP literature with the adoption of crowdsourcing and heterogeneous vehicles, which cover occasional drivers, in picking up and delivering customers’ requests. For a prior contribution, the problem settings of our work and those of Mühlbauer and Fontaine [19] are quite similar; however, we support the delivery activities of the second echelon by integrating a group of ODs. Although 2EVRPTW-OD was introduced by Yu et al. [17], they dealt with the covering locations while we handled the issue of transshipment nodes. Our latter contribution is addressed through the involvement of occasional drivers and city freighters for delivery purposes in our work’s second echelon.To our knowledge, existing studies on the crowd-shipping problem have also not yet dealt with 2EVRP. Both Huang and Ardiansyah [38] and Macrina et al. [40] handled VRP-CS with transfer locations, where regular drivers mainly serve both transshipment nodes and customers, and customers are serviced by both types of drivers: regular and occasional. Those works match our second-echelon work; however, the first-echelon activities are not analyzed in their works. Moreover, our work distinguishes between the depot-satellite delivery trucks and the satellite-transshipment node city freighters, which also have not been analyzed in previous studies. This is our addition to existing studies related to crowd-shipping problems. For more details, our paper refers to regular drivers as city freighters and first-echelon trucks for better discrimination in description and modelling. A table is given in the appendix to summarize and analyze the difference between our work and existing studies.


## 2.1. Two Echelon Vehicle Routing Problem

Crainic et al. [7] first introduced 2EVRP by considering the setting up of a two-tier distribution network. The model aims to mainly address the day-before-planning problem with several key decisions related to freight assignment to the satellites as well as the routing and scheduling of vehicles that operate in both first- and second-echelon networks. Perboli et al. [8] pioneered the method’s development by addressing the basic problem of 2EVRP, where time windows and satellite synchronizations are not considered. This problem is called a two-echelon capacitated vehicle routing problem (2E-CVRP). Recent work has focused on integrating trends and realistic scenarios into the 2EVRP structure. Wang et al. [9] considered the environmental impact of 2EVRP. In addition, Belgin et al. [10] developed an extension of 2EVRP by incorporating not only delivery but also pickup demand.Due to the rise in the utilization of electric vehicles, several works tackled similar networks by considering electric vehicles [11, 12]. Anderluh et al. [13] explicitly integrated emission minimization by formulating 2EVRP into a multiobjective problem in which both operational and emission costs are considered. In light of rapid technological development, several helpful innovative alternatives have become available, such as parcel lockers, crowd-shipping, delivery robots, and drones. These technology-enabled options led to the fruitful developments of 2EVRP [14–18]. Mühlbauer and Fontaine [19] adopted swap containers for solving 2E-CVRP, where the first echelon’s requests are fulfilled by a set of vans and the second echelon’s requests are serviced by a set of cargo bicycles or city freighters provided by the logistic companies. Various studies were collected and summarized in the review of 2EVRP variants by Sluijk et al. [20]. It can be seen that the existing literature has rarely considered heterogeneous vehicles, and so far, only Yu et al. [17] have introduced occasional drivers to 2EVRP with covering locations.


## 2.2. Location-Routing Problem

The location-routing problem (LRP) is an extension of VRP that involves two key decisions: the selection of depots and the routing of vehicles originating from a particular depot [21]. The network setting of LRP is similar to the second-echelon network of 2EVRP. The main differences between LRP and 2EVRP can be found in their objective functions and network structures. Various works have proposed solution methods such as heuristics and exact algorithms [22–26]. Schneider and Löffler [27] constructed a three-phase tree-based search algorithm for solving CLRP. Their numerical analysis showed the effectiveness of the algorithm in solving LRPs. New variants of LRP were also proposed to cope with emerging challenges and trends. Karaoglan et al. [28] set up a model and a heuristic to deal with LRP under simultaneous pickup and delivery to handle the existence of returned products from customers. Yu and Lin [29] tackled LRP under the condition that all operational fleets are supplied by a third-party logistics company. Nowadays, due to the efforts to reduce the number of operating vehicles within cities, LRP has been extended to two-tier distribution networks, namely, 2ELRP [4]. Several recent works have addressed the case of 2ELRP under specific real-world characteristics [30, 31]. Arnold and Sörensen [32] proposed a progressive filtering heuristic adopting the regret Clarke-Wright and knowledge-guided local search heuristic for solving 2ELRP and its special case, the single truck and trailer routing problem with satellite depots (STTRPSD), which have a configuration similar to the settings of 2EVRP.


## 2.3. Crowdshipping

The trend of providing crowd-shipping models and analyzing the positive impacts of considering this concept gives rise to a new VRP variant. Archetti et al. [33] pioneered this strand of research by introducing the vehicle routing problem with occasional drivers (VRPOD), where deliveries are conducted by two types of fleets: company vehicles and ODs. Macrina et al. [34] introduced three different characteristics to the VRPOD problem: time window constraints, multiple deliveries, and split delivery that can be performed by ODs. Chen et al. [35] proposed a multihop problem where parcels can be delivered using multiple consecutive tours to reach their destination within a designed time window and a maximum rate of detour compared to the drivers’ shortest paths.Based on the definition of crowdsourcing, Sampaio et al. [36] presented a system with a large number of drivers, a uniquely owned depot, and short working shifts when performing the pickup and delivery activities. These drivers do not have a desired destination, and they will return to the depot after delivering the requests. Moreover, they limit the number of transfers to once per request only, which distinguishes their work from the proposed study of Chen et al. [35]. These configurations have turned their problem into the classic pickup and delivery problem with time windows and transfer (PDPTW-T). Voigt and Kuhn [37] constructed their model based on the work of Sampaio et al. [36] with a modification in the shipping fleet team by categorizing them into occasional drivers and regular drivers. In their settings, the latter team will be the main force handling the delivery operation, and the former ones will accept the tasks based on their convenience. To sum up, recent works on the crowd-shipping system consider the existence of transshipment nodes so that ODs have more options for picking up the assigned demand [38–40]. The main advantage of a crowd-shipping system, as mentioned in all the aforementioned works, is savings on operational costs.


## 2.4. Transshipment Nodes

Transshipment nodes (TNs) are considered important factors when handling the vehicle routing problem with crowd shipping and are usually referred to as intermediate depots [37, 40]. The benefits TNs bring to the context of ODs have been proposed by various scholars. Macrina et al. [40] believed the adoption of these intermediate transfer points would gain interest and relax the barriers for more ODs to performing delivery tasks due to the advantages in the geography context, with a total savings cost of up to 31% compared to the nontransshipment nodes model in VRPOD. This rate will significantly increase in proportion to the number of available TNs and their capacities. In addition, by allowing the ODs to visit the depot, a saving of 6% in the consumed operation costs has been found in their study. Voigt and Kuhn [37] proposed a similar model to our work with the involvement of both regular drivers and ODs in their study and performed some experiments with several scenarios: VRPOD without TNs (T0), VRPOD with random TNs (TR), and VRPOD with fixed assigned TNs (T4). Though there is no significant difference in the outcomes when comparing T0 and TR in scenarios with 100 requests and 100 occasional drivers, other computational results show that the adoption of TNs usually has huge advantages in the consumed cost when compared with the reverse scenarios in their study. Yu et al. [41] also confirmed the benefit to total operation cost with the setup of TNs in the crowd shipping delivery system with four designed scenarios: alternative delivery options are applied/not applied, and TNs are applied/not applied. The mentioned works have proven the effectiveness of TNs in advancing savings in VRP in the crowd-shipping context and are worth considering when working with the crowd-shipping system.


## 2.5. Main Contributions

Our work contributes to the 2EVRP literature with the adoption of crowdsourcing and heterogeneous vehicles, which cover occasional drivers, in picking up and delivering customers’ requests. For a prior contribution, the problem settings of our work and those of Mühlbauer and Fontaine [19] are quite similar; however, we support the delivery activities of the second echelon by integrating a group of ODs. Although 2EVRPTW-OD was introduced by Yu et al. [17], they dealt with the covering locations while we handled the issue of transshipment nodes. Our latter contribution is addressed through the involvement of occasional drivers and city freighters for delivery purposes in our work’s second echelon.To our knowledge, existing studies on the crowd-shipping problem have also not yet dealt with 2EVRP. Both Huang and Ardiansyah [38] and Macrina et al. [40] handled VRP-CS with transfer locations, where regular drivers mainly serve both transshipment nodes and customers, and customers are serviced by both types of drivers: regular and occasional. Those works match our second-echelon work; however, the first-echelon activities are not analyzed in their works. Moreover, our work distinguishes between the depot-satellite delivery trucks and the satellite-transshipment node city freighters, which also have not been analyzed in previous studies. This is our addition to existing studies related to crowd-shipping problems. For more details, our paper refers to regular drivers as city freighters and first-echelon trucks for better discrimination in description and modelling. A table is given in the appendix to summarize and analyze the difference between our work and existing studies.


## 3. Problem Description

Similar to the classical 2EVRP addressed by Hemmelmayr et al. [24] and Breunig et al. [42], we model 2EVRP-TN-OD where the first-echelon network connects the central depot with satellites and the second-echelon network connects the satellites with customers. In 2EVRP-TN-OD, a set of TNs exists in the second-echelon network, serving as places for transferring goods from second-echelon vehicles (city freighters) to ODs. Figure 1 illustrates a possible solution for 2EVRP-TN-OD.Figure 1
An illustration of 2EVRP-TN-OD.In the first echelon of the distribution network in 2EVRP-TN-OD, there exists a set of homogeneous first-echelon trucks serving satellites. The demands of a satellite are the total demands of the customers assigned to the satellite. In the second echelon, there are two types of fleets for serving customers: a set of homogeneous city freighters and a set of heterogeneous ODs. Each OD has the willingness to serve one or more customers depending on his/her available capacity. The willingness is defined as a maximum extra distance that an OD willingly travels.Each customer has an amount of demand that must be fulfilled by one visit by a vehicle, either a city freighter or an OD. While no split delivery is allowed in the second echelon, 2EVRP-TN-OD adopts the split delivery of the first echelon of 2EVRP, allowing a particular satellite to be visited by more than one first-echelon truck. The first-echelon trucks start and finish at the depot, while city freighters start to visit customers or TNs from a particular satellite and end at the same satellite. Each OD starts at his/her origin, ends at his/her destination, and is able to collect demands either at a satellite or a TN before serving the assigned customers. In this work, our main objective is to minimize the total operation cost for delivery activities, including: (1) the travel costs of the first-echelon trucks, second-echelon city freighters, and ODs; and (2) the employment costs accounted for by ODs.An overview of our proposed work is illustrated in Figure1. There are 2 transshipment nodes, 3 satellites, 5 customer nodes, 1 vehicle serving the first echelon, 2 ODs and 3 city freighters fulfilling requests at the second echelon, and one depot. From the main depot, a large capacity vehicle serves two satellites (1 and 3), and then there are two city freighters that will come to the first satellite to pickup the requests, namely, freighters 1 and 2. At the same time, at the third satellite, an OD departs from origin node 1 and collects the request to deliver to the customer at node 7 before ending at destination node 1. In addition, another freighter 3 visits this satellite to pick up the requests and then deliver them to customers at Nodes 6, 10, 11, 15, and 18. Back at the first satellite, while city freighter 1 will normally deliver requests to customers at nodes 5, 9, 12, and 17, city freighter 2 is assigned to fulfill customers’ requests at nodes 8, 13, and 14 in addition to delivering the remaining requests to transshipment node 19 (green square) for later pickup and delivery activities performed by OD 2. This OD, who departs from original location 2, visits transshipment node 19 to collect requests and then continues delivering them to customers at nodes 4 and 16 on his travel route to the respected destination node 2. The other delivery activities performed by other ODs or freighters are similar to those performed by OD2 and freighter 1 in this illustrated example. To sum up, in this illustrated example, not all satellites and transshipment nodes are used for delivering requests to customers; one vehicle is used for servicing requests at satellites in the first echelon, while 2 occasional drivers and 3 city freighters perform the requests’ delivery to end-customers at the second echelon through the utilization of one transshipment node out of the two existing ones.


## 4. A Mixed-Integer Nonlinear Program Formulation

The problem is formally defined on a directed graphG=N,A, where N represents a set of vertices and A denotes a set of arcs. Here, N consists of a depot {0}, the set of satellites Ns, the set of customers Nc, the set of TNs Nr, and the set of origin and destination nodes of available ODs, NO and ND. There are two sets of arcs in A–that is, A1:=i,j|i,j∈NNc∪Nr∪∪NO∪ND,i≠j and A2:=A\A1, consecutively representing the arcs in the first- and second-echelons. For all arcs i,j in A, dij is defined as the distance from vertex i to vertex j, where tij=dij. In addition, several additional sets are defined; i.e., N1=0∪Ns, N2=Nc∪Nr∪NO∪ND, Np=Ns∪Nr, and Nv=Nc∪Nr, for the sake of simplicity of the formulation.Each customeri∈Nc has a nonnegative demand δi. A set of K1 homogeneous first-echelon trucks is available at the depot, each having a maximum capacity of Q1. In the second-echelon network, there exists a set of K2 homogeneous city freighters that can be positioned at any existing satellite and a set of KOD ODs. Each OD k∈KOD has his/her own origin nko stored in NO and his/her own destination nkD stored in ND. Each city freighter k∈K2 has a maximum capacity of Q2, while each OD k∈KOD has a maximum capacity of QkOD.LetΔk be the original distance between the origin and destination nodes of an OD k∈KOD. OD k is willing to serve the customers as long as the maximum extra-travelled distance is not greater than μΔk. The cost of employing an OD consists of two components: (1) a fixed cost H and (2) a distance-related variable cost ρOD. In addition, both first-echelon trucks and city freighters have distance-related variable costs cd. Finally, the following MINLP describes 2EVRP-TN-OD.


### 4.1. Sets

Ns Set of satellites.Nr Set of available TNs.Nc Set of customers.NO Set of origin nodes of ODs.ND Set of destination nodes of ODs.K1 Set of first-echelon trucks.K2 Set of city freighters.KOD Set of ODs.


### 4.2. Additional Sets

N1 Set of all first-echelon nodes, N1=0∪NsN2 Set of all second-echelon nodes, N2=Ns∪NC∪NrNp Set of all ODs’ pickup point nodes, Np=Ns∪NrNv Set of all nodes that are able to be visited by city freighters excluding the satellites, Nv=Nr∪Nc


### 4.3. Parameters

dij Travel distance from node i to node j, where i,j∈Aδi Demand of customer i, where i∈NcQiOD Maximum available capacity of OD i, where i∈KODΔi Distance between origin and destination nodes of OD i, where i∈KODQ1 Capacity of a first-echelon truck.Qs Capacity of a satellite.Q2 Capacity of a city freighter.Qr Capacity of a TN.H Fixed cost of employing an OD.ρOD Compensation per distance unit travelled by an OD.cd Variable cost per distance unit for either a first-echelon truck or a city freighter.μ A multiplier for calculating the extra-travelled distance of an OD.


### 4.4. Decision Variables

xijk Binary variables representing the selection of arc i,j to be traversed by first-echelon truck k∈K1; 1 if first-echelon truck k traverses through arc i,j and 0 otherwise.yijk Binary variables representing the selection of arc i,j to be traversed by city freighter k∈K2; 1 if city freighter k traverses through arc i,j and 0 otherwise.zijk Binary variables representing the selection of arc i,j to be traversed by OD k∈KOD; 1 if OD k traverses through arc i,j and 0 otherwise.αpk Binary variable indicating the assignment of OD k to pickup point p∈Np; 1 if OD k visits pickup point i and 0 otherwise.βik Binary variable representing whether customer i∈Nc is the last node visited by OD k∈KOD; 1 if OD k goes to the destination after visiting customer i and 0 otherwise.σkOD Total demands assigned to OD k∈KOD, σkOD≥0σk2 Total demands assigned to city freighter k∈K2, σk2≥0σk1 Total demands assigned to first-echelon truck k∈K1, σk1≥0liM Total demands assigned to TN i∈Nr, liM≥0liS Total demands assigned to satellite i∈NS, liS≥0qik1 Remaining loads carried by first-echelon truck k∈K1 when visiting node i∈NS, qik1≥0qik2 Remaining loads carried by city freighter k∈K2 when visiting node i∈Nv, qik2≥0qikOD Remaining loads carried by OD k∈KOD when visiting node i∈Nc, qikOD≥0


### 4.5. Objective Function

(1)MinZ=∑k∈K1∑i,j∈N1cddijxijk+∑k∈K2∑i,j∈N2cddijyijk+∑k∈KOD∑i∈nko∪N2j∈N2∪nkDρODdijzijk+∑k∈KOD∑p∈NpHαpk.Subject to(2)∑j∈N1x0jk≤1∀k∈K1,(3)∑j∈N1xj0k≤1∀k∈K1,(4)∑i∈N1xijk≤1∀j∈NS,∀k∈K1,(5)∑j∈N1xijk−∑j∈N1xjik=0i∈N1,∀k∈K1,(6)∑k∈K2∑j∈Nvypjkσk2+∑k∈KODαpkσkOD≤lpSp∈NS,(7)∑i∈N1∑j∈NSxijkljS≤σk1∀k∈K1,(8)σk1≤Q1∀k∈K1,(9)σk1≥qjk1+M1−x0jk∀j∈NS,∀k∈K1,(10)qik1−liS≤qjk1+M1−xijk∀i∈NS,∀j∈NS,∀k∈K1,(11)∑j∈Nvyijk≤∑l∈K1∑j∈N1xjili∈NS,∀k∈K2,(12)∑j∈NCzijk≤∑l∈K1∑j∈N1xjili∈NS,∀k∈KOD.The objective of 2EVRP-TN-OD in (1) is to minimize the total cost of the delivery operation. The first, second, and third terms consecutively address the travelling costs of first-echelon trucks, city freighters, and the employed ODs, while the last one represents the payments for employing these ODs. Constraints (2) to (12) focus on managing the characteristics of the first-echelon network. Constraints (2) to (4) ensure that for each first-echelon truck, every satellite can only be visited at most once. Constraint (5) represents the conservation-of-flow constraint of every first-echelon truck. Constraint (6) defines the total demands assigned to each satellite. Constraint (7) defines the total demands carried by each first-echelon truck, while constraint (8) limits the maximum number of demands that can be carried by each first-echelon truck. Constraints (9) and (10) track the remaining demands carried by each first-echelon truck while visiting a node. Constraint (11) ensures that a city freighter can originate from a satellite if and only if the satellite is served by first-echelon trucks. Similarly, constraint (12) guarantees that each OD can collect demands at a satellite if and only if the satellite is served by at least one first-echelon truck.(13)∑k∈K2∑i∈N2yijk+∑k∈KOD∑i∈N2zijk=1j∈NC,(14)αpk≤∑l∈K2∑i∈N2yiplp∈Nr,k∈KOD,(15)αpk≤∑l∈K2∑i∈N1xiplp∈NS,k∈KOD,(16)∑p∈Npαpk≤1k∈KOD,(17)∑i∈NCβik=∑p∈Npαpkk∈KOD,(18)∑i∈NCzpik=αpkp∈Np,k∈KOD,(19)∑i∈nkozipk=αpkp∈Np,k∈KOD,(20)∑j∈nkDzijk=βiki∈NC,k∈KOD,(21)∑j∈Np∪NCzjik−∑j∈NC∪nkDzijk=0πi∈NC,k∈KOD,(22)σkOD≤qjkOD+M1−zpjkj∈NC,p∈Np,k∈KOD,(23)qikOD−δi≤qjkOD+M1−zijki,j∈NC,k∈KOD,(24)∑i∈N2∑j∈NCzijkδj≤σkODk∈KOD,(25)σkOD≤QkODk∈KOD,(26)∑p∈Npdnkopαpk+∑i∈Np∪NC∑j∈NCdijzijk+∑j∈NCdjnkDβjk≤1+μΔkk∈KOD.Constraint (13) ensures that either an OD or a city freighter must serve each customer. Constraints (14) to (26) focus on managing the characteristics of ODs’ routes. In particular, constraints (14) to (15) define the assignment of each OD to a pickup point, either a satellite or a TN. Constraint (16) ensures that each OD can only visit one pickup point to take the assigned demands. Constraint (17) ensures that each OD needs to visit a customer before reaching his/her destination if the OD is assigned a demands. Constraints (18) to (21) ensure the flow-in and flow-out of routes for an OD. Constraints (22) and (23) track the remaining demands carried by an OD while visiting customer nodes. Constraint (24) calculates the total demands carried by an OD, and constraint (25) defines the maximum capacity of an OD. Lastly, constraint (26) ensures that the total distance travelled by an OD is not greater than the maximum willingness of the OD.(27)∑i∈NS∑j∈Nvyijk≤1k∈K2,(28)∑i∈NS∪Nvyijk−∑i∈NS∪Nvyjik=0j∈Nv,k∈K2,(29)∑k∈K2∑i∈N2yijk≤1j∈Nr,(30)∑k∈KODαpkσkOD≤lpMp∈Nr,(31)liM≤Qri∈Nr,(32)∑i∈N2∑j∈NCyijkδj+∑i∈N2∑j∈NryijkljM≤σk2k∈K2,(33)σk2≤Q2k∈K2,(34)σk2≤qjk2+M1−yijki∈NS,j∈Nv,k∈K2,(35)qik2−δi≤qjk2+M1−yijki∈NC,j∈Nv,k∈K2,(36)qik2−liM≤qjk2+M1−yijklimx⟶∞i∈Nr,j∈Nv,k∈K2.Constraints (27) to (36) focus on the route characteristics of the city freighters. Constraint (27) ensures that each city freighter can only originate from one satellite. Constraint (28) ensures the flow-in and flow-out of each city freighter. Constraint (29) guarantees that each TN can only be visited at most once. Constraint (30) defines the number of demands delivered to a TN by a city freighter. Constraint (31) limits the maximum allowable demands delivered to a TN. Constraint (32) defines the total demand carried by a city freighter from a satellite, while the maximum allowable demands that can be carried by a city freighter are limited by constraint (33). Moreover, constraints (34) to (36) track the remaining demands carried by a city freighter while visiting a node, either a customer node or a TN. Finally, the domains of the decision variables are given by constraints (37)–(49).(37)xijk∈0,1∀i∈N1,∀j∈N1,∀k∈K1,(38)yijk∈0,1∀i∈N2,∀j∈N2,∀k∈K2,(39)zijk∈0,1∀i∈nko∪N2,∀j∈N2∪nkD,∀k∈KOD,(40)αpk∈0,1∀p∈Np,∀k∈KOD,(41)βik∈0,1∀i∈Nc,∀k∈KOD,(42)σk1≥0∀k∈K1,(43)σk2≥0∀k∈K2,(44)σkOD≥0∀k∈KOD,(45)liM≥0∀i∈Nr,(46)liS≥0∀i∈NS,(47)qik1≥0∀i∈NS,∀k∈K1,(48)qik2≥0∀i∈Nv,∀k∈K2,(49)qikOD≥0∀i∈NC,∀k∈KOD.The formulation belongs to a mixed-integer nonlinear program due to the multiplication of two types of decision variables in several constraints, i.e., constraints (6), (7), (30), and (32).


## 4.1. Sets

Ns Set of satellites.Nr Set of available TNs.Nc Set of customers.NO Set of origin nodes of ODs.ND Set of destination nodes of ODs.K1 Set of first-echelon trucks.K2 Set of city freighters.KOD Set of ODs.


## 4.2. Additional Sets

N1 Set of all first-echelon nodes, N1=0∪NsN2 Set of all second-echelon nodes, N2=Ns∪NC∪NrNp Set of all ODs’ pickup point nodes, Np=Ns∪NrNv Set of all nodes that are able to be visited by city freighters excluding the satellites, Nv=Nr∪Nc


## 4.3. Parameters

dij Travel distance from node i to node j, where i,j∈Aδi Demand of customer i, where i∈NcQiOD Maximum available capacity of OD i, where i∈KODΔi Distance between origin and destination nodes of OD i, where i∈KODQ1 Capacity of a first-echelon truck.Qs Capacity of a satellite.Q2 Capacity of a city freighter.Qr Capacity of a TN.H Fixed cost of employing an OD.ρOD Compensation per distance unit travelled by an OD.cd Variable cost per distance unit for either a first-echelon truck or a city freighter.μ A multiplier for calculating the extra-travelled distance of an OD.


## 4.4. Decision Variables

xijk Binary variables representing the selection of arc i,j to be traversed by first-echelon truck k∈K1; 1 if first-echelon truck k traverses through arc i,j and 0 otherwise.yijk Binary variables representing the selection of arc i,j to be traversed by city freighter k∈K2; 1 if city freighter k traverses through arc i,j and 0 otherwise.zijk Binary variables representing the selection of arc i,j to be traversed by OD k∈KOD; 1 if OD k traverses through arc i,j and 0 otherwise.αpk Binary variable indicating the assignment of OD k to pickup point p∈Np; 1 if OD k visits pickup point i and 0 otherwise.βik Binary variable representing whether customer i∈Nc is the last node visited by OD k∈KOD; 1 if OD k goes to the destination after visiting customer i and 0 otherwise.σkOD Total demands assigned to OD k∈KOD, σkOD≥0σk2 Total demands assigned to city freighter k∈K2, σk2≥0σk1 Total demands assigned to first-echelon truck k∈K1, σk1≥0liM Total demands assigned to TN i∈Nr, liM≥0liS Total demands assigned to satellite i∈NS, liS≥0qik1 Remaining loads carried by first-echelon truck k∈K1 when visiting node i∈NS, qik1≥0qik2 Remaining loads carried by city freighter k∈K2 when visiting node i∈Nv, qik2≥0qikOD Remaining loads carried by OD k∈KOD when visiting node i∈Nc, qikOD≥0


## 4.5. Objective Function

(1)MinZ=∑k∈K1∑i,j∈N1cddijxijk+∑k∈K2∑i,j∈N2cddijyijk+∑k∈KOD∑i∈nko∪N2j∈N2∪nkDρODdijzijk+∑k∈KOD∑p∈NpHαpk.Subject to(2)∑j∈N1x0jk≤1∀k∈K1,(3)∑j∈N1xj0k≤1∀k∈K1,(4)∑i∈N1xijk≤1∀j∈NS,∀k∈K1,(5)∑j∈N1xijk−∑j∈N1xjik=0i∈N1,∀k∈K1,(6)∑k∈K2∑j∈Nvypjkσk2+∑k∈KODαpkσkOD≤lpSp∈NS,(7)∑i∈N1∑j∈NSxijkljS≤σk1∀k∈K1,(8)σk1≤Q1∀k∈K1,(9)σk1≥qjk1+M1−x0jk∀j∈NS,∀k∈K1,(10)qik1−liS≤qjk1+M1−xijk∀i∈NS,∀j∈NS,∀k∈K1,(11)∑j∈Nvyijk≤∑l∈K1∑j∈N1xjili∈NS,∀k∈K2,(12)∑j∈NCzijk≤∑l∈K1∑j∈N1xjili∈NS,∀k∈KOD.The objective of 2EVRP-TN-OD in (1) is to minimize the total cost of the delivery operation. The first, second, and third terms consecutively address the travelling costs of first-echelon trucks, city freighters, and the employed ODs, while the last one represents the payments for employing these ODs. Constraints (2) to (12) focus on managing the characteristics of the first-echelon network. Constraints (2) to (4) ensure that for each first-echelon truck, every satellite can only be visited at most once. Constraint (5) represents the conservation-of-flow constraint of every first-echelon truck. Constraint (6) defines the total demands assigned to each satellite. Constraint (7) defines the total demands carried by each first-echelon truck, while constraint (8) limits the maximum number of demands that can be carried by each first-echelon truck. Constraints (9) and (10) track the remaining demands carried by each first-echelon truck while visiting a node. Constraint (11) ensures that a city freighter can originate from a satellite if and only if the satellite is served by first-echelon trucks. Similarly, constraint (12) guarantees that each OD can collect demands at a satellite if and only if the satellite is served by at least one first-echelon truck.(13)∑k∈K2∑i∈N2yijk+∑k∈KOD∑i∈N2zijk=1j∈NC,(14)αpk≤∑l∈K2∑i∈N2yiplp∈Nr,k∈KOD,(15)αpk≤∑l∈K2∑i∈N1xiplp∈NS,k∈KOD,(16)∑p∈Npαpk≤1k∈KOD,(17)∑i∈NCβik=∑p∈Npαpkk∈KOD,(18)∑i∈NCzpik=αpkp∈Np,k∈KOD,(19)∑i∈nkozipk=αpkp∈Np,k∈KOD,(20)∑j∈nkDzijk=βiki∈NC,k∈KOD,(21)∑j∈Np∪NCzjik−∑j∈NC∪nkDzijk=0πi∈NC,k∈KOD,(22)σkOD≤qjkOD+M1−zpjkj∈NC,p∈Np,k∈KOD,(23)qikOD−δi≤qjkOD+M1−zijki,j∈NC,k∈KOD,(24)∑i∈N2∑j∈NCzijkδj≤σkODk∈KOD,(25)σkOD≤QkODk∈KOD,(26)∑p∈Npdnkopαpk+∑i∈Np∪NC∑j∈NCdijzijk+∑j∈NCdjnkDβjk≤1+μΔkk∈KOD.Constraint (13) ensures that either an OD or a city freighter must serve each customer. Constraints (14) to (26) focus on managing the characteristics of ODs’ routes. In particular, constraints (14) to (15) define the assignment of each OD to a pickup point, either a satellite or a TN. Constraint (16) ensures that each OD can only visit one pickup point to take the assigned demands. Constraint (17) ensures that each OD needs to visit a customer before reaching his/her destination if the OD is assigned a demands. Constraints (18) to (21) ensure the flow-in and flow-out of routes for an OD. Constraints (22) and (23) track the remaining demands carried by an OD while visiting customer nodes. Constraint (24) calculates the total demands carried by an OD, and constraint (25) defines the maximum capacity of an OD. Lastly, constraint (26) ensures that the total distance travelled by an OD is not greater than the maximum willingness of the OD.(27)∑i∈NS∑j∈Nvyijk≤1k∈K2,(28)∑i∈NS∪Nvyijk−∑i∈NS∪Nvyjik=0j∈Nv,k∈K2,(29)∑k∈K2∑i∈N2yijk≤1j∈Nr,(30)∑k∈KODαpkσkOD≤lpMp∈Nr,(31)liM≤Qri∈Nr,(32)∑i∈N2∑j∈NCyijkδj+∑i∈N2∑j∈NryijkljM≤σk2k∈K2,(33)σk2≤Q2k∈K2,(34)σk2≤qjk2+M1−yijki∈NS,j∈Nv,k∈K2,(35)qik2−δi≤qjk2+M1−yijki∈NC,j∈Nv,k∈K2,(36)qik2−liM≤qjk2+M1−yijklimx⟶∞i∈Nr,j∈Nv,k∈K2.Constraints (27) to (36) focus on the route characteristics of the city freighters. Constraint (27) ensures that each city freighter can only originate from one satellite. Constraint (28) ensures the flow-in and flow-out of each city freighter. Constraint (29) guarantees that each TN can only be visited at most once. Constraint (30) defines the number of demands delivered to a TN by a city freighter. Constraint (31) limits the maximum allowable demands delivered to a TN. Constraint (32) defines the total demand carried by a city freighter from a satellite, while the maximum allowable demands that can be carried by a city freighter are limited by constraint (33). Moreover, constraints (34) to (36) track the remaining demands carried by a city freighter while visiting a node, either a customer node or a TN. Finally, the domains of the decision variables are given by constraints (37)–(49).(37)xijk∈0,1∀i∈N1,∀j∈N1,∀k∈K1,(38)yijk∈0,1∀i∈N2,∀j∈N2,∀k∈K2,(39)zijk∈0,1∀i∈nko∪N2,∀j∈N2∪nkD,∀k∈KOD,(40)αpk∈0,1∀p∈Np,∀k∈KOD,(41)βik∈0,1∀i∈Nc,∀k∈KOD,(42)σk1≥0∀k∈K1,(43)σk2≥0∀k∈K2,(44)σkOD≥0∀k∈KOD,(45)liM≥0∀i∈Nr,(46)liS≥0∀i∈NS,(47)qik1≥0∀i∈NS,∀k∈K1,(48)qik2≥0∀i∈Nv,∀k∈K2,(49)qikOD≥0∀i∈NC,∀k∈KOD.The formulation belongs to a mixed-integer nonlinear program due to the multiplication of two types of decision variables in several constraints, i.e., constraints (6), (7), (30), and (32).


## 5. Proposed ALNS Algorithm for 2EVRP-TN-OD

This section presents the ALNS algorithm for solving 2EVRP-TN-OD. The following sections discuss the solution representation, the procedure for generating an initial solution, and the components of the proposed ALNS.


### 5.1. Adaptive Large Neighborhood Search

Due to the complexity of VRP variants, most of the aforementioned works opted for metaheuristics to solve the problems, which are summarized in Table1. Among them, a large neighborhood search-based algorithm is one of the most effective methods for solving VRPs [12, 14, 42]. ALNS was first introduced by Ropke and Pisinger [44] to deal with the pickup and delivery problem with time windows (PDPTW). This approach searches for a new solution using various removal as well as insertion mechanisms and the computed costs for reinserting the removed nodes from the current solutions. Originally, the removal operators included Shaw Removal, Random Removal, and Worst Removal, while the insertion ones were Greedy heuristic and Regret heuristic. They are selected based on rewards using the roulette wheel principle. A new solution will be accepted if it is better than the current one.Table 1
Summary of metaheuristics and heuristics used for solving various variants of 2EVRP, LRP, and crowd-shipping (CS) problems.
ReferenceProposed problemSolution approachArchetti et al. [33]CSIterated local searchChen et al. [35](∗∗)CSTime compatibility-based heuristic and time-expanded graph-based heuristic (heuristic approach)Macrina et al. [40]CSVariable neighborhood searchSampaio et al. [36](∗∗∗)CSAdaptive large neighborhood searchVoigt and Kuhn [37]CSAdaptive large neighborhood searchYu et al. [26]LRPSimulated annealingKaraoglan et al. [28]LRPTwo-phase simulated annealingHemmelmayr et al. [24]LRP (2EVRP based)Adaptive large neighborhood searchTing and Chen [25](∗)LRPMultiple ant colony optimizationYu and Lin [29]LRPSimulated annealingSchneider and Löffler [27]LRPTree-based search algorithm (heuristic approach)Pichka et al. [30]2ELRPSimulated annealingYu et al. [16]2ELRPSimulated annealingHuang and Ardiansyah [38]2ELRP and CSTabu searchArnold and Sörensen [32](∗∗∗∗)2ELRPProgressive filtering algorithm (heuristic approach)Wang et al. [9]2EVRPVariable neighborhood search and integer programmingBelgin et al. [10]2EVRPHybrid variable neighborhood Descent and local searchZhou et al. [18]2EVRPHybrid Genetic algorithmBreunig et al. [11]2EVRPLarge neighborhood searchJie et al. [12]2EVRPColumn generation-adaptive large neighborhood searchEnthoven et al. [14]2EVRPAdaptive large neighborhood searchKitjacharoenchai et al. [15]2EVRPLarge neighborhood search and drone truck route constructionYu et al. [31]2EVRPHybrid multistart metaheuristicAnderluh et al. [13]2EVRPCombined∈-constraints and large neighborhood search, and combined large neighborhood and the heuristic Rectangle Splitting methodMühlbauer and Fontaine [19]2EVRPParallelized large neighborhood searchYu et al. [17]2EVRP and CSAdaptive large neighborhood searchVoigt et al. [43](∗∗−∗∗)2EVRP and LRPHybrid adaptive large neighborhood searchThe works of the mentioned authors also involve the following problem: (∗): FLP and MDVRP (facility location problem and multidepot vehicle routing problem). (∗∗): The authors addressed CS through MDMPMP (Multidriver, multiparcel matching problem). (∗∗∗): The authors addressed CS through the vehicle settings of PDPTW-T (pickup and delivery problem with time windows and transfers). (∗∗∗∗): STTRPSD (single truck and trailer routing problem with satellite depots). (∗∗−∗∗): MDVRP (multidepot vehicle routing problem).Akpunar and Akpinar [45] proposed a hybrid adaptive large neighborhood (H-ALNS) search algorithm, which outperforms the existing ones in terms of computational time and solution quality. This is confirmed by Voigt et al. [43], who used H-ALNS to solve various instances of 2EVRP, LRP, and multidepot VRP (MDVRP). Despite the superiority in solving 2EVRP of HALNS, our work adopts ALNS to solve 2EVRP-TN-OD.Although the proposed ALNS of our work is partly similar to those of Hemmelmayr et al. [24] and Breunig et al. [42], they did not consider the crowd-shipping scenario in their works. Our work has a similar destroy mechanism to the “satellite-related removal mechanism” of theirs, but the destroyed nodes in our ALNS can be either satellites or transshipment nodes. When a TN is removed, a set of customers related to this node will also be removed from the second echelon, which is the uniqueness of our study. Moreover, our work’s Greedy insertion mechanism is slightly different from theirs due to the extra costs for servicing the TNs of the city freighters. The related insertion is also another difference between our work and that of the mentioned authors. To sum up, ALNS with destroy and repair mechanisms for solving the 2EVRP-TN-OD distinguishes our study from the works of Hemmelmayr et al. [24] and Breunig et al. [42].


### 5.2. Solution Representation

The solution representation of the 2EVRP-TN-OD problem consists of three parts: (1) routes of first-echelon trucks, (2) routes of city freighters, and (3) routes of ODs. First, letRk1=Rk10,Rk11,…,Rk1Rk1 be the route of first-echelon truck k, consisting of depot positioned at Rk10 and Rk1Rk1, and the sequence of visited satellites. Second, let Ri2=Ri20,Ri21,…,Ri2Ri2 denote the route of city freighter i. Herein, Ri20 represents the satellite from which city freighter i originates. Note that Ri20=Ri2Ri2 because city freighter i returns to the same satellite after visiting a sequence of nodes. The visited vertices by a city freighter can be customers and/or TNs. Third, let Ωj=Ωj0,Ωj1,…,ΩjΩj be the route of OD j consisting of two parts: (1) a satellite or a TN where the OD picks up the assigned demands Ωj0, and (2) the visited customers Ωj1,…,ΩjΩj.Figure2 illustrates the solution representation of Figure 1. For example, R11 represents the route of first-echelon truck 1, visiting satellites 1 and 3 before returning to the depot. R12 represents the route of city freighter 1, which originates from satellite 1, and performs a visit to a sequence of customers. Ω1 represents the route of OD 1 that visits transshipment node 19 and serves customers 4 and 16.Figure 2
The solution representation example of 2EVRP-TN-OD.


### 5.3. Construction Heuristic

The construction heuristic generates a feasible initial solution for ALNS. It first assigns customers to available ODs using the nearest neighborhood insertion based on the capacity of ODs. Each OD is assigned to the nearest collection point, either a satellite or a TN, to collect demands. Afterward, unassigned customers are assigned to their nearest satellite. From each utilized satellite, city freighters’ routes are then built to serve these customers and utilize TNs. Finally, routes of first-echelon trucks are constructed to serve the demands assigned to each satellite.


### 5.4. Proposed Adaptive Large Neighborhood Search Heuristic

Algorithm1 describes the general structure of the ALNS heuristic developed for solving 2EVRP-TN-OD. The construction algorithm described in Section 5.2 is executed to build an initial solution. The resulting solution is stored in Fc and copied to the working solution, Fw, and the best solution, F∗. ALNS begins with removing nodes from FW with a destroy operator selected by a roulette wheel mechanism. The number of removed nodes from the solution is randomly chosen between 1,θ. To accommodate the reconfiguration of both echelons, two conditions are expressed in Lines 5–8. The main purpose of the first condition in Lines 5–6 is to remove nodes visited in the second-echelon network, while the main purpose of the second condition in Lines 7–8 is not only for dealing with nodes in the second-echelon network but also for changing the configuration of satellites.After the removal phase is conducted,ChangeODPickupPoint is performed to change pickup points visited by ODs. The repair phase then starts by applying a selected repair operator to FW. The repair phase focuses on the second-echelon network. Thus, the first-echelon routes’ optimization procedure follows to reconstruct a set of feasible routes for first-echelon trucks, as described in Section 5.3.4.A local search procedure expressed in Lines 11–18 is performed for intensification purpose and therefore increases the possibility of obtaining a better solution. WhenFW is better (in this case, has a lower objective value) than FC, we replace FC with FW, as expressed in Lines 15–18. Furthermore, if FW successfully improves the best solution so far, then F∗ will be updated, and the selection probability of destroying and repairing operators is updated using UpdateOperatorScoreλr,λi, which only considers the performance of available operators after the previous score updates. Initially, every operator has a score of λ. In particular, an amount of score τ is added to an operator’s score whenever the operator can achieve a new best solution.The acceptance ofFW with worse quality is based on threshold acceptance; i.e., still accepting FW as long as the gap between the objective value of FW and FC is within a particular threshold value. This type of acceptance is expressed using a function, i.e., acceptFw,Fc, in Line 26. The algorithm performs a restart strategy when a certain number of nonimproving iterations, β, is reached. Finally, ALNS runs until the maximum time limit, Tmax, is reached.Algorithm 1: The proposed ALNS heuristic.
Input:OD,OI,λ,p,δ,θ,β,T,L and an instance of 2EVRP-TN-ODOutput:S∗(1)
Fc⟵InitialSolution(2)
F∗⟵Fw⟵Fc(3)
T⟵T0, i⟵0, nnoupdate⟵0(4)
q⟵RandIntδ,θ(5)
whileT≤Tmaxdo(6)
ifnnoupdate<Lthen(7)
Removeandreinsertqcustomersbasedonremovalandarepairoperatorλrandλiwherer∈ODandi∈OI(8)
else(9)
Pick−uppointremovalandarepairoperatorbasedonλiisSelected(10)
nnoupdate⟵0(11)
Fw⟵RepairChangeODPickupPointDestroyFc(12)
LocalSearchFw(13)
iffLocalSearchFw>fFwthen(14)
whilefLocalSearchFw>fFwdo(15)
LocalSearchFw(16)
iffLocalSearchFw≤fFwthen(17)
Fw⟵LocalSearchFw(18)
else(19)
Fw⟵LocalSearchFw(20)
iffFw<fFcthen(21)
Fc⟵Fw(22)
iffFc<fF∗then(23)
F∗⟵Fc(24)
λr,λi⟵UpdateOperatorScoreλr,λi(25)
nnoupdate⟵0(26)
else(27)
Fc⟵acceptFw,Fc(28)
ifnnoupdate≥βthen(29)
Fc⟵F∗(30)
nnoupdate⟵nnoupdate+1(31)
returnF∗


#### 5.4.1. Destroy Operators

The developed destroy operators are mainly adopted from Ropke and Pisinger [44] and Hemmelmayr et al. [24]. In the solution, there exist two types of nodes, i.e., customer nodes and TNs, in the second-echelon network. In this section, we thus refer to a node as either a customer node or a TN. The removal operators are specifically designed to address the existence of the aforementioned types of nodes. Whenever applicable, the number of customers and/or TNs to remove, q, is randomly chosen within the range of available customers and TNs δ,θ. The description of all destroy operators is provided as follows:(1) Pickup Point Removal. A pickup point in the solution space, either a TN or a satellite, is selected, and all customers associated with the point are put on the removal list. This operator only operates whenever nnoupdate reaches the maximum number of nonimprovements’ iterations L. After the implementation of this operator, nnoupdate is set to 0.(2) Random Removal. This operator will choose q nodes randomly and exclude them from the solution pool for reinserting purposes. If the selected removed node is a TN, then all associated customers with this TN are destroyed.(3) Worst Removal. A number of q nodes with the highest normalized removal gains, which are the absolute value of the subtraction between the original solution’s cost and the cost of the solution without the selected node, are removed using this operator. This value is then divided by the arcs’ average cost where the respected node belongs for normalization purposes. By doing this, candidates who have an extreme distance from other customers can limit the rate at which they are selected over iterations.(4) Node Neighborhood Removal. This operator randomly retrieves a node as a seed, and then q−1 nodes closest to this node are removed. If a removed node is a TN, then all related customers with the TN are removed.(5) Random Route Removal. A route in the solution pool is randomly selected and excluded from the pool along with its candidates. It is strictly prohibited for destroyed nodes using this operator to be reinserted back into their most recent satellite.(6) Route Redistribution. For each satellite, we remove a route containing one or more nodes that are actually located nearer to another satellites . The distances of nodes’ pairs are multiplied with a perturbation value pϵ [0.8, 1.2] for randomization purposes.


#### 5.4.2. Repair Operators

Previously removed customer nodes are stored in the removal list and are reinserted to generate a new solution from the existing ones in this section using the following operators:(1) Greedy Insertion. Customer nodes from the removal list are retrieved based on the ascending order of their insertion costs at all available locations in the searching space first. Those nodes are then ranked and reinserted at the location with the lowest cost among the currently available routes. A new route is also considered whenever applicable. An additional inserting cost is also accounted for by the city freighter to service unopened TNs that are launched for inserting customers to be served by unemployed ODs. Two other variants of greedy insertion are also proposed: Greedy insertion with noise, where the insertion cost is modified using the noise factor p∈0.8,1.2 for randomization purposes, and Greedy insertion with a prohibition, where nodes are not allowed to be reinserted into routes that they were recently excluded from.(2) Regret Insertion. Nodes are inserted based on their regret values, which are the absolute subtraction value between the insertion position of the node’s best cost and its second-best ones among existing routes. Similar to the greedy operator, the insertion costs also consider the scenario where a new route is opened. In addition, the regret value for each removed node holds its value until the current iteration finishes.(3) Related Insertion. This operator first randomly chooses a removed customer node as a seed and selects the following inserted candidates based on the distance between the seed and the remaining nodes in the removal pool in ascending order.


#### 5.4.3. Local Search

This study implements local search (LS) for the purpose of intensification. Five operators are considered in this procedure: insertion, interroute swap, intraroute swap, 2-opt, and 2-opt∗. For each local search operator, a solution is accepted when its objective value surpasses the current one. When a better solution is achieved, it is immediately implemented by the procedure, just prior to the termination of the search process, and the operator is restarted. This loop continues until no more improvements are achieved. In this case, the next local search operator using the same mechanism is implemented. Each move is explained as follows:(1) Insertion. One node is randomly excluded from its current route and then assigned to a different place within the same route that extracted it. If a TN is excluded for wrapping itself into a new position in the routing sequence of the second-level route, then there are no changes to the customers’ sequences served using this TN.(2) Inter-Route Swap. A pair of nodes is removed from their current route, and their origin position in the route sequence is switched with the other ones within the same route. If a TN is chosen for switching its position with another node, then no changes happen to all customers in the servicing sequence from this TN.(3) Intra-Route Swap. A pair of routes is chosen, and within themselves, a node is taken out and switched positions in the route sequence with the ones that are extracted from the other route. If a TN is chosen for swapping positions with nodes from other routes, then all customers in the servicing sequence using this TN are also delivered to the new route and are still served by this TN. It must be noted that the swapping method also checks the validity of the respected constraints when swapping nodes; hence, the demands of all TNs are determined, bringing no violations to the whole solving procedure.(4) 2-Opt. This option selects 2 consecutive pairs of nodes i,i+1 and j,j+1 in a random route through the solution pool, and then each candidate in a pair swaps with the respected candidate on the other pair. The two new pairs i,j and i+1,j+1 are then formed from the previous ones on the same route. Similar to the interroute swap, customers’ orders that are fulfilled by a selected TN are still handled by this TN.(5) 2-opt∗. This option is similar to the 2-opt, but instead of doing the swapping on the same route, a pair of random routes is selected for initiating the swapping process, and at each route, there is only one consecutive pair excluded from the existing route. Just like the intraroute swap, if a TN is chosen for this swapping operator, then its validity on the demand constraints is checked to ensure there is no violation of the route capacity, and all customers whose demands are delivered from this TN are still distributed using it.


#### 5.4.4. First-Echelon Routes’ Optimization

After the repair phase of the second-echelon network, a procedure for building first-echelon trucks’ routes is implemented to complete the solution. The procedure consists of two phases: the route construction phase and the route improvement phase. The first phase constructs first-echelon trucks’ routes following the algorithm proposed by Breunig et al. [42]. The second phase applies the improvement heuristic shown in Algorithm 2 to improve the solution obtained in Phase 1.LetF1EC and F1EW be the routes of the first-echelon network in FC and FW, respectively. For every iteration of Algorithm 2, a random number r is generated. One of the three local search operators (i.e., insertion, swap, and 2-opt) is then implemented. This procedure terminates when the number of iterations reaches μ1E. In this work, μ1E is set to 500 based on a trial-and-error experiment. In this procedure, there is no involvement with the routes that cross the TNs; thus, the cases when TNs are selected are not considered.Algorithm 2: The pseudocode of the first-echelon routes’ optimization.
(1)
fori=1toμ1Edo(2)
r⟵U0,1(3)
ifr<1/3then(4)
1E_InsertionF1EW(5)
else if1/3≤r<2/3then(6)
1E_SwapF1EW(7)
else(8)
1E_2optF1EW(9)
iffF1EW<fF1ECthen(10)
F1EC⟵F1EW(11)
end for(12)
returnF1EC


### 5.5. Our Contributions

Different from the existing 2EVRP, our 2EVRP-TN-DO included the pickup point operator for the removal process, in which the selected nodes can be either TNs or satellites, while Hemmelmayr et al. [24] addressed this as a satellite removal only. Our random and node neighborhood removals also include the retrieval of TNs, which is not mentioned in the work of the scholars mentioned.In the insertion phase, we also account for the cost of opening new TNs that are served by unemployed ODs in addition to the related insertion operator. This has not been suggested by the studies of Hemmelmayr et al. [24] and Breunig et al. [42]. Instead of implementing the local search procedure on the first level of existing work, we apply it to both first- and second-level route optimizations in our work.


## 5.1. Adaptive Large Neighborhood Search

Due to the complexity of VRP variants, most of the aforementioned works opted for metaheuristics to solve the problems, which are summarized in Table1. Among them, a large neighborhood search-based algorithm is one of the most effective methods for solving VRPs [12, 14, 42]. ALNS was first introduced by Ropke and Pisinger [44] to deal with the pickup and delivery problem with time windows (PDPTW). This approach searches for a new solution using various removal as well as insertion mechanisms and the computed costs for reinserting the removed nodes from the current solutions. Originally, the removal operators included Shaw Removal, Random Removal, and Worst Removal, while the insertion ones were Greedy heuristic and Regret heuristic. They are selected based on rewards using the roulette wheel principle. A new solution will be accepted if it is better than the current one.Table 1
Summary of metaheuristics and heuristics used for solving various variants of 2EVRP, LRP, and crowd-shipping (CS) problems.
ReferenceProposed problemSolution approachArchetti et al. [33]CSIterated local searchChen et al. [35](∗∗)CSTime compatibility-based heuristic and time-expanded graph-based heuristic (heuristic approach)Macrina et al. [40]CSVariable neighborhood searchSampaio et al. [36](∗∗∗)CSAdaptive large neighborhood searchVoigt and Kuhn [37]CSAdaptive large neighborhood searchYu et al. [26]LRPSimulated annealingKaraoglan et al. [28]LRPTwo-phase simulated annealingHemmelmayr et al. [24]LRP (2EVRP based)Adaptive large neighborhood searchTing and Chen [25](∗)LRPMultiple ant colony optimizationYu and Lin [29]LRPSimulated annealingSchneider and Löffler [27]LRPTree-based search algorithm (heuristic approach)Pichka et al. [30]2ELRPSimulated annealingYu et al. [16]2ELRPSimulated annealingHuang and Ardiansyah [38]2ELRP and CSTabu searchArnold and Sörensen [32](∗∗∗∗)2ELRPProgressive filtering algorithm (heuristic approach)Wang et al. [9]2EVRPVariable neighborhood search and integer programmingBelgin et al. [10]2EVRPHybrid variable neighborhood Descent and local searchZhou et al. [18]2EVRPHybrid Genetic algorithmBreunig et al. [11]2EVRPLarge neighborhood searchJie et al. [12]2EVRPColumn generation-adaptive large neighborhood searchEnthoven et al. [14]2EVRPAdaptive large neighborhood searchKitjacharoenchai et al. [15]2EVRPLarge neighborhood search and drone truck route constructionYu et al. [31]2EVRPHybrid multistart metaheuristicAnderluh et al. [13]2EVRPCombined∈-constraints and large neighborhood search, and combined large neighborhood and the heuristic Rectangle Splitting methodMühlbauer and Fontaine [19]2EVRPParallelized large neighborhood searchYu et al. [17]2EVRP and CSAdaptive large neighborhood searchVoigt et al. [43](∗∗−∗∗)2EVRP and LRPHybrid adaptive large neighborhood searchThe works of the mentioned authors also involve the following problem: (∗): FLP and MDVRP (facility location problem and multidepot vehicle routing problem). (∗∗): The authors addressed CS through MDMPMP (Multidriver, multiparcel matching problem). (∗∗∗): The authors addressed CS through the vehicle settings of PDPTW-T (pickup and delivery problem with time windows and transfers). (∗∗∗∗): STTRPSD (single truck and trailer routing problem with satellite depots). (∗∗−∗∗): MDVRP (multidepot vehicle routing problem).Akpunar and Akpinar [45] proposed a hybrid adaptive large neighborhood (H-ALNS) search algorithm, which outperforms the existing ones in terms of computational time and solution quality. This is confirmed by Voigt et al. [43], who used H-ALNS to solve various instances of 2EVRP, LRP, and multidepot VRP (MDVRP). Despite the superiority in solving 2EVRP of HALNS, our work adopts ALNS to solve 2EVRP-TN-OD.Although the proposed ALNS of our work is partly similar to those of Hemmelmayr et al. [24] and Breunig et al. [42], they did not consider the crowd-shipping scenario in their works. Our work has a similar destroy mechanism to the “satellite-related removal mechanism” of theirs, but the destroyed nodes in our ALNS can be either satellites or transshipment nodes. When a TN is removed, a set of customers related to this node will also be removed from the second echelon, which is the uniqueness of our study. Moreover, our work’s Greedy insertion mechanism is slightly different from theirs due to the extra costs for servicing the TNs of the city freighters. The related insertion is also another difference between our work and that of the mentioned authors. To sum up, ALNS with destroy and repair mechanisms for solving the 2EVRP-TN-OD distinguishes our study from the works of Hemmelmayr et al. [24] and Breunig et al. [42].


## 5.2. Solution Representation

The solution representation of the 2EVRP-TN-OD problem consists of three parts: (1) routes of first-echelon trucks, (2) routes of city freighters, and (3) routes of ODs. First, letRk1=Rk10,Rk11,…,Rk1Rk1 be the route of first-echelon truck k, consisting of depot positioned at Rk10 and Rk1Rk1, and the sequence of visited satellites. Second, let Ri2=Ri20,Ri21,…,Ri2Ri2 denote the route of city freighter i. Herein, Ri20 represents the satellite from which city freighter i originates. Note that Ri20=Ri2Ri2 because city freighter i returns to the same satellite after visiting a sequence of nodes. The visited vertices by a city freighter can be customers and/or TNs. Third, let Ωj=Ωj0,Ωj1,…,ΩjΩj be the route of OD j consisting of two parts: (1) a satellite or a TN where the OD picks up the assigned demands Ωj0, and (2) the visited customers Ωj1,…,ΩjΩj.Figure2 illustrates the solution representation of Figure 1. For example, R11 represents the route of first-echelon truck 1, visiting satellites 1 and 3 before returning to the depot. R12 represents the route of city freighter 1, which originates from satellite 1, and performs a visit to a sequence of customers. Ω1 represents the route of OD 1 that visits transshipment node 19 and serves customers 4 and 16.Figure 2
The solution representation example of 2EVRP-TN-OD.


## 5.3. Construction Heuristic

The construction heuristic generates a feasible initial solution for ALNS. It first assigns customers to available ODs using the nearest neighborhood insertion based on the capacity of ODs. Each OD is assigned to the nearest collection point, either a satellite or a TN, to collect demands. Afterward, unassigned customers are assigned to their nearest satellite. From each utilized satellite, city freighters’ routes are then built to serve these customers and utilize TNs. Finally, routes of first-echelon trucks are constructed to serve the demands assigned to each satellite.


## 5.4. Proposed Adaptive Large Neighborhood Search Heuristic

Algorithm1 describes the general structure of the ALNS heuristic developed for solving 2EVRP-TN-OD. The construction algorithm described in Section 5.2 is executed to build an initial solution. The resulting solution is stored in Fc and copied to the working solution, Fw, and the best solution, F∗. ALNS begins with removing nodes from FW with a destroy operator selected by a roulette wheel mechanism. The number of removed nodes from the solution is randomly chosen between 1,θ. To accommodate the reconfiguration of both echelons, two conditions are expressed in Lines 5–8. The main purpose of the first condition in Lines 5–6 is to remove nodes visited in the second-echelon network, while the main purpose of the second condition in Lines 7–8 is not only for dealing with nodes in the second-echelon network but also for changing the configuration of satellites.After the removal phase is conducted,ChangeODPickupPoint is performed to change pickup points visited by ODs. The repair phase then starts by applying a selected repair operator to FW. The repair phase focuses on the second-echelon network. Thus, the first-echelon routes’ optimization procedure follows to reconstruct a set of feasible routes for first-echelon trucks, as described in Section 5.3.4.A local search procedure expressed in Lines 11–18 is performed for intensification purpose and therefore increases the possibility of obtaining a better solution. WhenFW is better (in this case, has a lower objective value) than FC, we replace FC with FW, as expressed in Lines 15–18. Furthermore, if FW successfully improves the best solution so far, then F∗ will be updated, and the selection probability of destroying and repairing operators is updated using UpdateOperatorScoreλr,λi, which only considers the performance of available operators after the previous score updates. Initially, every operator has a score of λ. In particular, an amount of score τ is added to an operator’s score whenever the operator can achieve a new best solution.The acceptance ofFW with worse quality is based on threshold acceptance; i.e., still accepting FW as long as the gap between the objective value of FW and FC is within a particular threshold value. This type of acceptance is expressed using a function, i.e., acceptFw,Fc, in Line 26. The algorithm performs a restart strategy when a certain number of nonimproving iterations, β, is reached. Finally, ALNS runs until the maximum time limit, Tmax, is reached.Algorithm 1: The proposed ALNS heuristic.
Input:OD,OI,λ,p,δ,θ,β,T,L and an instance of 2EVRP-TN-ODOutput:S∗(1)
Fc⟵InitialSolution(2)
F∗⟵Fw⟵Fc(3)
T⟵T0, i⟵0, nnoupdate⟵0(4)
q⟵RandIntδ,θ(5)
whileT≤Tmaxdo(6)
ifnnoupdate<Lthen(7)
Removeandreinsertqcustomersbasedonremovalandarepairoperatorλrandλiwherer∈ODandi∈OI(8)
else(9)
Pick−uppointremovalandarepairoperatorbasedonλiisSelected(10)
nnoupdate⟵0(11)
Fw⟵RepairChangeODPickupPointDestroyFc(12)
LocalSearchFw(13)
iffLocalSearchFw>fFwthen(14)
whilefLocalSearchFw>fFwdo(15)
LocalSearchFw(16)
iffLocalSearchFw≤fFwthen(17)
Fw⟵LocalSearchFw(18)
else(19)
Fw⟵LocalSearchFw(20)
iffFw<fFcthen(21)
Fc⟵Fw(22)
iffFc<fF∗then(23)
F∗⟵Fc(24)
λr,λi⟵UpdateOperatorScoreλr,λi(25)
nnoupdate⟵0(26)
else(27)
Fc⟵acceptFw,Fc(28)
ifnnoupdate≥βthen(29)
Fc⟵F∗(30)
nnoupdate⟵nnoupdate+1(31)
returnF∗


### 5.4.1. Destroy Operators

The developed destroy operators are mainly adopted from Ropke and Pisinger [44] and Hemmelmayr et al. [24]. In the solution, there exist two types of nodes, i.e., customer nodes and TNs, in the second-echelon network. In this section, we thus refer to a node as either a customer node or a TN. The removal operators are specifically designed to address the existence of the aforementioned types of nodes. Whenever applicable, the number of customers and/or TNs to remove, q, is randomly chosen within the range of available customers and TNs δ,θ. The description of all destroy operators is provided as follows:(1) Pickup Point Removal. A pickup point in the solution space, either a TN or a satellite, is selected, and all customers associated with the point are put on the removal list. This operator only operates whenever nnoupdate reaches the maximum number of nonimprovements’ iterations L. After the implementation of this operator, nnoupdate is set to 0.(2) Random Removal. This operator will choose q nodes randomly and exclude them from the solution pool for reinserting purposes. If the selected removed node is a TN, then all associated customers with this TN are destroyed.(3) Worst Removal. A number of q nodes with the highest normalized removal gains, which are the absolute value of the subtraction between the original solution’s cost and the cost of the solution without the selected node, are removed using this operator. This value is then divided by the arcs’ average cost where the respected node belongs for normalization purposes. By doing this, candidates who have an extreme distance from other customers can limit the rate at which they are selected over iterations.(4) Node Neighborhood Removal. This operator randomly retrieves a node as a seed, and then q−1 nodes closest to this node are removed. If a removed node is a TN, then all related customers with the TN are removed.(5) Random Route Removal. A route in the solution pool is randomly selected and excluded from the pool along with its candidates. It is strictly prohibited for destroyed nodes using this operator to be reinserted back into their most recent satellite.(6) Route Redistribution. For each satellite, we remove a route containing one or more nodes that are actually located nearer to another satellites . The distances of nodes’ pairs are multiplied with a perturbation value pϵ [0.8, 1.2] for randomization purposes.


### 5.4.2. Repair Operators

Previously removed customer nodes are stored in the removal list and are reinserted to generate a new solution from the existing ones in this section using the following operators:(1) Greedy Insertion. Customer nodes from the removal list are retrieved based on the ascending order of their insertion costs at all available locations in the searching space first. Those nodes are then ranked and reinserted at the location with the lowest cost among the currently available routes. A new route is also considered whenever applicable. An additional inserting cost is also accounted for by the city freighter to service unopened TNs that are launched for inserting customers to be served by unemployed ODs. Two other variants of greedy insertion are also proposed: Greedy insertion with noise, where the insertion cost is modified using the noise factor p∈0.8,1.2 for randomization purposes, and Greedy insertion with a prohibition, where nodes are not allowed to be reinserted into routes that they were recently excluded from.(2) Regret Insertion. Nodes are inserted based on their regret values, which are the absolute subtraction value between the insertion position of the node’s best cost and its second-best ones among existing routes. Similar to the greedy operator, the insertion costs also consider the scenario where a new route is opened. In addition, the regret value for each removed node holds its value until the current iteration finishes.(3) Related Insertion. This operator first randomly chooses a removed customer node as a seed and selects the following inserted candidates based on the distance between the seed and the remaining nodes in the removal pool in ascending order.


### 5.4.3. Local Search

This study implements local search (LS) for the purpose of intensification. Five operators are considered in this procedure: insertion, interroute swap, intraroute swap, 2-opt, and 2-opt∗. For each local search operator, a solution is accepted when its objective value surpasses the current one. When a better solution is achieved, it is immediately implemented by the procedure, just prior to the termination of the search process, and the operator is restarted. This loop continues until no more improvements are achieved. In this case, the next local search operator using the same mechanism is implemented. Each move is explained as follows:(1) Insertion. One node is randomly excluded from its current route and then assigned to a different place within the same route that extracted it. If a TN is excluded for wrapping itself into a new position in the routing sequence of the second-level route, then there are no changes to the customers’ sequences served using this TN.(2) Inter-Route Swap. A pair of nodes is removed from their current route, and their origin position in the route sequence is switched with the other ones within the same route. If a TN is chosen for switching its position with another node, then no changes happen to all customers in the servicing sequence from this TN.(3) Intra-Route Swap. A pair of routes is chosen, and within themselves, a node is taken out and switched positions in the route sequence with the ones that are extracted from the other route. If a TN is chosen for swapping positions with nodes from other routes, then all customers in the servicing sequence using this TN are also delivered to the new route and are still served by this TN. It must be noted that the swapping method also checks the validity of the respected constraints when swapping nodes; hence, the demands of all TNs are determined, bringing no violations to the whole solving procedure.(4) 2-Opt. This option selects 2 consecutive pairs of nodes i,i+1 and j,j+1 in a random route through the solution pool, and then each candidate in a pair swaps with the respected candidate on the other pair. The two new pairs i,j and i+1,j+1 are then formed from the previous ones on the same route. Similar to the interroute swap, customers’ orders that are fulfilled by a selected TN are still handled by this TN.(5) 2-opt∗. This option is similar to the 2-opt, but instead of doing the swapping on the same route, a pair of random routes is selected for initiating the swapping process, and at each route, there is only one consecutive pair excluded from the existing route. Just like the intraroute swap, if a TN is chosen for this swapping operator, then its validity on the demand constraints is checked to ensure there is no violation of the route capacity, and all customers whose demands are delivered from this TN are still distributed using it.


### 5.4.4. First-Echelon Routes’ Optimization

After the repair phase of the second-echelon network, a procedure for building first-echelon trucks’ routes is implemented to complete the solution. The procedure consists of two phases: the route construction phase and the route improvement phase. The first phase constructs first-echelon trucks’ routes following the algorithm proposed by Breunig et al. [42]. The second phase applies the improvement heuristic shown in Algorithm 2 to improve the solution obtained in Phase 1.LetF1EC and F1EW be the routes of the first-echelon network in FC and FW, respectively. For every iteration of Algorithm 2, a random number r is generated. One of the three local search operators (i.e., insertion, swap, and 2-opt) is then implemented. This procedure terminates when the number of iterations reaches μ1E. In this work, μ1E is set to 500 based on a trial-and-error experiment. In this procedure, there is no involvement with the routes that cross the TNs; thus, the cases when TNs are selected are not considered.Algorithm 2: The pseudocode of the first-echelon routes’ optimization.
(1)
fori=1toμ1Edo(2)
r⟵U0,1(3)
ifr<1/3then(4)
1E_InsertionF1EW(5)
else if1/3≤r<2/3then(6)
1E_SwapF1EW(7)
else(8)
1E_2optF1EW(9)
iffF1EW<fF1ECthen(10)
F1EC⟵F1EW(11)
end for(12)
returnF1EC


## 5.4.1. Destroy Operators

The developed destroy operators are mainly adopted from Ropke and Pisinger [44] and Hemmelmayr et al. [24]. In the solution, there exist two types of nodes, i.e., customer nodes and TNs, in the second-echelon network. In this section, we thus refer to a node as either a customer node or a TN. The removal operators are specifically designed to address the existence of the aforementioned types of nodes. Whenever applicable, the number of customers and/or TNs to remove, q, is randomly chosen within the range of available customers and TNs δ,θ. The description of all destroy operators is provided as follows:(1) Pickup Point Removal. A pickup point in the solution space, either a TN or a satellite, is selected, and all customers associated with the point are put on the removal list. This operator only operates whenever nnoupdate reaches the maximum number of nonimprovements’ iterations L. After the implementation of this operator, nnoupdate is set to 0.(2) Random Removal. This operator will choose q nodes randomly and exclude them from the solution pool for reinserting purposes. If the selected removed node is a TN, then all associated customers with this TN are destroyed.(3) Worst Removal. A number of q nodes with the highest normalized removal gains, which are the absolute value of the subtraction between the original solution’s cost and the cost of the solution without the selected node, are removed using this operator. This value is then divided by the arcs’ average cost where the respected node belongs for normalization purposes. By doing this, candidates who have an extreme distance from other customers can limit the rate at which they are selected over iterations.(4) Node Neighborhood Removal. This operator randomly retrieves a node as a seed, and then q−1 nodes closest to this node are removed. If a removed node is a TN, then all related customers with the TN are removed.(5) Random Route Removal. A route in the solution pool is randomly selected and excluded from the pool along with its candidates. It is strictly prohibited for destroyed nodes using this operator to be reinserted back into their most recent satellite.(6) Route Redistribution. For each satellite, we remove a route containing one or more nodes that are actually located nearer to another satellites . The distances of nodes’ pairs are multiplied with a perturbation value pϵ [0.8, 1.2] for randomization purposes.


## 5.4.2. Repair Operators

Previously removed customer nodes are stored in the removal list and are reinserted to generate a new solution from the existing ones in this section using the following operators:(1) Greedy Insertion. Customer nodes from the removal list are retrieved based on the ascending order of their insertion costs at all available locations in the searching space first. Those nodes are then ranked and reinserted at the location with the lowest cost among the currently available routes. A new route is also considered whenever applicable. An additional inserting cost is also accounted for by the city freighter to service unopened TNs that are launched for inserting customers to be served by unemployed ODs. Two other variants of greedy insertion are also proposed: Greedy insertion with noise, where the insertion cost is modified using the noise factor p∈0.8,1.2 for randomization purposes, and Greedy insertion with a prohibition, where nodes are not allowed to be reinserted into routes that they were recently excluded from.(2) Regret Insertion. Nodes are inserted based on their regret values, which are the absolute subtraction value between the insertion position of the node’s best cost and its second-best ones among existing routes. Similar to the greedy operator, the insertion costs also consider the scenario where a new route is opened. In addition, the regret value for each removed node holds its value until the current iteration finishes.(3) Related Insertion. This operator first randomly chooses a removed customer node as a seed and selects the following inserted candidates based on the distance between the seed and the remaining nodes in the removal pool in ascending order.


## 5.4.3. Local Search

This study implements local search (LS) for the purpose of intensification. Five operators are considered in this procedure: insertion, interroute swap, intraroute swap, 2-opt, and 2-opt∗. For each local search operator, a solution is accepted when its objective value surpasses the current one. When a better solution is achieved, it is immediately implemented by the procedure, just prior to the termination of the search process, and the operator is restarted. This loop continues until no more improvements are achieved. In this case, the next local search operator using the same mechanism is implemented. Each move is explained as follows:(1) Insertion. One node is randomly excluded from its current route and then assigned to a different place within the same route that extracted it. If a TN is excluded for wrapping itself into a new position in the routing sequence of the second-level route, then there are no changes to the customers’ sequences served using this TN.(2) Inter-Route Swap. A pair of nodes is removed from their current route, and their origin position in the route sequence is switched with the other ones within the same route. If a TN is chosen for switching its position with another node, then no changes happen to all customers in the servicing sequence from this TN.(3) Intra-Route Swap. A pair of routes is chosen, and within themselves, a node is taken out and switched positions in the route sequence with the ones that are extracted from the other route. If a TN is chosen for swapping positions with nodes from other routes, then all customers in the servicing sequence using this TN are also delivered to the new route and are still served by this TN. It must be noted that the swapping method also checks the validity of the respected constraints when swapping nodes; hence, the demands of all TNs are determined, bringing no violations to the whole solving procedure.(4) 2-Opt. This option selects 2 consecutive pairs of nodes i,i+1 and j,j+1 in a random route through the solution pool, and then each candidate in a pair swaps with the respected candidate on the other pair. The two new pairs i,j and i+1,j+1 are then formed from the previous ones on the same route. Similar to the interroute swap, customers’ orders that are fulfilled by a selected TN are still handled by this TN.(5) 2-opt∗. This option is similar to the 2-opt, but instead of doing the swapping on the same route, a pair of random routes is selected for initiating the swapping process, and at each route, there is only one consecutive pair excluded from the existing route. Just like the intraroute swap, if a TN is chosen for this swapping operator, then its validity on the demand constraints is checked to ensure there is no violation of the route capacity, and all customers whose demands are delivered from this TN are still distributed using it.


## 5.4.4. First-Echelon Routes’ Optimization

After the repair phase of the second-echelon network, a procedure for building first-echelon trucks’ routes is implemented to complete the solution. The procedure consists of two phases: the route construction phase and the route improvement phase. The first phase constructs first-echelon trucks’ routes following the algorithm proposed by Breunig et al. [42]. The second phase applies the improvement heuristic shown in Algorithm 2 to improve the solution obtained in Phase 1.LetF1EC and F1EW be the routes of the first-echelon network in FC and FW, respectively. For every iteration of Algorithm 2, a random number r is generated. One of the three local search operators (i.e., insertion, swap, and 2-opt) is then implemented. This procedure terminates when the number of iterations reaches μ1E. In this work, μ1E is set to 500 based on a trial-and-error experiment. In this procedure, there is no involvement with the routes that cross the TNs; thus, the cases when TNs are selected are not considered.Algorithm 2: The pseudocode of the first-echelon routes’ optimization.
(1)
fori=1toμ1Edo(2)
r⟵U0,1(3)
ifr<1/3then(4)
1E_InsertionF1EW(5)
else if1/3≤r<2/3then(6)
1E_SwapF1EW(7)
else(8)
1E_2optF1EW(9)
iffF1EW<fF1ECthen(10)
F1EC⟵F1EW(11)
end for(12)
returnF1EC


## 5.5. Our Contributions

Different from the existing 2EVRP, our 2EVRP-TN-DO included the pickup point operator for the removal process, in which the selected nodes can be either TNs or satellites, while Hemmelmayr et al. [24] addressed this as a satellite removal only. Our random and node neighborhood removals also include the retrieval of TNs, which is not mentioned in the work of the scholars mentioned.In the insertion phase, we also account for the cost of opening new TNs that are served by unemployed ODs in addition to the related insertion operator. This has not been suggested by the studies of Hemmelmayr et al. [24] and Breunig et al. [42]. Instead of implementing the local search procedure on the first level of existing work, we apply it to both first- and second-level route optimizations in our work.


## 6. Computational Experiments

MINLP is solved using AMPL with the Gurobi solver, and the proposed ALNS is coded in Microsoft Visual Studio C++ 2019. Both of them were run on a computer with an Intel® Core i7-6700 CPU at 3.40 GHz and 8 GB of RAM under Windows 10 Professional.


### 6.1. Benchmark Instances

The 2EVRP-TN-OD datasets are generated based on the existing 2EVRP instances. There are four sets of benchmark instances. Set 1, originally proposed by Perboli et al. [8], consists of 1 depot, 2 satellites, and 12 customers. Sets 2 and 3 are larger instances, with 1 depot, 2–4 satellites, and 21–50 customers, respectively. A set of available ODs and a set of TNs are added to the 2EVRP instances. We follow Archetti et al. [33] by setting the number of ODs in the system to be equal to the number of customers. For a given instance with locations of customers xi,yi,∀i∈Nc, the origin and destination of ODs are randomly generated in the lower left-hand corner of the square mini∈Ncxi×0.75,mini∈Ncyi×0.75 and upper right-hand corner of the square maxi∈Ncxi×1.25,maxi∈Ncyi×1.25, respectively. Second, ODs’ capacity is generated under a uniform distribution of U0.05,0.25×Q2. Third, we make sure that every generated OD is able to at least serve one customer; otherwise, another OD is generated to replace the OD who is not able to serve any customers. The parameters of ODs, partly adopted from Archetti et al. [33], are provided as follows:(1)
The variable cost of an OD is set to 0.2 times the total travelled distance, and the fixed cost of an OD is set to 5.(2)
The maximum distance deviation that is willingly traversed by an OD is 0.5 times the distance between the origin and destination nodes of the OD.TNs are generated by following these rules. First, the number of TNs in the system is set to 3. This is reasonable since the number of satellites in these instances is 2 or 4. Second, for the location of TNs, we randomly select a position in the lower left-hand corner of the customer zonemini∈Ncxi,mini∈Ncyi and the upper right-hand corner of the customer zone maxi∈Ncxi,maxi∈Ncyi and locate a TN in the selected position. After that, we ensure that every TN is not too close to each other and the satellites by defining a minimum distance requirement, γ. Here, γ is set by the following equation (50). In case a newly generated TN violates the distance requirement, another location is generated using the aforementioned mechanism. This procedure repeats until all generated TNs fulfill the minimum distance requirement.(50)γ=sqrtmaxxi−minxi2+maxyi−minyi2×0.25.


### 6.2. Parameter Tuning

We follow the parameter tuning procedure for the ALNS heuristic developed by Hemmelmayr et al. [24] and Ropke and Pisinger [44]. Initially, the parameters’ values of the proposed ALNS follow the aforementioned references. By considering a set of candidate values for each parameter, the tuning procedure is conducted by only accepting the value leading to an improvement in terms of the average objective value over all selected testing instances. The aforementioned procedure is performed sequentially over all parameters being tuned. The final selected values for all parameters are presented in Table 2.Table 2
Final parameters in the proposed ALNS algorithm.
SymbolExplanationFinal valueλInitial score of destroy and repair operators50pPerturbation factor0.2δMinimum number of customers and/or TNs to be removed1θMaximum number of customers and/or TNs to be removed0.4 |N|βNumber of iterations without improvement before returning to the best-found solution1000LNumber of iterations without improvement before implementing the satellite removal100T2E−VRPTime budget for 2EVRP60 s (small instance), 900 s (large instance)T2E−VRP−TN−ODTime budget for 2EVRP-TN-OD60 s (small instance), 120 s (medium instance), 300 s (large instance)


### 6.3. Performance of the Proposed ALNS for Solving 2EVRP Instances

Three sets of 2EVRP instances from Perboli et al. [8] and Hemmelmayr et al. [24] are used to test the effectiveness of the proposed ALNS. We compare the results obtained by ALNS with three state-of-the-art algorithms proposed by Hemmelmayr et al. [24]; Breunig et al. [42]; and Enthoven et al. [14]. Tables 3–5 present the comparative Table 6 results. In these tables, columns “Best” and “CPU (s)” represent the best value obtained and the average computational time used by each method, respectively.Table 3
Comparative results for solving 2EVRP Set 2 instances
InstanceBKSHemmelmayr et al. [24]Breunig et al. [42]Enthoven et al. [14]Proposed ALNSAvg. gap to BKS (%)Best gap to BKS (%)BestCPU (s)BestCPU (s)BestCPU (s)Avg. 5BestCPU (s)Set 2aE-n22-k4-s6-17417.07417.0737417.0760417.0775417.07417.07600.000.00E-n22-k4-s8-14384.96384.9634384.9660384.96103384.96384.96600.000.00E-n22-k4-s9-19470.6470.635470.660470.664470.60470.60600.000.00E-n22-k4-s10-14371.5371.537371.560371.575371.50371.50600.000.00E-n22-k4-s11-12427.22427.2231427.2260427.22138427.22427.22600.000.00E-n22-k4-s12-16392.78392.7836392.7860392.7867392.78392.78600.000.00E-n33-k4-s1-9730.16730.1674730.1660730.16157730.16730.16600.000.00E-n33-k4-s2-13714.63714.6364714.6360714.6383714.63714.63600.000.00E-n33-k4-s3-17707.48707.4858707.4860707.4880707.48707.48600.000.00E-n33-k4-s4-5778.74778.7477778.7460778.7479778.74778.74600.000.00E-n33-k4-s7-25756.85756.8553756.8560756.8582756.85756.85600.000.00E-n33-k4-s14-22779.05779.0585779.0560779.0582779.05779.05600.000.00Set 2bE-n51-k5-s2-4-17–46530.76530.76154530.7660530.76145530.76530.76600.000.00E-n51-k5-s2-17597.49597.49100597.4960597.49140597.49597.49600.000.00E-n51-k5-s4-46530.76530.76173530.7660530.76211534.76530.76600.750.00E-n51-k5-s6-12554.81554.81149554.8160554.81199555.03554.81600.040.00E-n51-k5-s6-12-32–37531.92531.92150531.9260531.92141537.78531.92601.100.00E-n51-k5-s11-19581.64581.64182581.6460581.64131584.06581.64600.420.00E-n51-k5-s11-19-27–47527.63527.63147527.6360527.63137533.96527.63601.200.00E-n51-k5-s27-47538.22538.22136538.2260538.22237538.22538.22600.000.00E-n51-k5-s32-37552.28552.28141552.2860552.28263552.28552.28600.000.00Set 2cE-n51-k5-s2-4-17–46601.39601.3960601.39209606.75601.39600.890.00E-n51-k5-s2-17601.39601.3960601.39141602.66601.39600.210.00E-n51-k5-s4-46702.33702.3360702.33196706.78702.33600.630.00E-n51-k5-s6-12567.42567.4260567.42129569.36567.42600.340.00E-n51-k5-s6-12-32–37567.42567.4260567.42136569.97567.42600.450.00E-n51-k5-s11-19617.42617.4260617.42137621.18617.42600.610.00E-n51-k5-s11-19-27–47530.76530.7660530.76135532.10530.76600.250.00E-n51-k5-s27-47530.76530.7660530.76134530.76530.76600.000.00E-n51-k5-s32-37752.59752.5960752.6132757.74752.59600.680.00Average578.27565.5593578.2760578.27134579.76578.2760.000.250.00Table 4
Comparative results for solving 2EVRP set 3 instances.
InstanceBKSHemmelmayr et al. [24]Breunig et al. [42]Enthoven et al. [14]Proposed ALNSAvg. Gap to BKS (%)Best gap to BKS (%)BestCPU (s)BestCPU (s)BestCPU (s)Avg. 5BestCPU (s)Set 3aE-n22-k4-s13-14526.15526.1543526.1560526.1570532.25526.15601.160.00E-n22-k4-s13-16521.09521.0944521.0960521.0964521.96521.09600.170.00E-n22-k4-s13-17496.38496.3849496.3860496.3865511.59496.39603.060.00E-n22-k4-s14-19498.8498.843498.860498.863518.63498.80603.980.00E-n22-k4-s17-19512.81512.8126512.8160512.8172526.15512.81602.600.00E-n22-k4-s19-21520.42520.4234520.4260520.4262525.93520.42601.060.00E-n33-k4-s16-22672.17672.1776672.1760672.17114677.53672.17600.800.00E-n33-k4-s16-24666.02666.0277666.0260666.02106692.45666.02603.970.00E-n33-k4-s19-26680.36680.3684680.3660680.3679681.16680.37600.120.00E-n33-k4-s22-26680.37680.3777680.3760680.37136683.34680.37600.440.00E-n33-k4-s24-28670.43670.4388670.4360670.4377679.16670.43601.300.00E-n33-k4-s25-28650.58650.5863650.5860650.58112675.43650.58603.820.00Set 3bE-n51-k5-s12-18690.59690.59147690.5960690.59125693.62690.57600.440.00E-n51-k5-s12-41683.05683.05133683.0560683.05239690.74683.05601.130.00E-n51-k5-s12-43710.41710.41217710.4160710.41137710.83710.41600.060.00E-n51-k5-s39-41728.54728.54155728.5460728.54179736.08728.54601.030.00E-n51-k5-s40-41723.75723.75154723.7560723.75129736.46723.75601.760.00E-n51-k5-s40-43752.15752.15158752.1560752.15242763.2752.15601.470.00Set 3cE-n51-k5-s13-19560.73560.7360560.73119566.35560.73601.000.00E-n51-k5-s13-42564.45564.4560564.45174568.69564.45600.750.00E-n51-k5-s13-44564.45564.4560564.45121571.14564.45601.190.00E-n51-k5-s40-42746.31746.3160746.31122749.71746.31600.460.00E-n51-k5-s41-42771.56771.5660771.56203775.65771.56600.530.00E-n51-k5-s41-44802.91802.9160802.91169807.64802.91600.590.00Average641.44632.4592641.4460641.44124649.82641.4460.001.370.00Table 5
Comparative results for solving 2EVRP Set 5 instances.
InstanceBKSHemmelmayr et al. [24]Breunig et al. [42]Enthoven et al. [14]Proposed ALNSAvg. Gap to BKS (%)Best gap to BKS (%)BestCPU (s)BestCPU (s)BestCPU (s)Avg. 5BestCPU (s)100-5-11564.461578.44291564.469001605.595541583.751569.479001.230.32100-5-1b1108.621118.954761108.629001121.76991134.031112.889002.290.38100-5-21016.321016.324321016.329001022.76471024.361022.349000.790.59100-5-2b782.25784.06415782.25900782.25345803.73794.329002.751.54100-5-31045.291046.054181045.299001046.055651047.761045.299000.240.00100-5-3b828.54828.99391828.54900828.99372829.63828.549000.130.00100-10-11124.931133.173531125.5390011256431130.611127.069000.510.19100-10-1b916.25917.05397916.25900921.29402930.69924.589001.580.91100-10-2990.58997.424061002.159001010.724541014.791009.619002.441.92100-10-2b768.61770.7340774.11900775.72384795.48779.419003.501.41100-10-31043.251047.053521048.539001053.028101057.341049.489001.350.60100-1-3b850.92862.11391854.9900859.24682874.22871.199002.742.38200-10-11556.791597.198881556.799001570.8823751593.151572.309002.341.00200-10-1b1187.621225.96921187.629001195.128271253.271238.179005.534.26200-10-21365.741385.910721365.749001389.1113641413.511392.869003.501.99200-10-2b1002.851016.1410581002.859001002.638101055.701033.499005.273.06200-10-31787.731799.859161793.999001837.6210961833.301806.159002.551.03200-10-3b1197.91203.0512171197.99001219.928251266.691241.079005.743.60Average1118.811129.35591.381120.669001131.53769.671146.781134.35900.002.471.40BKS: best-known solution from the state-of-the-art algorithms. Avg. gap to BKS:ProposedALNSAvg.5−BestSolution/BestSolution×100%. Best gap to BKS: ProposedALNSBest−BestSolution/BestSolution×100%. Bold number indicates the same solution with BKS.Table 6
Literature on 2EVRP, LRP, and CS variants.
Literature1D/MD1E/2 EEICSTN/CVIRHVProblem nameHighlightsObjective functionSolution approachAnderluh et al. [13]MD2☒☐CV☐☐2eVRPSyn (two-echelon vehicle routing problem with vehicle synchronization)The existence of a gray zone, environmental factors such as greenhouse gas (GHG), disturbance, a second-echelon depot for small electric vehicles (small EVs), synchronization, and cloned satellites for multivisit purposesRouting cost in manner of time and distance, vehicle fixed cost, GHG, disturbanceCombinedϵ-constraints and large neighborhood search, and combined large neighborhood and the heuristic Rectangle Splitting methodArchetti et al. [33]1D1☐☒—☐☒VRPOD (vehicle routing problem with occasional drivers)Introduction of occasional driversRouting cost, compensation paid to occasional driversIterated local searchArnold and Sörensen [32]1D2☐☐—☐☐2ELRP and STTRPSD (two-echelon location-routing problem and the single truck and trailer routing problem with satellite depots)Proposed the single truck and trailer routing problem with satellite depots (STTRPSD) as the 2ELRP special case with configuration closed to 2EVRP, a new solution approach, limited capacity satellitesRouting cost, depot opening cost, and vehicle fixed costsProgressive filtering Algorithm(∗)Belgin et al. [10]MD2☐☐—☐☐2EVRPSPD (two-echelon vehicle routing problem with simultaneous pickup and delivery)Customers (and satellites) with pickup and delivery demand, limited capacity satellitesRouting costHybrid variable neighborhood Descent and local searchBreunig et al. [11]1D2☐☐—☐☐E2EVRP (electric two-echelon vehicle routing problem)Limited capacity satellites, small EVs, charging stations for small EVs, and limited energy capacity for small EvsRouting costs in the manner of time and energy, vehicle fixed costsLarge neighborhood searchChen et al. [35]MD1☐☒—☐☒MDMPMP (multidriver multiparcel matching problem)Safety constraints for securing parcels, synchronizations, ridesharing, relay visits at nodes.Routing cost, travel cost of ODs, parcel transferring cost of ODs, waiting time cost of ODs, and extra cost for ODs’ detours.Time compatibility-based heuristic and time-expanded graph-based heuristic(∗)Enthoven et al. [14]1D2☐☐CV☐☐2EVRP-CO (two-echelon vehicle routing problem with covering options)Have both satellites and covering locations, delivery optionsRouting cost and connection costAdaptive large neighborhood searchHemmelmayr et al. [24]1D2☐☐—☐☐LRP (based on 2EVRP) (location-routing problem as special case of the two-echelon vehicle routing problem)—Routing cost, vehicle fixed cost, and satellites opening costAdaptive large neighborhood searchHuang and Ardiansyah [38]1D2☐☒—☐☒2ELRP-TTRP (two-echelon location-routing problem with the truck and trailer routing problem)Relay pointsRouting cost, payments for hiring CSTabu searchJie et al. [12]1D2☐☐—☐☐2E-EVRP-BSS (two-echelon capacitated electric vehicle routing problem with battery swapping stations)Multiple type of EVs, charging (battery swapping) stations for EVs, limited energy capacity for EVsRouting cost, recharging cost, handling costColumn generation-adaptive large neighborhood searchKaraoglan et al. [28]1D1☐☐—☐☐LRPSPD (location-routing problem with simultaneous pickup and delivery)—Routing cost, depot opening cost, and vehicle usage costTwo-phase simulated annealingKitjacharoenchai et al. [15]1D2☐☐—☐☐2EVRPD (two-echelon vehicle routing problem with drones)Delivery truck carrying drones for deliveries at second-echelon, customers nodes are adopted as satellitesRouting cost (denoted as truck arrival time)Large neighborhood search and drone truck route constructionMacrina et al. [40]1D1☐☒TN☐☒VRPTW-TN-OD (vehicle routing problem with time windows, transshipment nodes, and occasional drivers)Capacitated transshipment nodes, synchronizationRouting cost, compensation cost for ODs activitiesVariable neighborhood searchMühlbauer and Fontaine [19]1D2☐☐—☐☐2E-CVRPSC (two-echelon capacitated vehicle routing problem with swap containers)Swap containersRouting costParallelized large neighborhood searchPichka et al. [30]1D2☐☐—☐☐2E-OLRP (two-echelon open location-routing problem)Propose the two-echelon open location-routing problem where vehicles will not return to depot and satellites after delivering parcels, limited capacity satellitesRouting cost, vehicle cost, and satellite open costSimulated annealingSampaio et al. [36]MD1☐☒TN☒☒PDPTW-T (pickup and delivery problem with time windows and transfers)Adopt the crowd-shipping characteristics in PDPTW-TRouting costAdaptive large neighborhood searchSchneider and Löffler [27]MD1☐☐—☐☐CLRP (capacitated location-routing problem)Proposed a tree-based search algorithm that is far more effective in dealing with CLRP compared to existing studiesRouting cost, depot cost, and vehicle costTree-based search Algorithm(∗)Ting and Chen [25]MD1☐☐—☒☐CLRP (capacitated location-routing problem)—Routing cost, depot opening cost, and vehicle costMultiple ant colony optimizationVoigt and Kuhn [37]MD1☐☒TN☐☒PDPTOD (pickup and delivery problem with transshipments and occasional drivers)Occasional drivers and transshipment points on multidepot scenariosRouting cost, compensation cost on detours, and stops for ODsAdaptive large neighborhood searchVoigt et al. [43]Various experimental design☐☐—☐☐Various experimental design (two-echelon VRP (2E-VRP), location-routing problem (LRP), and multidepot VRP (MDVRP))Proposed a hybrid ALNS that has higher impacts on solving 2EVRP, LRP, and MDVRP compared to existing studiesOperation costHybrid adaptive large neighborhood searchWang et al. [9]1D2☒☐—☐☐2E-CVRP-E (two-echelon capacitated vehicle routing problem with environmental considerations)Consider environmental impact on 2E-CVRPRouting cost, fuel cost, and satellites’ handling costVariable neighborhood search and integer programmingYu and Lin [29]MD1☐☐—☐☐OLRP (open location-routing problem)—Routing cost, depot opening cost, and vehicle costSimulated annealingYu et al. [26]MD1☐☐—☐☐LRP ((capacitated) location-routing problem)—Routing cost, depot opening cost, and vehicle costSimulated annealingYu et al. [16]MD1☐☐—☐☐2E-OLRP (two-echelon open location-routing problem)Considering third-party logistics (TPL) and loading-unloading zones (LUZs)Routing cost, satellites opening cost,and vehicle costSimulated annealingYu et al. [17]1D2☐☒CV☐☒2EVRPTW-CO-OD (two-echelon vehicle routing problem with time windows, covering options, and occasional drivers)Covering locations with occasional drivers in the 2EVRP, delivery optionsRouting costs, connection costs, and compensation for ODsAdaptive large neighborhood searchYu et al. [31]1D2☐☐—☐☐2EUDP-UV (two-echelon urban delivery problem with second level unmanned vehicles)Adopt the “mothership approach” scenario; first-echelon will not deliver to customers nodes; rendezvous nodesRouting cost, number of used large vehiclesHybrid multi-start metaheuristicZhou et al. [18]MD2☐☐—☐☐MD-TEVRP-DO (multidepot two-echelon vehicle routing problem with delivery options)Delivery options, multidepot 2EVRP, limited working time for vehicleRouting cost, connection cost, handling cost, and vehicle costHybrid Genetic algorithmThis work1D2☐☒TN☐☒2EVRP-TN-OD (two-echelon vehicle routing problem with transshipment nodes and occasional drivers)Transshipment nodes in the 2EVRP framework with heterogeneous vehicles, ALNS for 2EVRP-TNRouting costs, ODs employment costAdaptive large neighborhood searchAbbreviations: 1D, one depot; MD, multidepot; 1E, one-echelon; 2E, two-echelon; EI, environmental impact; CS, crowdsourcing; TN, transshipment nodes; CV, covering locations; IR, indirect related; HV, heterogeneous vehicles.Based on Tables3–5, the proposed ALNS provides results comparable to those obtained from the state-of-the-art algorithms. In particular, ALNS achieves exactly the same best solutions for Sets 2 and 3 instances, while the average gap between the best solutions obtained by our ALNS for Set 5 and BKS is 1.40%. The average gaps to BKS are 0.25% for Set 2, 1.37% for Set 3, and 2.47% for Set 5, respectively.


### 6.4. Performance of ALNS for Solving 2EVRP-TN-OD Small Instances

Table7 shows the results obtained by MINLP and the proposed ALNS. All of the instances of Set 1 are solved to optimality by the Gurobi solver, but it takes more than 5 hours on average to solve one small-scale instance. ALNS is run 5 times to solve each instance of Set 1. Based on Table 7, our ALNS can obtain results of equal quality to those produced by Gurobi’s solver. Moreover, both the best solution and average solution values are the same, showing the robustness of the proposed ALNS for solving small-scale instances. The computational time required by ALNS to obtain those results is significantly smaller. To sum up, the proposed ALNS outperforms the Gurobi solver in terms of computational time.Table 7
Comparison between Gurobi and the proposed ALNS for solving Set 1 instances of 2EVRP-TN-OD.
InstanceGurobiProposed ALNSAvg. Gap (%)Best gap (%)CostCPU (s)Avg. 5Best 5Avg. CPU (s)ODE-n13-s2-m2-1319.748,643319.7319.73000ODE-n13-s2-m2-2291.197,692291.19291.193000ODE-n13-s2-m2-3372.7419,379372.74372.743000ODE-n13-s2-m2-4321.5811,940321.58321.583000ODE-n13-s2-m2-5324.226,099324.2324.23000ODE-n13-s2-m2-6273.833,489273.83273.833000ODE-n13-s2-m2-7315.132,445315.13315.133000ODE-n13-s2-m2-8300.545,933300.54300.543000ODE-n13-s2-m2-9303.838,101303.83303.833000ODE-n13-s2-m2-10334.2858,091334.27334.273000Average315.719,181315.7315.73000


### 6.5. Performance of ALNS for Solving 2EVRP-TN-OD Larger Instances

Tables8 and 9 show the results of solving the larger instances; i.e., Sets 2 and 3 of 2EVRP-TN-OD instances. For solving each instance of Sets 2 and 3, we utilize both Gurobi and the proposed ALNS. We set the maximum running time to 3600 seconds for Gurobi. ALNS is run5 times. We provide BKS from solving the related 2EVRP instance for the purpose of analyzing the benefit of considering ODs and TNs.Table 8
Comparison between Gurobi and the proposed ALNS for solving Set 2 instances of 2EVRP-TN-OD.
InstanceBKS of 2EVRPGurobiProposed ALNSGap+Gap∗BestAverageCPU (s)Set 2aE-n22-k4-s6-17417.07329.64328.18328.2160−21.31−0.44E-n22-k4-s8-14384.96316.96316.96316.9660−17.660.00E-n22-k4-s9-19470.6372.77351.51351.6960−25.31−5.70E-n22-k4-s10-14371.5277.09277.09277.0960−25.410.00E-n22-k4-s11-12427.22367.62350.58350.5860−17.94−4.64E-n22-k4-s12-16392.78311.88306.24306.5860−22.03−1.81E-n33-k4-s1-9730.16606.01511.93513.07120−29.89−15.52E-n33-k4-s2-13714.63763.04510.77513.83120−28.53−33.06E-n33-k4-s3-17707.48—512.48512.95120−27.56—E-n33-k4-s4-5778.74761.61555.73557.81120−28.64−27.03E-n33-k4-s7-25756.85612.65536.78536.92120−29.08−12.38E-n33-k4-s14-22779.05713.55573.87583.34120−26.34−19.58Set 2bE-n51-k5-s2-4-17-46530.76—334.54350.99300−36.97—E-n51-k5-s2-17597.49—389.54401.84300−34.80—E-n51-k5-s4-46530.76—358.97365.06300−32.37—E-n51-k5-s6-12554.81—347.34360.41300−37.39—E-n51-k5-s6-12-32-37531.92—336.03341.66300−36.83—E-n51-k5-s11-19581.64—398.27413.28300−31.53—E-n51-k5-s11-19-27-47527.63—341.57363.55300−35.26—E-n51-k5-s27-47538.22—349.85353.78300−35.00—E-n51-k5-s32-37552.28—371.96373.12300−32.65—Set 2cE-n51-k5-s2-4-17-46601.39—485.8491.75300−19.22—E-n51-k5-s2-17601.39—448.58452.99300−25.41—E-n51-k5-s4-46702.33—534.21540.66300−23.94—E-n51-k5-s6-12567.42—415.9421.48300−26.70—E-n51-k5-s6-12-32-37567.42—463.34473.26300−18.34—E-n51-k5-s11-19617.42—452.49461.79300−26.71—E-n51-k5-s11-19-27-47530.76—421.11428300−20.66—E-n51-k5-s27-47530.76—440.76444.43300−16.96—E-n51-k5-s32-37752.59—670.74683.62300−10.88—Average423.10429.02−26.71Gap+=ALNSBest−BKS of2E−VRP/BKS of2E−VRP×100%,Gap∗=ALNSBest−Gurobi/ALNSBest×100%.Table 9
Comparison between Gurobi and the proposed ALNS for solving Set 3 instances of 2EVRP-TN-OD.
InstanceBKS of 2EVRPGurobiProposed ALNSGap+Gap∗BestAverageCPU (s)Set 3aE-n22-k4-s13-14526.15556.14457.51457.6560−13.05−17.74E-n22-k4-s13-16521.09442.88430.08430.0860−17.47−2.89E-n22-k4-s13-17496.38485.54438.44438.4460−11.67−9.70E-n22-k4-s14-19498.8444.79411.3411.360−17.54−7.53E-n22-k4-s17-19512.81491.22443.2443.3760−13.57−9.78E-n22-k4-s19-21520.42490.14448.91448.9160−13.74−8.41E-n33-k4-s16-22672.17744.65469.01469.88120−30.22−37.02E-n33-k4-s16-24666.02720.02478.78487.85120−28.11−33.50E-n33-k4-s19-26680.37503.14450.44450.44120−33.79−10.47E-n33-k4-s22-26680.37—462.8463.22120−31.98—E-n33-k4-s24-28670.43—494.86495.55120−26.19—E-n33-k4-s25-28650.58726.24453.11455.65120−30.35−37.61Set 3bE-n51-k5-s12-18690.59—494.18506.34300−28.44—E-n51-k5-s12-41683.05—506.84520.78300−25.80—E-n51-k5-s12-43710.41—538.31548.97300−24.23—E-n51-k5-s39-41728.54—606.12611.77300−16.80—E-n51-k5-s40-41723.75—633.6640.05300−12.46—E-n51-k5-s40-43752.15—631.78643.05300−16.00—Set 3cE-n51-k5-s13-19560.73—359.91372.4300−35.81—E-n51-k5-s13-42564.45—421.6425.68300−25.31—E-n51-k5-s13-44564.45—461.58468.5300−18.22—E-n51-k5-s40-42746.31—594.25603.82300−20.37—E-n51-k5-s41-42771.56—682.28698.77300−11.57—E-n51-k5-s41-44802.91—664.06693.47300−17.29—Average498.08507.75−21.50Based on both Tables8 and 9, Gurobi can only provide 11 and 10 feasible solutions for Sets 2 and 3, respectively. The results show that ALNS outperforms Gurobi in terms of both solution quality and computational time. While the time limit of Gurobi is 1 hour, the time limit of ALNS is set to 60–300 seconds, depending on the size of a particular instance.The impacts of considering ODs and TNs are analyzed by comparing the BKS of 2EVRP and the best solution obtained by solving the related 2EVRP-TN-OD instance. Based on Tables8 and 9, the existence of ODs and TNs on average reduces the total operational costs up to 26.71% and 21.50% for Sets 2 and 3, respectively. The higher cost resulting from 2EVRP occurs due to the longer distances that should be travelled by city freighters in order to serve all customers. In 2EVRP-TN-OD, some customers are more beneficially served by OD.


### 6.6. Sensitivity Analyses

We further elaborate on the impact of varying several characteristics of ODs on operational costs. We also analyze the impact of considering TNs in the system. The considered characteristics of ODs are described as follows:(1)
The variable costρOD and the fixed cost H of utilizing an OD.(2)
The number of available ODs,KOD.(3)
The delivery flexibility of an OD,μ.When the crowd-shipping system is implemented, the fixed and variable costs of an OD surely influence the total operational costs. Figures3(a) and 3(b) show that changes in both fixed and variable costs lead to similar behavior. When the fixed/variable cost of an OD increases, the operational cost also increases while the number of employed ODs decreases.Figure 3
Sensitivity analyses of the cost of ODs. (a) Sensitivity analyses of the variable cost of ODs. (b) Sensitivity analyses of the fixed cost of ODs.
(a)(b)Based on the observation in Figure3, it is better to keep the fixed and variable costs as low as possible. However, from a practical perspective, fixed and variable costs will be one of the concerns when one implements such a system. Although the lower fixed and variable costs of an OD lead to lower operational costs, they also impact the number of available ODs. If the values are too low, then the number of available ODs will likely decrease.Figure4 shows the impact of varying the number of available ODs on the total operational costs. The lower the number of available ODs is, the higher the operational cost is. Based on the sensitivity analyses presented in both Figures 3 and 4, anyone who plans to implement the system needs to carefully design the fixed and variable costs of employing ODs so that the costs should be reasonable enough to attract ODs and to gain operational cost savings.Figure 4
Sensitivity analyses of the number of ODs.Figure5 shows the sensitivity analyses of varying the flexibility of ODs. When the flexibility of ODs increases, the operational cost drops while the utilization of ODs increases. Similar behaviors are exhibited by varying the number of available ODs, as shown in Figure 4. However, both the marginal reduction of operational cost and the marginal improvement of the utilization of ODs are lower when the value of flexibility of ODs reaches 0.5 and the number of available ODs equals the number of customers. These phenomena indicate that the company that plans to implement this system needs to carefully design a compensation scheme to attract a reasonable number of ODs who offer a reasonable value of flexibility in order to leverage the crowd-shipping system.Figure 5
Sensitivity analyses of the flexibility of ODs.The impact of TNs on operational costs is presented in Table10. When TNs are introduced to the delivery system, there is an average cost reduction of 1.67%.Table 10
Comparison of the average operational costs of 2EVRP-TN-OD and 2EVRP-OD solutions.
Average operational costWithout TN515.11With TN506.52Difference (%)1.67


## 6.1. Benchmark Instances

The 2EVRP-TN-OD datasets are generated based on the existing 2EVRP instances. There are four sets of benchmark instances. Set 1, originally proposed by Perboli et al. [8], consists of 1 depot, 2 satellites, and 12 customers. Sets 2 and 3 are larger instances, with 1 depot, 2–4 satellites, and 21–50 customers, respectively. A set of available ODs and a set of TNs are added to the 2EVRP instances. We follow Archetti et al. [33] by setting the number of ODs in the system to be equal to the number of customers. For a given instance with locations of customers xi,yi,∀i∈Nc, the origin and destination of ODs are randomly generated in the lower left-hand corner of the square mini∈Ncxi×0.75,mini∈Ncyi×0.75 and upper right-hand corner of the square maxi∈Ncxi×1.25,maxi∈Ncyi×1.25, respectively. Second, ODs’ capacity is generated under a uniform distribution of U0.05,0.25×Q2. Third, we make sure that every generated OD is able to at least serve one customer; otherwise, another OD is generated to replace the OD who is not able to serve any customers. The parameters of ODs, partly adopted from Archetti et al. [33], are provided as follows:(1)
The variable cost of an OD is set to 0.2 times the total travelled distance, and the fixed cost of an OD is set to 5.(2)
The maximum distance deviation that is willingly traversed by an OD is 0.5 times the distance between the origin and destination nodes of the OD.TNs are generated by following these rules. First, the number of TNs in the system is set to 3. This is reasonable since the number of satellites in these instances is 2 or 4. Second, for the location of TNs, we randomly select a position in the lower left-hand corner of the customer zonemini∈Ncxi,mini∈Ncyi and the upper right-hand corner of the customer zone maxi∈Ncxi,maxi∈Ncyi and locate a TN in the selected position. After that, we ensure that every TN is not too close to each other and the satellites by defining a minimum distance requirement, γ. Here, γ is set by the following equation (50). In case a newly generated TN violates the distance requirement, another location is generated using the aforementioned mechanism. This procedure repeats until all generated TNs fulfill the minimum distance requirement.(50)γ=sqrtmaxxi−minxi2+maxyi−minyi2×0.25.


## 6.2. Parameter Tuning

We follow the parameter tuning procedure for the ALNS heuristic developed by Hemmelmayr et al. [24] and Ropke and Pisinger [44]. Initially, the parameters’ values of the proposed ALNS follow the aforementioned references. By considering a set of candidate values for each parameter, the tuning procedure is conducted by only accepting the value leading to an improvement in terms of the average objective value over all selected testing instances. The aforementioned procedure is performed sequentially over all parameters being tuned. The final selected values for all parameters are presented in Table 2.Table 2
Final parameters in the proposed ALNS algorithm.
SymbolExplanationFinal valueλInitial score of destroy and repair operators50pPerturbation factor0.2δMinimum number of customers and/or TNs to be removed1θMaximum number of customers and/or TNs to be removed0.4 |N|βNumber of iterations without improvement before returning to the best-found solution1000LNumber of iterations without improvement before implementing the satellite removal100T2E−VRPTime budget for 2EVRP60 s (small instance), 900 s (large instance)T2E−VRP−TN−ODTime budget for 2EVRP-TN-OD60 s (small instance), 120 s (medium instance), 300 s (large instance)


## 6.3. Performance of the Proposed ALNS for Solving 2EVRP Instances

Three sets of 2EVRP instances from Perboli et al. [8] and Hemmelmayr et al. [24] are used to test the effectiveness of the proposed ALNS. We compare the results obtained by ALNS with three state-of-the-art algorithms proposed by Hemmelmayr et al. [24]; Breunig et al. [42]; and Enthoven et al. [14]. Tables 3–5 present the comparative Table 6 results. In these tables, columns “Best” and “CPU (s)” represent the best value obtained and the average computational time used by each method, respectively.Table 3
Comparative results for solving 2EVRP Set 2 instances
InstanceBKSHemmelmayr et al. [24]Breunig et al. [42]Enthoven et al. [14]Proposed ALNSAvg. gap to BKS (%)Best gap to BKS (%)BestCPU (s)BestCPU (s)BestCPU (s)Avg. 5BestCPU (s)Set 2aE-n22-k4-s6-17417.07417.0737417.0760417.0775417.07417.07600.000.00E-n22-k4-s8-14384.96384.9634384.9660384.96103384.96384.96600.000.00E-n22-k4-s9-19470.6470.635470.660470.664470.60470.60600.000.00E-n22-k4-s10-14371.5371.537371.560371.575371.50371.50600.000.00E-n22-k4-s11-12427.22427.2231427.2260427.22138427.22427.22600.000.00E-n22-k4-s12-16392.78392.7836392.7860392.7867392.78392.78600.000.00E-n33-k4-s1-9730.16730.1674730.1660730.16157730.16730.16600.000.00E-n33-k4-s2-13714.63714.6364714.6360714.6383714.63714.63600.000.00E-n33-k4-s3-17707.48707.4858707.4860707.4880707.48707.48600.000.00E-n33-k4-s4-5778.74778.7477778.7460778.7479778.74778.74600.000.00E-n33-k4-s7-25756.85756.8553756.8560756.8582756.85756.85600.000.00E-n33-k4-s14-22779.05779.0585779.0560779.0582779.05779.05600.000.00Set 2bE-n51-k5-s2-4-17–46530.76530.76154530.7660530.76145530.76530.76600.000.00E-n51-k5-s2-17597.49597.49100597.4960597.49140597.49597.49600.000.00E-n51-k5-s4-46530.76530.76173530.7660530.76211534.76530.76600.750.00E-n51-k5-s6-12554.81554.81149554.8160554.81199555.03554.81600.040.00E-n51-k5-s6-12-32–37531.92531.92150531.9260531.92141537.78531.92601.100.00E-n51-k5-s11-19581.64581.64182581.6460581.64131584.06581.64600.420.00E-n51-k5-s11-19-27–47527.63527.63147527.6360527.63137533.96527.63601.200.00E-n51-k5-s27-47538.22538.22136538.2260538.22237538.22538.22600.000.00E-n51-k5-s32-37552.28552.28141552.2860552.28263552.28552.28600.000.00Set 2cE-n51-k5-s2-4-17–46601.39601.3960601.39209606.75601.39600.890.00E-n51-k5-s2-17601.39601.3960601.39141602.66601.39600.210.00E-n51-k5-s4-46702.33702.3360702.33196706.78702.33600.630.00E-n51-k5-s6-12567.42567.4260567.42129569.36567.42600.340.00E-n51-k5-s6-12-32–37567.42567.4260567.42136569.97567.42600.450.00E-n51-k5-s11-19617.42617.4260617.42137621.18617.42600.610.00E-n51-k5-s11-19-27–47530.76530.7660530.76135532.10530.76600.250.00E-n51-k5-s27-47530.76530.7660530.76134530.76530.76600.000.00E-n51-k5-s32-37752.59752.5960752.6132757.74752.59600.680.00Average578.27565.5593578.2760578.27134579.76578.2760.000.250.00Table 4
Comparative results for solving 2EVRP set 3 instances.
InstanceBKSHemmelmayr et al. [24]Breunig et al. [42]Enthoven et al. [14]Proposed ALNSAvg. Gap to BKS (%)Best gap to BKS (%)BestCPU (s)BestCPU (s)BestCPU (s)Avg. 5BestCPU (s)Set 3aE-n22-k4-s13-14526.15526.1543526.1560526.1570532.25526.15601.160.00E-n22-k4-s13-16521.09521.0944521.0960521.0964521.96521.09600.170.00E-n22-k4-s13-17496.38496.3849496.3860496.3865511.59496.39603.060.00E-n22-k4-s14-19498.8498.843498.860498.863518.63498.80603.980.00E-n22-k4-s17-19512.81512.8126512.8160512.8172526.15512.81602.600.00E-n22-k4-s19-21520.42520.4234520.4260520.4262525.93520.42601.060.00E-n33-k4-s16-22672.17672.1776672.1760672.17114677.53672.17600.800.00E-n33-k4-s16-24666.02666.0277666.0260666.02106692.45666.02603.970.00E-n33-k4-s19-26680.36680.3684680.3660680.3679681.16680.37600.120.00E-n33-k4-s22-26680.37680.3777680.3760680.37136683.34680.37600.440.00E-n33-k4-s24-28670.43670.4388670.4360670.4377679.16670.43601.300.00E-n33-k4-s25-28650.58650.5863650.5860650.58112675.43650.58603.820.00Set 3bE-n51-k5-s12-18690.59690.59147690.5960690.59125693.62690.57600.440.00E-n51-k5-s12-41683.05683.05133683.0560683.05239690.74683.05601.130.00E-n51-k5-s12-43710.41710.41217710.4160710.41137710.83710.41600.060.00E-n51-k5-s39-41728.54728.54155728.5460728.54179736.08728.54601.030.00E-n51-k5-s40-41723.75723.75154723.7560723.75129736.46723.75601.760.00E-n51-k5-s40-43752.15752.15158752.1560752.15242763.2752.15601.470.00Set 3cE-n51-k5-s13-19560.73560.7360560.73119566.35560.73601.000.00E-n51-k5-s13-42564.45564.4560564.45174568.69564.45600.750.00E-n51-k5-s13-44564.45564.4560564.45121571.14564.45601.190.00E-n51-k5-s40-42746.31746.3160746.31122749.71746.31600.460.00E-n51-k5-s41-42771.56771.5660771.56203775.65771.56600.530.00E-n51-k5-s41-44802.91802.9160802.91169807.64802.91600.590.00Average641.44632.4592641.4460641.44124649.82641.4460.001.370.00Table 5
Comparative results for solving 2EVRP Set 5 instances.
InstanceBKSHemmelmayr et al. [24]Breunig et al. [42]Enthoven et al. [14]Proposed ALNSAvg. Gap to BKS (%)Best gap to BKS (%)BestCPU (s)BestCPU (s)BestCPU (s)Avg. 5BestCPU (s)100-5-11564.461578.44291564.469001605.595541583.751569.479001.230.32100-5-1b1108.621118.954761108.629001121.76991134.031112.889002.290.38100-5-21016.321016.324321016.329001022.76471024.361022.349000.790.59100-5-2b782.25784.06415782.25900782.25345803.73794.329002.751.54100-5-31045.291046.054181045.299001046.055651047.761045.299000.240.00100-5-3b828.54828.99391828.54900828.99372829.63828.549000.130.00100-10-11124.931133.173531125.5390011256431130.611127.069000.510.19100-10-1b916.25917.05397916.25900921.29402930.69924.589001.580.91100-10-2990.58997.424061002.159001010.724541014.791009.619002.441.92100-10-2b768.61770.7340774.11900775.72384795.48779.419003.501.41100-10-31043.251047.053521048.539001053.028101057.341049.489001.350.60100-1-3b850.92862.11391854.9900859.24682874.22871.199002.742.38200-10-11556.791597.198881556.799001570.8823751593.151572.309002.341.00200-10-1b1187.621225.96921187.629001195.128271253.271238.179005.534.26200-10-21365.741385.910721365.749001389.1113641413.511392.869003.501.99200-10-2b1002.851016.1410581002.859001002.638101055.701033.499005.273.06200-10-31787.731799.859161793.999001837.6210961833.301806.159002.551.03200-10-3b1197.91203.0512171197.99001219.928251266.691241.079005.743.60Average1118.811129.35591.381120.669001131.53769.671146.781134.35900.002.471.40BKS: best-known solution from the state-of-the-art algorithms. Avg. gap to BKS:ProposedALNSAvg.5−BestSolution/BestSolution×100%. Best gap to BKS: ProposedALNSBest−BestSolution/BestSolution×100%. Bold number indicates the same solution with BKS.Table 6
Literature on 2EVRP, LRP, and CS variants.
Literature1D/MD1E/2 EEICSTN/CVIRHVProblem nameHighlightsObjective functionSolution approachAnderluh et al. [13]MD2☒☐CV☐☐2eVRPSyn (two-echelon vehicle routing problem with vehicle synchronization)The existence of a gray zone, environmental factors such as greenhouse gas (GHG), disturbance, a second-echelon depot for small electric vehicles (small EVs), synchronization, and cloned satellites for multivisit purposesRouting cost in manner of time and distance, vehicle fixed cost, GHG, disturbanceCombinedϵ-constraints and large neighborhood search, and combined large neighborhood and the heuristic Rectangle Splitting methodArchetti et al. [33]1D1☐☒—☐☒VRPOD (vehicle routing problem with occasional drivers)Introduction of occasional driversRouting cost, compensation paid to occasional driversIterated local searchArnold and Sörensen [32]1D2☐☐—☐☐2ELRP and STTRPSD (two-echelon location-routing problem and the single truck and trailer routing problem with satellite depots)Proposed the single truck and trailer routing problem with satellite depots (STTRPSD) as the 2ELRP special case with configuration closed to 2EVRP, a new solution approach, limited capacity satellitesRouting cost, depot opening cost, and vehicle fixed costsProgressive filtering Algorithm(∗)Belgin et al. [10]MD2☐☐—☐☐2EVRPSPD (two-echelon vehicle routing problem with simultaneous pickup and delivery)Customers (and satellites) with pickup and delivery demand, limited capacity satellitesRouting costHybrid variable neighborhood Descent and local searchBreunig et al. [11]1D2☐☐—☐☐E2EVRP (electric two-echelon vehicle routing problem)Limited capacity satellites, small EVs, charging stations for small EVs, and limited energy capacity for small EvsRouting costs in the manner of time and energy, vehicle fixed costsLarge neighborhood searchChen et al. [35]MD1☐☒—☐☒MDMPMP (multidriver multiparcel matching problem)Safety constraints for securing parcels, synchronizations, ridesharing, relay visits at nodes.Routing cost, travel cost of ODs, parcel transferring cost of ODs, waiting time cost of ODs, and extra cost for ODs’ detours.Time compatibility-based heuristic and time-expanded graph-based heuristic(∗)Enthoven et al. [14]1D2☐☐CV☐☐2EVRP-CO (two-echelon vehicle routing problem with covering options)Have both satellites and covering locations, delivery optionsRouting cost and connection costAdaptive large neighborhood searchHemmelmayr et al. [24]1D2☐☐—☐☐LRP (based on 2EVRP) (location-routing problem as special case of the two-echelon vehicle routing problem)—Routing cost, vehicle fixed cost, and satellites opening costAdaptive large neighborhood searchHuang and Ardiansyah [38]1D2☐☒—☐☒2ELRP-TTRP (two-echelon location-routing problem with the truck and trailer routing problem)Relay pointsRouting cost, payments for hiring CSTabu searchJie et al. [12]1D2☐☐—☐☐2E-EVRP-BSS (two-echelon capacitated electric vehicle routing problem with battery swapping stations)Multiple type of EVs, charging (battery swapping) stations for EVs, limited energy capacity for EVsRouting cost, recharging cost, handling costColumn generation-adaptive large neighborhood searchKaraoglan et al. [28]1D1☐☐—☐☐LRPSPD (location-routing problem with simultaneous pickup and delivery)—Routing cost, depot opening cost, and vehicle usage costTwo-phase simulated annealingKitjacharoenchai et al. [15]1D2☐☐—☐☐2EVRPD (two-echelon vehicle routing problem with drones)Delivery truck carrying drones for deliveries at second-echelon, customers nodes are adopted as satellitesRouting cost (denoted as truck arrival time)Large neighborhood search and drone truck route constructionMacrina et al. [40]1D1☐☒TN☐☒VRPTW-TN-OD (vehicle routing problem with time windows, transshipment nodes, and occasional drivers)Capacitated transshipment nodes, synchronizationRouting cost, compensation cost for ODs activitiesVariable neighborhood searchMühlbauer and Fontaine [19]1D2☐☐—☐☐2E-CVRPSC (two-echelon capacitated vehicle routing problem with swap containers)Swap containersRouting costParallelized large neighborhood searchPichka et al. [30]1D2☐☐—☐☐2E-OLRP (two-echelon open location-routing problem)Propose the two-echelon open location-routing problem where vehicles will not return to depot and satellites after delivering parcels, limited capacity satellitesRouting cost, vehicle cost, and satellite open costSimulated annealingSampaio et al. [36]MD1☐☒TN☒☒PDPTW-T (pickup and delivery problem with time windows and transfers)Adopt the crowd-shipping characteristics in PDPTW-TRouting costAdaptive large neighborhood searchSchneider and Löffler [27]MD1☐☐—☐☐CLRP (capacitated location-routing problem)Proposed a tree-based search algorithm that is far more effective in dealing with CLRP compared to existing studiesRouting cost, depot cost, and vehicle costTree-based search Algorithm(∗)Ting and Chen [25]MD1☐☐—☒☐CLRP (capacitated location-routing problem)—Routing cost, depot opening cost, and vehicle costMultiple ant colony optimizationVoigt and Kuhn [37]MD1☐☒TN☐☒PDPTOD (pickup and delivery problem with transshipments and occasional drivers)Occasional drivers and transshipment points on multidepot scenariosRouting cost, compensation cost on detours, and stops for ODsAdaptive large neighborhood searchVoigt et al. [43]Various experimental design☐☐—☐☐Various experimental design (two-echelon VRP (2E-VRP), location-routing problem (LRP), and multidepot VRP (MDVRP))Proposed a hybrid ALNS that has higher impacts on solving 2EVRP, LRP, and MDVRP compared to existing studiesOperation costHybrid adaptive large neighborhood searchWang et al. [9]1D2☒☐—☐☐2E-CVRP-E (two-echelon capacitated vehicle routing problem with environmental considerations)Consider environmental impact on 2E-CVRPRouting cost, fuel cost, and satellites’ handling costVariable neighborhood search and integer programmingYu and Lin [29]MD1☐☐—☐☐OLRP (open location-routing problem)—Routing cost, depot opening cost, and vehicle costSimulated annealingYu et al. [26]MD1☐☐—☐☐LRP ((capacitated) location-routing problem)—Routing cost, depot opening cost, and vehicle costSimulated annealingYu et al. [16]MD1☐☐—☐☐2E-OLRP (two-echelon open location-routing problem)Considering third-party logistics (TPL) and loading-unloading zones (LUZs)Routing cost, satellites opening cost,and vehicle costSimulated annealingYu et al. [17]1D2☐☒CV☐☒2EVRPTW-CO-OD (two-echelon vehicle routing problem with time windows, covering options, and occasional drivers)Covering locations with occasional drivers in the 2EVRP, delivery optionsRouting costs, connection costs, and compensation for ODsAdaptive large neighborhood searchYu et al. [31]1D2☐☐—☐☐2EUDP-UV (two-echelon urban delivery problem with second level unmanned vehicles)Adopt the “mothership approach” scenario; first-echelon will not deliver to customers nodes; rendezvous nodesRouting cost, number of used large vehiclesHybrid multi-start metaheuristicZhou et al. [18]MD2☐☐—☐☐MD-TEVRP-DO (multidepot two-echelon vehicle routing problem with delivery options)Delivery options, multidepot 2EVRP, limited working time for vehicleRouting cost, connection cost, handling cost, and vehicle costHybrid Genetic algorithmThis work1D2☐☒TN☐☒2EVRP-TN-OD (two-echelon vehicle routing problem with transshipment nodes and occasional drivers)Transshipment nodes in the 2EVRP framework with heterogeneous vehicles, ALNS for 2EVRP-TNRouting costs, ODs employment costAdaptive large neighborhood searchAbbreviations: 1D, one depot; MD, multidepot; 1E, one-echelon; 2E, two-echelon; EI, environmental impact; CS, crowdsourcing; TN, transshipment nodes; CV, covering locations; IR, indirect related; HV, heterogeneous vehicles.Based on Tables3–5, the proposed ALNS provides results comparable to those obtained from the state-of-the-art algorithms. In particular, ALNS achieves exactly the same best solutions for Sets 2 and 3 instances, while the average gap between the best solutions obtained by our ALNS for Set 5 and BKS is 1.40%. The average gaps to BKS are 0.25% for Set 2, 1.37% for Set 3, and 2.47% for Set 5, respectively.


## 6.4. Performance of ALNS for Solving 2EVRP-TN-OD Small Instances

Table7 shows the results obtained by MINLP and the proposed ALNS. All of the instances of Set 1 are solved to optimality by the Gurobi solver, but it takes more than 5 hours on average to solve one small-scale instance. ALNS is run 5 times to solve each instance of Set 1. Based on Table 7, our ALNS can obtain results of equal quality to those produced by Gurobi’s solver. Moreover, both the best solution and average solution values are the same, showing the robustness of the proposed ALNS for solving small-scale instances. The computational time required by ALNS to obtain those results is significantly smaller. To sum up, the proposed ALNS outperforms the Gurobi solver in terms of computational time.Table 7
Comparison between Gurobi and the proposed ALNS for solving Set 1 instances of 2EVRP-TN-OD.
InstanceGurobiProposed ALNSAvg. Gap (%)Best gap (%)CostCPU (s)Avg. 5Best 5Avg. CPU (s)ODE-n13-s2-m2-1319.748,643319.7319.73000ODE-n13-s2-m2-2291.197,692291.19291.193000ODE-n13-s2-m2-3372.7419,379372.74372.743000ODE-n13-s2-m2-4321.5811,940321.58321.583000ODE-n13-s2-m2-5324.226,099324.2324.23000ODE-n13-s2-m2-6273.833,489273.83273.833000ODE-n13-s2-m2-7315.132,445315.13315.133000ODE-n13-s2-m2-8300.545,933300.54300.543000ODE-n13-s2-m2-9303.838,101303.83303.833000ODE-n13-s2-m2-10334.2858,091334.27334.273000Average315.719,181315.7315.73000


## 6.5. Performance of ALNS for Solving 2EVRP-TN-OD Larger Instances

Tables8 and 9 show the results of solving the larger instances; i.e., Sets 2 and 3 of 2EVRP-TN-OD instances. For solving each instance of Sets 2 and 3, we utilize both Gurobi and the proposed ALNS. We set the maximum running time to 3600 seconds for Gurobi. ALNS is run5 times. We provide BKS from solving the related 2EVRP instance for the purpose of analyzing the benefit of considering ODs and TNs.Table 8
Comparison between Gurobi and the proposed ALNS for solving Set 2 instances of 2EVRP-TN-OD.
InstanceBKS of 2EVRPGurobiProposed ALNSGap+Gap∗BestAverageCPU (s)Set 2aE-n22-k4-s6-17417.07329.64328.18328.2160−21.31−0.44E-n22-k4-s8-14384.96316.96316.96316.9660−17.660.00E-n22-k4-s9-19470.6372.77351.51351.6960−25.31−5.70E-n22-k4-s10-14371.5277.09277.09277.0960−25.410.00E-n22-k4-s11-12427.22367.62350.58350.5860−17.94−4.64E-n22-k4-s12-16392.78311.88306.24306.5860−22.03−1.81E-n33-k4-s1-9730.16606.01511.93513.07120−29.89−15.52E-n33-k4-s2-13714.63763.04510.77513.83120−28.53−33.06E-n33-k4-s3-17707.48—512.48512.95120−27.56—E-n33-k4-s4-5778.74761.61555.73557.81120−28.64−27.03E-n33-k4-s7-25756.85612.65536.78536.92120−29.08−12.38E-n33-k4-s14-22779.05713.55573.87583.34120−26.34−19.58Set 2bE-n51-k5-s2-4-17-46530.76—334.54350.99300−36.97—E-n51-k5-s2-17597.49—389.54401.84300−34.80—E-n51-k5-s4-46530.76—358.97365.06300−32.37—E-n51-k5-s6-12554.81—347.34360.41300−37.39—E-n51-k5-s6-12-32-37531.92—336.03341.66300−36.83—E-n51-k5-s11-19581.64—398.27413.28300−31.53—E-n51-k5-s11-19-27-47527.63—341.57363.55300−35.26—E-n51-k5-s27-47538.22—349.85353.78300−35.00—E-n51-k5-s32-37552.28—371.96373.12300−32.65—Set 2cE-n51-k5-s2-4-17-46601.39—485.8491.75300−19.22—E-n51-k5-s2-17601.39—448.58452.99300−25.41—E-n51-k5-s4-46702.33—534.21540.66300−23.94—E-n51-k5-s6-12567.42—415.9421.48300−26.70—E-n51-k5-s6-12-32-37567.42—463.34473.26300−18.34—E-n51-k5-s11-19617.42—452.49461.79300−26.71—E-n51-k5-s11-19-27-47530.76—421.11428300−20.66—E-n51-k5-s27-47530.76—440.76444.43300−16.96—E-n51-k5-s32-37752.59—670.74683.62300−10.88—Average423.10429.02−26.71Gap+=ALNSBest−BKS of2E−VRP/BKS of2E−VRP×100%,Gap∗=ALNSBest−Gurobi/ALNSBest×100%.Table 9
Comparison between Gurobi and the proposed ALNS for solving Set 3 instances of 2EVRP-TN-OD.
InstanceBKS of 2EVRPGurobiProposed ALNSGap+Gap∗BestAverageCPU (s)Set 3aE-n22-k4-s13-14526.15556.14457.51457.6560−13.05−17.74E-n22-k4-s13-16521.09442.88430.08430.0860−17.47−2.89E-n22-k4-s13-17496.38485.54438.44438.4460−11.67−9.70E-n22-k4-s14-19498.8444.79411.3411.360−17.54−7.53E-n22-k4-s17-19512.81491.22443.2443.3760−13.57−9.78E-n22-k4-s19-21520.42490.14448.91448.9160−13.74−8.41E-n33-k4-s16-22672.17744.65469.01469.88120−30.22−37.02E-n33-k4-s16-24666.02720.02478.78487.85120−28.11−33.50E-n33-k4-s19-26680.37503.14450.44450.44120−33.79−10.47E-n33-k4-s22-26680.37—462.8463.22120−31.98—E-n33-k4-s24-28670.43—494.86495.55120−26.19—E-n33-k4-s25-28650.58726.24453.11455.65120−30.35−37.61Set 3bE-n51-k5-s12-18690.59—494.18506.34300−28.44—E-n51-k5-s12-41683.05—506.84520.78300−25.80—E-n51-k5-s12-43710.41—538.31548.97300−24.23—E-n51-k5-s39-41728.54—606.12611.77300−16.80—E-n51-k5-s40-41723.75—633.6640.05300−12.46—E-n51-k5-s40-43752.15—631.78643.05300−16.00—Set 3cE-n51-k5-s13-19560.73—359.91372.4300−35.81—E-n51-k5-s13-42564.45—421.6425.68300−25.31—E-n51-k5-s13-44564.45—461.58468.5300−18.22—E-n51-k5-s40-42746.31—594.25603.82300−20.37—E-n51-k5-s41-42771.56—682.28698.77300−11.57—E-n51-k5-s41-44802.91—664.06693.47300−17.29—Average498.08507.75−21.50Based on both Tables8 and 9, Gurobi can only provide 11 and 10 feasible solutions for Sets 2 and 3, respectively. The results show that ALNS outperforms Gurobi in terms of both solution quality and computational time. While the time limit of Gurobi is 1 hour, the time limit of ALNS is set to 60–300 seconds, depending on the size of a particular instance.The impacts of considering ODs and TNs are analyzed by comparing the BKS of 2EVRP and the best solution obtained by solving the related 2EVRP-TN-OD instance. Based on Tables8 and 9, the existence of ODs and TNs on average reduces the total operational costs up to 26.71% and 21.50% for Sets 2 and 3, respectively. The higher cost resulting from 2EVRP occurs due to the longer distances that should be travelled by city freighters in order to serve all customers. In 2EVRP-TN-OD, some customers are more beneficially served by OD.


## 6.6. Sensitivity Analyses

We further elaborate on the impact of varying several characteristics of ODs on operational costs. We also analyze the impact of considering TNs in the system. The considered characteristics of ODs are described as follows:(1)
The variable costρOD and the fixed cost H of utilizing an OD.(2)
The number of available ODs,KOD.(3)
The delivery flexibility of an OD,μ.When the crowd-shipping system is implemented, the fixed and variable costs of an OD surely influence the total operational costs. Figures3(a) and 3(b) show that changes in both fixed and variable costs lead to similar behavior. When the fixed/variable cost of an OD increases, the operational cost also increases while the number of employed ODs decreases.Figure 3
Sensitivity analyses of the cost of ODs. (a) Sensitivity analyses of the variable cost of ODs. (b) Sensitivity analyses of the fixed cost of ODs.
(a)(b)Based on the observation in Figure3, it is better to keep the fixed and variable costs as low as possible. However, from a practical perspective, fixed and variable costs will be one of the concerns when one implements such a system. Although the lower fixed and variable costs of an OD lead to lower operational costs, they also impact the number of available ODs. If the values are too low, then the number of available ODs will likely decrease.Figure4 shows the impact of varying the number of available ODs on the total operational costs. The lower the number of available ODs is, the higher the operational cost is. Based on the sensitivity analyses presented in both Figures 3 and 4, anyone who plans to implement the system needs to carefully design the fixed and variable costs of employing ODs so that the costs should be reasonable enough to attract ODs and to gain operational cost savings.Figure 4
Sensitivity analyses of the number of ODs.Figure5 shows the sensitivity analyses of varying the flexibility of ODs. When the flexibility of ODs increases, the operational cost drops while the utilization of ODs increases. Similar behaviors are exhibited by varying the number of available ODs, as shown in Figure 4. However, both the marginal reduction of operational cost and the marginal improvement of the utilization of ODs are lower when the value of flexibility of ODs reaches 0.5 and the number of available ODs equals the number of customers. These phenomena indicate that the company that plans to implement this system needs to carefully design a compensation scheme to attract a reasonable number of ODs who offer a reasonable value of flexibility in order to leverage the crowd-shipping system.Figure 5
Sensitivity analyses of the flexibility of ODs.The impact of TNs on operational costs is presented in Table10. When TNs are introduced to the delivery system, there is an average cost reduction of 1.67%.Table 10
Comparison of the average operational costs of 2EVRP-TN-OD and 2EVRP-OD solutions.
Average operational costWithout TN515.11With TN506.52Difference (%)1.67


## 7. Conclusion

This research introduces a new variant of 2EVRP called 2EVRP-TN-OD by considering the existence of transshipment nodes (TNs) and occasional drivers (ODs). The available ODs are heterogeneous in terms of capacity. In addition, a set of TNs is introduced to leverage the ODs so that ODs may collect the assigned demands at either satellites or TNs. We developed MINLP and ALNS to deal with the problem. The experiments indicate that MINLP took a significantly large amount of computational time to solve small-scale instances; i.e., 19,181 seconds. The proposed ALNS combines the problem-specific destroy and repair operators with a local search procedure for the purpose of intensification. Our developed ALNS outperforms the performance of MINLP in terms of solution quality over a significantly shorter computational time. In addition, we benchmark our ALNS with the state-of-the-art algorithms for solving the 2EVRP benchmark instances. The results show that our ALNS provides a comparable result to BKS; i.e., 0.0% gap for both Sets 2 and 3 and 1.4% gap for Set 5.The results show that considering TNs and ODs in a two-echelon distribution system is beneficial in terms of cost savings. By utilizing the developed 2EVRP-TN-OD instances, we compare the obtained results with the BKS of the associated 2EVRP instances. The results show that considering TNs and ODs can result in cost reductions of up to 26.71% and 21.50% for Sets 2 and 3, respectively.To give a deeper understanding of how ODs affect the operational cost and decisions in a two-echelon distribution system, sensitivity analyses regarding three OD parameters are provided: (1) the cost of utilizing ODs, (2) the number of available ODs, and (3) the flexibility of ODs. As expected, the cost reduction can be further achieved by considering the TNs. The marginal improvements from adopting TNs in our research can be explained by various reasons (e.g., the generated instances for experimental designs); however, the benefit of adopting TNs has validated the correctness of our work compared to existing studies. In addition, other benefits of introducing TNs, as mentioned in Section 2.4, can be considered for future work.Several interesting extensions can be deliberated due to the limitations of this study, resulting in a range of research opportunities. They include (1) designing various compensation schemes for ODs, (2) heterogeneous occasional drivers, (3) a multiobjective problem to analyze the impact of considering ODs toward both operational costs and environment-related objectives, (4) time windows of customers and ODs, and (5) the dynamic version to approach the realistic case of a system where only partial information is provided at the beginning of a planning period.



---

*Source: 5603956-2022-12-22.xml*