{'fraction_non_alphanumeric': 0.03443577356620835, 'fraction_numerical': 0.0063173541434411, 'mean_word_length': 4.606481481481482, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
STD-Trees: Spatio-temporal Deformable Trees for Multirotors Kinodynamic Planning

In constrained solution spaces with a huge number of homotopy classes, stand-alone sampling-based kinodynamic planners suffer low efficiency in convergence. Local optimization is integrated to alleviate this problem. In this paper, we propose to thrive the trajectory tree growing by optimizing the tree in the forms of deformation units, and each unit contains one tree node and all the edges connecting it. The deformation proceeds both spatially and temporally by optimizing the node state and edge time durations efficiently. The unit only changes the tree locally yet improves the overall quality of a corresponding sub-tree. Further, variants to deform different tree parts considering the computation burden and optimizing level are studied and compared, all showing much faster convergence. The proposed deformation is compatible with different RRT-based kinodynamic planning methods, and numerical experiments show that integrating the spatio-temporal deformation greatly accelerates the convergence and outperforms the spatial-only deformation.


I. INTRODUCTION
Kinodynamic planning [1] respects directly the system kinematics and the dynamical constraints in the process of globally searching for an optimal trajectory in exploring the entire constrained solution space and bring an optimal trajectory. Although much progress has been made, it remains challenging to globally plan a near-optimal trajectory efficiently in large-scale environments with complex obstacles.
For multirotor systems of high dimensional states, sampling-based variants [2]- [4] are popular by holding the anytime property and being asymptotically optimal. They connect state samples in the discretized state space and grow trajectory trees to explore the constrained solution space. By introducing some randomness, they are naturally suitable to search into spaces near different local minima and potentially bring better solutions as more computation resources are available. Although the randomness necessities optimality, it impedes the convergence since the probability to sample exactly near the optimal solution is extremely low, especially in high dimensional solution spaces where a huge number of homotopy classes exist.
We embrace local optimization to alleviate this problem. The optimization is performed by locally deform the trajectory tree and keep it well grown, both spatially and temporally, as shown in Fig. 1. Tree deforming is studied by Hauer and Tsiotras [5] in the geometric path planning problem. They reposition the locations of nodes to minimize (a) With ST Deformation.
(b) Without deformation. Fig. 1: Current tree and best trajectory after adding a same number of nodes. Planning with deformation acquires a better-organized tree structure with overall higher quality and gets a smoother trajectory.
the overall tree cost and improve the convergence. However, they do not account for trajectory planning. Focused on kinodynamic planning problems where an extra time dimension needs to be considered, we show that deforming the trajectory tree in additional the time dimension, that is, spatiotemporally (ST), can further facilitate the convergence. The difficulties, however, appear in designing practical objectives and attaining efficiency. Therefore, we propose to deform the tree in forms of deformation units to balance the optimization level and the computation burden.
In this work, we focus on accelerating the convergence in sampling-based kinodynamic planning by optimizing specific node states and the time durations of some related edges. We summarize the main contributions as follows: 1) We propose a novel sampling-based kinodynamic planner for multirotors. It thrives the trajectory tree growing by optimizing the tree in forms of deformation unit that contains one tree node and the edges connecting to it. One deformation unit only constitutes a small part of the entire tree, yet deforming it improves the overall quality of the sub-tree rooted at the deformed node. 2) We propose to efficiently refine the deformation unit in not only shapes but also the time dimension. It is compatible with different RRT-based kinodynamic planners, and integrating it shows a faster convergence and generates much better final trajectories. 3) We propose several variants to perform the deformation by different levels of computation burden. Numerical results show that all the variants converge quicker than the baseline and overtake some search-based planners.
II. RELATED WORK Basically, different sampling-based algorithms differ in two procedures, the exploration of solution space and the exploitation of information brought by new samples, resulting in different convergence. The original RRT algorithm [6] simply connects a state node to the nearest node in the current tree and is very fast in exploration to get solutions. It is, however, proved to converge to a suboptimal solution by Karaman and Frazzoli in [7] where the RGG and RRT* algorithms are introduced. They stand out in the process of choosing parent nodes and rewire local tree structure according to the best cost-from-start values and can converge to the optimal with probability one. The RRT # [8] algorithm resolves the under-exploitation in RRT* and improves convergence by further rewiring in cascade to propagate the new information to other parts of the tree and maintain a well-grown tree of consistent nodes. We adopt the RRT # scheme to balance the exploration-exploitation tradeoff in our kinodynamic algorithm framework. To acquire higher convergence rate, recent works focus on improving sampling strategies by sampling directly in the Informed Set [9], the Relevant Region Set [10], or the Local Subset of GuILD [11]. However, these designs depend on the L2-Norm objective. For applications with varied desired objectives, analytic forms of the sampling set can not be obtained.
Combining local optimization in the global searching process is studied by many for geometric shortest path planning. RABIT* [12] adopts CHOMP [13] to find feasible alternative connections for edges in collision and thus accelerate the search. However, the convergence rate highly depends on shrinking the L2-Norm based informed set in their BIT* [14]-like batch sampling and searching process. The rather heavy optimization can also hinder exploring different homotopy classes. DRRT [5] does not require any informed sampling set. It minimizes the overall tree cost by optimizing the locations of some existed nodes in the tree. A well-grown tree is maintained to cover the search space in each iteration, and thus paths obtained are of high quality and the convergence is improved. For kinodynamic planning, INSAT [15] combines trajectory optimization into discretized graph search. However, the solving time is intractable because of high-dimensional graph search, and long trajectories are being optimized many times in each node expanding step. Although weighted A* [16,17] is used to accelerate searching, it eventually generates a sub-optimal solution since the admissibility criterion is relaxed. In our combining scheme, the sampling-based search globally reasons for the optimal solution with asymptotic optimality, and the optimization only deforms limited edges of the tree.

