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

: 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 deﬁning the main characteristics of trucks and drones, a sectorization of the delivery area is ﬁrst 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


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.
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 2 reviews the related work. Section 3 presents the description of our problem. Section 4 establishes the set of assumptions that have been considered. Section 5 remarks the main contributions of this research. Section 5 presents 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 6 contains 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 7 contains 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 8 presents several conclusions.

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][7][8][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][14][15][16], and a generalization of vehicle-routing problems to include drones (VRP-D) [17][18][19][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 literatureexcept for only a few works, see, e.g., [25][26][27][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 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.
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.

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 n d and a given number of replacement fully charge batteries n b . As we will see later, for efficiency reasons, the ration between the number of replacement batteries n b and the number of drones n d 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 ith 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 + S 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.  Table 1 summarizes 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 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.  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 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].

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

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

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 + S i with speed v d must be less or equal to the time spent by the truck to travel the distance S with speed v c . Moreover, since the distance S 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-. Consider two concentric circumferences, with radius R and r, R > r and a chord C of the same length at both circumferences. As illustrated in Figure 4a, the sagitta of the inner circumference f will be larger than the sagitta of the outer circumference f . Thus, since the curvature is greater in the inner circumference, the length of the arc corresponding to the chord C is larger in the inner circumference. Then, if we consider two arcs of the same length at two concentric circumferences, the length of the chord in the inner one is shorter than in the outer circumference. A mathematical demonstration is easy if we consider R = 2r or R = 3r. Suppose that R = 2r, see Figure 4b.  R θ = r θ = S The length of chords d and d can be expressed as: Demonstrating that d is greater than d is equivalent to demonstrate that d /d < 1.
which is <1 except for θ = 0 and θ = 2π, situations that do not make any sense in our case.
A similar reasoning can be followed for the case R = 3r using the formula of the triple angle. Obtaining a general expression for R = ar with a > 1 is not easy. For interested readers, Table A1 in the Appendix A resumes the values of d /d for angles varying from 15 to 360 degrees and for relations r/R from 1/2 to 9/10. To summarize, for a given truck trip length S the length of the chord increases with the radius, the total length of the drone flight x + z (Figure 4b) also increases with the radius.
Thus, to determine the maximum trajectory length of the truck, the circular outer crown must be chosen. - Determining the width of circular crowns and the circular truck trajectory length Taking as reference the outside circular crown, the next step consists of determining the width of circular crowns W, the number of crowns N and the maximum circular truck length S to ensure the synchronization of both vehicles for the worst case.
The time required by the drone to travel the distance W/2 + S 1 (which at the worst case must be equal to the drone flight endurance E -measured in min-) must be lees or equal to the time spent by the truck to travel a distance S, but, for a precise synchronization, thinking in moving the truck as quick as possible and taking into account that we are considering the worst case and that, as later reported, due to a final rounding procedure the length of truck trips will be later shortened (see the end of this section), we can consider: where, according to Figure 5, Since d 1 is the base of an equilateral triangle (see Figure 5), we can obtain d 1 as half of the triangle side opposite to the angle θ 1 : then, using the sin theorem, , the next expression holds: By substituting d1 as a function of W and θ 1 , we have: By other hand, if a truck transports the same number of products as drones (each drone with only one battery-no replacements-) the delivery area A covered by a truck at circular crown 1 can be expressed as a function of the angle θ 1 , R and W (later expressed as a function of the customer density): Figure 5. Drone trajectories at two concentric circumferences.
As explained in Section 3, each truck is carrying P products, n d drones, and n b fully charged batteries. If n d = n b (n = 1) each drone can perform only one trip (one delivery per truck trip). Ifn = 2 each drone can make two trips per truck trip, and then, the total length traveled by the truck will be 2S. Then, the total delivery area will be 2A. Obviously, for efficiency reasons, the number of products to be delivered by a truck must be equal to the number of replacement batteries; otherwise, some products could not be delivered (there are not enough replenishment batteries), or some of the carried batteries will not be used.
For a similar reason, it does not make to carry several replacement batteries lower than a multiple of the number of drones.
Thus, if the truck carriesn batteries and P products, the delivery area, the customer density and the number of products carried are related by: The system of Equations (2) and (3): allows us to numerically determine W and θ 1 , as well as S, for a given set of values of parameters P, E, ρ,n and R.
Please note that the value of S = θ 1 R − W 2 , corresponding to the outer circular crown, is obtained for an endurance exactly equal to E, but as we will see next, in the majority of cases S is diminished to ensure an integer number of truck delivery cycles, thus guaranteeing that the drone autonomy is enough to perform the worst-case delivery. - Determining number of circular crowns and rounding to integer values The number of circular crowns can be computed as: which turns in valuesŴ = R N ≤ W, and the number of truck delivery cycles for the ith circular crown will be: which ensures values of truck circular trip lengthsŜ i ≤ S for each circular crown.

Assigning Trucks to Crowns and Equilibrating the Distance Traveled by Trucks-Model A
After computing the number of circular crowns N and the number of cycles of trucks per crown C i , assuming a fleet size of T trucks, the assignment of trucks to delivery cycles can be obtained by solving the next MIP model, where x ij is an integer variable representing the number of cycles of circular crown i that are assigned to truck j = {1, . . . , T}. Please note that the length of a cycle corresponding to the circular crown ith, l i can be obtained as: where the last addend collects twice the distance from the circular crown central line to the depot.

Min L (9)
subject to : x ij ∈ N + The objective function (9) minimizes the maximum length traveled by any truck. Constraints (10) 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 (11) compute the total length traveled by truck j as a function of variables x ij . Finally, the constraints (12) ensure that all the cycles of each circular crown are assigned to trucks. As a result, the distance traveled by all the trucks is equilibrated.

Key Performance Indicators for the Combined Truck-Drone Delivery System
In this section, we develop two performance key indicators that later will be used for comparing width the truck-only delivery to customers.
Moreover, at the endpoint of each truck circular movement (inside each cycle), the truck must wait for the return of all the launched drones. The worst-case corresponds to a drone launched just when the truck reaches the endpoint of its trajectory and the customer is just located at the border of the circular crown. In this case, the distance traveled by the drone is W and the waiting time W v d . Summarizing, the truck delivery time for a cycle at circular crown i is: Using the optimal values of variables x * ij of Model A ((9)-(13)), the total delivery time spent by truck j is obtained as: and the maximum system delivery time, denoted by DT, corresponds to the maximum of these values: To compare the results of the combined truck-drone delivery system with those obtained using only trucks, two additional key performance indicators are defined. Let lg be the number of fuel liters consumed by a truck to travel 1 km at speed v c , bc the number of kWh needed for a full charge of one battery, κ a factor expressing the number of CO 2 (Kg) emitted by a truck when burning 1 L of fuel and φ the average number of CO 2 (Kg) emitted to the atmosphere when producing 1 kW of electricity.
The CO 2 emissions of trucks T CO 2 are directly related to the length of the traveled distance: where both T * and L * j are the optimal values of Model A ((9)-(13)). For computing the CO 2 emissions due to the drones, represented by D CO 2 , we will consider the worst case, supposing that the batteries are full discharged at each trip.
Then, the total pessimistic CO 2 emissions are:

The Truck-Only Delivery System
To test the efficiency of the combined truck-drone delivery system, two models with the same structure are proposed, but delivering products only with the use of trucks. We consider trucks to have the same capacity (number of products) as in the combined delivery system. The delivery area considered is also the same. Furthermore, to make the comparison as objective as possible, the distribution zone will be divided into circular crowns as in the combined model, respecting the amplitude W of the previously calculated circular crowns, but adapting the length of each circular truck trip. Please note that for given values of customer density and truck capacity, the number of delivery zones is the same if we use drones or not, which is determined by In this case, the dynamics of the distribution process will be as follows: The trucks leave the depot loaded with the products and travel up to some point of the central circumference of the corresponding circular crown. From there, each truck can perform to types of movements: a circular movement on the circular route (as in the combined truck-drone system) or radial movements towards customers to deliver the corresponding parcel. After delivering all the products, the trucks will return to the depot to be loaded again before moving on to a new route. Figure 6 shows a graphical representation of the truck-only delivery system. Importantly, as explained in the assumptions, we are interested (as in the combined case) in defining a constant set of daily main routes the trucks must follow independently of the exact customer location. i.e., the trucks follow a fixed main route and then, every day, deliver goods to different customers located inside each zone at a priori unknown positions. In what follows, we derive the length of circular movements and the angle of each sector zone. Since each truck is carrying P products, the delivery area covered by a truck at circular crown ith must be equal to P/ρ.
Then, the angle θ i corresponding to the delivery area can be obtained as: The length of the trip performed by each truck at circular crown ith will be: Thus, the number of cycles corresponding to the circular crown ith is obtained by rounding to the upper nearest integer the quotient between the length of the central circumference and the length of the trip l i .
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 θ i is also updated.
The number of products carried by each truck is now: And the total length of each delivery cycle at circular crown ith will be: Where the first term takes into account the average distance corresponding to the radial movements from the central line of circular crown ith 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.

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 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 ofT trucks, the assignment of trucks to cycles can be obtained by solving the next MIP model, where x ij 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 : x ij ∈ N + (32) The objective function (28) minimizes the maximum length traveled by any truck, L j , 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 x ij . 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 j is:

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 x ij represent the number of delivery cycles for the circular crown ith assigned to vehicle jth. Since the number of vehicles must be minimized, for each vehicle, a binary variable δ j , taking value 1 if truck jth is used for the delivery of goods, is included. Denoting byĈ 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,Ĉ = ∑ N i=1 C i , the MIP model is formulated as follows: subject to : 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, H j 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 jth as a function of variables x ij . Constraints (37) express the total delivery time of truck jth. Constraints (38) ensure that all the cycles of each circular crown are assigned to trucks. Constraints (39) are used to relate variables x ij and δ j . If truck jth 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.

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.  40)), as appropriate, the total system delivery time DT will be: where NT =T in the case of Model B1 and NT = ∑Ĉ j=1 δ * j when using Model B2.
(b) Co 2 emissions Using the optimal values of variables x * ij in the Model B1 ((28)-(32)) or Model B2 ((34)-(40)), as appropriate, the total system CO 2 emissions are: where NT =T in the case of Model B1 and NT = ∑Ĉ j=1 δ * j when using Model B2. Table 2 summarizes 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 (v max d , v min d , v max c , v min c ), respectively, and the number of trucks. The columns A1, . . . , A13 reflect the different values of those parameters. The first three experiments A1 to A3 vary the density of customers (represented in boldface), maintaining constant values for the rest of the parameters. Experiments A4 to A7 vary the number of available trucks. Finally, the experiments A8 to A13 vary simultaneously the number of products carried per truck and the number of available trucks.

