Motion- and Communication-Planning of Unmanned Aerial Vehicles in Delay Tolerant Network using Mixed-Integer Linear Programming

Large amounts of data are typically generated in applications such as surveillance of power lines and railways, inspection of gas pipes, and security surveillance. In the latter application it is a necessity that the data is transmitted to the control centre “on-the-fly” for analysis. Also missions related to other applications would greatly benefit from near real-time analysis and operator interaction based on captured data. This is the motivation behind this paper on coarse offline motionand communication-planning for cooperating Unmanned Aerial Vehicles (UAVs). A Mixed-Integer Linear Programming (MILP) problem is defined in order to solve the surveillance mission. To efficiently transmit the data back to the base station the vehicles are allowed to store data for later transmission and transmit via other vehicles, in addition to direct transmission. The paths obtained by solving the optimization problem are analyzed using a realistic radio propagation path loss simulator. If the radio propagation path loss exceeds the maximum design criterion the optimization problem is solved again with a stricter communication constraint, and the procedure is continued in an iterative manner until the criterion is met. The proposed algorithm is supported by simulations showing the resulting paths and communication topologies for different choices of delay tolerance.


Introduction
Among the first to consider multi-vehicle path planning using MILP, were Schouwenaars et al. (2001), Earl and D'Andrea (2002) and Richards and How (2002).Trajectory generation for autonomous vehicles was considered in Schouwenaars et al. (2001), while a robotic ball game was the application of Earl and D'Andrea (2002).Since then MILP has been used extensively for path planning problems, in particular for UAVs, both for single and multi-vehicle systems, see Richards and How (2002), Ma and Miller (2006), Shengxiang and Hailong (2008), Kuwata and How (2011), Grøtli and Johansen (2012b) and Grøtli and Johansen (2012c).In Schouwenaars et al. (2006) connectivity constrained trajectory planning for autonomous helicopters through cluttered environments was studied.Line-of-sight connectivity between a leader helicopter and base station was maintained by coordinating the motion of intermediate helicopters.Other motion planning techniques for improving network or communication properties, which do not necessarily involve solving an optimization prob-lem, have been treated in Spanos and Murray (2005), Dixon and Frew (2005) and Dixon and Frew (2007) among many others.
In on of our previous papers, Grøtli and Johansen (2012b), we described how paths for multiple UAVs can be planned to create a communication chain between two base stations, also known as chaining, Dixon and Frew (2007).In Grøtli and Johansen (2012c) paths were planned for a multi-task mission similar to that presented in this paper, and where the mission objective only could be met if one of the vehicles were used as a relay node.By relaying, we mean that additional nodes are used to receive a transmission from a source and retransmit it to a destination, Dixon and Frew (2007), Frew and Brown (2008).Task assignment for multiple UAVs has been considered in several earlier papers, see for instance Kim et al. (2007), Alighanbari et al. (2003) and Kingston and Schumacher (2005), but the focus of these papers has mainly been on scheduling or selection of straight line paths.
The results from Grøtli and Johansen (2012c) are generalized in this paper in order to allow also for ferrying in a Delay Tolerant Network (DTN).The first notion of message ferrying was developed by Zhao and Ammar (2003).Ferrying means that a mobile node physically stores and carries data from one location to another Frew et al. (2006), Frew and Brown (2008).This has the advantage of extended range and possibly increased total data transmission rate compared to direct communication, see Figure 1.Henkel and Brown (2006) studied route design for aerial data ferrying nodes.The route designs were evaluated analytically for relative comparison based on node velocity, data rate, and buffer size.In this paper ferrying of data is incorporated into the optimization problem.This makes the optimization problem much more complex, as it is necessary for instance to keep track of which node stores the data, how much data is stored on each node and how long data has been in the network.The reason for this is physical constraints on the storage capacity of each node and desired limitation on the time taken from when the data is collected by the UAV until it is received at the base station.On the upside, ferrying means increased flexibility with respect to coordination between vehicles compared to relaying, and it also means that the surveillance range of a mission can be extended beyond the line-of-sight distance of a chain of relaying vehicles.
As in Grøtli and Johansen (2012b) and Grøtli and Johansen (2012c) we analyze the planned paths in a radio propagation simulator to get a more accurate prediction of the radio path loss between the vehicles, and between the vehicles and the base station.In the MILP formulation the ability to communicate at a certain data rate depends on the distance between the nodes.For this distance to give a realistic picture of the communication properties, the paths found by solving the MILP problem are analyzed at every time step using the radio path loss simulator SPLAT!Maglicane (2010(Accessed August 18, 2010), Wright (2011(Accessed June 29, 2011).SPLAT! uses digital elevation data to calculate field strength and path loss based on the Longley-Rice Irregular Terrain Model Longley and Rice (1968).If the path loss estimate calculated by SPLAT! is too high to maintain communication at the desired rate, the communication constraints of the MILP problem are tightened.This means that the maximal distance where communication is assumed feasible is reduced by a certain value.The process is then repeated until paths are found in which communication can be maintained at a predefined criterion during the servicing of the tasks.We emphasize that we have proposed an offline algorithm, which provides coarse motion-and communication-planning for the UAVs.In Beard and McLain (2012) they discriminate between two different approaches to motion planning: "deliberative motion planning, where explicit paths and trajectories are computed based on global world knowledge, and reactive motion planning, which uses behavioral methods to react to local sensor information".Our method falls in the first category as it is designed for preplanning of trajectories, which then can serve as an input towards reactive planning methods.For instance the online re-planning methods described in Grancharova et al. (2012Grancharova et al. ( , 2014) ) are well suited for improving robustness towards inaccuracies in the model and uncertainties that become apparent in real-time.Many other methods have been proposed in order to solve the motion planning problem.For an extensive overview over path-and motion-planning techniques the interested reader is referred to LaValle (2006) or Tsourdos et al. (2010).In our opinion MILP is well suited for complex missions as outlined here.Particular benefits of MILP are that logics, approximations of nonlinear functions, and non-convexity can easily be handled, Richards and How (2005); Bemporad and Morari (1999); Williams (1999).Another important important property is the ability to define hard constraints.Other path-and motion-planning methods do not easily handle hard inequality constraints.Instead, dissatisfied constraints are penalized in the cost function.This makes these methods hard to tune, because of the many (often competing) objectives that must be weighted appropriately.Also, many optimization methods are not easily modified to allow for coordination between multiple vehicles.Examples of popular methods used for path-or motion planning include rapidly exploring random trees, LaValle and Kuffner (2001); Kuwata et al. (2009), particle swarm optimization Kennedy and Eberhart (1995); Ho et al. (2013); Saska et al. (2006), A andD Hart et al. (1968); Stentz (1994); Likachev et al. (2005) or variants of these.The performance of these methods relies heavily on the choice of a good heuristic potential.Also the application to the specific mission in this paper is challenging because of the above-mentioned reasons.Their primary strengths are within real-time local path planning, a category where many methods are shown to outperform MILP, Thunberg et al. (2008).Finally, the fact that the MILP provides a globally optimal solution makes it useful for benchmarking other optimization methods, i.e how far away is the solution of other sub-optimal methods.

Contribution
The main benefit of the proposed approach compared to earlier works, lies in the flexibility obtained by allowing for direct transmission, relaying and ferrying of data in the same framework.To our knowledge, we present a novel approach to an optimization problem by combining motion-and data-transmission planning while incorporating important constraints on vehicle motion (visiting of waypoints, vehicle speed, height above ground) and communication parameters (bandwidth, buffer size, path loss, delay, etc.).Although there are several authors who consider a combined path-and communication-planning problem, the ability to communicate is typically approximated as a function of the relative distance between vehicles.Even though this is also true for the MILP problem proposed in this paper, we emphasize that all paths will eventually have to satisfy a radio propagation path loss requirement in the simulator SPLAT! as described above.
In Grøtli and Johansen (2012b) we dealt with motion planning in order to create a relay chain between two stationary ground-or surface nodes.In comparison, in this paper we have presented a solution to a surveillance mission planning problem.In order to achieve the main objective of surveillance of all tasks as efficiently as possible, the means of communication (direct, relaying, ferrying) during the mission is chosen based on important communication parameters such as buffer size, delay, and bandwidth.We point out that this paper is an extended version of our conference paper, Grøtli and Johansen (2012a).The latter paper lacked a detailed description of many of the constraints of the MILP optimization problem due to space constraints.These are now included.In addition, we have extended the simulations section to cover both when the UAVs have the possibility to store data for later transmission, and when the transmission has to be immediate.working can be used not only to extend the surveillance range, but also to increase the total data transmission rate.In the dark red areas the UAV can communicate with the base stations at a high transmission rate, whereas the transmission rate is poor in the light red areas.By bringing the UAV into range of the dark red areas, the total data transmission rate can be substantially increased even when taking into account the UAV's flight time.To quote Tanenbaum (2003): "Never underestimate the bandwidth of a station wagon full of tapes hurtling down the highway".The process of storing data on a mobile node and physically bringing it from one location to another is called ferrying.

MILP Problem formulation
The mission objective is to perform data acquisition along sequences of waypoints by UAVs, and real-time transmission of sensor data back to the base station.The waypoints are not necessarily within communication range for direct transmission, and we will therefore allow for one or more UAVs to function as relay links.We will assume that the tasks consist of segments (e.g.road, railway or power line segments), which can be described by one or more waypoints.The interest in real-time transmission of sensor data while observing these segments, is to allow for a human operator to intervene if something irregular is found in the sensor data.

Notation
Abbreviations: • MILP -Mixed Integer Linear Programming • ECEF -Earth Centered, Earth Fixed • ENU -East-North-Up • TIN -Triangulated Irregular Network Constants: • n p -total number of vehicles refer, depending on the situation, either to the physical vehicles or the communication nodes onboard these vehicles.In addition we consider one communication node at the base station, which is given the index {n p + 1}, such that the index set of all communication nodes is P • p pi := (p 1pi , p 2pi , p 3pi ) := (x pi , y pi , z pi ) -position vector of vehicle p and time step i along the axes of a local East-North-Up (ENU) coordinate reference frame, see Figure 3 • v pi := (v 1pi , v 2pi , v 3pi ) -velocity of vehicle p at time step i along the axes of a local ENU coordinate reference frame and vehicle q ∈ P np p , in at least one of the direction of the ENU frame is greater than the minimum separation distance • w acc pi := (w acc 1pi , w acc 2pi , w acc 3pi ) -vector of slack variables, used to penalize the acceleration in each direction of the ENU frame • χ pqi := (x pi − x qi , y pi − y qi , z pi − z qi ) -vector of relative distance between node p ∈ P np+1 1 and node q ∈ P np+1 1 along each of the ENU coordinate directions at time step i ∈ I

Vehicle model
In the planning problem we assume that the p th UAV is described by the discrete time model , where n p is the number of UAVs, ∆t is the sample time, and p pi and v pi are vectors with elements being positions and velocities along the orthogonal axes of a local East-North-Up (ENU) coordinate reference frame.The simplicity of the model makes it suitable for the optimization problem described in this paper, where we focus on coarse path planning.The relationship between the ENU and the Earth-Centered Earth-Fixed (ECEF) frames is explained in Figure 3.The ECEF frame, is a coordinate system where its origin is at the center of Earth, and with axes X, Y , Z, rotating with the angular velocity of Earth.The ENU coordinate frame is a local geodetic coordinate system whose tangent plane is fitted to the geodetic reference ellipsoid at some convenient point for local measurements.The x axis points towards East, the y axis points towards North and the z axis completes the right-handed orthogonal frame by pointing away from the Earth perpendicular to the reference ellipsoid.The origin of the ENU frame is typically represented in terms of the reference ellipsoidal parameters longitude l, and geodetic latitude µ.

Velocity constraints
We approximate V pi in a similar manner as in Chaudhry et al. (2004), here in the three-dimensional case as in Grøtli and Johansen (2012b), by introducing the constraints: , where b vel pikl are binary optimization variables, the unit vector with θ k := 2πk/D vel , φ l := 2πl/D vel , k ∈ {1, . . ., D vel /2}, l ∈ {1, . . ., D vel } and the discretization level D vel is some constant even integer greater or equal to four.By (2), V pi is constrained to be larger than the scalar projection of v pi onto any of the unit vectors ξ kl .To illustrate this constraint, consider the simpler constraint v pi ξ kl ≤ 1 in the horizontal plane in Figure 2, for some particular p ∈ P np 1 , some particular i ∈ I N −1 0 , ∀k ∈ {1, . . ., D vel = 6}, and l = 1.The dashed arrows illustrate the unit vectors, ξ k1 , k ∈ {1, . . ., D vel = 6}.These constraints ensure that the projection of the velocity v pi onto the unit vectors ξ k1 (the scalar product of v pi and ξ k1 ) are less or equal to 1. Hence, the velocity vector should be within the red area in the horizontal plane, which is a polygonal approximation to the unit circle.The constant D vel is required to be an even integer greater or equal to four in order to keep symmetry about both axes in the horizontal plane.In the constraints in (2) the right-hand side is substituted with the speed V pi , which in turn is constrained by V pi ≤ V p , a constant maximum allowed speed for vehicle p.
In (3) V p is required to be less than α vel v pi ξ kl , or the corresponding binary variable b vel pikl must be zero.To prevent (3) to be trivially satisfied with all b vel pikl being zero, we add the additional constraints (4).The accuracy of the approximation depends of course on the discretization level D vel (see Figure 4), but also on α vel , a constant slightly greater than one.The closer to one α vel is, the better is the approximation; however, taking it too close may have a negative impact on the computation time of the MILP problem, Chaudhry et al. (2004).The constants M vel pkl should be chosen sufficiently large.
Remark 1 Here, and in the rest of the document, we mean by sufficiently large (or sufficiently small) in this context that the constant should be chosen large (small) enough to maintain the original logical implication the constraint is meant to realize.Consider for instance the constraint f (x) ≤ 0, where f : R n → R is linear, and x ∈ X for a given bounded set X .Then a sufficiently large (small) constant can be chosen as Bemporad and Morari (1999).Although M (m) in theory could be taken to be arbitrarily large (small), this is not recommended for computational efficiency, Williams (1999).
In YALMIP Löfberg (2008), a MATLAB toolbox for implementing optimization problems, logic implications can be expressed instead of big-M formulations such as (3), and YALMIP will automatically derive big-M coefficients by analyzing the constraints on the other variables involved in the expression.As the speed of the vehicles is approximated by (2), ( 3) and (4), we simply use that , where V p and V p are the minimum and maximum velocities, respectively, of vehicle p.If the binary variable b wp piw is true, this means that vehicle p is visiting waypoint w at time step i.More on this implication and visiting of waypoints is postponed and considered in Section 3.8.Equation (6) will constrain the velocity of the vehicles to zero when they have arrived at the waypoint w ∈ W nt , which is the final waypoint of the mission.

Acceleration cost
To avoid fluctuations in the speed, we introduce the following cost function similar to the one proposed in Schouwenaars et al. (2001), with the additional constraints where w acc pi := (w acc 1pi , w acc 2pi , w acc 3pi ) and r p ∈ R 3 ≥0 is a nonnegative weighting vector.The motivation behind ( 7) is to penalize the absolute value of acceleration in each direction of the ENU frame.Also, to avoid a piecewise linear cost function, we have introduced slack variables w acc jpi .

Position constraints
There are typically restrictions on where UAVs are allowed to fly.This may for instance be air space used for other air traffic, air space over a certain altitude, or air space over populated areas.In addition the operator might want to avoid flying into regions with bad weather, outside the area where the operator is able to communicate with UAVs or in case of military applications: areas with enemies and enemy radars.If the region the UAVs are required to stay within is convex (e.g. a rectangular box), the constraints may simply be written , where x, y, z and x, y, z are the constant lower and upper bounds, respectively, on the state vector in the east, north and up directions.More generally, unions of convex sets can be implemented.

Anti-collision constraints
To avoid collision between vehicles we will implement the method of Schouwenaars et al. (2001).Let the position of vehicle p and vehicle q at time step i be given by (x pi , y pi , z pi ) and (x qi , y qi , z qi ), respectively.The constraints on their relative position are then given as , where d x , d y and d z are the safety distances in the east, north and up directions, respectively.These safety distances represents the separation required to still maintain the ability to perform avoidance maneuvers.The binary variables b col pqil , ensure that there is a minimum separation distance between the vehicles in at least one of the directions of the ENU frame.The constant M col 1 should be taken sufficiently large, see Remark 1, for instance M col 1 > x − x + d x with x, x as in ( 10), and correspondingly for M col 2 and M col 3 .

Anti-grounding constraints
As in Ma and Miller (2006) we will represent the terrain as a triangulated irregular network (TIN).Terrain avoidance constraints in MILP form are given in Shengxiang and Hailong (2008), and will be used here.T TIN non-overlapping triangles with m TIN vertices ) are used to represent the piecewise affine terrain surface.The point strictly below vehicle p at time step i is given by (x pi , y pi , h pi ), and satisfy is the set of indices of triangles that have a common vertex P l .An example of the P 1 enumeration of triangles, vertices and the corresponding sets D l is shown in Figure 5. Equations ( 17)-( 21) describe the position strictly below the vehicle in terms of its barycentric coordinates, as illustrated in Figure 6.By ( 22), the binary variables b TIN pit are forced to be true if and only if vehicle p flies over triangle t at time step i.By ( 23), the variables λ pil corresponding to vertices which are not adjacent to the particular triangle are set to zero.Finally, the terrain avoidance constraint can be expressed as , where d TIN is the minimum vertical distance from the UAVs to the ground.The TIN is generated from the elevation data by incremental Delaunay triangulation.This reduces the complexity of the problem, since only a subset of the available data is used in the MILP formulation.

Task assignment
We will assume that there are n t tasks, and that each task t ∈ T nt 1 is comprised of a set of waypoint indices.Let W t denote the set of the indices of the waypoints which belong to task t ∈ T nt 1 .A special meaning is given to the final task W nt .It contains only one waypoint, which is located above the designated landing site for the UAVs.Take-off and landing is not considered in this paper, but assumed to be handled separately.We require that W r W s = ∅, ∀r, s ∈ T nt 1 , r = s, i.e. that each waypoint belongs to one and only one task.A waypoint characterized by the ENU coordinates (x wp 1w , y wp 2w , z wp 3w ) is considered to be visited if a UAV is flying through a cube containing the waypoint.More precisely, we assume each waypoint to be a cube with sides of length 2d wp , and require that Richards x C be the position vectors of the blue, red and green vertices of the triangle, respectively.The position of any point x P on the triangle is given by A , and λ C = A C A .A A , A B , A C are the areas of the blue, red and green subtriangles, respectively, and A is the total area of the triangle such that and How (2002) , where M wp pw1 . . .M wp pw6 are chosen sufficiently large.This way, the binary variable b wp piw is true then vehicle p flies through waypoint w at time step i.We will require that each waypoint of tasks W t , ∀t ∈ T nt−1 1 , is visited once and once only.Mathematically this is formulated by the equality constraint By assigning all vehicles the final task, all vehicles will return to the landing site before the end of the mission.This assignment is ensured by the constraints , w ∈ W nt , which means that once b wp piw has become true for some time step i, it will remain true for the rest of the horizon.Together with the implication presented in ( 6) and ( 25), this also means that vehicle p will remain at the final waypoint (landing site) once it has arrived there.Given that b piw is true, the implication in (25) constrains the position of vehicle p to the final waypoint, whereas Equation ( 6) constrains the velocity to zero.These constraints indirectly force all tasks T nt−1 1 to be executed before the the last vehicle returns to the landing site.Since each waypoint w ∈ W 1 ∪ . . .∪ W nt−1 , is visited only once, the time steps elapsed before a waypoint is visited are given by , where θ w is a variable we have introduced in our optimization problem.Furthermore, we require that the waypoints within the same task are visited in a specific order.Let T + represent those tasks with more than one waypoint.Then, W t \W Last t , t ∈ T + contains the indices of all the waypoints of task t, except the last one.Then the visiting order can be achieved by requiring that ∀w ∈ W t \W Last t , t ∈ T + .The number of time steps elapsed before vehicle p returns to the landing site is given by θ finish p , if we use the constraints ∀p ∈ P np 1 , i ∈ I N 0 , w ∈ W nt where M finish is a constant chosen sufficiently large, see Remark 1, for instance as M finish := N .Recall that for each vehicle p, b wp piw may be true for many consecutive time steps i, so we cannot use the same approach as in (34) to find the time elapsed before vehicle p arrives at waypoint w ∈ W nt .Instead with the upper-and lower-bounds on θ finish p given by ( 36) and (37), respectively, θ finish p will represent the exact time step of arrival at the final waypoint.Since we want to minimize the overall mission timethe time elapsed until the last vehicle arrives at the final waypoint -we introduce the variable η finish and require that ∀p ∈ P np 1 , and set our objective to minimize the cost function where γ finish is a positive scalar.Equation ( 38) can be satisfied for any sufficiently large η finish .However, by minimizing η finish in the cost function, we achieve the desired effect, which is to minimize the overall mission time.We do not want vehicles to arrive at the final waypoint simultaneously, as this may cause the UAVs to collide.Therefore, we also require a temporal separation between the arrival at the final waypoint, that is, , where t separation ∈ N is the number of time steps separating the UAVs at the arrival of the final waypoint.Thus far, there is nothing restricting multiple vehicles each accomplishing parts of a task.As this may be undesirable, for instance because we want video to be recorded continuously between waypoints of a task, we introduce an additional binary variable b task pt which is true if and only if task t is served by vehicle p.This is achieved by imposing the constraints , where n Wt is the number of waypoints of task t.Still, there is a possibility that a UAV will switch back and forth between different tasks.This behavior is allowed, but as it is shown in the simulations, it is more beneficial to accomplish one task at the time.This is due to the demanding communication constraints during the accomplishment of a task, which we will impose in Section 3.10.

Data gathering
We introduce the binary variable b sensor pi which is true if and only if vehicle p is serving a task at time step i.This implication can be achieved by using the auxiliary variable λ pit ∈ [0, 1] and imposing the constraints , and finally and W Last t represent the set with index to the first and the last waypoint, respectively, of task t.The inequalities in (43) force λ pit to be less or equal to one for every time step i during the servicing of task t, and zero otherwise.Equation ( 44) makes i∈I N 1 λ pit , for a specific vehicle p and a specific task t, to be equal to the number of time steps elapsed from when the vehicle visited the first waypoint of the task and until it visited the last (and where we have assumed that the servicing of the task ends at the end of the step).Together, ( 43) and ( 44) constraint λ pit to be one at the time steps i at which vehicle p is servicing task t and zero otherwise.Finally, (45) will give the variables b sensor pi the desired property.

Data flow for delayed transmission
In the following we will sometimes commonly refer to the base station (where the antenna for the operator user interface is located) and the UAVs as nodes.The communication network will be modeled as a buffered flow network, a directed graph where each edge represents a limited transmission capacity, but where the nodes have the ability to store data.We assume that while vehicle p is servicing a task -that is for those time steps i the binary variable b sensor pi is true -the rate at which sensor data is gathered for later (or immediate) transmission back to base station, is given by c sensor .In our setup, the base station is the only sink, whereas the UAVs act as sources during servicing of a task.We assume that the bandwidth required for transmission is substantially larger during the execution of a task, and ignore the possible need for communication during transit between tasks.Immediate transmission of data back to the base station requires that the vehicles and therefore also the tasks lie within direct communication distance of the base station, or at least the communication distance of a chain of multiple UAVs.As this limits the surveillance area, we also allow for ferrying.This may increase the surveillance area.It could potentially also increase the data transfer rate, as the vehicles can ferry the data into a range where high bandwidth transmission is possible.The drawback is added complexity, and of course the additional delay between sensing and receiving of data.In our task assignment scenario, we require that the collected data is forwarded within t delay time steps.This allows a human operator to analyze the data, and possibly command the UAVs to service the task again if the data shows something of interest.
We will assume that each vehicle could possibly service a task at any time instant over the horizon N , and that the servicing of a task would mean that the vehicle has collected some piece of data to be transmitted back to the base station.Such piece of data will from now on be referred to as a message, and a new message is created every time step a vehicle is servicing a task.These messages can be divided into even smaller pieces and each piece could possibly be routed differently back to the base station.It will therefore be important to label when each piece is collected in order to constrain the total time elapsed before the whole message is received by the base station.Motivated by Jain et al. (2004), we introduce the optimization variables m pisj to represent the amount of the message with source node s created at time step j which is stored at node p at time step i.Similarly, c pqisj is the transmission rate from node p to node q at time step i of the message with source node s created at time step j.These variables can only take on nonnegative values, and we introduce the constraints m pisj ≥ 0 , (46) The subscripts p and q can take on the value n p + 1, which refers to the base station.Notice that the subscript i of the the variables m pisj and c pqisj start at j in ( 46) and (47) since no message or transmission of message can exist before the message is created.The flow equations relating the introduced variables are given by , and where c p(np+1)isj denotes the data rate at which the message with source s created at time step j is transmitted from vehicle p to the base station (denoted by subscript n p + 1) at time step i.Equations ( 48) and ( 49) can be thought of as representing the conservation of data.In particular, if b sensor pj is true vehicle p is servicing a task at time step j and the amount of data gathered (message size) is given by ∆tc sensor .The amount of data immediately transmitted to other vehicles is given by the second term on the right-hand side of ( 48), and the left-hand side represents the amount that will be stored on vehicle p for the next time step.After the data has been gathered the amount of a specific message stored at node p will remain unchanged, unless parts of this message are transmitted to or received from other nodes, see (49).We assume that each vehicle is equipped with a buffer or a hard drive, and that the buffer can store a limited amount of data.The data, which is the sum of all messages on the node, should be less than the buffer size, h p , that is If a task is not serviced by a vehicle at a specific time instant, then the message size will be zero for the whole horizon.We achieve this with the constraints where the constant M msg := ∆tc sensor is the maximum message size.The delay criterion is ensured by the constraints ∀s ∈ P np 1 , j ∈ I N 1 .These constraints mean that no part of a message will remain on any node p ∈ P np 1 , t delay time steps after it was created on a node s ∈ P np 1 .Furthermore, it also enforces that no data is left on any of the vehicle nodes at the end of the optimization horizon.If instantaneous transmission is required, the constant t delay should be set to 0. Since the base station is only receiving data, we require that To reflect the fact that transmission is only possible when the different nodes are within each others' communication range, we also add to our optimization problem the constraints , where the constant C max is the maximum data rate and bcon pqi is a binary variable which is true if and only if vehicle q is within communication distance of vehicle p at time step i.The constraints required to give bcon pqi this property, are introduced in the section to follow.We also want to bound the collective incoming and outgoing data rate at each time step, and require that The constants C max out and C max in may be different from C max in Equation ( 54), and they may be different for each UAV although we consider them all to be equal for all UAVs in the simulations in Section 4. Depending on the communication equipment, each node may for instance be communicating with multiple other nodes over different bands simultaneously.
Remark 2 Notice that we have chosen to consider the raw flow between nodes in the above formulation.The advantage of considering the raw flow instead of the net flow, is that the raw flow can correctly model asymmetric flow in a duplex communication network.Raw network flows are characterized by nonnegative upper and lower capacity constraints (as in ( 47) and ( 53)-( 56)), and conservation of flow (as in ( 48) and ( 49)).
Remark 3 In the above formulation we have not penalized the total amount of data sent in the network.This means that information can transmitted back and forth between nodes, or in loops when there are three or more communication nodes.This behavior does not influence the feasibility of the MILP problem, but the planned path may unnecessarily be considered infeasible in SPLAT! with respect to the radio path loss design criterion introduced in the next section.This behavior can for instance be avoided by penalizing in the objective function the amount of data transmitted in the network.In our case the problem will not be dominant since we penalize the total time taken from data is gathered and until it is received by the base station.

Connectivity constraints
When formulating the MILP, we will assume that the ability of node p to successfully transmit data at a specified rate to node q, at some time instance i, depends on whether the relative distance between the two nodes is below a certain threshold.This threshold would typically depend on the antenna gains of the receiver and transmitter node, surrounding terrain, data rate, radio frequency band, etc.We stress that R con pqi is not necessarily equal to R con qpi , that is, the threshold depends on the direction of communication.Instead of requiring that node q is within a sphere of radius R con pqi of node p, we require that node q is within a polyhedron that approximates the sphere.The approximation is formed by taking the inner product of the vector χ pqi := (x pi − x qi , y pi − y qi , z pi − z qi ) and ξ kl , where ξ kl was defined in (5) with θ k := 2πk/D con and φ l := 2πl/D con and D con is some constant even integer greater or equal to four.As already pointed out, we introduce binary indicator variables bcon pqi such that bcon ∀p, q ∈ P np+1 1 , i ∈ I N 1 , k ∈ {1, . . ., D con /2}, l ∈ {1, . . ., D con }, that is, the indicator variable bcon pqi is true if and only if, node p can directly transmit to vehicle q, where p, q = n p +1, denote the base station.The logical statement in (57) can be achieved by the introduction of a number of additional optimization variables, but out of brevity we will simply refer to Bemporad and Morari (1999) for this procedure.
Initially R con pqi should be chosen large enough to give a feasible solution to the MILP problem.The solution will then be analyzed in SPLAT!, meaning that the radio propagation path loss between all nodes will be calculated at every time step.We denote the maximum allowed path loss between vehicle p and q to maintain the desired bandwidth during transmission by L max pq .If, for any time step i, the path loss calculated by SPLAT! is below L max pq while at the same time c pqi = 0, then R con pqi will be reduced and we solve the MILP problem again.

Mission set-up
We will consider three different scenarios all including two UAVs.Scenarios 1 and 2 show simulations for the proposed method when delay tolerant communication is not allowed and allowed, respectively.Scenario 0 does not include the rather complex communication constraints we have proposed in this paper, and serves as a justification for why this is necessary in order to satisfy the communication constraints.The tasks and their corresponding waypoints for all scenarios are listed in Table 1.The initial positions of the vehicles with respect to the local ENU frame are given in Table 2, whereas the base station is located at (3750,1450,175) .We assume that a human operator will bring the UAVs up to the initial position, and down from the final waypoint, as this is a common approach in practice.
Table 3 shows the parameters used in the MILP problem.The default solver parameters of the IBM CPLEX Optimizer 12.6 were used, except for epgap which was set to 1.0 × 10 −2 .The algorithm was run on a HP EliteBook 8540w, with Intel Core i7 CPU Q720 @1.6 GHz, 16 GB RAM and a Windows 7, 64bit operating system.Furthermore, we used MATLAB   Wright (2011(Accessed June 29, 2011), and the parameters used are shown in Table 4. SPLAT! is mainly intended for ground-based antennas.However, it is possible to adjust the antenna height, and we have therefore set the antenna height equal to the altitude of the vehicles when analyzing the radio path loss.The point-to-point radio path loss is calculated using the Irregular Terrain Model (ITM), Longley and Rice (1968), with 60% clearance of the first Fresnel zone.It should be noticed, however, that there are some constraints to the parameters used in the ITM.For instance should the antenna heights be less than 1 km above ground level, and the relative distance should be between 1 km and 1000 km.We take a conservative approach and use the maximum of the path loss calculated using ITM and the free space model.Due to the additional constraint of d TIN = 100 m safety distance to the ground, the attenuation due to terrain shielding is likely to be small for short relative distances, and the free space path loss should be a good approximation to the actual path loss.We set the maximum allowed path loss L max pq = L max = 98 for any p, q ∈ P 3 1 , p = q, although it can be different between different nodes, and does not even have to be symmetric between a pair of nodes.We initially set R con pqi = 750 ∀p, q ∈ P 3 1 , p = q, i ∈ I N 1 .How to choose the initial value of R con pqi can for instance be guided by using the area prediction mode in SPLAT!, which gives a prediction of the regional coverage from the base station.If the calculated path loss at time  6 Polarization (0 = Horizontal, 1 = Vertical) 0 Fraction of situations (50% of locations) 0.5 Fraction of time (50% of the time) 0.5 step i is below L max while at the same time c pqi = 0 or c qpi = 0, then R con qpi = R con pqi = R con pqi − r con , so as to maintain the symmetry of R con in this particular example.We use r con = 150.The horizon is 160 s for both Scenarios 1 and 2. With a discretization step of 5 s, we find that N = 32.

Scenario 0: Two UAVs, without communication constraints
An interesting question is whether or not the solution obtained by iterating between the optimization problem and the simulator SPLAT! can justify the increased complexity.We have therefore solved the pure task assignment problem without taking into account any of the communication constraints, in order to be able to compare the performance.The cost to be minimized is given by where J acc and J finish are given in ( 7) and (39), subject to the constraints (1)-( 4), ( 6), ( 8)-( 38), ( 40)-(42).Consider therefore Figure 7 which shows a top view of the planned path, and Figure 8 which shows the calculated path losses.During the servicing of a task we have assumed that the UAV will send sensed information on the route with the least maximum path loss, either direct or by relaying through the other UAV.Even by chosing the best possible data transmission route, we can see that the path loss calculated using SPLAT!(solid black) violates the maximum path loss of 98 dB (dotted black) during communication (shaded areas).This shows that the introduced iterative interplay between optimization and simulation is necessary in order to enforce the design requirements with respect to path loss.By comparing Table 1, Figure 8 and Figure 7, it can be noted that UAV 1 communicates directly with the base station during the servicing of task W 2 .UAV 2 also communicates directly with the base station during the servicing of task W 3 .UAV 2 then services task W 4 , first by relaying information to base station via UAV 1, then by direct communication.Finally, task W 1 is serviced by UAV 1, by relaying information to the base station.The UAVs fly in almost straight lines between waypoints.This is partly because there are no communication path loss constraints enforcing synchronization of their motion during relaying.Further penalizing of the acceleration constraints could have resulted in planned trajectories of more accuracy with respect to vehicle dynamics.
Figure 8: Scenario 0: The calculated path loss using SPLAT!(solid black), maximum allowed path loss for communication 98 dB (dashed black) and time steps when high bandwidth communication is required (shaded grey area).Notice that, the communication path loss requirement of 98 dB is violated at several instants.
Figure 9 shows a top view of the planned path for different iterations.The black circles represent the waypoints, and the yellow circle is the initial location of the UAVs, which East-North location coincides with the location of the base station.In particular, take notice that the trajectories planned at the first iteration are rather straight lines, while after the last iteration they are optimized in order to satisfy the communication path loss constraints during relaying.Figure 10 shows the calculated path loss using SPLAT!(solid black), maximum allowed path loss for communication 98 dB (dashed black), communication data rate (solid blue) and time steps when high bandwidth communication is required (shaded grey area).Comparing this figure with Table 1 and Figure 9b we see that UAV 1 and UAV 2 serve tasks W 2 and W 3 , respectively, and transmit the sensor data directly to the base station.
Then UAV 2 serves task W 4 by relaying data by UAV 1 to the base station, before UAV 1 uses UAV2 as a relay node when servicing task W 1 .Figure 11 which shows the relative distance between the nodes at the final iteration.Also here the shaded areas represents time steps the nodes are communicating.The relative distance (solid lines) should preferably be below the MILP maximum relative distance (dashed lines) at these time steps.This may not always be the case, since we used an over-approximation of the Euclidean norm of the relative between the nodes in (57).Notice that the dashed lines are no longer constantly equal to 750, but have been reduced at some time steps during the iterations.Figure 12 shows the altitude of the vehicles at the final iteration.
Figure 10: Scenario 1: The calculated path loss using SPLAT!(solid black), maximum allowed path loss for communication 98 dB (dashed black), communication data rate (solid blue) and time steps when high bandwidth communication is required (shaded grey area) to accomplish the task at the 3rd and final iteration.Notice that the communication path loss requirements of maximum 98 dB between all nodes are now satisfied.
Figures 13 to 16 correspond to Figures 9 to 11 in Scenario 1.In addition, the amount of data stored at each node is shown in Figure 17.For easier comparison between buffered and transmitted data, we have used the unit bit for stored data.Comparing Figure 9 with 13 it can be noticed that the UAVs can fly in slightly straighter paths between the tasks when delayed networking is allowed.This means that the total mission time is also shorter.The UAVs will still not fly in straight lines for two reasons: 1) We are minimizing the total time of the mission, and not the time spent by the individual UAVs.Hence, if one UAV will have to take a longer route than the other, the UAV with the shortest route will not be hurried to the final waypoint, but may rather serve as a relay vehicle.This may, or may not be a desirable feature depending on whether it is most important with short flight time or getting the sensor data transferred back to the base station as soon as possible.2) For better comparison, Figure 11: Scenario 1: The relative distance between nodes (solid black), and the MILP relative distance constraint for communication (dashed black) at the 3rd and final iteration.The shaded grey areas are time steps when high bandwidth communication is required, and therefore represents time instances when the solid line should preferably be at or below the dashed line.The reason why this is not the case here, is due to the fact that the solid line in the figure represents the actual relative distance, where as in the MILP problem we have used an over-approximation of the relative distance.
The violation is irrelevant for the objective of our algorithm, which is to ensure that the communication path loss requirements are satisfied.
we are using the same initial conditions as in the previous simulation scenario.This means that the maximal distance between the UAVs is bounded by R con which was set to 750 meters.As a result, seen from Figure 15 the UAVs will typically relay sensor information when possible, and delay the transmission only when necessary.Taking a close look at Figure 15 and comparing it with Table 1, Figure 13b and Figure 17 we see that first UAV 1 and UAV 2 serve tasks W 2 and W 3 , respectively, and transmit the sensor data directly to the base station.Then UAV 2 starts collecting data at the waypoints of task W 4 , even if it is not possible to com-Figure 12: Scenario 1: Altitude of the vehicles at the final iteration.The shaded region is the terrain height, and the dashed lines are the upper and lower bounds on the altitude of the vehicles.The actual altitude of the vehicles is depicted with the continuous line during the mission, and even after they reach the final waypoint.When the final waypoint is reached they are expected to be landed manually by a pilot.
municate with the base station neither directly nor by relaying data through UAV 1. Instead it will have to store the data onboard until it gets into communication range of UAV 1.It will then start relaying both the sensed and the stored data through UAV 1, hence the increased data transmission rate seen in Figure 15.Finally, task W 1 is then served by UAV 1, using UAV 2 as a relay node.From Figure 15 alone, there is no apparent reason why the transmission should be interrupted since the communication path loss criteria is satisfied between all nodes.Taking a close look at Figure 11 reveals that the relative distance constraint between UAV 1 and UAV 2 has been tightened at the time instant, and the transmission will therefore need to be interrupted.The sensed data is temporarily stored on UAV 1, as seen from Figure 17.