III. PROBLEM STATEMENT
The multirotors kinodynamic planning problem is to find a trajectory that minimizes some cost and satisfies some constraints. Following [3,4,18,19], the cost minimized is a trade-off between time and energy, which suits many applications. The constraints imposed include obstacle avoidance, system dynamics, kinematics and start and goal constraints. According to the multirotor systems' differential flatness property [20], we can use a linear model of chain integrator x new ← Sampling(E) 5: One Solution Found. 10: end if 11: U ← SelectDeformationUnits(x n , L) 12: DeformInOrder(U) 13: RewireInCascade(T , x new ) 14: end while 15: return T to represent its dynamics with four flat outputs p x , p y , p z (position in each axis), ψ (yaw), and their derivatives being the state variables. Thus, the trajectory planning problem is formulated as a Linear Quadratic Minimum Time (LQMT) problem in the area of optimal control. For each axis (yaw is excluded in this paper), it is formulated as follows: where τ is the time duration of the trajectory, ρ the tradeoff weight to penalize time against energy, x start the initial state, x goal the goal state, and G 0 the obstacle avoidance and higher derivative limit constraints.

IV. KINODYNAMIC DEFORMABLE TREES
The nonlinear constraints make the entire solution space highly non-convex. Directly applying nonlinear programming will most definitely fall into local minima. We adopt a sampling-based approach as described in Alg. 1 to address the problem. Basically, the algorithm is based on kinodynamic RRT* [3] with features of RRT # [8] and with an extra step to deform the shape and the time duration of the current trajectory tree (line 12) in each iteration. The searching framework breaks the problem down into smaller subproblems by incrementally adding state samples to probe the solution space. Each subproblem is of the same form as Equ.1 with different boundary states, and we get a trajectory segment by solving one.
In each iteration, after a state node is sampled, we seek its best parent node in the BackwardNear node set X backward according to a cost-from-start value estimated by solving the LQMT subproblems. If a qualified best parent is found, the sampled node and a new edge are added to the tree. It then tries to connect to the goal node for a possible solution. After that, the core part, the tree deformation is performed by optimizing some trajectory edges, both the profile and the time duration. This step reduces the overall tree cost without adding more nodes. Since spending more time on exploitation means less exploration, the deformation is only activated after finding a first feasible solution. Finally, since the structure of the implicit graph is changed while adding nodes, we rewire some nodes for potential improvement of their cost-from-start values. The rewire performs in cascade to propagate the new information brought by the new node to other parts of the entire tree.
The final trajectory tree consists of many trajectory segments as its edges, and the optimal trajectory between start and goal consists of a chain of trajectory segments between a series of successive states.