Definition of Experiments for the Combined Truck-Drone Delivery System
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]. Table 3 contains 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 A1 to A13. 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 A2 and A3, 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 A4 to A7 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 A8 to A13 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.

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 4 and 5 contain 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 B1.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 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 5 contains 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.

Results of the Experiments for the Truck-Only Delivery System
We first analyze the results of Model B1. Table 6 summarizes 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/km 2 the truck-only delivery system always needs 23 trucks, which supposes an average increment of 43.7%. Figure 7 shows a graphical comparison of these results. Table 7 presents 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 CO 2 emissions reaches a 115%. Figure 8 shows a graphical comparison of these results.   Since the advantage of the combined truck-drone delivery system increases with the density of customers, see Table 9 for 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 C1 to C12), 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.67e (supposing a truck life of 8 years and a cost per truck of 60,000e), a driver cost per hour of 20e, a per drone hourly cost of 0.089e (supposing a drone cost of 1000e and an useful life of 4 years), regarding the batteries, a cost of 200e 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.8e/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 batteriesn is set as "number o f products/5" for each experiment. Table 10 shows the result of the three models A, B1 and B2. 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.

Conclusions
We have discussed the design of a circular hybrid truck-drone last-mile delivery for an urban context, where the delivery area is characterized by a certain density of customers. The need for synchronization between drones and trucks influences the definition of the design of the system. The circular delivery system divides the city into a set of circular crowns whose width depends on the drones endurance and other system parameters, including capacity of trucks and average speeds of drones and trucks.
After completing the definition of the delivery area, a first MIP truck-drone delivery model is proposed to assign trucks to delivery cycles while equilibrating the length traveled by trucks, thus minimizing the system delivery time (Section 6.2). Two additional truckonly MIP models have been developed to compare the results obtained against a delivery procedure using only trucks. The first one (see Section 7.1) minimizes the truck-only system delivery time, fixing or not the number of vehicles, as appropriate for comparisons, and the second one (see Section 7.2) minimizes the number of used trucks, considering as a goal the delivery times obtained after solving the truck-drone MIP model.
Several sets of experiments have been carried out. First, to compare the truck-drone and the truck-only (minimizing delivery time) models, 13 experiments which imply the resolution of 13 non-linear systems of equations (to determine the width and number of the circular crowns as well as the speeds of vehicles) and 39 MIP models (13 of each type). Second, to compare the truck-drone and the truck-only (minimizing the number of trucks) models, 18 experiments are carried out, solving nine truck-drone problems to minimize the delivery time with a fixed fleet of trucks and nine truck-only problems minimizing the fleet of trucks. Finally, a set of ten experiments are developed to compare the truck-drone and the truck-only systems in the case of low truck capacities and low customer density, solving ten non-linear systems of equations for determining the main design area parameters plus 30 MIP models, ten of each type. This set of experiments also proposes a daily cost assessment for both delivery systems.
In all the experiments, the truck-drone combined delivery system outperforms the results of the truck-only system for the different defined key performance indicators: daily system delivery time and CO 2 emissions. As expected, the advantage of the combined system increases with the density of customers, which respond to the motivation of this research. According to the last set of experiments, the truck-only delivery system results more expensive even for low densities of customers/km 2 and relatively low capacity of trucks.
In a whole, we can state that our methodology is appropriate for setting a repetitive (day-to-day) delivery process on truck-drone usage to service the huge last-mile delivery demand arisen in dense urban areas. Although the number of delivery orders and the location of customers vary every day, our assumption is that a great number of customers are spread over a broad service area and that the density of customers remains approximately constant. The applied continuous approximation technique allows us to derive analytical formulas to determine the best area sectorization (expressing the division of the circular area as a function of the truck capacity along with the speed and endurance of the drones).
Furthermore, a mixed-integer programming model allows decision makers to define the set of daily routes for a fleet of homogeneous trucks, each of which is assigned to service a certain number of sectors, carrying a certain number of drones. Aside from the typical limited capacity for carrying product exhibited by each truck, we also include other real-life operational assumptions; a set of replacement batteries is available for each drone at each truck, the trucks start and end their routes at a central depot where they are filled with new products and new sets of batteries before being ready to serve again the crowns.
Arguably, our continuum approximation approach leads to average rather than a detailed description of results. However, it gives us the possibility to study these complex problems in a simplified way to obtain a general understanding of the performance and convenience of the system. We derive a sort of quantitative measures such as delivery times, truck CO 2 emissions, and total length of delivery routes, which demonstrated that the combined truck-drone system outperforms the truck-only system.
Although several legal/regulatory questions (concerning the use of drones in civil areas) and other issues related to the final customer delivery infrastructure remain, the analysis developed in this work points out, in a quantitative way, the advantage of using UAVs to perform last-mile delivery activities in broad delivery areas.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.