Dynamic Motion Planning for Aerial Surveillance on a FixedWing UAV
Abstract
We present an efficient path planning algorithm for an Unmanned Aerial Vehicle surveying a cluttered urban landscape. A special emphasis is on maximizing area surveyed while adhering to constraints of the UAV and partially known and updating environment. A Voronoi bias is introduced in the probabilistic roadmap building phase to identify certain critical milestones for maximal surveillance of the search space. A kinematically feasible but coarse tour connecting these milestones is generated by the global path planner. A local path planner then generates smooth motion primitives between consecutive nodes of the global path based on UAV as a Dubins vehicle and taking into account any impending obstacles. A Markov Decision Process (MDP) models the control policy for the UAV and determines the optimal action to be undertaken for evading the obstacles in the vicinity with minimal deviation from current path. The efficacy of the proposed algorithm is evaluated in an updating simulation environment with dynamic and static obstacles.
1 Introduction
In the last decade, Unmanned Aerial Vehicles (UAVs) have gained importance in a myriad of military and civilian applications. They are especially useful in missions pertaining to surveillance, rescue and exploration. In these tasks UAVs are often required to operate in unfavourable conditions where a complete description of operating environment may not be available.A fully autonomous UAV eliminate the need to plan the mission beforehand in an unknown , unpredictable environment, it is capable of handling any new situation which is thrown upon it spontaneously. Consider the following scenario:An urban setting where a UAV is released in an unknown airspace for the purpose of surveillance. The aim of the UAV is to maximize mapped area with constraints such as minimum turning radius and avoiding collision with other planes/UAVs operating in vicinity, buildings , other dynamic and unknown obstacles.In such a challenging scenario fully autonomous path planning plays a vital role in ensuring success of the mission specially for fixedwing UAVs which lack hovering capability and hence need to maintain a minimum airspeed to maintain lift.
Several techniques have been proposed for 3D path planning like Rapidlyexploring Random Graph[1],Visibility Graphs[2],PRM[3],RRT*[4],Artificial Potential Field. All these techniques though effective,are computationally expensive and need to be recomputed as soon as the environment is changed.The above mentioned techniques suffer the curse of dimensionality and a combinatorial explosion takes place with increase in search space.
This paper proposes a 3D Dynamic path planning approach which overcome the aforementioned limitations and facilitates conduction of fully autonomous missions. The rest of the paper is organized as follows In Section II:Scenario describing the environment,obstacles and coordinate space is presented.Section III comprises of proposed approach giving details of global and local path generation along with the decision process for optimal control policy generation.The simulation results and analysis for the proposed approach are presented in Section IV.We conclude with final remarks and discuss future works in Section V.
2 Scenario
A fixedwing autonomous UAV is to conduct surveillance of a prespecified area(Search Space) within a given Geofence.Considering the given scenario, the aim of our approach is to:

Maximize the area surveyed in minimum flight time.