V. TREE EDGE REPRESENTATION
As seen, solving the LQMT subproblems makes up a fundamental part of the searching process and will be called many times in each iteration. Therefore, it needs to be solved fast. Suppose the integrator model is of s th -order, which leads the state variable to be x(t) = [p(t), · · · , p (s−1) (t)] T and the control to be u(t) = p (s) (t). According to [3,21], if all the inequality constraints are temporally disregarded, the optimal control law for the unconstrained LQMT problems are parameterized by polynomials of degree s − 1, and then each state variable in x(t) can be obtained by integration up to s times. Thus, we have p(t) = c T β(t), t ∈ [0, T ], where c ∈ R 2s is the polynomial coefficient vector, T the time duration, and β(t) = (1, t, t 2 , · · · , t 2s−1 ) T the natural basis.
This leads to a conclusion that given time durations T and boundary conditions d = [x T (t)| t=0 , x T (t)| t=T ] T ∈ R 2s , the optimal solutions for these specific LQMT problems are fully determined since no freedom is left. As claimed by our previous work [22], a smooth bijection exists between the polynomial coefficient vector c and the boundary condition vector d, that is, an analytic transformation exist between the two kinds of descriptions of polynomials {c, T } and {d, T }: where A f (T ) and A b (T ) are forward and backward mapping matrices whose entries are analytically determined. The cost of one edge can then be calculated by {c, T } as or by {d, T } as both analytically. The previously ignored constraints are checked afterwards. If satisfied, a tree edge is grown, and the cost provides an estimation of the optimal transition cost between two state samples in the constrained solution space. The two kinds of representation of trajectory edges support the following tree deformation process.

VI. TRAJECTORY TREE DEFORMATION
By the end of each iteration of the tree growing process, some new information of the solution space is collected in the tree. However, the cost-from-start value of some nodes may not well estimate the actual value due to insufficient sampling, presenting that the tree does not grow well, especially in the time dimension. A deforming of the edges can improve its quality, as depicted in Fig. 2.

A. Deformation Unit
Given current tree structure T , each node n ∈ T contains the following information: • x n : the corresponding state of n, • p n : the parent node of n in T , • T n : the time duration of the edge from p n to n, • c n : the coefficient vector of the edge from p n to n, • C n : the children node set of n in T , • g n : the cost-from-start value of n following T . If we reposition one of the nodes and fix all others, meanwhile keeping the connection between nodes unchanged, then only the edges that connect to the repositioned node will be affected. Thus, we set a deformation unit by the collection of one node, the edge from its parent, and the edges to its children. Within a unit, the deformation is performed by optimizing the node state and the time duration of these edges, while all other nodes' state and edge time durations are fixed. We denote the element collection in a deformation unit as {x n , T n } = {x n , T n , T i , |i ∈ C n }, where n is the node being deformed, and T n = {T n , T i , |i ∈ C n } contains the time durations of all the edges connecting to n. By setting the node state as decision variables instead of edge coefficients, the equality constraints imposed by the continuity requirement in-between tree nodes are implicitly eliminated, and the number of decision variables is reduced.

