# Truck and Unmanned Vehicle Routing Problem with Time Windows: A Satellite Synchronization Perspective

Authors: Hongqi Li; Jiamin Zhao; Zhuopeng Zhan
Journal: Journal of Advanced Transportation
 (2022)

---

## Abstract

We consider an important feature of satellite synchronization in the practical scenario of using unmanned vehicles (UVs) carried by trucks for “last-meter” delivery and introduce the truck and UV routing problem with time windows (TUVRP-TW) for optimizing the routes of a homogeneous fleet of truck-UV combinations. A UV that has been dispatched from its truck must be picked up by the same truck or must return by itself to the depot. Customers with time windows are classified into two types: truck-UV customers (TUCs) and UV customers (UCs). The TUCs where trucks dispatch or pick up the carried UVs are regarded as satellites. Fleet coordination and satellite synchronization are essential for modelling the TUVRP-TW. We classify satellite synchronization into inner-satellite synchronization and intersatellite synchronization. The inner-satellite synchronization generally considered in the literature focuses on synchronization operations at the same satellite. Intersatellite synchronization, which focuses on synchronization operations at various satellites, allows UVs to not return to the dispatched locations, if necessary. In the mixed-integer linear programming model of the TUVRP-TW, both binary variables for identifying the appointed satellites and continuous variables for time continuity constraints are introduced to ensure the interaction between truck routes and UV routes. A hybrid algorithm based on a greedy randomized adaptive search procedure (GRASP) and a variable neighborhood search (VNS) is provided. Based on generated instances and benchmark instances, computational experiments are conducted to evaluate the performance of the intersatellite synchronization, the performance of the developed formulation, and the applicability of the hybrid algorithm.


---

## Body

## 1. Introduction

