Design and Assessment of an Urban Circular Combined Truck–Drone Delivery System Using Continuum Approximation Models and Integer Programming

Citation: Canca, D.; Navarro-Carmona, B.; Andrade-Pineda, J.L. Design and Assessment of an Urban Circular Combined Truck–Drone Delivery System Using Continuum Approximation Models and Integer Programming. Sustainability 2022,14, 13459. https://doi.org/10.3390/ su142013459 Academic Editor: Tamás Bányai Received: 23 September 2022 Accepted: 17 October 2022 Published: 18 October 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). sustainability Article Design and Assessment of an Urban Circular Combined Truck–Drone Delivery System Using Continuum Approximation Models and Integer Programming David Canca 1,* , Belén Navarro-Carmona 1and José Luis Andrade-Pineda 2 1Department of Industrial Engineering and Management Science, School of Engineering, Universidad de Sevilla, Av. de los Descubrimentos s/n, 41092 Seville, Spain 2Robotics, Vision and Control Group, School of Engineering, Universidad de Sevilla, Av. de los Descubrimentos s/n, 41092 Seville, Spain *Correspondence: [email protected] Abstract: The analysis of tandem truck–drone delivery systems has recently attracted the attention of the research community, mainly focused on extending classical operational research problems such as the multiple traveling salesperson or the vehicle-routing problem. In this paper, we explore the design of an urban massive combined delivery system using a continuum approximation (CA) method for a circular city characterized by a certain density of customers. Starting from a set of parameters defining the main characteristics of trucks and drones, a sectorization of the delivery area is first determined. Then, for a given truck capacity, the optimal number of trucks is obtained considering different scenarios using three integer programming models. We propose several performance indicators to compare the tandem approach with the alternative solely truck delivery system. Keywords: transportation; last-mile delivery; truck–drone delivery; continuum approximation; mixed-integer programming 1. Introduction The new express-delivery services demanded around e-commerce sourced by citizens in urban areas are transforming the operational model of courier firms. Where, in the past, they might focus on visiting a set of streets with minimal travel time constrained to energy-feasible routes, presently the goal is connected to the definition of delivery plans that alleviate the disruptions caused by the current fleet of electric-powered delivery vehicles (avoiding congestion and transport externalities) vans (others have forbidden access to the city center). This paper explores the inclusion of drones for the delivery of packages, thus minimizing the need for vans in the city center. Although the hybrid drone–truck problem has been largely studied in recent years, the focus has generally been put on operational-level decisions—specifically, finding the optimal path for the truck and the optimal schedule for the drones—for serving a discrete number of locations, while balancing the delivery cost and time. However, in our research, the density of the concerned e-commerce customers is so high that we prefer to aim the hybrid truck–drone delivery system at the definition of competitive long-term routes for several vans in the city center (namely the members of the fleet following a fixed set of routes, no matter where the specific daily customers were). Take, as an example, the case of New York City, where the average daily density of customers reaches 1834 customers/km 2 —see the article published in the New York Times “1.5 Million Packages a Day: The Internet Brings Chaos to N.Y. Streets”, 28 October 2019, which refers to population (8,336,817 inhabitants) and average population (10,194 inhabitants/km 2 ) data [ 1 ]. Hence, our decision-making process on the last-mile delivery in large cities considers that customers are distributed throughout the city, regardless of the specific locations. Sustainability 2022,14, 13459. https://doi.org/10.3390/su142013459 https://www.mdpi.com/journal/sustainability Sustainability 2022,14, 13459 2 of 30 In addressing these strategic/tactical decisions (size of the fleet, delivery area assigned to each van), we use continuum approximation (CA) methods, replacing the numerical solution of problems by analytical techniques which allow us to analyze general problems with imprecise data [ 2 , 3 ]. Under this approach, the problems under study can be formulated using a relatively small number of design parameters, and the effect of these parameters on the obtained results can be analyzed. Importantly, the CA approach allows us to maintain the tractability of the problem, focusing our attention in the key issues and to obtain reasonable solutions with as little information as possible [2,4]. According to [ 5 ], the use of CA models that incorporate the most relevant cost components and capture the most important cost trade-offs can create accurate and useful representations of logistics systems. As in many other truck–drone delivery systems, we assume that drones provide service to customers and then make a return trip to a truck or van that is itself moving in its assigned area. Moreover, the truck acts as a mobile depot from which the drones pick up a package and where it receives replenished batteries when necessary. In short, we have selected a continuum approximation approach wherein the inclusion of in-route synchronization between truck and drones definitely influenced the hybrid delivery service design. Specifically, we analyze the worst-case scenario concerning drone autonomy, although this consideration overestimates the real delivery capacity of the combined system; otherwise, the analytic expressions that allow the determination of the system design would become intractable and scenario-dependent, refusing effectiveness to the general evaluation we intend. Regarding the need for synchronizing trucks and drone movements, we first outline the parameters which determine the partition of the city into delivery areas. We then propose a mixed-integer programming (MIP) model to assign trucks to delivery cycles while minimizing the total delivery time for a given number of trucks. Finally, two MIP formulations are proposed to optimize the truck-only system, devising two performance measures to compare both systems. Summarizing, in this research, we propose a methodology using both CA methods and mixed-integer programming to determine the system design, i.e., determine the number of sectors in which the study area must be divided and assign ground vehicles to delivery routes. The use of a continuum approach allows us to determine the number and dimension of the delivery zones (depending on the endurance and speed of drones and the capacity of trucks), whereas the MIP models are used to help assess the system performance by assigning ground vehicles to the delivery zones. Importantly, our purpose is not the exact determination of all the details involved in a massive delivery system (which definitively will depend on the specific city structure) , but obtaining a set of quantitative measures to demonstrate the benefits of using a combined truck–drone system when the specific locations of customers are not fixed. The remainder of the paper is structured as follows. Section 2reviews the related work. Section 3presents the description of our problem. Section 4establishes the set of assumptions that have been considered. Section 5remarks the main contributions of this research. Section 5presents the equations to determine the truck–drone system design and proposes a MIP formulation to obtain, as a function of the model parameters, the minimum delivery time for the combined system. Section 6contains two MIP formulations for the truck-only delivery system. The first minimizes the delivery time for a given number of trucks. The second minimizes the number of vehicles for a given delivery time. Both formulations are necessary to make a fair comparison against the combined delivery system. Section 7contains the definition and the solutions, respectively, of a set of experiments designed to assess the effect due to the change in the different design parameters and to report the comparison between the combined and the truck-only delivery systems. Finally, Section 8presents several conclusions. Sustainability 2022,14, 13459 3 of 30 2. Related Work The use of drones for logistics operations is increasingly being adopted. The flexibility to follow free routes makes them very attractive for last-mile delivery of goods, surveillance activities, medical samples transportation and drug delivery in inaccessible places, among others. Regarding the last-mile delivery problem studied in this paper, there are also shortcomings to be considered: the flight restrictions in urban areas usually imposed by authorities, the limited payload they can carry, and the short traveling radius (flight time bounded by lifetime of batteries). However, authorities are opening their regulation, and the advent of commercially massive drone-based services is coming, and the technology is ready to cope with the other two cons. The battery endurance limitation can be partially solved by combining the drones with ground vehicles (generally called trucks) that act as mobile depots to pick up packages up and where drones can take advantage of battery replacements. Last-mile delivery is by far the most important research stream on hybrid truck–drone logistic problems, since the seminal work by [ 6 ]. Some authors give an active role to ground vehicles, assuming that they can also deliver the products to the customers (see, e.g., [ 6 – 9 ]) while, in other works, it is assumed that the truck only performs a secondary role, interacting solely with drones (see, e.g., [ 10 , 11 ]). If the truck launches the drones and waits at the same position until their return, no synchronization problems arise, but, in general, if the truck launches drones and moves to a different location, the routing plan must deal with synchronization issues, since, in general, one of the vehicles must wait to the arrival of the other at rendezvous points. The first approximation to a combined delivery system using trucks and drones is due to [ 6 ]. These authors proposed two types of combined truck–drone delivery problem: the flying sidekick traveling salesman problem (FSTSP), in which the truck must visit some fixed customers while the rest of the parcels are delivered by a drone carried by the truck, and the parallel drone scheduling traveling salesman problem (PDSTSP) , in which a significant proportion of customers are located within the drone flight range from a depot. The drones depart from the depot, deliver parcels to customers, and return to the depot while, at the same time, in a parallel fashion, the truck visits several customers, delivering products without interacting with the drones. For a relatively similar problem, [ 12 ] consider that a single truck carries the drones to some locations from which drones fly to serve customers, returning to the truck after the delivery has been completed. The combined truck–drone delivery situations analyzed have evolved from the initial works to a set of more complicated problems, mostly focused on extending classical routing problems—a variant of the traveling salesman problem with drones (TSP-D) [6,7,13–16] , and a generalization of vehicle-routing problems to include drones (VRP-D) [ 17 – 20 ]. Noticeably, most of the existing literature approaches how to serve a set of discrete locations through a mathematical optimization formulation for the considered problem and a metaheuristic-type solving procedure. For example, ref. [ 21 ] proposed a simulated annealing (SA) heuristic to solve the FSTSP. Ref. [ 14 ] used a dynamic programming (DP) approach for the FSTSP. Ref. [ 22 ] proposed a greedy randomized adaptive search procedure to solve an FSTSP. Ref. [ 13 ] proposed a route-first cluster-second heuristic method for the same problem. Recently, ref. [ 15 ] extended the FSTSP problem to include the existence of nonfly zones and the effect of payloads on the drone energy consumption and proposed a two-phase constructive and search heuristic algorithm to deal with real-world problems. For a comprehensive recent review of the current research on drone-based logistics, the reader is referred to the works by [23,24]. According to the classification of drone-base logistics system by [ 24 ], our motivating problem settles into the category “multimodal synchronized”. As in the seminal work by [ 6 ], we consider the commonly accepted hypothesis in the last-mile delivery literature— except for only a few works, see, e.g., [ 25 – 28 ]—for which only one single location is visited on each drone trip. However, concerning the movement until the synchronization meetings, we consider that once the truck launches a drone it does not stay at the launch Sustainability 2022,14, 13459 4 of 30 location to retrieve the drone. Although the common assumption in this research stream has been that synchronizations only occur while the truck stops at one of the discrete locations, ref. [ 19 ] have recently considered the possibility of launching and retrieving drones along a route—i.e., at discrete locations other than the customer locations—and have reported experimental results for their VNS/Tabu search heuristic applied on up to 50 customer instances: (they claimed the potential of a reduced completion time of the delivery service and a higher use of drones). Furthermore, some other authors [ 3 , 29 , 30 ] adopted a continuous approach to allow drones to meet the truck at any point of the route between two locations. Ref. [ 30 ] compare two delivery models after dividing the service region into hexagonal areas. The first model allows only drones to serve customers, while the second model allows both trucks and drones to serve customers. In their CA approach, mathematical cost formulations are obtained as a function of the area size, density of customer locations, and the number of customers that can be served by drones. The authors show that once the density of customer locations increases, it is more cost effective to allow both trucks and drones to serve customers. In [ 29 ], drones are launched from the truck to simultaneously serve customers while the truck moves to the next customer location. A cost function is provided which can be used to assess the optimal number of trucks and drones and the number of drones in each truck. The results demonstrate that this model can produce considerable cost-savings compared to a truck-only model, especially in rural areas and when multiple drones per truck are used. In addition to the operating cost of the drone and truck, the relative stop cost for the drone and truck, and the spatial density of customers are other major factors that affect the cost performance of this model. In contrast to them, our study is concerned only with supporting the role of trucks and on addressing the last-mile service in highly populated urban areas. Finally, the work by [ 3 ] is more directly connected to our motivating case study; they assume that a truck travels in continuous planar space and dispatches the drones to deliver packages to the customers. These customers can only be served by drones, while the truck serves as a moving hub for the drones. A continuous approximation method is used to approximate the length of the truck route while a heuristic is provided to assign and schedule the drone trips from the truck to the customers. The main conclusion of their analysis indicates that the total savings in delivery time of this model is proportional to the square root of the truck-to-drone speed ratio. CA models have been applied in many fields. For example, ref. [ 31 ] made decisions about the production-distribution system design problem by developing CA models representing spatial distributions of cost and customer demand. Ref. [ 32 ] demonstrated that CA models can represent the exact model with an accepted gap and that cost differences are insensitive to gradual demand variations. Ref. [ 33 ] estimated the distance of a vehicle-routing problem using CA methods to solve the facility location problem. Ref. [ 34 ] developed CA models for the location of stores in a retailer network. Concerning the application of CA models to transportation and logistics problems, the review of [ 35 ] summarized the first scientific contributions in the field of freight distribution problems, stressing important principles and key results from continuum approximation models. The more recent review by [ 36 ]) included many studies that develop CA models for transportation, distribution, and logistics problems, identifying current research gaps. Ref. [ 37 ] developed an analytical model for designing a hybrid grid network that combines flexible demand-responsive services with a fixed route service. The objective of the model was to determine the optimal number of zones, while each zone is served by several on-demand vehicles. In the field of location analysis, ref. [ 38 ] presented a game-theoretical model based on a continuum approximation (CA) scheme to determine the location of the service facility under spatial competition and facility disruption risks. They first analyzed the existence of Nash equilibria in a symmetric two-company competition case and then developed a leader-follower Stackelberg competition model to derive the optimal facility location design when one of the companies has the first-move advantage over its competitor. The authors solved both models and devised closed-form analytical solutions for special cases. Sustainability 2022,14, 13459 5 of 30 Concerning the design of public transportation systems, ref. [ 39 ] developed a methodology for designing bus networks for cities where travel demand varies gradually over space. Compared to homogeneous bus-route grids, the resulting heterogeneous route configurations can reduce user and operating agency costs. Much of the savings are due to the lower access costs that users experiment when high-demand neighborhoods are served by local grids with closely spaced routes and stops. Ref. [ 40 ] formulated two CA models for designing city-wide transit systems at a minimum cost. The paper shows how to solve these CA optimization problems numerically. As main results, the optimal headways and spacings in the periphery increase with the distance from the center and, at the boundary between the central district and the periphery, both the optimal service frequency and the line spacing for radial lines decrease abruptly in the outbound direction. The results suggested that the proposed CA procedures can be used to design transit systems over real street networks when they are not too different from the ideal and that the resulting costs might usually be very close to those predicted. Ref. [ 41 ] focused on the planning and design of ring-radial rail transit systems. An analytical model to find the optimal number of radial lines in a city is first introduced. The passenger route choice for different rail networks is analyzed considering a per passenger total travel time cost objective. The authors considered different cases, radial lines only, ring lines only, or a combination of radial and ring lines. By assessing changes in passenger costs, this study shows the potential benefit of introducing a ring line in the design. In analyzing the problem of a demand adaptive paired-line hybrid transit system, ref. [ 12 ] studied the differences in operating the demand adaptive service along circular or radial transit lines. A continuum approximation approach is employed to develop the optimal design problem, which is formulated as a mixed-integer model. A set of numerical experiments is performed to compare various cost components corresponding to the optimal design of the two systems, and a discrete-event simulation model is developed to validate the analysis. Ref. [ 42 ] proposed a CA modeling framework to optimize a hybrid transit grid network in the central district of a city and a hub-and-spoke structure in the periphery. Two CA models are formulated incorporating a local route service and a short-turn strategy, respectively. The configuration of the transit network is optimized considering both the cost of passengers and the cost of the agency. Ref. [ 43 ] incorporated spatial heterogeneity into the optimal design of paired-line hybrid transit systems, with the aim of striking a better balance between accessibility and efficiency. The authors proposed a simple trip production and distribution model to differentiate the central business district (CBD) of a city from its periphery. Closely related to our research, returning to the combined truck–drone delivery problem and assuming that the delivery area is characterized by a continuous distribution of customers, ref. [ 3 ] used a CA method to approximate the length of the truck route while a heuristic is provided to assign and schedule drone trips from the truck to the customers. Ref. [ 29 ] proposed a CA model to obtain important design parameters such as the optimal number of drones per truck, and compare the delivery cost against a truck-only method for different scenarios. The authors develop a set of formulas to calculate the expected delivery cost and time. Using this methodology. Ref. [ 30 ] explored the economics of combined delivery systems using trucks and drones. Distances and costs are approximated as simple functions using CA methods. The main research contributions of this paper are the development of cost models and the analysis of delivery activities focusing on the trade-off between the major components of the problem. 3. Problem Description Consider a circular area (let us say a city) of radius R measured in km, which is characterized by a density of customers ρ measured in customers/km 2 . A depot is located at the center of the area from which a certain number of trucks T depart following radial trajectories. The fleet of trucks is homogeneous, and each truck has a certain capacity P of products and carries a certain number of drones nd and a given number of replacement Sustainability 2022,14, 13459 6 of 30 fully charge batteries nb . As we will see later, for efficiency reasons, the ratio ˆ n between the number of replacement batteries nb and the number of drones nd must be an integer greater than or equal to 1. Each truck will move along a radial street until a certain position, from which it will start a circular concentric movement (no matter about the direction, that is, the whole system can follow a clockwise or a counterclockwise direction or even both directions in different parts of the city, according to the design of the road network). At this stage, different customer orders–inside a certain area A around the truck trajectory–are served by drones, one by one, each drone carrying one order, until the truck is empty–all the carried products have been delivered–or the drones cannot do more trips–the number of fully charged replacement batteries is zero. Then each truck returns to the depot following a radial street to be loaded again and to replenish new fully charged batteries. The described system divides the city into a set of concentric circular crowns so that the trucks’ trajectories are located at the central circumference of each circular crown. The delivery area of each drone cycle A is then a sector of a circular crown. At this moment, the number of circular crowns N , the stopping position of each truck, the width of the delivery area around the truck trajectory W and the length of the truck trip over the crown S are unknown and dependent on the truck and drone characteristics as well as on the density of customers. Figure 1a depicts a complete truck cycle whereas Figure 1b shows to trucks operating at different circular crowns. Please note that the number of circular crowns { 1, . . . , N} are numbered from outside the city to the central depot (the inner circle is numbered as N and is a priori unknown). Regarding drone movements, Figure 2a depicts the way a drone is launched from the truck for delivery purposes. The drone follows a trajectory composed of two movements; the first part of the trip is a radial flight from the truck to the customer location. After delivering the parcel, the drone returns to the truck following a straight line. Please note that during the drone flight the truck continues its circular trips over the central circumference of the corresponding circular crown. In the worst case, see Figure 2b, considering the i th circular crown, the customer will be located at the border of the circular crown. Supposing that the drone is launched at the first delivery instant, the total distance traversed by the drone will be W/ 2 +S0 i and the maximum length of the truck trip S (which depends on the truck speed) must be enough to receive the drone. At this point, if the truck has a replacement charged battery, the battery of the incoming drone can be replaced, and a new customer order can be dispatched. When necessary, the truck can stop at the end of its complete trajectory, before returning to the depot, and wait for the return of the latest launched drones. (a) (b) Figure 1. Truck trajectories. (a) Truck cycle. (b) Area covered by two trucks at different cycles. Table 1summarizes the notation used for the different parameters involved in the determination of the geometry of the combined truck–drone delivery system. A similar approach could be followed for a radial delivery system (where trucks move in radial Sustainability 2022,14, 13459 7 of 30 directions and drones fly perpendicularly to trucks); however, the radial structure gives rise to a variable trapezoidal delivery area which changes with the distance to the center of the city, thus complicating the development of the analytical expressions needed to define the main system structure as well as the practical operation, since the delivery areas are not similar in the radial direction. (a) (b) Figure 2. Drone trajectories. (a) Drone launching scheme. (b) Detail of the worst-case drone flight. Table 1. Parameters and variables for determining the structure of the hybrid delivery system. Parameters RRadius of the delivery area (km) ρCustomers density (customers/km2) vmin dLower bound for the drone speed (km/h) vmax dUpper bound for the drone speed (km/h) vmin cLower bound for the truck speed (km/h) vmax cUpper bound for the truck speed (km/h) EDrone battery endurance (min) PNumber of products carried by the truck ndNumber of drones carried by the truck nbNumber of batteries carried by the truck ˆ nRatio that represents the relationship between the number of batteries and the number of drones Variables WWidth of each circular crown (km) SLength of the circular truck trip from the start of the circular movement to the reception of a drone launched at the start of the trip (km) ADelivery area for a truck with as many batteries as customer orders (km2) NNumber of circular crowns S0 iLength of the drone movement from a customer located at the border of the circular crown ith to the truck (km) θiAngle corresponding to the distance Sat circular crown ith (radians) vdDrone speed (km/h) vcTruck speed (km/h) Please note that our purpose is not the exact determination of all the details involved in a massive delivery system, which definitively will depend on the specific city structure; in contrast, as in other CA analysis, our goal is to obtain a set of general measures which allow us to demonstrate the benefits of using a combined truck–drone system when the Sustainability 2022,14, 13459 8 of 30 specific locations of customers are not fixed and to obtain quantitative measures of such advantage, even when several simplifications are made, in a similar way of [12,41,42]. 4. Assumptions • We are interested in designing a set of daily delivery routes to serve a large circular area divided into a set of zones characterized by a certain density of customers. The specific locations of customers inside a zone will vary every day, but the trucks follow a fixed pattern in such a way that they cover the entire city. • For simplicity, we consider that the drone battery endurance is constant for the range of drone speeds considered in the experiments, between 30 and 45 km/h. (The inclusion of a law relating the endurance to the drone speeds complicates the resolution of the models and does not imply significant changes to the results. Furthermore, for our purposes, battery endurance could be pessimistically estimated considering the upper speed value of the speed range, providing us with a worst-case value that could be used as data in our experiments). • We assume that all the batteries are in perfect state, that they can be charged to the same level, and that their endurance is the same for all the replacement batteries. • We suppose that the road network capillarity (the density of roads in the network) is large enough to accommodate the horizontal and radial movements of trucks, i.e., there is always a very close street that allows for the movement of trucks in both circular and radial directions as obtained after solving our model. Please note that similar assumptions have been previously considered in previous research, see, for instance, [12] or [40], and it is reasonable when dealing with CA models that involve routing vehicles in extensive areas. • Without loss of generality, we consider a constant time of truck loading of 30 min for all the experiments, i.e., the time needed to charge new products–to be delivered in the next delivery cycle–and fully charged batteries (in the case of the combined truck–drone delivery system). This parameter value is used in both combined and truck-only designs, not affecting the general comparison of their performance. • In both delivery systems, truck–drone and truck-only, we assume that the delivery time, once reached the customer location, is negligible. • To compute the truck emissions of CO 2 we suppose that the per liter fuel emission is constant and independent of the truck speed for the given range of truck speeds, between 30 and 60 km/h. Again, if convenient, the worst value could be used in the case of the combined truck–drone delivery system if a worst-case analysis is considered. 5. Contributions The main contributions of this research are as follows: • We analyze the performance of a combined truck–drone delivery system for extensive areas. In contrast to the usual MIP methodology to determine routes for a given set of customer locations, we consider here that the customers are spread throughout the city, regardless of specific locations, and consequently we have selected a continuum approximation approach to deal with this problem. • We take care of the synchronization issues between trucks and drones, which definitively influence the design of the system. To this end, for synchronization purposes, we will analyze the worst-case scenario concerning the drone autonomy, even if this consideration overestimates the real delivery capacity of the combined system; otherwise, the analytic expressions that allow the determination of the system design become intractable and scenario-dependent, refusing effectiveness to the general evaluation we intend. • We alternatively analyze the entire delivery process using only trucks. To this end, we consider that the only-truck system follows the same circular design as in the Sustainability 2022,14, 13459 9 of 30 combined system, but redefining the parameters to the best values for delivering using only trucks. • We propose a methodology using both continuum approximation methods and mixedinteger programming to determine the system design. The continuum approach allows us to determine the number and dimension of the delivery zones. The MIP models are used to help assess the system performance by assigning the ground vehicle to delivery zones. • We perform a set of numerical experiments to obtain key performance indicators of both delivery systems for different values of the main parameters. 6. The Combined Truck–Drone Delivery System 6.1. Determining the Width and the Number of Circular Crowns To determine the main parameters of the combined truck–drone delivery system; the width and the number of circular crowns needed to deliver items to customers, we shall first analyze the maximum length S of a truck movement over the central line of a circular crown. As expected, the dimensions of the delivery zones will depend on the density of customers and the endurance of the drone battery; moreover, the need for synchronization between the drones and the trucks will play a very important role in determining the number of delivery zones. Remember that the circular crowns are numbered from outside to the depot (located in the city center). – Selecting the worst-case scenario Since, according to the Figure 2b, which represents the worst-case scenario, a drone can be launched from the truck following a radial movement until reaching the outside border of the corresponding circular crown and then it returns to the truck following a straight trajectory, the length of the truck movement S has to be calculated to ensure the drone endurance allows it to reach the ground vehicle. Then, the time needed for the drone to cover a distance of W/ 2 +S0 i with speed vd must be less or equal to the time spent by the truck to travel the distance S with speed vc . Moreover, since the distance S0 i differs from a circular crown to another, see Figure 3, to ensure the truck–drone synchronization, the calculations must be performed in the worst case, which, as depicted in what follows, corresponds to the circular crown number 1 –the most outside one–. Figure 3. Drone trajectories at two consecutive circular crowns. Sustainability 2022,14, 13459 16 of 30 Since each truck is carrying P products, the delivery area covered by a truck at circular crown ith must be equal to P/ρ. P ρ=θi 2h(R−(i−1)W)2−(R−iW)2i= =θi 2h2RW +W2(1−2i)i(20) Then, the angle θicorresponding to the delivery area can be obtained as: θi=2P ρ[2RW +W2(1−2i)](21) The length of the trip performed by each truck at circular crown ith will be: li=θiR−(2i−1)W 2(22) Thus, the number of cycles corresponding to the circular crown i th is obtained by rounding to the upper nearest integer the quotient between the length of the central circumference and the length of the trip li. Ci=    2πR−(2i−1)W 2 li    (23) Once computed the number of delivery cycles for each circular crown, the length of the delivery trip per cycle is adjusted, obtaining a set of cycles with the same delivery trip length. The angle of the delivery area θiis also updated. ¯ li=θiR−(2i−1)W 2 Ci (24) ¯ θi=¯ liR−(2i−1)W 2(25) The number of products carried by each truck is now: P=ρ¯ θi 2h2RW +W2(1−2i)i(26) And the total length of each delivery cycle at circular crown ith will be: lci=2PW 4+2R−(2i−1)W 2+¯ li(27) Where the first term takes into account the average distance corresponding to the radial movements from the central line of circular crown i th to the customers and the returning movement to the central line. The second term corresponds to the length of movements from the depot to the central line of the circular crown, and the final movement to the depot. The last term includes the circular truck movement inside the delivery area. 7.1. Assigning Trucks to Crowns and Equilibrating the Distance Traveled by Trucks–Model B1 After explaining the truck-only delivery system and determining the angle, the total trip delivery length, and the cycle length (circular movement length), in what follows, we propose a MIP model with the objective of equilibrating the distance traveled by each truck, which is equivalent to equilibrate the delivery time per truck. Please note that both Sustainability 2022,14, 13459 17 of 30 the number of circular crowns N and the width of each circular crown W were obtained after solving the combined truck–drone delivery model in Section 6.2. Assuming a fixed fleet size of ˆ T trucks, the assignment of trucks to cycles can be obtained by solving the next MIP model, where xij is an integer variable that represents the number of cycles of circular crown i that are assigned to truck j={ 1, . . . , T} . Please note that in this model, to perform a reasonable comparison, the total delivery length (or total delivery time) will be obtained for a given number of vehicles, which corresponds to the fleet size calculated for the combined truck–drone delivery system. Min L (28) subject to : Lj≤L j =1, . . . , ˆ T(29) Lj= N ∑ i=1 lcixij j=1, . . . , ˆ T(30) Ci= ˆ T ∑ j=1 xij i=1, . . . , N(31) xij ∈N+(32) The objective function (28) minimizes the maximum length traveled by any truck, Lj , j== 1, . . . , ˆ T . Constraints (29) bound the length assigned to truck j with the variable L , which is minimized in the objective function, thus equilibrating the distance traveled by the trucks. Constraints (30) compute the total length traveled by truck j as a function of variables xij . Finally, constraints (31) ensure that all the cycles of each circular crown are assigned to trucks. As a result, the distance traveled by trucks is equilibrated. Using the optimal values of variables x∗ ij of Model B1 ( (28) – (32) ), the total delivery time spent by truck jis: Tj= N ∑ i=1 lci vc x∗ ij (33) 7.2. Assigning Trucks to Crowns and Minimizing the Number of Required Trucks–Model B2 Instead of fixing the number of trucks, an interesting alternative consists of fixing the delivery time, using the value obtained in the combined truck–drone delivery system, and thus determining the required number of trucks. From this point of view, we are interested in obtaining the size of the delivery fleet needed to reach a delivery time as good as the one determined in the combined delivery system. However, depending on the values of the parameters affecting the solution, it is possible that even with a high number of vehicles the system delivery time cannot be less than or equal to the optimal value obtained in the combined system. To practically avoid infeasibility, a positive excess variable is included to relax the delivery time per truck. Again, integer variables xij represent the number of delivery cycles for the circular crown i th assigned to vehicle j th. Since the number of vehicles must be minimized, for each vehicle, a binary variable δj , taking value 1 if truck j th is used for the delivery of goods, is included. Denoting by ˆ C the total number of cycles in the full delivery area, which can be easily obtained as the sum of the number of cycles of all the circular crowns, ˆ C=∑N i=1Ci, the MIP model is formulated as follows: Sustainability 2022,14, 13459 18 of 30 Min ˆ C ∑ j=1 δj+α ˆ C ∑ j=1 Hj(34) subject to : Tj≤t∗ j+Hjj=1, . . . , ˆ C(35) Lj= N ∑ i=1 lcixij j=1, . . . , ˆ C(36) Tj=Lj vc j=1, . . . , ˆ C(37) Ci= ˆ C ∑ j=1 xij i=1, . . . , N(38) N ∑ i=1 xij ≤ˆ Cδjj=1, . . . , ˆ C(39) xij ∈N+,δj∈ {0, 1},Hj≥0 (40) The objective function (34) minimizes the number of trucks and the sum of excess variables, which is weighted by a parameter α to measure the importance of this term in the minimization. Constraints (35) bound the delivery time of truck j to the value obtained after solving the combined truck–drone system Model A ( (9) – (13) ). If possible, Hj variables should take values equal to zero; otherwise, it is not possible to complete the delivery of parcels spending a system delivery time as in the combined system. Constraints (36) compute the total length traveled by truck j th as a function of variables xij . Constraints (37) express the total delivery time of truck j th. Constraints (38) ensure that all the cycles of each circular crown are assigned to trucks. Constraints (39) are used to relate variables xij and δj . If truck j th does not take part in the delivery of goods, none of the area cycles can be assigned to it. As a result of the model, the distance traveled by trucks and the truck delivery times are equilibrated. 7.3. Key Performance Indicators for the Truck-Only Delivery System In a similar way to the combined truck–drone delivery system, two key performance indicators are formulated. (a) System delivery time Using the optimal values of variables x∗ ij in the Model B1 ( (28) – (32) ) or in Model B2 ((34)–(40)), as appropriate, the total system delivery time DT will be: DT =max j=1,...,NT{Tj}, (41) where NT =ˆ Tin the case of Model B1 and NT =∑ˆ C j=1δ∗ jwhen using Model B2. (b) Co2emissions Using the optimal values of variables x∗ ij in the Model B1 ( (28) – (32) ) or Model B2 ((34)–(40)), as appropriate, the total system CO2emissions are: Co2=lg κ ˆ NT ∑ j=1 L∗ j, (42) where NT =ˆ Tin the case of Model B1 and NT =∑ˆ C j=1δ∗ jwhen using Model B2. Sustainability 2022,14, 13459 19 of 30 8. Experimental Analysis 8.1. Definition of Experiments for the Combined Truck–Drone Delivery System Table 2summarizes the definition of the experiments designed to assess the performance of the combined truck–drone delivery system. The first column contains the parameters of Model A: customer density ( ρ ), number of products ( P ), drone endurance ( E ), ratio between the number of replacement batteries and drones ( ˆ n ), upper and lower bounds on drone speed and truck speed ( vmax d , vmin d , vmax c , vmin c ), respectively, and the number of trucks. The columns A 1, . . . , A 13 reflect the different values of those parameters. The first three experiments A 1 to A 3 vary the density of customers (represented in boldface), maintaining constant values for the rest of the parameters. Experiments A 4 to A 7 vary the number of available trucks. Finally, the experiments A 8 to A 13 vary simultaneously the number of products carried per truck and the number of available trucks. Each experiment must be solved in two stages: First, it is necessary to solve the initial system of non-linear Equations (2) and (3) to determine the value of the width of each circular crown W , the drones and truck speeds, and the length of the truck delivery trip for the outside circular crown. This system of equations is solved to optimality using the dynamic link libraries (DLLs) of the non-linear global solver LINGO v18 [ 44 ]. Second, the Model A MIP instance corresponding to each experiment is coded in Python and solved on an Intel i7-6500U CPU processor with 2.5 GHz and 8 GB of RAM using the Python application program interface (API) of Gurobi 9.1 [45]. 8.2. Results of the Experiments for the Combined Truck–Drone Delivery System Table 3contains the results of the 13 experiments for the combined truck–drone delivery system. The first four rows provide the results of the non-linear system of Equations (2) and (3) for each experiment. Specifically, the drone speed, the truck speed, the number of circular crowns, and the width of each circular crown. The three last rows show the results of solving the trucks assignment problem Model A, providing the system delivery time, the CO 2 emissions, and the number of used trucks. Each column shows the results of the corresponding experiment from A 1 to A 13. Concerning the first three experiments, as expected, the system delivery time, the CO 2 emissions and the number of used trucks increase when the density of customers increases. In experiments A 2 and A 3, corresponding to densities of 5 and 1 customer/km 2 , only 15 and 6 trucks of the initial fleet of 20 trucks are used. The increment in the values of both the system delivery time and the CO 2 emissions is non-linear; the higher the density is, the higher the increment of both variables. Experiments A 4 to A 7 measure the increment of the system delivery time with respect to the number of used trucks. Delivery time decreases as the size of the fleet increases. Again, the decrement in the delivery time is non-linear with the number of used vehicles. The decrease is higher for lower values of the available number of trucks. Since the number and width of the circular crowns remain constant, there is no change concerning CO 2 emissions. Experiments A 8 to A 13 measure the changes in the system delivery time and the CO 2 emissions when modifying both the number of available vehicles and the number of products carried by each truck. 8.3. Definition of Experiments for the Truck-Only Delivery System Concerning the truck-only delivery system, two different sets of experiments have been carried out. Tables 4and 5contain the parameters values for models B1 and B2, respectively. The first set corresponds to the model with the objective of minimizing delivery time. As in the previous case, the different experiments B 1.1 to B1.13 are organized by columns. The first three columns correspond to variations in the customer density. Experiments from B1.4 to B1.7 measure the effect of the variation in the number of available trucks. The rest of the experiments, from B1.8 to B1.13 correspond to the simultaneous variation in the number of products carried by each truck and the number of available trucks. Please note that to facilitate the comparison, for each experiment, the rows corresponding to the parameters Sustainability 2022,14, 13459 20 of 30 truck speed, circular crown width and the number of circular crowns (whose values are written in boldface) are fixed to the values obtained after solving Model A. The second set of experiments, Table 5, corresponds to the parameters values for Model B2. In this case, since the objective consists of minimizing the size of the fleet, the experiments of the second block (experiments to measure the variation of delivery time and emissions with the number of trucks) have been removed. Then, Table 5contains only two blocks of experiments, the first set, B2.1 to B2.3 measuring the impact of varying the customer density and the second one, columns B2.8 to B2.13 measuring the effect of changing the number of products carried by the truck. Please note that in all these experiments, the rows corresponding to the parameters truck speed, circular crowns, number of circular crowns, and delivery system time (whose values are written in boldface) are fixed to the values obtained after solving Model A. As previously mentioned, instances of Model B2 may be infeasible when imposing an upper limit to the system delivery time. For this reason, the upper bound is relaxed, and the optimization tries to respect as much as possible the imposed delivery time while minimizing the fleet size. 8.4. Results of the Experiments for the Truck-Only Delivery System We first analyze the results of Model B1. Table 6summarizes the results of the experiments presented in Table 4. The first row shows the system delivery time, the second row contains the obtained values for CO 2 emissions, and the third row shows the required number of trucks. A first look allows us to conclude that the truck-only delivery system produces higher delivery times in all cases. Concerning the variation in the density of customers, the average increase in the system delivery time reaches a value of 119.18%. When measuring the effect due to the available number of trucks, again the truck-only delivery system obtains worse results concerning the delivery time, obtaining an average increase of 112.77%. The average increase in the delivery time is about 37% for the third block of experiments, those related to the variation in the number of carried products. Similar results are obtained when focusing on the CO 2 emissions, with even higher percentages, 171.14%, 204.79% and 84.28% for the three blocks of experiments. In general, the truck-only delivery system requires more trucks than the combined truck–drone delivery system, especially when the density of customers is high or when the number of available trucks in the combined system is low. For example, in the second block of experiments, with a density of 10 customers/km2the truck-only delivery system always needs 23 trucks, which supposes an average increment of 43.7%. Figure 7shows a graphical comparison of these results. Table 7presents the results of experiments for Model B2. In this situation, the value of the parameter α in (34) has been set to 0.1, thus prioritizing the minimization of the number of trucks. When compared with the combined truck–drone system, the average increase of the system delivery time for the first three experiments (B2.1 to B2.3) is about 341% and of 140% for the rest of experiments (B2.8 to B2.13). As in previous experiments, the delivery time increases with the density and when the number of products carried per truck decreases. Since the objective of Model B2 is the minimization of the number of used trucks, an increase in the system delivery time was expected. On the contrary, the number of trucks has been significantly reduced compared to Model A, where trucks were taken as parameters. To make a fair comparison, we run again Model A but fix the number of trucks to those obtained in the optimal solution of Model B2. The result of these new experiments is reported in Table 8, where each column, denoted as A2 instead of B2, shows the values of the system delivery time, CO 2 emissions, and the number of trucks used in the experiment. Again, the combined truck–drone delivery system outperforms the truck-only system in both the system delivery time and the CO 2 emissions. Concerning the former, the average reduction of total delivery time is of 57.7%, whereas the average reduction in CO2emissions reaches a 115%. Figure 8shows a graphical comparison of these results. Sustainability 2022,14, 13459 21 of 30 Table 2. Experiments for the combined Truck–Drone delivery system. Combined Truck–Drone Delivery System. Model A Density (ρ) Fleet Size Number of Products Parameters Units A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 A13 Density (ρ) customers/km 210 5 1 10 10 10 10 5 5 5 5 5 5 Number of products (P) 50 50 50 50 50 50 50 30 20 10 30 20 10 Drone battery Endurance (E) h 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 Rate Batteries/Drones (ˆ n) 2 2 2 2 2 2 2 2 2 2 2 2 2 Radius (R) km 6 6 6 6 6 6 6 6 6 6 6 6 6 vmax dkm/h 45 45 45 45 45 45 45 45 45 45 45 45 45 vmin dkm/h 30 30 30 30 30 30 30 30 30 30 30 30 30 vmax ckm/h 60 60 60 60 60 60 60 60 60 60 60 60 60 vmin ckm/h 30 30 30 30 30 30 30 30 30 30 30 30 30 Number of trucks 20 20 20 10 15 20 30 20 20 20 10 10 10 Table 3. Results of the experiments for the combined Truck–Drone delivery system. Combined Truck–Drone Delivery System. Model A Density (ρ) Fleet Size Number of Products Variables Units A1 A2 A3 A4 A5 A6 A7 A8 A9 A10 A11 A12 A13 Drone speed km/h 31.28 34.37 37.43 31.29 31.29 31.29 31.29 30.55 31.24 30.07 30.55 31.24 30.07 Truck speed km/h 31.25 33.33 31.25 31.25 31.25 31.25 31.25 30 31.67 30.83 30 31.67 30.83 Number of circular crowns 15 8 3 15 15 15 15 12 19 37 12 19 37 Crown width 0.40 0.75 2.00 0.40 0.40 0.40 0.4 0.50 0.32 0.16 0.50 0.32 0.16 Delivery time h 2.06 1.20 1.18 3.34 2.25 2.06 1.41 1.98 2.13 4.03 3.08 4.14 7.91 CO2emissions kg 351.16 171.40 65.97 351.16 351.16 351.16 351.16 297.69 440.57 825.39 297.69 440.57 825.39 Number of used trucks 20 15 6 10 15 20 30 20 20 20 10 10 10 Sustainability 2022,14, 13459 22 of 30 Table 4. Experiments for the Truck-only delivery system. Equilibrating truck delivery length. Truck-Only Delivery System. Model B1 Density (ρ) Fleet Size Number of Products Parameter Units B1.1 B1.2 B1.3 B1.4 B1.5 B1.6 B1.7 B1.8 B1.9 B1.10 B1.11 B1.12 B1.13 Density (ρ)products/km210 5 1 10 10 10 10 5 5 5 5 5 5 Number of products (P) 50 50 50 50 50 50 50 30 20 10 30 20 10 Radius (R) km 6 6 6 6 6 6 6 6 6 6 6 6 6 Truck speed (vc) km/h 31.25 33.33 31.25 31.25 31.25 31.25 31.25 30.00 31.67 30.83 30.00 31.67 30.83 Circular crown width (W) km 0.40 0.75 2.00 0.40 0.40 0.40 0.4 0.50 0.32 0.16 0.50 0.32 0.16 Number of circular crowns (N)15 8 3 15 15 15 15 12 19 37 12 19 37 Number of trucks 20 20 20 10 15 20 30 20 20 20 10 10 10 Table 5. Experiments for the Truck-only delivery system. Minimizing the number of trucks for given delivery times. Truck-Only Delivery System. Model B2 Density (ρ) Number of Products Parameter Units B2.1 B2.2 B2.3 B2.8 B2.9 B2.10 B2.11 B2.12 B2.13 Density (ρ) products/km 21051555555 Number of products (P) 50 50 50 30 20 10 30 20 10 Radius (R) km 6 6 6 6 6 6 6 6 6 Truck speed (vc) km/h 31.25 33.33 31.25 30.00 31.67 30.83 30.00 31.67 30.83 Circular crown width (W) km 0.40 0.75 2.00 0.50 0.32 0.16 0.50 0.32 0.16 Number of circular crowns (N)15 8 3 12 19 37 12 19 37 System delivery time h 1.7683 1.1931 0.9299 1.9868 2.1349 4.0401 3.0883 4.1458 7.9128 Sustainability 2022,14, 13459 23 of 30 Since the advantage of the combined truck–drone delivery system increases with the density of customers, see Table 9for a comparison of models A and B1 considering densities varying from 50 to 100 customers/km 2 with 100 and 200 products per truck (experiments C 1 to C 12), a final set of experiments with low customer densities and low truck capacities (number of products per truck) has been carried out. For this set, a daily delivery cost evaluation is also performed. To this end, we have considered a per hour truck cost of 2.67 e (supposing a truck life of 8 years and a cost per truck of 60,000 e ), a driver cost per hour of 20 e , a per drone hourly cost of 0.089 e (supposing a drone cost of 1000 e and an useful life of 4 years), regarding the batteries, a cost of 200 e per unit has been considered, supposing a battery useful life of 6 months. Additional data include the gas oil consumption of trucks, fixed to 10 L/100 km with a cost of 0.8 e /L. The three models have been solved for two density values, 5 and three customers per km 2 , by varying for each case the number of products carried from 30 to 10, as depicted in Table 10, where the experiments are denoted as D1, . . . , D10. Concerning the combined truck–drone delivery system, we consider that 5 drones are carried by each truck in all the experiments, while the number of replacement batteries ˆ n is set as “ number o f products/ 5 ” for each experiment. Table 10 shows the result of the three models A , B 1 and B 2. As reported, the daily delivery cost of the combined truck–drone delivery system outperforms the truck-only one in all the experiments, (see also Figure 9a). Concerning delivery times, the models A and B1 obtain quite similar results (see Figure 9b), since Model B1 also minimizes the delivery time, but contrarily increases the number of trucks, producing the most expensive results. (a) (b) Figure 7. Comparison of results of models A and B1. (a) System delivery time. (b) CO2emissions. (a) (b) Figure 8. Comparison of results of models B2 and A. (a) System delivery time. (b) CO2emissions. Sustainability 2022,14, 13459 24 of 30 Table 6. Results for the Truck-only delivery system. Equilibrating truck delivery length. Truck-Only Delivery System. Model B1 Density (ρ) Fleet Size Number of Products Variables Units B1.1 B1.2 B1.3 B1.4 B1.5 B1.6 B1.7 B1.8 B1.9 B1.10 B1.11 B1.12 B1.13 Delivery time h 3.85 3.05 2.57 6.86 4.74 3.85 3.51 2.68 3.14 4.86 4.55 5.72 9.63 CO2emissions kg 1070.32 540.47 127.53 1070.32 1070.32 1070.32 1070.32 641.36 792.05 1190.90 641.36 792.05 1190.90 Fleet size 23 14 3 23 23 23 23 20 18 22 17 14 13 Table 7. Results for the Truck-only delivery system. Minimizing the number of required trucks. Truck-Only Delivery System. Model B2 Density (ρ) Number of Products Variables Units B2.1 B2.2 B2.3 B2.8 B2.9 B2.10 B2.11 B2.12 B2.13 Delivery time h 8.6 5.8 5.0 7.00 7.00 9.70 7.00 7.00 10.00 CO2emissions km/h 1088.76 558.65 138.16 654.95 787.11 1201.32 654.95 787.11 1201.32 Fleet size 8 6 2 7 8 10 7 8 10 Table 8. Results for the combined Truck–Drone delivery system after fixing the values of fleet sizes to the previously obtained for Model B2. Combined Truck–Drone Delivery System. Model A with Fixed Values on the Number of Trucks Density (ρ) Number of Products Variables Units A2.1 A2.2 A2.3 A2.8 A2.9 A2.10 A2.11 A2.12 A2.13 Delivery time h 4.16 2.85 3.22 4.19 5.18 7.90 4.19 5.18 7.90 CO2emissions km/h 350.67 171.16 65.92 297.44 440.33 825.15 297.44 440.33 825.15 Fleet Size 8 6 2 7 8 10 7 8 10 Sustainability 2022,14, 13459 25 of 30 Table 9. Comparison between the combined and the truck-only delivery systems for high-demand scenarios. Minimizing delivery time. Parameters Units C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 Density (ρ) products/km 250 60 70 80 90 100 50 60 70 80 90 100 Number of products (P) 100 100 100 100 100 100 200 200 200 200 200 200 Drone battery Endurance (E) h 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 Rate Batteries/Drones (ˆ n) 5 5 5 5 5 5 10 10 10 10 10 10 Radius (R) km 6 6 6 6 6 6 6 6 6 6 6 6 Data vmax dkm/h 45 45 45 45 45 45 45 45 45 45 45 45 vmin dkm/h 30 30 30 30 30 30 30 30 30 30 30 30 vmax ckm/h 60 60 60 60 60 60 60 60 60 60 60 60 vmin ckm/h 30 30 30 30 30 30 30 30 30 30 30 30 Number of trucks 20 20 20 20 20 20 20 20 20 20 20 20 Variables Combined truck–drone system Drones speed km/h 30.14 30.09 30.06 30.03 30.01 30.18 30.14 30.09 30.06 30.03 30.01 30.18 Trucks speed km/h 31.25 31.25 31.25 31.25 31.25 31.46 31.25 31.25 31.25 31.25 31.25 31.46 Number of circular crowns 75 90 105 120 136 151 75 90 105 120 136 151 Model Crown width km 0.08 0.07 0.06 0.05 0.04 0.04 0.08 0.07 0.06 0.05 0.04 0.04 A Delivery time h 5.89 7.07 8.25 9.43 10.71 11.80 4.89 5.85 6.82 7.79 8.83 9.75 CO2emissions kg 1385.87 1663.04 1940.21 2217.38 2516.52 2788.45 1221.95 1466.34 1710.73 1955.12 2215.80 2460.19 Number of used trucks 20 20 20 20 20 20 20 20 20 20 20 20 Variables Truck-only system Delivery time h 14.95 18.48 21.56 23.91 25.44 29.91 13.92 17.25 20.12 22.28 23.60 27.86 Model CO2emissions kg 5093.31 6334.48 7390.27 8149.46 8562.02 10,253.03 4929.40 6137.79 7160.80 7887.20 8266.54 9924.77 B1 Number of used trucks 20 20 20 20 20 20 20 20 20 20 20 20