Discussion
As pointed out in the introduction the intention of this paper is to provide a coarse motion-and communication-plan for a complex surveillance mission.In addition, an online re-planner might be used to account for model uncertainties and effects that become apparent in real-time.We will now discuss some of the uncertainties that should be captured by such a re-planner.Here, a simple point mass model has been used for the UAV, whereas the maneuverability of the physical UAVs will highly depend on its design parameters such as weight, engine power and aerody- namics including lift-and drag forces.We have also ignored some of the environmental effects, and in particularly wind will affect the lift and drag of the UAVs.
In addition the communication properties can be heavily attenuated by rain at certain frequencies, and although SPLAT! is a relatively advanced radio propagation model there will be local effects that are not captured.The sampling based nature of the problem, and in particular when analyzing the radio propagation path loss in the simulator, means that the results may be quite different between samples.Furthermore, the Doppler effect is not taken into account, which may have some effect for high velocity flights.A general recommendation would be to use conservative constraints in the offline path planner, in particular for the vehicle speed and the radio propagation constraints.This means that there would be some robustness to changes in radio propagation path loss between samples, and that vehicle speed can be decreased or increased in order to compensate for wind-effects and lack of abil-Figure 14: Scenario 2: Altitude of the vehicles at the final iteration.The shaded region is the terrain height, and the dashed lines are the upper and lower bounds on the altitude of the vehicles.The actual altitude of the vehicles is depicted with the solid line.When arriving at the final waypoint, the UAV will be brought down by remote control.
ity to perform planned maneuvers.Finally, a minor drawback of the method was illustrated in Section 4.4 where the relative distance constraints of the MILP problem may cause interruption in the communication even when the actual communication path loss is below the design criterion.In some cases this could be remedied by only tightening the relative distance constraint by a little amount at each iteration.This will not work in general as the order in which the tasks are served may change from iteration to iteration and therefore the necessity of transmission at a certain time instant will also change.
The motion-and communication-planning algorithm described in this paper is only intended for offline use, and computational complexity may therefore be of less importance than for real-time algorithms.It is however important to realize that low computation time is desirable also for offline methods.An example is the inspections of critical infrastructure with a known failure, where it is important to locate the failure as soon as possible for it to be dealt with.Reducing the computational complexity has not been emphasized in this work, but in this section we will discuss the many means available for decreasing the computation time.For instance, the code could be made more efficient by estimating bounds on the optimization variables in advance such that tighter constraints could be provided to the solver; bounds on the costs of the optimization problem could be estimated and added as constraints; the required accuracy of the optimal solution could be reduced by adjusting the solver parameters (it should be related to the accuracy of the model pro- ing SPLAT! (solid blue), maximum allowed path loss for communication 98 dB (dashed blue), communication data rate (solid blue) and time steps when high bandwidth communication is required (shaded grey area) to accomplish the task at the 3rd and final iteration.
vided); a non-optimal solution to the problem could be provided at start-up, Vitus et al. (2008); Kuwata and How (2011); the computation could be (decomposed and) run in parallel over several computers, Gunnerud and Foss (2010); Kuwata and How (2011); the computational complexity could be reduced by reducing the number of binary variables, Vielma et al. (2010); Prodan et al. (2012).In general it takes longer to find a feasible solution when we allow for ferrying due to the additional variables we have introduced.However, ferrying means that it is possible to find a feasible solution with lower maximum allowed path loss and/or for tasks placed further away from the base station.Also, it does not require coordination between vehicles, and can therefore be considered a more robust alternative to relaying.One of the benefits of using MILP is that it gives a global optimal solution, and is therefore suitable for the benchmarking of other, possibly sub-optimal methods.One such method is Particle Swarm Optimization (PSO), Kennedy and Eberhart (1995), which has successfully been applied for optimal path planning for UAVs in e.g Sujit and Beard (2009), Ho et al. (2013) Figure 16: Scenario 2: The relative distance between nodes (solid black), and the MILP relative distance constraint for communication (dashed black) at the 3rd and final iteration.See also caption text for Figure 11.and Ho et al. (2015).However, it turned out that PSO is not very well suited for an implementation of an optimization problem similar to that of Section 3. By just requiring the UAVs to visit all the waypoints (and disregard any communication constraints), the resulting paths using PSO were comparable to the one achieved with MILP, but with computation time exceeding the computation time of MILP with communication constraints.The dimension of the search space was in this case n p × N .By including the communication constraints in the PSO the dimension of the search space would increase to n 2 p × (n p + 1) 2 × N 3 which leads to very high computational complexity.The key problems of using PSO for this scenario were identified as: a high-dimensional search space; PSO lacks the ability to implement hard inequality constraints, which resulted in a large number of cost function terms that were difficult to weigh in case of competing objectives (e.g minimizing the mission time, while constraining the velocity of the UAVs); PSO lacks the ability to implement hard equality constraints, which resulted in a problem of conserving the data flow in and out of network nodes.Notice that the increments are typically of 10 Mbits, which is equivalent to c sensor ∆t, i.e. the rate at which data is collected multiplied with the sample time.