Through developing and adopting new technologies to improve operations, logistics enterprises can offer high-quality express services to customers. Great advancements in the area of unmanned vehicle (UV) technologies are bringing many new possibilities for the utilization of UVs in the logistics sectors. A number of large enterprises are investing in the delivery modes that take advantage of the utilization of UVs. At a small scale, some enterprises have tested or attempted “last-meter” small parcel deliveries through using UVs.We strive to develop an innovative methodological approach for tackling the practical scenarios that logistics and supply-chain enterprises are presently facing in the course of using UVs that are carried by trucks for “last-meter” delivery. Many enterprises are now racing toward delivery by UVs. The project “Vans & Robots,” which was launched at Mercedes-Benz Vans in the Future Transportation unit, aims at developing the “mothership approach,” which is expected to be an ideal solution for ground transport of parcels via the cooperation between vans and carried UVs (Figure1(a)). The UVs are developed as ground-based autonomous delivery robots for parcel deliveries. Compared with a single UV for “last-meter” delivery, the “mothership approach” realizes more efficient and more flexible delivery operations: on an ideal route, the van drops off one or more UVs with parcels at each stop. Each UV rolls autonomously through the last meters of the delivery route. When the delivery has been completed, the UV finds its way back to the mothership (the van) by itself. Returned UVs can be reloaded for the next delivery missions. Although many technical issues must be overcome, Daimler is not alone in developing delivery by UV. For example, JD e-business employs UVs for campus delivery in some cities in China (Figure 1(b)).Figure 1
UV examples that are developed by companies. (a) The “mothership approach” in the project “Vans & Robots,” which was launched at Mercedes-Benz Vans in the Future Transportation unit (https://www.daimler.com/innovation/specials/future-transportation-vans/parcel-carrier-2-0.html). (b) A UV that was developed by JD e-business (http://tech.qq.com/a/20160901/039363.htm).
(a)(b)We are aware of the operational challenges in optimizing delivery routes for truck-UV combinations. In the delivery application of pairing one UV with a truck, a truck departs from the depot carrying one UV and cargoes. At a specified customer, the truck drops off the carried UV with the assigned cargoes. En route, the UV requires no intervention from the truck; however, the UV must travel along a predetermined route to return to the paired truck, which may continuously make deliveries and move to a new customer. From the two-echelon network perspective, truck routes that originate at the depot are on one echelon and UV routes that originate at customers or the depot are on the other echelon. Specified customers for UV dropping or return are regarded as satellites. At a satellite for UV dropping, the arrival time of the truck should be earlier than or equal to the dispatch time of the carried UV. At a satellite for UV return, the time difference between the arrival time of the truck and that of the UV should be considered. To address the interdependence, we introduce intersatellite synchronization, which focuses on synchronization operations at satellites.The main contributions of the paper are as follows. First, we introduce and define the truck and UV routing problem with time windows (TUVRP-TW) that caters to the delivery applicability of a fleet of truck-UV combinations with respect to customer time windows. Second, we propose a mixed-integer linear programming model, in which intersatellite synchronization constraints are specially developed. Third, we provide a hybrid algorithm based on a greedy randomized adaptive search procedure (GRASP) and a variable neighborhood search (VNS). Fourth, we experimentally evaluate the intersatellite synchronization performance, the effectiveness of the developed TUVRP-TW mathematical formulation, and the applicability of the hybrid algorithm for large-scale instances.The remainder of this paper is organized as follows. Section2 reviews the literature on routing problems that involve the satellite synchronization in specified and simplified situations. Section 3 defines the TUVRP-TW and introduces the mixed-integer linear programming model, which is expressed in the vehicle flow formulation. Section 4 presents a hybrid algorithm. Section 5 describes the generated instances, the exact and heuristic results. We present the conclusions of this work in Section 6.


## 2. Literature Review

Almost all the mathematical methods for solving routing problems involving satellite synchronization are based on studies on the two-echelon vehicle routing problem (2E-VRP), the truck and trailer routing problem (TTRP), and the two-echelon location routing problem (2E-LRP). According to Crainic and Sgalambro [1], the echelon interactions, such as the fleet coordination and satellite synchronization operations, are the main challenges in modelling the routing variants. In the 2E-VRP, the TTRP, and the 2E-LRP literature, satellite synchronization constraints require that first-echelon vehicles and second-echelon vehicles meet at the same satellite. The satellite synchronization, which is called inner-satellite synchronization in this paper, focuses on synchronization operations at the same satellite. Despite the wide application backgrounds of routing problems involving satellite synchronization, many important realistic features remain to be considered. The scenario in which an autonomous delivery robot or a drone works in collaboration with a truck to distribute parcels brings new and interesting topics of routing problems involving satellite synchronization.In the literature, many forms of the travelling salesman problem (TSP) variants (e.g., the FSTSP or the TSPD) are introduced considering truck-drone delivery modes. For an overview of these variants, we refer to Chung et al. [2]; Macrina et al. [3]; Rojas Viloria et al. [4]; and Li et al. [5]. The FSTSP or the TSPD has some basic characters, for example, one truck paired with one drone is used. The truck and drone depart from and return to the depot exactly once. The drone can launch from and land on the truck only while the truck is parked at a node. One drone route is assumed to cover exactly one customer. Murray and Chu [6] introduced the FSTSP. de Freitas and Penna [7, 8], Dell’Amico et al. [9], and Jeong et al. [10] provided various methods for the FSTSP and variants. Agatz et al. [11] introduced the TSPD. Es Yurek and Ozmutlu [12]; Ha et al. [13, 14]; and Marinelli et al. [15] introduced the TSPD variants and solution methods. Poikonen and Golden [16] and Gonzalez-R et al. [17] relaxed the constraint of one drone route covering exactly one customer. Murray and Raj [18]; Moshref-Javadi et al. [19, 20]; Salama and Srinivas [21]; and Dell’Amico et al. [22] relaxed the constraint of using one truck paired with one drone. Poikonen and Golden [23] introduced the k-multi-visit drone routing problem that considered a tandem between a truck and k drones. Each drone could launch from the truck with one or more packages to deliver to customers. There is literature that investigated the vehicle routing problem (VRP) variants with drones, called the VRPD. Wang et al. [24] introduced the VRPD. Schermer et al. [25]; Sacramento et al. [26]; and Euchi and Sadok [27] proposed suitable methods for the VRPD and variants. Di Puglia Pugliese and Guerriero [28] and Das et al. [29] introduced the VRPD with time windows. Kitjacharoenchai et al. [30] and Li et al. [31] studied VRPD variants from a two-echelon perspective.Especially, Boysen et al. [32] introduced scheduling procedures for a truck-based robot delivery, which aims at minimizing the weighted number of late customer deliveries after the announced delivery time. Customers were exclusively serviced by robots. A truck loaded the shipments for a set of customers at a central depot, and a fixed part of the truck’s loading capacity was reserved to also load autonomous robots on board. Once the truck reached a drop-off point, one or multiple robots were loaded with shipments and launched to autonomously deliver goods to customers. The truck exclusively acted as a mobile satellite, and it did not access customer locations to also service customers directly. A multi-start local search procedure was proposed. A tailor-made local search procedure was provided for determining an optimal assignment of customers to drop-off points and robot depots for a predetermined truck route. Two instances with 16 and 86 nodes were randomly generated. The computational study revealed that the truck fleet could considerably be reduced if autonomous robots supported the delivery process.According to the literature, intersatellite synchronization constraints have been imposed in specified and simplified scenarios. Methodological innovation is necessary for formulating the routing variant with intersatellite synchronization constraints to cater to more general and practical scenarios that involve, e.g., customer time windows and a fleet of truck-UV combinations. In many promotional videos and other material distributed by companies developing UV or drone delivery, a customer often must be present to receive delivery. The use of time windows and intersatellite synchronization is a key point of distinction of the introduced variant in this paper. In Table1, we compare the literature with this paper to highlight the main differences. Although the body of literature on truck-drone routing variants is presently relatively large, time windows are seldom considered by the majority of the truck-drone routing variants. The variants introduced by Di Puglia Pugliese and Guerriero [28] and Das et al. [29] involve time windows, while exactly one customer is covered by each drone route. Li et al. [31] studied a variant with time windows from a two-echelon perspective, while intersatellite synchronization is forbidden.Table 1
Studies that are related to intersatellite synchronization.
ReferencesTruckTruck and UVUV capacityTime windowObjective and modelMurray and Chu [6]; de Freitas and Penna [7, 8]; Dell’Amico et al. [9]; Agatz et al. [11]; Es Yurek and Ozmutlu [12]; Ha et al. [13, 14]; Jeong et al. [10]; and Marinelli et al. [15]11 : 1OneMin makespan, time, or costs; MILPPoikonen and Golden [16] and Gonzalez-R et al. [17]11 : 1ManyMin makespan; MIPPoikonen and Golden [23]11 :nManyMin makespan; MILPMurray and Raj [18]; Moshref-Javadi et al. [19, 20]; Salama and Srinivas [21]; and Dell’Amico et al. [22]11 :nOneMin makespan, time, or costs; MILPSacramento et al. [26] and Euchi and Sadok [27]m1 : 1OneMin time or costs; MILPDi Puglia Pugliese and Guerriero [28] and Das et al. [29]m1 : 1One✓Min distance or costs; MIPWang et al. [24] and Schermer et al. [25]m1 :nOneMin makespan; MILPKitjacharoenchai et al. [30]m1 :nManyMin makespan; MILPLi et al. [31]m1 :nMany✓Min costs; MILPThis paperm1 : 1Many✓Min costs; MILPNote. “✓” denotes that the factor is considered. “MIP” denotes mixed-integer programming. “MILP” denotes mixed-integer linear programming.


## 3. The TUVRP-TW Definition and Mathematical Formulation

### 3.1. Problem Definition

From the methodological innovation perspective, the TUVRP-TW utilizes intersatellite synchronization to realize more efficient and more flexible use of UVs, and customer time windows are respected. The intersatellite synchronization permits the meeting satellites for a truck and its UV to vary, which is expected to provide more possibilities in determining the optimized routes. From the perspective of model constructions, both binary variables for identifying the appointed satellites and continuous variables for time continuity constraints ensure the interaction between truck routes and UV routes.The TUVRP-TW network includes one depot and a specified number of customers. Customers are classified into two types: truck-UV customers (denoted as TUCs) and UV customers (denoted as UCs). Each TUC can be served by one truck-UV combination or one UV, and each UC can only be served by one UV. The customer demand with time windows is known in advance and cannot be divided. At TUCs, a UV may be dispatched from the paired truck to run by itself to visit UCs or TUCs. The TUCs where trucks drop off or pick up carried UVs are regarded as satellites. Each truck can drop off and pick up the carried UV only at TUCs. The truck can continue to visit TUCs after the carried UV has been dispatched and the truck may pick up the UV at another TUC, if possible. The UV or the paired truck that arrives first at the pick-up TUC should wait for the other vehicle.We use an illustrative example in Figure2 to show the route forms in the TUVRP-TW. Other characters of the TUVRP-TW, e.g., time windows, are omitted from Figure 2, for simplicity. Figure 2 involves a network of one depot and eleven customers. The eleven customers are classified into six TUCs and five UCs. Among the six TUCs, TUC1, TUC2, and TUC5 are used as satellites. Routes travelled by trucks or truck-UV combinations are represented by solid lines and routes travelled by UVs are represented by dotted lines. In the route travelled on by a UV, the UV (which is paired with truck C) with cargoes departs by itself from the depot, visits UC4 and UC5 in order, and returns to the depot. In the route travelled by a truck-UV combination, the truck-UV combination (of truck A and UV a) with cargoes departs from the depot to TUC1. At TUC1, UV a is dispatched to serve UC1 and UC2. After serving TUC1, truck A leaves and travels to serve TUC2. The truck and the UV meet at TUC2. The truck-UV combination leaves from TUC2, serves TUC3, and returns to the depot. In the route travelled by another truck-UV combination, the truck-UV combination (of truck B and UV b) with cargoes departs from the depot to TUC5. At TUC5, UV b is dispatched to serve UC3. Then, UV b returns by itself to the depot. After serving TUC5, truck B leaves and travels to serve TUC4 and TUC6. Truck B leaves from TUC6 and returns to the depot.Figure 2
Route examples of the TUVRP-TW.The following are some of the key components that constitute the TUVRP-TW.First, the TUVRP-TW is defined in a directed complete graph, in which the nodes represent the depot and customers, and the arcs represent the possible moves of trucks or UVs. At the depot, a homogeneous fleet of truck-UV combinations is available. The depot is the only source of cargoes. UVs can be used for direct delivery from the depot.Second, the truck-UV combination includes one truck that is carrying one UV. Each truck that is in use departs from the depot carrying one UV and cargoes for the intended customers. The UVs are autonomous vehicles for cargo delivery. A truck and its UV work in a pair mode. The pair mode means a UV that has been dispatched from its truck must be picked up by the paired truck or return by itself to the depot. A UV is not allowed to depart itself and then be picked up by its truck at a satellite, while the opposite case is allowed.Third, routes are classified into three types: (i) UV routes (UV routes, which are travelled by UVs for direct delivery, originate and terminate at the depot); (ii) truck routes (truck routes are travelled by trucks to visit TUCs); and (iii) combination routes (combination routes consist of main tours and subtours). Each combination route includes one main tour and one or several subtours and is travelled by one truck-UV combination. The main tour, which includes TUCs and the depot, is travelled by the truck-UV combination or the truck. The subtour is travelled by the dispatched UV. The origin of each subtour is a satellite, namely, a specified TUC on the main tour. The destination of each subtour is a satellite (namely, a TUC on the main tour) or the depot.Fourth, the objective of the TUVRP-TW is to minimize the integrated distance, which includes the total travelling distance of the vehicles and the transformed distance of the probable waiting times of vehicles at customers. The probable waiting time is that the vehicle waits for the time window to open at a customer.To formulate the TUVRP-TW, the following assumptions are considered. First, the UV that is carried by each truck does not occupy the truck’s capacity. Second, a running range limitation related to battery-powered UVs is not considered. Third, there are enough vehicles available at the depot.


### 3.2. Mathematical Formulation

We develop a mixed-integer linear programming model for the TUVRP-TW. The mixed-integer linear programming model, which includes the objective function and constraints, is expressed by the vehicle flow formulation.The TUVRP-TW is presented on a directed complete graphG=V,A that includes vertices and arcs. The vertex set is denoted as V=0∪Vtuc∪Vuc=0,1,2,...,ntuc,...,n, in which vertex 0 is the depot and vertex subset Vc=V\0 is the customer set. The customer set consists of the TUC subset (Vtuc=1,2,...,ntuc) and the UC subset (Vuc=ntuc+1,...,n). The origins and the destinations of routes or subtours that are travelled by UVs should be the depot or satellites. We introduce the vertex set (denoted as Vtu=0,1,2,...,ntuc) of all vertices that can act simultaneously as the origin and the destination of UV routes or subtours. The arc set is denoted as A=i,j|i,j∈V. dij denotes the travelling distance of arc (i, j) (i,j∈V). We introduce set K to denote the truck-UV combination set available at the depot initially, and K=n.The weight of the cargoes that are required by customeri (i∈Vc) is denoted as qi. ei,li denotes the time window of customer i, where ei and li represent the earliest and latest service starting times, respectively, of customer i. si is the service time at customer i.The capacities of each truck and of each UV are denoted asQt and Qu, respectively. The total amount of cargoes that is needed by the intended customers in a route that is travelled by the truck-UV combination should not exceed Qt + Qu, and the capacity of each truck-UV combination is Qtu=Qt+Qu. vt and vu denote the average velocities of a truck-UV combination and of a UV. The waiting cost of a vehicle for a unit of time at a customer is τ.Several types of binary variables and continuous variables are introduced.xkijt is a binary variable that is 1 only if truck k (k∈K) covers arc (i, j) (i,j∈Vtu). xkrciju is a binary variable that is 1 only if UV k that is carried by truck k is dispatched at vertex r (r∈Vtuc) and is picked up at vertex c by covering arc (i, j) (i,j∈V) on a subtour. dtkt and dtku denote the departure time of truck k and UV k from the depot. atkit (i∈Vtu) and atkiu (i∈V) denote the arrival time of truck k and UV k at vertex i. wtkit (i∈Vtuc) and wtkiu (i∈Vc) denote the possible waiting time of truck k and UV k at customer i. The continuous variables for time are nonnegative. In various constraints, M denotes a sufficiently large positive integer.


#### 3.2.1. Objective Function

The objective function consists of four parts: the travelling distance of trucks, the transformed distance of possible waiting times of trucks at customers, the travelling distance of UVs, and the transformed distance of possible waiting times of UVs at customers. The objective function is expressed as(1)min∑k∈K∑i∈Vtu∑j∈Vtudij⋅xkijt+∑k∈K∑i∈Vcwtkit⋅vt⋅τ+∑k∈K∑r∈Vtu∑c∈Vtu∑i∈V∑j∈Vdij⋅xkrciju+∑k∈K∑i∈Vcwtkiu⋅vu⋅τ.


#### 3.2.2. Constraints

We divide the constraints into four categories: constraints on the vehicle capacity, constraints on the time continuity, constraints on the vehicle visiting vertices, and constraints on the intersatellite synchronization.Constraint (2) restricts the capacity of each UV in UV routes for direct delivery or in subtours by UVs that have been dispatched from trucks at satellites. Constraint (3) restricts the capacity of each truck-UV combination in combination routes or truck routes.(2)∑i∈Vc,i≠r,i≠c∑j∈Vqixkrciju≤Qu,k∈K,r∈Vtu,c∈Vtu,(3)∑i∈Vtuc∑j∈Vtuqixkijt+∑r∈Vtuc∑c∈Vtu∑i∈Vc,i≠r,i≠c∑j∈Vqixkrciju≤Qtu,k∈K.Constraints on the time continuity aim to decide the arrival, leaving, and waiting time for vehicles in use. According to constraints (4) and (5), the truck arrival time at the first TUC is determined by the truck departure time from the depot and the running time on the arc that connects the depot and the first TUC. According to constraints (6) and (7), the truck arrival time at a specified TUC or at the depot is determined by the truck arrival time at previous vertex, the possible waiting time for time window opening, the service time, and the running time on the arc. Constraints (8) and (9) ensure that the truck visits TUCs by respecting the TUC time windows. According to constraints (10) and (11), the UV arrival time at the first customer is determined by the UV departure time from the depot and the running time on the arc that connects the depot and the first customer. According to constraints (12) and (13), the UV arrival time at a specified vertex is determined by the UV arrival time at previous vertex, the possible waiting time for time window opening, the service time, and the UV running time on the arc. Constraints (14) and (15) ensure that the dispatched UV visits customers by respecting the customer time windows.(4)dtkt+d0ivt+M⋅1−xk0it≥atkit,k∈K,i∈Vtuc,(5)dtkt+d0ivt−M⋅1−xk0it≤atkit,k∈K,i∈Vtuc,(6)atkit+wtkit+si+dijvt+M⋅1−xkijt≥atkjt,k∈K,i∈Vtuc,j∈Vtu,(7)atkit+wtkit+si+dijvt−M⋅1−xkijt≤atkjt,k∈K,i∈Vtuc,j∈Vtu,(8)atkit+wtkit+M⋅1−∑j∈Vtuxkijt≥ei,k∈K,i∈Vtuc,(9)atkit+wtkit−M⋅1−∑j∈Vtuxkijt≤li,k∈K,i∈Vtuc,(10)dtku+d0ivu+M⋅1−xk000iu≥atkiu,k∈K,i∈Vc,(11)dtku+d0ivu−M⋅1−xk000iu≤atkiu,k∈K,i∈Vc,(12)atkiu+wtkiu+si+dijvu+M⋅1−xkrciju≥atkju,k∈K,r∈Vtu,c∈Vtu,i∈V,i≠r,j∈V,(13)atkiu+wtkiu+si+dijvu−M⋅1−xkrciju≤atkju,k∈K,r∈Vtu,c∈Vtu,i∈V,i≠r,j∈V,(14)atkju+wtkju+M⋅1−∑i∈Vxkrciju≥ej,k∈K,r∈Vtu,c∈Vtu,j∈Vc,j≠r,j≠c,(15)atkju+wtkju−M⋅1−∑i∈Vxkrciju≤lj,k∈K,r∈Vtu,c∈Vtu,j∈Vc,j≠r,j≠c.Constraints on the vehicle visiting vertices ensure the limit on the times of vehicle visiting a node and the vehicle retention. Constraints (16)–(18) ensure that a TUC is visited exactly once by a truck or a UV when the TUC is not being used as the satellite or a TUC is visited more than once if it is being used as the satellite. Constraint (19) guarantees that each UC is visited exactly once by one UV. Constraint (20) ensures that for a truck, the number of arrivals at vertex j equals the number of departures from vertex j. Constraint (21) guarantees that for a UV, the number of arrivals at a customer vertex equals the number of departures from the customer vertex. According to constraint (22), each dispatched UV leaves once from the origin of a subtour and returns once to the destination of the subtour. Constraint (23) requires that a UV departs by itself or with the paired truck no more than once from the depot.(16)∑k∈K∑i∈Vtuxkijt+∑k∈K∑r∈Vtu∑c∈Vtu∑i∈Vxkrciju≥1,j∈Vtuc,(17)∑k∈K∑i∈Vtuxtijt≤1,j∈Vtuc,(18)∑k∈K∑i∈Vtu∑r∈Vtu∑c∈Vtuxkrciju≤1,j∈Vtuc,(19)∑k∈K∑r∈Vtu∑c∈Vtu∑i∈Vxkrciju=1,j∈Vuc,(20)∑i∈Vtuxkijt=∑a∈Vtuxkjat,k∈K,j∈Vtu,(21)∑i∈Vxkrciju=∑i∈Vxkrcjiu,k∈K,r∈Vtu,c∈Vtu,j∈V,j≠r,j≠c,(22)∑i∈Vxkrcriu=∑j∈Vxkrcjcu,k∈K,r∈Vtu,c∈Vtu,i,j≠r,i,j≠c,(23)∑j∈Vtucxk0jt+∑j∈Vcxk000ju≤1,k∈K.Constraints on the intersatellite synchronization ensure the connection between the main tour and its subtour(s). The intersatellite synchronization permits the origin and the destination of a subtour to differ. It is assumed that the truck arrives once at and leaves once from the subtour origin, and the truck arrives once and leaves once from the subtour destination. The dispatched UV leaves once from the subtour origin and arrives once at the subtour destination. After dispatching the carried UV and finishing the service at a satellite, the truck continues to visit other TUCs to complete each intended service. Only if the arrival time of the dispatched UV at a satellite is between the arrival time of the paired truck at the satellite and the departure time of the paired truck UV from the satellite can the truck pick up its UV. According to constraint (24), a TUC should be included in a main tour if the TUC is used as the subtour origin. According to constraint (25), the subtour destination should be included in a main tour. From the temporal perspective, constraints (26)–(28) ensure that a UV and the paired truck meet at the appointed time. According to constraint (26), the arrival time of a UV at the subtour destination is not later than the departure time of the paired truck from the subtour destination. If the depot is used as the subtour destination, the appointed time for a UV meeting the paired truck is not restricted. According to constraints (27) and (28), the arrival time of a dispatched UV at the first customer in a subtour is determined by the arrival time of the truck-UV combination at the satellite and the UV running time on the arc.(24)∑j∈Vc∑c∈Vtuxkrcrju≤∑i∈Vtuxkirt,k∈K,r∈Vtuc,(25)∑i∈Vc∑r∈Vtucxkrcicu≤∑i∈Vtuxkict,k∈K,c∈Vtu,(26)atkct+wtkct+sc+M⋅1−∑r∈Vtuc∑i∈Vcxkrcicu≥atkcu,k∈K,c∈Vtuc,(27)atkrt+drivu+M⋅1−xkrcriu≥atkiu,k∈K,r∈Vtuc,c∈Vtu,i∈Vc,(28)atkrt+drivu−M⋅1−xkrcriu≤atkiu,k∈K,r∈Vtuc,c∈Vtu,i∈Vc.The results of the computational experiments on small-scale instances demonstrate that the TUVRP-TW formulation can be solved directly by commercial IP solver software (CPLEX in this paper). The exact solutions of small-scale instances that are found by CPLEX, which can be referenced in calibrating the heuristic, support the validity of the TUVRP-TW formulation. Instances are seldom solved directly via CPLEX in practical applications because practical instances are of typically large scale and too much computer memory and computation time are required by computers. To find satisfactory solutions for large-scale instances, we propose a hybrid algorithm based on a GRASP and a VNS. On small-scale instances, the exact solutions that are found using CPLEX and the satisfactory solutions that are output by the hybrid algorithm are compared to evaluate the performance of the heuristic. Large-scale instances are solved using the proposed heuristic.


### 3.3. A Variant without Intersatellite Synchronization

Considering that truck-UV combinations are used for “last-meter” delivery in the “mothership approach” of truck and UV delivery modes, we introduce the TUVRP-TW. The intersatellite synchronization in the TUVRP-TW allows a truck and its UV to meet at different satellites. Due to the involvement of the intersatellite synchronization, the TUVRP-TW differs from the 2E-VRP variants in the literature. In the 2E-VRP variant literature, satellite synchronization constraints require that first-echelon vehicles and second-echelon vehicles meet at the same satellite. The inner-satellite synchronization permits synchronization operations at the same satellite.We also introduce and formulate a routing variant with inner-satellite synchronization, named the TUVRP-TW-B. In the TUVRP-TW-B, intersatellite synchronization is not permitted. The TUVRP-TW-B formulation is presented in the Appendix.


## 3.1. Problem Definition

From the methodological innovation perspective, the TUVRP-TW utilizes intersatellite synchronization to realize more efficient and more flexible use of UVs, and customer time windows are respected. The intersatellite synchronization permits the meeting satellites for a truck and its UV to vary, which is expected to provide more possibilities in determining the optimized routes. From the perspective of model constructions, both binary variables for identifying the appointed satellites and continuous variables for time continuity constraints ensure the interaction between truck routes and UV routes.The TUVRP-TW network includes one depot and a specified number of customers. Customers are classified into two types: truck-UV customers (denoted as TUCs) and UV customers (denoted as UCs). Each TUC can be served by one truck-UV combination or one UV, and each UC can only be served by one UV. The customer demand with time windows is known in advance and cannot be divided. At TUCs, a UV may be dispatched from the paired truck to run by itself to visit UCs or TUCs. The TUCs where trucks drop off or pick up carried UVs are regarded as satellites. Each truck can drop off and pick up the carried UV only at TUCs. The truck can continue to visit TUCs after the carried UV has been dispatched and the truck may pick up the UV at another TUC, if possible. The UV or the paired truck that arrives first at the pick-up TUC should wait for the other vehicle.We use an illustrative example in Figure2 to show the route forms in the TUVRP-TW. Other characters of the TUVRP-TW, e.g., time windows, are omitted from Figure 2, for simplicity. Figure 2 involves a network of one depot and eleven customers. The eleven customers are classified into six TUCs and five UCs. Among the six TUCs, TUC1, TUC2, and TUC5 are used as satellites. Routes travelled by trucks or truck-UV combinations are represented by solid lines and routes travelled by UVs are represented by dotted lines. In the route travelled on by a UV, the UV (which is paired with truck C) with cargoes departs by itself from the depot, visits UC4 and UC5 in order, and returns to the depot. In the route travelled by a truck-UV combination, the truck-UV combination (of truck A and UV a) with cargoes departs from the depot to TUC1. At TUC1, UV a is dispatched to serve UC1 and UC2. After serving TUC1, truck A leaves and travels to serve TUC2. The truck and the UV meet at TUC2. The truck-UV combination leaves from TUC2, serves TUC3, and returns to the depot. In the route travelled by another truck-UV combination, the truck-UV combination (of truck B and UV b) with cargoes departs from the depot to TUC5. At TUC5, UV b is dispatched to serve UC3. Then, UV b returns by itself to the depot. After serving TUC5, truck B leaves and travels to serve TUC4 and TUC6. Truck B leaves from TUC6 and returns to the depot.Figure 2
Route examples of the TUVRP-TW.The following are some of the key components that constitute the TUVRP-TW.First, the TUVRP-TW is defined in a directed complete graph, in which the nodes represent the depot and customers, and the arcs represent the possible moves of trucks or UVs. At the depot, a homogeneous fleet of truck-UV combinations is available. The depot is the only source of cargoes. UVs can be used for direct delivery from the depot.Second, the truck-UV combination includes one truck that is carrying one UV. Each truck that is in use departs from the depot carrying one UV and cargoes for the intended customers. The UVs are autonomous vehicles for cargo delivery. A truck and its UV work in a pair mode. The pair mode means a UV that has been dispatched from its truck must be picked up by the paired truck or return by itself to the depot. A UV is not allowed to depart itself and then be picked up by its truck at a satellite, while the opposite case is allowed.Third, routes are classified into three types: (i) UV routes (UV routes, which are travelled by UVs for direct delivery, originate and terminate at the depot); (ii) truck routes (truck routes are travelled by trucks to visit TUCs); and (iii) combination routes (combination routes consist of main tours and subtours). Each combination route includes one main tour and one or several subtours and is travelled by one truck-UV combination. The main tour, which includes TUCs and the depot, is travelled by the truck-UV combination or the truck. The subtour is travelled by the dispatched UV. The origin of each subtour is a satellite, namely, a specified TUC on the main tour. The destination of each subtour is a satellite (namely, a TUC on the main tour) or the depot.Fourth, the objective of the TUVRP-TW is to minimize the integrated distance, which includes the total travelling distance of the vehicles and the transformed distance of the probable waiting times of vehicles at customers. The probable waiting time is that the vehicle waits for the time window to open at a customer.To formulate the TUVRP-TW, the following assumptions are considered. First, the UV that is carried by each truck does not occupy the truck’s capacity. Second, a running range limitation related to battery-powered UVs is not considered. Third, there are enough vehicles available at the depot.


## 3.2. Mathematical Formulation

We develop a mixed-integer linear programming model for the TUVRP-TW. The mixed-integer linear programming model, which includes the objective function and constraints, is expressed by the vehicle flow formulation.The TUVRP-TW is presented on a directed complete graphG=V,A that includes vertices and arcs. The vertex set is denoted as V=0∪Vtuc∪Vuc=0,1,2,...,ntuc,...,n, in which vertex 0 is the depot and vertex subset Vc=V\0 is the customer set. The customer set consists of the TUC subset (Vtuc=1,2,...,ntuc) and the UC subset (Vuc=ntuc+1,...,n). The origins and the destinations of routes or subtours that are travelled by UVs should be the depot or satellites. We introduce the vertex set (denoted as Vtu=0,1,2,...,ntuc) of all vertices that can act simultaneously as the origin and the destination of UV routes or subtours. The arc set is denoted as A=i,j|i,j∈V. dij denotes the travelling distance of arc (i, j) (i,j∈V). We introduce set K to denote the truck-UV combination set available at the depot initially, and K=n.The weight of the cargoes that are required by customeri (i∈Vc) is denoted as qi. ei,li denotes the time window of customer i, where ei and li represent the earliest and latest service starting times, respectively, of customer i. si is the service time at customer i.The capacities of each truck and of each UV are denoted asQt and Qu, respectively. The total amount of cargoes that is needed by the intended customers in a route that is travelled by the truck-UV combination should not exceed Qt + Qu, and the capacity of each truck-UV combination is Qtu=Qt+Qu. vt and vu denote the average velocities of a truck-UV combination and of a UV. The waiting cost of a vehicle for a unit of time at a customer is τ.Several types of binary variables and continuous variables are introduced.xkijt is a binary variable that is 1 only if truck k (k∈K) covers arc (i, j) (i,j∈Vtu). xkrciju is a binary variable that is 1 only if UV k that is carried by truck k is dispatched at vertex r (r∈Vtuc) and is picked up at vertex c by covering arc (i, j) (i,j∈V) on a subtour. dtkt and dtku denote the departure time of truck k and UV k from the depot. atkit (i∈Vtu) and atkiu (i∈V) denote the arrival time of truck k and UV k at vertex i. wtkit (i∈Vtuc) and wtkiu (i∈Vc) denote the possible waiting time of truck k and UV k at customer i. The continuous variables for time are nonnegative. In various constraints, M denotes a sufficiently large positive integer.


### 3.2.1. Objective Function

The objective function consists of four parts: the travelling distance of trucks, the transformed distance of possible waiting times of trucks at customers, the travelling distance of UVs, and the transformed distance of possible waiting times of UVs at customers. The objective function is expressed as(1)min∑k∈K∑i∈Vtu∑j∈Vtudij⋅xkijt+∑k∈K∑i∈Vcwtkit⋅vt⋅τ+∑k∈K∑r∈Vtu∑c∈Vtu∑i∈V∑j∈Vdij⋅xkrciju+∑k∈K∑i∈Vcwtkiu⋅vu⋅τ.


### 3.2.2. Constraints

We divide the constraints into four categories: constraints on the vehicle capacity, constraints on the time continuity, constraints on the vehicle visiting vertices, and constraints on the intersatellite synchronization.Constraint (2) restricts the capacity of each UV in UV routes for direct delivery or in subtours by UVs that have been dispatched from trucks at satellites. Constraint (3) restricts the capacity of each truck-UV combination in combination routes or truck routes.(2)∑i∈Vc,i≠r,i≠c∑j∈Vqixkrciju≤Qu,k∈K,r∈Vtu,c∈Vtu,(3)∑i∈Vtuc∑j∈Vtuqixkijt+∑r∈Vtuc∑c∈Vtu∑i∈Vc,i≠r,i≠c∑j∈Vqixkrciju≤Qtu,k∈K.Constraints on the time continuity aim to decide the arrival, leaving, and waiting time for vehicles in use. According to constraints (4) and (5), the truck arrival time at the first TUC is determined by the truck departure time from the depot and the running time on the arc that connects the depot and the first TUC. According to constraints (6) and (7), the truck arrival time at a specified TUC or at the depot is determined by the truck arrival time at previous vertex, the possible waiting time for time window opening, the service time, and the running time on the arc. Constraints (8) and (9) ensure that the truck visits TUCs by respecting the TUC time windows. According to constraints (10) and (11), the UV arrival time at the first customer is determined by the UV departure time from the depot and the running time on the arc that connects the depot and the first customer. According to constraints (12) and (13), the UV arrival time at a specified vertex is determined by the UV arrival time at previous vertex, the possible waiting time for time window opening, the service time, and the UV running time on the arc. Constraints (14) and (15) ensure that the dispatched UV visits customers by respecting the customer time windows.(4)dtkt+d0ivt+M⋅1−xk0it≥atkit,k∈K,i∈Vtuc,(5)dtkt+d0ivt−M⋅1−xk0it≤atkit,k∈K,i∈Vtuc,(6)atkit+wtkit+si+dijvt+M⋅1−xkijt≥atkjt,k∈K,i∈Vtuc,j∈Vtu,(7)atkit+wtkit+si+dijvt−M⋅1−xkijt≤atkjt,k∈K,i∈Vtuc,j∈Vtu,(8)atkit+wtkit+M⋅1−∑j∈Vtuxkijt≥ei,k∈K,i∈Vtuc,(9)atkit+wtkit−M⋅1−∑j∈Vtuxkijt≤li,k∈K,i∈Vtuc,(10)dtku+d0ivu+M⋅1−xk000iu≥atkiu,k∈K,i∈Vc,(11)dtku+d0ivu−M⋅1−xk000iu≤atkiu,k∈K,i∈Vc,(12)atkiu+wtkiu+si+dijvu+M⋅1−xkrciju≥atkju,k∈K,r∈Vtu,c∈Vtu,i∈V,i≠r,j∈V,(13)atkiu+wtkiu+si+dijvu−M⋅1−xkrciju≤atkju,k∈K,r∈Vtu,c∈Vtu,i∈V,i≠r,j∈V,(14)atkju+wtkju+M⋅1−∑i∈Vxkrciju≥ej,k∈K,r∈Vtu,c∈Vtu,j∈Vc,j≠r,j≠c,(15)atkju+wtkju−M⋅1−∑i∈Vxkrciju≤lj,k∈K,r∈Vtu,c∈Vtu,j∈Vc,j≠r,j≠c.Constraints on the vehicle visiting vertices ensure the limit on the times of vehicle visiting a node and the vehicle retention. Constraints (16)–(18) ensure that a TUC is visited exactly once by a truck or a UV when the TUC is not being used as the satellite or a TUC is visited more than once if it is being used as the satellite. Constraint (19) guarantees that each UC is visited exactly once by one UV. Constraint (20) ensures that for a truck, the number of arrivals at vertex j equals the number of departures from vertex j. Constraint (21) guarantees that for a UV, the number of arrivals at a customer vertex equals the number of departures from the customer vertex. According to constraint (22), each dispatched UV leaves once from the origin of a subtour and returns once to the destination of the subtour. Constraint (23) requires that a UV departs by itself or with the paired truck no more than once from the depot.(16)∑k∈K∑i∈Vtuxkijt+∑k∈K∑r∈Vtu∑c∈Vtu∑i∈Vxkrciju≥1,j∈Vtuc,(17)∑k∈K∑i∈Vtuxtijt≤1,j∈Vtuc,(18)∑k∈K∑i∈Vtu∑r∈Vtu∑c∈Vtuxkrciju≤1,j∈Vtuc,(19)∑k∈K∑r∈Vtu∑c∈Vtu∑i∈Vxkrciju=1,j∈Vuc,(20)∑i∈Vtuxkijt=∑a∈Vtuxkjat,k∈K,j∈Vtu,(21)∑i∈Vxkrciju=∑i∈Vxkrcjiu,k∈K,r∈Vtu,c∈Vtu,j∈V,j≠r,j≠c,(22)∑i∈Vxkrcriu=∑j∈Vxkrcjcu,k∈K,r∈Vtu,c∈Vtu,i,j≠r,i,j≠c,(23)∑j∈Vtucxk0jt+∑j∈Vcxk000ju≤1,k∈K.Constraints on the intersatellite synchronization ensure the connection between the main tour and its subtour(s). The intersatellite synchronization permits the origin and the destination of a subtour to differ. It is assumed that the truck arrives once at and leaves once from the subtour origin, and the truck arrives once and leaves once from the subtour destination. The dispatched UV leaves once from the subtour origin and arrives once at the subtour destination. After dispatching the carried UV and finishing the service at a satellite, the truck continues to visit other TUCs to complete each intended service. Only if the arrival time of the dispatched UV at a satellite is between the arrival time of the paired truck at the satellite and the departure time of the paired truck UV from the satellite can the truck pick up its UV. According to constraint (24), a TUC should be included in a main tour if the TUC is used as the subtour origin. According to constraint (25), the subtour destination should be included in a main tour. From the temporal perspective, constraints (26)–(28) ensure that a UV and the paired truck meet at the appointed time. According to constraint (26), the arrival time of a UV at the subtour destination is not later than the departure time of the paired truck from the subtour destination. If the depot is used as the subtour destination, the appointed time for a UV meeting the paired truck is not restricted. According to constraints (27) and (28), the arrival time of a dispatched UV at the first customer in a subtour is determined by the arrival time of the truck-UV combination at the satellite and the UV running time on the arc.(24)∑j∈Vc∑c∈Vtuxkrcrju≤∑i∈Vtuxkirt,k∈K,r∈Vtuc,(25)∑i∈Vc∑r∈Vtucxkrcicu≤∑i∈Vtuxkict,k∈K,c∈Vtu,(26)atkct+wtkct+sc+M⋅1−∑r∈Vtuc∑i∈Vcxkrcicu≥atkcu,k∈K,c∈Vtuc,(27)atkrt+drivu+M⋅1−xkrcriu≥atkiu,k∈K,r∈Vtuc,c∈Vtu,i∈Vc,(28)atkrt+drivu−M⋅1−xkrcriu≤atkiu,k∈K,r∈Vtuc,c∈Vtu,i∈Vc.The results of the computational experiments on small-scale instances demonstrate that the TUVRP-TW formulation can be solved directly by commercial IP solver software (CPLEX in this paper). The exact solutions of small-scale instances that are found by CPLEX, which can be referenced in calibrating the heuristic, support the validity of the TUVRP-TW formulation. Instances are seldom solved directly via CPLEX in practical applications because practical instances are of typically large scale and too much computer memory and computation time are required by computers. To find satisfactory solutions for large-scale instances, we propose a hybrid algorithm based on a GRASP and a VNS. On small-scale instances, the exact solutions that are found using CPLEX and the satisfactory solutions that are output by the hybrid algorithm are compared to evaluate the performance of the heuristic. Large-scale instances are solved using the proposed heuristic.


## 3.2.1. Objective Function

The objective function consists of four parts: the travelling distance of trucks, the transformed distance of possible waiting times of trucks at customers, the travelling distance of UVs, and the transformed distance of possible waiting times of UVs at customers. The objective function is expressed as(1)min∑k∈K∑i∈Vtu∑j∈Vtudij⋅xkijt+∑k∈K∑i∈Vcwtkit⋅vt⋅τ+∑k∈K∑r∈Vtu∑c∈Vtu∑i∈V∑j∈Vdij⋅xkrciju+∑k∈K∑i∈Vcwtkiu⋅vu⋅τ.


## 3.2.2. Constraints

We divide the constraints into four categories: constraints on the vehicle capacity, constraints on the time continuity, constraints on the vehicle visiting vertices, and constraints on the intersatellite synchronization.Constraint (2) restricts the capacity of each UV in UV routes for direct delivery or in subtours by UVs that have been dispatched from trucks at satellites. Constraint (3) restricts the capacity of each truck-UV combination in combination routes or truck routes.(2)∑i∈Vc,i≠r,i≠c∑j∈Vqixkrciju≤Qu,k∈K,r∈Vtu,c∈Vtu,(3)∑i∈Vtuc∑j∈Vtuqixkijt+∑r∈Vtuc∑c∈Vtu∑i∈Vc,i≠r,i≠c∑j∈Vqixkrciju≤Qtu,k∈K.Constraints on the time continuity aim to decide the arrival, leaving, and waiting time for vehicles in use. According to constraints (4) and (5), the truck arrival time at the first TUC is determined by the truck departure time from the depot and the running time on the arc that connects the depot and the first TUC. According to constraints (6) and (7), the truck arrival time at a specified TUC or at the depot is determined by the truck arrival time at previous vertex, the possible waiting time for time window opening, the service time, and the running time on the arc. Constraints (8) and (9) ensure that the truck visits TUCs by respecting the TUC time windows. According to constraints (10) and (11), the UV arrival time at the first customer is determined by the UV departure time from the depot and the running time on the arc that connects the depot and the first customer. According to constraints (12) and (13), the UV arrival time at a specified vertex is determined by the UV arrival time at previous vertex, the possible waiting time for time window opening, the service time, and the UV running time on the arc. Constraints (14) and (15) ensure that the dispatched UV visits customers by respecting the customer time windows.(4)dtkt+d0ivt+M⋅1−xk0it≥atkit,k∈K,i∈Vtuc,(5)dtkt+d0ivt−M⋅1−xk0it≤atkit,k∈K,i∈Vtuc,(6)atkit+wtkit+si+dijvt+M⋅1−xkijt≥atkjt,k∈K,i∈Vtuc,j∈Vtu,(7)atkit+wtkit+si+dijvt−M⋅1−xkijt≤atkjt,k∈K,i∈Vtuc,j∈Vtu,(8)atkit+wtkit+M⋅1−∑j∈Vtuxkijt≥ei,k∈K,i∈Vtuc,(9)atkit+wtkit−M⋅1−∑j∈Vtuxkijt≤li,k∈K,i∈Vtuc,(10)dtku+d0ivu+M⋅1−xk000iu≥atkiu,k∈K,i∈Vc,(11)dtku+d0ivu−M⋅1−xk000iu≤atkiu,k∈K,i∈Vc,(12)atkiu+wtkiu+si+dijvu+M⋅1−xkrciju≥atkju,k∈K,r∈Vtu,c∈Vtu,i∈V,i≠r,j∈V,(13)atkiu+wtkiu+si+dijvu−M⋅1−xkrciju≤atkju,k∈K,r∈Vtu,c∈Vtu,i∈V,i≠r,j∈V,(14)atkju+wtkju+M⋅1−∑i∈Vxkrciju≥ej,k∈K,r∈Vtu,c∈Vtu,j∈Vc,j≠r,j≠c,(15)atkju+wtkju−M⋅1−∑i∈Vxkrciju≤lj,k∈K,r∈Vtu,c∈Vtu,j∈Vc,j≠r,j≠c.Constraints on the vehicle visiting vertices ensure the limit on the times of vehicle visiting a node and the vehicle retention. Constraints (16)–(18) ensure that a TUC is visited exactly once by a truck or a UV when the TUC is not being used as the satellite or a TUC is visited more than once if it is being used as the satellite. Constraint (19) guarantees that each UC is visited exactly once by one UV. Constraint (20) ensures that for a truck, the number of arrivals at vertex j equals the number of departures from vertex j. Constraint (21) guarantees that for a UV, the number of arrivals at a customer vertex equals the number of departures from the customer vertex. According to constraint (22), each dispatched UV leaves once from the origin of a subtour and returns once to the destination of the subtour. Constraint (23) requires that a UV departs by itself or with the paired truck no more than once from the depot.(16)∑k∈K∑i∈Vtuxkijt+∑k∈K∑r∈Vtu∑c∈Vtu∑i∈Vxkrciju≥1,j∈Vtuc,(17)∑k∈K∑i∈Vtuxtijt≤1,j∈Vtuc,(18)∑k∈K∑i∈Vtu∑r∈Vtu∑c∈Vtuxkrciju≤1,j∈Vtuc,(19)∑k∈K∑r∈Vtu∑c∈Vtu∑i∈Vxkrciju=1,j∈Vuc,(20)∑i∈Vtuxkijt=∑a∈Vtuxkjat,k∈K,j∈Vtu,(21)∑i∈Vxkrciju=∑i∈Vxkrcjiu,k∈K,r∈Vtu,c∈Vtu,j∈V,j≠r,j≠c,(22)∑i∈Vxkrcriu=∑j∈Vxkrcjcu,k∈K,r∈Vtu,c∈Vtu,i,j≠r,i,j≠c,(23)∑j∈Vtucxk0jt+∑j∈Vcxk000ju≤1,k∈K.Constraints on the intersatellite synchronization ensure the connection between the main tour and its subtour(s). The intersatellite synchronization permits the origin and the destination of a subtour to differ. It is assumed that the truck arrives once at and leaves once from the subtour origin, and the truck arrives once and leaves once from the subtour destination. The dispatched UV leaves once from the subtour origin and arrives once at the subtour destination. After dispatching the carried UV and finishing the service at a satellite, the truck continues to visit other TUCs to complete each intended service. Only if the arrival time of the dispatched UV at a satellite is between the arrival time of the paired truck at the satellite and the departure time of the paired truck UV from the satellite can the truck pick up its UV. According to constraint (24), a TUC should be included in a main tour if the TUC is used as the subtour origin. According to constraint (25), the subtour destination should be included in a main tour. From the temporal perspective, constraints (26)–(28) ensure that a UV and the paired truck meet at the appointed time. According to constraint (26), the arrival time of a UV at the subtour destination is not later than the departure time of the paired truck from the subtour destination. If the depot is used as the subtour destination, the appointed time for a UV meeting the paired truck is not restricted. According to constraints (27) and (28), the arrival time of a dispatched UV at the first customer in a subtour is determined by the arrival time of the truck-UV combination at the satellite and the UV running time on the arc.(24)∑j∈Vc∑c∈Vtuxkrcrju≤∑i∈Vtuxkirt,k∈K,r∈Vtuc,(25)∑i∈Vc∑r∈Vtucxkrcicu≤∑i∈Vtuxkict,k∈K,c∈Vtu,(26)atkct+wtkct+sc+M⋅1−∑r∈Vtuc∑i∈Vcxkrcicu≥atkcu,k∈K,c∈Vtuc,(27)atkrt+drivu+M⋅1−xkrcriu≥atkiu,k∈K,r∈Vtuc,c∈Vtu,i∈Vc,(28)atkrt+drivu−M⋅1−xkrcriu≤atkiu,k∈K,r∈Vtuc,c∈Vtu,i∈Vc.The results of the computational experiments on small-scale instances demonstrate that the TUVRP-TW formulation can be solved directly by commercial IP solver software (CPLEX in this paper). The exact solutions of small-scale instances that are found by CPLEX, which can be referenced in calibrating the heuristic, support the validity of the TUVRP-TW formulation. Instances are seldom solved directly via CPLEX in practical applications because practical instances are of typically large scale and too much computer memory and computation time are required by computers. To find satisfactory solutions for large-scale instances, we propose a hybrid algorithm based on a GRASP and a VNS. On small-scale instances, the exact solutions that are found using CPLEX and the satisfactory solutions that are output by the hybrid algorithm are compared to evaluate the performance of the heuristic. Large-scale instances are solved using the proposed heuristic.


## 3.3. A Variant without Intersatellite Synchronization

Considering that truck-UV combinations are used for “last-meter” delivery in the “mothership approach” of truck and UV delivery modes, we introduce the TUVRP-TW. The intersatellite synchronization in the TUVRP-TW allows a truck and its UV to meet at different satellites. Due to the involvement of the intersatellite synchronization, the TUVRP-TW differs from the 2E-VRP variants in the literature. In the 2E-VRP variant literature, satellite synchronization constraints require that first-echelon vehicles and second-echelon vehicles meet at the same satellite. The inner-satellite synchronization permits synchronization operations at the same satellite.We also introduce and formulate a routing variant with inner-satellite synchronization, named the TUVRP-TW-B. In the TUVRP-TW-B, intersatellite synchronization is not permitted. The TUVRP-TW-B formulation is presented in the Appendix.


## 4. A GRASP and the VNS

The TUVRP-TW is complex because of the intersatellite synchronization and time continuity constraints. Since only small instances of the TUVRP-TW could be solved directly by using the mathematical formulation in Section3, we provide a hybrid algorithm based on a GRASP and a VNS. The hybrid algorithm incorporates two powerful features, the effective construction and improving ability of the GRASP, and the flexibility of the VNS to explore different search spaces.The GRASP consists of iterations made up from successive constructions of a greedy randomized solution and subsequent iterative improvements of it. The greedy randomized solutions are generated by adding elements to the solution set, and the list of elements is ranked by a greedy function according to the quality of the solution they will achieve. To obtain variability in the candidate set of greedy solutions, well-ranked candidate elements are placed in a restricted candidate list (RCL) and are randomly chosen when constructing the solution.The VNS, proposed by Hansen and Mladenović [33], is an effective metaheuristic for solving combinatorial optimization problems. The VNS can exploit systematically the neighborhood change, both in descent to local minima and in escape from the valleys that contain them [34]. The basic idea is that a local optimum defined by one neighborhood structure is not necessarily the local optimum of another neighborhood structure, and thus the search can systematically traverse different search spaces which are defined by different neighborhood structures.


### 4.1. Constructing an Initial Solution by a GRASP

There are three types of routes in the TUVRP-TW solution: UV routes, truck routes, and combination routes. From the perspective of node insertion in solution construction operation, a TUC may be inserted into a truck route, a UV route, a main tour of a combination route, or a subtour of a combination route; a UC may be inserted into a UV route or a subtour of a combination route. In order to construct routes of an initial solution, the GRASP adopts different adaptive probabilistic mechanisms for TUCs and UCs. Considering that the types of TUC insertion operations are more than those of UC insertion operations, a reactive mechanism is adopted for TUC insertion operations, and a random plus greedy mechanism is adopted for UC insertion operations. The following ways to insert a TUC are adopted, i.e., constructing a new UV route, inserting into a constructed UV route, constructing a new truck route, inserting into a constructed truck route, inserting into a constructed subtour, and constructing a new subtour. The following ways to insert a UC are adopted, i.e., constructing a new UV route, inserting into a constructed UV route, constructing a new subtour, and inserting into a constructed subtour. The insertion cost of a customer is estimated as the travelling distance and the transformed distance of possible waiting times of vehicles at customers.The capacity constraints are firstly used to decide the feasibility of the insertion operation. To handle the time windows, (i) when constructing a new truck route, constructing a new UV route or inserting into a constructed UV route is adopted, the departure time of the route is initialized as the opening time of the depot, and then the arrival time and vehicle waiting time at each customer in the route are estimated. The departure time of the route is postponed by the method provided in Li et al. [35], so as to decrease vehicle waiting time at each customer as much as possible, while customer time windows are guaranteed. (ii) When the way of constructing a new subtour or inserting into a constructed subtour is adopted, the opening time of the depot is used as the earliest departure time of the corresponding main tour, and the latest departure time of the corresponding main tour is estimated as the postponed time by the method provided in Li et al. [35]. The earliest and latest departure times of the main tour are referenced to estimate the departure time of the subtour, and the result is not accepted if customer time windows are not ensured or if the UV covering the subtour cannot arrive at the subtour destination before its truck leaves. (iii) When the way of inserting into a constructed truck route is adopted, the handling method that is same as that used in the way of constructing a new truck route is used if there is no subtour involved. If there is one or several subtours corresponding to the constructed truck route, the handling method is the same as that used in the way of constructing a new subtour.To handle the number of truck-UV pairs, when the way of constructing a new subtour or inserting into a constructed subtour is adopted, at a TUC, only if there is a UV carried by a truck, the TUC can be used as the origin of a subtour; at a TUC, only if there is no UV on a truck, the TUC can be used as the destination of a subtour.We name the insertion of customeri by a way for customer insertion as temp_in and the insertion cost change of inserting customer i or not as temp_obj. The best insertion of customer i and its cost change are named best_in (i) and best_obj (i), respectively. For each customer insertion way, every location of each route is enumerated, for customer i. If a temp_in results in a solution that cannot ensure the TUVRP-TW constraints, the corresponding temp_obj is estimated as a large enough number. If a temp_in leads to “temp_obj < best_obj (i),” best_in (i) and best_obj (i) are updated. If temp_in leads to “temp_obj = best_obj (i),” best_in (i) and best_obj (i) are updated according to the probability of 50%.The GRASP process is shown in Figure3. The GRASP includes the following main steps.Step 1. Initialize the solution (denoted as initial_S) and the unserved customer set (denoted as CustomerList).Step 2. To insert customer i (i ∈ CustomerList), each feasible position is referred to estimate the insertion cost, and the minimum of insertion cost of customer i is identified and denoted as ce (i), and ce (i) = best_obj (i).Step 3. Rank the unserved customers by minimum of insertion costs and introduce cemin ← min{ce (i)|i∈CustomerList}, and cemax ← max{ce (i)|i∈CustomerList}.Step 4. Define the RCL set, RCL ← {i|i∈CustomerList, and ce (i) ≤ cemin + α·(cemax − cemin)}, where α ∈ [0, 1].Step 5. Randomly select a customer in the RCL set for customer insertion and renew initial_S and CustomerList.Figure 3
The GRASP process.Repeat Steps 2–5 untilCustomerList is empty.A feasible solution is constructed iteratively, one customer at a time. At each iteration, all unserved customers are ordered according to the insertion cost, and then the customers are selected to join the RCL with respect to an adaptive or greedy function. The customer to be inserted is selected randomly from the RCL set. The RCL set is updated at each iteration to reflect the changes brought on by the insertion of the previous customer, which indicates the adaptive mechanism. The probabilistic character is indicated by randomly selecting one of the customers in the RCL set.For UC insertion, a random plus greedy mechanism is adopted. Instead of combining greediness and randomness at each iteration, the random plus greedy scheme applies randomness during the firstp iterations to produce a random partial solution, and the parameter α = 1. Then, one or more pure greedy iterations construct the feasible solution, and the parameter α = 0. The balance between greediness and randomness can be controlled by changing the value of the parameter p. In the paper, the value of the parameter p is 50% of the number of UCs in the TUVRP-TW network.For TUC insertion, a reactive mechanism is adopted. The parameterα that is involved in defining the RCL set is not fixed but is selected at each iteration from a set of discrete values. The value of the parameter α is selected by referring to the solution quality found along the previous iterations. Let the set of discrete values of α be {0.0, 0.1, 0.2, 0.3, …, 0.8, 0.9, 1.0}. We denote the selection times of the i-th value for α as usetimes (i). We denote the incumbent solution at each iteration with the i-th value of α as aobj (i) and the average of objective values of all solutions found by using the i-th value for α as aq (i). Let aq (i) = aobj (i)/usetimes (i). The probabilities ap (i) associated with the selection of each value are initialized as ap (i) = 1/11. usetimes (i) and aq (i) are initialized as 0. When the i-th value of α is selected, usetime (i) is increased by 1, and aq (i) is increased by the insertion cost of the selected customer. The selection probability of the i-th value of α is periodically reevaluated by taking ap (i) = aq (i)/∑aq (i). The probabilities associated with the more appropriate values will increase when they are reevaluated.


### 4.2. VNS

#### 4.2.1. The VNS Process

A VNS procedure is used to improve the initial solution (i.e.,initial_S). We set the number of external iterations as Nout and the number of internal iterations as Nin. We define the optimization process as Improving and the disturbance process as Shaking. Within the Shaking process, a poor solution may be accepted, provided that its deviation does not exceed a given threshold (denoted as gap). Nout determines the number of shaking process, and Nin determines the number of internal iterations after the solution perturbation. The VNS combines deterministic and stochastic changes of neighborhood. The deterministic part is represented by a local search heuristic, and the stochastic phase is represented by the random selection of one point in one of the neighborhoods (Villegas et al., 2011).At the beginning of the external iterations,initial_S obtained by the GRASP is regarded as the incumbent solution (denoted as Solution) and the best solution (denoted as best_S). At each external iteration, Solution is perturbed by Shaking. Based on the selected operators of Shaking, Solution is changed, and a solution (denoted as shaking_S), which has the most disturbance within gap, is retained. After the Shaking, Improving is performed, which is iterated Nin times. If the attained result is worse than Shaking_S, Shaking_S is retained; if the attained result is better than Shaking_S, the better solution is retained. After each iteration, it will be judged whether the incumbent solution is optimized. If so, the best solution attained so far will be updated. We denote a finite set of preselected neighborhood structures with Nk1 (k1 = 1, …, kmax1) in shaking and Nk2 (k2 = 1, …, kmax2) in Improving. We introduce Nki (x) (i = 1, 2) to denote the set of solutions in the k1-th or k2-th neighborhood of x, where x is a feasible solution.The process of the whole VNS for solving the TUVRP-TW is shown in Figure4 and is described in pseudo-code in Figure 5.Figure 4
The VNS process.Figure 5
The VNS algorithm.InImproving and Shaking, various types of local search operators are used to search neighborhoods. The objective function of the TUVRP-TW is used as the estimation criterion for substituting the newly generated solution for the incumbent solution. Operators are designed based on relocation operators and interchange operators, which are typically used in the search phase for the VRP with time windows (VRPTW) in the literature. Operators for combination routes probably change the time continuity due to the intersatellite synchronization and the varying subtour origins and destinations.InShaking, seven operators are adopted, that is, a relocation operator for truck routes or main tours of combination routes, a relocation operator for UV routes, a relocation operator for subtours, and four interchange operators for two randomly selected UV routes (denoted as (1, 0)PUR, (1, 1)PUR, (1, 2)PUR, and (2, 0)PUR).Eighteen operators are adopted inImproving, that is, a relocation operator for truck routes and main tours of combination routes, a relocation operator for UV routes, a relocation operator for subtours, four interchange operators for two randomly-selected UV routes (denoted as (1, 0)PUR, (1, 1)PUR, (1, 2)PUR and (2, 0)PUR), four interchange operators for two randomly-selected truck routes or main tours (denoted as (1, 0)PTR, (1, 1)PTR, (1, 2)PTR and (2, 0)PTR), two interchange operators for two randomly-selected truck route (or main tour) and UV route (denoted as (1,0)PTR_PUR, (0,1)PTR_PUR, (1,1)PTR_PUR), two interchange operators for a randomly-selected subtour and a randomly-selected UV route (denoted as (1, 0)PUR_SUB and (1, 1)PUR_SUB), and two interchange operators for a randomly-selected subtour and a randomly-selected truck route(or main tour) (denoted as (1, 0)PTR_SUB and (1, 1)PTR_SUB).


#### 4.2.2. Local Search Operators

The operators for truck routes or UV routes are similar to those usually adopted in the search phase for the routing variants involving time windows in the literature (e.g., [35]). The operators involving combination routes are relatively complex because of the intersatellite synchronization and time constraints. We take several routes shown in Figures 6–9 as examples to explain several special operators. In Figures 6–9, TUC is denoted as the box with the serial number of a customer. UC is denoted as the circle with the serial number of a customer.Figure 6
An example of the relocate operator for subtours.Figure 7
Examples of the interchange operators for one main tour (or truck route) and one UV route.
(a)(b)(c)Figure 8
Examples of the interchange operators for one subtour and one UV route.
(a)(b)Figure 9
Examples of the interchange operators for one main tour (or truck route) and one subtour.
(a)(b)Relocation operators aim to relocate a vertex of a randomly selected truck route, UV route, main tour, or subtour of a combination route. Figure6 shows an illustrative example of the relocate operator for subtours. A combination route including the main tour “0-1-2-3-4-5-6-0” and the subtour “2-7-8-6” is randomly selected. A randomly selected vertex (“7”) in the subtour is relocated, which results in a new feasible subtour “2-8-7-5” with the destination of the subtour being adjusted at vertex “5.” It is worth noting that the relocate operator for subtours probably leads to the changing time continuity for the intersatellite synchronization and the subtour destination varying.Interchange operators for two randomly selected routes aim at selecting vertices in each of the two routes for exchange. Figure7(a) gives an example of the (1, 0)PTR_PUR operator. A TUC “2” in the main tour of the combination route is randomly selected and inserted into the UV route, on condition of respecting the TUVRP-TW constraints. Figure 7(b) gives an example of the (0, 1)PTR_PUR operator. A TUC “8” in the UV route is randomly selected and inserted into the main tour of the combination route, on condition of respecting the TUVRP-TW constraints. The main tour “0-2-1-3-6-7-4-0” is adjusted as “0-2-1-3-6-7-4-8-0,” and the subtour “3-11-5-10-0” is transformed into “3-11-5-10-8.” Figure 7(c) gives an example of the (1, 1)PTR_PUR operator. A TUC “1” in the main tour of the combination route is randomly selected, and a TUC “8” included in the UV route is randomly selected. The two selected TUCs are exchanged with each other.Figure8(a) gives an example of the (1, 0)PUR_SUB operator. A customer “12” in the UV route is randomly selected and assigned as the customer of TUC “2” included in the main tour of the combination route, on condition of respecting the TUVRP-TW constraints. A new subtour “2-12-1” is constructed. The UV route “0-12-8-9-13-0” is transformed into “0-8-9-13-0”. Figure 8(b) gives an example of the (1, 1)PUR_SUB operator. A customer “11” included in a subtour of the combination route is randomly selected, and a customer “8” included in the UV route is randomly selected. The two selected customers are exchanged with each other.Figure9(a) gives an example of the (1, 1)PTR_SUB operator. In the main tour, a TUC “1” that is used as the origin of a subtour is randomly selected, and a TUC “4” included in the subtour is randomly selected. The TUCs “1” and “4” are exchanged with each other. Respecting the TUVRP-TW constraints, the main tour “0-1-2-3-7-8-11-12-0” is transformed into “0-4-2-3-7-8-11-12-0,” and the subtour “1-4-5-3” is transformed into “4-1-5-7.” The TUCs “4” and “7” are renewed as the origin and destination of the subtour, respectively. Figure 9(b) gives an example of the (1, 0)PTR_SUB operator. A TUC “7” in the main tour is randomly selected and inserted into a randomly selected subtour, on the condition that the TUVRP-TW constraints are ensured.


## 4.1. Constructing an Initial Solution by a GRASP

There are three types of routes in the TUVRP-TW solution: UV routes, truck routes, and combination routes. From the perspective of node insertion in solution construction operation, a TUC may be inserted into a truck route, a UV route, a main tour of a combination route, or a subtour of a combination route; a UC may be inserted into a UV route or a subtour of a combination route. In order to construct routes of an initial solution, the GRASP adopts different adaptive probabilistic mechanisms for TUCs and UCs. Considering that the types of TUC insertion operations are more than those of UC insertion operations, a reactive mechanism is adopted for TUC insertion operations, and a random plus greedy mechanism is adopted for UC insertion operations. The following ways to insert a TUC are adopted, i.e., constructing a new UV route, inserting into a constructed UV route, constructing a new truck route, inserting into a constructed truck route, inserting into a constructed subtour, and constructing a new subtour. The following ways to insert a UC are adopted, i.e., constructing a new UV route, inserting into a constructed UV route, constructing a new subtour, and inserting into a constructed subtour. The insertion cost of a customer is estimated as the travelling distance and the transformed distance of possible waiting times of vehicles at customers.The capacity constraints are firstly used to decide the feasibility of the insertion operation. To handle the time windows, (i) when constructing a new truck route, constructing a new UV route or inserting into a constructed UV route is adopted, the departure time of the route is initialized as the opening time of the depot, and then the arrival time and vehicle waiting time at each customer in the route are estimated. The departure time of the route is postponed by the method provided in Li et al. [35], so as to decrease vehicle waiting time at each customer as much as possible, while customer time windows are guaranteed. (ii) When the way of constructing a new subtour or inserting into a constructed subtour is adopted, the opening time of the depot is used as the earliest departure time of the corresponding main tour, and the latest departure time of the corresponding main tour is estimated as the postponed time by the method provided in Li et al. [35]. The earliest and latest departure times of the main tour are referenced to estimate the departure time of the subtour, and the result is not accepted if customer time windows are not ensured or if the UV covering the subtour cannot arrive at the subtour destination before its truck leaves. (iii) When the way of inserting into a constructed truck route is adopted, the handling method that is same as that used in the way of constructing a new truck route is used if there is no subtour involved. If there is one or several subtours corresponding to the constructed truck route, the handling method is the same as that used in the way of constructing a new subtour.To handle the number of truck-UV pairs, when the way of constructing a new subtour or inserting into a constructed subtour is adopted, at a TUC, only if there is a UV carried by a truck, the TUC can be used as the origin of a subtour; at a TUC, only if there is no UV on a truck, the TUC can be used as the destination of a subtour.We name the insertion of customeri by a way for customer insertion as temp_in and the insertion cost change of inserting customer i or not as temp_obj. The best insertion of customer i and its cost change are named best_in (i) and best_obj (i), respectively. For each customer insertion way, every location of each route is enumerated, for customer i. If a temp_in results in a solution that cannot ensure the TUVRP-TW constraints, the corresponding temp_obj is estimated as a large enough number. If a temp_in leads to “temp_obj < best_obj (i),” best_in (i) and best_obj (i) are updated. If temp_in leads to “temp_obj = best_obj (i),” best_in (i) and best_obj (i) are updated according to the probability of 50%.The GRASP process is shown in Figure3. The GRASP includes the following main steps.Step 1. Initialize the solution (denoted as initial_S) and the unserved customer set (denoted as CustomerList).Step 2. To insert customer i (i ∈ CustomerList), each feasible position is referred to estimate the insertion cost, and the minimum of insertion cost of customer i is identified and denoted as ce (i), and ce (i) = best_obj (i).Step 3. Rank the unserved customers by minimum of insertion costs and introduce cemin ← min{ce (i)|i∈CustomerList}, and cemax ← max{ce (i)|i∈CustomerList}.Step 4. Define the RCL set, RCL ← {i|i∈CustomerList, and ce (i) ≤ cemin + α·(cemax − cemin)}, where α ∈ [0, 1].Step 5. Randomly select a customer in the RCL set for customer insertion and renew initial_S and CustomerList.Figure 3
The GRASP process.Repeat Steps 2–5 untilCustomerList is empty.A feasible solution is constructed iteratively, one customer at a time. At each iteration, all unserved customers are ordered according to the insertion cost, and then the customers are selected to join the RCL with respect to an adaptive or greedy function. The customer to be inserted is selected randomly from the RCL set. The RCL set is updated at each iteration to reflect the changes brought on by the insertion of the previous customer, which indicates the adaptive mechanism. The probabilistic character is indicated by randomly selecting one of the customers in the RCL set.For UC insertion, a random plus greedy mechanism is adopted. Instead of combining greediness and randomness at each iteration, the random plus greedy scheme applies randomness during the firstp iterations to produce a random partial solution, and the parameter α = 1. Then, one or more pure greedy iterations construct the feasible solution, and the parameter α = 0. The balance between greediness and randomness can be controlled by changing the value of the parameter p. In the paper, the value of the parameter p is 50% of the number of UCs in the TUVRP-TW network.For TUC insertion, a reactive mechanism is adopted. The parameterα that is involved in defining the RCL set is not fixed but is selected at each iteration from a set of discrete values. The value of the parameter α is selected by referring to the solution quality found along the previous iterations. Let the set of discrete values of α be {0.0, 0.1, 0.2, 0.3, …, 0.8, 0.9, 1.0}. We denote the selection times of the i-th value for α as usetimes (i). We denote the incumbent solution at each iteration with the i-th value of α as aobj (i) and the average of objective values of all solutions found by using the i-th value for α as aq (i). Let aq (i) = aobj (i)/usetimes (i). The probabilities ap (i) associated with the selection of each value are initialized as ap (i) = 1/11. usetimes (i) and aq (i) are initialized as 0. When the i-th value of α is selected, usetime (i) is increased by 1, and aq (i) is increased by the insertion cost of the selected customer. The selection probability of the i-th value of α is periodically reevaluated by taking ap (i) = aq (i)/∑aq (i). The probabilities associated with the more appropriate values will increase when they are reevaluated.


## 4.2. VNS

### 4.2.1. The VNS Process

A VNS procedure is used to improve the initial solution (i.e.,initial_S). We set the number of external iterations as Nout and the number of internal iterations as Nin. We define the optimization process as Improving and the disturbance process as Shaking. Within the Shaking process, a poor solution may be accepted, provided that its deviation does not exceed a given threshold (denoted as gap). Nout determines the number of shaking process, and Nin determines the number of internal iterations after the solution perturbation. The VNS combines deterministic and stochastic changes of neighborhood. The deterministic part is represented by a local search heuristic, and the stochastic phase is represented by the random selection of one point in one of the neighborhoods (Villegas et al., 2011).At the beginning of the external iterations,initial_S obtained by the GRASP is regarded as the incumbent solution (denoted as Solution) and the best solution (denoted as best_S). At each external iteration, Solution is perturbed by Shaking. Based on the selected operators of Shaking, Solution is changed, and a solution (denoted as shaking_S), which has the most disturbance within gap, is retained. After the Shaking, Improving is performed, which is iterated Nin times. If the attained result is worse than Shaking_S, Shaking_S is retained; if the attained result is better than Shaking_S, the better solution is retained. After each iteration, it will be judged whether the incumbent solution is optimized. If so, the best solution attained so far will be updated. We denote a finite set of preselected neighborhood structures with Nk1 (k1 = 1, …, kmax1) in shaking and Nk2 (k2 = 1, …, kmax2) in Improving. We introduce Nki (x) (i = 1, 2) to denote the set of solutions in the k1-th or k2-th neighborhood of x, where x is a feasible solution.The process of the whole VNS for solving the TUVRP-TW is shown in Figure4 and is described in pseudo-code in Figure 5.Figure 4
The VNS process.Figure 5
The VNS algorithm.InImproving and Shaking, various types of local search operators are used to search neighborhoods. The objective function of the TUVRP-TW is used as the estimation criterion for substituting the newly generated solution for the incumbent solution. Operators are designed based on relocation operators and interchange operators, which are typically used in the search phase for the VRP with time windows (VRPTW) in the literature. Operators for combination routes probably change the time continuity due to the intersatellite synchronization and the varying subtour origins and destinations.InShaking, seven operators are adopted, that is, a relocation operator for truck routes or main tours of combination routes, a relocation operator for UV routes, a relocation operator for subtours, and four interchange operators for two randomly selected UV routes (denoted as (1, 0)PUR, (1, 1)PUR, (1, 2)PUR, and (2, 0)PUR).Eighteen operators are adopted inImproving, that is, a relocation operator for truck routes and main tours of combination routes, a relocation operator for UV routes, a relocation operator for subtours, four interchange operators for two randomly-selected UV routes (denoted as (1, 0)PUR, (1, 1)PUR, (1, 2)PUR and (2, 0)PUR), four interchange operators for two randomly-selected truck routes or main tours (denoted as (1, 0)PTR, (1, 1)PTR, (1, 2)PTR and (2, 0)PTR), two interchange operators for two randomly-selected truck route (or main tour) and UV route (denoted as (1,0)PTR_PUR, (0,1)PTR_PUR, (1,1)PTR_PUR), two interchange operators for a randomly-selected subtour and a randomly-selected UV route (denoted as (1, 0)PUR_SUB and (1, 1)PUR_SUB), and two interchange operators for a randomly-selected subtour and a randomly-selected truck route(or main tour) (denoted as (1, 0)PTR_SUB and (1, 1)PTR_SUB).


### 4.2.2. Local Search Operators

The operators for truck routes or UV routes are similar to those usually adopted in the search phase for the routing variants involving time windows in the literature (e.g., [35]). The operators involving combination routes are relatively complex because of the intersatellite synchronization and time constraints. We take several routes shown in Figures 6–9 as examples to explain several special operators. In Figures 6–9, TUC is denoted as the box with the serial number of a customer. UC is denoted as the circle with the serial number of a customer.Figure 6
An example of the relocate operator for subtours.Figure 7
Examples of the interchange operators for one main tour (or truck route) and one UV route.
(a)(b)(c)Figure 8
Examples of the interchange operators for one subtour and one UV route.
(a)(b)Figure 9
Examples of the interchange operators for one main tour (or truck route) and one subtour.
(a)(b)Relocation operators aim to relocate a vertex of a randomly selected truck route, UV route, main tour, or subtour of a combination route. Figure6 shows an illustrative example of the relocate operator for subtours. A combination route including the main tour “0-1-2-3-4-5-6-0” and the subtour “2-7-8-6” is randomly selected. A randomly selected vertex (“7”) in the subtour is relocated, which results in a new feasible subtour “2-8-7-5” with the destination of the subtour being adjusted at vertex “5.” It is worth noting that the relocate operator for subtours probably leads to the changing time continuity for the intersatellite synchronization and the subtour destination varying.Interchange operators for two randomly selected routes aim at selecting vertices in each of the two routes for exchange. Figure7(a) gives an example of the (1, 0)PTR_PUR operator. A TUC “2” in the main tour of the combination route is randomly selected and inserted into the UV route, on condition of respecting the TUVRP-TW constraints. Figure 7(b) gives an example of the (0, 1)PTR_PUR operator. A TUC “8” in the UV route is randomly selected and inserted into the main tour of the combination route, on condition of respecting the TUVRP-TW constraints. The main tour “0-2-1-3-6-7-4-0” is adjusted as “0-2-1-3-6-7-4-8-0,” and the subtour “3-11-5-10-0” is transformed into “3-11-5-10-8.” Figure 7(c) gives an example of the (1, 1)PTR_PUR operator. A TUC “1” in the main tour of the combination route is randomly selected, and a TUC “8” included in the UV route is randomly selected. The two selected TUCs are exchanged with each other.Figure8(a) gives an example of the (1, 0)PUR_SUB operator. A customer “12” in the UV route is randomly selected and assigned as the customer of TUC “2” included in the main tour of the combination route, on condition of respecting the TUVRP-TW constraints. A new subtour “2-12-1” is constructed. The UV route “0-12-8-9-13-0” is transformed into “0-8-9-13-0”. Figure 8(b) gives an example of the (1, 1)PUR_SUB operator. A customer “11” included in a subtour of the combination route is randomly selected, and a customer “8” included in the UV route is randomly selected. The two selected customers are exchanged with each other.Figure9(a) gives an example of the (1, 1)PTR_SUB operator. In the main tour, a TUC “1” that is used as the origin of a subtour is randomly selected, and a TUC “4” included in the subtour is randomly selected. The TUCs “1” and “4” are exchanged with each other. Respecting the TUVRP-TW constraints, the main tour “0-1-2-3-7-8-11-12-0” is transformed into “0-4-2-3-7-8-11-12-0,” and the subtour “1-4-5-3” is transformed into “4-1-5-7.” The TUCs “4” and “7” are renewed as the origin and destination of the subtour, respectively. Figure 9(b) gives an example of the (1, 0)PTR_SUB operator. A TUC “7” in the main tour is randomly selected and inserted into a randomly selected subtour, on the condition that the TUVRP-TW constraints are ensured.


## 4.2.1. The VNS Process

A VNS procedure is used to improve the initial solution (i.e.,initial_S). We set the number of external iterations as Nout and the number of internal iterations as Nin. We define the optimization process as Improving and the disturbance process as Shaking. Within the Shaking process, a poor solution may be accepted, provided that its deviation does not exceed a given threshold (denoted as gap). Nout determines the number of shaking process, and Nin determines the number of internal iterations after the solution perturbation. The VNS combines deterministic and stochastic changes of neighborhood. The deterministic part is represented by a local search heuristic, and the stochastic phase is represented by the random selection of one point in one of the neighborhoods (Villegas et al., 2011).At the beginning of the external iterations,initial_S obtained by the GRASP is regarded as the incumbent solution (denoted as Solution) and the best solution (denoted as best_S). At each external iteration, Solution is perturbed by Shaking. Based on the selected operators of Shaking, Solution is changed, and a solution (denoted as shaking_S), which has the most disturbance within gap, is retained. After the Shaking, Improving is performed, which is iterated Nin times. If the attained result is worse than Shaking_S, Shaking_S is retained; if the attained result is better than Shaking_S, the better solution is retained. After each iteration, it will be judged whether the incumbent solution is optimized. If so, the best solution attained so far will be updated. We denote a finite set of preselected neighborhood structures with Nk1 (k1 = 1, …, kmax1) in shaking and Nk2 (k2 = 1, …, kmax2) in Improving. We introduce Nki (x) (i = 1, 2) to denote the set of solutions in the k1-th or k2-th neighborhood of x, where x is a feasible solution.The process of the whole VNS for solving the TUVRP-TW is shown in Figure4 and is described in pseudo-code in Figure 5.Figure 4
The VNS process.Figure 5
The VNS algorithm.InImproving and Shaking, various types of local search operators are used to search neighborhoods. The objective function of the TUVRP-TW is used as the estimation criterion for substituting the newly generated solution for the incumbent solution. Operators are designed based on relocation operators and interchange operators, which are typically used in the search phase for the VRP with time windows (VRPTW) in the literature. Operators for combination routes probably change the time continuity due to the intersatellite synchronization and the varying subtour origins and destinations.InShaking, seven operators are adopted, that is, a relocation operator for truck routes or main tours of combination routes, a relocation operator for UV routes, a relocation operator for subtours, and four interchange operators for two randomly selected UV routes (denoted as (1, 0)PUR, (1, 1)PUR, (1, 2)PUR, and (2, 0)PUR).Eighteen operators are adopted inImproving, that is, a relocation operator for truck routes and main tours of combination routes, a relocation operator for UV routes, a relocation operator for subtours, four interchange operators for two randomly-selected UV routes (denoted as (1, 0)PUR, (1, 1)PUR, (1, 2)PUR and (2, 0)PUR), four interchange operators for two randomly-selected truck routes or main tours (denoted as (1, 0)PTR, (1, 1)PTR, (1, 2)PTR and (2, 0)PTR), two interchange operators for two randomly-selected truck route (or main tour) and UV route (denoted as (1,0)PTR_PUR, (0,1)PTR_PUR, (1,1)PTR_PUR), two interchange operators for a randomly-selected subtour and a randomly-selected UV route (denoted as (1, 0)PUR_SUB and (1, 1)PUR_SUB), and two interchange operators for a randomly-selected subtour and a randomly-selected truck route(or main tour) (denoted as (1, 0)PTR_SUB and (1, 1)PTR_SUB).


## 4.2.2. Local Search Operators

The operators for truck routes or UV routes are similar to those usually adopted in the search phase for the routing variants involving time windows in the literature (e.g., [35]). The operators involving combination routes are relatively complex because of the intersatellite synchronization and time constraints. We take several routes shown in Figures 6–9 as examples to explain several special operators. In Figures 6–9, TUC is denoted as the box with the serial number of a customer. UC is denoted as the circle with the serial number of a customer.Figure 6
An example of the relocate operator for subtours.Figure 7
Examples of the interchange operators for one main tour (or truck route) and one UV route.
(a)(b)(c)Figure 8
Examples of the interchange operators for one subtour and one UV route.
(a)(b)Figure 9
Examples of the interchange operators for one main tour (or truck route) and one subtour.
(a)(b)Relocation operators aim to relocate a vertex of a randomly selected truck route, UV route, main tour, or subtour of a combination route. Figure6 shows an illustrative example of the relocate operator for subtours. A combination route including the main tour “0-1-2-3-4-5-6-0” and the subtour “2-7-8-6” is randomly selected. A randomly selected vertex (“7”) in the subtour is relocated, which results in a new feasible subtour “2-8-7-5” with the destination of the subtour being adjusted at vertex “5.” It is worth noting that the relocate operator for subtours probably leads to the changing time continuity for the intersatellite synchronization and the subtour destination varying.Interchange operators for two randomly selected routes aim at selecting vertices in each of the two routes for exchange. Figure7(a) gives an example of the (1, 0)PTR_PUR operator. A TUC “2” in the main tour of the combination route is randomly selected and inserted into the UV route, on condition of respecting the TUVRP-TW constraints. Figure 7(b) gives an example of the (0, 1)PTR_PUR operator. A TUC “8” in the UV route is randomly selected and inserted into the main tour of the combination route, on condition of respecting the TUVRP-TW constraints. The main tour “0-2-1-3-6-7-4-0” is adjusted as “0-2-1-3-6-7-4-8-0,” and the subtour “3-11-5-10-0” is transformed into “3-11-5-10-8.” Figure 7(c) gives an example of the (1, 1)PTR_PUR operator. A TUC “1” in the main tour of the combination route is randomly selected, and a TUC “8” included in the UV route is randomly selected. The two selected TUCs are exchanged with each other.Figure8(a) gives an example of the (1, 0)PUR_SUB operator. A customer “12” in the UV route is randomly selected and assigned as the customer of TUC “2” included in the main tour of the combination route, on condition of respecting the TUVRP-TW constraints. A new subtour “2-12-1” is constructed. The UV route “0-12-8-9-13-0” is transformed into “0-8-9-13-0”. Figure 8(b) gives an example of the (1, 1)PUR_SUB operator. A customer “11” included in a subtour of the combination route is randomly selected, and a customer “8” included in the UV route is randomly selected. The two selected customers are exchanged with each other.Figure9(a) gives an example of the (1, 1)PTR_SUB operator. In the main tour, a TUC “1” that is used as the origin of a subtour is randomly selected, and a TUC “4” included in the subtour is randomly selected. The TUCs “1” and “4” are exchanged with each other. Respecting the TUVRP-TW constraints, the main tour “0-1-2-3-7-8-11-12-0” is transformed into “0-4-2-3-7-8-11-12-0,” and the subtour “1-4-5-3” is transformed into “4-1-5-7.” The TUCs “4” and “7” are renewed as the origin and destination of the subtour, respectively. Figure 9(b) gives an example of the (1, 0)PTR_SUB operator. A TUC “7” in the main tour is randomly selected and inserted into a randomly selected subtour, on the condition that the TUVRP-TW constraints are ensured.


## 5. Computational Experiments

We conduct computational experiments to evaluate the performance of the intersatellite synchronization and to evaluate the effectiveness of the TUVRP-TW mathematical formulation and the applicability of the proposed heuristic for large-scale instances. The heuristic that is proposed in Section4 is coded in C++. We use CPLEX 12.9 to solve directly the TUVRP-TW formulation and the TUVRP-TW-B formulation. The codes for the heuristic and the exact algorithm by CPLEX 12.9 are run on a computer with a Windows 10 operating system configured with an Intel(R) Core(TM) i5-10400F CPU @ 2.90-GHz processor with 16.00 GB memory.This section is organized as follows. Section5.1 describes the generated instances. Section 5.2 shows the intersatellite synchronization performance by referring to results of small-scale instances obtained by solving directly the TUVRP-TW formulation and the TUVRP-TW-B formulation. Section 5.3 shows heuristic results of small-scale instances obtained by the algorithm proposed in Section 4. Section 5.4 shows the performance of the TUVRP-TW heuristic for benchmark instances. Sections 5.5 and 5.6 show heuristic results of large-scale instances obtained by the proposed algorithm and the impact of customer time windows, respectively.


### 5.1. Test Instances

We select VRPTW test problems of Solomon [36] and convert them into two types of TUVRP-TW instances. Various scales of instances are generated for two types of computational experiments. In the first computational experiment, we use 21 small-scale instances to evaluate the intersatellite synchronization performance by solving directly the TUVRP-TW formulation and the TUVRP-TW-B formulation. In the second computational experiment, we use 21 small-scale instances and 27 large-scale instances to evaluate the performance of the TUVRP-TW formulation and the applicability of the proposed heuristic.We derive the test instances from VRPTW benchmark problems of Solomon [36] with 100 customers. The Solomon [36] benchmark instances are classified into three categories: the R type, which consists of uniformly distributed customers; the C type, which consists of clustered customers; and the RC type, which is a mix of the R and C types. Each of the categories consists of two sets: sets R1, C1, and RC1 have narrow scheduling horizons, whereas sets R2, C2, and RC2 have wide scheduling horizons. We select the nine instances with 100 customers from set C1 for generating test instances for the two types of computational experiments.Considering the computational experiment requirements for the TUVRP-TW and respecting the customer locations and the demand with time windows of the Solomon [36] benchmark instances, we modify the VRPTW benchmark instances by adjusting the number of customers and classifying the customers into TUCs and UCs. Each of the instances that has the same number of customers is converted into three types of test instances via a method that is similar to the method that Chao [37] used to generate the TTRP instances. In a VRPTW benchmark instance, the distance between each customer i and its nearest neighboring customer is calculated and denoted by Ai. In the first type of TUVRP-TW instances, 25% of the customers with the smallest Ai values are specified as UCs. This percentage is increased to 50% or 75% in the second type or the third type of TUVRP-TW instances, respectively.The 21 small-scale instances that are used to evaluate the intersatellite synchronization performance are generated by selecting randomly a specified number of customers from the Solomon [36] instances with 100 customers from set C1. In an instance, the number of customers (denoted as cn) varies from 6 to 12 with an increment of 1. For instances that have the same number of customers, the UC percentage varies from 25% to 75% with an increment of 25%. The small-scale instances are classified according to cn and the UC percentage. Each small-scale instance is denoted by C1-cn-cp, where cp is the UC percentage. The 27 large-scale instances are generated by classifying 100 customers from set C1 according to the customer types of the Solomon [36] instances. Each instance with 100 customers from set C1 is converted into three types of large-scale instances via an approach that is similar to the method that was used to generate the small-scale instances. Each large-scale instance is denoted by “Solomon (1987) instance No.”-cp. The parameters involved in the TUVRP-TW formulation and algorithm are evaluated in Table 2.Table 2
Values of involved parameters.
ParametersValue (unit)Small-scale instanceLarge-scale instanceQt100 (kg)100 (kg)Qu100 (kg)100 (kg)vt50 (km/h)50 (km/h)vu50 (km/h)50 (km/h)τ11Nout200500Nin10002000gap3070


### 5.2. Intersatellite Synchronization Performance

Computational experiments on small-scale instances are conducted to evaluate the performance of intersatellite synchronization via a comparison of exact solutions for small-scale instances. For each small-scale instance, the total computation time for CPLEX 12.9 is set to 14400 s. CPLEX 12.9 runs with default settings until it has found an exact solution or the computation time has been exhausted. The computational results of the 21 small-scale instances are presented in Table3. Column 1 lists the test instance. Columns 2–4 list the exact results of the TUVRP-TW formulation. Columns 2 and 3 list the results of the objective (ObjE) and the numbers of main tours, subtours, and UV routes (denoted as “m/s/U”). Column 4 lists the computation time (seconds). Columns 5 and 7 list the exact results of the TUVRP-TW-B formulation. Columns 5 and 6 list the results of the objective (Obj) and the numbers of main tours, subtours, and UV routes (denoted as “m/s/U”). Column 7 lists the computation time (seconds).Table 3
Exact solutions for small-scale instances.
InstanceTUVRP-TW formulationTUVRP-TW-B formulationObjEm/s/UTime (s)Objm/s/UTime (s)C1-6-2546.20/0/13546.20/0/14C1-6-5046.20/0/1546.20/0/1<1C1-6-7546.20/0/1246.20/0/1<1C1-7-2570.61/1/067382.20/0/2214C1-7-5072.71/1/03582.20/0/28C1-7-7572.71/1/01182.20/0/23C1-8-2567.21/1/027681.51/0/160C1-8-5067.21/1/03381.50/0/25C1-8-7581.50/0/2481.50/0/2<1C1-9-25115.91/1/0867126.20/0/2143C1-9-50126.20/0/2138126.20/0/27C1-9-75126.20/0/211126.20/0/23C1-10-25106.41/1/014400∗121.91/0/13967C1-10-50121.90/0/214400∗121.90/0/2654C1-10-75121.90/0/2195121.90/0/239C1-11-25207.30/0/314400∗207.30/0/32234C1-11-50207.30/0/310413207.30/0/3304C1-11-75207.30/0/3670207.30/0/357C1-12-25208.91/0/214400∗208.91/0/21656C1-12-50221.80/0/41036221.80/0/493C1-12-75221.80/0/437221.80/0/410Note. ∗ denotes that the computation time is exhausted.The exact solutions for small-scale instances are affected by the instance data. We summarize several conclusions as follows.(i)
Intersatellite synchronization enables more flexible delivery operations of vehicles, which reduces the integrated distance in the TUVRP-TW model. Of the 17 instances with exact TUVRP-TW solutions, 6 instances have smaller objective values in the TUVRP-TW formulation than in the TUVRP-TW-B formulation. The other 11 instances have the same objective values in the TUVRP-TW formulation and in the TUVRP-TW-B formulation. Two instances, namely, C1-8-25 and C1-9-25, are selected as examples for demonstrating flexible delivery operations of vehicles. For instance C1-8-25, the optimal solution of the TUVRP-TW formulation includes one main tour, namely, “0-4-0,” and one subtour, namely, “4-7-5-6-8-3-1-2-0,” while the optimal solution of the TUVRP-TW-B formulation includes one truck route, namely, “0-2-0,” and one UV route, namely, “0-4-7-5-6-8-3-1-0.” For instance C1-9-25, the optimal solution of the TUVRP-TW formulation includes one main tour, namely, “0-7-6-3-0,” and one subtour, namely, “3-8-4-5-9-2-1-0,” while the optimal solution of the TUVRP-TW-B formulation includes two UV routes, namely, “0-7-6-8-3-0” and “0-4-5-9-2-1-0.” Hence, the combination of the main tour and the subtour reduces the integrated distance, compared with pure UV routes or pure truck routes.(ii)
A lower UC percentage for an instance corresponds to more TUCs being included in the instance. The inclusion of more TUCs in an instance corresponds to the availability of more candidate satellites for constructing subtours, which leads to a higher computational workload in the solution process. For instances that have the same number of customers, the computation time varies substantially with the UC percentage of the instance.(iii)
The results of the computational experiments on small-scale instances demonstrate that the TUVRP-TW formulation can be directly solved by CPLEX 12.9. However, CPLEX 12.9 fails to solve directly some small-scale instances that include more customers. Various strategies (e.g., the decomposition method) for improving the performance of the exact method by CPLEX 12.9 may be useful. However, instances are seldom solved directly by CPLEX 12.9 for large-scale instances because too much computer memory and computation time are required. It is necessary to evaluate the performance of the heuristic.


### 5.3. Heuristic Results of Small-Scale Instances

The performance of the proposed heuristic, which is measured in terms of the quality of the heuristic solutions, is assessed by comparing heuristic solutions with exact solutions for small-scale instances. Table4 lists the heuristic results for the 21 small-scale instances, which include the results of the initial solution (columns 2 and 3) and the best solution that were obtained by the heuristic (columns 4 and 5). Column 1 lists the instance. Column 2 lists the objective value (denoted as ObjI) of the initial solution. Column 3 lists the numbers of main tours, subtours, and UV routes (denoted as “m/s/U”) in the initial solution. Column 4 lists the objective value (denoted as ObjO) of the heuristic. Column 5 lists the numbers of main tours, subtours, and UV routes in the heuristic solution. Column 6 lists the computation time (seconds). Column 7 lists the percentage gap between the objective value of the initial solution (ObjI) and the objective value of the heuristic solution (ObjO), which is estimated as Gap1 = 100% × (ObjI − ObjO)/ObjI. Column 8 lists the percentage gap between the objective value of the heuristic solution and the objective value of the exact solution (denoted as ObjE), which is estimated as Gap2 = 100% × (ObjO − ObjE)/ObjO.Table 4
Heuristic solutions for 21 small-scale instances.
InstanceInitial solutionOptimized solutionTime (s)Gap1 (%)Gap2 (%)Objlm/s/UObjOm/s/UC1-6-2546.30/0/146.20/0/11.60.20.0C1-6-5046.30/0/146.20/0/11.50.20.0C1-6-7546.20/0/146.20/0/12.40.00.0C1-7-2585.50/0/270.61/1/08.317.40.0C1-7-5086.60/0/272.71/1/04.816.10.0C1-7-7577.61/1/072.71/1/05.36.30.0C1-8-2567.71/1/067.21/1/01.30.70.0C1-8-5086.00/0/267.21/1/07.621.90.0C1-8-75118.70/0/281.50/0/20.731.30.0C1-9-25214.31/1/2115.91/1/06.145.90.0C1-9-50191.91/0/2126.20/0/22.634.20.0C1-9-75165.70/0/2126.20/0/22.423.80.0C1-10-25173.11/0/2106.41/1/09.838.50.0#C1-10-50148.01/1/1121.90/0/22.917.60.0#C1-10-75150.01/1/1121.90/0/23.218.70.0C1-11-25276.11/0/3207.30/0/35.524.90.0#C1-11-50307.31/0/3207.30/0/33.232.50.0C1-11-75247.21/0/2207.30/0/33.116.10.0C1-12-25265.10/0/4208.91/0/28.221.20.0#C1-12-50316.20/0/4221.80/0/43.629.90.0C1-12-75222.100/4221.80/0/44.50.10.0Note. # denotes that Gap2 is estimated by referring to a solution obtained by CPLEX 12.9 in the limited computation time.According to the computational results, the proposed heuristic has a general and adaptable structure. For the 21 small-scale instances, the heuristic does not require instance-related iteration settings.According to the computation time and the Gap2 values for the small-scale instances, the performance of the heuristic is satisfactory. Compared with the exact method by CPLEX 12.9, the heuristic can output the solution quickly, and the heuristic is relatively stable. Based on the heuristic running time of the heuristic solution, the heuristic can solve each of the 21 small-scale instances in less than 10 s, and the average computation time is 4.2 s. In terms of the solution quality, the heuristic performs well. Of the 17 small-scale instances with exact solutions, Gap2 is 0.0% in all the 17 cases; hence, these instances can be solved optimally via the proposed heuristic. The gap (Gap1) between the objective value of the initial solution and that of the heuristic solution is relatively large. Of all the 21 small-scale instances, Gap1 exceeds 20.0% in 10 cases. The largest value of Gap1 is 45.9% and the average value of Gap1 is 18.9%. Gap1 is less than 1.0% in 5 cases. Although the GRASP may obtain a satisfactory solution in some cases, the VNS is necessary and important for the proposed heuristic.


### 5.4. Performance of the TUVRP-TW Heuristic for Benchmark Instances

The TUVRP-TW is fundamentally a new variant. From the problem definition perspective, the TUVRP-TW situation is similar to the situation of the TTRP with time windows (TTRPTW). The combination vehicle used in the TUVRP-TW involves pairing one UV with one truck. Both trucks and UVs are autonomous and can serve assigned customers. The TTRPTW involves pairing one trailer with one truck, while a trailer is nonautonomous and cannot serve customers without its truck pulling.We adopt the TUVRP-TW heuristic to solve the TTRPTW benchmark instances. Because vehicles (i.e., truck-trailer combination and truck) involved in the TTRPTW are obviously different from vehicles (i.e., truck with/without carried UV and UV) involved in the TUVRP-TW, we assume that truck and UV have the same capacities and velocities as those of truck-trailer combination and truck and adjust the objective function of the TUVRP-TW. The results attained by the proposed heuristic and the known results of the TTRPTW benchmark instances are shown in Table5. Column 1 shows the instances. Column 2 shows the best known solution value (BKS), which is attained in Lin et al. [39]. Columns 3 and 4 show the best results reported by Parragh and Cordeau [38] and their deviations from the BKS. Columns 5 and 6 show the best objective values encountered in Parragh and Cordeau [38] during all parameter tuning tests and their respective deviations from the best known results. Columns 7 and 8 show the best results reported by the TUVRP-TW heuristic and their deviations from the BKS. We introduce Gap5 to denote the percentage gap between columns 3 and 7. We introduce Gap6 to denote the percentage gap between columns 5 and 7.Table 5
The TUVRP-TW heuristic results for TTRPTW benchmark instances.
InstanceBKSParragh and Cordeau [38]TUVRP-TW algorithm in the paperBest (of 5)Dev (%)Best seenDev (%)BestDev (%)Gap5 (%)Gap6 (%)C101-25971.82905.02−6.87903.84−7.001015.214.512.212.3C101-501106.901004.39−9.261004.39−9.261120.411.211.611.6C101-751154.801068.92−7.441058.85−8.311179.842.210.411.4C201-25671.31671.370.01671.370.01713.606.36.36.3C201-50713.06709.07−0.56709.07−0.56711.46−0.20.30.3C201-75711.45711.460.00711.460.00711.460.00.00.0R101-251651.161644.64−0.391644.64−0.391633.37−1.1−0.7−0.7R101-501644.641644.640.001644.640.001637.06−0.5−0.5−0.5R101-751644.641644.640.001644.640.001642.81−0.1−0.1−0.1R201-251166.641152.25−1.231147.80−1.611148.20−1.6−0.40.0R201-501148.171149.790.141147.80−0.031154.790.60.40.6R201-751147.801147.880.011147.800.001156.310.70.70.7RC101-251717.011654.39−3.651654.39−3.651708.13−0.53.23.2RC101-501770.401724.23−2.611710.32−3.391759.35−0.62.02.9RC101-751784.061779.61−0.251779.45−0.261794.560.60.80.8RC201-251277.051265.56−0.901265.56−0.901284.390.61.51.5RC201-501273.691266.11−0.601265.56−0.641274.420.10.70.7RC201-751266.111266.110.001265.98−0.011272.790.50.50.5Average1267.821245.001243.201273.23To the best of our knowledge, the best results of TTRPTW benchmark instances are reported by Parragh and Cordeau [38]. Parragh and Cordeau [38] proposed a branch-and-price algorithm for the TTRPTW, using problem specific enhancements in the pricing scheme and alternative lower bound computations. They tailored an adaptive large neighborhood search algorithm to the TTRPTW in order to obtain good initial columns. We introduce Gap5 and Gap6 to compare results attained by the TUVRP-TW heuristic with the best results in the literature. For the 18 benchmark instances, the TUVRP-TW heuristic can attain new best solutions for four benchmark instances and the same solution for one benchmark instance. The smallest Gap5 is −0.7%, and the average Gap5 is 2.7%. The smallest Gap6 is −0.7%, and the average Gap6 is 2.9%. Generally speaking, the TUVRP-TW heuristic can obtain average gaps to best known solutions of less than 2.9%. We investigate the detailed results attained by the TUVRP-TW heuristic and find that instance data affect obviously the attained results. Some TTRPTW instances (e.g., C101-25, C101-50, and C101-75) may not be suitable for the TUVRP-TW situation.


### 5.5. Heuristic Applicability for Large-Scale Instances

The conclusions that are obtained by comparing heuristic solutions with exact solutions for small-scale instances and the TUVRP-TW heuristic results for TTRPTW benchmark instances support the satisfactory performance of the proposed heuristic. Several large-scale instances are solved by the heuristic in the computational experiments. Table6 lists the heuristic results for the 27 large-scale instances. For each large-scale instance in the computational experiments, the heuristic is repeated ten times. The objective value, the number of tours or routes, and the heuristic running time of the best solution among the ten repetitions are reported. Column 1 lists the instance. Column 2 lists the objective value of the initial solution. Column 3 lists the numbers of main tours, subtours, and UV routes in the initial solution. Column 4 lists the objective value of the heuristic. Column 5 lists the numbers of main tours, subtours, and UV routes in the heuristic solution. Column 6 lists the computation time for finding the best solution. Column 7 lists the percentage gap between the objective value of the initial solution and the objective value of the heuristic solution, namely, Gap1.Table 6
Heuristic solutions of 27 large-scale instances.
InstanceInitial solutionOptimized solutionTime (s)Gap1 (%)ObjIm/s/UObjOm/s/UC101-253642.46/0/142802.83/1/162533.423.0C101-504020.45/2/182807.30/0/221211.830.2C101-753857.45/3/192846.10/0/201139.526.2C102-253110.06/1/122235.55/1/112309.128.1C102-503438.34/2/162236.01/0/181385.135.0C102-753472.06/2/192303.50/0/201342.333.7C103-252913.36/3/111638.36/1/72463.143.8C103-502807.04/0/151763.93/1/141715.637.2C103-752975.74/0/171851.10/0/191217.337.8C104-252576.55/2/111250.07/2/51557.251.5C104-502705.23/5/141424.94/2/111787.247.3C104-752721.62/1/161447.65/1/151433.046.8C105-253048.58/4/92067.33/1/142965.732.2C105-503438.45/2/122142.13/1/142740.537.7C105-753518.23/3/152027.60/0/191193.942.4C106-253268.63/0/151904.14/1/122070.441.7C106-503329.75/3/101983.33/3/161770.140.4C106-753352.63/1/172078.81/1/171775.838.0C107-252890.74/1/111771.85/3/82674.238.7C107-502895.64/1/111770.71/1/182280.638.8C107-753075.05/3/141833.11/1/181344.240.4C108-252646.86/1/81511.76/1/82664.542.9C108-503227.93/0/141622.83/1/132023.349.7C108-752951.63/2/161655.41/1/171399.543.9C109-252454.66/2/71242.34/2/83179.649.4C109-502726.93/0/131448.32/1/151790.146.9C109-752797.03/2/141596.82/2/151856.242.9The computation time of the heuristic for each of the 27 large-scale instance does not exceed 0.9 h, and the average computation time is 0.5 h. Hence, the proposed heuristic can solve one large-scale instance in an acceptable time for practices.The initial solution that is constructed by the GRASP is feasible for enterprise operations; however, an initial solution that is constructed using the GRASP is likely to fall into a local optimum. The gap (Gap1) between the objective value of the initial solution and that of the heuristic solution is large. For the 27 large-scale instances, Gap1 is more than 40.0% in 14 cases. The largest value of Gap1 is 51.5%, and the average value of Gap1 is 39.5%. The VNS is highly important for finding a locally optimal solution. For the 27 large-scale instances, the average value of Gap1 is larger than the average value of Gap1 for the 21 small-scale instances. The scale of each large-scale instance is larger than that of any small-scale instance, and the initial solution of a large-scale instance may have more neighborhoods than a small-scale instance.Comparing the initial solution with the final solution that is provided by the heuristic, various optimization strategies can be identified. For example, the number of main tours or subtours of the initial solution and the number of main tours or subtours of the final solution differ substantially. The number of main tours of the initial solution, the number subtours of the initial solution, the number of main tours of the final solution, and the number subtours of the final solution are denoted as NMI, NSI, NMO, and NSO, respectively. For 17 of the 27 large-scale instances, the percentage gap between NMI and NMO, namely, Gap3 = 100% × (NMl − NMO)/NMI, exceeds 33.0%. The average value of Gap3 equals 33.9%. For 11 of the 27 large-scale instances, the percentage gap between NSI and NSO, namely, Gap4 = 100% × (NSI − NSO)/NSI, exceeds 31.0%. The average value of Gap4 equals 31.8%. Hence, additional savings may be realized by optimizing the main tours and subtours. For enterprises, focusing on optimizing the main tours and subtours is more desirable.The results of the computational trials on the 27 large-scale instances demonstrate that the instance data and the estimated values of various parameters (e.g., the vehicle velocity and the vehicle capacity) substantially affect the heuristic solutions. We summarize the management insights as follows, which are obtained by respecting the limitations of the instances and the computational experiments. (i) Catering to the “mothership approach” in the expected delivery practice, the TUVRP-TW includes intersatellite synchronization, which is expected to result in more efficient and more flexible delivery operations. Aiming at decreasing the integrated distance, it is not always suitable to adopt the “mothership approach.” Of the 27 large-scale instances, six instances have no subtours in the heuristic solutions. The heuristic results of all 27 large-scale instances include 73 main tours and 29 subtours; hence, the approach of combining the main tour with one or several subtours is not always widely adopted in these instances. Considering the comparison of the exact solutions for small-scale instances that are obtained by directly solving the TUVRP-TW formulation and the TUVRP-TW-B formulation, we find that the “mothership approach” is suitable for pursuing savings in various practical scenarios. Both practitioners and researchers are suggested to identify suitable practical scenarios to employ truck-UV combinations for delivery. (ii) Using UVs for direct delivery from the depot seems necessary for reducing the integrated distance. For each instance, the number of UV routes exceeds the number of main tours or the number of subtours. Hence, the direct delivery strategy is necessary and important in implementing the “mothership approach.”


### 5.6. Impact of Customer Time Windows

In the computational experiments on large-scale instances, sensitivity analysis is conducted to evaluate the impact of time windows on the large-scale instances. By changing the width of customer time windows, the large-scale instances designed in Section5.1, in which the time window of customer i is ei,li and the width of the time window of customer i is denoted as witw (witw=li−ei), are converted to another two types of instances. First, each large-scale instance is denoted by Cbr10J-p. The only difference between instances Cbr10J-p and C10J-p is the width of customer time windows. The time window of customer i in instances Cbr10J-p is ei−0.5⋅witw,li+0.5⋅witw. Second, each large-scale instance is denoted by Cna10J-p. The only difference between instances Cna10J-p and C10J-p is the width of customer time windows. The time window of customer i in instances Cna10J-p is ei+0.25⋅witw,li−0.25⋅witw. The proposed heuristic is used to solve instances Cbr10J-p and Cna10J-p. For each instance, the proposed heuristic is run 10 times, and the objective value and the computation time of the best solution of the 10 results are reported in Table 7.Table 7
Heuristic results of large-scale instances with changed time windows.
InstanceObjOTime (s)Gap3InstanceObjOTime (s)Gap4Cbr101−252101.42814.0−25.0Cna101−253300.92357.017.8Cbr101−502144.51358.9−23.6Cna101−503380.71233.820.4Cbr101−752205.91215.8−22.5Cna101−753520.31295.423.7Cbr102−251750.82663.1−21.7Cna102−252663.52737.619.1Cbr102−501742.91808.0−22.1Cna102−502949.72406.731.9Cbr102−751796.21660.3−22.0Cna102−752743.81301.419.1Cbr103−251445.92940.9−11.7Cna103−252002.02138.322.2Cbr103−501575.32221.5−10.7Cna103−502069.52260.417.3Cbr103−751657.51246.5−10.5Cna103−752200.01264.618.9Cbr104−251193.92336.4−4.5Cna104−251398.72187.911.9Cbr104−501356.32380.4−4.8Cna104−501556.73175.79.2Cbr104−751403.71370.1−3.0Cna104−751707.01414.117.9Cbr105−251520.52763.8−26.5Cna105−252774.72420.734.2Cbr105−501570.92040.4−26.7Cna105−502901.92318.635.5Cbr105−751685.12322.3−16.9Cna105−752761.61591.436.2Cbr106−251428.82444.8−25.0Cna106−252462.31722.329.3Cbr106−501604.62260.6−19.1Cna106−502050.31222.63.4Cbr106−751610.91430.5−22.5Cna106−752587.11342.124.4Cbr107−251300.13063.9−26.6Cna107−252158.52688.621.8Cbr107−501474.52073.7−16.7Cna107−502291.41684.729.4Cbr107−751561.71294.7−14.8Cna107−752230.11277.721.7Cbr108−251160.33752.1−23.2Cna108−251923.22614.627.2Cbr108−501398.33265.9−13.8Cna108−502201.82195.135.7Cbr108−751426.41468.5−13.8Cna108−752033.01233.022.8Cbr109−251018.43337.0−18.0Cna109−251541.72825.224.1Cbr109−501191.32453.6−17.7Cna109−501645.31811.613.6Cbr109−751353.61807.3−15.2Cna109−751723.51297.07.9To compare the objective values of solutions of instances Cbr10J-p and C10J-p, we introduce a percentage gap, Gap3, where Gap3 = 100% × (the objective value of solution of instance “Cbr10J-p” − the objective value of solution of instance “C10J-p”)/the objective value of solution of instance “C10J-p.” To compare the objective values of solutions of instances Cna10J-p and C10J-p, we introduce a percentage gap, Gap4, where Gap4 = 100% × (the objective value of solution of instance “Cna10J-p” − the objective value of solution of instance “C10J-p”)/the objective value of solution of instance “C10J-p.” The average, largest, and the smallest Gap3 values are −17.7%, −3.0%, and −26.7%, respectively. The average, largest, and the smallest Gap4 values are 22.1%, 36.2%, and 3.4%, respectively. It indicates that wide time windows may lead to the decrease of the integrated distance while narrow time windows may lead to the increase of the integrated distance. In terms of the heuristic computation time, the proposed heuristic is able to solve each of the 27 Cbr10J-p instances in an average computational time of 2215 s and is able to solve each of the 27 Cna10J-p instances in an average computational time of 1927 s. Wide time windows may enlarge the neighborhood, so that more computation time is spent by the heuristic to find satisfactory solutions. Narrow time windows may cut the neighborhood, so that less computation time is spent by the heuristic to find satisfactory solutions.


## 5.1. Test Instances

We select VRPTW test problems of Solomon [36] and convert them into two types of TUVRP-TW instances. Various scales of instances are generated for two types of computational experiments. In the first computational experiment, we use 21 small-scale instances to evaluate the intersatellite synchronization performance by solving directly the TUVRP-TW formulation and the TUVRP-TW-B formulation. In the second computational experiment, we use 21 small-scale instances and 27 large-scale instances to evaluate the performance of the TUVRP-TW formulation and the applicability of the proposed heuristic.We derive the test instances from VRPTW benchmark problems of Solomon [36] with 100 customers. The Solomon [36] benchmark instances are classified into three categories: the R type, which consists of uniformly distributed customers; the C type, which consists of clustered customers; and the RC type, which is a mix of the R and C types. Each of the categories consists of two sets: sets R1, C1, and RC1 have narrow scheduling horizons, whereas sets R2, C2, and RC2 have wide scheduling horizons. We select the nine instances with 100 customers from set C1 for generating test instances for the two types of computational experiments.Considering the computational experiment requirements for the TUVRP-TW and respecting the customer locations and the demand with time windows of the Solomon [36] benchmark instances, we modify the VRPTW benchmark instances by adjusting the number of customers and classifying the customers into TUCs and UCs. Each of the instances that has the same number of customers is converted into three types of test instances via a method that is similar to the method that Chao [37] used to generate the TTRP instances. In a VRPTW benchmark instance, the distance between each customer i and its nearest neighboring customer is calculated and denoted by Ai. In the first type of TUVRP-TW instances, 25% of the customers with the smallest Ai values are specified as UCs. This percentage is increased to 50% or 75% in the second type or the third type of TUVRP-TW instances, respectively.The 21 small-scale instances that are used to evaluate the intersatellite synchronization performance are generated by selecting randomly a specified number of customers from the Solomon [36] instances with 100 customers from set C1. In an instance, the number of customers (denoted as cn) varies from 6 to 12 with an increment of 1. For instances that have the same number of customers, the UC percentage varies from 25% to 75% with an increment of 25%. The small-scale instances are classified according to cn and the UC percentage. Each small-scale instance is denoted by C1-cn-cp, where cp is the UC percentage. The 27 large-scale instances are generated by classifying 100 customers from set C1 according to the customer types of the Solomon [36] instances. Each instance with 100 customers from set C1 is converted into three types of large-scale instances via an approach that is similar to the method that was used to generate the small-scale instances. Each large-scale instance is denoted by “Solomon (1987) instance No.”-cp. The parameters involved in the TUVRP-TW formulation and algorithm are evaluated in Table 2.Table 2
Values of involved parameters.
ParametersValue (unit)Small-scale instanceLarge-scale instanceQt100 (kg)100 (kg)Qu100 (kg)100 (kg)vt50 (km/h)50 (km/h)vu50 (km/h)50 (km/h)τ11Nout200500Nin10002000gap3070


## 5.2. Intersatellite Synchronization Performance

Computational experiments on small-scale instances are conducted to evaluate the performance of intersatellite synchronization via a comparison of exact solutions for small-scale instances. For each small-scale instance, the total computation time for CPLEX 12.9 is set to 14400 s. CPLEX 12.9 runs with default settings until it has found an exact solution or the computation time has been exhausted. The computational results of the 21 small-scale instances are presented in Table3. Column 1 lists the test instance. Columns 2–4 list the exact results of the TUVRP-TW formulation. Columns 2 and 3 list the results of the objective (ObjE) and the numbers of main tours, subtours, and UV routes (denoted as “m/s/U”). Column 4 lists the computation time (seconds). Columns 5 and 7 list the exact results of the TUVRP-TW-B formulation. Columns 5 and 6 list the results of the objective (Obj) and the numbers of main tours, subtours, and UV routes (denoted as “m/s/U”). Column 7 lists the computation time (seconds).Table 3
Exact solutions for small-scale instances.
InstanceTUVRP-TW formulationTUVRP-TW-B formulationObjEm/s/UTime (s)Objm/s/UTime (s)C1-6-2546.20/0/13546.20/0/14C1-6-5046.20/0/1546.20/0/1<1C1-6-7546.20/0/1246.20/0/1<1C1-7-2570.61/1/067382.20/0/2214C1-7-5072.71/1/03582.20/0/28C1-7-7572.71/1/01182.20/0/23C1-8-2567.21/1/027681.51/0/160C1-8-5067.21/1/03381.50/0/25C1-8-7581.50/0/2481.50/0/2<1C1-9-25115.91/1/0867126.20/0/2143C1-9-50126.20/0/2138126.20/0/27C1-9-75126.20/0/211126.20/0/23C1-10-25106.41/1/014400∗121.91/0/13967C1-10-50121.90/0/214400∗121.90/0/2654C1-10-75121.90/0/2195121.90/0/239C1-11-25207.30/0/314400∗207.30/0/32234C1-11-50207.30/0/310413207.30/0/3304C1-11-75207.30/0/3670207.30/0/357C1-12-25208.91/0/214400∗208.91/0/21656C1-12-50221.80/0/41036221.80/0/493C1-12-75221.80/0/437221.80/0/410Note. ∗ denotes that the computation time is exhausted.The exact solutions for small-scale instances are affected by the instance data. We summarize several conclusions as follows.(i)
Intersatellite synchronization enables more flexible delivery operations of vehicles, which reduces the integrated distance in the TUVRP-TW model. Of the 17 instances with exact TUVRP-TW solutions, 6 instances have smaller objective values in the TUVRP-TW formulation than in the TUVRP-TW-B formulation. The other 11 instances have the same objective values in the TUVRP-TW formulation and in the TUVRP-TW-B formulation. Two instances, namely, C1-8-25 and C1-9-25, are selected as examples for demonstrating flexible delivery operations of vehicles. For instance C1-8-25, the optimal solution of the TUVRP-TW formulation includes one main tour, namely, “0-4-0,” and one subtour, namely, “4-7-5-6-8-3-1-2-0,” while the optimal solution of the TUVRP-TW-B formulation includes one truck route, namely, “0-2-0,” and one UV route, namely, “0-4-7-5-6-8-3-1-0.” For instance C1-9-25, the optimal solution of the TUVRP-TW formulation includes one main tour, namely, “0-7-6-3-0,” and one subtour, namely, “3-8-4-5-9-2-1-0,” while the optimal solution of the TUVRP-TW-B formulation includes two UV routes, namely, “0-7-6-8-3-0” and “0-4-5-9-2-1-0.” Hence, the combination of the main tour and the subtour reduces the integrated distance, compared with pure UV routes or pure truck routes.(ii)
A lower UC percentage for an instance corresponds to more TUCs being included in the instance. The inclusion of more TUCs in an instance corresponds to the availability of more candidate satellites for constructing subtours, which leads to a higher computational workload in the solution process. For instances that have the same number of customers, the computation time varies substantially with the UC percentage of the instance.(iii)
The results of the computational experiments on small-scale instances demonstrate that the TUVRP-TW formulation can be directly solved by CPLEX 12.9. However, CPLEX 12.9 fails to solve directly some small-scale instances that include more customers. Various strategies (e.g., the decomposition method) for improving the performance of the exact method by CPLEX 12.9 may be useful. However, instances are seldom solved directly by CPLEX 12.9 for large-scale instances because too much computer memory and computation time are required. It is necessary to evaluate the performance of the heuristic.


## 5.3. Heuristic Results of Small-Scale Instances

The performance of the proposed heuristic, which is measured in terms of the quality of the heuristic solutions, is assessed by comparing heuristic solutions with exact solutions for small-scale instances. Table4 lists the heuristic results for the 21 small-scale instances, which include the results of the initial solution (columns 2 and 3) and the best solution that were obtained by the heuristic (columns 4 and 5). Column 1 lists the instance. Column 2 lists the objective value (denoted as ObjI) of the initial solution. Column 3 lists the numbers of main tours, subtours, and UV routes (denoted as “m/s/U”) in the initial solution. Column 4 lists the objective value (denoted as ObjO) of the heuristic. Column 5 lists the numbers of main tours, subtours, and UV routes in the heuristic solution. Column 6 lists the computation time (seconds). Column 7 lists the percentage gap between the objective value of the initial solution (ObjI) and the objective value of the heuristic solution (ObjO), which is estimated as Gap1 = 100% × (ObjI − ObjO)/ObjI. Column 8 lists the percentage gap between the objective value of the heuristic solution and the objective value of the exact solution (denoted as ObjE), which is estimated as Gap2 = 100% × (ObjO − ObjE)/ObjO.Table 4
Heuristic solutions for 21 small-scale instances.
InstanceInitial solutionOptimized solutionTime (s)Gap1 (%)Gap2 (%)Objlm/s/UObjOm/s/UC1-6-2546.30/0/146.20/0/11.60.20.0C1-6-5046.30/0/146.20/0/11.50.20.0C1-6-7546.20/0/146.20/0/12.40.00.0C1-7-2585.50/0/270.61/1/08.317.40.0C1-7-5086.60/0/272.71/1/04.816.10.0C1-7-7577.61/1/072.71/1/05.36.30.0C1-8-2567.71/1/067.21/1/01.30.70.0C1-8-5086.00/0/267.21/1/07.621.90.0C1-8-75118.70/0/281.50/0/20.731.30.0C1-9-25214.31/1/2115.91/1/06.145.90.0C1-9-50191.91/0/2126.20/0/22.634.20.0C1-9-75165.70/0/2126.20/0/22.423.80.0C1-10-25173.11/0/2106.41/1/09.838.50.0#C1-10-50148.01/1/1121.90/0/22.917.60.0#C1-10-75150.01/1/1121.90/0/23.218.70.0C1-11-25276.11/0/3207.30/0/35.524.90.0#C1-11-50307.31/0/3207.30/0/33.232.50.0C1-11-75247.21/0/2207.30/0/33.116.10.0C1-12-25265.10/0/4208.91/0/28.221.20.0#C1-12-50316.20/0/4221.80/0/43.629.90.0C1-12-75222.100/4221.80/0/44.50.10.0Note. # denotes that Gap2 is estimated by referring to a solution obtained by CPLEX 12.9 in the limited computation time.According to the computational results, the proposed heuristic has a general and adaptable structure. For the 21 small-scale instances, the heuristic does not require instance-related iteration settings.According to the computation time and the Gap2 values for the small-scale instances, the performance of the heuristic is satisfactory. Compared with the exact method by CPLEX 12.9, the heuristic can output the solution quickly, and the heuristic is relatively stable. Based on the heuristic running time of the heuristic solution, the heuristic can solve each of the 21 small-scale instances in less than 10 s, and the average computation time is 4.2 s. In terms of the solution quality, the heuristic performs well. Of the 17 small-scale instances with exact solutions, Gap2 is 0.0% in all the 17 cases; hence, these instances can be solved optimally via the proposed heuristic. The gap (Gap1) between the objective value of the initial solution and that of the heuristic solution is relatively large. Of all the 21 small-scale instances, Gap1 exceeds 20.0% in 10 cases. The largest value of Gap1 is 45.9% and the average value of Gap1 is 18.9%. Gap1 is less than 1.0% in 5 cases. Although the GRASP may obtain a satisfactory solution in some cases, the VNS is necessary and important for the proposed heuristic.


## 5.4. Performance of the TUVRP-TW Heuristic for Benchmark Instances

The TUVRP-TW is fundamentally a new variant. From the problem definition perspective, the TUVRP-TW situation is similar to the situation of the TTRP with time windows (TTRPTW). The combination vehicle used in the TUVRP-TW involves pairing one UV with one truck. Both trucks and UVs are autonomous and can serve assigned customers. The TTRPTW involves pairing one trailer with one truck, while a trailer is nonautonomous and cannot serve customers without its truck pulling.We adopt the TUVRP-TW heuristic to solve the TTRPTW benchmark instances. Because vehicles (i.e., truck-trailer combination and truck) involved in the TTRPTW are obviously different from vehicles (i.e., truck with/without carried UV and UV) involved in the TUVRP-TW, we assume that truck and UV have the same capacities and velocities as those of truck-trailer combination and truck and adjust the objective function of the TUVRP-TW. The results attained by the proposed heuristic and the known results of the TTRPTW benchmark instances are shown in Table5. Column 1 shows the instances. Column 2 shows the best known solution value (BKS), which is attained in Lin et al. [39]. Columns 3 and 4 show the best results reported by Parragh and Cordeau [38] and their deviations from the BKS. Columns 5 and 6 show the best objective values encountered in Parragh and Cordeau [38] during all parameter tuning tests and their respective deviations from the best known results. Columns 7 and 8 show the best results reported by the TUVRP-TW heuristic and their deviations from the BKS. We introduce Gap5 to denote the percentage gap between columns 3 and 7. We introduce Gap6 to denote the percentage gap between columns 5 and 7.Table 5
The TUVRP-TW heuristic results for TTRPTW benchmark instances.
InstanceBKSParragh and Cordeau [38]TUVRP-TW algorithm in the paperBest (of 5)Dev (%)Best seenDev (%)BestDev (%)Gap5 (%)Gap6 (%)C101-25971.82905.02−6.87903.84−7.001015.214.512.212.3C101-501106.901004.39−9.261004.39−9.261120.411.211.611.6C101-751154.801068.92−7.441058.85−8.311179.842.210.411.4C201-25671.31671.370.01671.370.01713.606.36.36.3C201-50713.06709.07−0.56709.07−0.56711.46−0.20.30.3C201-75711.45711.460.00711.460.00711.460.00.00.0R101-251651.161644.64−0.391644.64−0.391633.37−1.1−0.7−0.7R101-501644.641644.640.001644.640.001637.06−0.5−0.5−0.5R101-751644.641644.640.001644.640.001642.81−0.1−0.1−0.1R201-251166.641152.25−1.231147.80−1.611148.20−1.6−0.40.0R201-501148.171149.790.141147.80−0.031154.790.60.40.6R201-751147.801147.880.011147.800.001156.310.70.70.7RC101-251717.011654.39−3.651654.39−3.651708.13−0.53.23.2RC101-501770.401724.23−2.611710.32−3.391759.35−0.62.02.9RC101-751784.061779.61−0.251779.45−0.261794.560.60.80.8RC201-251277.051265.56−0.901265.56−0.901284.390.61.51.5RC201-501273.691266.11−0.601265.56−0.641274.420.10.70.7RC201-751266.111266.110.001265.98−0.011272.790.50.50.5Average1267.821245.001243.201273.23To the best of our knowledge, the best results of TTRPTW benchmark instances are reported by Parragh and Cordeau [38]. Parragh and Cordeau [38] proposed a branch-and-price algorithm for the TTRPTW, using problem specific enhancements in the pricing scheme and alternative lower bound computations. They tailored an adaptive large neighborhood search algorithm to the TTRPTW in order to obtain good initial columns. We introduce Gap5 and Gap6 to compare results attained by the TUVRP-TW heuristic with the best results in the literature. For the 18 benchmark instances, the TUVRP-TW heuristic can attain new best solutions for four benchmark instances and the same solution for one benchmark instance. The smallest Gap5 is −0.7%, and the average Gap5 is 2.7%. The smallest Gap6 is −0.7%, and the average Gap6 is 2.9%. Generally speaking, the TUVRP-TW heuristic can obtain average gaps to best known solutions of less than 2.9%. We investigate the detailed results attained by the TUVRP-TW heuristic and find that instance data affect obviously the attained results. Some TTRPTW instances (e.g., C101-25, C101-50, and C101-75) may not be suitable for the TUVRP-TW situation.


## 5.5. Heuristic Applicability for Large-Scale Instances

The conclusions that are obtained by comparing heuristic solutions with exact solutions for small-scale instances and the TUVRP-TW heuristic results for TTRPTW benchmark instances support the satisfactory performance of the proposed heuristic. Several large-scale instances are solved by the heuristic in the computational experiments. Table6 lists the heuristic results for the 27 large-scale instances. For each large-scale instance in the computational experiments, the heuristic is repeated ten times. The objective value, the number of tours or routes, and the heuristic running time of the best solution among the ten repetitions are reported. Column 1 lists the instance. Column 2 lists the objective value of the initial solution. Column 3 lists the numbers of main tours, subtours, and UV routes in the initial solution. Column 4 lists the objective value of the heuristic. Column 5 lists the numbers of main tours, subtours, and UV routes in the heuristic solution. Column 6 lists the computation time for finding the best solution. Column 7 lists the percentage gap between the objective value of the initial solution and the objective value of the heuristic solution, namely, Gap1.Table 6
Heuristic solutions of 27 large-scale instances.
InstanceInitial solutionOptimized solutionTime (s)Gap1 (%)ObjIm/s/UObjOm/s/UC101-253642.46/0/142802.83/1/162533.423.0C101-504020.45/2/182807.30/0/221211.830.2C101-753857.45/3/192846.10/0/201139.526.2C102-253110.06/1/122235.55/1/112309.128.1C102-503438.34/2/162236.01/0/181385.135.0C102-753472.06/2/192303.50/0/201342.333.7C103-252913.36/3/111638.36/1/72463.143.8C103-502807.04/0/151763.93/1/141715.637.2C103-752975.74/0/171851.10/0/191217.337.8C104-252576.55/2/111250.07/2/51557.251.5C104-502705.23/5/141424.94/2/111787.247.3C104-752721.62/1/161447.65/1/151433.046.8C105-253048.58/4/92067.33/1/142965.732.2C105-503438.45/2/122142.13/1/142740.537.7C105-753518.23/3/152027.60/0/191193.942.4C106-253268.63/0/151904.14/1/122070.441.7C106-503329.75/3/101983.33/3/161770.140.4C106-753352.63/1/172078.81/1/171775.838.0C107-252890.74/1/111771.85/3/82674.238.7C107-502895.64/1/111770.71/1/182280.638.8C107-753075.05/3/141833.11/1/181344.240.4C108-252646.86/1/81511.76/1/82664.542.9C108-503227.93/0/141622.83/1/132023.349.7C108-752951.63/2/161655.41/1/171399.543.9C109-252454.66/2/71242.34/2/83179.649.4C109-502726.93/0/131448.32/1/151790.146.9C109-752797.03/2/141596.82/2/151856.242.9The computation time of the heuristic for each of the 27 large-scale instance does not exceed 0.9 h, and the average computation time is 0.5 h. Hence, the proposed heuristic can solve one large-scale instance in an acceptable time for practices.The initial solution that is constructed by the GRASP is feasible for enterprise operations; however, an initial solution that is constructed using the GRASP is likely to fall into a local optimum. The gap (Gap1) between the objective value of the initial solution and that of the heuristic solution is large. For the 27 large-scale instances, Gap1 is more than 40.0% in 14 cases. The largest value of Gap1 is 51.5%, and the average value of Gap1 is 39.5%. The VNS is highly important for finding a locally optimal solution. For the 27 large-scale instances, the average value of Gap1 is larger than the average value of Gap1 for the 21 small-scale instances. The scale of each large-scale instance is larger than that of any small-scale instance, and the initial solution of a large-scale instance may have more neighborhoods than a small-scale instance.Comparing the initial solution with the final solution that is provided by the heuristic, various optimization strategies can be identified. For example, the number of main tours or subtours of the initial solution and the number of main tours or subtours of the final solution differ substantially. The number of main tours of the initial solution, the number subtours of the initial solution, the number of main tours of the final solution, and the number subtours of the final solution are denoted as NMI, NSI, NMO, and NSO, respectively. For 17 of the 27 large-scale instances, the percentage gap between NMI and NMO, namely, Gap3 = 100% × (NMl − NMO)/NMI, exceeds 33.0%. The average value of Gap3 equals 33.9%. For 11 of the 27 large-scale instances, the percentage gap between NSI and NSO, namely, Gap4 = 100% × (NSI − NSO)/NSI, exceeds 31.0%. The average value of Gap4 equals 31.8%. Hence, additional savings may be realized by optimizing the main tours and subtours. For enterprises, focusing on optimizing the main tours and subtours is more desirable.The results of the computational trials on the 27 large-scale instances demonstrate that the instance data and the estimated values of various parameters (e.g., the vehicle velocity and the vehicle capacity) substantially affect the heuristic solutions. We summarize the management insights as follows, which are obtained by respecting the limitations of the instances and the computational experiments. (i) Catering to the “mothership approach” in the expected delivery practice, the TUVRP-TW includes intersatellite synchronization, which is expected to result in more efficient and more flexible delivery operations. Aiming at decreasing the integrated distance, it is not always suitable to adopt the “mothership approach.” Of the 27 large-scale instances, six instances have no subtours in the heuristic solutions. The heuristic results of all 27 large-scale instances include 73 main tours and 29 subtours; hence, the approach of combining the main tour with one or several subtours is not always widely adopted in these instances. Considering the comparison of the exact solutions for small-scale instances that are obtained by directly solving the TUVRP-TW formulation and the TUVRP-TW-B formulation, we find that the “mothership approach” is suitable for pursuing savings in various practical scenarios. Both practitioners and researchers are suggested to identify suitable practical scenarios to employ truck-UV combinations for delivery. (ii) Using UVs for direct delivery from the depot seems necessary for reducing the integrated distance. For each instance, the number of UV routes exceeds the number of main tours or the number of subtours. Hence, the direct delivery strategy is necessary and important in implementing the “mothership approach.”


## 5.6. Impact of Customer Time Windows

In the computational experiments on large-scale instances, sensitivity analysis is conducted to evaluate the impact of time windows on the large-scale instances. By changing the width of customer time windows, the large-scale instances designed in Section5.1, in which the time window of customer i is ei,li and the width of the time window of customer i is denoted as witw (witw=li−ei), are converted to another two types of instances. First, each large-scale instance is denoted by Cbr10J-p. The only difference between instances Cbr10J-p and C10J-p is the width of customer time windows. The time window of customer i in instances Cbr10J-p is ei−0.5⋅witw,li+0.5⋅witw. Second, each large-scale instance is denoted by Cna10J-p. The only difference between instances Cna10J-p and C10J-p is the width of customer time windows. The time window of customer i in instances Cna10J-p is ei+0.25⋅witw,li−0.25⋅witw. The proposed heuristic is used to solve instances Cbr10J-p and Cna10J-p. For each instance, the proposed heuristic is run 10 times, and the objective value and the computation time of the best solution of the 10 results are reported in Table 7.Table 7
Heuristic results of large-scale instances with changed time windows.
InstanceObjOTime (s)Gap3InstanceObjOTime (s)Gap4Cbr101−252101.42814.0−25.0Cna101−253300.92357.017.8Cbr101−502144.51358.9−23.6Cna101−503380.71233.820.4Cbr101−752205.91215.8−22.5Cna101−753520.31295.423.7Cbr102−251750.82663.1−21.7Cna102−252663.52737.619.1Cbr102−501742.91808.0−22.1Cna102−502949.72406.731.9Cbr102−751796.21660.3−22.0Cna102−752743.81301.419.1Cbr103−251445.92940.9−11.7Cna103−252002.02138.322.2Cbr103−501575.32221.5−10.7Cna103−502069.52260.417.3Cbr103−751657.51246.5−10.5Cna103−752200.01264.618.9Cbr104−251193.92336.4−4.5Cna104−251398.72187.911.9Cbr104−501356.32380.4−4.8Cna104−501556.73175.79.2Cbr104−751403.71370.1−3.0Cna104−751707.01414.117.9Cbr105−251520.52763.8−26.5Cna105−252774.72420.734.2Cbr105−501570.92040.4−26.7Cna105−502901.92318.635.5Cbr105−751685.12322.3−16.9Cna105−752761.61591.436.2Cbr106−251428.82444.8−25.0Cna106−252462.31722.329.3Cbr106−501604.62260.6−19.1Cna106−502050.31222.63.4Cbr106−751610.91430.5−22.5Cna106−752587.11342.124.4Cbr107−251300.13063.9−26.6Cna107−252158.52688.621.8Cbr107−501474.52073.7−16.7Cna107−502291.41684.729.4Cbr107−751561.71294.7−14.8Cna107−752230.11277.721.7Cbr108−251160.33752.1−23.2Cna108−251923.22614.627.2Cbr108−501398.33265.9−13.8Cna108−502201.82195.135.7Cbr108−751426.41468.5−13.8Cna108−752033.01233.022.8Cbr109−251018.43337.0−18.0Cna109−251541.72825.224.1Cbr109−501191.32453.6−17.7Cna109−501645.31811.613.6Cbr109−751353.61807.3−15.2Cna109−751723.51297.07.9To compare the objective values of solutions of instances Cbr10J-p and C10J-p, we introduce a percentage gap, Gap3, where Gap3 = 100% × (the objective value of solution of instance “Cbr10J-p” − the objective value of solution of instance “C10J-p”)/the objective value of solution of instance “C10J-p.” To compare the objective values of solutions of instances Cna10J-p and C10J-p, we introduce a percentage gap, Gap4, where Gap4 = 100% × (the objective value of solution of instance “Cna10J-p” − the objective value of solution of instance “C10J-p”)/the objective value of solution of instance “C10J-p.” The average, largest, and the smallest Gap3 values are −17.7%, −3.0%, and −26.7%, respectively. The average, largest, and the smallest Gap4 values are 22.1%, 36.2%, and 3.4%, respectively. It indicates that wide time windows may lead to the decrease of the integrated distance while narrow time windows may lead to the increase of the integrated distance. In terms of the heuristic computation time, the proposed heuristic is able to solve each of the 27 Cbr10J-p instances in an average computational time of 2215 s and is able to solve each of the 27 Cna10J-p instances in an average computational time of 1927 s. Wide time windows may enlarge the neighborhood, so that more computation time is spent by the heuristic to find satisfactory solutions. Narrow time windows may cut the neighborhood, so that less computation time is spent by the heuristic to find satisfactory solutions.


## 6. Conclusions

The satellite synchronization that is considered in the 2E-VRP, the TTRP, and the 2E-LRP models in the literature, which is called inner-satellite synchronization in this paper, focuses on synchronization operations at a single satellite. We consider an important realistic feature in a practical scenario that logistics and supply-chain enterprises are presently facing in using UVs that are carried by trucks for “last-meter” delivery, and we propose the TUVRP-TW as a suitable methodology for overcoming the operational challenges in optimizing delivery routes for truck-UV combinations. We introduce intersatellite synchronization, which focuses on synchronization operations at multiple satellites. Using intersatellite synchronization, UVs need not return to the dispatched locations.The TUVRP-TW involves a homogeneous fleet of truck-UV combinations for deliveries. UVs may be used for direct delivery from the depot. A UV that has been dispatched by its truck must be picked up by the same truck or must return by itself to the depot. Customers are classified into two types: TUCs and UCs. The TUCs where trucks drop off or pick up carried UVs are regarded as satellites. The TUVRP-TW involves truck-UV combinations or trucks delivering cargoes from the depot to TUCs, involves carried UVs being dispatched at satellites to visit TUCs or UCs, and involves UVs delivering cargoes directly from the depot. The TUVRP-TW objective is to minimize the integrated distance. To formulate the TUVRP-TW, both binary variables for identifying the appointed satellites and continuous variables for time continuity constraints are introduced to ensure the interaction between truck routes and UV routes.We propose a mixed-integer linear programming model that can be solved directly using CPLEX. Compared with all types of models of the 2E-VRP, TTRP, and 2E-LRP, intersatellite synchronization constraints provide an innovative approach for formulating routing problems with satellite synchronization. A hybrid algorithm based on a GRASP and a VNS is proposed. We use computational experiments to evaluate the performances of the intersatellite synchronization and of the TUVRP-TW formulation and the applicability of the heuristic. Through a comparison of the exact solutions for small-scale instances, the satisfactory performance of the intersatellite synchronization is demonstrated. The results on small-scale and large-scale instances demonstrate that the TUVRP-TW formulation and the heuristic are effective. Besides, for 18 TTRPTW benchmark instances, the TUVRP-TW heuristic can attain new best solutions for four benchmark instances and the same solution for one benchmark instance.For future research, extensions of the TUVRP-TW by considering, e.g., one truck carrying several UVs, are expected. In addition, more effective heuristics for the TUVRP-TW should be developed.



---

*Source: 6599089-2022-05-12.xml*