University of Southern California Dissertations and Theses

Package delivery with trucks and UAVs

(USC Thesis Other)

Package delivery with trucks and UAVs

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Package Delivery With Trucks And UAVs by Siyuan Song A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Industrial Engineering) August 2018 Copyright 2018 Siyuan Song Acknowledgements Like many others, my road to Ph.D is full of doubt and hesitation. I can't remember how many times during my eight years, the idea of giving up will occur. But because of people around me, I come to where I'm now. Without any of them, my Ph.D life might end up stopping at somewhere in the past. I would rst like to give thanks to my adviser, Professor John Carlsson. As what I always tell my friends, the wisest decision I have made in recent years is to be John's student. During my Ph.D years, John not only fully supports me nancially, but also always impresses me with his innovative ideas, brilliant solutions and beautiful writings. I have learned a lot from him. Every time I'm getting stuck and feeling depressed, meeting with him will bring out hopes. Besides these, John is always supportive like a friend in many aspects. My Ph.D years become much easier with John and I really appreciate it. I would also express my great gratitude towards my committee members Professor Ketan Savla and Phebe Vayanos. They are also in my committee of qualifying exam, their time and eorts to support and advice me at those crucial steps towards my Ph.D are of great value to me. Lots of other professors have also given me important helps. Especially for Professor Shouhong Wang and Elizabeth Houseworth at Indiana University, Professor William Cooper at University of Minnesota, Professor Sheldon Ross and Meisam Razaviyayn at University of Southern California, I give my sincere thanks. My current and former lab-mates, Fan Jia, Mehdi Behroozi, Ying Li, Xiangfei Meng, Ye Wang, Jiachuan Chen and Bo Jones, have also given me great supports. Although with lots of ii challenges in courses, research, job hunting and etc., it's good to have colleagues like you facing those together. Last but not the least, thanks to my girlfriend Yuexin Xu who for so many years, has accom- panied me when I'm humble, encouraged me when I'm depressed and rejoiced with me when I'm joyful. Also thanks to my parents who are always supportive behind me. iii Table of Contents Acknowledgements ii List Of Tables vi List Of Figures vii Abstract ix Chapter 1: Introduction 1 1.1 Related work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.1 Traveling salesman problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.2 Generalized travelling salesman problem . . . . . . . . . . . . . . . . . . . . 7 1.1.3 UAV in logistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.1.4 Vehicle routing problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.1.5 Backbone network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2 Notational conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Chapter 2: Preliminaries 13 Chapter 3: Continuous Horse y Problem 16 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2 Integer programming formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.3 Asymptotic analysis under uniform distribution . . . . . . . . . . . . . . . . . . . . 19 3.3.1 Asymptotic lower bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.3.2 Asymptotic upper bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.4 Asymptotic analysis under general density function . . . . . . . . . . . . . . . . . . 28 3.5 Cost as objective function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.6 Computational results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.6.1 (r) for time of completion as objective . . . . . . . . . . . . . . . . . . . . 38 3.6.2 (r) for total cost as objective . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.6.3 Experiments using road network distances . . . . . . . . . . . . . . . . . . . 45 Chapter 4: Horse y Problem with Multiple UAVs 48 4.1 Theoretical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.2 Computational Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.2.1 Experiments using Euclidean distances . . . . . . . . . . . . . . . . . . . . . 52 Chapter 5: Horse y Problem with Multiple Objectives 56 5.1 Theoretical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 5.2 Computational Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 iv Chapter 6: Horse y Problem with Battery Life Constraints 63 6.1 Theoretical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 6.2 Computational Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 Chapter 7: Horse y Problem with Truck Unreachable Regions 71 7.1 Theoretical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 7.2 Computational Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 7.2.1 Time of completion estimation . . . . . . . . . . . . . . . . . . . . . . . . . 76 Chapter 8: Discrete Horse y Problem 80 8.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 8.1.1 Asymmetrical traveling salesman problem formulation . . . . . . . . . . . . 81 8.2 Theoretical Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 8.3 Computational Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 8.3.1 Algorithm 1: Directly Solving AGTSP . . . . . . . . . . . . . . . . . . . . . 84 8.3.2 Algorithm 2 : solving AGTSP with launching sites pre-selection . . . . . . 86 8.3.3 Algorithm 3: adjusting solutions of continous horse y problem . . . . . . . 86 8.3.4 Optimal value with respect to dierent K . . . . . . . . . . . . . . . . . . . 87 Chapter 9: Future Work 90 9.1 Horse y problem with charging time . . . . . . . . . . . . . . . . . . . . . . . . . . 90 9.2 Uncertainty of trac time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 9.3 Better heuristics for horse y problem . . . . . . . . . . . . . . . . . . . . . . . . . . 92 9.4 Multiple robots problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 9.5 Pareto frontier for multiple objectives . . . . . . . . . . . . . . . . . . . . . . . . . 93 Reference List 95 Appendix A Proof of Lemma 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 v List Of Tables 3.1 (r) values with dierent n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.1 Lower bound of (r;k) for dierent k . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.2 Estimate of ( 1 ;k) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.1 Optimal values comparisons with dierent rations for multiple objectives problem . 61 7.1 Algorithms comparison with and without projections . . . . . . . . . . . . . . . . . 76 7.2 Comparisons between numerical results and theoretical bounds . . . . . . . . . . . 78 vi List Of Figures 1.1 Comparisons between dierent delivery systems . . . . . . . . . . . . . . . . . . . . 3 2.1 Illustration of neighborhoods for lemma 2 . . . . . . . . . . . . . . . . . . . . . . . 15 3.1 Solutions to horse y problem with dierent UAV speed . . . . . . . . . . . . . . . 17 3.2 An example of griding an area R . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.3 An example of lower bound for horse y problem . . . . . . . . . . . . . . . . . . . 22 3.4 Two lower bounds of horse y problem . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.5 An example of the approximation algorithm . . . . . . . . . . . . . . . . . . . . . . 27 3.6 Example of creating sub-tour from solution of horse y problem . . . . . . . . . . . 31 3.7 Cases for subdividing horse y problem . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.8 Horse y problem comparison between waiting time and cost as objective function . 34 3.9 Partition horse y problem with cost as objective . . . . . . . . . . . . . . . . . . . 37 3.10 Exact solutions for horse y problem of small size . . . . . . . . . . . . . . . . . . . 39 3.11 Heuristic algorithms for solving horse y problem . . . . . . . . . . . . . . . . . . . 40 3.12 (r) values with respect to dierent r . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.13 (r) values with dierent n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.14 (r) values with respect to dierent r (cost version) . . . . . . . . . . . . . . . . . 44 3.15 Experiments of road network (Pasadena) . . . . . . . . . . . . . . . . . . . . . . . . 46 3.16 Pasadena results histogram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.1 Examples of solution for problem with multiple UAVs . . . . . . . . . . . . . . . . 54 vii 5.1 Computational results for two objectives . . . . . . . . . . . . . . . . . . . . . . . . 60 6.1 Comparisons between problems with and without battery life constriants . . . . . 66 6.2 Optimal value changes with respect to battery life constriants . . . . . . . . . . . . 68 6.3 Optimal value changes with respect to dierent parameters for horse y problem with battery life constriants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 7.1 An example of horse y problem with areas that are unreachable for trucks . . . . . 72 7.2 Upperbound scheme illustration with truck unreachable regions . . . . . . . . . . . 74 7.3 Examples of solutions for problems with truck unreachable regions . . . . . . . . . 77 7.4 Estimation for problems with truck unreachable regions . . . . . . . . . . . . . . . 79 8.1 Illustration of extra time caused by moving launching sites . . . . . . . . . . . . . . 83 8.2 Example of solution of discrete horse y problem by solving AGTSP directly . . . . 85 8.3 Comparisons between algorithms solving discrete horse y problem . . . . . . . . . 88 8.4 Optimal values of discrete horse y problem with respect to dierent launching sites 89 9.1 Pareto frontier example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 A.1 Re ecting the polygonal chain s about the line joining v 1 and v 2 . . . . . . . . . . . 101 viii Abstract Unmanned aerial vehicl (UAV) has seen its development recently in lots of areas, both in military and commercial areas. This paper is on its applicatoin of package delivery where UAVs, provides service to customers while making return trips to a truck that is itself moving. In other words, a UAV picks up a package from the truck (which continues on its route), and after delivering the package, the UAV returns to the truck to pick up the next package. In this paper, we rst formulate this problem called horse y routing problem and then give both asymptotic analysis and heuristic algortihms. Then some variants of it are considered. Those problems are all considering some extra factors under some circumstances in practical uses. Chapter 4 considers the problem with more than one UAV in the system and chapter 5 works on a problem with both cost and time of completion as objectives. Chapter 6 considers the fact that nowadays UAV only has limited battery and handles with the horse y routing problem with battery life constriants. Chapter 7 is on problem with region restriction on truck, it happens when taking the adavantage of UAV that reach to some rural areas without roads. Finally, chapter 8 solves a problem that launching sites are restricted to some pre-dened locations. ix Chapter 1 Introduction If a salesman, starting from his home city, is to visit exactly once each city on a given list and then return home, it is plausible for him to select the order in which he visits the cities so that the total of the distance traveled in his tour is as small as possible (gure 1.1a). This is traveling salesman problem (TSP) and It was rst considered mathematically in the 1930s by Merrill Flood who was looking to solve a school bus routing problem [65]. Although with simplicity, it is one of the NP-hard problems and is extremely dicult to solve. The problem is practically useful in logistics where deliveries are made from depot to customs. And this kind of delivery is called last-mile delivery. Last-mile delivery is one of the most important problems in transportation science. Transport- ing goods via freight rail networks and container ships is often a rather ecient and cost-eective manner of shipping. However, when goods arrive at a high-capacity freight station or port, they must then be transported to their nal destination. This last leg of the supply chain is often less ecient, comprising up to 28% of the total cost to move goods. This has become known as the "last mile problem" [87, 11]. The last mile problem can also include the challenge of making deliveries in urban areas where retail stores, restaurants, and other merchants in a central business district often contribute to congestion and safety problems [83, 11]. 1 One of the most talked-about developments in transportation and logistics in recent years has been the potential use of unmanned aerial vehicles (UAVs), or \drones", for transporting packages, food, medicine, and other goods. The most famous proof-of-concept of such a service is the \Amazon Prime Air" system, which was introduced in late 2013 and has since undergone several iterations [18]. Other similar systems include Google's \Project Wing" [58], DHL's \Parcelcopter" [21], and a joint eort between the Swiss Post, Swiss Worldcargo (the air freight division of Swiss International Air Lines), and the California-based startup Matternet [32]. As described in [56], there are many reasons to be optimistic about the role that UAVs can play in next-generation logistics systems: [Demand for last-mile delivery] is likely to increase as e-commerce volumes grow. ... UAVs could provide major relief for inner cities, taking trac o the roads and into the skies. So far, payloads are limited but a network of UAVs could nevertheless support rst and last-mile logistics networks. ... This urban rst and last-mile use case is probably the most tangible and spectacular in the logistics industry. UAVs are also already being used in many related industries including energy [90], agriculture and forestry [3], environmental protection [70], and emergency response [69]. From a transportation scientist's perspective, many of the benets of a UAV-based delivery system are obvious: UAVs have a low per-mile cost, can operate without human intervention, and can travel at high speeds while being unaected by road trac. The shortcomings of such a system are equally apparent: UAVs have an extremely low carrying capacity and short travelling radius, both of which necessitate frequent returns to a central depot. Thus, as suggested in Figures 1.1a and 1.1b, a conventional truck delivery system benets from an economy of scale and suers from high per-mile costs, whereas a UAV delivery system benets from low per-mile costs and lacks an economy of scale. 2 (a) Truck only (b) UAV only (c) Truck + UAV Figure 1.1: (a) shows a travelling salesman tour of a set of client destinations and a central depot, which corresponds to the optimal provision of service by a delivery scheme that uses only a truck. Similarly, (b) shows a collection of direct trips between the same set of points and the central depot, which corresponds to the optimal provision of service by a delivery scheme that uses only UAVs. Finally, (c) shows the solution to a problem in which we have a truck that follows a route that is much shorter than that of (a), as indicated by the thick tour, as well as a single UAV that alternates between visiting the client destinations and and the truck, as indicated by the dashed lines. The routes shown correspond to the case where the UAV is about four times as fast as the truck; it is also worth mentioning that the optimal sequence of the client destinations is not necessarily the same as that induced by the TSP tour which is shown in (a), and indeed one can readily observe that there are several places where the two sequences dier. The purpose of this paper is to determine the eciency of a hybrid approach in which a UAV (or a eet of UAVs) provides service to customers while making return trips to a truck that is itself moving, as illustrated in Figure 1.1c. In other words, a UAV picks up a package from the truck (which continues on its route), and after delivering the package, the UAV returns to the truck to pick up the next package. Although the hardware for such systems has already been constructed { one particular implementation is called the \HorseFly" and was developed jointly with AMP Electric Vehicles and the University of Cincinnati [96] { it is not yet understood to what extent such an approach can actually provide a signicantly improved quality of service. Our goal is to describe precisely how much improvement can be realized in a mathematically sound way. Obviously, the \horse y" problem is extremely dicult to solve to optimality because it is a generalization of the travelling salesman problem (TSP) that requires the consideration of the locations where the truck and the UAV meet. Thus, we approach this problem from asymptotic 3 aspect where analyze the behavior of optimality when number of customs goes to innity. This paradigm is common in transportation researches [91, 92]. Such analysis exists for many combina- torial problems, such as the traveling salesman problem [16, 43], facility location [51, 57, 77], and any subadditive Euclidean functional such as a minimum spanning tree, Steiner tree, or matching [82, 91, 92]. We primarily focus on the time to completion of all of the deliveries, and our results are also easily applicable to other cost functions such as the net energy expenditure over both vehicles. In addition to the eciency analysis for horse y system, some aspects that may be big concerns in practice are also discussed including how UAV's battery life can aect the eciency and the strategy to balance the total cost and time of completion in such system. 1.1 Related work 1.1.1 Traveling salesman problem TSP attracts lots of researchers and tons of papers were produced from dierent perspectives. We mention some of them that either intrigues our ideas or relates with some materials in this paper which include Mathematical formulation Approximation algorithms Probabilistic Analysis Heuristics and empirical Analysis The most famous integer programming formulation for TSP was put forward by Miller,Tucker & Zemlin [22] based on an additional loop elimination constraints. We will use this formulation for exact solutions of TSP in later sections and thus state it here 4 minimize n X i=1 n X j=1 c ij x ij (1.1) X i x ij = 18j (1.2) X j x ij = 18i (1.3) u i u j +nx ij n 1; 2i;jn (1.4) u i 2Z; x ij 2f0; 1g (1.5) where c ij is the cost between two customs location and equations (4) are the loop elimination constraints. TSP is denitely in the category of NP-hard. To make it worse, when the c ij above can be arbitrary, there will be no approximation algorithm ifP6=NP [65]. However, when triangular inequality is satised and the problem becomes metric TSP, the algorithm based on minimum spanning tree can give guarantee of at most twice of the optimal value. The best approximation algorithm of guarantee was from Christodes [28] which gives at most optimal value multiplied by 3=2 with polynomial time complexity. It combines the minimum spanning tree with a solution of another problem, minimum-weight perfect matching. The problem we are working on is Euclidean TSP where distance is dened by ` 2 norm. Thus it's a special case of metric TSP. But even for the Euclidean TSP, it's NP-hard [76]. The Held- Karp heuristic is conjectured to have an approximation ratio 4=3 (some results of Goemans [49] support this conjecture) but the best upperbound known is 3=2 (Shmoys and Williamson [34]). Another aspect is to nd Polynomial Time Approximation Scheme (PTAS) which is a aglortihm for each xed c> 1, can approximate the problem within a factor 1 + 1=c. And the running time could depend upon c which has to be polynomial in the input size. PTASs are known for very few problems; two important ones are Subset-Sum (Ibarra and Kim [73]) and Bin-Packing (Fernandez de la Vega and Lueker [95]). 5 Arora, Lund and etc. [84] showed that ifP6=NP, then metric TSP and many other problems do not have a PTAS. However, Arora [13] showed that Euclidean TSP has a PTAS. For every xed c > 1, a randomized version of this algorithm computes a (1 + 1=c)-approximation to the optimal tour in O(n(logn) O (c)) time. Researches on heuristic algorithms of TSP are enormous. Most of them with best performance are based on local search where an initial tour is rst constructed and then tour improvements methods are implemented. For example, Nearest-neighbor and insertion algorithms were discussed in [72, 4] for initial tour construction step. For tour improvements step, 2-opt and 3-opt where two or three lines were swapped to reach a better tour are the most simple algorithms. Lin-Kernigan heuristic [66] where the number of edges swapped are dynamically determined at each step are the best known algorithm in tour improvement. It will be used repetitively in later sections. In addition to strict improvement process, simulated annealing [60], derived from an analogy with a material annealing process used in mechanics, allow the objective to deteriorate in order to avoid local minimum. It's also one of the popular algorithms in solving TSP. One variant of the TSP that is conceptually related to the horse y problem is the covering salesman problem (CSP) as introduced by Current and Schilling [33]. The CSP aims to nd the shortest tour of a subset of given nodes such that every node that is not on the tour is within a predened covering distance of a node on the tour. Current and Schilling [33] propose a heuristic solution procedure that is based on the set covering problem. Several generalizations and special cases of the CSP have been studied in the literature, see e.g. Gulczynski et al. [36] and Golden et al. [20]. Similar to the CSP, in the horse y problem, the truck does not have to visit all nodes. However, it is not enough to pass nodes that are not visited by the truck within a certain distance threshold, but drone visits have to be scheduled for each such node, and truck and drone need to be synchronized. 6 1.1.2 Generalized travelling salesman problem We also need to solve generalized travelling salesman problem (GTSP) when trying to handle with discrete horse y problem. Thus some work on GTSP will be introudced here. The GTSP is an extension of the (TSP) where the set of cities is partitioned into clusters, and the salesman has to visit every cluster exactly once. The GTSP coincides with the TSP when all clusters are singletons. Like TSP, if the distance matrix C = (c ij ) is symmetric which means c ij =c ji , then it's symmetric GTSP otherwise called asymmetic GTSP. The generalized TSP is a reasonably well-known combinatorial problem and was rst put for- ward by Srivastava, Kumar and etc. [85]. There are exact algorithms such as branch- and-bound and branch-and-cut described in [68]. Laporte and Nobert [45] developed an exact algorithm for GTSP by formulating an integer programming and nding the shortest Hamiltonian circuit through some clusters of nodes. Since the problem is an extension of TSP, a natural way to solve the problem is to transform GTSP to TSP. Transformation of a GTSP into TSP was rst introduced by Lien and Ma [99], in which the number of nodes of the transformed TSP was rela- tively very large which are more than three times the number of nodes in the associated GTSP. Later, Dimitrijevic and Saric [39] developed another transformation to decrease the size of the corr esponding TSP. In their method, the number of nodes of the TSP was twice the number of nodes of the original GTSP. Behzad and Modarres put foward another method in [17] where the number of nodes of TSP is equal to that of the original GTSP. Heuristic TSP solvers are often considered inappropriate for GTSP ([63, 59]). However, Hels- gaun shows that this need not be the case if the state-of-the-art heuristic TSP solver LKH is used [55]. And we are using this method when solving GTSP in later chapters. 7 1.1.3 UAV in logistics The use of UAVs in logistics is a prospect that is very much in its infancy, and as such, there is little research on the economics of such systems. The recent paper by Murry and Chu [71] describes mixed integer programming formulations and heuristic algorithms for solving "Flying Sidekick Traveling Salesman Problem" (FSTSP) in which a truck and a drone perform service in tandem; they are able to nd optimal solutions for up to 20 customer locations and heuristic solutions for up to 75 customers locations. Ponza [80] extended the work of Murray and Chu [71] to solve the FSTSP in his master's thesis, proposing an enhancement to the MIP model and solving the problem by a heuristic method based on Simulated Annealing. Wang and ec al. derives a number of worst-case results for the vehicle routing problem with drones, in which a eet of trucks equipped with drones deliver packages to customers [79]. Agatz, Bouman and ec al. put forward the traveling salesman problem with drones (TSP-D) where drone and truck also together deliver packages and formulate it as an MIP model and develop several fast route rst-cluster second heuristics based on local search and dynamic programming. They also proved worst-case approximation ratios for the heuristics and test their performance by comparing the solutions to the optimal solutions for small instances [10]. All problems above focus on quality level (time). Ha and et al. put forward a problem where the objective is to minimize the total transportation cost with formulation and two algorithms followed. [81] The most concrete discussion of the economic desirability thus far can be found in [37], which analyzes per-mile costs based on the estimated energy consumption of UAVs and describes several feasibility issues in implementation; the author concludes that \From a cost perspective, the numbers do not look unreasonable". Within the operations research community, extensive research has been conducted on the use of UAVs for military applications; for example, [42] uses robust optimization to plan reconnaissance missions subject to uncertainty in fuel usage between locations and weather conditions, and [62] uses scenario-based optimization to decide the number of UAVs 8 to deploy; once deployed, more information about the environment is learned, at which point the route must be determined. Another important application is in healthcare for timely delivery of urgently needed medications, blood and vaccines [88, 8]. 1.1.4 Vehicle routing problem Our problem can essentially be thought of as an intermodal instance of the Vehicle Routing Problem (VRP) [46, 93]. Since it was rst proposed by Dantzig and Ramser [38] where a simple matching-based heuristic for its solution was given, tons of papers have been devoted to the exact and approximate solution. Christodes, Mingozzi and Toth in 1981 published two papers on exact algorithms of VRP in Networks [29] and Mathematical Programming [30]. The rst algorithm is based on dynamic programming with state-space relaxation and the second algorithm proposed two mathematical formulations making use of q- paths and k-shortest spanning trees. Laporte, Desrochers and Nobert, Laporte et al. [64] proposed the rst cutting plane approach for a VRP based on the solution of linear relaxation of an integer model. In addition to that, Fukasawa and etc [44] formulated the problem as vehicle ow or commodity ow variables and are solved by branch-and-cut. Baldacci et al. [15] formulated the problem as a set partitioning problem to which some valid inequalities are added. For heuristic algorthms, there are several heuristics based on a variety of simple principles including savings, geographical and proximity, customer matchings, of which Clarke and Wright savings heuristic [31] has proven to perform well of its speed, simplicity and reasonably good accuracy. The development of research on metaheuristic algortihms also has eects on VRP. These include some tabu search-based approaches ([48]). The metaheuristic algorithms proven to work well always perform a wide and deep search of the solution including using adaptive large neighbourhood search (Pisinger and Ropke [78]), or combining genetic search with local search, as in the hybrid genetic algorithm proposed by Vidal et al. ([94]). There are also much reseach on variants of this problem such as the Capacitated VRP (CVRP), in which a homogeneous eet of vehicles is available and the only constraint is the vehicle capacity, 9 or the VRP with Time Windows (VRPTW), where customers may be served within a specied time interval and the schedule of the vehicle trips needs to be determined. Because we seek to optimally coordinate two classes of vehicles that have diametrically opposing strengths and weaknesses, our problem is particularly related to instances of VRP in which vehicle heterogeneity plays an important role, such as [47]. A closely related problem to ours is described in [67], which gives an integer programming formulation and a heuristic algorithm for solving a routing problem in which a truck can carry a eet of \foot couriers" on a single- or multi-route assignment, and the goal is to coordinate resources between the truck and the couriers eectively. The horse y problem also shares some aspects with the truck and trailer routing problem (TTRP). In this problem, a vehicle composed of a truck with a detachable trailer serves the demand of a set of customers reachable by truck and trailer or only by the truck without the trailer. Several heuristics based on the cluster rst-route second principle have been proposed to solve the TTRP, including Tabu search [86] and simulated annealing [89]. Furthermore, some approaches for exact solutions based on branch-and-cut [41] and branch-and-price [40] have been studied. The main dierence between the TTRP and the horse y problem is that both the truck and the drone can separately serve customer locations in horse y problem,while the trailer cannot serve customers without the truck in the TTRP. 1.1.5 Backbone network One of the basic phenomena that is of interest to us is the trade-o between eciency in trans- portation along a backbone network (in this case, the route of the truck) versus direct trips between locations (in this case, the direct trips taken by UAVs); this is arguably one of the funda- mental dichotomies in transportation and logistics[24, 25]. In this sense, our problem of interest is philosophically similar is [23], which asks whether small local retail stores are preferable to \big-box" retailers, with [97], which estimates the changes in net CO 2 emissions that result by 10 introducing grocery delivery services, and with [98], which computes the optimal layout of a set of facility locations that are themselves connected with a backbone network. An important problem in backbone network is hub location. Perhaps, Goldman [50] is the rst paper addressing the network hub location problem. However, O'Kelly [75] presented the rst recognized mathematical formulation for a hub location problem by studying airline passenger networks. His formulation is referred to the single allocation p-hub median problem. Given n demand nodes, ow between origin-destination pairs and the required number of hubs (p), the objective is to minimize the total transportation cost (time, distance, etc.) to serve the given set of ows. A generalization of this problem is multiple allocation where each demand center can receive and send ow through more than one hub. Campbell [27] was the rst to formulate the multiple allocation p-hub median problem as a linear integer programming problem. Our problem is closer to continuous or planar hub location problem where hubs can be located anywhere in a region other than demand nodes [26, 74]. Aykin and Brown [14] examine the problems involving interacting new facilities on a plane and on a sphere and an algorithm is given to solve the problems. Another problem on continuous approximation models where distribution of the demand rather than precise information on every demand point is given. These continuous approximations allow simple, yet plausible models to be developed, which often yield closed form solutions (e.g. minimal cost expressions [9]. The main goal of this approach is to obtain reasonable solutions with as little information as possible, and to gain a clear understanding of the trade- os. The monograph of Daganzo [35] is the rst presentation of the continuous approximation method in a more pedagogical manner. It addresses the dierent logistics costs and leads the reader through the development of a number of continuous approximation models. 11 1.2 Notational conventions Our notational conventions are as follows: we assume that there are n customer locations in the Euclidean plane to be visited with a truck and a UAV whose speeds are 0 and 1 respectively, with 0 < 1 and costs per unit distance arec 0 andc 1 withc 0 >c 1 . These customers are assumed to follow an absolutely continuous probability distributionf, which is dened on a compact planar regionR. We useL to denote a loop inR (representing the truck's tour), we use Loop(R) to denote the set of all loopsL inR whose length is well-dened, and we let d(x;L) denote the distance between point x to loopL; that is, d(x;L) = min x 0 2L kxx 0 k; wherekk is the usual Euclidean distance. The set of permutations off1;:::;ng will be written as S n , with a particular permutation written as 2 S n ; because these permutations always correspond to a tour that our truck takes, we will adopt the convention that (n + 1) =(1) for brevity. The (closed) -neighborhood of a (compact) setS in the plane will be written asN (S), which is to say, N (S) =fx2R 2 : min x 0 2S kxx 0 kg: We say that a functiong(x) satisesg(x)h(x) asx!1 if lim x!1 g(x)=h(x) = 1. Finally, we use the expression 1 x2S to denote the indicator function for membership in setS. 12 Chapter 2 Preliminaries In this chapter, we will introduce some core lemmas and theorems that will be used in some of our results. The following classical theorem, originally stated in [16] and further developed in [91, 92], relates the length of a TSP tour of a sequence of points with the distribution from which they were sampled: Theorem 1 (BHH Theorem). Suppose that X =fX 1 ;X 2 ;:::g is a sequence of random points i.i.d. according to an absolutely continuous probability density function f dened on a compact planar regionR. Then with probability one, the length TSP(X) of the optimal travelling salesman tour through X satises lim N!1 TSP(X) p N = ZZ R p f(x) dx where is a constant. It is additionally known that 0:6250 0:9204 and estimated that 0:7124; see [12, 16]. The following geometric result will be helpful to us later in relating the workload of the truck and the workload of the UAV; given a loopL (which represents the tour taken by the truck), it is useful to know how much area is within a given distance ofL: 13 Lemma 2. LetL be a loop in the plane with length` and let> 0. The area of an-neighborhood ofL is at most Area(N (L)) 8 > > > < > > > : 2` if ` 2 2 +` + ` 2 4 otherwise; which is tight whenL is a circle, wherebyN (L) is either an annulus (if `=(2)) or a disk (if >`=(2)). Proof. Assume without loss of generality thatL forms the boundary of a convex region and thatL is piecewise linear, i.e. polygonal (see the appendix for justications of both of these assumptions). The -neighborhood ofL has an \inner" portion R in and an \outer" portion R out as shown in Figure 2.1a, and we can also see from Figure 2.1a that the area of the \outer" portion R out is always exactly 2 +`. It is also obvious that the outer perimeter of R out is exactly ` + 2. The area of the \inner" portionR in is a little more complicated to bound; rst, for any 0 , we letL 0 denote the closed curve insideL consisting of points that are exactly 0 away from their nearest point inL (thus, our originalL would simply be written asL 0 under this notation). It is of course possible thatL 0 =; for suciently large 0 . Since R in is simply the union of all curvesL 0 over 0 2 [0;], the coarea formula [61] says that the area of R in is obtained by simply integrating the length ofL 0 from 0 = 0 to 0 =: Area(R in ) = Z 0 length(L 0) d 0 : Let R 0 out denote the \outer" portion of the 0 -neighborhood ofL 0, as shown in Figure 2.1b. We then haveR 0 out R in , which by convexity implies that the outer perimeter of R 0 out is less than or equal to the length ofL, i.e. `. However, by the same reasoning as our calculation of the outer 14 (a) (b) Figure 2.1: (a) shows the neighborhoods R in and R out ; note that R out is the union of the shaded rectangles and the hatched circular sectors, and therefore Area(R out ) = 2 +` and perimeter(R out ) =` + 2. (b) shows an inner loopL 0 and illustrates thatR 0 out R in for all 0 (which characterize R 0 out ). perimeter of R out , we also see that the outer perimeter of R 0 out is exactly length(L 0) + 2 0 , and therefore ` = length(L) perimeter(R 0 out ) = length(L 0) + 2 0 =) length(L 0) ` 2 0 provided thatL 0 exists, i.e. that length(L 0)> 0. Thus, we see that Area(R in ) = Z 0 length(L 0) d 0 Z 0 maxf` 2 0 ; 0g d 0 = 8 > > > < > > > : ` 2 if ` 2 ` 2 4 otherwise from which the desired result follows, since Area(N (L)) = Area(R in ) + Area(R out ). 15 Chapter 3 Continuous Horse y Problem 3.1 Introduction We begin by formally dening the continuous horse y routing problem in which we coordinate a truck and a single UAV: Denition 3. [Horse y routing with one truck and one UAV] Let p 1 ;:::;p n be a collection of points in the plane and let 0 ; 1 > 0 denote the speeds of a truck and a UAV respectively, with 0 < 1 . The optimal horse y tour of p 1 ;:::;p n is the solution to the optimization problem minimize x1;:::;xn;2Sn n X i=1 max 1 0 kx (i) x (i+1) k; 1 1 kx (i) p (i) k +kp (i) x (i+1) k (3.1) where S n 3 is the set of all permutations of the setf1;:::;ng, with the added convention that (n + 1) =(1) for notational convenience. Figure 3.1 shows an example of a solution to the above problem. Each variablex i corresponds to the \launch site" at which the truck releases the UAV to visit customer point p i ; the rst term in the maxf;g expression simply corresponds to the amount of time needed for the truck to move from one launch site to the next, whereas the second term corresponds to the amount of time needed for the UAV to leave its launch site, arrive at its customer point, and return to rendezvous 16 x 6 x 4 x 1 x 3 x 5 x 2 p 4 p 6 p 1 p 3 p 5 p 2 (a) 1 = 0 = 3=2 p 4 p 6 p 1 p 3 p 5 p 2 x 6 x 4 x 1 x 3 x 5 x 2 (b) 1 = 0 = 3 Figure 3.1: Solutions to an instance of the horse y routing problem, with points p 1 ;:::;p 6 as shown and with 1 = 0 = 3=2 and 1 = 0 = 3 respectively. Dening x 1 arbitrarily as a starting point, we have =f1; 3; 5; 2; 6; 4g in both routes. The fact that 1 = 0 = 3=2 in (a) is re ected by the observation that each pair of dashed line segments to and from a customer point p (i) (corresponding to a UAV's trip that starts with the truck, visits the customer, and returns to the truck) has length equal to 3=2 that of the thickened line segment fromx (i) tox (i+1) . Similarly, in (b), each pair of dashed line segments has length equal to 3 times that of the corresponding thickened line segment. with the truck at the next launch site. It is worth noting that, if one xes the permutation , the remaining optimization problem over variables x i is convex. This problem, although simple in its form, illustrates the core idea behind the delivery system with truck and UAVs described in chapter 1. All our problems later stated will be extension of the it in some certain aspect. For example, the problem considering battery life of UAVs is adding the limit on the distance that UAV can travel. The problem with multiple UAVs allow more than one UAV to travel along with UAVs. The real circumstances may be complicated, however, understanding the solution of continuous horse y problem will help understand the delivery system with truck and UAVs. Moreover, the techniques used in this problem can be easily extended to more complicated problems. Thus in this problem, we will explore this problem in many aspects. 17 3.2 Integer programming formulation Denition 3 of horse y problem is based on a permutation which is not a common setting when solving through computer programs. Formulating it by an integer programming instead will give us benets of solving it by commonly used methods like branching and bounding. Following the above notation, assumep 1 ;p 2 p n as a collection of points in the plane. Let c ij to be an integer variable of valuesf0; 1g with value 1 if the ith points UAV visits is p j and 0 otherwise where i = 1; 2; 3;n and j = 1; 2; 3n. Let t i to be the time that UAV and truck uses to visit the ith custom where i = 1; 2;n. We can get the following formulation Denition 4. [Integer programming formulation of horse y problem] The horse y problem with one truck and one drone can be formulated as follows: minimize n X i=1 t i (3.2) t i kx i x i+1 k 0 8i (3.3) t i M(1c ij ) + kx i p j k +kx i+1 p j k 1 8i;j (3.4) x n+1 =x 1 (3.5) X i c ij = 18j (3.6) X j c ij = 18i (3.7) c ij 2f0; 1g; x i 2R (3.8) (3.3) - (3.4) ensures that t i is the maximum time of truck and UAV travel. Notice that M is a very large number and can be made the double largest distance between two points in region 18 R divided by 1 . It makes sure that (3.4) doesn't cause any eect if c ij = 0. (3.6) - (3.7) ensures each point is visited once. One thing to mention here is that the right hand sides of (3.3) and (3.4) are both convex functions and thus the relaxation of the integer programming formulation is a convex optimization problem when the regionR can be written as some convex constraints. 3.3 Asymptotic analysis under uniform distribution In this section , we are going to derive upper and lower bounds for problem (3.1), under the assumption that the customers are distributed uniformly within a RegionR. Our work is mainly on truck with only one UAV in previous sections, but for upper and lower bounds, our result can be easily generalized to truck with multiple drones. Without loss of generality, we will assumeR having unit area in later sections. We begin with a key observation regarding the cost of a horse y tour in the asymptotic limit as the number of customer points p i becomes large: Theorem 5 (Limit of problem (3.1)). Let parameters 0 and 1 be xed and let T n denote the optimal cost of a horse y tour (i.e. problem (3.1)) of uniformly distributed points p 1 ;:::;p n in the unit square. Then we have lim T n p n =( 0 ; 1 ) with probability one as n!1. Moreover, we have ( 0 ; 1 ) = 1 1 ( 0 = 1 ; 1). Proof. For xed values of 0 and 1 , a horse y tour is a monotone subadditive Euclidean func- tional, and hence the desired result is an immediate consequence of well-established geometric probability theory; see [91]. The second observation is an obvious consequence of linear scal- ing. 19 Theorem 5 is appealing because it partially establishes the limiting behavior of a horse y tour, although it does not tell us anything about how the parameters 0 and 1 aect the optimal solution. Thus, this section is devoted to nding bounds for ( 0 ; 1 ). By linear scaling, as described above, we can assume without loss of generality that 1 = 1, and thereby concern ourselves only with the ratio r = 0 = 1 , which satises 0<r < 1. We therefore seek bounds on r :=(r; 1). An important technique on analyzing the lower and upper bound of the horse y routing problem is on partition. The following lemma, used in later analysis for both lower and upper bounds lays the foundation. Lemma 6. PutR with contour of length ` ctr in a griding system of identical squares with length (see Firgure 3.2 for an example), assume n = 1= 2 , thenR contains at least n 2 p 2` ctr = squares and is contained in a region of n + 2 p 2` ctr = squares. Proof. Consider two kinds of squares : the squares totally contained inR and the squares having some parts both inside and outside ofR. Assume the number of rst kind squares is n 1 and the second kind with number n 2 . It's obvious that n 1 nn 1 +n 2 . As all squares in the grid are disjoint and points in second kind squares has distance of at most p 2 toR. Thus from lemma 2, when we assume is small enough, 2 p 2` ctr n 2 2 orn 2 2 p 2` ctr =. Thus the two kinds of squares together will from a region with n 1 +n 2 n +n 2 n + 2 p 2` ctr = squares will contain R. As n 1 +n 2 n and n 2 2 p 2` ctr =, n 1 n 2 p 2` ctr =. Obviously, if customs are distributed uniformly in a larger region, it takes more time to visit times. That is to say, if regionsR 1 R 2 , the expected optimal time for horse y problem will have T (R 1 ) T (R 2 ). Lemma 6 tells us that when the squares in the grid is small enough, we can think of any regionR composed of 1= 2 squares. For the horse y problem, if we assume = 1= p n, we can think of any regionR as a region with n squares where n is the number of points to be visited. 20 R Figure 3.2: An example of griding an areaR, shaded area is the squares that contain both inside and outside area of R, here we have n 1 = 16;n 2 = 27 which are the number of squares totally contained inside R and number of squares of shaded area respectively. Let N to be the number of squares that contains some points infp 1 ;p 2 ;p n g. Its obvious thatNn. Actually,E(N) =E(N 1 )+E(N 2 )+E(N n ) whereN i = 1 if there are some points in ith square and 0 otherwise. The probability of a square to be empty is (1 1 n ) n . Thus when n!1, the approximation result E(N) (1 1 e )n holds almost surely where e is the Euler's number. This is an important property to reach a lower bound because obviously, if the tour can visit all points, each square containing some points will also be touched. As all squares are disjoint and having some areas thus can't be clustered together like points, it will lead us to some lower bounds. 3.3.1 Asymptotic lower bound In order to nd a lower bound of r , we will study a natural lower bound of problem (3.1): Lemma 7. Consider a set of input parameters 0 ; 1 ;fp i g to a horse y problem. The solution to the following innite-dimensional optimization problem is a lower bound of the horse y tour: minimize L max ( 1 0 length(L); 2 1 n X i=1 dist(p i ;L)) ) ; (3.9) 21 ` (a) ` (b) Figure 3.3: (a) is an example of horse y problem and (b) shows a possible solution of transformed problem from original horse y problem, the objective function of the problem is composed of the length of` and sum of distance from points to` This is also an example for horse y problem with objective function as cost. where the optimization variable consists of all loopsL with well-dened arc length and dist(p i ;L) is the distance from p i toL. Proof. AssumeL be the loop truck travels and its length be`, assume that drone travels distance ` d . The following inequalities about the time till completion T will hold T max( ` 0 ; 2` d 1 ) max( ` 0 ; 2 P n i=1 dist(pi;L) 1 ) The last inequality is because the fact that drone has to visit some p i starting from and coming back to some locations onL . Thus the 3.9 is a lower bound. Instead of nding n points and a permutation in (3.1), we are nding a loopL. It seems that we are upgrading a nite dimensional problem to an innite dimensional problem thus making it more complicated. But it's actually easier when coming to lower bound analysis. An example of the problem is shown in Figure 3.3. 22 The problem 3.9 is to make trade-o between truck and drone's travel distance. The longer L is, the shorter distance UAV needs to travel to reach L, the following lemmas will give a quantitative result on this trade-o. Lemma 8. LetL be a loop in the plane and letS be a region with area A. Then RR S d(x;L) dx satises ZZ S d(x;L) dx 8 > > > < > > > : 2A 3=2 3 p A` 2 + ` 3 12 2 if ` p A A 2 4` otherwise: (3.10) and equivalently the average distance is 1 A ZZ S d(x;L) dx 8 > > > < > > > : 2 p A 3 p ` 2 + ` 3 12A 2 if ` p A A 4` otherwise: Proof. We can assume without loss of generality that A = 1 because, if we apply the transforma- tion A7!cA and `7! p c` for any c> 0, then the right-hand side of the above is simply scaled by c 3=2 . Thus, the quantity of interest RR S d(x;L) dx is simply equal to the expected distance between a pointX uniformly sampled inS andL, i.e. Ed(X;L). Recall that for any non-negative random variable Z on the real line, we have E(Z) = Z 1 0 1F (z)dz; where F is the cumulative distribution function of Z (this is a simple consequence of integration by parts [52]). If we set Z =d(X;L), where X is uniformly sampled inL, then Lemma 2 tells us that F (z) = Area(N z (L)\S) Area(N z (L)) 8 > > > < > > > : 2z` if z ` 2 z 2 +z` + ` 2 4 otherwise: 23 Suppose that `> p as in the second case of the desired inequality. We then have ZZ S d(x;L) dx = Ed(X;L) = Z 1 0 1F (z) dz Z 1 2` 0 1F (z) dz Z 1 2` 0 1 2z` dz = 1 4` : On the other hand, if ` p as in the rst case of the desired inequality, we have `=(2) 1= p `=(2), whence ZZ S d(x;L) dx = Z 1 0 1F (z) dz Z ` 2 0 1F (z) dz + Z 1 p ` 2 ` 2 1F (z) dz Z ` 2 0 1 2z` dz + Z 1 p ` 2 ` 2 1 (z 2 +z` + ` 2 4 ) dz = 2 3 p ` 2 + ` 3 12 2 which completes the proof. Lemma 9. LetL be a loop and letG denote a grid consisting of N squares with side length . Then we have X g2G d(g;L) 8 > > > < > > > : N 3=2 4 p p 2N if ` p N N 2 2 4` p 2N otherwise: Proof. If we have N squares with area 2 , then we substitute the above lemma withS =G and A =N 2 . The average distance to a square diers from the minimum distance to that square by at most p 2, thus the sum of distances is at most X g2G d(g;L) 8 > > > < > > > : 2N 3=2 3 p `N 2 + ` 3 12 2 2 if ` p N N 2 2 4` otherwise: p 2N 24 it is easy to see that this is decreasing, so let's dramatically weaken the upper component by simply substituting ` = p N in the rst case: X g2G d(g;L) 8 > > > < > > > : N 3=2 4 p p 2N if ` p N N 2 2 4` p 2N otherwise: The remainder of the proof is simple: the above expression is obviously continuous and de- creasing in `, thus we want to minimize maxf 1 r `;f(`)g, which occurs at the intersection of the two components, which gives an objective value of opt = 8 > > > < > > > : N 3=2 4 p p 2N if N 4 p 2N r 4 p N 2 p 2 + 1=r p 2 otherwise. We now sample n points, set = 1= p n (or a ceiling or whatever), and the law of large numbers says that the number of empty squares N satises N=n! := (1 1=e) with probability one. Because of this, we see that N 4 p 2N r!1 with probability one and so we can disregard the rst component of the piecewise expression. Substituting N 7! n and 7! 1= p n in the second component gives a lower bound of 2 r 2r + 1 r p 2 ! p n 0:32 r 2r + 1 r 0:45 ! p n: This lemma gives lower bound r > 2 p 2r + 1 p 2r p n= p r. However, it works well when r 0 but much worse even than TSP tour when r 1, it's because in the process of creating a lower bound tour, UAV's travel distance is underestimated too much when 0 1 . We will next claim another lower bound which is simple but works well in the case that 0 1 . 25 Theorem 10 (Lower bound 2). r tsp Proof. Assume ` to be the length of truck's drone and ` d to be the distance that drone travels. Time till completion is max( ` r ;` d ). The tour UAV travels form a tour visiting all points, thus is a candidate for TSP, this is also true for truck. Thus both ` d and ` are less than tsp p n which means r tsp . Remark. The lower bound only works for the case where UAV leaves immediately from truck after it meets with truck at some location. Thus the lower bound won't work for cases discussed later where cost or battery life is considered. However, we always have the fact that the tours of truck and drone combined together a solution for TSP. It means when solving minimize max( ` r ;` d ) subject to ` +` d tsp p nS It gives a lower bound of r tsp 1+r which is slightly worse. For simplicity and a form that is in accordance with upper bound later, we will write the lower bound as r = (r)= p r. We dene it this way because the rst lower bound shows us that the speed r blows up is 1= p r. Combine the the two lower bounds together, we will get (r) max tsp p r; 2 ( p 2r + 1 p 2r) with (r) 0:2103 (See gure 3.4). 3.3.2 Asymptotic upper bound Now we are deducing an upper bound for problem (3.1) by giving an algorithm and analyze the performance of it. We rst partition the regionR by K small squares (notice that it's dierent from our previous partition where we assume the number of small squares is equal to number of customs n). Truck will take a tour of visiting center of each square once, drone will serve all customs in the square by traveling back and forth from customs to the truck. Truck will wait until all customes were visited. An example is shown in the gure 3.5. Theorem 11 (Upper Bound). r 2 p 2c p r 26 Figure 3.4: The plot of two lower bounds for with respect tor, the second lower bound dominates whenr is close to 0 then rst lower bound dominates. r which dened as the maximum of the two bounds reaches its minimum when the two bounds become the same at r = 0:0877 with value 0:2103. When r is 0, (r) has value 0:3161 and when r = 1, (r) = tsp . Figure 3.5: An example of the approximation algorithm, customs are distributed in a unit square and the region is divided into 9 regions with each region a square. The bold tour is what truck travels. It touches each center of square. The line froms points to the tour are the tour drone travels. Truck will wait at the center of square until all customs were visited. 27 Proof. Following the scheme above, the distance truck travels is always p K (some constant dif- ference existis depending on the tour we used to visit all centers of each small squares, but as K goes to innity when n is large, this number is negligible). For visiting each custom, the average distance for drone will be 2c= p K where c 0:383, thus for time of completion T = p K r + 2cn q 1 K When K = 2cnr, we will have the upper bound for time of 2 p 2cn p r and thus r 2 p 2c p r . Before heading to result of general density function, let's rst wrap up the results we have now. First, from theorem 5, lim Tn p n = r exists, and combining the results from lower and upper bound, we can actually write it as r = (r)= p r, with 0:2103 (r) 1:7504. Then we can state a BHH similar result under uniform distribution. Theorem 12 (BHH for horse y problem under uniform distribution). Suppose thatX =fX 1 ;X 2 ; :::X N g is a sequence of random points i.i.d. uniformly distributed within a compact planar region R 2 with area A. Then with probability one, the optimal time T N of horse y problem with one UAV through X satises lim N!1 T N p N = (r) p 0 1 p A where (r) is only related with r = 0 = 1 and 0:2103(r) 1:7504. 0 , 1 are speed of truck and UAV respectively. 3.4 Asymptotic analysis under general density function Now, we are going to show the asymptotic results for the general absolutely continuous density function f with compact support. The proof is based on the fact in probability theory that any absolutely continuous density function can be approximated by step functions. It was applied in proof of BHH theorem in [92]. Before heading to the proof, we will rst state variant of an 28 important lemma used in proving BHH theorem. It is going to play a role later in the proof of asymptotic results for horse y problem. Lemma 13. Without loss of generality, consider the regionR with unit area that can be partitioned into m 2 squaresfQ i g with each square having each side of length 1 m . Assume X 1 ;X 2 ;X n are points inR. The the following inequality holds length(TSP (X 1 ;X 2 ;X n )) m 2 X 1 length(TSP (fX j :X j 2Q i g)) + p 5m (3.11) m 2 X 1 length(TSP (fX j :X j 2Q i g))length(TSP (X 1 ;X 2 ;X n )) + 10 p 2m + 8m (3.12) WhereTSP (X 1 ;X 2 ;X n ) represents the TSP tour visitingX 1 ;X 2 ;X n andlength(TSP (X 1 ; X 2 ;X n )) the length of it. Proof. The proof of (3.11) is obvious. Each TSP (fX j :X j 2Q i g connects the points inside Q i . Connecting all the Q i together will give us a tour visiting all points X 1 ;X 2 ;X n which is a feasible solution to TSP problems which is larger than the optimal solution TSP (X 1 ;X 2 ;X n ). Connecting adjacent Q i will cost at most p 5 m . Thus (3.11) holds. Proof of (3.12) needs some extra work. Consider the TSP tour inside a specicQ i orTSP (X 1 ;X 2 ; X n )\Q i . To prove (3.12), we need to make a tour visiting all points inside each Q i out of the TSP tour. Obviously, if the truck enters and exits Q i just once, keeping the original tour without changes will nish the work. However truck may enter and exist multiple times as shown in gure 3.6a. Notate the boundaries of Q i to be @Q i . It consists of four sides of square Q i and use @Q ij to represent them where j = 1; 2; 3; 4. Considering @Q ij \TSP (X 1 ;X 2 ;X n ) which are the entering and exiting points of truck . If the number of points is even, we are good now and can go to next step. Otherwise, take any point out to make the number even. For convenience, we label points of TSP (X 1 ;X 2 ;X n )\Q i in clockwise order for each j. We then need to do two steps of preprocessing. 29 First, connect the last point of @Q ij \TSP (X 1 ;X 2 ;X n ) with the rst point of @Q i;j+1 \ TSP (X 1 ;X 2 ;X n ) for each j = 1; 2; 3; 4 (j + 1 = 1 if j = 4). Second, assume @Q ij \ TSP (X 1 ;X 2 ;X n ) consists of 2d points (one point is excluded if there are odd number of points). We double all the edges connecting the (2i 1)th and (2i)th points where i = 1; 2;d. We are doing these steps because when assuming all X i s are not identical and not on @Q i , there is only one edge passing points in @Q i;j+1 \TSP (X 1 ;X 2 ;X n ) in TSP (X 1 ;X 2 ;X n ). And thus after the above preprocessing, the degrees of those points are all even which guarantees a tour using those edges exactly once. We have some extra points (at most of number 4) excluded before. Connecting all those points together and use two edges to connect one of the points to any point in@Q ij \TSP (X 1 ;X 2 ;X n ), this will make a tour of pointsfX j : X j 2 Q i g which will be larger than TSP (fX j : X j 2 Q i g. An example is shown in 3.6b. At step one, we add four lines conencting points in each side of square which costs at most 4 p 2 m . The second step costs at most double the length of the sides which is 8m in total. The last step adding extra points in will add at most 6 lines with costs at most 6 p 2 m . Thus inequality (3.12) holds. The above Lemma shows that by partition a region into small squares and nd a TSP tour for points inside of each square, we get both a lower and upper bound for optimal solution of TSP tour. Based on this lemma and its proof, we are giving a similar result for horse y problem. Assume HF (X 1 ;X 2 ;X n ) to be optimal tour for horse y problem with only one UAV for visiting customs at X 1 ;X 2 ;X n And time(HF (X 1 ;X 2 ;X n )) the time of completion for the tour. The following lemma will hold 30 (a) (b) Figure 3.6: (a) An example of possible TSP (X 1 ;X 2 ;X n )\Q i which is the part of a TSP tour inside a small square. (b) nal graph after we taking steps of connecting each side of the square, doubling every other edge connecting points on the sides and adding extra points. The bold line means the line is used twice in the nal tour. We can see that each point on the sides of square has even degree and thus can form a tour that will take each edge just once. Lemma 14. For a regionR with unit area that can be partitioned intom 2 squaresfQ i g with each square having side length 1=m. AssumeX 1 ;X 2 ;X n are points inR. The following inequalities hold time(HF (X 1 ;X 2 ;X n )) m 2 X 1 time(HF (fX j :X j 2Q i g)) + p 5m 0 (3.13) m 2 X 1 time(HF (fX j :X j 2Q i g)) 2time(HF (X 1 ;X 2 ;X n )) + 10 p 2m + 8m 0 (3.14) Proof. (3.13) is easy to proof. When we nish serving customs in each square, truck and UAV travel together to the next square. This will make a candidate to the original horse y problem. For (3.14), we want to create horse y tour for eachfX j :X j 2Q i g out ofHF (X 1 ;X 2 ;X n ). First we recreate the tour according to methods shown in gure 3.7 which will guarantee points inside the square be visited and at the end, truck will hold its UAV at the loading location for next delivery if it's inside square or at some points on the square sides. Then follow the same proof in 3.12, we will be able to get a tour connecting all tours with dierence 10 p 2m+8m 0 . The only thing left to be proved is the new tour we created in 3.7 is within two multiples of the time in original horse y tour. In order to do that, we just need to calculate how much time are 31 used in each square for tour visiting a specic custom, sayp i . As gure 3.7c takes no time and we can ignore it, the tour visiting p i can be calculated in new tour with following combinations: (a) + (f), (b) + (b) + (f), (b) + (e), (d). It's obvious that in each case, the time added together in all cases will not exceed two multiples of the time in the original tour. This completes our proof. If the density function has the format of step function f(x) = P m 2 k=1 a k I Q k (x). Since inside each Q i , customs are uniformly distributed and thus from Theorem 5, we have time(HF (fX j : X j 2 Q i g) 1 m p a k m 2 n= p 0 1 . Combining lemma 14 will give us asymptotic time of T n = p n r RR R 2 p f(x)dx 2T n = p n. This step is exactly the same as the proof of BHH theorem, details can be found in [92]. It means that the results we got in section 3.3 can be extended to general density function f in a weaker version as follows Theorem 15 (BHH for horse y problem). Suppose that X =fX 1 ;X 2 ;:::g is a sequence of random points i.i.d. according to an absolutely continuous probability density function f dened on a compact planar regionR 2 . Then with probability one, the optimal time T N of horse y problem with one UAV through X satises (r) 2 p 0 1 ZZ R 2 p f(x) dx lim inf N!1 T N p N lim sup N!1 T N p N (r) p 0 1 ZZ R 2 p f(x) dx where (r) is dened the same as that in theorem 12. 3.5 Cost as objective function In previous sections, we all focus on time of completion as objective function. However, sometimes cost may be also one of the considerations. Especially in our later chapter where we consider both cost and time of completion as our objectives. So it's good to mention the problem with cost as 32 p i l1 l2 (a) p i l1 l2 (b) p i l1 l2 (c) p1 l1 l2 (d) p i l1 l2 Meeting point (e) l1 l2 p i Truck waits here (f) Figure 3.7: The six pictures above show how we create a truck tour out of the original horse y tour in the six dierent cases. We use dashed line to represent the tour drone travels, black lines for truck tours in the horse y tour and bold line for the new truck tour. In (a) and (d), truck travels along the same tour, in (b), truck only takes also the same tour as the original horse y tour but only the part within the square. For (c), we can just ignore this case. For (e), truck takes the tour within square, and when it hits any side of the square, it travels toward the drone along square sides until it gets the drone. For (f), truck stays at the intersection point of square and the line chosen between p i l 1 and p i l 2 with shorter length, drone will leave to visit p i and then comes back. Those tours will guarantee truck at the end reaching the loading location l 2 with drone loaded for the next delivery and also served custom p i if it's inside square. 33 (a) Time (b) Cost Figure 3.8: Seven points are generated inside unit square, (a) is the solution for horse y with time as objective function , (b) is the horse y problem with cost as objective function. The dierence is obvious from the intrinsic nature of the problems. For time as objective function, UAV will travel from one launching site to custom and then back to the next launching site. When cost is objective function, UAV will travel back to the same launching site as it's the shorter from triangular inequality. objective function. Most results are not very dierent from what we had before. Thus we mention them here instead of giving an additional chapter. Here we use c 0 to represent the unit cost of truck and c 1 to represent the unit cost of UAV. Notice here that c 0 > c 1 as truck will generally have higher cost than UAV. An illustration of dierence of solutions between problems with cost and time of completion as objective function can be seen in Figure 3.8. The problem can be derived easily following the integer programming formulation of horse y problem with time of completion as objective function. The only dierences are the constraints parts: minimize n X i=1 cost i (3.15) cost i M(1c ij ) +c 0 kx i x i+1 k + 2c 1 kx i p j k8i;j (3.16) x n+1 =x 1 (3.17) 34 X i c ij = 18j (3.18) X j c ij = 18i (3.19) c ij 2f0; 1g; x i 2R (3.20) The problem can be solved directly when the size is small. And to solve the problem of large size with heuristic algorithms described later, we need the following formation similar with denition 3: Denition 16. [Horse y routing with one truck and one UAV] Let p 1 ;:::;p n be a collection of points in the plane and let c 0 ;c 1 > 0 denote the speeds of a truck and a UAV respectively, with c 0 c 1 . The optimal horse y tour for cost of p 1 ;:::;p n is the solution to the optimization problem minimize x1;:::;xn;2Sn n X i=1 c 0 kx (i) x (i+1) k + 2c 1 kx (i) p (i) k (3.21) where S n 3 is the set of all permutations of the setf1;:::;ng, with the added convention that (n + 1) =(1) for notational convenience. We will claim the some theorems directly with brief proof if there are dierences from the problem minimizing waiting time. First, under uniform distribution, recall that in the proof of our previous problem with waiting time as objective function, we nd lower bound by converting the lower bound to the horse y problem with cost as objective function, thus the lower bound also holds. For the upper bound, the solution we created is also a feasible solution for this problem, thus the following lemma holds automatically. Theorem 17 (BHH for horse y problem under uniform distribution (cost version)). Suppose that X =fX 1 ;X 2 ;:::X N g is a sequence of random points i.i.d. uniformly distributed within a 35 compact planar regionR 2 with areaA. Then with probability one, the optimal cost C N of horse y problem with one UAV through X satises lim N!1 C N p N =(r) p c 0 c 1 p A where(r) is only related withr =c 1 =c 0 and 0:2103(r) 1:7504. c 0 ,c 1 are unit cost of truck and UAV respectively. Remarks. There are some remarks we need to mention here. First, although we used the same notation (r) both in theorem 12 and theorem 17. Actually, they are dierent. Later, we will use cost and time to distinguish them when necessary. Second, the ratio r in this problem is UAV unit cost/truck unit cost while in the problem with waiting time as objective it's truck speed/UAV speed. It's obvious because we need to make sure thatr is between 0 and 1. Another thing to mention is that from the proof of theorem 12, we have the time (r) cost (r) 2 time (r) for any specic 0r 1. We will develop more in our next section on computational results. For problem under general distribution, the result turns out to be cleaner due to the super- additivity property of the solution. As the proof has some dierences, we spend some time showing the lemma and theorem related with it. Lemma 18. For a regionR with unit area that can be partitioned intom 2 squaresfQ i g with each square having side length 1=m. AssumeX 1 ;X 2 ;X n are points inR. The following inequalities hold cost(HF (X 1 ;X 2 ;X n )) m 2 X 1 cost(HF (fX j :X j 2Q i g)) + p 5m 0 (3.22) m 2 X 1 cost(HF (fX j :X j 2Q i g))cost(HF (X 1 ;X 2 ;X n )) + 10 p 2m + 8m 0 (3.23) Proof. The only dierence here from lemma 14 is that we don't have the coecient 2 in the rst term of equation 3.23. The reason is that when considering the waiting time, truck and UAV 36 Figure 3.9: Left hand side is a feasible solution for the horse y problem with cost as objective, six points need to be visited. Bold lines are the tour that truck travels and dotted lines as tour of UAV, same as before. The right hand side is an example of how we partition the problem into smaller problems in the four sub-squares. What we do is that we make truck's traveling path a tour rst, if UAV's launching sites are not in the current sub-square, we will set launching sites to be the intersection of UAV's path with the side of the sub-square, as the bottom left sub-square in the graph above. travels together and when we do partition, the split may be uneven and cause some additional time as the time of truck and UAV have to be counted separately. However, when we are handling with the cost, the costs of truck and UAV should be considered separately anyway and thus partition doesn't incur any additional cost, it's illustrated in gure 3.9. Combining the techniques we used in proof of lemma 14, it's easy to see that the lemma holds. After getting the above lemma, following the exact the same discussion between lemma 14 and theorem 15, we will get the BHH theorem we mentioned for cost objective: Theorem 19 (BHH for horse y problem (cost version)). Suppose that X =fX 1 ;X 2 ;:::g is a sequence of random points i.i.d. according to an absolutely continuous probability density function f dened on a compact planar regionR 2 . Then with probability one, the optimal time C N of horse y problem with one UAV through X satises lim sup N!1 C N p N =(r) p c 0 c 1 ZZ R 2 p f(x) dx 37 where (r) is dened the same as that in theorem 17. 3.6 Computational results Computational results of continuous problem will be considered in this section. As the problem can be solved by directly from their integer programming formulations, we tried problems of small sizes and one result is shown in gure 6.1 which includes the cases for both waiting time and cost as objectives. Use Gurobi Solver to solve the problems of size 10 takes hours for horse y and horse y-BL. In order to compare the solutions of optimal solution for TSP with that for horse y problems for larger size, we will use the following heuristic algorithm: 1. run Lin-Kernighan heuristic [66, 53, 54] to get a near-optimal solution of TSP tour and use this order as an input for horse y problem. 2. solve it through a corresponding convex optimization problem. Figure 3.11 is an example for solving 100 points. This is the basic structure of the heuristics. But there will be some variations for horse y problems under dierent scenarios. We will explain in details when the corresponding problem appears. There are two benets of this algorithm, rst the Lin-Kernighan heuristic is fast enough for required experiments, and an order making a good solution for TSP will also make the that of horse y problem not bad. In later literature, this heuristic will used repeatedly for numerical experiments. 3.6.1 (r) for time of completion as objective Theorem 15 shows that the time of completion for horse y problem is around (r) p n p 01 asymptot- ically under uniform distribution. (r) varies according to dierent ratio between 0 and 1 , 38 (a) TSP (b) Horse y (c) Horse y-cost Figure 3.10: Eight points are uniformly generated from unit square, (a) is the tour of TSP which gives 3:228 time units, (b) is from solving the integer programming formulations before respectively for horse y with waiting time and cost as objective function under this situation. Here, truck is set to be half the speed of drone and unit cost is double of that of UAV. Horse y problem gives optimal time 1:873 and optimal cost 3:017. It satises the relationship we claimed before that optimal time will be less than optimal cost but larger than half of optimal cost given the ratios r are the same (here 0:5 for both cases). Although both numbers are less than that of TSP problem. 39 (a) TSP (b) Horse y (c) Horse y-cost Figure 3.11: 100 points are generated randomly from a unit square, (a) is the near optimal solution of TSP by running Lin-Kernighan heuristic. The total length is 7:9346. The ratio 7:9346= p 100 = 0:79346. (b) is the near optimal solution of horse y problem with the same sample points. It is reached by taking the solution of (a) as default order of points to be visited and solve the optimization problem in denition 3 with this order per-dened by this order. (c) is the solution for problem with cost as objective function by using the same method for getting (b). 40 although it can be bounded from above and below some constants. A natural question for com- putational experiments is on good estimates for (r) as the lower and upper bound given above still have a large gap and behavior of (r) when 0 ; 1 vary. For simulation, n samples are rst drawn uniformly in a unit square, then Lin-Kernighan heuristic will be run for a near optimal solution for TSP, it will be used to solve equation (3.1) which is convex and can be solved fast, this process will be repeated three times and the average will be taken to undermine the randomness. The cases for n = 500 and 100 are shown in gure 3.12. When r is very small even compared to 1 n , the drone only delivery time n 1 is better than horse y system of (r) p n p 01 , under this situation, the y value in the simulation is almost p n p r which denitely decreases as ratio decreases. A guess is (r) doesn't have the initial increasing part and is actually a inverted bell curve. To test our guess, we do the simulation on xed ratio r = 0:01 but changing sample size and result is shown in gure 3.13. (r) becomes stable after n becomes larger. The increasing part of the curve when ratio goes close to 1 is because of the usability of drone tends not to be constrained as its speed is almost the same as the truck. TSP tour will be the optimal solution when r is exactly 1. An important thing is the minimum value that can be served as a lower bound. Our simulation gives the ratio 0:911 whenr is around 0:65, with the fact TSP 0:7124, the value is about 0:6490. This value is with only small variance in our simulations and can be seen from gure3.12. Last but not least, the maximum value of (r) should be reached at a point close to 0, if our assumption of monotonicity is correct, this value is lim r!0 (r). For dierent sample sizes, golden-section search is implemented on r to nd the peak value of (r)= tsp , some results are shown in the table below: One thing need to mention here is that our simulation shows large variance for those values, but it's obvious that when n is less than 1500, the peak value of (r)= tsp is increasing. However after that, we are not sure whether it's still increasing or vibrating because of large variance. We adopt 1:3737 or = 1:3737 tsp 0:9786 as an estimate to it. 41 n r (r)= tsp 500 0:1627 1:2185 1000 0:1485 1:2640 1500 01137 1:3031 2000 0:1346 1:3327 2500 0:0900 1:4040 3000 0:0873 1:4136 3500 0:0907 1:3870 4000 0:1063 1:3949 4500 0:0814 1:3607 5000 0:0770 1:3737 Table 3.1: (r) values with dierent n (a) n=500 (b) n=1500 Figure 3.12: (a) and (b) are both the simulation results for nding (r) where n = 500; 1500 respectively, x-axis is p r and y-axis is the ratio between horse y and tsp solution multiplied by q 1 0 which is about (r)= tsp , both the pictures have the similar shape where y value goes to a peak, then decreases to a bottom and nally goes back to 1, the dierence is that when n = 1500, the shape is much sharper and reaches a higher peak at a lower r. 42 Figure 3.13: When x the ratior = 0:01 or p 0 = 1 = 0:1, the y-value which is about(r)= tsp increases fast when the sample size is less than 1500. Increasing the sample size even more, the value vibrates between 1:3 and 1:4. It illustrate the fact that the initial increase in gure 3.12, specically when p 0 = 1 = 0:1, is caused by the comparability of 0 = 1 to n and may disappear when sample size becomes larger. With the above discussions and simulation, we give a bold guess on the following three prop- erties of function (r) without rigorous proofs: 1. (r) is an inverted bell curve 2. (r) reaches minimum value at around 0:65 with (0:65) 0:6490 3. (r) has an upper bound at lim r!0 (r) 0:9786 Remark. Although the value of(r) is still not known accurately near 0, the value of is around 1:15 tsp in the extreme case that the speed of UAV is regarded as ten times speed of truck where p r 0:316. In practice, r should be much smaller than that thus the value of (r) is quite at. Using the coecient in our simulations, for example(r) = 1:05 tsp 0:7455, for time estimation will not be bad. 3.6.2 (r) for total cost as objective In this section, we are doing the same thing for the problem with cost as objective. We can then compare the two (r) so that we can know how much improvement can be made by using 43 (a) n=500 (b) n=1500 Figure 3.14: (a) and (b) are both the simulation results for nding (r) where n = 500; 1500 respectively, x-axis is p r and y-axis is the ratio between horse y and tsp solution multiplied by q c0 c1 which is about (r)= tsp , both the pictures have the similar shape where y value goes to a peak, then decreases 1, the dierence is that when n = 1500, the shape is much sharper and reaches a higher peak at a lower r. UAVs for delivery. In addition to that, although we are claiming the problem as horse y problem with total cost as objective, there is another way to view this problem. It can be regarded as the problem requiring truck to stop when UAV is in use. This also happens a lot especially we need manually control UAVs. Comparing these two functions can give us a sense of how much improvement can we make to automate UAVs. Simulation results are shown in gure 3.14. One similar part with the previous problem with time as objective is the initial increasing part. And as discussed before, that part should not exist and is caused by the fact that r is large compared to n. But it's dierent as the value of (r) decreases to 1 directly without the inverted bell curve. And the variation is much larger from 0:72 to more than 1:5. Some interesting results will thus come out from simulation here by comparing the two (r). For convenience, we will use time and cost here. First, when ratio r is close to 1, the two values will be both equal to 1. However, when r is not signicantly large, let's say r = 0:5, using UAV doesn't give much improvements. Because (r)= TSP 1:35 and p 1=r 1:4. Less than 5% improvement is achieved. But for time of completion, more than 40% improvement is 44 gained. When r is very small, say r = 0:01. time = TSP 1:4 and cost = TSP 2:2 which also gives around 50% improvements. Combining these together, we are concluding that if we are automating UAV in the system, it will give us around 50% savings in time. Cost will not be aected if we are not considering cost of labor work. 3.6.3 Experiments using road network distances In our third experiment, we use the results of the preceding section to predict the changes in service time when UAVs are introduced to a truck that visits a collection of n destination points in a map of Pasadena, California, with n2f25; 50; 100g. Thesen destination points are sampled uniformly from the centers of all 1734 US census blocks that belong to Pasadena, as shown in Figure 3.15b, which were obtained from [19]. Next, we compute a TSP tour of those n points with respect to the road network distance, as shown in Figure3.15c; the length and duration of this TSP tour is estimated using the Google Maps Directions API [2], the Google Distance Matrix API [5], and the Concorde TSP Solver [1]. Finally, we apply the same two-stage heuristic from the previous section to compute a horse y tour of those same points, as shown in Figure 3.15d The "adjusted" truck speed 0 is generally around 20 km hr in these experiments, and the UAV speed 1 satises 1 2f30 km hr ; 40 km hr ; 50 km hr ; 60 km hr g. Based on our preceding analysis, we expect that the two service times should satisfy Service time with UAVs = Service time without UAVs p k 1 = 0 (3.24) where there is k = 1 UAV and so we use = 1:00 (as determined from the previous section). Figure 3.16a shows these "predicted" service times compared with the actual service times for 10 individual experiments for n2f25; 50; 100g and 1 2f30 km hr ; 40 km hr ; 50 km hr ; 60 km hr g. The "quality" of our approximation, so to speak, is s therefore obtained by dividing the actual time by the predicted times (which we hope is close to 1); 45 (a) (b) (c) (d) Figure 3.15: Figure (a) shows a map of Pasadena, California, and (b) shows the initial dataset consisting of the centers of all 1734 US census blocks located in Pasadena, California. Figure (c) shows the TSP tour with respect to the road network of n = 25 points sampled from those centers, and (d) shows the horse y tour computed using the two-stage heuristic. 46 1 2 3 4 5 6 7 8 9 10 0 0.5 1 1.5 2 2.5 1 2 3 4 5 6 7 8 9 10 0 1 2 3 1 2 3 4 5 6 7 8 9 10 0 2 4 6 1 2 3 4 5 6 7 8 9 10 0 0.5 1 1.5 2 2.5 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 1 2 3 4 5 6 7 8 9 10 0 2 4 6 1 2 3 4 5 6 7 8 9 10 0 1 2 3 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 1 2 3 4 5 6 7 8 9 10 0 2 4 6 1 2 3 4 5 6 7 8 9 10 0 0.5 1 1.5 2 2.5 1 2 3 4 5 6 7 8 9 10 0 1 2 3 4 1 2 3 4 5 6 7 8 9 10 0 2 4 6 w/o UAV prediction actual trial # time (hrs) (a) Figure 3.16: Each of the 12 charts above corresponds to an experiment performed for a particular number of customersn and UAV speed 1 . Each chart contains 10 triplets, corresponding to 10 individual experiments. The three bars (black, gray, and white) in each triplet are interpreted as follows: the black bar simply indicates the amount of time, in hours, that a single truck takes to visit the n points. The gray bar represents the predicted amount of time that the truck and the UAV will take in a horse y tour, as determined by the right-hand side of 3.24. The white bar indicates the actual amount of time that the truck and the UAV take in performing their horse y tour. 47 Chapter 4 Horse y Problem with Multiple UAVs In our previous chapter, the delivery system contains only one truck and one UAV. However, there is no technical barrier that prevents us from using two UAVs in our system if we can implement the system with one UAV. And recently, Mercedes made a new drone truck with two drones[7]. Truck with more UAVs will denitely come out in near future to further improve the eciency of delivery system for solving the last-mile problem. In this section, we are handling with this system, following the approach to solve horse y rout- ing problem with one UAV, we will rst give theoretical results which are asymptotic like before, then we will introduce heuristic algorithms for solving this problem which will be a generalization of the previous problem with one UAV. Finally, some simulation results will be given. In general, most of the results we got before can't be extended directly without some mod- ications. However, basic ideas can be applied here so getting to the results for multiple UAVs don't take a lot of extra work. We will emphasize on the dierences in this chapter when coming to proofs. In this chapter, we assume our system now has k drones and as before, each drone is with capacity of one package but can leave truck to delver packages. And each UAV can travel inde- pendently with each other and truck. 48 4.1 Theoretical results The rst result for existence of limit comes directly with the fact that the sub-additivity property of the problem with multiple drones, following the notation we used before, the theorem follows: Theorem 20 (Existence of limit under uniform distribution). Time of completionT n is monotone and sub-additive, thus the following result holds from [91] lim T n p n =(r;k) The next step is to bound(r;k), remember that in our previous problem with only one UAV, we have two lower bounds and thus for the multiple UAVs problem, we are also giving two lower bounds. Lemma 21 (Lower Bound 1). (r;k) 2k q 2r+k r p 2 (0:32 q 2r+k r 0:45)=k Proof. The proof follows the proof of lemma 8 and the discussion after it. In our previous proof, after getting the lower bound in lemma 8, we minimize max(`=r;f(`)) where the rst term is the time for truck and second as the time of UAVs. Here as we have k UAVs and the best we can achieve is max(`=r;f(`)=k) if the time each UAV is coordinated well. We can imagine that r is substituted with r=k this time and it's not to reach the lower bound: 2k r 2r +k r p 2 ! p n (0:32 r 2r +k r 0:45) p n=k and this completes the proof. The second lower bound comes by assuming that we are visiting all the points with tsp tour and all UAVs will be able to nish 1=k parts at the same time which means that the total time will be tsp p n k if we assume the UAV has speed 1 as before. Thus we will have the following lower bound: 49 k lower bound of (r;k) 1 0:2111 2 0:2111 3 0:2111 0:2111 12 0:2124 20 0:2315 1 0:32 Table 4.1: Lower bound of (r;k) for dierent k Lemma 22 (Lower Bound 2). (r;k) tsp k Before heading to upper bound, we will see how the lower bound performs by combining the two bounds together. Notice that as we will show (r;k) = p rk so we will compare the value of (r;k) p rk. Some results are shown in table 4.1. When k is very large or r is very small, the rst lower bound will dominate. Our results show that when k 12, we just need the rst lower bound. And before that, as the intersection part always have the same y value 0:2111, we will thus have the xed lower bound for (r;k). Then we will claim the upper bound: Lemma 23 (Upper Bound). (r;k) 2 p p 2(1+ p 2) p kr Proof. Doing the same for the partition algorithm, but now, if we assume the ith square has n i customs, an upper bound of the time drones travel inside the ith square will bed ni k e p 2S 1 p K wheredne the smallest integer that is larger or equal to n. Thus an upper bound for the time till completion T in this scheme will be: T = p KS 0 + K X i=1 d n i k e p 2S 1 p K p KS 0 + K X i=1 ( n i k + 1) p 2S 1 p K = p KS 0 + p 2KS 1 + n p 2S k 1 p K 50 (1 + p 2) p KS 0 + n p 2S k 1 p K Now when K = p 2n0 (1+ p 2)k1 , we will have the upper bound 2 p p 2(1+ p 2) p S p k01 . Combine all the discussions above will give the following theorem which is a weaker version of theorem 15 Theorem 24. Suppose that X =fX 1 ;X 2 ;:::g is a sequence of random points i.i.d. according to uniform distribution dened on a compact planar regionR 2 . Then with probability one, the optimal time T N of horse y problem with k UAV through X satises lim N!1 T N p N = (r;k) p k 0 1 p A when N is large enough where 0:2111(r;k) 2 q p 2(1 + p 2) and A is the area ofR 2 . This theorem is dicult to be extended to the case of arbitrary function because lemma 14 is not easy to use here. Thus we are approaching the result under general density function in a dierent way. The basic idea is to set up a relationship with the problem with cost as objectives. Theorem 25 (BHH for horse y problem). Suppose that X =fX 1 ;X 2 ;:::g is a sequence of random points i.i.d. according to an absolutely continuous probability density function f dened on a compact planar regionR 2 . Then with probability one, the optimal time T N of horse y problem with k UAVs through X satises (r) (k + 1) p 0 1 ZZ R 2 p f(x) dx lim inf N!1 T N p N lim sup N!1 T N p N (r;k) p k 0 1 ZZ R 2 p f(x) dx where (r) is dened the same as that in theorem 12 and (r;k) dened in theorem 24. Proof. The second inequality comes from the sub additivity of the problem which is discussed before when we are proving theorem 12. For the proof of the rst inequality, let's assume ` to be 51 the length of tourL that truck travels and ` i to be the length of tourL i that ith UAV travels. Then it's obvious that the time of completion T max(`= 0 ;` 1 = 1 ` k = 1 ). However,L and allL i is a solution with cost larger than the optimal cost which means (k + 1)T optimal cost (r) p 0 1 ZZ R 2 p f(x) dx where the second inequality comes from theorem 19. This completes our proof for the rst inequality. Comment This theorem actually shows that the system with k UAVs are equivalent to the system of a truck with speed p k 0 1 . This result gives us how ecient the system is in general which is increasing by the speed of p k as we increase the number of UAVs. One extra information we get here is that(r;k)(r; 1) p k=(k+1) as it's obvious that(r) =(r; 1), so we can actually use (r; 1) for the lower bound rather than using the theoretical lower bounds we deduced which may sometimes be better. 4.2 Computational Results In this section we are conducting computational experiment. It will be done in the unit square with uniformly distributed demand, which we use to estimate the coecient (k;r) introduced in previous section. 4.2.1 Experiments using Euclidean distances In our rst experiment, we sample n = 500 points uniformly in the unit square, we x a truck speed of 0 = 1, and we allow the UAV speed 1 to vary with 1 2f1:5; 2; 2:5; 3; 5g. Notice that the method used for solving truck with one UAV can't be used directly here as we also need to 52 decide points are going to be visited by which UAV. Because of the diculty of solving a given problem instance to optimality, we apply the following intuitive heuristic rule: Stage 1: { Assume initially that points p 1 ; ;p n are ordered according to their optimal TSP tour.The targets of the k UAVs are initially set to the points p 1 ; ;p k . { The truck always drives towards the next point that is not currently a target of any of the UAVs (for example, initially, the truck drives towards p k+1 ). { Each UAV ies to its target and returns to the truck, at which point its new target becomes the next point that is not currently a target of any of the UAVs. Stage 2: { After recording the sequence of events that occurred in Stage 1, we have a list of the points that each UAV visited, which is written as P i , for i2f1; ;kg. { Solve the following convex optimization problem, which nds an optimal horse y tour for xed ordered assignments P i : minimize x1;;xn+1;t1;;tn+1 t n+1 s:t: t r t q + (kx q p q k +kp q x r )= 1 8(q;r) :q =P i j ;r =P i j1 t i t i1 +kx i x i1 k= 0 8i2f2; ;n + 1g t 1 = 0 x 1 =x n+1 Figure 4.1 shows two of the solutions that are obtained when using this procedure. Table 4.2 shows he estimates of the coecient (r;k) as computed by performing the above procedure 50 53 (a) (b) Figure 4.1: Two dierent horse y tours with k = 1 and 1 = 1:5 (a) and withk = 3 and 1 = 3 (b); in the later, the dierent colors of the points p i (black, gray, and white) indicate which of the three UAVs visits that point. 1 nk 1 2 3 4 1.5 1.02 0.88 0.84 0.80 2 1.00 0.93 0.86 0.78 3 0.95 0.89 0.85 0.74 5 1.02 0.92 0.83 0.80 Table 4.2: Estimate of ( 1 ;k), the proportionality term introduced in theorem 24 54 times for each pair of 1 andk (hence, 50 4 4 trials in total). Each estimate ^ is obtained by computing the TSP tour of the sampled pointsp 1 p 5 00 as well as the horse y tour that results from performing the heuristic above, and then setting ^ = Service time without UAVs Service time with UAVs p k 1 = 0 the values shown in the table are simply the averages of the 50 dierent trials. Note that, for xedk, the estimates of are more or less consistent for varying 1 , and that seems to decrease in terms of k. The results shown are for the most part consistent with our prediction that 0:2111 (r;k) 2 q p 2(1 + p 2). Thus, by taking averages over the columns of table 4.2, we estimate that for k = 1, we have = 1:00, for k = 2, = 0:91, for k = 3, we have = 0:84, and for k = 5, we have = 0:78. 55 Chapter 5 Horse y Problem with Multiple Objectives In previous sections, we only considered time of completion as our only objective. However, many other factors may be incorporated at the same time when making decisions in reality. As a prospect way of solving the last mile delivery problem, horse y system has at least two benets over the traditional delivery with trucks only system: cost eciency and faster delivery. As shown in theorem 15, the horse y system with one UAV can be imagined as a new vehicle with speed p 0 1 when optimal solution is achieved. A similar result is true if objective function substituted by total cost as shown in theorem 19. It means that when solving the horse y problem with total cost as objective function, optimal value will have the approximation of p c 0 c 1 RR R 2 p f(x) dx if c 0 represents the cost of truck for traveling unit distance and c 1 represent that for UAV. When making decisions, we can only optimize one of them. In this section, we claim a way to combine the two objective functions and state a fact that when one of the objective functions is optimized, the other one will also reach a value that is not bad. It makes sure that the two benets of using UAV can be taken simultaneously. In addition, a upper bound will be given on how bad the other objective can be. Lastly, we estimate the performance by of a good under the fact r close to 1. In order to distinguish the dierence between total cost and time of completion, we assume the in theorem 15 to be t and the corresponding one for total cost to be c . It's not dicult to see when ratio of cost is inverse of that of time, t c 2 t . 56 5.1 Theoretical Results A natural way of solving the problem will be making a new objective function by combining the two functions together with a coecient 0 1 making trade-o between total cost and time of completion: minimize x1;:::;xn;2Sn p n P n i=1 c 0 kx (i) x (i+1) k + 2c 1 kx (i) p (i) k p c 0 c 1 + (1) p n p 0 1 n X i=1 max 1 0 kx (i) x (i+1) k; 1 1 kx (i) p (i) k +kp (i) x (i+1) k (5.1) Changing can balance the weights between total cost and time of completion. And if n is large enough, the both terms in the objective function are constants from theorem 15 and 19. As the optimal solutions generally can't be achieved at the same time, it's obvious a lower bound for the above solution will be c + (1) t . Next, we are going to get a asymptotic upper bound based on theorem 15. As shown in Figure 3.3, when we adopt the optimal solution where total cost reaches the optimal value and its asymptotic value to be c p c 0 c 1 RR R 2 p f(x) dx, truck travels at most c p c0c1 RR R 2 p f(x)dx c0 and drone with distance at most c p c0c1 RR R 2 p f(x)dx c1 , the time of completion is at most ( p c1 p c00 + p c0 p c11 ) c RR R 2 p f(x) dx. When the optimal solution for time of completion is taken, the same discussion shows the total cost of at most ( p 0c1 p 1 + p 1c1 p 0 ) t RR R 2 p f(x) dx. Do the similar thing for time of completion when using the optimal solution of cost will lead to the following theorem Theorem 26. Assume r = p c11 p c00 , the optimal solution of problem with combined objective func- tion shown above has a lower bound of c + (1) t and upper bound of min( c + (1)(r + 1 r ) c ;(r + 1 r ) t + (1) t ). This theorem shows that when optimal time is reached, cost will be within (r + 1 r ) t c multiples of optimal cost, while when optimal solution for cost is taken, time of completion will be within (r + 1 r ) c t multiples of optimal time. Notice c ; t can be bounded by the values in theorem 15 thus can be regarded as constant. Let's rst analyze what it tells us when optimal value for 57 time is used. In this case, the time of completion is theoretically t p 01 , the total cost is at most (r + 1 r ) t p c 0 c 1 . If r is a number close to 1, it will deteriorate cost within some constant factor r + 1 r . Otherwise if r is large which means 0 << 1 , the term is almost q 1 0 t c 1 which can be much larger than the truck case. Lastly, r is small where c 1 << c 0 , cost is almost q 0 1 t c 0 . In this case, cost can be far away from optimal value while is better than truck case with a factor q 0 1 . As to decision making, with the discussion above, the bound r is a constant close to 1, we can choose to be any number and both objective function will be within factor of (r + 1 r ). When r is much larger than 1, we prefer the value of to be close 0. In this way, optimal cost is reached and time of completion can also be improved by a subtracting factor of q c0 c1 where compared with truck delivery. And when r is much smaller than 1, setting to 1 will give us optimal cost and time improved by factor of q 0 1 . The discussion here gives us some insights on delivery with horse y system. First, the horse y system can improve both delivery time and cost when chosen appropriately. Second, if we focus only on time of completion and keep on improving speed of UAVs regardless of its cost, the system may behave much worse than truck delivery only in cost. Third, it's better to improve the bottleneck part of UAVs even if you focus on the other objective function. For example, if we want to optimize time, however UAV already has 1 >> 0 . Keeping on improving the speed of 1 will dramatically deteriorate performance of cost and improving q 1 0 under 1 >> 0 is dicult . However, if we can lower the cost of UAV, keeping cost at the same level can improve the time of completion correlated with q c1 c0 , it's easier to work on. Lastly, we estimate the performance of a good by assuming that it will make P n i=1 c 0 kx (i) x (i+1) k + 2c 1 kx (i) p (i) k c p c 0 c 1 p 0 1 t n X i=1 max 1 0 kx (i) x (i+1) k; 1 1 kx (i) p (i) k +kp (i) x (i+1) k 58 It's a reasonable assumption as when is close to 1, the term on left hand side is close to 1 and right hand side term is larger than 1, it's reversed when is close to 0, noticing the upper bound of theorem 26min( c +(1)(r + 1 r ) c ;(r + 1 r ) t +(1) t ) can be written as format ofmin( c +(1)(r + 1 r 1) c ;(r + 1 r 1) t + t ) . The rst term is decreasing and the second one is increasing, thus the maximum value is reached when = t+(r+ 1 r )c (t+c)(r+ 1 r 1) and the value is (r+ 1 r +1)ct c+t . Thus if there is some that makes the two terms in the objective function for 5.1 are close, assume the optimal value is opt r . We will have t opt r (r+ 1 r +1)ct c+t or opt r (r+ 1 r +1)c c+t . It means that the solution it gives are pretty close to both time of completion and cost almost equal to the optimal value they can reach at most. In summary, we claim that from theoretical perspective, when r is much larger than 1, it's reasonable to make objective function as cost. When r is much smaller than 1, take time as objective function is a good choice. And 5.1 can be used with choice from 5.1 if the above choices are not enough. 5.2 Computational Results In our previous section, some theoretical work considering both time of completion and total cost shows the fact that under the scenario that r = p c11 p c00 is neither too small or too large, solution of optimization under one objective function will not deteriorate the other too bad. However, there are two questions of importance in practice remaining unanswered. First, the performance of optimal solutions for 5.1 on both objective functions as changes. Second, how we can choose a reasonably well which will make a trade-o between time of completion and cost. In this section, we are trying to answer the two questions by solving the optimization problem 5.1 with the points visiting order from the method as before. One thing worth mentioning here is that although we are solving problem 5.1 directly, the optimal value itself is not what interests us as the value can't be enough to measure whether the solution is good or not. Thus we rst 59 (a) r = 0.2 (b) r = 5 (c) r = 1 Figure 5.1: For each graph of the above, problem 5.1 is solved eleven times for each = 0; 0:1 ; 1:0. y-axis is time(x;)=opt-time and x-axis is cost(x;)=opt-cost. (a) is the case that 1 0 = 2 and c1 c0 = 0:1, (b) is where 1 0 = 10 and c1 c0 = 0:5 and in (c), 1 0 = 2 and c1 c0 = 0:5. solve the problem for its optimal solution and compare it with the respective optimal solutions for optimizing objective function as only time of completion and total cost. Mathematically, if the optimal value of problem in denition 3 to be opt-time and that of denition 16 to be opt-cost. Solving 5.1 for the optimal solution x; and calculate the respective objective functions with it for time and cost, say time(x;) and cost(x;), the values of time(x;)=opt-time and cost(x;)=opt-cost are what we consider here. Simulation results are shown in gure 5.1. The graphs are all having the shape of a convex curve. However, the curve is closer to the most desirable point (1; 1) when r = 1. It means that the optimal solutions of cost and time are closer than the cases that r is much less or larger than 60 1 0 c1 c0 time(x;)=opt-time cost(x;)=opt-cost 2 0:04 0:4659 1:2866 1:2811 25 0:5 0:4698 1:6059 1:6019 25 0:04 0:4094 1:0633 1:0629 10 0:2 0:4211 1:1329 1:1326 Table 5.1: Optimal values comparisons with dierent rations for multiple objectives problem 1. Another thing is that points are evenly distributed when r = 1 while are more biased toward top when r is large and toward bottom when r is small. It other words, if we want to get the point closest to (1; 1) in the curve, needs to be small when r is small and large when it's large. This is because when r is large, time has a more rapid change thus doesn't need to be large to be signicant, it's discussed theoretically in previous section. One thing that is not talked about is the convexity of the graph. The benet of this property is that the choosing a good will be much better than the theoretical bound (r+ 1 r +1)c c+t deduced. From simulation results, the curvature of the graph is larger when r is further away from 1. It means the its actual performance may be much better than theoretical upper bound even r is not close to 1. Actually, from 5.1b, a good choice of can give a factor of 1:3 for both time and cost even if r = 5. Thus we propose a method that makes binary search on until the two ratios are close enough, in order to see its capacity, we tested it on the following three cases under n = 1000, simulation results show the table From the table, rst we claim that the ratio r is highly related with results of the method, with r close to 1 , both time and cost can reach a satisfying level. When r is either too large or small, it may not be able to obtain both objective optimal. For example, in the worst case where 1 0 = 25 and c1 c0 = 0:5. Cost is even worse than the truck delivery by a factor of 1:6019= p 2 1:14, although time having a subtracting factor of 5=1:6059 3:11. Lastly, to verify our claims on always improving bottleneck, we compare the second column in the table with the last column in the table. Imagine the original UAV has speed ten times of truck 61 while unit traveling cost of only half truck. Two options are at our hands, either improve speed by a factor of 2:5 or decrease the traveling cost by a factor of 2:5. Compared with truck only delivery, the rst option can make our system 3:11 times faster while 1:14 times more expensive. For the second option, it will make the system 2:79 times faster and 1:97 times cheaper. In practice, the second option is denitely more desirable, in accordance with our discussion before. 62 Chapter 6 Horse y Problem with Battery Life Constraints In previous sections, the horse y problem where truck and UAVs combined together is discussed. However we assume UAV can travel a distance as long as it wants. This is not the case practically as UAV is limited by a lot of constraints in reality. First is the battery issue where commercial UAVs nowadays usually only have minutes of battery life which is short. Second is because of the regulations in United States that UAVs can only travelled within sight. The last one is that some cases, people need to control the UAVs and thus may not be possible if UAV is far away. Above all, it's reasonable to think that UAVs will have the travel distance limit in practice. A question raised then is how much the optimal solution will be aected when UAV is limited to travel no more than some distance d. It's obvious that if the number of customers is larger, UAV will tend to travel a shorter distance. So it's from common sense that if d is small enough, the optimal value will be the same as that of original horse y problem. And after the threshold, the optimal value will increase. In this section we are discussing something related with this. 6.1 Theoretical Results This problem can be formulated easily as an integer programming problem following the method for denition 4 by adding one constraint for battery life. One thing to mention here is that 2n 63 points representing loading and unloading locations are needed here while just n points are used before because in this case, UAV has to be carried by truck if it has too short battery life minimize n X i=1 (t 2i1 +t 2i ) (6.1) t 2i1 kx 2i1 x 2i k 0 8i (6.2) t 2i1 M(1c ij ) + kx 2i1 p j k +kx 2i p j k 1 8i;j (6.3) t 2i kx 2i x 2i+1 k 0 8i;j (6.4) M(1c ij ) +kx 2i1 p j k +kx 2i p j kdf8i;j (6.5) x 2n+1 =x 1 (6.6) X i c ij = 18j (6.7) X j c ij = 18i (6.8) c ij 2f0; 1g; x i 2R (6.9) Recall that in the case where we don't consider the limit of UAV's travel distance, the rst lower bound is from solving an optimization problem in the proof of theorem 12. Here, to modify it to t into our current problem considering battery life of UAVs, an additional constraint (the last constraint in bold font below) is added and solving the following problem will give us the lower bound. minimize ` max( ` r ;f(`)) subject to 2(d + p 2 p n )` 1 From the proof of lemma 9, optimal value is reached at ` = 1 2 r( q 2 + 1 r p 2) p n, thus the lower bound will still hold if d> ( q 2 + 1 r )= p n. 64 The second lower bound in discussion following theorem 10 is from minimize `i max( ` r ;` d ) subject to ` +` d tsp p nS Here there is an additional constraint ` d d p n as UAV can travel at most d q 1 n for visiting one custom, and the bound won't be aected a lot (from tsp p n to tsp 1+r p n) if d 1 1+r tsp (or d 1 2 p r tsp ). Lastly, the problem of nding an upper bound becomes minimize K p K r + 2cn q S K subject to p K p 2n d The upper bound won't change compared with the problem without constraining battery life if d p c r . Remark. The discussion above doesn't directly shows that when d is large enough, the optimal value won't be aected. But it illustrates the fact that whend satises the above condition, all the discussions in problem without the constraint will hold. In other words, no matter what values if 1 ; 0 are, the dierence between optimal solution with and without battery life constraint will be within some constant factors. If we let the critical value of d to bed c , the discussion above shows d c q 1 nr . In practice, it means that improving the capacity of battery won't aect a lot when it passes some threshold value (the value that can guarantee UAV to travel more than distance d c ). We will analyze relationship between d and time of completion (instead of upper and lower bound of that) and the critical value in the next section through numerical experiments. 6.2 Computational Results In our previous section, we analyzed lower and upper bound related with distance limit of drone d. However, the direct relation between optimal value which is the time of completion and d was 65 (a) TSP (b) Horse y (c) Horse y-BL Figure 6.1: Eight points are uniformly generated from unit square, (a) is the tour of TSP which gives 2:689 time units, (b) and (c) are both from solving the integer programming formulations before respectively for horse y and horse y-BL under this situation. Here, truck is set to be half the speed of drone and for horse y-BL the distance that drone can travel is no more than q 1 n0 = 0:5. Horse y problem gives optimal value 1:792 and that for horse y-BL is 2:146. 66 not analyzed. In this section, numerical experiments will be on it. It will give us quantitative results on how optimal value will be aected when ability of UAV is limited. The algorithm we are working on follows the procedures as that in the previous subsection, draw sample points, run Lin-Kernighan algorithm and solve a convex optimization problem. The only dierence is the last convex optimization problem where there is an additional constraint and 2n variables are needed for loading and unloading points. minimize x1;:::;x2n;2Sn n X i=1 max 1 0 kx 2i1 x 2i k; 1 1 kx 2i1 p (i) k +kp (i) x 2i k + 1 0 kx 2i x 2i+1 k (6.10) s:t kx 2i1 p (i) k +kp (i) x 2i kd (6.11) where is the xed permutation representing order of customs being visited and x 1 =x 2n+1 . An example is shown in gure 6.2. It is obvious that optimality of the problem is decreasing when d increases. But the gure also shows that after some threshold value of d, the optimal value keeps almost the same. We dene d c to be the smallest value that will make the optimal value within one percent of that for problem without battery life constraint. In practice, when making trade-o between battery life capacity and optimal value, d c can be taken as a choice. It will keep the optimal value within one percent of the best can be reached and from discussion, d c won't be very large in order to reach this goal. We have discussed how d will aect our upper and lower bounds in previous section which shows that when d is larger than some value, the same upper and lower bounds still hold. And this value is proportional to q 1 n0 . Thus we guessd c will have the same property and will analyze it through simulations. d c is found in our simulations by binary search and it terminates when the interval remaining is of length less than 0:001 and simulation results are shown in gure 6.3. One thing to mention here is that we are doing linear regression ond c with respect to 1 p nr , but is theoretically incorrect 67 Figure 6.2: The gure shows the change of optimal value for horse y problem when the limit on the distance that UAV can travel changes where n = 1000 and 1 = 0 = 2. The points are plotted at d = 0; 0:2 q 1 n0 ; 0:4 q 1 n0 ; 3 q 1 n0 . The optimal value decreases extremely fast at rst and then the graph becomes at. as it's obvious that when n!1 orr = 1,d c ! 0. Andd c should be bounded by p 2 2 . And some of those properties can't be satised. Despite of that, our linear regression has its meanings. Figure 6.3a shows some linearity with respect with 1 p n when n!1 which gives us a reasonable way to bound d c by c q 1 n0 and from gure 6.3b, its linearity also guarantee the upper bound although with negative intercept. In summary, the simulation ensures in some sense that c q 1 n0 is a reasonable upper bound and the number from simulation is 2:2579. In order to verify the correctness of the upper bound 2:2579 q 1 n0 . We randomly generate n from 20 to 500 and 1 n0 from 1:0 to 10:0 and ndd c to compare with the upper bound. The result is shown below 68 (a) dc v.s. 1= p n (b) dc v.s. p 1 = 0 . Figure 6.3: (a) is the values of d c with respect to dierent 1= p n where n = 100; 200; 2000 under xed r = 0:5. (b) is d c with respect to p 1 = 0 under n = 1000. Both graphs show some linearity. Run linear regressions for data in (a) gives d c = 0:0014 + 2:0469 1 p n or d c = 0:0014+1:4474 1 p nr . Do the same thing for (b) givesd c =0:0575+0:0714 1 p r ord c =0:0574+ 2:2579 1 p nr 69 n 1 0 d c upper bound percentage dierence 391 7:6882 0:2448 0:3166 29:33% 313 7:6603 0:3025 0:3532 16:77% 57 8:7734 0:6810 0:8858 30:08% 460 6:3243 0:1793 0:2647 47:64% 140 6:4661 0:4571 0:4852 6:16% 79 4:0232 0:3597 0:5095 41:67% The upper bound given above are correct for all the six samples generated although with some gaps. The gap is due to the randomness of samples, for example when n = 1, d c can vary from 0 to p 2 2 , and also the inaccuracy of d c represented as function of c q 1 n0 . From the examples, it shows that the upper bound is robust and with reasonable gap. If more accurate numbers are needed, simulations for d c can be done for the specic cases needed. 70 Chapter 7 Horse y Problem with Truck Unreachable Regions An important benet of drone is its broader reachable areas than truck. Sometimes, we may restrict the region that truck can travel in some scenarios. For example, drones can make deliveries to rural areas where infrastructure like highway is not covered. Companies like JD.com in China are taking the benet to boost consumption in those areas by leveraging delivery service [6]. Another scenario is the post-disaster monitoring in areas that are either with potential dangers to humans or becoming unreachable as infrastructure destroyed in the disasters. For example, drones were used to test radiation levels at the site of the 2011 nuclear disaster site in Fukushima [8]. In summary, this assumption is true whenever there are some restrictions on the region that truck can travel. Unlike the problems discussed before where UAVs are limited by battery life or charging time, truck is here limited within some reachable ares while customs can appear in areas that are unreachable. It's a generalization of horse y problem where truck can reach anywhere and also of pure drone delivery where truck is limited at the depot. In addition to 0 and 1 , size of reachable area and even location of it may aect optimal solution. Take the problem in gure 7.1 as an example, the small areas located at the center is far better than it being placed at any corner if customs are uniformly distributed. 71 Figure 7.1: An example of horse y problem with areas that are unreachable for trucks. In the graph, grey area can't be reached by truck while four customs are all located at unreachable areas. We will show our theoretical results later that for the problem, the cost asymptotically only relates with the average distance from the area unreachable by truck to the distance reachable and area portion of area that are reachable, here we assume the average distance to be d avg and the percentage of area reachable by truck is p (which means the area of the region pA if A is the area of the whole region) in this chapter. 7.1 Theoretical results It's obvious that if we restrict truck at a xed location, the truck can be seen as a depot and we have a pure drone system and if we don't have any restrictions, the system will be the horse y system we mentioned in our previous chapters. But asymptotically, we will show that our current system is actually more like a pure drone system. We will give lower and upper bounds and then asymptotic results which are not hard to obtain. Before heading to that, we do some discussions. Notice that from results we had before, to visit the points inside the truck reachable area, the time of completion will be only(r) p n, however, when n is large, the number of customs outside reachable area are proportional to n. Those customs 72 will be visited by UAVs with average distance d avg . Thus it's obvious that asymptotically, time of completion will be proportional to n if some restrictions were given to truck. The results and related proofs are as follows: Lemma 27. [Lower Bound] Suppose that X =fX 1 ;X 2 ;:::X N g is a sequence of random points i.i.d. uniformly distributed within a compact planar regionR 2 with area A. Assume that only regionsR in with area pA are reachable by truck, and the average distance from theR in to re- maining regions is d avg . Then with probability one, the optimal time T N of horse y problem with one UAV through X satises T N (r) p 0 1 p p NA + 2n(1p)Nd avg = 1 +o( p N) where (r) is dened same as theorem 15. Proof. The proof of the theorem is quite straightforward, we handle customs inside and outside R in separately. Notice that from theorem 15, if we are only visiting points withinR in , the optimal time will be (r) p 01 p pNA. For the points outsideR in , as truck can't travel outsideR in , UAV will travel at least d(x;R in ) where x is the point to be visited. Thus those points will take at least 2n(1p)NAd avg = 1 which completes our proof. Lemma 28. [Upper Bound] Following the notation in the above lemma, T N satises T N (r) p 0 1 p p NA + 2n(1p)Nd avg = 1 + (` +d Rin )= 0 +o( p N) where c is dened in theorem 11, ` is the contour of regionR in and d Rin is the longest distance of any two points inR in . Proof. Without loss of generality, we can assume the area of whole region is 1. As we are nding upper bound this time, we can't directly separate point insideR in and outside. Thus we are 73 Rin R Figure 7.2: This is an illustration of the scheme we adopted for upper bound. First all points insideR in will be visited with the optimal solution for the horse y problem. After all those points visited, a shortest tour to the contour ofR in will be used and then truck will travel along the contour and each point outside will be visited by UAV from the closest point inR in . following the idea of proving theorem 11 where we give a scheme that will visit all points. The basic idea is that we rst visit all points insideR in with the optimal solution for horse y problem and then choose the shortest path to the contour ofR in (see gure 7.2) and travel the contour. UAV will visit all points outsideR in from the closest point inR in . In addition to time needed to visit all points insideR in , truck travels additional ` +d Rin distance where ` is the contour of regionR in . It's obvious for time of completion T : T N (r) p 0 1 p p NA + 2n(1p)Nd avg = 1 + (` +d Rin )= 0 +o( p N) which completes our proof. With the above two lemmas, we will be able to get the asymptotic results similar to theorem 12 as follows: 74 Theorem 29 (BHH for horse y problem with truck unreachable area under uniform distribution). Following the notation in lemma 27, with probability one, the optimal time T N of horse y problem with one UAV through X satises lim N!1 T N = (r) p 0 1 p p NA + 2n(1p)Nd avg = 1 +o( p N) where (r) is only related with r = 0 = 1 and 0:2103(r) 1:7504. 0 , 1 are speed of truck and UAV respectively. For this theorem we assume the points are uniformly distributed, the lower bound can be easily extended to the case of uniform distribution as we have similar results for general density function in theorem 15. We will thus have a similar result under general density function: Theorem 30 (BHH for horse y problem with truck unreachable area under uniform distribution). Suppose that X =fX 1 ;X 2 ;:::X N g is a sequence of random points i.i.d. distributed with an absolutely continuous function f within a compact planar regionR 2 with area A. Assume that only regionsR in with area pA are reachable by truck. Then with probability one, the optimal time T N of horse y problem with one UAV through X satises (r) p 0 1 p p N Z Z R p f(x)dx + 2N(1p) Z Z RnRin f(x)d(x;R in )dx= 1 +o( p N) limT N 2(r) p 0 1 p p N Z Z R p f(x)dx + 2N(1p) Z Z RnRin f(x)d(x;R in )dx= 1 +o( p N) where (r) is only related with r = 0 = 1 and 0:2103(r) 1:7504. 0 , 1 are speed of truck and UAV respectively. 75 center area n without projection with projection 0.5, 0.5 0.25 100 36.8240 35.9396 0.25, 0.25 0.25 100 63.9317 63.4013 0.5, 0.5 0.01 100 96.7562 96.6030 0.5, 0.5 0.25 300 104.7312 103.4857 Table 7.1: Algorithms comparison with and without projections 7.2 Computational Results In this section, we are going to solve the problem computationally. The basic structure we are using to solve it is quite similar with our previous one as the only dierences here is that now the truck location will have to be withinR in instead of the whole regionR. We only need to add those constraints when solving it. In other words, we are also rst solving traveling salesman problem to get an order then use it as an input for the second optimization problem. One change we make here is that before solving TSP using Lin-Kernign algorithm, a projection for all points outsideR in to the closest points inR in is implemented. The intuition is that as points outsideR in will be visited by UAVs, but the benet of visiting two adjacent points will be from truck. In our current settings, as truck can't visit points outsideR in , points are viewed as close by truck when their projections are close. Some results comparing two algorithms are shown in table 7.1. Points are generated randomly in a unit square andR in will be a square inside. In general, the results with projection are slightly better than that without projection, although only with 1% improvements. There are some extra information about values in the table and in general they are the same as explained in the theoretical parts. We will do further discussion later. Here we show some graphs on how the solution looks like in gure 7.3. 7.2.1 Time of completion estimation From gure 7.3, we can get some intuitions on how to estimate the time of completion on the problem with truck unreachable area. One important thing worth mentioning is that for points 76 (a) (b) (c) Figure 7.3: 100 points are generated randomly in the unit square, in (a) truck can only travel at the center square with side length 0:5, for (b), truck can also only travel at a square with side length 0:5 but located at the bottom left, (c) also restricts truck at a center square but with side length 0:8. It's obvious that truck travels a large amount of time at the boundary of the restricted area. 77 lower bound numerical result upper bound 12.1781 12.9145 14.6781 14.0402 11.0821 16.5402 4.8385 5.5349 8.8385 6.7446 7.4419 10.7446 57.7243 57.4218 60.2243 67.0348 69.0748 69.5348 16.0874 16.6979 20.0874 25.6178 25.1652 29.6178 Table 7.2: Comparisons between numerical results and theoretical bounds outside the restricted area, the solution is almost taking the shortest tour to the restricted area. Combining the theoretical results in the previous section, we are thinking of making estimation by lower bound: T =(r)p p N= p 0 1 + 2N(1p)d avg = 1 As we are doing simulation in unit square and will restrict truck also to a square inside, thus the upper bound will be: T =(r)p p N= p 0 1 + 2N(1p)d avg = 1 + 5` where ` is the side length of restricted area. We will test the case where r = 0:5 which means UAV has twice the speed of truck, and (r) 0:9 tsp 0:648. Numerical results are shown in table 7.2 and gure 7.4, most of the cases will be within error of 10% and around 5%. And upper bound are larger than the real number in all the cases. Notice that the upper bound and lower bound only have constant dierence, so it will be more accurate when the number of points is larger. In this simulation, we only test cases where n is between 100 and 500. But even at this range, the estimation is satisfying. Thus using the lower bound to estimate the time of completion will be a good estimation in general. 78 Figure 7.4: Eight cases are generated with dierent area, location of restricted area for truck and number of customs. Left bar is the theoretical lower bound and middle bar is the solution of our algorithm introduced, the right bar is the theoretical upper bound. 79 Chapter 8 Discrete Horse y Problem In previous chapters, horse y problem is tackled with theoretical bounds and a heuristic algo- rithms. However, all our problems assume that launching points for loading and unloading drones from the truck can be anywhere in a region. This is not always true especially nowadays there are many regulations and limits on drones. For companies like JD.com, drone deliveries are carried out by launching drones at some xed locations and packages are manually loaded to drones. Mathematically, it means that possible launching sites are predened points or areas rather than continuous and can be any point in the region as discussed before. And in this paper, we are handling with this discretized horse y problem. One thing that we need to mention is that this problem is dierent from the problem where we limit the region that travel can travel. Here we truck can still travels to any where it wants, but launching sites are restricted to some selected points, although it may also exert some restrictions on truck as it has to visit between launching sites. After the discussions in later sections, we will see that those two problems are very dierent both in theoretical and computational results. 8.1 Problem Formulation First we are formulating the problem, this formulation won't be used directly for our later results, but it gives us a sense of the diculty of the problem. 80 Assumep 1 ;p 2 p K are points (customs) to be visited andq 1 ;q 2 q N are possible launching points, the problem is then to nd a permutation of 1; 2;K determining order of points to be visited and a one to one mappingf from domain 1; 2;K to range 1; 2;N determining the launching points to be used. Let 0 ; 1 to be speed of truck and drone respectively, the discretized horse y problem can be dened as follows: minimize ;f n X i=1 max 1 0 kq f(i) q f(i+1) k; 1 1 kq f(i) p (i) k +kp (i) q f(i+1) k (8.1) where we assume f(n + 1) =f(1) and (n + 1) =(1) for simplicity of notation. One thing to mention here is that one launching site can be used repeatedly, an extreme case is that only one launching site exists and our system delegates to the pure drone delivery. The problem is obviously more dicult as an additional map f is needed even with known. The common notation is to use N as the number of customers and K as the number of launching sites. We are using it reversely here for the reason of the formulation later. Formulating this way, it will be dicult to solve other than enumeration. Like before, a natural way of thinking is to formulate it as an integer programming problem. However, we will see that it's not necessary here as we are going to transform the problem to a asymmetric generalized traveling salesman problem (AGTSP). 8.1.1 Asymmetrical traveling salesman problem formulation To transform the problem to a specic case of AGTSP. First, each of launching pointsq 1 ;q 2 q N can be viewed as the same point with K replicas. In other words, q i becomes x ij ;j = 1; 2;K which represents that customer j will be visited from q i . And imagine that points x ik ;i = 1; 2; ;N are one cluster for any specic k. The distance between x ij and x st is dened as d(x ij ;x st ) = max 1 0 kq i q s k; 1 1 (kq i p j k +kp s q j k) 81 which is the time of truck traveling from q i to q s with drone visiting p j , thus our problem now is to nd a Generalized TSP tour that will visit exactly one point in each of the cluster. That is to say, discretized horse y problem is equivalent to a generalized ATSP with K clusters and each cluster contains N points. Notice here that K is actually the number of customers and N is the number of launching points. After this formulation, we can then use the heuristic algorithms solving AGTSP and transform back to our original problem. We will use one algorithm that experimented to work well in computational section later. 8.2 Theoretical Analysis A natural question is how the number of launching points will aect the optimal value of the problem. And as a follow up and an important application is where should we place the launching points. And this section is on theoretical analysis of the problem that will also answer the above questions to some extend. In order to do theoretical analysis, custom points p 1 ;p 2 p K are now assumed uniformly distributed in a compact regionR with unit area. Launching sites q 1 ;q 2 ;q N are variables to be determined to minimize the expected time to serve customs. In other words, dene the optimal value of 8.1 to beoptTime(p 1 ;p 2 ;p K ;q 1 ;q 2 ;q N ), we will have the following problem regarding choosing the locations fo launching sites: minimize q1;q2;q N E p1;p2;p K opt(p 1 ;p 2 ;p K ;q 1 ;q 2 ;q N ) (8.2) lower and upper bounds will be given. For simplicity, dene the optimal value for 8.2 to be opt(K;N). 82 length d Both truck and UAV travel at most 2d extra distance Figure 8.1: This is an illustration of moving launching site a distance of d. The the box point is the new launching site, both truck and UAV will travel at most 2d extra distance from comparing the new and old tour. Thus the extra time will be at most 2d=(1= 0 + 1= 1 ). The rst lower bound can be achieved by relaxing launching sites to be anywhere in the region, this is exactly the horse y problem discussed before. Thus when K!1, we will have opt(K;N)= p K! 1 p 01 (r) where r is the ratio between 0 and 1 . This lower bound is not good when N is small as launching sites can hardly be relaxed to be any place in the region. In this case, we are only considering the time drone travels. Notice that the best case is when each launching site is in charge of the area of a circle with it as a center and radius r satisfying r 2 N = 1, the time drone travels at this case is opt(K;N) 4K 3 p N1 . Theorem 31 (Lower Bound). opt(K;N)= p Kmax( 1 p 01 (r); 4 p K 3 p N1 ) Upper bound can be reached by giving some schemes for visiting the points given. The method based on the horse y tour in continuous case will give Theorem 32 (Upper Bound). opt(K;N)= p K 1 p 01 (r) + ( 1 0 + 1 1 ) p 2 p K p N . Proof. The scheme we are using is to adjust the launching sites from solving continuos horse y problem to the closest possible launching sites. This will move the launching site at most 1= p 2N. And it will incur at most extra time of q 2 N (1= 0 +1= 1 )(illustrated in gure 8.1) which compeltes our proof. 83 In summary, with the above lower and upper bounds, we will have the following asymptotic results for three cases: (i)N = (K), discretized horse y problem have the same result as continuous one withopt(K;N) = p K! 1 p 01 (r). (ii) N = o(K), opt(K;N) will have the optimal value of (K= p N) which will deteriorate the ( p K) in continuous case. (iii) N = (K), opt(K;N) has the order of ( p K) but may cause a increase of converging constant. 8.3 Computational Results In this section, we are going to introduce three heuristic algorithms that can solve the discrete horse y problem. Each algorithm has it own benets in some cases. Some algorithms are closer to optimal solution but takes huge amounts of computational and meomry resources, while some algorithm takes less resources but are with sligtly worse results. 8.3.1 Algorithm 1: Directly Solving AGTSP The rst algorithm we are implementing is to use the formulation we descibed before and solve the problem directly with Lin-Kernighan-Helsgaun Algorithm which proves to work well within 5% errors for most cases. In summary, the problem contains two parts: stage 1: Formulate the problem as AGTSP withK clusters and each cluster containsN points where K is the number of points to be visited andN is the number of launching points, calculating the corresponding distance matrix. 84 Figure 8.2: This is an example of solution of discrete horse y problem by solving AGTSP directly. We generated 20 random points in a unit square to represent the points to be visited and another 20 random points, both uniformly distributed, to represent possible launching sites. The launching sites are circled points and visiting points are black dots in the picture. It's obvious that unlike continous horse y problem, truck stop at some launching sites as no other better launching sites are nearby. statge 2: Solve the problem formulated above with Lin-Kernighan-Helsgaun Algorithm and transform the solution back to original problem. An example is shown in Figure 8.2. This problem solves AGTSP of size NK at the second stage. So it works ne if we don't have many launching sites. But in reality, that may not be true. For example, in real delivery, it may possible that each small area may have some possilbe lauching sites, which means K has at least the same scale as N. In that case, even if N = 100, we have a AGTSP of size at least 10000 which takes huge amount of time and memory. It happens because each launching site becomes a candidate for every customer point. But from intuition, that may not be true. Some launching sites may not be even be chosen by any chance and some others may be a choice for only small amount of points. We will introduce another algorithm based on this intuition. 85 8.3.2 Algorithm 2 : solving AGTSP with launching sites pre-selection We rst brie y discuss the distance that UAV will travel on average. Notice that from theorem 12, the total time of completion will be (r) p n= 0 1 and as UAV can inhabit on the truck, UAV will travel at most (r) p n=r or on average a distance of (r)= p rn. Now if we assume that the launching sites are uniformly distributed, the neigbhorhood of (r)= p rn contains about (r) 2 k=rn clusters. It means that only aboutck=rn launching sites are needed for a specic point. So we are claiming the following heuristic algorithm with a parameter k representing the number of launching sites to use. stage 1: For each of the customer points, choose k nearest launching sites that are candidates of the point. Formulate the problem as AGTSP withK clusters and each cluster containsk points where K is the number of points to be visited. Each point is a pair of the customer point and the launching site for that point selected before. The distance between two points is dened as before. statge 2: Solve the problem formulated above with Lin-Kernighan-Helsgaun Algorithm and transform the solution back to original problem. It's obvious that the benet of the algorithm is that the complexity doesn't increase with the number of launching sites if k is selected as a constant. Before heading to our simulation results, we will rst head to another algorithm and compare the results from the two algorithms. 8.3.3 Algorithm 3: adjusting solutions of continous horse y problem The last algorithm we put forward is based on the solution of continuos horse y problem. However, it's obvious that some launching sites may not be the possible launching sites for this problem. A 86 basic idea is to adjust those launching sites to cloeset avialble one. This is also the idea we use to prove the upper bound of horse y problem. It's not dicult to see that if the number of lauching sites is larger, the result will be better as this approximation is more closer to optimal solution. The algorithm is also a two step algortihm as follows: Stage 1: Solve the problem as continuos horse y problem disregarding the required launching sites. (This step contains two steps as shown before, rst to solve a TSP using Lin-Kernighn algorthm and then a optimization problem). Stage 2: Adjust the lauching site to the cloest possible launching sites Notice that from the proof of upper bound, this method will have extra time of at most ( 1 0 + 1 1 ) 2K p N thus larger N will give out a better results. However, for the method before, our discussion shows that larger N will yeilds worse result because more launching sites need to be considered. To compare the two algorithms numerically, we generate 50 points uniformly distributed to be visited. And then compare the performance of the two algorithms by changing the number of launching sites. The launching sites are also generated uniformly. Results are shown in gure 8.3. The result is the same as we discussed before theoretically. Algorithm 2 performs better when N is small and algorithm 3 becomes better as N increases. Following our discussion above that the algorithm works well when N ckrn, here we have k = 4, r = 0:5, n = 50, and 250 is the number that algorithm 3 becomes better, c 2:5 in this case. 8.3.4 Optimal value with respect to dierent K We rasied a question about how the number of launching sitesK will aect the optimal value at the begining of this chapter. The theoretical part shows that the error is at most 2(1= 0 +1= 1 )K= p N. 87 Figure 8.3: Comparisons between algorithms solving discrete horse y problem. A case of 50 points is rst generated, then we generate laucnhing sites uniformly with dierent numbers in each test, the algorithm 2 and algorithm 3 described in this section are used to solve each test case where the parameter k for algorithm 2 is chosen to be 4 In the graph, x-axis is the the number of launching sites. The star point at the bottom right is the value for solving continuos horse y problem. Comparing with the optimal value of horse y problem which is (r) p 0 1 K. The percentage of increase is at most 2( 0 + 1 )=((r) p 0 1 ) p K=N. This result is enough for asymptotic result but with a huge gap in practice. Thus we are going to do it numerically for the increase percentage. As shown in our previous chapter that algorithm 3 works well under the case that N is large, thus we are going to do the experiments in this section with this algorithm. And in accord with theoretical results, p K=N will be used for the value of x-axis. Results are shown in Figure 8.4, the shapes of the graphs are almost the same and with the relative error about 3% when p K=N = 15. And relative error will be under 5% when p K=N = 10 for both case. It means that if K=N 100, the performance will be close to the continous one. 88 (a) K=100 (b) K=500 Figure 8.4: The two graphs show how the relative error of discrete horse y problem varies with respect the number of launching sites N, in (a) the number of points to be visited is 100 while for (b), the number is 500. The x-axis is the value of p N=K. 89 Chapter 9 Future Work 9.1 Horse y problem with charging time In previous chapter, problem related with battery life is considered. However, charging time is not put into the time. The problem without charging time as a factor can be taken as the circumstance that we change another battery once UAV is running out of battery. In practice, this may cause extra work of controlling battery life for our drone truck vehicle. Or we can change battery every time UAV comes back which makes truck carry a lot more batteries. Thus making UAV recharge on truck may be a good option especially when it can be charged very fast. Then the problem on how to schedule the delivery with charging time. A way to consider it is every time UAV comes back to truck, it has to recharge the battery it consumed before leaving again. It means the UAV is on full battery every time it leaves truck. To consider the charging time, we just need to add the following constraint to the problem with battery life constraints M(1c ij ) + (kx 2i1 p j k +kx 2i p j k)=ct 2i where c is the charging speed. The optimal solution of this problem can have the situation that truck will have to stay in some location for recharge. 90 Theoretically, an upper bound for time of completion can be found if we assume the speed of UAV is slowed down by a factor of 1 1+c . The lower bound needs some extra work which is not apparent. Computational results may also need to analyze the aect of charging time and whether a result like that in battery life where when d is beyond some value, the aect will be small. 9.2 Uncertainty of trac time In reality, trac time is not deterministic as the settings of our problems. Two questions may need to be answered in regarding with uncertainty of trac time. First, the uncertainty of trac aects truck more but has much smaller impact on travelling time of UAV. It's thus possible that the uncertainty in horse y system is smaller than truck only delivery. How much will horse y system reduce the variance compared to truck only delivery on dierent number of customs n and UAV speed 1 will be an interesting research topic. Second, deciding a delivery strategy under uncertainty is much more dicult. If a decision should be made ahead of time where we only know a distribution of trac time of truck, it will be a stochastic optimization problem where the minimum expected time of completion is needed. This scenario often happens when delivery time needs to be predicted ahead of time ( commence like amazon for example). On the other hand, when coming to practical delivery, it's more a dynamic problem where trac changes over time and we only know the current trac situation and have to make a decision with that and the decision may vary over time. For example, if the trac is extremely bad now, we tend to have a tour letting UAV travel more while in the condition of good trac, a tour with less UAV usage will be preferred. Both problems are intrinsically very dierent and much more dicult than the problems we are having before. 91 9.3 Better heuristics for horse y problem Finally, better heuristics is another important area that is worth to research. The heuristic we used in section 3.6 is not even a local minimum solution. For TSP, heuristic is always a two step algorithm where an initial solution is found then tour improvement step is done. Under this framework, the heuristic for horse y problem should be three-step algorithms. 1. Find an initial order of customs to be visited, 2. take a tour improvement algorithm for better solution, 3. get a solution for horse y problem out of the order. Step 2 and step 3 can be mingled. The heuristic we used before actually take only the rst and second step. Thus possible improvements can be made based on the current solution. One possible way of tour improving is based on the fact that when the loading locations for UAV is xed, some improvements can be made by swapping orders of some customs. It's just almost the same as the 2-opt, 3-opt and Lin Kernighan heuristic in TSP. And if some improvements are made, the new order can change the optimal loading locations again by solving a convex optimization problem. This procedure can be repeatedly done until no further improvements can be made. How much can it improve the solution we currently have still remains as a question and is worth to research. 9.4 Multiple robots problem In almost all of the previous parts are handling with horse y system with UAVs. But in United States, many regulations on UAVs have limited the usage, for example, there are limit on the height and locations that drone can y. Another possible way to get the benet for horse y problem is to use robots for delivery. Although robots will be able to inhabit on trucks and move independently like UAVs in horse y system, it's obvious that robots may travel slower than trucks. 92 Theoretically, we assume in most of our results that truck travels slower than UAVs, here, this assumption will no longer be held. Some results may be changed under this situation. This will be an interesting problem both as an extension of the current horse y system and as a new multiple robots problem. 9.5 Pareto frontier for multiple objectives We already discussed the horse y problem with two objective functions total cost and time of completion. Some results are obtained. An important result is that when solving optimization problem with linear combination of two objectives for new objective. We got some points and it forms a convex curve as shown in gure 9.1. When solving problems of multiple objectives, a common way is to nd all Pareto frontiers where no other feasible solution will give both better cost and time of completion. It's obvious that all the points on the curve that can be obtained by solving problems with linear combination of objectives as new objective will be Pareto frontiers. However, as the problem is a mixed integer programming problem and thus the space of feasible set is not continuous. There may be some Pareto frontiers as shown in gure 9.1. Either to prove that there is no other Pareto frontiers or to nd them if they exist will be another possible direction. 93 Pareto frontier not obtained Figure 9.1: Points on the curve can be obtained by solving linear combination of two objectives, however, the point shown in the picture is also Pareto frontier but can't be achieved by solving optimization problems. 94 Reference List [1] Concorde TSP Solver. http://www:math:uwaterloo:ca/tsp/concorde:html. [Online; ac- cessed 22-October-2015]. [2] Google Maps Directions API. https://developers:google:com/maps/documentation/ directions/. [Online; accessed 22-October-2015]. [3] U.S. Forest Service UAS policy. http://www:fs:fed:us/science-technology/fire/ unmanned-aircraft-systems. [4] An analysis of several heuristics for the traveling sales- man problem. SIAM Journal on Computing, pages 563{581, 1977. [5] The Google Distance Matrix API. https://developers:google:com/maps/documentation/ distancematrix/intro, 2015. [6] In china, drone delivery promises to boost consumption ... especially in rural areas. APR 2017. https://www:forbes:com/sites/sarahsu/2017/04/13/in-china-drone-delivery- promises-to-boost-consumption-especially-in-rural-areas/#1baf285768fe. [7] Mercedes made a crazy van with built-in drones and robot arms to deliver the packages of to- morrow. Jan 2017. https://qz:com/879605/mercedes-made-a-crazy-van-with-built- in-drones-and-robot-arms-to-deliver-the-packages-of-tomorrow/,. [8] E. Cychosz J. P. Dingens M. Duey K. Koebel S. Siriphanh M. Yurika Tulen H. Watanabe T. Juskauskas J. Holland A. Choi-Fitzpatrick, D. Chavarria and L. Almquist. Up in the air: A global estimate of non-violent drone use 2009-2015. Joan B. Kroc School of Peace Studies at Digital@USanDiego, University of SanDiego, 2016. [9] J.F. Campbell A. Langevin, P. Mbaraga. Continuous approximation models in freight dis- tribution: An overview. Transportation Research, 30B:163{188, 1996. [10] Niels Agatz, Paul Bouman, and Marie Schmidt. Optimization approaches for the traveling salesman problem with drone. ERIM Report Series Reference No. ERS-2015-011-LIS, June 2016. [11] Brigitte Allen. Improving freight eciency within the 'last mile': A case study of wellington's central business district. 2012. [12] D. Applegate, W. Cook, D. S. Johnson, and N. J. A. Sloane. Using large-scale computation to estimate the Beardwood-Halton-Hammersley TSP constant. Presentation at 42 SBPO, 2010. [13] S. Arora. Polynomial time approximation schemes for euclidean traveling salesman and other geometric problems. JACM, 45(5):1998, 1998. 95 [14] Brown G.F. Aykin, T. Interacting new facilities and location-allocation problem. Trans- portation Science, 26:212{222, 1992. [15] Mingozzi A Baldacci R, Christodes N. An exact algorithm for the vehicle routing problem based on the set partitioning formulation with additional cuts. Math Program Ser A, 115:351{ 385, 2008. [16] J. Beardwood, J. H. Halton, and J. M. Hammersley. The shortest path through many points. Mathematical Proceedings of the Cambridge Philosophical Society, 55(4):299{327, 1959. [17] A. Behzad and M. Modarres. New ecient transformation of the gen- eralized traveling salesman problem into traveling salesman problem. Proc. 15th Int. Conf. Syst. Eng., 2002. [18] G. Bensinger. Amazon's drones for deliveries. Wall Street Journal, Dec 2013. http://blogs:wsj:com/digits/2013/12/01/amazons-bezos-shows-flying-drones- for-package-delivery/. [19] E. Brackett. Amazon seeks OK for testing of 50-mph drones. USA Today, Jul 2014. http: //www:usatoday:com/story/money/business/2014/07/10/amazon-drones/12505605/. [20] S Raghavan Majid Salari Bruce Golden, Zahra Naji-Azimi and Paolo Toth. The generalized covering salesman problem. INFORMS Journal on Computing, 24:534{553, 2012. [21] V. Bryan. Drone delivery: DHL 'parcelcopter' ies to German isle. Reuters, Sep 2014. http://www:reuters:com/article/2014/09/24/us-deutsche-post-drones- idUSKCN0HJ1ED20140924. [22] R. A. Zemlin C. E. Miller, A. W. Tucker. Integer programming formulations and traveling salesman problems. J. ACM, pages 326{329, 1960. [23] G. P. Cachon. Retail store density and the cost of greenhouse gas emissions. Management Science, 60(8):1907{1925, 2014. [24] J. F. Campbell, A. T. Ernst, and M. Krishnamoorthy. Hub arc location problems: part I { introduction and results. Management Science, 51(10):1540{1555, 2005. [25] J. F. Campbell and M. E. O'Kelly. Twenty-ve years of hub location research. Transportation Science, 46(2):153{169, 2012. [26] J.F. Campbell. Locating transportation terminals to serve an expanding demand. Trans- portation Research B, 24B:173{192, 1990. [27] J.F. Campbell. Location and allocation for distribution systems with transshipments and transportation economies of scale. Annals of Operations Research, 40:77{99, 1992. [28] Nicos Christodes. Worst-case analysis of a new heuristic for the travelling salesman problem. Graduate School of Industrial Administration, CMU, page 388, 1976. [29] Toth P Christodes N, Mingozzi A. Exact algorithms for the vehicle routing problem based on the spanning tree and shortest path relaxations. Math Program, 20:255282, 1981. [30] Toth P Christodes N, Mingozzi A. State-space relaxation procedures for the computation of bounds to routing problems. Networks, 11:145{164, 1981. [31] Wright JV Clarke G. Scheduling of vehicles from a central depot to a number of delivery points. Oper Res, 12:568{581, 1964. 96 [32] V. Culpan. Watch how Swiss Post is delivering with drones. Wired, Jul 2015. http: //www:wired:co:uk/news/archive/2015-07/09/swiss-delivery-drones. [33] John R. Current and David A. Schilling. The covering salesman problem. Transportation science, 23:208{213, 1989. [34] D. P. Williamson. D. B. Shmoys. Analyzing the held-karp tsp bound: A monotonicity prop- erty wtih application. Information Processing Letters, 35:281{285, 1990. [35] Carlos F. Daganzo. Logistics Systems Analysis. Springer Berlin Heidelberg. [36] Jerey W Heath Damon J Gulczynski and Carter C Price. The close enough traveling salesman problem: A discussion of several heuristics. In Perspectives in Operations Research, pages 271{283, 2006. [37] R. D'Andrea. Guest editorial: Can drones deliver? Automation Science and Engineering, IEEE Transactions on, 11(3):647{648, 2014. [38] G. Dantzig and J. Ramser. The truck dispatching problem. Management Science, 6:80{91, 1959. [39] V. Dimitrijevic and Z. Saric. Ecient transformation of the generalized traveling salesman problem into the traveling salesman problem on digraphs. Information Sci, 102:65{110, 1997. [40] Michael Drexl. Branch-and-price and heuristic column generation for the generalized truck- and-trailer routing problem. Revista de Metodos Cuantitativos para la Economia y la Em- presa, pages 5{38, 2011. [41] Michael Drexl. Branch-and-cut algorithms for the vehicle routing problem with trailers and transshipments. Networks, 63:119{133, 2014. [42] L. Evers, T. Dollevoet, A. I. Barros, and H. Monsuur. Robust UAV mission planning. Annals of Operations Research, 222(1):293{315, 2014. [43] L. Few. The shortest path and the shortest road through n points. Mathematika, 2:141{144, 1955. [44] Lysgaard J Poggi de Aragao M Reis M Uchoa E Werneck RF Fukasawa R, Longo H. Robust branch-and-cut-and-price for the capacitated vehicle routing problem. Math Program Ser. A, 106:491{511, 2006. [45] Y. Nobert G. Laporte and S. Taillefer. Solving a family of multi-depot vehicle routing and location-routing problems. Transp. Sci, 22:161{172, 1988. [46] M. Gendreau, A. Hertz, and G. Laporte. A tabu search heuristic for the vehicle routing problem. Management Science, 40(10):1276{1290, 1994. [47] M. Gendreau, G. Laporte, C. Musaraganyi, and E. D. Taillard. A tabu search heuristic for the heterogeneous eet vehicle routing problem. Computers & Operations Research, 26(12):1153{ 1173, 1999. [48] Laporte G Gendreau M, Hertz A. A tabu search heuristic for the vehicle routing problem. Manag Sci, 40:1276{1290, 1994. [49] M. Goemans. Worst-case comparison of valid inequalities for the tsp. Mathematical Pro- gramming, 69:335{349, 1995. 97 [50] A.J. Goldman. Optimal location for centers in a network. Transportation Science, 3:352{360, 1969. [51] M. Haimovich and Thomas L. Magnanti. Extremum properties of hexagonal partitioning and the uniform distribution in Euclidean location. SIAM J. Discrete Math., 1:50{64, 1988. [52] B. Hajek. Random Processes for Engineers. Cambridge University Press, 2015. [53] K. Helsgaun. An eective implementation of the lin-kernighan traveling salesman heuristic. European Journal of Operational Research, (1):106{130, 2000. [54] K. Helsgaun. General k-opt submoves for the lin-kernighan tsp heuristic. Mathematical Programming Computation, 2009. [55] K. Helsgaun. Solving the equality generalized traveling salesman problem using the link- ernighanhelsgaun algorithm. Mathematical Programming Computation, 7(3):269{287, 2015. [56] M. Heutger and M. K uckelhaus. Unmanned aerial vehicle in logistics: a DHL perspective on implications and use cases for the logistics industry. Technical report, DHL Trend Research, 2014. http://www:dhl:com/content/dam/downloads/g0/about us/logistics insights/ DHL TrendReport UAV:pdf. [57] D. S. Hochbaum. When are NP-hard location problems easy? Annals of Operations Research, 1:201{214, 1984. [58] M. Kanellos. Google working on drones too. Forbes, Aug 2014. http://www:forbes:com/ sites/michaelkanellos/2014/08/29/google-working-on-drones-too/. [59] G. Karapetyan. D., Gutin. Ecient local search algorithms for known and new neighborhoods for the generalized traveling salesman problem. Eur. J. Oper. Res, 219(2):234{251, 2012. [60] Gelatt Jr. C.D. Kirkpatrick, S. and M.P Vecchi. Optimization by simulated annealing. Sci- ence, pages 671{680, 1983. [61] S.G. Krantz and H.R. Parks. Geometric Integration Theory. Cornerstones Series. Birkh auser, 2008. [62] M. Kress and J. O. Royset. Aerial search optimization model (ASOM) for UAVs in special operations. Military Operations Research, 13(1):23{33, 2008. [63] Semet F Laporte, G. Computational evaluation of a transformation procedure for the sym- metric generalized traveling salesman problem. INFOR, 37(2):114{120, 1999. [64] Nobert Y Laporte G, Desrochers M. Two exact algorithms for the distance-constrained vehicle routing problem. Networks, 14:161{172, 1984. [65] J.; Rinnooy Kan A.; Lawler, E.; Lenstra and Shmoy. D. The Traveling Salesman Problem: A Guided Tour of Combinatorial Optimization. New York: Wiley, 1985. [66] B. W. Lin, Shen; Kernighan. An eective heuristic algorithm for the traveling-salesman problem. Operations Research, (2):498{516, 1973. [67] C. K. Y. Lin. A vehicle routing problem with pickup and delivery time windows, and co- ordination of transportable resources. Computers & Operations Research, 38(11):1596{1609, 2011. [68] J. J. S. Gonz alez M. Fischetti and P. Toth. A branch-and-cut algorithm for the symmetric generalized traveling salesman problem. Operations Research, 45(3):378{394, 1997. 98 [69] A. Madrigal. Inside the drone missions to Fukushima. The Atlantic, Apr 2011. http://www:theatlantic:com/technology/archive/2011/04/inside-the-drone- missions-to-fukushima/237981/. [70] C. McCall. Sea Shepherd says drones nd, photograph Japan's whaling eet. Reuters, Dec 2011. http://www:reuters:com/article/2011/12/25/us-australia-japan- whaling-idUSTRE7BO01K20111225. [71] C. C. Murray and A. G. Chu. The ying sidekick traveling salesman problem: Optimization of drone-assisted parcel delivery. Transportation Research Part C: Emerging Technologies, 54:86{109, 2015. [72] Love R. Norback, J. Geometric approaches to solving the traveling salesman problem. Man- agement Science, pages 1208{1223, 1977. [73] C. E. Kim O. H. Ibarra. Fast approximation algorithms for the knapsack and sum of subsets problems. JACM, 22(4):463{468, 1975. [74] M.E. O'Kelly. The location of interacting hub facilities. Transportation Science, 20:92{105, 1986a. [75] M.E. O'Kelly. A quadratic integer program for the location of interacting hub facilities. European Journal of Operational Research, 32:393{404, 1987. [76] C. H. Papadimitriou. Euclidean tsp is np-complete. Theoretical Computer Science, 4:237{ 244, 1977. [77] C. H. Papadimitriou. Worst-case and probabilistic analysis of a geometric location problem. SIAM Journal on Computing, 10:542, 1981. [78] Ropke S Pisinger D. A general heuristic for vehicle routing problems. Comput Oper Res, 34:2403{2435, 2007. [79] Xingyin Wang; Stefan Poikonen and Bruce Golden. The vehicle routing problem with drones: several worst-case results. Optimization Letters, pages 1{19, 2016. [80] A. Ponza. Optimization of drone-assisted parcel delivery. [81] Quang Dung Pham Minh Hoang Ha Quang Minh Ha, Yves Deville. On the min-cost traveling salesman problem with drone. ICTEAM/INGI/EPL, 2016. Technical Report. [82] C. Redmond and J. E. Yukich. Limit theorems and rates of convergence for euclidean func- tionals. The Annals of Applied Probability, 4(4):pp. 1057{1073, 1994. [83] Jean-Paul; Claude Comtois; Brian Slack Rodrigue. The "Last Mile" in Freight Distribution. Routledge, 2009. [84] R. Motwani M. Sudan S. Arora, C. Lund and M. Szegedy. Proof verication and intractability of approximation problems. In Proc. 33rd IEEE Symp. on Foundations of Computer Science, pages 13{22, 1992. [85] R. C. Garg S. S. Srivastava, S. Kumar and P. Sen. Generalized traveling salesman problem through n sets of nodes. CORS journal, 7:97{101, 1969. [86] Stephan Scheuerer. A tabu search heuristic for the truck and trailer routing problem. Com- puters & Operations Research, 33:894{909, 2006. 99 [87] Martia Scott. Improving freight movement in delaware central business districts. Institute for Public Administration, University of Delaware, November 2009. [88] Scott C Scott, J. Drone delivery models for healthcare. Proceedings of the 50th Hawaii International Conference on System Sciences, 2017. [89] F Yu Vincent Shih-Wei Lin and Shuo-Yan Chou. Solving the truck and trailer routing prob- lem based on a simulated annealing heuristic. Computers & Operations Research, 36:1683{ 1692, 2009. [90] R. Smith. Utilities turn to drones to inspect power lines and pipelines. Wall Street Jour- nal, May 2015. http://www:wsj:com/articles/utilities-turn-to-drones-to-inspect- power-lines-and-pipelines-1430881491. [91] J. M. Steele. Subadditive Euclidean functionals and nonlinear growth in geometric probabil- ity. The Annals of Probability, 9(3):pp. 365{376, 1981. [92] J.M. Steele. Probability Theory and Combinatorial Optimization. CBMS-NSF Regional Con- ference Series in Applied Mathematics. Society for Industrial and Applied Mathematics, 1987. [93] P. Toth and D. Vigo. The vehicle routing problem. Society for Industrial and Applied Mathematics, 2001. [94] Gendreau M Prins C Vidal T, Crainic TG. A unied solution framework for multi-attribute vehicle routing problems. Tech Rep, CIRRELT, 2012. [95] G. S. Lueker W. Fernandez de la Vega. Bin packing can be solved within 1 + in linear time. Combinatorica, 1(4):349{355, 1981. [96] M. Wohlsen. The next big thing you missed: Amazon's delivery drones could work { they just need trucks. Wired, Jun 2014. http://www:wired:com/2014/06/the-next-big-thing- you-missed-delivery-drones-launched-from-trucks-are-the-future-of-shipping/. [97] E. Wygonik and A. Goodchild. Evaluating the ecacy of shared-use vehicles for reducing greenhouse gas emissions: a US case study of grocery delivery. In Journal of the Transporta- tion Research Forum, volume 51, 2012. [98] W. Xie and Y. Ouyang. Optimal layout of transshipment facility locations on an innite homogeneous plane. Transportation Research Part B: Methodological, 75:74{88, 2015. [99] E. Ma Y. Lien and B. W. Wah. Transformation of the generalized traveling salesman problem into the standard traveling salesman problem. Information Sci., 74:177{189, 1993. 100 Appendix A Proof of Lemma 2 The assumption thatL is polygonal is a valid one because the standard denition of the length of L is simply the limit of a discretization; if we say thatL is the image of the continuous mapping : [0; 1]!R 2 , then the length ofL is dened as length(L) = sup 0=t0<t1<<t M =1 M1 X i=1 k (t i ) (t i1 )k where the supremum is taken over all possible partitions of [0; 1] and M is unbounded. The assumption thatL forms the boundary of a convex region is also straightforward: at the very least, we can certainly assume thatL is simple, that is, thatL does not intersect itself (since one can always \un-cross" a pair of intersecting edges ofL in an obvious way). IfL is not the boundary of a convex region, then (since we have assumed thatL is polygonal) there must exist a pair of vertices v 1 ;v 2 ofL such that v 1 and v 2 are adjacent on the convex hull Conv(L) ofL, but not adjacent onL itself; see Figure A.1a. If we let s denote the component ofL that lies between v 1 and v 2 , then it is obvious that we can re ect s about the line joining v 1 and v 2 , to obtain a new curveL 0 with the same length asL, as shown in Figure A.1b. It is then an entirely straightforward argument to verify that, for any , we have Area(N (L 0 )) Area(N (L)), which completes the proof. v 1 v 2 s (a) v 1 v 2 s (b) Figure A.1: Re ecting the polygonal chain s about the line joining v 1 and v 2 . 101