Avoid collision with any obstacle, be it stationary or dynamic.
The UAV has a constraint on its turning radius , climbrate and minimum airspeed.The path planning problem is finding the optimal path in the search space given the constraints on UAV and environment.
2.1 Environment
We assume that the UAV is GPS enabled and carries some specialized sensors capable of accurately detecting the location of the obstacles. The search space is assumed to be partially known for the PRM and the approach is designed keeping in mind that any kind of obstacle can turn up at any coordinates on any given time.
Static Obstacles
Static or Stationary obstacles are assumed to be a solid cylinder. The coordinates of the centre of the cylinders along with its radius and height is updated whenever there is a change in the environment.
Dynamic Obstacles
Dynamic or Moving obstacles are assumed to be solid spheres. Similar to static obstacles the coordinates of sphere’s centre, radius and velocity are synced with the changes in the environment.The trajectory of the dynamic obstacles are randomly selected from a pool of precomputed paths so as to mimic real life aerial vehicles that may operate in the vicinity of the UAV.
2.2 Coordinates Transformation
WGS 84 to ECEF coordinates
Most navigation systems are based on World Geodetic System including GPS. WGS is modelled after the standard coordinate system for Earth, a standard ellipsoid reference surface for raw altitude data, and a gravitational equipotential surface that defines the nominal sea level. We need to convert to a Local Tangent Plane Coordinate System for the ease of mathematical computation so as to make the problem retractable. The order of approximation involved in the conversion is of the magnitude of for the small distances involved for the scenario of our paper. Geodetic coordinates(latitude , longitude , height h) can be converted to ECEF coordinates using the equations:
(1) 
(2) 
(3) 
where
(4) 
here and are semimajor axis and the first numerical eccentricity of the ellipsoid respectively. The values of which can be found in the definition of WGS 84[5].
ECEF to WGS 84
We need to transform the ECEF coordinates back to geodetic coordinate system so as to facilitate the waypoint planning on the onboard autopilot system.
(5) 
(6) 
(7) 
where auxiliary values are:
(8) 
(9) 
3 Proposed Approach
3.1 Overview
Two phase motion[6] planning is a widely adopted strategy for reducing the planning complexities of highorder dynamic systems[7][8]. At the first stage, a coarse but feasible path is generated taking into consideration the kinematics of the system and the known static elements of the environment. System design constraints and dynamic elements in the environment may be ignored in this phase. The primary objective of this phase is to confine the search space to the vicinity of the generated path and thus reduce the complexity of the planning task. At the second stage, a smoother more finely attuned path is generated in response to changes in the environment and dynamics of the system. The two stages compensate each other in terms of planning complexity, risk associated with the path and type of obstacles handled. Two phase planning combines the merits of both the methods and facilitates realtime operation in a cluttered environment.
The major components of our system along with the flow of control are represented in system architecture shown in the Fig.1. Global planner works in a 3D occupancy grid with boundaries demarcated by the geofence. The global planner takes starting location, geofence and known static elements of environment as input and generates a tour of search space with objective of maximizing area coverage while minimizing collision risk. The output of the global planner is a sequence of waypoints representing the tour. Certain waypoints in the search space are marked as milestones to avoid being trapped in local minima by virtue of heuristic search algorithm and to ensure maximum area coverage (details in later sections).Local planner works alongside the global planner and caters to moving obstacles and unforeseen stationary obstacles. Local planner checks for potential collision on the current path and assesses the risk associated with it. In case current path is deemed risky, the path is modified with a sequence of dynamic waypoints to avert the danger. The trajectory of the modified path is based on smooth Dubin curves to ensure a dynamically feasible and steerable path. A decision process determines the type of path adjustment to be undertaken by the global or local planner and accordingly waypoints maybe modified. If both original and modified paths are deemed risky, the current milestone is marked unsafe and local planner computes a route to the nearest unvisited safe milestone. In absence of such a milestone UAV is guided to last known safe milestone. When each milestone has been visited the UAV returns to the starting location and lands. A detailed explanation for each of the phases follows in the subsequent sections.
3.2 Global Planner
As mentioned earlier, the objective of global planner is to compute a kinematically feasible tour of the search space which maximizes area covered and minimizes collision risk. Initially, a discretised 3D occupancy grid is computed from the search space with each cell value indicating the risk associated with it. The grid is then sampled randomly for safe cells. A probabilistic road map (PRM) is built using the computed samples. The samples are clustered into safe regions and the centroid of every safe region is considered a milestone. Starting from the initial location, a greedy approach is used to form a tour going through all the milestones. The optimal path between individual milestones is determined by A* algorithm. The path so formed is bound to be kinematically viable as for each node in the path, only successors following the constraints are considered for path expansion.
Building Occupancy Grid
The search space within the confines of the geofence is represented as a discretised 3D grid. Each cell of the grid is associated with a collision risk involved with the cell. A global occupancy score for a cell is determined from the obstacle field distribution in the search space. Both static and dynamic obstacles contribute to this obstacle field. The global occupancy score for a cell is expressed as
(10) 
where is the count of known obstacles, is the Euclidian distance of the cell centred at (x,y,z) from the obstacle and is the radius of the obstacle. This score provides an accurate estimation to the proximity of the obstacle but leads to a coarse distribution of values in the occupancy grid around the obstacles. Thus a local score is also considered alongside the global score to facilitate a evener distribution of values in the grid.
(11) 
The local score, denoted as for a cell centred at (x,y,z), is given as the mean of values of all the adjacent cells in the grid.The collision risk for the cell is taken as the maximum of its local and global occupancy score.
(12) 
Voronoi Tessellation of search space
To ensure maximum coverage of search space for surveillance, it is necessary to ensure that the path planning algorithm does not remain confined to any one region as a result of being trapped in local minima which have been a problem with many planners including potential fieldbased methods. Thus there exists a need for division of target space (free space) into disjoint sub regions which must each be covered to facilitate maximum coverage.Voronoi tessellation of search space can provide the required partitions given the seed points. For a typical UAV mission, the number of seed points can be estimated to be roughly equal to the number of identified critical tasks. For a discretized configuration space as considered, voronoi tessellation can be achieved by means of KMedoids Clustering[8] with being set to the number of milestones.The objective function minimized during the iterations is given in Eq.13.
(13) 
where is the number of milestones, is medoid of the cluster, is a sample data point and is the cluster.
Only a randomly chosen subset of safe points in target space are considered for clustering due to the high runtime complexity associated with the procedure.A sample cell is considered safe if the associated collision risk score is less than the collision risk threshold .An adequate number of sample points are chosen to be a good indication of the overall distribution of large target spaces within the search space.
For a given scenario a good heuristic for estimating the number of milestones can be taken as
(14) 
where is the area of geofence, is the anticipated cruise altitude for the UAV for in the mission and is the field of view of the onboard camera.
Building Probabilistic Road Map
The subset of points chosen from target space are taken as vertices for the Probabilistic Road Map construction. KNearest Neighbours (KNN)[10] algorithm is used for determining the vertices adjacent to a given vertex and Euclidian distance is chosen as the comparison metric. The value of K has to be high enough to maximise the probability of having at least one kinematically feasible path when the graph is traversed in realtime but at the same time it is required to keep the value of K as moderate as possible to curtail the runtime complexity. We found that in a typical urban city scenario the value of K varies between 40 and 120 for optimal results.
The initial order in which milestone way points are visited is determined from the nearest neighbour algorithm (greedy algorithm) which lets the UAV choose the nearest unvisited milestone as its next milestone.This technique has been chosen for approximating the shortest tour as it gives a path 25% longer than the optimal path on average. Also, Tassiulas[11] has shown that a tight upper bound of exists for the nearest neighbour tour. Here N is the number of cities and d is the optimal path distance for the Euclidian travelling salesman problem wherein the triangle inequality holds.
Kinematically feasible Path
Kinematic constraints of minimum turning radius and maximum climb rate for the UAV in current state are enforced during generation of motion primitives when the PRM is traversed using A* algorithm. At any given node, only those adjacent nodes are considered for shortest path which satisfies the kinematic constraints with respect to the current state of the UAV. This modification not only reduces the complexity of the graph but also leads to a smoother path for traversing.
Thus the node expansion from the current node depends not only on its distance from the current node and occupancy grid value but also on its spatial orientation from it.The heuristic chosen for selection of next best node from the list of adjacent nodes is given by
(15) 
(16) 
(17) 
(18) 
where is the Euclidian path cost, is the number of times the node has already been explored, is evaluation function for node expansion and for imposing a high penality on revisitation. Revisitation penalty is imposed on an already visited node to ensure that preference is given to new unsearched areas rather than remaining confined to local minima and lose gatherable intelligence.
Pseudocode
Algorithm 1 gives the pseudo code of the approach described in previous sections: {algorithm}[h!] \SetKwInOutInputInput \SetKwInOutOutputOutput \InputGeofence, Initial location and nature of Obstacles, UAV parameters,Source Point \OutputGloabal Path waypoints get_range Convert Geodatic coordinates to ECEF \ForEachpoint in space_matrix \eIfpoint is outside geofence or point is inside obstacle
sample_list Randomly Sample minimum required points points whose score initialize \ForEachpoint in sample_list inititalize node Add node to PRM \ForEachnode in PRM \ForEachnode in neighbour_list \IfCollision_Free(Extend(node,neighbour_list)) initialize Global Path \ForEachnode in milestone
In lines 12 we discretize the search space based on resolution distance along all the three axes.In lines 412 collision risk score is computed for all the cells in the discretized grid.In line 13 random samples are acquired for which the collision risk score is less then the collision risk threshold().In line 14 acquired samples are clustered and centroids of each cluster are set as milestones.In line 15, starting from the source point we rearrange the milestones according to the euclidean distance in a greedy manner such that nearest milestone to source point is taken first then the milestone nearest to the previous milestone which has not already been visited comes next till all milestones have been connected.In lines 1628 PRM is built as described in Section 3.2.3.In lines 2933 the path segments between the consecutive milestones are computed using with heuristic given in Eq. 15.Individual path segments are appended in order to obtain the final global path.
The computational complexities involved with each of the phases in the global planner are given as follows:

Initial occupancy grid initialisation and collision risk score assignment collectively have a time complexity of order of .

Voronoi tessellation of discretized search space on the randomly sampled safe points has an associated time complexity of order of .

Probabilistic Roadmap building phase by means of KMedoid Clustering with max heap optimisation runs in . This phase along with run on chosen heuristic are the dominating factors during the execution of the algorithm.
Here, K = number of obstacles, N=discrete division along one dimmension of search space, m= number of milestones/clusters, n= count of sampled safe points, t= number of iterations, k = count of nearest neighbours.
3.3 Local Planner
The path generated in the grid planning phase though feasible has its limitations. Continuously updating the global path by querying the grid planner in response to changes in the environment is a computationally expensive process and also the path generated by the global planner is jagged. To remedy these limitations and provide a higher resolution of control in motion planning a local planner is built atop the global planner. The motion primitives generated between two configurations represented by successive nodes of the global path ensure that vehicle dynamics strictly adhere.
The local planner proposed generates motion primitives based on shortest Dubins path between the two configurations. The trajectory between these configurations is decomposed onto two orthogonal planes and for each of these 2D Dubins curves are calculated taking into account any obstacles that maybe present in between these configurations.
Dubins Motion Primitives
A node in the global path is defined as follows
(19) 
where the first four components ( describe the location and desired heading vector, represents the target airspeed and is the collision risk score associated with the given node. Given any two consecutive nodes in the global path i.e. (and ( the connecting path is constrained by minimum turning radius, initial and final bearings. The 3D path connecting the two nodes can be decomposed into two orthogonal planes as shown in Fig.7.
Bearing and minimum turning radius constrains UAV’s movement in the horizontal plane while max climb rate and direction of climb(ascent/descent) affect UAV’s propagation in the vertical plane. 2D Dubins[12] curve is created for each of these orthogonal planes. Depending on the distance of separation between the two nodes, either one of CSC or CCC trajectories[Fig.8] is chosen for the shortest path. The benefit of choosing Dubins curves to join the nodes is twofold as the initial and final headings are automatically aligned and it renders the shortest path given the radius curvature constraints.
Obstacle Path Adjustment
If an obstacle lies on the connecting path between the two consecutive nodes, then the Dubins path is modified to steer the UAV clear of obstruction with minimum deviation from the determined path. Depending on nature (static/dynamic) and radius of obstacle different path modification strategies are adopted. If the radius of the obstacle is less than the minimum turning radius of UAV, then its radius is taken to be equal to for computation of viable curves in resultant Dubins path.
Static Obstacle For a static obstacle multiple paths exist that can be generated for safe propagation of UAV, to ensure minimum deviation, path passing through tangentially closest point on the periphery of the obstacle is chosen, given that the change in steering angle for such a point is feasible.Assume that the UAV is initially at and going towards next node at ,then the initial bearing and heading of the UAV is given as and where
(20) 
(21) 
Also consider a static obstacle that appears at with radius and height .The obstacle can be dodged in two ways either by changing the current yaw or pitch angle.
The change in yaw angle() is given as
(22) 
where is expressed as
(23) 
here are xy coordinates of the tangentially closest point on the periphery of the obstacle required to dodge the obstacle and can be given as
(24) 
(25) 
where
(26) 
and
(27) 
In Eq.26 is minimum safety clearance. denotes the euclidean distance between and . The zaxis is taken as the axis of rotation for yaw and hence there is no change in z coordinate.
Similarly the change in pitch angle is given as
(28) 
where is expressed as
(29) 
here are yz coordinates of the tangentially closest point on the periphery of the obstacle required to dodge the obstacle and can be given as
(30) 
(31) 
The xaxis is taken as the axis of rotation for pitch and hence there is no change in x coordinate.The appropriate way to dodge the obstacle by changing pitch ot yaw is determined by wheather and are within UAV’s constraints. If both are feasible the smaller change in angle requirement is followed.
Dynamic Obstacle Consider the inital and target nodes as given in the previous section.Also consider a dynamic obstacle that appears at and velocity .A similar procedure as followed to dodge the static obstacle is followed here with the difference being that are not fixed and hence a sequence of dodge points are generated until the UAV steer clears of the obstacle.
Again the UAV can dodge by either changing the yaw or pitch to generate the next dodge point.At any instant the computation for the next dodge point as the result of the change in yaw angle is similar to that of the static obstacle and is given by Eq. 22,23,24,25.
The change required in pitch angle(heading) is given as in Eq. 28 but here
(32) 
where
(33) 
(34) 
and
(35) 
Here is the dodge point in the sequence, is the previous dodge point in the yz plane and l is the euclidean distance between and .The dodge point generation is terminated when the bearing from the current dodge point to the target point() is no longer intercepted by the obstacle.For more details on and refer to Appendix.
L1 Navigation Controller
A nonlinear navigation controller based on[13] was used to guide the UAV along curved segments of the trajectory. The controller not only works like a PD controller for straight paths but also contains an anticipatory control element for minimizing drifts when following curved paths. The guidance logic works by selecting a reference point at a distance on the desired trajectory in front of the UAV. Then a lateral acceleration command is generated which gives it an adaptive capability to face external disturbances like winds. The expression for lateral acceleration command is given by
(36) 
Where is the velocity of UAV, is defined from current position to reference point, is the angle between and. Details regarding the calculation of the expression and error estimation are omitted due to space constraints and can be found in [13].
3.4 Markov Decision Process
Three sets of actions can be adopted to prevent a collision with obstacles. A control policy is needed for determining the efficacy of each action in a given state. The control policy determines whether variation in velocity, local path adjustment, global path amendment or a combination is required to evade an impending obstacle or tackle changes in the environment. The UAV control problem and decision process is modeled as a finite state Markov Decision Process[14].
We made use of the Policy Search algorithm as given in [15] to learn the critical parameters of our control policy. For the reasons of brevity, only a short description of each of the learned parameters is given in Fig. 9
4 Results and Experiments
We validate our algorithm on the sig rascal model provided by JSBSim as a simulation platform for a fixed wing UAV. The onboard autopilot, control system and sensors are simulated using the ArduPilot platform with APM Planner as the ground control station. A custom C# .Net application was built for controlling and configuring the obstacles and Caesium was used for depicting the final scenarios in realtime. The kinematic constraints for the considered system are given in Table 1.
Parameters  Values 

Turning Radius  
ClimbRate  
Cruise Speed  
Max. ClimbAngle  
Max. Roll Angle  
Max. Pitch Angle  
Min. Pitch Angle  
NAV L1 Damping  
NAV L1 Period  
NAV L1 X TRAC I 
We conducted over 44 full system simulations on a workstation equipped with an octa core intel i7, 2.50 Ghz CPU.In any typical mission there were 36 stationary obstacles of varying radii spread across the search space and 13 dynamic obstacles traversing along predetermined paths. The criterion involved in evaluating our system were the fraction of total search area (enclosed within the geofence boundary) surveyed, length of the path and distance of separation of UAV from the obstacles within the search space.
From Fig.10 and Fig.11 we conclude that the minimum distance of separation of UAV from obstacles decreases as collision risk threshold is increased from 0.2 to 0.8 at the same time the surveyed region increases from 72% to 94.7%.Hence an appropriate collision risk threshold has to be estimated experimentally for the desired UAV.Ideally the value lies between 0.4  0.7 depending on the mission profile.
The choice of number of milestone waypoints for mission also affects the expanse of search space surveyed. Fig. 12(a) shows the increase in the area covered with increase in number of milestones. From Fig.12(b) it is also noticed that the length of the tour also increases substantially.From Fig.13 it is noticed that there exists a tradeoff between area surveyed and path length with respect to milestone count. In our simulations we found that for the geofence considered, there is a consistent increase in area surveyed with gradual increase in path length till 14 milestones, beyond which the path length increases without any appreciable growth in area surveyed. This is consistent with our heuristic as proposed in Eq.14 where for an area of dimensions 834x577 the optimal milestone count should be 14 for a mission in which the onboard camera has fov and anticipated average altitude is 180.
CSpace Dim.  A*  A*PRM  Speedup 

25x30x27  27685  3109  8.90 
40x30x30  31084  5942  5.23 
50x50x50  52031  6232  8.95 
Precomputation of a probabilistic roadmap from sample space greatly reduces planning time. Each node in PRM is extended by connecting to KNearest Neighbours and pruning those extensions which do not adhere to the dynamic and kinematic constraints of the system, then A* is used to find the optimal path between the milestone nodes in PRM. In comparison applying A* over the complete search space takes much longer to find a similar path with no added benefit. Table 2 shows the speedup in computation when a precomputed PRM is used instead of searching through the entire search space. On an average, the global path planning phase runs 7.69 times faster due to this optimization.
Fig.14, 15and 16 illustrate the results of simulations, in which a consistent wind blowing at 2 m/s with the bearing of 120° from North was considered. Nonlinear L1 controller with the aforementioned parameters was used as the waypoint controller to ensure that the UAV does not drift considerably from its trajectory in presence of wind. For all the illustrations the boundary of geofence was demarcated at points A (28.754812 N, 77.115204 E), B (28.754812 N, 77.115204 E), C (28.755188 N, 77.119689 E) and D (28.748698 N,77.120161 E). The UAV after takeoff enters the search space at coordinates at point (28.752088 N,77.116211 E). Fig.14 demonstrates the working of the global grid planner in the presence of obstacles , four stationary obstacles of radii 25,35,28 and 30 are present at locations (28.7536640 N,77.1160412 E), (28.7522719 N,77.1180367 E),(28.7547551,77.1193027) and (28.7503155 N,77.1195173 E). The global planner generates a kinematically feasible path represented by the yellow coarse path as shown in Fig.14(a).The local planner generates interconnecting motion primitives which result in a smooth dynamically feasible path which is traversed by the UAV as shown in Fig.14(b).
Fig.15 depicts the response of the local planner when an additional fifth obstacle appears unexpectedly at (28.7502026 N,77.1157837 E) of radius 32.The corrected path is shown in Fig.15(b).We see that an additional dodde point has been added to the current path so as to avert the collision.
Fig.16 illustrates the case where in addition to previous obstacles a dynamic obstacle of radius 30 approaches the path of the UAV with velocity 15 at distance of 105.A path alteration was made by the local planner to steer the UAV away safely from the oncoming obstacle. The UAV was able to stay clear of the obstacles and maintain a minimum separation of more than 3 which was taken as the safety clearance.
5 Conclusion
In this study, a dynamic motion planning approach has been presented for a fixed wing UAV surveying urban environment. The specialty of the proposed approach lies in its utility to work in partially known and updating environment.The computation involved in generating feasible trajectories for the UAV through the search space is greatly reduced by the use of precomputed PRM.Voronoi Tessellation of search space and identification of key waypoints in the form of milestones leads to efficient mapping of the region to be surveyed.The changes in the environment are handled effectively by the decision process in the form of local or global planner.The application of 3D Dubins curve leads to smooth and dynamically feasible trajectories at runtime.The efficacy of the proposed approach is demonstrated in the complex simulation scenarios in highly cluttered environments comprising of both static and dynamic elements.In future, we plan on extending the capability of the proposed system by integrating an array of sensors to impart perception ability to our UAV platform.
Appendix A Derivation for dodge point coordinates
Equation of the circle centered at is given as:
(37) 
Let the Equation of the tangent from point A be of the form
(38) 
then,after substituting and Eq. 37 can be written as:
(39) 
Since line AD is a tangent to the circle, therefore both the roots of Eq. 39 are equal and hence the discriminant is zero. Using this we obtain the following quadratic relation in m:
(40) 
where
(41) 
(42) 
Solving the quadratic Eq. 40 we get the roots as:
(43) 
where m is the possible values of the slope of line AD. The slope of line CD which is perpendicular to line AD is given by . Thus coordinates of the dodge points in the yaw plane are obtained as in Eq. 24 and 25 using parametric form equation of a line. Similar steps can be followed for pitch plane to obtain coordinates of dodge points.
References
 Karaman S, Frazzoli E. Incremental samplingbased algorithms for optimal motion planning. arXiv preprint arXiv:1005.0416, 2010
 Geraerts, Roland. "Planning short paths with clearance using explicit corridors."Robotics and Automation (ICRA), 2010 IEEE International Conference on. IEEE, 2010.
 Karaman S, Frazzoli E. Samplingbased algorithms for optimal motion planning[J]. The International Journal of Robotics Research, 2011, 30(7):846894.
 DuToit R, Holt M, Lyle M, et al. UAV Collision Avoidance Using RRT* and LOS Maximization Technical Report# CSSE1203[J].
 Anonymous,"World geodetic system 1984 (wgs84)its definition and relationships with local geoderic systems," Defense Mapping Agency,Tech. Rep. NIMA TR 8350.2, 2000.
 M. Hwangbo, J. Kuffner and T. Kanade, "Efficient Twophase 3D Motion Planning for Small Fixedwing UAVs," Proceedings 2007 IEEE International Conference on Robotics and Automation, Roma, 2007, pp. 10351041. doi: 10.1109/ROBOT.2007.363121
 J. E. Bobrow, "Optimal robot path planning using the minimumtime criterion," IEEE Trans. on Robotics and Automation, vol. 4, no. 4, pp.443450, 1988.
 Y. Kuwata and J. How, "Three dimensional receding horizon control for uavs," AIAA Guidance, Navigation, and Control Conference and Exhibit,Aug. 2004.
 Jin X., Han J.,"KMedoids Clustering",Encyclopedia of Machine Learning, 2010, Springer US, Boston, MA, pp. 564565, doi:10.1007/9780387301648_426
 ltman, N. S."An introduction to kernel and nearestneighbor nonparametric regression". The American Statistician, 1992,pp. 175?185. doi:10.1080/00031305.1992.10475879
 L. Tassiulas, "Worst case length of nearest neighbor tours for the Euclidean traveling salesman problem," SIAM J. Discrete Math ., Vol. 10, No. 2, pp 171179, May 1997.
 L.E.Dubins, On curves of minimal length with a constraint on average curvature and with prescribed initial and terminal positions and tangents,? American Journal of Mathematics, vol. 79, no. 3, pp. 497516, jul 1957.
 S. Park, J. Deyst, and J. P. How, "A New Nonlinear Guidance Logic for Trajectory Tracking," Proceedings of the AIAA Guidance, Navigation and Control Conference, Aug 2004. AIAA20044900.
 Martin L. Puterman,"Markov Decision Processes: Discrete Stochastic Dynamic Programming (1st ed.)",1994, John Wiley & Sons, Inc., New York, NY, USA.
 Andrew Y. Ng, Michael I. Jordan, "Policy Search Method for Large MDPs and POMDPs", CoRR, 2013, Vol. abs/1301.3878