B. Objective Design
When optimizing the deformation unit, we need an appropriate objective. It should follow that after the deformation, the overall tree quality is improved and it stands a better chance to find a better solution or the optimal solution. As revealed by Hauer et al. [5], the sum of the cost-from-start values of a fixed number of sample nodes stands for an estimation of the optimal value function for the solution space covered by the selected sample set. Considering also that optimizing the deformation unit will and only will change the cost-from-start values of node n and all its descendant nodes, leaving all other nodes unaffected. We therefore set the objective as where D n is the set of all n's descendant nodes. Optionally, we can assign a weight w ∈ [0, 1] for each node suggesting the possibility of it constituting the final best trajectory. A heuristic estimation of the cost-to-go can be an alternative for this weight. The objective suggests an optimization of the sub-tree rooted at node n. By growing of the tree, we have where c(x i , x j ) is the edge cost connecting node i to node j, analytically estimated by Equ. 4. We can then use the weighted sum of all the edge costs in the tree to calculate the objective. The weight of an edge (i, j) is assigned by counting the total number of paths that start from the start node to any other node in the sub-tree using this edge. The objective thus becomes where T n is the sub-tree rooted at node n, d n equals 1 + nb des(n) with nb des(n) the number of descendants of node n, and C a constant indicating the sum of all other edges' cost in the tree, which is irrelevant to the elements in the deformation unit. Therefore, the objective depends only on the edges in the deformation unit yet expresses the quality of a sub-tree that estimates the value function of a part of the solution space.

C. Unconstrained Formulation
In the process of tree growing, we ignore the constraints and check afterwards to achieve faster exploration of the entire solution space. When deforming, it is wished to focus more on the local space covered by the deformation unit and exploit the information already gathered in the tree. The constraints on obstacle avoidance and dynamical feasibility should now be considered.
Denote G(p [s] (t)) ∈ R s+1 0 the functional-type constraints that consider obstacle avoidance and dynamical limitations such as the velocity, acceleration and other higher order derivative constraints up to order s. To eliminate these inequality constraints, we construct penalty functions and turn them into soft ones. Since they consist of an infinite number of inequality constraints and can hardly be directly handled, we transform them into finite-dimensional ones via integral of constraint violations, and the integral is estimated by weighted sum of the sampled penalty functions. Note that in practice, these constraints are decoupled between edges, that is, G(p [s] (t)) with t ∈ [0, T i ] are solely determined by c i and T i . Thus, for one trajectory edge parameterized by {c i , T i }, we compute the penalty function as where k i is the sample number on this edge, X ∈ R s+1 ≥0 is a vector of penalty weights for each entry of G, (ω 0 , ω 1 , · · · , ω ki−1 , ω ki ) = (1/2, 1, · · · , 1, 1/2) are the quadrature coefficients following the trapezoidal rule, max[ ] is an element-wise comparison, and t = j/k i · T i is the sampled time stamp. Integrating the penalty into the objective by each edge and the unconstrained optimization of one deformation unit is formulated as min xn,Tn i∈{n} Cn with {x n , T n } being the decision variables.

D. Spatio-temporal Optimization
For one edge {c i , T i }, i ∈ {n} C n in the deformation unit, we derive the gradient of the decoupled objective w.r.t {x n , T n } by chain rule as follows: where is the Hadamard product, Sign[ ] is element-wise indication of the sign of a number, and A xx b (T ) contains block elements of A b (T ) partitioned by s. The gradients of an edge w.r.t. the time durations of other edges are all 0.
The conditional computation implies different cases that x n being tail or head state of an edge.
The only parts undefined in the above equations are ∂G/∂c i and ∂G/∂t. For obstacle avoidance, we wish the edge positions have certain clearance and thus define where r > 0 is the preferred distance away from obstacles, and F(p(t)) computes a minimum distance to the closest obstacle given a position in the edge, which can be precomputed with a distance field built incrementally [23,24] or in batch [25]. For dynamical feasibility, we limit the amplitude of higher-order derivatives and define where m k is the maximum allowed value of the k th derivative. The corresponding gradients are where ∂F/∂p(t) is computed with interpolation of the distance field. The analytic transformation between {c i , T i } and {x pi , x i , T i } provides efficiency for calculating the objective and the gradients since no matrix inversion is required. With the gradients at hand, numerical optimization is applied then to solve this optimization problem. Considering that the interpolation of the distance fields F introduces nonsmoothness, Newton-type methods like BFGS can be inaccurate or fail sometimes. We thus adopt bundle methods to address it. In this work, the Limited Memory Bundle Method (LMBM) [26] is used.