Abstract (if available)

Abstract Unmanned aerial vehicle (UAV) has seen its development recently in lots of areas, both in military and commercial areas. This paper is on its application of package delivery where UAVs, provides service to customers while making return trips to a truck that is itself moving. In other words, a UAV picks up a package from the truck (which continues on its route), and after delivering the package, the UAV returns to the truck to pick up the next package. In this paper, we first formulate this problem called horsefly routing problem and then give both asymptotic analysis and heuristic algorithms. ❧ Then some variants of it are considered. Those problems are all considering some extra factors under some circumstances in practical uses. Chapter 4 considers the problem with more than one UAV in the system and chapter 5 works on a problem with both cost and time of completion as objectives. Chapter 6 considers the fact that nowadays UAV only has limited battery and handles with the horsefly routing problem with battery life constraints. Chapter 7 is on problem with region restriction on truck, it happens when taking the advantage of UAV that reach to some rural areas without roads. Finally, chapter 8 solves a problem that launching sites are restricted to some pre-defined locations.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Dynamic programming-based algorithms and heuristics for routing problems

PDF

Routing problems for fuel efficient vehicles

PDF

The warehouse traveling salesman problem and its application

PDF

Train routing and timetabling algorithms for general networks

PDF

A continuous approximation model for the parallel drone scheduling traveling salesman problem

