Recently, our community has witnessed the development of planning methods for quadrotors. Spline-based methods, which decompose the spatial and temporal parameters of a planning problem and focus on optimizing its spatial part, are widely applied for real-time applications [1, 2, 3].
Although spline-based methods can efficiently and accurately generate energy-optimal solutions for online usage, they usually omit temporal planning for simplicity. A typical spatial-temporal joint planning problem has high nonlinearity and nonconvexity coming from its objective and constraints. Since temporal planning has underlying coupling with spatial parameters and implicit gradients, the spatial-temporal joint optimization cannot be solved by general nonlinear programming (NLP) in real-time. Even though existing methods can provide online motion planning without temporal planning, they are often too conservative to be used for autonomous flights with high aggressiveness. To bridge this gap, we propose a framework to split the spatial and temporal aspects of a trajectory optimization problem, then solve them alternately. With our method, we can obtain the energy-time joint optimal trajectory in milliseconds.
The proposed method is based on the naturally alternating structure of the spatial-temporal trajectory optimization and designed to have guaranteed optimality and efficiency. Our method is motivated by the fact that a polynomial trajectory with an odd order can be uniquely determined by its endpoint derivatives, i.e., the boundary condition, and the time duration. For a piecewise polynomial, once all boundary conditions are fixed, each piece of the trajectory depends only on its time duration, which can be optimized separately. By utilizing the widely-adopted linear-quadratic objective of the optimization, the optimal time durations can be updated efficiently. Moreover, inspired by , the closed-form solution is adopted to update derivatives on waypoints. Based on the above observations, the joint optimization can be efficiently processed by an alternating minimization (AM) procedure . With our method, a large-scale joint optimization can be done in a few milliseconds to generate optimal trajectory with the best time allocation.
To the best knowledge of us, the proposed method is the first one that generates trajectories for a quadrotor with the spatial-temporal optimality, in such a limited time. Summarizing our contributions in this work:
An unconstrained alternating minimization algorithm is proposed to generate spatial-temporal optimal trajectories efficiently, with a proven non-asymptotic rate of the global/local convergence.
A computationally efficient feasibility check method is designed for a wide range of constraints in our method.
A constrained alternating minimization algorithm is constructed to optimize feasible trajectories in a recursive fashion, with global convergence verified.
The proposed method is integrated into an autonomous quadrotor system then evaluated by real-world experiments as well as extensive benchmarks. The source code is released for the reference of the community.
In what follows, we discuss related literature in Sec. II. Preliminaries of this paper are given in Sec. III. The proposed spatial-temporal trajectory generation method for unconstrained and constrained planning cases are detailed in Sec. IV and V, respectively. Experiments and benchmarks are given in Sec. VI. The paper is concluded in Sec. VII.
Ii Related Work
For quadrotor planning, polynomial splines have long been used for trajectory parametrization since , because of their flexibility and analytical convenience. In , the minimization of squared derivatives is used as the objective of quadratic programming (QP), which can be solved efficiently and accurately. Based on this formulation, intensive works have been recently proposed. A method for obtaining a closed-form solution of the above QP program is proposed in , where a safe geometric path guides the generation of the trajectory to ensure its safety. By recursively adding intermediate waypoints to the path, a safe trajectory is finally generated after solving the minimum-snap problem several times. In [6, 7, 8]
, safe and dynamically feasible trajectories are online generated within a safe flight corridor, which excludes all obstacles in complex environments. However, in these methods, the time allocation of the piece-wise trajectory is pre-defined or online adjusted by heuristics. Although these heuristics are cheap to compute, the trajectories generated are often far less optimal and over-conservative, making them incapable of high-speed flights.
To address the time allocation problem, Mellinger et al. 
compute the projected gradient with respect to durations on a hyperplane where the sum is fixed. They optimize time allocation through backtracking gradient descent. Temporal scaling is applied on the solution until dynamical feasibility is achieved. Both the finite difference and scaling used in this method are expensive operations when the number of trajectory pieces is large. Liu et al. propose a way to calculate the proper scaling factor such that a single scaling operation suffices, while it only applies for rest-to-rest trajectories. Sun et al. 
formulate the problem as a two-level optimization. They estimate the projected gradient analytically through the dual solution of the low-level QP, which improves accuracy compared with the numerical gradient. Nonetheless, the projected gradient is still inconvenient to compute and inefficient for nonlinear optimization. To avoid this situation, Richter et al. use total duration as a regularization term, thus making each duration an independent variable. The time allocation is optimized through gradient descent, while actuator constraints are also fulfilled by scaling. However, the optimal ratio of time allocation under the constrained case may differ a lot from the unconstrained case. Consequently, scaling can ruin a trajectory where constraints are violated on a very short piece. Burri et al.  optimize the squared total duration instead. They soften all constraints by penalizing them in the objective and optimize durations through an NLP solver, while the pieces number is limited.
Gradient-based direct optimization is not satisfactory when computational overhead or solution quality is critical. To address this, Gao et al.  propose a method which decouples geometrical and temporal information of a trajectory. They use QP to generate a spatial trajectory with guaranteed safety in a virtual domain. A temporal trajectory is then optimized through second-order conic programming (SOCP) based on direct collocation , which maps time to the virtual domain. Their method achieves near real-time. The drawback is that the spatial trajectory can restrict the aggressiveness in consequent optimization, resulting in slow motions. In , they improve the above method by alternating minimization between coefficients of the two-layer parametrization. Time durations of the optimized temporal trajectory is used to re-generate a spatial trajectory, on which the optimization is applied again. Although high-quality trajectories can be obtained by this process, several rounds of QP and SOCP make it unapplicable to use online. Moreover, this method lacks theoretical convergence analysis and only works in the rest-to-rest case. Almeida et al. 
propose a machine learning method to train a supervised neural network offline, thus online application is able to refine good initial guesses in real-time. However, the neural network has to be trained case by case.
In this paper, we adopt the time regularized objective  as well as the alternating minimization framework . For efficiency, the spatial and temporal parameters are alternately updated in separate phases. Each minimization phase exploits the objective structure and is solved algebraically, making it free from gradient estimation and step-size choosing. To handle various constraints, we also design a simple yet solid feasibility checker. The proposed framework is able to generate aggressive trajectories at extremely high frequency and not limited to the rest-to-rest case.
Differential flatness of quadrotor dynamics is validated by Mellinger et al. . It means the trajectory planning for a quadrotor can be done solely in the translational space. The kinodynamic feasibility is implicitly transformed into smoothness of the trajectory. Then, actuator constraints can be enforced by restricting norms on high-order derivatives.
In this paper, we employ the piece-wise polynomial trajectory, with each piece denoted as an -order polynomial:
where is the coefficient matrix, is the duration and is a basis function.
It is worth noting that we only consider odd-order polynomial trajectories. Since is odd, the mapping is bijective between the coefficient matrix and the boundary condition. To further explain this, consider derivatives of up to order:
we have where
We denote and by and , respectively. The boundary condition of this polynomial is described by the tuple . The following mapping holds:
where is the mapping matrix. is a non-singular square matrix only if is an odd number. Otherwise, becomes over-determined, which means for any given and , there may not exist polynomial whose satisfies (4).
Moreover, the inverse can be implemented with zero overhead when is an odd number. Burri et al.  explore the structure of and find that can be computed more efficiently through its Schur-Complement, which only involves submatrix inverse. We take things one step further. Actually, all entries of are power functions of , thus Gaussian-Elimination is applied to get its analytic form. Consequently, time-consuming operation is no longer needed when is computed online. To achieve this, we pre-compute matrices offline, where :
Following mapping matrices are computed online:
Therefore, provided with an odd order, we show the practical equivalence between the tuple and in the sense of polynomial representation.
Consequently, we consider an -piece trajectory parametrized by time allocation as well as boundary conditions of all pieces, as shown in Fig. 2. The trajectory is defined by
where lies in the -th piece and is a boundary condition of the -th piece. This definition implicitly involves order continuity at boundaries of each piece. Normally, some entries in are fixed, such as the position of waypoints . We split into two parts, the fixed part which is viewed as constant, and the free part which is to be optimized. Then, the whole trajectory can be fully determined by .
Iv Spatial-Temporal Trajectory Optimization: Unconstrained Case
In this section, we describe our method for jointly optimizing spatial and temporal parameters of a trajectory, for the Unconstrained Case, where no constraint is considered.
Iv-a Optimization Objective
We use the time regularized quadratic cost over the whole trajectory, as the objective of the optimization:
where and are the lowest and the highest order of derivative to be penalized respectively, is the weight of the -order derivative and is the weight of time regularization. The weight adjusts the aggressiveness of the trajectory , which allows total duration varies adaptively. For now, we consider the unconstrained optimization:
where free boundary conditions and durations are decision variables. is used for brevity.
The cost for the -th piece can be calculated as
in which is a symmetric matrix  consisting of high powers of , and is trace operation which only sums up diagonal costs produced in three dimensions. The overall objective can be formulated as
where is the direct sum of its diagonal blocks, and is a permutation matrix. We make sure that the setting for is legal by assuming that the -sublevel set of for any finite is bounded and only consists of positive time allocation. For example, consecutive repeating waypoints with identical boundary conditions fixed in are not allowed.
Iv-B Unconstrained Trajectory Optimization
To optimize Eq. (9), we propose an AM-based algorithm, as shown in Alg. 1. Initially, is solved for the provided . After that, the minimization of the objective function is done through a two-phase process, in which only one of and is optimized while the other is fixed.
In the first phase, the sub-problem
is solved for each . We employ the unconstrained QP formulation by Richter et al. , and briefly introduce it here. The matrix is partitioned as
then the solution is be obtained analytically through
For efficiency, we solve the sparse linear system
through Sparse LU Factorization to get since and are both sparse.
In the second phase, the sub-problem
is solved for each . In this phase, the scale of sub-problem can be reduced into each piece. Due to our representation of trajectory, once is fixed, the boundary conditions isolate each entry in from the others. Therefore, can be optimized individually to get all entries of . As for the -th piece, its cost in (8) is indeed a rational function of . We show the structure of and omit the deduction for brevity:
where and are orders of numerator and denominator, respectively. The coefficient is determined by . Due to positiveness of , we have as or . Therefore, the minimizer exists for
To find all candidates, we compute the derivative of (16):
The minimum exists in solutions of , which can be calculated through any modern univariate polynomial real-roots solver . In this paper, we utilize the Continued Fraction method  to isolate all positive roots of any high order () polynomial. The second phase is completed by updating every entry in .
Iv-C Convergence Analysis
Alg. 1 is globally convergent. Moreover, it is faster than conventional gradient descent used in time allocation refinement, under no assumption on convexity.
Consider the process described in Algorithm 1. Provided with any , the following inequality always holds for the -th iteration:
where and are both constant, is Frobenius norm. It shows the process globally converges to a stationary point with non-asymptotic sublinear rate .
See  for details. ∎
Thm. 1 shows that our algorithm shares the same global convergence rate as that of gradient descent with the best step-size chosen in each iteration . The best step-size is practically unavailable. In contrast, our algorithm does not involve any step-size choosing in its iterations. Sub-problems in Eq. (11) and Eq. (15) both are solved exactly and efficiently due to their algebraic convenience. Therefore, Alg. 1 is faster than gradient-based methods in practice.
Another key advantage of our algorithm is its capability of escaping from a local minimum in the time optimization. Watching Eq. (9), despite is convex in as proved in , it is a rational function which can have multiple local minima in . Therefore, a case may occur where the initial time allocation falls into one of these local minima instead of the global minimum in . Under this situation, naturally, the global minimum in time allocation cannot be attained by gradient-based methods. However, with our method, all local minima are compared directly. Thus, the situation can be avoided.
It is worth noting that, here the global optimality is not guaranteed because our algorithm still exploits local structure of the problem. Although convergence to a stationary point is ensured, we argue that strict saddle points are theoretically and numerically unstable for our first-order AM method . Moreover, when the stationary point is a strict local minimum, we show that the convergence rate is faster.
Let denote any strict local minimum of to which Alg. 1 converges. There exist and , such that
for all , where .
See  for details. ∎
Thm. 2 shows that a faster convergence rate can be attained for a strict local minimum than the general case in Thm. 1. Note that it is possible to accelerate our method to attain the optimal rate of first-order methods or use second-order methods to achieve better performance. However, we still employ the first-order AM process because of its simplicity in implementation and its good performance when the trajectory is far from optimum.
The non-asymptotic property implies that the convergence is bounded strictly by the rate, rather than approximated. This property, along with the monotone decrease of the objective, shows guaranteed progress in each iteration, while gradient-based methods may try bad step-size thus making no/negative progress.
V Spatial-Temporal Trajectory Optimization: Constrained Case
In this section, we present our method to incorporate safety and dynamical feasibility constraints into our optimization process. To begin with, we introduce a computationally efficient feasibility check method that applies to a wide range of constraints. Then this method is used in a constrained trajectory optimization process.
V-a Computationally Efficient Feasibility Check
A piece of the trajectory is denoted by
Constraint over should be satisfied by the -order derivative of the piece. It is required that has the form of a multivariate polynomial:
where is the highest degree. Many kinds of constraints can be expressed by , such as the safe distance constraint to keep away from an obstacle located at :
or maximum speed constraint:
Provided with any piece , we check whether constraint is fulfilled for all . We define which is indeed a polynomial of . The procedure is as follows: Firstly, check the sign of and . Then, If both two endpoints satisfy the constraint, we have to make sure the constraint is not violated inside the interval . Instead of locating all extrema of and checking their values, we only need to check the existence of root of in the interval. If the equation has any root in , then is infeasible. Fortunately, it is convenient for a polynomial to achieve this leveraging Sturm’s Theory . Now that neither nor is the root of , we compute the Sturm sequence by
where is remainder in the Euclidean division of by . When becomes constant, we stop expanding this sequence. Let denote the number of sign variations of Sturm sequence at , in which zero values should be ignored. Then the number of distinct roots inside equals . Here the feasibility check is done for . Sometimes, constraint has the form . In practice, it can be equally handled by checking , where is a small positive real number. What’s more, non-polynomial constraint can also be efficiently checked through its Taylor series within acceptable approximation error.
V-B Constrained Trajectory Optimization
For the Constrained Case, we enforce constraints on norms of derivatives of the trajectory:
Generally, the constraint does not have to be like (23). If only a constraint is representable in (20) and its feasible solution can be constructed, then it can be handled in our optimization. With a slight abuse of notation, we use to denote that is feasible. is used to denote that the -th piece is feasible. We say is tight by meaning that, at least one constraint is tight at a on the -th piece.
The constrained version of our method is shown in Alg. 2. An initial feasible trajectory can be constructed from conservative time allocation. The spatial-temporal parameters are then recovered from the trajectory, which is used in the consequent two-phase constrained minimization.
In the first phase, is fixed. An illustration is provided in Fig. 3(a), where the unconstrained minimum is obtained as is done in Alg. 1. The trajectory may not be feasible. Since the feasibility of is ensured in last iteration, a line search is done as
where is the convex combination of and . Convexity of implies the convexity of . Moreover, is a feasible solution, while is the unconstrained global minimum. We simply take if it is feasible. If not, a bisection procedure is done on the interval . In this procedure, the feasibility check method described in Sec. V-A is employed to shrink the interval. The procedure stops at an acceptable interval length, with taking the feasible lower bound. After that, we update by . Meanwhile, a set is maintained to store indices of tightened pieces.
In the second phase, is fixed. An illustration is given in Fig. 3(b). Each entry in is updated by solving
As stated in Alg. 1, all extrema of can be computed exactly. However, the constrained minimum may not exist in them. It can be any at which some constraints are exactly tightened. When infeasible extremum exists, must be located between any infeasible extremum and the neighboring feasible one or . A bisection procedure with feasibility check suffices to compute . After that, we compare on those feasible extrema together with .
When all iterations are done, the set indicates pieces stuck by active constraints. If is not empty, we recursively apply Alg. 2 on split sub-trajectories, while boundary conditions of pieces indexed by should be totally fixed. Finally, is updated and returned. The recursive process is essential, since it ensures that pieces with no room for optimization do not prevent other pieces to decrease the objective. Alg. 2 is globally convergent to a solution set where constraints are tight or local minimum is attained, which can be checked by Zangwill’s theorem .
Vi-a Comparison of Feasibility Check Methods
Firstly, we compare our feasibility check method with Mueller’s recursive bound check , Burri’s analytical extrema check , as well as the widely used sampling-based check. In each case, 1000 trajectory pieces are randomly generated along with velocity constraints to estimate average time consumption. As is shown in Fig. 4, our method outperforms all other methods in computation speed because of its resolution independence and scalability with higher polynomial orders. The recursive check and sampling-based check may have false positives under rough temporal resolution. The efficiency of analytical check and recursive check deteriorates with higher polynomial orders, because both of them involve roots finding which has closed-form solution only for low order (). In comparison, our method is able to do a solid feasibility check within .
Vi-B Benchmark for Trajectory Optimization Methods
Secondly, we conduct the benchmark comparison of our trajectory optimization method against state-of-the-art methods. The benchmark is done as follows: We generate a sequence of waypoints by random walk, of which the step is uniformly distributed overfor each axis. The maximum speed and acceleration rates are set to and , respectively. The derivatives on the first and the last waypoints are set to zero. The objective function is set as , , and . For a given number of pieces, each method is applied to 1000 sequences of waypoints. The optimization process stops until the relative decrease of objective is less than . The cost is then normalized by the cost of Alg. 1. An illustration in Fig. 5 shows that the minimized cost in unconstrained case can be used as a baseline. All comparisons are conducted on an Intel Core i7-8700 CPU under Linux environment.
We compare our method with Richter’s method  and Mellinger’s method . Richter’s method optimizes trajectory derivatives on waypoints through an unconstrained QP while time allocation is adjusted by gradient descent and scaling. To use it in constrained case, we soften the constraints by penalizing them in objective function as suggested in  and directly optimize the time allocation through NLopt . Mellinger’s method optimizes the time allocation with total duration fixed first, using backtracking gradient descent (BGD). Then dynamical feasibility is ensured through time scaling, of which the ratio is properly calculated using Liu’s method .
, the performance of different methods are provided. Dashed lines indicate standard deviation.
As is shown in Fig. 6(b), our Alg. 2 has the fastest speed and the lowest cost when constraints are taken into consideration. Our method is capable of computing trajectories with pieces within , i.e., at least. However, both benchmarked methods fail to accomplish real-time computing for trajectories with more pieces. Moreover, our Alg. 2 always obtains better trajectories in terms of the cost function, while benchmarked methods cannot fully utilize the capability of system dynamics.
Vi-C Aggressive Flight Experiment
To validate the performance of our method in real-world applications, we deploy it on a self-developed compact quadrotor platform. The proposed method is implemented with C++ 11, and all tasks are conducted on an onboard computer with Intel Core i7-8550U CPU. The pose of our quadrotor is obtained through a robust visual-inertial state estimator . Besides, no external positioning system nor offboard computing is used. A geometric controller is employed for trajectory tracking control .
The experiment is conducted in a complex indoor scene, which is shown in Fig. 1. A globally consistent map of the scene is pre-built, from which some compulsory waypoints are selected offline. At the beginning of each flight, an optimal path is produced by RRT* , which starts from an initial position and passes all compulsory waypoints. Our method generates an optimal trajectory online based on this path within milliseconds. Immediately, the quadrotor starts its aggressive flight. Different from parameters used in benchmark, we set , and . The aggressive flight along with the generated trajectory is shown in Fig. 7. More details are included in the attached video.
In this paper, we propose an efficient trajectory generation method for quadrotor aggressive flight, which has guaranteed convergence and feasibility. Benchmarks for components in our method show its superior computation speed, trajectory quality as well as scalability against state-of-the-art methods. Aggressive flight experiments in limited space with dense obstacles validate the practical performance of our method. Currently, in the proposed framework, positions of waypoints are fixed during optimization. However, the method is underlying compatible with waypoints as part of decision variables. Our feasibility checker also supports various safety constraints. In the future, we plan to apply and improve our method in time-critical large-scale motion planning scenarios where complex spatial constraints exist.
-  D. Mellinger and V. Kumar, “Minimum snap trajectory generation and control for quadrotors,” in Proc. of the IEEE Intl. Conf. on Robot. and Autom., Shanghai, China, May 2011.
-  C. Richter, A. Bry, and N. Roy, “Polynomial trajectory planning for aggressive quadrotor flight in dense indoor environments,” in Proc. of the Intl. Sym. of Robot. Research, Singapore, Dec 2013.
-  M. Burri, H. Oleynikova, M. Achtelik, and R. Siegwart, “Real-time visual-inertial mapping, re-localization and planning onboard mavs in unknown environments,” in Proc. of the IEEE/RSJ Intl. Conf. on Intell. Robots and Syst., Hamburg, Germany, Sep 2015.
-  D. P. Bertsekas, Dynamic Programming and Optimal Control, 1995.
-  A. Beck, First-Order Methods in Optimization, 2017.
-  S. Liu, M. Watterson, K. Mohta, K. Sun, S. Bhattacharya, C. J. Taylor, and V. Kumar, “Planning dynamically feasible trajectories for quadrotors using safe flight corridors in 3-D complex environments,” IEEE Robotics and Automation Letters (RA-L), pp. 1688–1695, 2017.
-  F. Gao, W. Wu, Y. Lin, and S. Shen, “Online safe trajectory generation for quadrotors using fast marching method and Bernstein basis polynomial,” in Proc. of the IEEE Intl. Conf. on Robot. and Autom., Brisbane, Australia, May 2018.
-  L. Campos-Macías, D. Gómez-Gutiérrez, R. Aldana-López, R. de la Guardia, and J. I. Parra-Vilchis, “A hybrid method for online trajectory planning of mobile robots in cluttered environments,” IEEE Robotics and Automation Letters (RA-L), vol. 2, pp. 935–942, 2017.
-  W. Sun, G. Tang, and K. Hauser, “Fast UAV trajectory optimization using bilevel optimization with analytical gradients,” ArXiv, vol. abs/1811.10753, 2018.
-  A. Bry, C. Richter, A. Bachrach, and N. Roy, “Aggressive flight of fixed-wing and quadrotor aircraft in dense indoor environments,” The International Journal of Robotics Research, vol. 34, pp. 1002–969, 2015.
-  F. Gao, W. Wu, J. Pan, B. Zhou, and S. Shen, “Optimal time allocation for quadrotor trajectory generation,” in Proc. of the IEEE/RSJ Intl. Conf. on Intell. Robots and Syst., Madrid, Spain, Oct 2018.
-  D. Verscheure, B. Demeulenaere, J. Swevers, J. De Schutter, and M. Diehl, “Time-optimal path tracking for robots: A convex optimization approach,” IEEE Transactions on Automatic Control, vol. 54, no. 10, pp. 2318–2327, 2009.
-  F. Gao, L. Wang, B. Zhou, L. Han, J. Pan, and S. Shen, “Teach-repeat-replan: A complete and robust system for aggressive flight in complex environments,” ArXiv, vol. abs/1907.00520, 2019.
-  M. M. de Almeida, R. Moghe, and M. R. Akella, “Real-time minimum snap trajectory generation for quadcopters: Algorithm speed-up through machine learning,” in Proc. of the IEEE Intl. Conf. on Robot. and Autom., Montreal, Canada, May 2019.
-  M. Sagraloff and K. Mehlhorn, “Computing real roots of real polynomials,” Journal of Symbolic Computation, vol. 73, pp. 46–86, 2013.
-  E. P. Tsigaridas and I. Z. Emiris, “Univariate polynomial real root isolation: Continued fractions revisited,” in Proc. of the Conf. on Ann. Euro. Sym., Zwitserland, Zurich, Sep 2006.
-  Z. Wang, X. Zhou, C. Xu, and F. Gao, “Detailed proofs of alternating minimization based trajectory generation for quadrotor aggressive flight,” https://arxiv.org/abs/2002.09254.
-  Y. Nesterov, Lectures on Convex Optimization, 2018.
-  J. D. Lee, I. Panageas, G. Piliouras, M. Simchowitz, M. I. Jordan, and B. Recht, “First-order methods almost always avoid strict saddle points,” Mathematical Programming, vol. 176, pp. 311–337, 2019.
-  S. Basu, R. Pollack, and M.-F. Roy, Algorithms in Real Algebraic Geometry, 2003.
-  M. W. Mueller, M. Hehn, and R. D’Andrea, “A computationally efficient motion primitive for quadrocopter trajectory generation,” IEEE Transactions on Robotics, vol. 31, pp. 1294–1310, 2015.
-  W. I. Zangwill, Nonlinear programming: A Unified Approach, 1972.
-  S. G. Johnson, “The NLopt nonlinear-optimization package.”
-  T. Qin, P. Li, and S. Shen, “VINS-Mono: A robust and versatile monocular visual-inertial state estimator,” IEEE Transactions on Robotics, vol. 34, no. 4, pp. 1004–1020, 2018.
-  T. Lee, M. Leok, and N. H. McClamroch, “Geometric tracking control of a quadrotor UAV on SE(3),” in Proc. of the IEEE Control and Decision Conf., Atlanta, Georgia, USA, Dec 2010.
-  S. Karaman and E. Frazzoli, “Sampling-based algorithms for optimal motion planning,” The International Journal of Robotics Research, vol. 30, pp. 846 – 894, 2011.