E. Deformation Variants
The optimization of one deformation unit can explore its surrounding local solution space more thoroughly while other parts of the tree are left unchanged. After each iteration of Alg. 1 though, the current tree topology in different parts of the tree can be changed by adding new nodes and by RewireInCascade. Applying deformation to more tree parts can improve the overall tree quality but may at the cost of more computation loads. Weighing the optimization levels and the computation burden, we propose four kinds of variants. As depicted in Fig. 3, after one node new is newly added to the tree, denoting its parent node in the tree by n, the tree deformation is performed as one of the following variants.
1) NODE: optimizes only one deformation unit which contains n and the edges connecting to it. 2) TRUNK: optimizes several deformation units, and the units are selected by following parent pointers from n up to a child of start. The deformation is then performed in an order from the child node to n. This is the strategy DRRT [5] adopts. 3) BRANCH: optimizes every nodes and edges in the sub-tree rooted at n. The order follows breadth-first search starting from n. All the leaf nodes are excluded. 4) TREE: optimizes every node and edge in the entire tree except the start node and all the leaf nodes. The order follows breadth-first search starting from the children nodes of start.

A. Experiment Settings
For numerical comparisons, we set s = 3, which means a third-order integrator is used to model our multirotor system. The weight of time ρ is set 100. The dynamical limitations are set as 5m/s for velocity, 7m/s 2 for acceleration, and 15m/s 3 for jerk. For collision checking, the maps are prebuilt and transformed into occupancy grids of 0.1m resolution with obstacles inflated by 0.2m. All the numerical experiments are conducted with a desktop computer with a 3.4GHz Intel i7-6700 processor.

B. Variants Comparison
In this experiment, the proposed planning method w/o deformation or with one of the four deformation variants introduced in Sec. VI-E are benchmarked with each other and with search-based kinodynamic planners. For search-based ones, the weighted version of A* [17] is used to accelerate the search. Three different weights are used to amplify the heuristic: 1.7, 2.3, and 2.8. Each variant runs for 100 trials with a 3 second time budget. The environment and final trajectories of one trial are shown in Fig. 4. As Fig. 6 shows, the proposed methods find the first solution within milliseconds and then quickly converge. Compared to the baseline (w/o), the convergence rate all improves remarkably by planning with any variants of the proposed spatio-temporal deformation. Among the variants, BRANCH slightly beats TREE and much outperforms NODE and TRUNK. We think it is because a new sample brings potential improvements mostly on the sub-tree while other tree parts are less likely influenced, and thus deforming just the sub-tree is a good balance on computation cost and tree quality improvement. The search-based ones with different weights, however, take much longer time to get a solution with higher cost.

C. Deformation Comparison
For three sampling-based kinodynamic planning methods, kRRT, kRRT*, and kRRT # , we compare the proposed spatiotemporal deformation with deforming only spatially and   with no deformation, denoting a suffix of -ST, -S, and no suffix, respectively. For each deformation, variant BRANCH is adopted. Each method runs for 100 trials with a 7 second time budget. The environment and final trajectories of one trial are shown in Fig. 5. Fig. 7 shows that for all the RRTbased methods, both deformations accelerate the convergence evidently. The proposed spatio-temporal deformation further considerably outperforms the spatial-only deforming .
VIII. CONCLUSION In this paper, we propose a sampling-based kinodynamic planning method for multirotors combining local optimiza- tion by deforming some selected tree edges spatially and temporally. The optimization is performed efficiently within deformation units which contain the state of one node and the time durations of the edges connecting it. Though only limited parts are deformed, the overall quality of a sub-tree is improved and a well-grown tree is maintained after each node is added to the tree. By deforming the state of some selected nodes and related edge durations without adding more nodes, the convergence is improved. Benchmark results show that integrating the proposed deformation achieves a much faster convergence rate. We open source our code for the reference of the community.