PDF

Train scheduling and routing under dynamic headway control

PDF

Models and algorithms for pricing and routing in ride-sharing

PDF

Routing for ridesharing

PDF

Asymptotic analysis of the generalized traveling salesman problem and its application

PDF

Models and algorithms for the freight movement problem in drayage operations

PDF

Applications of explicit enumeration schemes in combinatorial optimization

PDF

Cost-sharing mechanism design for freight consolidation

PDF

Applications of Wasserstein distance in distributionally robust optimization

PDF

On the interplay between stochastic programming, non-parametric statistics, and nonconvex optimization

PDF

Continuous approximation formulas for cumulative routing optimization problems

PDF

New approaches for routing courier delivery services

PDF

Stochastic games with expected-value constraints

PDF

Continuous approximation formulas for location and hybrid routing/location problems

PDF

Continuous approximation for selection routing problems

PDF

Computational geometric partitioning for vehicle routing

Asset Metadata

Creator Song, Siyuan (author)

Core Title Package delivery with trucks and UAVs

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Industrial and Systems Engineering

Publication Date 06/27/2018

Defense Date 04/23/2018

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag logistics,OAI-PMH Harvest,package delivery,traveling salesman problem,UAV,vehicle routing

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Carlsson, John Gunnar (committee chair), Savla, Ketan (committee member), Vayanos, Phebe (committee member)

Creator Email siyuanso@usc.edu,ssydyc@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-505659

Unique identifier UC11268054

Identifier etd-SongSiyuan-6351.pdf (filename),usctheses-c40-505659 (legacy record id)

Legacy Identifier etd-SongSiyuan-6351.pdf

Dmrecord 505659

Document Type Dissertation

Format application/pdf (imt)

Rights Song, Siyuan

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA