Real-time Trajectory Generation for Quadrotors using
B-spline based Non-uniform Kinodynamic Search
In this paper, we propose a time-efficient approach to generate safe, smooth and dynamically feasible trajectories for quadrotors in the obstacle-cluttered environment. By using the uniform B-spline to represent trajectories, we transform the trajectory planning to a graph-search problem of B-spline control points in discretized space. Highly strict convex hull property of B-spline is derived to guarantee the dynamical feasibility of the entire trajectory. A novel non-uniform kinodynamic search strategy is adopted, and the step length is dynamically adjusted during the search process according to the Euclidean signed distance field (ESDF), making the trajectory achieve reasonable time-allocation and be away from obstacles. Non-static initial and goal states are allowed, therefore it can be used for online local replanning as well as global planning. Extensive simulation and hardware experiments show that our method achieves higher performance compared with the state-of-the-art method.
Many practical applications of quadrotors require safe navigation in cluttered environments, such as search, rescue, and exploration. Motion planning is the foundation of safe navigation, which generates trajectories that are followed by the motion controller.
Safety, smoothness and good time-allocation are important aspects of evaluating the quality of the trajectory. Smooth and dynamically feasible trajectory allows it to be well tracked by the quadrotor without significant control errors. A trajectory that maintains a safe distance from obstacles as much as possible is not only conducive to avoiding unexpected collision, but also facilitates observation of the environment due to less obstruction by obstacles. Time-allocation, also called time parameterization, affects the trajectory velocity and acceleration in different segments, which in turn affects the safety and efficiency of the execution of the trajectory. And it is also considered as one of the main factors leading to sub-optimality of the trajectory.
In this paper, these aspects are considered simultaneously. Uniform B-spline is utilized to represent trajectory and convert the trajectory generation problem into control point placement. A time-efficient non-uniform kinodynamic search method based on Euclidean signed distance field (ESDF) is proposed, ensuring the safety and reasonable time-allocation, and achieving optimality in time duration and control cost. Main contributions of this paper are summarized as follows:
A uniform B-spline trajectory representation for quadrotors. A highly strict convex hull property is derived for a nonconservative dynamical feasibility check.
A time-efficient non-uniform kinodynamic search algorithm based on ESDF, which generates smooth, safe and reasonably time-allocated trajectories.
Implementation on simulator and hardware, for analyzing and verifying the performance. The source code of this work will be released on https://github.com/tlb9551/BNUKsearch.
Ii Related Works
Trajectory generation for quadrotor can be transformed into the generation of the time-parameterized curve due to its differentially flat . When collision avoidance is considered in cluttered environments, curves such as B-splines , Bézier  and piecewise polynomials [4, 5] are used to represent shape-constrained trajectories.
Many existing planning methods take a two-step pipeline, i.e. a collision-free path is planned at first, then the smoothness and time-allocation of the relevant trajectory are optimized based on the shape of the path. At the front-end, sampling-based [6, 7] and searching-based [8, 9] methods are used to plan a collision-free path. In the back-end, gradient-based methods  and several other methods [3, 11, 12] are employed to guarantee the smoothness and dynamical feasibility. However, these methods separate the trajectory shape and trajectory parameterization, are susceptible to some problems. For example, global optimal trajectory or even feasible trajectory may not be inside the homology class of path which is generated by the front-end methods without dynamic consideration.
Some methods consider both the trajectory shape and the dynamics in parallel. Usenko et al.  used an optimization-based method that integrated the smoothness, dynamic constraints and obstacle-avoidance into the cost function of optimization but faced a low success fraction. Liu et al.  proposed a primitive-based search method, the dynamical feasibility was guaranteed but met the problem of time-efficiency when high-order control inputs were required.
The time-allocation of trajectory is another significant factor that affects the quality of trajectories. Ding et al.  refined the control points of B-spline while the total duration remains fixed. Fernbach et al.  sampled the duration of segments until the related quadratic program is solved, but it may be not suitable to generate trajectories with many segments. Gao et al.  generated time-indexed paths according to the density of obstacle, but the final trajectories do not strictly follow the heuristics after optimization. In this paper, we simultaneously deal with kinodynamic search and time-allocation of the trajectory, which shows noticeable benefits.
Iii B-spline based Trajectory Representation
Since the snap of the trajectory of quadrotor should be continuous , quintic uniform B-spline is selected to represent the trajectory of quadrotors. The B-spline is smooth and has some useful properties, which are suitable for our non-uniform kinodynamic search process.
Iii-a Local Control Property and Kinodynamic Search
Search-based trajectory generation methods need to evaluate the cost of candidate trajectory increments frequently, and the local properties of B-spline curves are tailored to this. The value of a B-spline curve of degree can be evaluated by the following equation:
where are the control points corresponding to and are the blending functions which can be computed by the De Boor-cox recursive formula .
Note that for uniform B-spline, the time interval between is fixed to be , and the is no-zero only in the interval of , which means that a curve span in only depends on local control points , it is called local control property. Therefore, a quintic B-spline curve with control points can be divided into spans (Fig. 2 (a)).
We define as a uniform time representation, as a normalized parameter corresponding to the -th span of the curve , where is the integer part of . Therefore, the B-spline curve can be evaluated by the following matrix representation as mentioned in :
where is the basis vector, is a constant basis matrix of quintic uniform B-spline and is the vector of local control points of -th span.
The derivatives of curve with respect to time can be also expressed as matrix representation. And the derivative curves are also B-spline curves with a lower degree, whose control points are a constant linear combination of higher degree B-spline’s control points.
The integral over squared time derivatives of the -th span, which expresses the control cost of quadrotor trajectory , can be also locally computed in closed form during the search process:
where is a constant matrix and can be computed in advance. The local control property and the matrix representation of enable incremental calculation of the cost during the search process.
One understanding of Eq. 2 is that a span of quintic B-spline can be regarded as weighted combinations of six local control points, and the weight is the time-varying basis, . Fig. 2 (b) shows the weight of each control point that contributes to a span when the uniform time parameter . The weights of and are much higher than other control points, so the span will be close to the connection line between and as shown in Fig. 2 (a). As a result, we can easily confirm that the trajectory mostly has no collision with obstacles as long as the connection line of successive control points does not pass through obstacles. Note that it is not a sufficient condition, and we will introduce the convex hull property to build a sufficient condition in the next subsection.
Iii-B Strict Convex Hull Property and Feasibility Checking
The dynamical feasibility of trajectory guarantees that trajectory can be well tracked by the motion controller while some additional tasks are performed well . As an important property of B-spline, the convex hull property can be used to conservatively evaluate the spans. Existing work is faced with limitations due to its conservativeness, but the novel strict convex hull property we derived here can be used to evaluate the spans of B-spline in a nonconservative way, so as to assess more accuracy in feasibility checking.
The convex hull property of B-spline means that a span of B-spline lies within the convex hull of its local control points , where is the degree of B-spline (see Fig. 2 (c)). Since the derivatives of a B-spline curve are also B-spline curves with a lower degree , they share the same property: local control property and convex hull property. Thus the dynamical feasibility and collision-free constraints can both be expressed as bounds on the placement of control points.
In the previous work , the convex hull property of B-spline is directly used to check the dynamical feasibility, a span is judged to be dynamically infeasible unless all the local control points of its derivative curve are within the border of the dynamically feasible domain, i.e. and . However, since the convex hull of uniform B-spline control points does not tightly wrap the span, this method has some limitations of conservativeness as shown in Fig. 2 (c).
To solve this problem, we try to find a linear transformation of local control points , so that the convex hull formed by more tightly wrap the span. We noticed that the Bézier curve also has convex hull property and the Bézier curve is guaranteed to pass through the first and the last control points, so we transform the local span of B-spline into a Bézier curve representation to achieve tighter convex hull(Fig. 2 (d)). Given:
where is the constant basis matrix of Bézier curve, and are the local control points of Bézier curve. Thus we can use the strict convex hull property for dynamically feasible checking as follows.
Dynamical feasibility checking: firstly, the control points of the derivatives curve of trajectory are evaluated. Then we computed the corresponding Bézier control points with Eq. 5. Finally, if all the Bézier control points are in the dynamically feasible domain, i.e. and , the span of B-spline trajectory is judged to be dynamically feasible (Fig. 2 (d)).
Iii-C Non-static Initial and Goal State
In order to keep the trajectory continuous up to a higher derivative of the quadrotor when the replanning happens in flying, the non-static initial state should be allowed. Given the non-static state of the quadrotor, position and its derivative, as the beginning of trajectory, the first five control points can be obtained by solving the following equation set:
and the 6-th control point is arbitrary.
The goal control points can also be calculated by a similar approach. In our implementation, the goal state is set as a waypoint with specified velocity. Therefore, only two equations related to position and velocity, confirm the last two control points of the whole trajectory.
Iv Non-uniform Kinodynamic Search
As discussed in the previous section, the whole B-spline trajectory can be expressed as a set of B-spline spans. Successive spans share some same control points, i.e. the last control points of span are the first control points of . The increase in control points leads to an increase in the spans and total trajectory duration. Therefore, the trajectory planning is transformed as a problem to find a set of B-spline control points, which generates a collision-free, dynamically feasible and cost-optimized B-spline trajectory between start and goal state. Such a problem can be treated as a graph optimization problem and solved by our search-based algorithm. For convenience, we call it BNUK (B-spline based Non-uniform Kinodynamic) search (Alg. 1).
Iv-a Algorithm Overview
The graph structure is chosen as follows. Vertices is the cell center of voxels in discretized 3-D space, indicating control points, and edges is their connections. Each vertex is connected to its neighbors in 26 directions like the usual A* search algorithm in 3-D space but obeys the rule of choosing specific step length between neighbors, which will be exhaustively discussed in the next subsection.
Thus the kinodynamic search problem is to find an ordered set of control points from the graph , which satisfy the dynamical feasibility, and minimizes the cost of trajectory as given below,
which is a weight of control cost and trajectory duration. The cost can be evaluated incrementally due to the local control property.
Iv-B Non-uniform Node Expansion
Node expansion refers to finding neighboring nodes of the current node and is implemented in function in Alg. 1.
Different from traditional path finding as A* search, the output of BNUK search, i.e. a set of control points is not required to be connected in 3-D discrete space. Therefore, we propose a novel non-uniform node expanding strategy. The step length of the node expansion, i.e. the spacing between two control points, is dynamically adjusted during the search process.
As discussed in Sec. III, the spacing length between control points will affect the derivatives of the trajectory due to the fixed time interval . We argue that it is reasonable to reduce the velocity of quadrotors in the vicinity of the obstacles for security consideration, because the actual trajectory is always deviated from the planned trajectory due to the existence of control errors. Similarly, the velocity should be increased to take full advantage of quadrotor’s dynamics when away from obstacles. Therefore the spacing length between control points needs to be set in reference to some quantity describing the distance from the obstacle.
The Euclidean signed distance field (ESDF) is employed as the reference. The cells of ESDF store the Euclidean distance to the nearest obstacle. We use the following function to mapping distance to step length of node expansion:
where is a parameter related to the radius of the quadrotor. This function is limited under , so it does not make the spacing between control points too large to cause the trajectory exceeding the dynamical velocity limit. And step length is proportional to the distance and is equal to 0 when the distance is under , meaning that the node is always expanding in the collision-free area. This avoids time-consuming obstacle inflation and also improves safety.
The non-uniform node expansion benefits the search process in three aspects. Firstly, safety is improved because the step length is always shorter than the distance to the nearest obstacle, no control point will be set in non-free or dangerous area. Secondly, since the trajectory duration is a part of cost function, equal to , expanding to a node with large step length is less costly than a node with small step length. Thus the trajectory will be naturally away from obstacles to achieve less trajectory duration. Thirdly, trajectory is reasonably time-allocated, slowly in obstacle-dense area and fast in obstacle-sparse areas.
Iv-C Dynamical Feasibility Checking and Heuristic Function
The dynamical checking and cost evaluating both need the value of local control points, thus we utilize the function to get the control points of the local span. Then dynamical checking and cost evaluating are done by the method discussed in Sec.III. As for the collision-free checking, it is time-consuming and unnecessary so we leave it out of the search process.
The heuristic function is useful to speed up the search process due to the reduction of expanded node number, it is an estimation of cost from the current node to the goal:
where denotes the diagonal distance. This heuristic function is admissible since it underestimates the optimal cost from to , and is efficient in speeding up the search process.
V-a Comparison with State-of-the-Art Method
In this section, we compared our BNUK search with a state-of-the-art method proposed by Gao et al.  in two ways, trajectory generation on random maps and navigation simulation on large maps.
|Gao’s method ||12.39||11.94||68.56||1.04||0.29||0.54||1.63|
on 50 Random Maps
Gao’s method uses the fast marching in the velocity field transformed from ESDF to find a time-indexed path, the trajectories represented by piecewise Bézier curves are generated by solving a convex quadratic program. We compared our method with Gao’s method because both methods use the obstacle density information in EDSF to generate trajectories with reasonable time-allocation, but in different ways and trajectory representation, resulting in different performance in time-efficiency and trajectory quality.
We generate 50 random maps and plan trajectories from the center to side and corner with two methods, simulating local replanning that often occurs in autonomous navigation. The map size is and resolution is . The maximum velocity and acceleration are set as and . As the result shown in Tab. I, both methods generate dynamically feasible trajectories and behave similarly in terms of trajectory length and duration. The mean and minimum distance from obstacles of our methods are larger than Gao’s method, meaning that the trajectory generated by our methods is safer. In addition, our method is much more time-efficient, and the average calculation time is less than a quarter of theirs.
We also simulate the navigation in several large random forest maps with two methods. RHC (Receding Horizon Control)  is used as the same strategy. Replanning happens when the quadrotor arrives at the border of prior sensing horizon or when a collision is detected. The map size is and resolution is , the local sensing map size is . We set the same maximum velocity and acceleration as and , same start and goal are set so that the trajectory length is longer than . As shown in Fig. 4 (d), our trajectory (red curve) is safer than Gao’s (green curve). Tab. II shows the average data of two methods, the trajectory length is similar but our method reaches a lower duration. Our method is much more time-efficient. The replanning times of Gao’s method are higher than ours because it encountered failures during the trajectory optimization process.
|Gao’s method ||51.75||56.94||194.27||26.3|
V-B Flight Experiment
Flight experiment is done to demonstrate the collision avoidance and real-time calculation of our method. We use a quadrotor based on the Pixhawk flight controller, and an UP2 embedded board is utilized as the onboard computer. A motion capture system is used to locate the quadrotor and obstacles. The map size of our indoor flight is , and the resolution is . The maximum velocity and acceleration are set as and . We employ the controller in  to track the trajectory generated by our method. Fig. 6 shows the snapshots of our experiment. The time of trajectory generation is . The max velocity, acceleration, and jerk of trajectory are , and , thus the whole trajectory is dynamically feasible and smooth. The minimum distance to the obstacle is while no obstacle inflation is applied, which means the trajectory generated by our method is quite safe. We recommend readers to watch the attached video for more details.
In this paper, we propose a time-efficient search-based approach to generate safe, smooth and dynamically feasible trajectories for quadrotors. Benefiting from ESDF-based non-uniform search, our algorithm generates trajectories away from obstacles and has reasonable time-allocation. The highly strict convex hull property we derived makes the trajectory dynamically feasible and not too conservative. We have also enabled non-zero initial and termination states to be used so that algorithms can be applied in different scenarios. Extensions to this work may include developing informative heuristics for our search process and exploring the application of the non-uniform rational B-spline curve in trajectory planning.
-  D. Mellinger and V. Kumar, “Minimum snap trajectory generation and control for quadrotors,” in Proc. IEEE Int. Conf. Robotics and Automation, May 2011, pp. 2520–2525.
-  E. Koyuncu and G. Inalhan, “A probabilistic b-spline motion planning algorithm for unmanned helicopters flying in dense 3D environments,” in Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems, Sept. 2008, pp. 815–821.
-  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. IEEE Int. Conf. Robotics and Automation (ICRA), May 2018, pp. 344–351.
-  C. Richter, A. Bry, and N. Roy, “Polynomial trajectory planning for aggressive quadrotor flight in dense indoor environments,” in Robotics Research. Springer, 2016, pp. 649–666.
-  L. Zhao, H. Wang, and W. Chen, “Optimal trajectory planning for manipulators with flexible curved links,” in International Conference on Intelligent Autonomous Systems. Springer, 2016, pp. 1013–1025.
-  S. Karaman and E. Frazzoli, “Sampling-based algorithms for optimal motion planning,” The international journal of robotics research, vol. 30, no. 7, pp. 846–894, 2011.
-  Z. Wei, W. Chen, H. Wang, and J. Wang, “Manipulator motion planning using flexible obstacle avoidance based on model learning,” International Journal of Advanced Robotic Systems, vol. 14, no. 3, p. 1729881417703930, 2017.
-  M. Likhachev, G. J. Gordon, and S. Thrun, “Ara*: Anytime a* with provable bounds on sub-optimality,” in Advances in neural information processing systems, 2004, pp. 767–774.
-  Q. H. Do, L. Han, H. Tehrani Nik Nejad, and S. Mita, “Safe path planning among multi obstacles,” in Proc. IEEE Intelligent Vehicles Symp. (IV), June 2011, pp. 332–338.
-  N. Ratliff, M. Zucker, J. A. Bagnell, and S. Srinivasa, “Chomp: Gradient optimization techniques for efficient motion planning,” in Proc. IEEE Int. Conf. Robotics and Automation, May 2009, pp. 489–494.
-  J. Chen, T. Liu, and S. Shen, “Online generation of collision-free trajectories for quadrotor flight in unknown cluttered environments,” in 2016 IEEE International Conference on Robotics and Automation (ICRA). IEEE, 2016, pp. 1476–1483.
-  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, vol. 2, no. 3, pp. 1688–1695, July 2017.
-  W. Ding, W. Gao, K. Wang, and S. Shen, “Trajectory replanning for quadrotors using kinodynamic search and elastic optimization,” in Proc. IEEE Int. Conf. Robotics and Automation (ICRA), May 2018, pp. 7595–7602.
-  V. Usenko, L. von Stumberg, A. Pangercic, and D. Cremers, “Real-time trajectory replanning for mavs using uniform b-splines and a 3D circular buffer,” in Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems (IROS), Sept. 2017, pp. 215–222.
-  S. Liu, N. Atanasov, K. Mohta, and V. Kumar, “Search-based motion planning for quadrotors using linear quadratic minimum time control,” in 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). IEEE, 2017, pp. 2872–2879.
-  P. Fernbach, S. Tonneau, and M. Taïx, “Croc: Convex resolution of centroidal dynamics trajectories to provide a feasibility criterion for the multi contact planning problem,” in Proc. IEEE/RSJ Int. Conf. Intelligent Robots and Systems (IROS), Oct. 2018, pp. 1–9.
-  L. Piegl and W. Tiller, The NURBS book. Springer Science & Business Media, 2012.
-  K. Qin, “General matrix representations for b-splines,” The Visual Computer, vol. 16, no. 3, pp. 177–186, 2000.
-  H. Wang, Y. Lai, and W. Chen, “The time optimal trajectory planning with limitation of operating task for the tokamak inspecting manipulator,” Fusion Engineering and Design, vol. 113, pp. 57–65, 2016.
-  J. Bellingham, A. Richards, and J. P. How, “Receding horizon control of autonomous aerial vehicles,” in Proc. American Control Conf. (IEEE Cat. No.CH37301), vol. 5, May 2002, pp. 3741–3746 vol.5.
-  T. Lee, M. Leok, and N. H. McClamroch, “Geometric tracking control of a quadrotor UAV on se(3),” in Proc. 49th IEEE Conf. Decision and Control (CDC), Dec. 2010, pp. 5420–5425.