Conclusion
In this paper we have addressed the problem of task assignment and path-and communication-planning for multiple UAVs.In the MILP formulation a general communication topology is allowed such that sensor data can be transferred efficiently back to the base station by direct transmission, ferrying or relaying.Tolerating delays introduces flexibility into the motion-and communication-planning problem, and means that the range between tasks and base station can be extended, and/or that the maximum allowed path loss can be reduced.We have improved the accuracy of the communication properties by using a path loss simulator to analyze the paths generated from the optimization problem.

Figure 1 :
Figure 1: This figure illustrates how delay tolerant net-working can be used not only to extend the surveillance range, but also to increase the total data transmission rate.In the dark red areas the UAV can communicate with the base stations at a high transmission rate, whereas the transmission rate is poor in the light red areas.By bringing the UAV into range of the dark red areas, the total data transmission rate can be substantially increased even when taking into account the UAV's flight time.To quoteTanenbaum (2003): "Never underestimate the bandwidth of a station wagon full of tapes hurtling down the highway".The process of storing data on a mobile node and physically bringing it from one location to another is called ferrying.
p + 1} • I b a = {a, a + 1, . . ., b} a, b ∈ Z, -sample time index • D l -set of indices to triangles that have a common vertex P l , where l ∈ {1, . . ., m l } • T b a = {a, a+1, . . ., b} a, b ∈ Z -set of task indices • W t -indices of the waypoints belonging to task t, where t ∈ T n l 1 .Every waypoint belongs to one and only one task, that is W r ∩W s = ∅ ∀r, s ∈ T n l 1 and r = s • T + -tasks with more than one waypoint • W first t , W last t -index to the first, respectively, the last waypoints of task t ∈ T + Optimization variables:

Figure 3 :
Figure 3: The Earth-Centered Earth-Fixed (ECEF) frame with axes X, Y , Z, and the East-North-Up (ENU) with axes x, y, z.The origin of ENU frame is typically represented in terms of the reference ellipsoidal parameters longitude l, and geodetic latitude µ.

Figure 5 :
Figure 5: In 5a a top view of a TIN with vertices P l , l ∈ {1, . . ., 7} is presented, and 5b shows the corresponding sets D l of indices to triangles which have P l as a common vertex.

Figure 7 :
Figure 7: Scenario 0: Top view of the planned path.The black circles represent the waypoints, and the yellow circle is the initial location of the UAVs where also the base station is located.The red circles represent UAV 1, and the green circles represent UAV 2. The contours in the background represent the height of the terrain.

Figure 9 : 1 :
Figure 9: 1: Top view of the planned path at (a) 1st iteration, and (b) 3rd and final iteration.See also caption text for Figure 7.

Figure 13 :
Figure 13: Scenario 2: Top view of the planned path at (a) 1st iteration, and (b) 3rd and final iteration.See also caption text for Figure 7.

Figure 15 :
Figure 15: Scenario 2: The calculated path loss us-ing SPLAT! (solid blue), maximum allowed path loss for communication 98 dB (dashed blue), communication data rate (solid blue) and time steps when high bandwidth communication is required (shaded grey area) to accomplish the task at the 3rd and final iteration.

Figure 17 :
Figure 17: Scenario 2: The amount of data stored at the nodes at the 3rd and final iteration.Notice that the increments are typically of 10 Mbits, which is equivalent to c sensor ∆t, i.e. the rate at which data is collected multiplied with the sample time.

Table 1 :
Tasks, and coordinates of their corresponding waypoints of Scenarios 1 and 2

Table 2 :
Initial position of UAVs