---
abstract: 'In this paper, we propose a complete and robust motion planning system for the aggressive flight of autonomous quadrotors. The proposed method is built upon on a classical teach-and-repeat framework, which is widely adopted in infrastructure inspection, aerial transportation, and search-and-rescue. For these applications, human’s intention is essential to decide the topological structure of the flight trajectory of the drone. However, poor teaching trajectories and changing environments prevent a simple teach-and-repeat system from being applied flexibly and robustly. In this paper, instead of commanding the drone to precisely follow a teaching trajectory, we propose a method to automatically convert a human-piloted trajectory, which can be arbitrarily jerky, to a topologically equivalent one. The generated trajectory is guaranteed to be smooth, safe, and kinodynamically feasible, with a human preferable aggressiveness. Also, to avoid unmapped or dynamic obstacles during flights, a sliding-windowed local perception and re-planning method are introduced to our system, to generate safe local trajectories onboard. We name our system as *teach-repeat-replan*. It can capture users’ intention of a flight mission, convert an arbitrarily jerky teaching path to a smooth repeating trajectory, and generate safe local re-plans to avoid unmapped or moving obstacles. The proposed planning system is integrated into a complete autonomous quadrotor with global and local perception and localization sub-modules. Our system is validated by performing aggressive flights in challenging indoor/outdoor environments. We release all components in our quadrotor system as open-source ros-packages[^1].'
author:
- 'Fei Gao, Luqi Wang$^*$, Boyu Zhou$^*$, Luxin Han, Jie Pan, and Shaojie Shen[^2] [^3]'
bibliography:
- 'tro2019fei.bib'
title: |
    Teach-Repeat-Replan:\
    A Complete and Robust System for Aggressive Flight in Complex Environments
---

Aerial Systems: Applications, Motion and Path Planning, Collision Avoidance, Autonomous Vehicle Navigation.

Introduction
============

the development of autonomy in aerial robots, Micro Aerial Vehicle (MAV) has been more and more involved in our daily life. Among all applications emerged in recent years, quadrotor teach-and-repeat has shown significant potentials in aerial videography, industrial inspection, and human-robot interaction. In this paper, we investigate and answer the problem of what is the best way to incorporate a human’s intention in autonomous and aggressive flight, and what is a flexible, robust and complete aerial teach-and-repeat system.

There is a massive market for consumer drones nowadays. However, we observe that most of the operators of consumer drones are not professional pilots and would struggle in generating their ideal trajectory for a long time. In some scenarios, such as the drone racing or aerial filming, a beginner-level pilot is impossible to control the drone to finish the race safely or take an aerial video smoothly unless months of training. Also, there is considerable demand in applying drones to repetitive industrial inspections or search-and-rescue missions, where human provides a preferable routine. In these situations, demonstrating a desirable trajectory and letting the drone to repeat it is a common wish. However, the taught trajectory generated by an unskilled pilot is usually incredibly hard or dynamically infeasible to repeat, especially in some cluttered environments. Moreover, most of the vision-based teach-and-repeat applications [@fei2019ral], [@fehr2018visual], [@furgale2010visual], such as our previous work [@fei2019ral], are sensitive to changing environments. In [@fei2019ral], even the environment changes very slightly, the global map has to be rebuilt, and the teaching has to be redone.

Based on these observations, in this paper, instead of asking the drone to follow the human-piloted trajectory exactly, we only require the human operator to provide a rough trajectory with an expected topological structure. Such a human’s teaching trajectory can be arbitrarily slow or jerky, but it captures the rough route the drone is expected to fly. Our system then autonomously converts this poor teaching trajectory to a topological equivalent and energy-time efficient one with an expected aggressiveness. Moreover, during the repeating flight, our system locally observes environmental changes and re-plans sliding-windowed safe trajectories to avoid unmapped or moving obstacles. In this way, our system can deal with changing environments. Our proposed system extends the classical robotics teach-and-repeat framework and is named as *teach-repeat-replan*. It is complete, flexible, and robust.

In our proposed system, the surrounding environment is reconstructed by onboard sensors. Then the user’s demonstrated trajectory is recorded by virtually controlling the drone in the map with a joystick or remote controller. Afterward, we find a flight corridor that preserves the topological structure of the teaching trajectory. The global planning is decoupled as spatial and temporal planning sub-problems. Having the flight corridor, an energy-optimal spatial trajectory which is guaranteed to be safe, and a time-optimal temporal trajectory which is guaranteed to be physically feasible, are iteratively generated. In repeating, while the quadrotor is tracking the global spatial-temporal trajectory, a computationally efficient local map [@han2019fiesta] is fused onboard by stereo cameras. Based on local observations, our proposed system uses a sliding-window fast re-planning method [@boyu2019ral] to avoid possible collisions. The re-planning module utilizes gradient information to locally wrap the global trajectory to generate safe and kinodynamic feasible local plans against unmapped or moving obstacles.

The concept of generating optimal topology-equivalent trajectories for quadrotor teach-and-repeat was first proposed in our previous research [@fei2019ral]. In [@fei2019ral], once the repeating trajectory is generated, the drone executes it without any other considerations. In that work [@fei2019ral], the environment must remain intact during the repeating, and the localization of the drone is assumed to be perfect. These requirements are certainly not guaranteed in practice, therefore, prevent the system from being applied widely. In this paper, we extend the framework of the classical teach-and-repeat and propose several new contributions to make our system complete, robust, and flexible. Contributions of this paper are listed as:

1.  We advance our flight corridor generation method. The flight corridor we use now provides much more optimization freedom compared to our previous work [@fei2019ral]. The improvement of the flight corridor facilitates the generation of more efficient and smooth global trajectories. Moreover, we propose methods to accelerate the corridor generation on both CPU and GPU.

2.  We introduce our previous works on online mapping [@han2019fiesta] and re-planning [@boyu2019ral] into our system, to improve the robustness against errors of global maps, drifts of localization, and environmental changes and moving obstacles.

3.  We present a whole set of experiments and comparisons in various scenarios to validate our system.

4.  We release all components in the system as open-source packages, which include local/global planning, perception, and localization, and onboard controller.

![image](sys_architecture){width="1.8\columnwidth"}

![The hardware setting of our autonomous drone system. \[fig:sys\_hardware\] ](drone){width="0.9\columnwidth"}

In what follows, we discuss related literature in Sect. \[sec:related\_work\] and introduce our system in Sect. \[sec:system\_overview\]. Our methods for finding a flight corridor consisting of large convex polyhedrons, and spatial-temporal trajectory optimization are detailed in Sect. \[sec:corridor\_generation\] and Sect. \[sec:trajectory\_optimization\], respectively. The local re-planning is introduced in Sect. \[sec:online\_local\_replanning\]. Experimental and benchmarked results are given in Sect. \[sec:results\]. The paper is concluded in Sect. \[sec:conclusion\].

Related Works {#sec:related_work}
=============

**Robotics teach-and-repeat:** Many robotics teach-and-repeat works, especially for mobile robots, have been published in recent years. Most of them focus on improving the accuracy or robustness in repeating/following the path by operators, which is fundamentally different from our motivation. A lidar-based teach-and-repeat system is proposed in [@sprunk2013lidar], where laser scans are used to localize the ground vehicle against its taught path driven by the user. Furgale et al. [@furgale2014there] [@krusi2015lighting] also develop a lidar-based ground robot, which is specially designed for repeating long-term motions in highly dynamic environments. This system equips a local motion planner which samples local trajectory to avoid dynamic elements during route following. A map maintenance module is used to identify moving objects and estimate their velocities. An iterative learning controller is proposed in [@ostafew2013visual], to reduce the tracking error during the repeating of the robot. This controller can compensate disturbances such as unmodelled terrains and environmental changes by learning a feedforward control policy. Vision-based teaching-and-repeat systems are also proposed in several works, such as the visual localization used by the rover in [@furgale2010visual]. In this work, the authors build a manifold map during the teaching and then use it for localization in the repeating, In [@paton2016bridging], a multi-experience localization algorithm is proposed to address the issue of environmental changes. The ground robot is localized robustly against several past experiences. In [@paton2015s] and [@berczi2016s], to further improve the accuracy and robustness in localization, illumination and terrain appearances are considered in their proposed visual navigation system used for teach-and-repeat, Compared to ground teach-and-repeat works, research on aerial teach-and-repeat is few. In [@fehr2018visual], a vision-based drone is used to inspect infrastructure repetitively. In the teaching phase, the desired trajectory is demonstrated by the operator, and some keyframes in the visual SLAM are recorded as checkpoints. While repeating, local trajectories are generated to connect those checkpoints by using the minimum-snap polynomials [@MelKum1105]. To function properly, in this work, the teaching trajectory itself must be smooth, and the environments must have no changes during the whole repeating. In contrast, our proposed method can convert an arbitrarily poor path to a safe and efficient trajectory with expected flying aggressiveness. Also, our system is flexible. Since it records the teaching path by virtually controlling the drone in simulation, a manually piloted teaching process is not necessary. Finally, our proposed system is robust to environmental changing and can even avoid moving obstacles.

**Quadrotor trajectory planning:** Trajectory optimization is essential in generating safe and executable repeating trajectories from poor teaching ones. The minimum-snap trajectory optimization is proposed by Mellinger [@MelKum1105], where piecewise polynomials are used to represent the quadrotor trajectory and are optimized by quadratic programs (QP). A method for solving a closed-form solution of the minimum snap is proposed in [@RicBryRoy1312]. In this work, a safe geometric path is firstly found to guide the generation of the trajectory. By adding intermediate waypoints to the path iteratively, a safe trajectory is finally generated after solving the minimum-snap problem several times. Our previous works [@fei2018icra] [@fei2016ssrr] [@fei2018jfr] carve a flight corridor consisting of simple convex shapes (sphere, cube) in a complex environment. The flight corridor constructed by a series of axis-aligned cubes or spheres can be extracted very fast on occupancy map or Kd-tree. Then we use the flight corridor and physical limits to constrain a piecewise Bézier curve, to generate a guaranteed safe and kinodynamic feasible trajectory, Other works are proposed to find general convex polyhedrons for constraining the trajectory. In [@liu2017ral], a piecewise linear path is used to guide and initialize the polyhedron generation. In [@deits2015computing], by assuming all obstacles are convex, SDP and QP are iteratively solved to find the maximum polyhedron seeded at a random coordinate in 3-D space. Gradient information in maps is also valuable for local trajectory optimization. In CHOMP [@ratliff2009chomp], the trajectory optimization problem is formulated as a nonlinear optimization over the penalty of safety and smoothness. In [@oleynikova2016continuous], [@lin2018autonomous] and [@helen2019system], gradient-based methods are combined with piecewise polynomials for local planning of quadrotors. In this paper, we also utilize gradient-based optimization for local re-planning.

Time optimization or so-called time parametrization is used to optimize the time profile of a trajectory, given the physical limits of a robot. Methods can be divided as direct methods [@choset2005principles] and indirect methods[@roberts2016generating]. Direct methods generate an optimal spatial-temporal trajectory directly in the configuration space. For indirect methods, a trajectory independent of time is firstly generated, the relationship between time and the trajectory is optimized by an additional optimization process. The method in [@jamieson2016near] finds a mapping function between the time and the trajectory, which is done by recursively adding key points into the function, and squeeze out the infeasibility of the time profile. This method obtains an optimal local solution and is computationally expensive. [@roberts2016generating] also proposes a mapping function, which maps time to a virtual parametrization of the trajectory. The mapping function is then optimized under a complicated nonlinear formulation. However, the global optimality is not guaranteed, and a feasible initial solution is necessary to bootstrap the optimization. Convex optimization [@verscheure2009time] and numerical integration [@pham2014general] are two typical methods of robotics time optimal path parametrization (TOPP) problem. Although numerical integration [@pham2014general] [@pham2018new] has shown superior performance in computing efficiency, convex optimization [@verscheure2009time] has the advantage of adding regularization terms other than total time into its objective function. This specialty suits well for our application where the user defines the expected aggressiveness of the drone, and sometimes the drone may not be expected to fly as fast as possible. As for the efficiency, since we do temporal optimization off-line before the repeating, computing time is not critical.

System Overview {#sec:system_overview}
===============

System Architecture
-------------------

The overall software and hardware architecture of our quadrotor system are shown in Fig. \[fig:sys\_architecture\] and \[fig:sys\_hardware\]. The global mapping, flight corridor generator, and global spatial-temporal planner are done on an off-board computer. Other online processings are running onboard on the drone during the flight. Before teaching, the global map is built by onboard sensors. During teaching, a flight corridor is generated by inflating the teaching trajectory. Then the spatial and temporal trajectories are optimized iteratively within the flight corridor under a coordinate descent scheme [@wright2015coordinate]. The local planner using gradient-based optimization is running onboard to avoid unexpected obstacles observed in the repeating flights. For trajectory tracking, we use a geometric controller [@lee2010]. And the attitude is stabilized by the autopilot.

Globally Consistent Localization and Mapping {#subsec:localization_mapping}
--------------------------------------------

We use VINS [@qin2018vins], a robust visual-inertial odometry (VIO) framework, to localize the drone. Moreover, the loop closure detection and global pose graph optimization are used in our system, to globally correct the pose estimation. The global mapping is done by fusing depth measurements from the stereo cameras with the pose estimation. By using our previous research on deformable map [@wang2019surfel], our global mapping module maintains a series of sub-maps with each attached to a keyframe. In this way, the map is attached to the pose graph and is therefore globally driftless. During the mapping, when a loop closure is detected, keyframes in the global pose graph are corrected, and all sub-maps are deformed accordingly. The global pose graph optimization is also activated during the repeating. When loop closure is detected, the pose of the drone is corrected accordingly to eliminate the drift.

Global Spatial-Temporal Planning {#subsec:coordinate_descent}
--------------------------------

For an extremely poor teaching trajectory, both the geometric shape and time profile of it is far from optimal and therefore useless, or even harmful for conducting optimization. However, the topological information of the teaching trajectory is essential since it reflects the human’s intention. To preserve the topological information, we group the free space around the teaching trajectory to form a flight corridor (Sect. \[sec:corridor\_generation\]). The corridor contains the teaching trajectory within it, shares the same topological structure, and provides large freedom for optimization. It’s hard to concurrently optimize a trajectory spatially and temporally in the flight corridor. However, generating a safe spatial trajectory given a fixed time allocation (Sect. \[subsec:space\_optimization\]) and optimizing the time profile of a fixed spatial trajectory (Sect. \[subsec:time\_optimization\]) are both conquerable. Therefore, we iteratively optimize the trajectory in the space-time joint solution space by designing a coordinate descent [@wright2015coordinate] framework. An objective with weighting energy and time duration is defined for optimization. We firstly generate a spatial trajectory whose energy is minimized, then we use the temporal optimization to obtain the optimal time profile of it. The optimal time profile is used to parametrize a trajectory again for spatial optimization. The spatial-temporal optimizations are done iteratively until the total cost cannot be reduced any more.

Local Collision Avoidance {#subsec:local_collision_avoidance}
-------------------------

In practice, the accumulated drift of VIO is unavoidable, and the recall rate of loop closure is unstable. Although we have built a dense global map, when the drift is significant and not corrected by loop detection in time, the quadrotor may have collisions with obstacles. Moreover, the environment may change or contain moving obstacles. Our previous work [@fei2019ral] has to re-build the map when changes happen and can not deal with dynamic obstacles. To resolve the above issues, we integrate our previous local map fusion module [@han2019fiesta] into our system to detect collisions locally and serves the local trajectory optimization. Also, we propose a sliding-window local replanning method based on our previous research on quadrotor local planning [@boyu2019ral], to avoid collisions on the flight.

In the repeating phase, the drone controls its yaw angle to face its flying direction and build a local map by stereo cameras. We consistently check the local trajectory within a replanning time horizon. If collisions along the local trajectory are reported, replanning is triggered to wrap the trajectory out of obstacles by gradient-based optimization [@boyu2019ral].

![An illustration of free space captured by an axis-aligned cube and a general polyhedron. Obstacles are shown in dashed lines. The blue curve is the teaching trajectory of humans. The red triangle is the seed for finding local free space. The axis-aligned cube and a corresponding general convex polyhedron are shown in yellow and green, respectively. \[fig:cube\_polygon\_compare\]](cube_polygon_compare){width="0.8\columnwidth"}

Flight Corridor Generation {#sec:corridor_generation}
==========================

As stated in Sect. \[subsec:coordinate\_descent\], the first step of our global planning is to build a flight corridor around the teaching trajectory for spatial-temporal trajectory optimization. In our previous work [@fei2019ral], the flight corridor is constructed by finding a series of axis-aligned cubes, which may sacrifice much space, especially in a highly nonconvex environment, as is shown in Fig \[fig:cube\_polygon\_compare\]. A more illustrative comparison is shown in Fig. \[fig:poly\_cube\_rviz\], where the convex polyhedron captures much more free space than the simple cube. Using simple axis-aligned cubes significantly limit the solution space of trajectory optimization, which may result in a poor solution. What’s more, in situations where the free space is very limited, such as flying through a very narrow circle, a cube-based corridor [@fei2019ral] may even fail to cover all teaching trajectory and result in no solutions existing in the corridor. Therefore, to utilize the free space more sufficiently and adapt to even extremely cluttered maps, we propose a method to generate general, free, large convex polyhedrons.

Since the human’s teaching trajectory may be arbitrarily jerky, we cannot assume there is a piecewise linear path to initiate the polyhedron generation, as in [@liu2017ral]. Also, we make no requirements on the convexity of obstacles in the map as in [@deits2015computing]. Our method is based on convex set clustering, which is similar to [@blochliger2018topomap], but is different and advanced at:

1.  We make no assumption on the growing directions of convex clusters and generate completely collision-free polyhedrons based on our dense occupancy map.

2.  We introduce several careful engineering considerations which significantly speed-up the clustering.

3.  We fully utilize the parallel structure of this algorithm and accelerate it over an order of magnitude in GPUs.

4.  We introduce a complete pipeline from building the convex polyhedron clusters to establishing constraints in trajectory optimization.

Convex Cluster Inflation {#subsec:convex_polyhedron_inflation}
------------------------

The core algorithm for the construction of the flight corridor is to find the largest convex free polyhedron at a given coordinate. In this paper, we use an occupancy grid map to represent the environment. Each polyhedron in the flight corridor is the convex hull of a voxel set, which is convex and contains only free voxels. The voxel set is found by clustering as many free voxels as possible around an arbitrary seed voxel. In this paper, we name the voxel set as *convex cluster*, and the process of finding such a set as *convex cluster inflation*. Our method for finding such a *convex cluster* is based on the definition of the convex set:

*Definition*: A set $\mathcal{S}$ in a vector space over $\mathcal{R}$ is called a convex set if the line segment joining any pair of points of $\mathcal{S}$ lies entirely in $\mathcal{S}$. [@lay2007convex].

The pipeline for iteratively inflating such a cluster while preserving convexity is stated in Alg. \[alg:polyhedron\_inflation\]. Our method operates on a 3D occupancy map $\mathcal{M}$ where voxels are labeled as *obstacle* or *free*. Three voxel sets are maintained in the algorithm. $\mathcal{C}$ stands for the targeting convex voxel cluster. $\mathcal{C}^+$ is the set of voxels that are tentative to be added to $\mathcal{C}$ in this iteration. And $\mathcal{C}^*$ contains newly added voxels which preserve the convexity. The cluster inflation starts by adding the seed voxel $p$ to $\mathcal{C}$, and adding all neighboring voxels of $p$ to $\mathcal{C}^+$. In an iteration, each voxel $p^+$ in $\mathcal{C}^+$ is checked whether it can preserve convexity using the function CHECK\_CONVEXITY($p^+, \mathcal{C}$, $\mathcal{M}$). This function, as is shown in Alg. \[alg:convex\_check\], casts rays from $p^+$ to each existing voxels in $\mathcal{C}$. According to the definition of the convex set, $\mathcal{M}$ with $p^+$ is convex if and only if all rays are collision-free. Based on this criteria, qualified voxels are considered as *active* voxels and are added into $\mathcal{C}$ and $\mathcal{C}^*$. Neighboring voxels of all *active* voxels $p^*$ are traversed and collected by the function GET\_NEIGHBORS($\mathcal{C}^*$) for next iteration. The inflation ends when $\mathcal{C}^+$ becomes empty, which implies no additional voxels can be added into $\mathcal{C}$. Fig. \[fig:cluster\_ray\_cast\] also illustrates the procedure of the *convex cluster inflation*.

**Notation**: Seed Voxel: $p^s$, Grid Map: $\mathcal{M}$, Convex Cluster: $\mathcal{C}$,\
Candidate Voxel Set: $\mathcal{C}^+$, Active Voxel Set: $\mathcal{C}^*$ **Input**: $p^s$, $\mathcal{M}$ $\mathcal{C} \:\:\: \leftarrow \{p^s\}$ $\mathcal{C}^* \:   \leftarrow \varnothing$ $\mathcal{C}^+ \leftarrow$ GET\_NEIGHBORS($\mathcal{C}$) $\mathcal{C} \:\:\: \leftarrow \mathcal{C} \:\:\cup p^+ $ $\mathcal{C}^* \: \leftarrow \mathcal{C}^* \cup p^+$ $\mathcal{C}^+ \leftarrow \varnothing$ $\mathcal{C}^+ \leftarrow \mathcal{C}^+ \cup $ GET\_NEIGHBORS($\mathcal{C}^*$) $\mathcal{C}^* \leftarrow \varnothing$ **Output**: $\mathcal{C}$

**Notation**: Ray Cast: $l$ $l \leftarrow$ CAST\_RAY($p^+$, $p$) False True

Having a *convex cluster* consists of a number of voxels, we convert it to the algebraic representation of a polyhedron for following spatial trajectory optimization. Quick hull algorithm [@barber1996quickhull] is adopted here for quickly finding the convex hull of all clustered voxels. The convex hull is in vertex representation (V-representation) $\{V_0, V_1, ..., V_m \}$, and is then converted to its equivalent hyperplane representation (H-representation) by using double description method [@fukuda1995double]. The H-representation of a 3-D polyhedron is a set of affine functions: $$\begin{aligned}
\label{eq:h_representation}
&a_0^x \cdot \textbf{x} + a_0^y \cdot \textbf{y} + a_0^z \cdot \textbf{z} \leq k_0, \nonumber \\
&a_1^x \cdot \textbf{x} + a_1^y \cdot \textbf{y} + a_1^z \cdot \textbf{z} \leq k_1, \nonumber \\
&\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\: \vdots  \\
&a_n^x \cdot \textbf{x} + a_n^y \cdot \textbf{y} + a_n^z \cdot \textbf{z} \leq k_n, \nonumber\end{aligned}$$ where $\{a_n^x, a_n^y, a_n^z \}$ is the normal vector of the 3-D hyperplane and $k_n$ is a constant.

CPU Acceleration {#subsec:cpu_acceleration}
----------------

As is shown in Alg. \[alg:polyhedron\_inflation\], determining whether a voxel preserves the convexity needs to conduct ray-casting to all existing voxels in the *convex cluster*. Iterating with all voxels and rays makes this algorithm impossible to run in real-time, especially when the occupancy grid map has a fine resolution. To make the polyhedron generated in real-time, we take careful engineering considerations on the implementations and propose some critical techniques to significantly increase the overall efficiency.

### Polyhedron Initialization

We initialize each convex cluster as an axis-aligned cube using our previous method [@fei2018icra], which can be done very fast since only index query ($\mathcal{O}$(1)) operations are necessary. After inflating the cube to its maximum volume, as in Fig. \[fig:cube\_polygon\_compare\], we switch to the convex clustering to further group convex free space around the cube.

The proposed polyhedron initialization may result in a final polyhedron different from the one which is clustered from scratch. This is because an axis-aligned cube only inflates in $x, y, z$ directions while a *convex cluster* grows in all possible directions (26-connections in a 3D grid map). However, this initialization process is reasonable. Our purpose is not making each polyhedron optimal but capturing as much as possible free space than a simple cube cannot. In practice, the initialization provides a fast discovery of nearby space which is easy to group, and does not prevent the following *convex cluster inflation* to refine the polyhedron and find sizeable free space. In Sect. \[subsubsec:compare\_corridor\], we show that the initialization process significantly improves the computing efficiency with only a neglectable sacrifice on the volume of the final polyhedron.

### Early Termination

We label all voxels in the cluster as *inner* voxels which inside the *convex cluster*, and *outer* voxels which on the boundary of the *convex cluster*. When traversing a ray from a candidate voxel to a voxel in the *convex cluster*, we early terminate the ray casting when it arrives at a voxel labeled as *inner*.

*Theorem 1*: \[Theo: early\_termination\] The early termination at *inner* voxels is sufficient for checking convexity.

*Proof*: According to the definition of convex set, a ray connecting an *inner* voxel to any other voxel in the *convex cluster* lies entirely in the *convex cluster*. Hence, the extension line of an *inner* voxel must lie inside the *convex cluster* and therefore it must pass the convexity check.

### Voxel Selection

To further reduce the number of voxels that need to be cast rays, given a candidate voxel, only *outer* voxels are used to check its convexity.

*Theorem 2*: \[Theo: voxel\_selection\] Using *outer* voxels of a *convex cluster* is sufficient for checking convexity.

*Proof*: Obviously, the *convex cluster* is a closed set with *outer* voxels at its boundary. The candidate voxel is outside this set. Therefore, casting a ray from any *inner* voxel to the candidate voxel must pass one of the *outer* voxels. According to *Theorem 1*, checking convexity of this ray can terminate after the ray passes an *outer* voxels, which means for a candidate voxel, checking rays cast to *outer* voxels is sufficient.

By introducing above techniques, the proposed *convex cluster inflation* can work in real time for a mediate grid resolution ($0.2m$) on CPUs. The efficacy of these techniques is numerically validated in Sect. \[subsubsec:compare\_corridor\].

GPU Acceleration {#subsec:gpu_acceleration}
----------------

We propose a parallel computing scheme that significantly speeds up the inflation by one order of magnitude where a GPU is available. As is shown in Sec. \[subsec:convex\_polyhedron\_inflation\], when the *convex cluster* discovers a new neighboring voxel, sequentially traversing and checking all rays is naturally parallelizable. With the help of many core GPUs, we can cast rays and check collisions parallelly. Moreover, to fully utilize the massively parallel capability of a GPU, reduce the serialize operations, and minimize the data transferring between CPU and GPU, we examine all potential voxels of the cluster parallelly in one iteration. Instead of discovering a new voxel and checking its rays, we find all neighbors of the active set $\mathcal{C}^*$ and check their rays all in parallel. The detailed procedure is presented in Alg. \[alg:parallel\_polyhedron\_inflation\], where GET\_NEIGHBORS($\mathcal{C}$) collects all neighbors of a set of voxels, and PARA\_CHECK\_CONVEXITY($\mathcal{C}^+, \mathcal{C}$, $\mathcal{M}$) checks the convexity of all candidate voxels parallelly in GPUs. Note that in the serialized version of the proposed method, the voxel discovered earlier may prevent later ones from being clustered, as is illustrated in Fig. \[fig:cluster\_ray\_cast\]. However, in the parallel clustering, all voxels examined at the same time may add conflicting voxels to the cluster. Therefore, we introduce an additional variable, $r$ to record sequential information of voxels. As shown in Alg. \[alg:parallel\_convex\_check\], the kernel function is running on the GPU per block. It checks the ray cast from every candidate voxel in $\mathcal{C}^+$ to a cluster voxel in $\mathcal{C}$ and to each other candidate voxel which has a prior index. After that, the function CHECK\_RESULTS($r$) selects all qualified voxels and adds them into $\mathcal{C}$. Firstly, candidate voxels that have collisions with $\mathcal{C}$ are directly excluded. Then, candidate voxels having collisions with other candidates that have already been added into $\mathcal{C}$ are excluded. In this way, we finally get the same results as in the serialized version of the clustering. The efficacy of the parallel computing is shown in Sect. \[subsubsec:compare\_corridor\].

**Notation**: Parallel Raycasting Result: $r$ $\mathcal{C} \:\:\: \leftarrow \{ p^s \}$ $\mathcal{C}^* \:   \leftarrow \varnothing$ $\mathcal{C}^+      \leftarrow$ GET\_NEIGHBORS($\mathcal{C}$)

/\* *GPU data uploads* \*/ $r \:     \leftarrow$ PARA\_CHECK\_CONVEXITY($\mathcal{C}^+, \mathcal{C}$, $\mathcal{M}$) /\* *GPU data downloads* \*/ $\mathcal{C}^* \: \leftarrow$ CHECK\_RESULTS($r$) $\mathcal{C}   \:\:\: \leftarrow \mathcal{C} \cup \mathcal{C}^*$ $\mathcal{C}^+        \leftarrow$ GET\_NEIGHBORS($\mathcal{C}^*$) **Output**: $\mathcal{C}$

$r[i].status \leftarrow True$ /\* *Kernel function starts* \*/ $l \leftarrow$ CAST\_RAY($p^+_i$, $p$) $r[i].status \leftarrow False$ /\* *Kernel function ends* \*/ /\* *Kernel function starts* \*/ $l \leftarrow$ CAST\_RAY($p^+_i$, $p^+_j$) $r[i].status \leftarrow Pending$ $r[i].list.push\_back(j)$ /\* *Kernel function ends* \*/ $r$

\[marker\] $\mathcal{C}^* \leftarrow \mathcal{C}^* \cup p^+$ $\mathcal{C}^* \leftarrow \mathcal{C}^* \cup p^+$ $\mathcal{C}^*$

Corridor Generation and Loop Elimination {#subsec:loop_elimination}
----------------------------------------

Since the trajectory provided by a user may behave arbitrarily jerky and contain local loops, we introduce a specially designed mechanism to elliminate unnecessary loops, i.e., repeatable polyhedrons. The exclusion of repeatable polyhedrons is essential since in following trajectory optimization (Sect. \[sec:trajectory\_optimization\]), each polyhedron is assigned with one piece of the trajectory. Repeatable polyhedrons would result in an optimized trajectory loops as the user does, which is obviously not efficient. The pipeline of the corridor generation is shown in Alg. \[alg:local\_loop\] and Fig. \[fig:corridor\_pipeline\]. At the beginning of the teaching, the flight corridor is initialized by finding the maximum polyhedron around the position of the drone. Then as the human pilots the drone to move, we keep checking the drone’s position. If it goes outside the last polyhedron ($\mathcal{G}$\[-1\]), we further check whether the drone discovers new free space or not. If the drone is contained within the second last polyhedron ($\mathcal{G}$\[-2\]), we can determine that the teaching trajectory has a loop, as shown in Fig. \[fig:corridor3\]. Then, the last polyhedron in the corridor is regarded as repeatable and is therefore popped out from the corridor. Otherwise, as shown in Fig. \[fig:corridor4\], the drone is piloted to discover new space. Then a new polyhedron $\mathcal{P}$ is inflated and added to the tail of the corridor. The corridor generation is terminated when the teaching finish. The final flight corridor shares the same topological structure with the teaching trajectory since no obstacles are included in the corridor. And it has no unnecessary loops.

**Notation**: Flight Corridor $\mathcal{G}$, Drone Position $p$, Convex Polyhedron $\mathcal{P}$ Initialize : $\mathcal{P} \leftarrow$ CONVEX\_INFLATION($p, \mathcal{M}$) $\mathcal{G}$.push\_back($\mathcal{P}$) $p \leftarrow$ UPDATE\_POSE() $\mathcal{G}$.pop\_back() $\mathcal{P}$ = CONVEX\_INFLATION($p, \mathcal{M}$) $\mathcal{G}$.push\_back($\mathcal{P}$) $\mathcal{G}$

Spatial-Temporal Global Trajecotry Optimization {#sec:trajectory_optimization}
===============================================

Spatial Trajecotry Optimization {#subsec:space_optimization}
-------------------------------

For the spatial optimization, we use the Bernstein basis to represent the trajectory as a piecewise Bézier curve, since it can be easily constrained in the flight corridor by enforcing constraints on control points. An $i^{th}$-order Bernstein basis is: $$b_{n}^{i}(t) = \binom{n}{i} \cdot t^i \cdot (1-t)^{n-i},$$ where $n$ is the degree of the basis, $\binom{n}{i}$ is the binomial coefficient and $t$ is the variable parameterizing the trajectory. An $N$-piece piecewise Bézier curve is written as: $$\label{eq:spatial_curve_n_piece}
 \textit{f}_{\mu}\textit{(t)} =  
 \begin{cases} 
    \sum_{i=0}^n c_{\mu, 1}^ib_{n}^{i}(t / T_1), & t\in[0, T_1], \\
    \sum_{i=0}^n c_{\mu, 2}^ib_{n}^{i}(t / T_2), & t\in[0, T_2], \\
      \:\:\: \:\:\:\:\:\:\:\:\:\vdots      &\:\:\:\:\:\:\:\:\:\vdots \\
    \sum_{i=0}^n c_{\mu, N}^ib_{n}^{i}(t / T_N), &t\in[0, T_N].
   \end{cases}$$ For the $m^{th}$ piece of the curve, $c_{\mu. m}^i$ is the $i^{th}$ control point, and $T_m$ is the time duration. The spatial trajectory is generated in $x, y, z$ dimensions, and $\mu \in {x, y, z}$. $\mu$ is omitted in following derivation for brevity. In this equation, $t$ is scaled by $T_m$ since a standard Bézier curve is defined on $[0,1]$.

Follow the formulation in minimum-snap [@MelKum1105], the squared jerk is minimized in this paper. Since the $3^{rd}$ order derivative of a curve corresponds to the angular velocity, the minimization of jerks alleviates the rotation and therefore facilitates visual tracking. The objective of the piecewise curve is: $$J = \sum_{\mu}^{x,y,z} \sum_{m=1}^N \int_{0}^{T_m} \left(\frac{d^3f_{\mu, m}(t)}{dt^3}\right)^2\, dt.$$ which is in a quadratic form denoted as $\mathbf{c}^T \mathbf{Q} \mathbf{c}$. Here $\mathbf{c}$ is composited by all control points in *x, y, z* dimensions. $\mathbf{Q}$ is a semi-definite Hessian matrix.

For a Bézier curve, its higher order derivatives can be represented by linear combinations of corresponding lower-order control points. For the $1^{st}$ and $2^{nd}$ order derivatives of the $m^{th}$ piece of the curve in Eq. \[eq:spatial\_curve\_n\_piece\], we have: $$\begin{aligned}
\label{eq:bezier_v_m}
& f^\prime_m(t)          = \sum_{i=0}^{n-1} n (c_m^{i+1} - c_m^{i}) b_{n-1}^{i}(\frac{t}{T_m}),  \\
& f^{\prime \prime}_m(t) = \sum_{i=0}^{n-2} n (n - 1) (c_m^{i+2} - 2 c_m^{i+1} + c_m^{i}) b_{n-2}^{i}(\frac{t}{T_m}). \nonumber\end{aligned}$$

### Boundary Constraints

The trajectory has the boundary constraints on the initial state ($p^0, v^0, a^0$) and the final state ($p^f, v^f, a^f$) of the quadrotor. Since a Bézier curve always passes the first and last control points, we enforce the boundary constraints by directly setting equality constraints on corresponding control points in each dimension: $$\begin{aligned}
\label{eq:boundary_constraints}
&c_0^{0} = p^0, \nonumber \\
&c_N^{n}  = p^f, \nonumber \\
&n (c_0^{1} - c_0^{0}) = v^0,  \\
&n (c_N^{n} - c_N^{n-1}) = v^f, \nonumber \\
&n (n - 1) (c_0^{2} - 2 c_0^{1} + c_0^{0}) = a^0, \nonumber \\
&n (n - 1) (c_N^{n} - 2 c_N^{n-1} + c_N^{n-2}) = a^f. \nonumber \end{aligned}$$

### Continuity Constraints

For ensuring smoothness, the minimum-jerk trajectory must be continuous for derivatives up to $2^{nd}$-order at all connecting points on the piecewise trajectory. The continuity constraints are enforced by setting equality constraints between corresponding control points of two consecutive curves. For the $j^{th}$ and $(j+1)^{th}$ pieces of the curve, we can write the equation in each dimension as: $$\begin{aligned}
\label{eq:continunous_constraints}
&  c_j^{n} = c_{j+1}^{0}, \nonumber \\
& (c_j^{n} - c_j^{n-1}) / T_j = (c_{j+1}^{1} - c_{j+1}^{0}) / T_{j+1},  \\
& (c_j^{n} - 2 c_j^{n-1} + c_j^{n-2}) / T^{2}_j = (c_{j+1}^{2} - 2 c_{j+1}^{1} + c_{j+1}^{0}) / T^{2}_{j+1}, \nonumber \end{aligned}$$

### Safety Constraints

The safety of the trajectory is guaranteed by enforcing each piece of the curve to be inside the corresponding polyhedron. Thanks to the convex hull property, an entire Bézier curve is confined within the convex hull formed by all its control points. Therefore we constrain control points using hyperplane functions obtained in Eq. \[eq:h\_representation\]. For the $i^{th}$ control point $c_{j, x}^{i}, c_{j, y}^{i}, c_{j, z}^{i}$ of the $j^{th}$ piece of the trajectory in $x, y, z$ dimensions, constraints are: $$\begin{aligned}
\label{eq:safety_constraints}
&a_0^x \cdot c_{j, x}^{i} + a_0^y \cdot c_{j, y}^{i} + a_0^z \cdot c_{j, z}^{i} \leq k_0, \nonumber \\
&a_1^x \cdot c_{j, x}^{i} + a_1^y \cdot c_{j, y}^{i} + a_1^z \cdot c_{j, z}^{i} \leq k_1, \nonumber \\
&\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\: \vdots  \\
&a_n^x \cdot c_{j, x}^{i} + a_n^y \cdot c_{j, y}^{i} + a_n^z \cdot c_{j, z}^{i} \leq k_n, \nonumber\end{aligned}$$

Constraints in Equs. \[eq:boundary\_constraints\] and \[eq:continunous\_constraints\] are affine equality constraints (*$\mathbf{A}_{eq}\mathbf{c} = \mathbf{b}_{eq}$*) and Eq. \[eq:safety\_constraints\] is in affine in-equality formulation (*$\mathbf{A}_{ie}\mathbf{c} \leq \mathbf{b}_{ie}$*). Finally, the spatial trajectory optimization problem is formulated as a QP as follows: $$\begin{aligned}
\label{eq:spatial_qp_program}
 \text{min}  \:\:\:\:\:\: &\mathbf{c}^T \mathbf{Q} \mathbf{c} \nonumber \\
 \text{s.t.} \:\:\:\:\:\: &\mathbf{A}_{eq} \mathbf{c}   = \mathbf{b}_{eq}, \\
&\mathbf{A}_{ie}\mathbf{c} \leq \mathbf{b}_{ie}.  \nonumber\end{aligned}$$ Unlike our previous works on corridor constrained trajectory [@fei2018icra; @fei2018jfr], here the kinodynamic feasibility (velocity and acceleration) is not guaranteed by adding higher-order constraints into this program, but by temporal optimization (Sect. \[subsec:time\_optimization\]). For a rest-to-rest trajectory, the program in Eq. \[eq:spatial\_qp\_program\] is always mathematically feasible.

Temporal Trajectory Optimization {#subsec:time_optimization}
--------------------------------

![The effect of the temporal optimization. $t$ and $\tau$ are the time profile of the spatial trajectory before and after optimization. \[fig:t\_tau\_timeline\]](t_tau_timeline){width="0.99\columnwidth"}

In spatial optimization, a corridor-constrained spatial trajectory is generated given a fixed time allocation. To optimize the trajectory temporally, we design a re-timing function $\{t(\tau): \tau \rightarrow t \}$ to map the original time variable $t$ to a variable $\tau$. The relation between $\tau$ and $t$ is shown in Fig. \[fig:t\_tau\_timeline\]. In this paper, the re-timing function $t(\tau)$ is named as the temporal trajectory, and finding the optimal $t(\tau)$ is called the temporal optimization. For the N-piece spatial curve defined in Equ. \[eq:spatial\_curve\_n\_piece\], we write $t(\tau)$ as a corresponding N-piece formulation: $$\label{eq:piecewise_tau_function}
t(\tau) =  
 \begin{cases} 
      t_1(\tau), &t_1(0) = 0, t_1(\mathcal{T}^*_1) = T_1, t_1 \in [0, T_1] \\
      t_2(\tau), &t_2(0) = 0, t_2(\mathcal{T}^*_2) = T_2, t_2 \in [0, T_2]   \\
      \:\:\: \vdots     & \:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\vdots \\
      t_N(\tau), &t_N(0) = 0, t_N(\mathcal{T}^*_N) = T_N, t_N \in [0, T_N]
   \end{cases}$$ where $T_1, T_2, ... T_N$ are original time durations of the spatial curve $f_{\mu}(t)$, and $\mathcal{T}^*_1, \mathcal{T}^*_2, ... \mathcal{T}^*_N$ are time durations after temporal optimization. Since physically time only increases, $t(\tau)$ is a monotonically increasing function. Therefore we have $\dot{t}(\tau) \geq 0$. For clarity, in what follows, we use $c^\prime = dc/dt$ to denote taking derivatives with respect to $t$, and $\dot{c} = dc/d\tau$ for taking derivatives with respect to $\tau$. By substituting $t$ with $t(\tau)$ in $f_{\mu}(t)$ and taking derivatives with chain rule, we can write the velocity as: $$\label{eq:velocity}
\dot{f}(t(\tau))  = f^\prime(t) \cdot \dot{t},\nonumber \\$$ \[eq:kinodynamic\] and acceleration as: $$\ddot{f}(t(\tau)) = f^\prime(t) \cdot \ddot{t} + f^{\prime\prime}(t) \cdot \dot{t}^2.$$ The velocity and acceleration are also piecewise functions.

Minimum-Time Formulation {#subsec:minimum_time_form}
------------------------

### Objective {#subsubsec:objective}

The total time $\mathcal{T}$ of the temporal trajectory can be written as: $$\label{eq:objective_raw}
\mathcal{T} =  \int_{0}^{\mathcal{T}} 1 d\tau = \sum_{m = 1}^{N} \int_{0}^{T_m} \frac{1}{\dot{t_m}} dt,$$ considering $\dot{t} = dt/d\tau$. We can introduce a regularization term that penalizes the changing rate of $t$, to trade-off the minimization of time and control extremeness, or so-called motion aggressiveness, in our final temporal trajectory. The objective function is then written as: $$\label{eq:objective}
\mathcal{J} =  \sum_{m = 1}^{N} \int_{1}^{T_m} \Big( \frac{1}{\dot{t_m}} + \rho \cdot \ddot{t_m}^2 \Big) dt,$$ where $\rho$ is a weight of the aggressiveness. By setting a larger $\rho$ we can obtain more gentle motions in the temporal trajectory. If $\rho = 0$, the temporal optimization is solved for generating motions as fast as possible. The motions generated with a large $\rho$ can be viewed in our previous work [@fei2019ral].

Following the direct transcription method in [@verscheure2009time], $\alpha(t)$ and $\beta(t)$ are introduced as two additional piecewise functions: $$\begin{aligned}
\label{eq:a_b_def}
\alpha_m(t) = \ddot{t}_m, \:\:\: \beta_m(t) = \dot{t}^2_m. \:\:\:\:\: m = 1,2,...,N.\end{aligned}$$ According to the relationship between $\ddot{t}_m$ and $\dot{t}$, we can have: $$\begin{aligned}
\beta_m(t) \geq 0, \:\:\: \beta_m^\prime(t) = 2\cdot \alpha_m(t). \label{eq:basic_con_2}\end{aligned}$$ Then the objective function in Equ. \[eq:objective\] is re-formulated as: $$\label{eq:objective_ab}
\mathcal{J} =  \sum_{m = 1}^{N} \int_{0}^{T_m} \Big( \frac{1}{\sqrt{\beta_m(t)}} + \rho \cdot \alpha_m(t)^2 \Big) dt,$$

### Constraints {#subsubsec:temporal_constraints}

The continuities of $t(\tau)$ are enforced by setting constraints between every two consecutive pieces of it. In each dimension $\mu$ $\in {x, y, z}$, we have: $$\begin{aligned}
&f_{\mu, m}^\prime(T_m) \cdot \sqrt{\beta_m(T_m)} = f_{\mu, m+1}^\prime(0) \cdot \sqrt{\beta_{m+1}(0)},\\
&f_{\mu, m}^\prime(T_m) \cdot \alpha_m(T_m) + f_{\mu, m}^{\prime\prime}(T_m) \cdot \beta_m(T_m) \nonumber \\
= &f_{\mu, m+1}^\prime(0) \cdot \alpha_{m+1}(0) + f_{\mu, m+1}^{\prime\prime}(0) \cdot \beta_{m+1}(0)\end{aligned}$$ Then, to satisfy the initial and the final velocity and acceleration $a_{0}, v_0, a_f, v_f$, we set boundary constraints: $$\begin{aligned}
& f_{\mu, 1}^\prime(0) \cdot \sqrt{\beta_1(0)} &= v_0, \\
& f_{\mu, N}^\prime(T_N) \cdot \sqrt{\beta_N(T_N)} &= v_f, \\
& f_{\mu, 1}^\prime(0) \cdot \alpha_1(0) + f_{\mu, 1}^{\prime\prime}(0) \cdot \beta_1(0) &= a_0, \\
& f_{\mu, N}^\prime(T_N) \cdot \alpha_N(T_N) + f_{\mu, N}^{\prime\prime}(T_N) \cdot \beta_N(T_N) &= a_f,\end{aligned}$$ Finally, kinodynamic feasibility constraints are set as: $$\begin{aligned}
& - v_{max} \leq f_{\mu, m}^\prime(t) \cdot \sqrt{\beta_m(t)} \leq v_{max}, \\
& - a_{max} \leq f_{\mu, m}^\prime(t) \cdot \alpha_m(t) + f_{\mu, m}^{\prime\prime}(t) \cdot \beta_m(t) \leq a_{max},\end{aligned}$$ where $v_{max}$ and $a_{max}$ are the physical limits of the drone.

### SOCP Re-formulation {#subsubsec:socp_reform}

The above optimization problem has convex objective and constraints and is, therefore, a convex program. To make it easily solvable, for each piece of the trajectory, $t_m \in [0, T_m]$ is discretized to $t_m^0, t_m^1, ... t_m^{K_m}$ according to a given resolution $\delta t$. $K_m = \lceil T_m / \delta t\rceil + 1$. Then, $\alpha_m(t)$ becomes piecewise constant at each discretization point. According to Equ. \[eq:basic\_con\_2\], $\beta_m(t)$ is piecewise linear. In this way, $\alpha_m(t)$ and $\beta_m(t)$ are modeled by a series of discrete variables $\alpha_m^k$ and $\beta_m^k$, where $\beta_m^k$ is evaluated at $t_m^k$ and $\alpha_m^k$ is evaluated at $(t_m^k + t_m^{k+1}) / 2$.

By applying the above discretization, the objective in Eq. \[eq:objective\_ab\] is derived as: $$\begin{aligned}
\label{eq:objective_dis}
\mathcal{J} = \sum_{m = 1}^{N} \sum_{k = 0}^{K_i-1} \bigg( \frac{2}
  { \sqrt{\beta_{m}^{k+1}} + \sqrt{ \beta_m^{k} } } + \rho \cdot (\alpha_m^k)^2 \bigg) \cdot \delta t,  \end{aligned}$$ which is mathematically equivalent to the affine formulation: $$\label{eq:objective_equ}
\sum_{m = 1}^{N} \sum_{k = 0}^{K_i-1} \Big( 2  \cdot \gamma_m^k + \rho \cdot (\alpha_i^k)^2 \Big) \cdot \delta t,$$ by introducing $\gamma_m^k$ and $$\label{eq:slack_equ_1}
\frac{1}{\sqrt{\beta_m^{k+1}} + \sqrt{ \beta_m^{k}}} \leq \gamma_m^k, \:\: k = 0,...K_i-1; m = 1, ... N$$ as slack variables and additional constraints.

Eq. \[eq:slack\_equ\_1\] is further derived to a quadratic form: $$\begin{aligned}
\frac{1}{\zeta_m^{k+1} + \zeta_m^k} &\leq \gamma_m^k, \:\:\:\:\: k = 0,...K_i-1; m = 1, ... N. \label{eq:slack_equ_2}  \\
\zeta_m^k &\leq \sqrt{\beta_m^k}, \:\: k = 0,...K_i;  \:\:\:\:\:\: m = 1, ... N. \label{eq:slack_equ_3}\end{aligned}$$ by introducing $\zeta_m^k$ as slack variables .

Eq. \[eq:slack\_equ\_2\] can be formulated as a standard rotated quadratic cone: $$\label{eq:objective_r_cone}
2 \cdot \gamma_m^k \cdot \Big( \zeta_m^{k+1} + \zeta_m^k \Big) \geq \sqrt{2}^2,$$ which is denoted as $$\label{eq:objective_r_cone_q}
(\gamma_m^k, \zeta_m^{k+1} + \zeta_m^k, \sqrt{2}) \in Q_r^3.$$ Also, Equ. \[eq:slack\_equ\_3\] can be written as a standard (non-rotated) quadratic cone: $$\label{eq:objective_cones}
\big( \beta_m^k + 1 \big) ^2 \geq \big( \beta_m^k - 1 \big) ^2 + \big( 2 \cdot \zeta_m^k \big)^2,$$ and is denoted as $$\label{eq:objective_cones_q}
(\beta_m^k + 1, \beta_m^k - 1, 2\zeta_m^k) \in Q^3.$$ Finally, a slack variable $s$ is introduced to transform the objective in Equ. \[eq:objective\_equ\] to an affine function: $$\label{eq:objective_equ_2}
 \sum_{m = 1}^{N} \sum_{k = 0}^{K_i-1} (2 \cdot \gamma_m^k + \rho \cdot s) \cdot \delta t,$$ with a rotated quadratic cone: $$\label{eq:objective_cone}
2 \cdot s \cdot 1 \geq \sum_{m=1}^{N}\sum_{k=0}^{K_i-1}(\alpha_m^k)^2,$$ i.e. $$\label{eq:objective_cone_q_r}
(s, 1, \boldsymbol{\alpha}) \in Q_r^{2+\sum_{m=1}^{N}(K_i)},$$ where $\boldsymbol{\alpha}$ contains $\alpha_m^k$ in all pieces of the trajectory.

Also, the discretization is applied to $\alpha_m(t)$ and $\beta_m(t)$ in constraints listed in Sect. \[subsubsec:temporal\_constraints\]. Details are omitted for brevity. After that, we re-formulate these constraints as affine equality and in-equality functions. Besides, although we assume $\alpha_k$ is piecewise constant, we bound the changing rate of $\alpha_k$ considering the response time of the actuators of our quadrotor. We also write this changing rate constraint in an affine form: $$\label{eq:a_d_constraints}
-\delta \alpha \leq (\alpha_m^k - \alpha_m^{k-1}) / \delta t \leq -\delta \alpha,$$ where $\delta \alpha$ (not jerk) is a pre-defined bound of the changing rate of acceleration. Since the difference of $\tau$ between $t_m^k$ and $t_m^{k-1}$ cannot be determined during the optimization, we only bound the changing rate of $\alpha_k$ in $t$ domain.

The temporal optimization problem in Sect. \[subsec:minimum\_time\_form\] is formulated as a standard Second Order Cone Program (SOCP): $$\begin{aligned}
\label{eq:socp_final_form}
\text{min}  \:\:\:\:\:\: & \mathbf{h}^T \boldsymbol{\gamma} + \rho \cdot s, \nonumber \\
\text{s.t.} \:\:\:\:\:\: & \mathbf{A}_{eq} \cdot \mathbf{x} = \mathbf{b}_{eq}, \nonumber \\
& \mathbf{A}_{ie} \cdot \mathbf{x} \leq \mathbf{b}_{ie}, \nonumber \\
&(s, 1, \boldsymbol{\alpha}) \in Q_r^{2+\sum_{m=1}^{N}(K_i)},  m = 1, ... N.  \\
&(\gamma_m^k, \zeta_m^{k+1} + \zeta_m^k, \sqrt{2}) \in Q_r^3, k = 0,..., K_i-1,\nonumber \\
&(\beta_m^k + 1, \beta_m^k - 1, 2\zeta_m^k) \in Q^3, k = 0,..., K_i. \nonumber\end{aligned}$$

Here $\boldsymbol{\gamma}$ and $\mathbf{x}$ consist of all $\gamma^k$ and $\alpha^k, \beta^k, \zeta^k, \gamma^k$. $\delta t$ is the resolution of discretization of the problem. The effect of different $\rho$ and $\delta t$ to the temporal trajectory and a more detailed derivation of the SOCP can be viewed in [@fei2018iros].

In our *teach-repeat-replan* system, since the global repeating trajectory always has static initial and final states, Equ. \[eq:socp\_final\_form\] is always mathematically feasible regardless of the solution of spatial optimization. Because a feasible solution of the optimization program can always be found by infinitely enlarging the time. Combined with the fact that the spatial optimization also always has a solution (Sect. \[subsec:space\_optimization\]), once a flight corridor is given, a spatial-temporal trajectory must exist.

Online Local Re-planning {#sec:online_local_replanning}
========================

In our previous work [@fei2019ral], once the global planning finished, the drone would execute the trajectory without other considerations. This strategy is based on assumptions that 1) the map of the environment is perfectly built and remains intact; 2) globally consistent pose estimation is provided. We use a VIO system with loop closure to correct local pose drifts, and our dense map is globally deformed according to the global pose graph. However, the first assumption does not always hold, especially when new obstacles suddenly appear or the environment changes. As for the second assumption, our global pose estimation relies on the loop closure detection, which also does not guarantee an extremely high recall rate. In situations where there are significant pose drifts but no timely loop closure corrections, the drone may have collisions with obstacles, as in Fig. \[fig:drift\_crash\].

Local Re-planning Framework
---------------------------

To address above issues fundamentally, we propose a local re-planning framework which reactively wraps the global trajectory to avoid unmodeled obstacles. A sliding local map is maintained onboard, where obstacles are fused, and an ESDF (Euclidean Signed Distance Field) [@felzenszwalb2012distance] is updated accordingly. Note that the dense global map is attached to the global pose graph but the local map introduced here is associated with the local VIO frame and sliding with the drone.

![An illustration of colliding with obstacles when there are significant pose drifts but no timely loop closure corrections. Obstacles are depicted in the global frame. The flight path of the drone in the VIO frame is shown in the red curve. But the actual trajectory in the global frame is the blue curve, which collides with obstacles on the global map. \[fig:drift\_crash\] ](drift_crash){width="0.9\columnwidth"}

### ESDF Mapping

![\[fig:esdf\_example\] The local occupancy map its corresponding ESDF map visualized at a given height of $0.6m$. ](esdf_local_rviz){width="0.9\columnwidth"}

We adopt our previous work FIESTA [@han2019fiesta], which is an advanced incremental ESDF [@schouten2010incremental] mapping framework, to build the local map for online re-planning. FIESTA fuses the depth information into a voxel-hashed occupancy map [@klingensmith2015chisel] and updates the distance value of voxels as few as possible using a breadth-first search (BFS) framework. It is lightweight, efficient, and produces near-optimal results. Details can be checked in [@han2019fiesta]. The ESDF is necessary for the following gradient-based trajectory wrapping. An example of a local occupancy map and its corresponding ESDF map are shown in Fig. \[fig:esdf\_example\]. Note, in our system the range of the local map is decided by the range of current depth observation.

### Sliding Window Re-planning

Due to the limited onboard sensing range and computing resources, it is impossible and unnecessary to conduct global re-planning. In this article, we maintain a temporal sliding window over the global trajectory and conduct local re-planning within it. As is shown in Fig. \[fig:replan\_traj\], when obstacles are observed to block the trajectory in the sliding window, a re-planed trajectory is generated to avoid obstacles, and rejoin the global trajectory afterward.

Gradient-Based B-spline Optimization
------------------------------------

### B-spline Trajectory Formulation

A B-spline is a piecewise polynomial function defined by a series of control points $ \{ \mathbf{Q}_{0},\mathbf{Q}_{1}, \cdots, \mathbf{Q}_{N}  \} $ and knot vector $ [ t_{0}, t_{1}, \cdots, t_{m} ] $. For a $p$-degree B-spline, we have $ m = N+p+1 $. Following the matrix representation of the De Boor–Cox formula [@de1971subroutine], the value of a B-spline can be evaluated as: $$\label{equ:matrix}
    P(u) = 
    \left[
        1, u, \dots, u^{p} 
    \right] 
    \cdot
    \mathbf{M}_{p+1}
    \cdot
    \left[
        \mathbf{Q}_{i-p}, \mathbf{Q}_{i-p+1}, \dots, \mathbf{Q}_{i}
    \right]^T$$ here $ \mathbf{M}_{p+1} $ is a constant matrix depends only on $ p $. And $ u = (t - t_{i})/ (t_{i+1} - t_{i})$, for $ t \in [t_{i}, t_{i+1}) $.

### B-spline Initialization

We initialize the local trajectory optimization by re-parameterizing the trajectory in the re-planning horizon as a uniform B-spline. The reason we use uniform B-spline is that it has a simple mathematical formula that is easy to evaluate in the following optimization. For a uniform B-splines, each knot span $ \Delta t_{i} = t_{i+1} - t_{i} $ has identical value $ \Delta t $. The local trajectory is first discretized to a set of points according to a given $\Delta t$. Then these points are fitted to a uniform B-spline by solving a min-least square problem.

Note that, a $ p $ degree uniform B-spline is naturally $ p-1 $ order continuous between consecutive spans. Therefore, there is no need to enforce continuity constraints in the following optimization explicitly. Besides, for a $ p $ degree B-spline trajectory defined by $ N+1 $ control points, the first and last $ p $ control points are fixed due to the continuous requirement of the starting and ending states of the local trajectory.

### Elastic Band Optimization

The basic requirements of the re-planed B-spline are three-folds: smoothness, safety, and dynamical feasibility. We define the smoothness cost $ J_{s} $ using a jerk-penalized elastic band cost function[@quinlan1993elastic; @zhu2015convex]: $$\begin{aligned}
\label{equ:elastic_cost}
& J_{s} = \nonumber \\ 
    &\sum\limits_{i=1}^{N-1} \Vert \underbrace{(\mathbf{Q}_{i+2}-2\mathbf{Q}_{i+1}+\mathbf{Q}_{i})}_{\mathbf{F}_{i+1,i}} \ - \ \underbrace{(\mathbf{Q}_{i+1}-2\mathbf{Q}_{i}+\mathbf{Q}_{i-1})}_{\mathbf{F}_{i-1,i}} \Vert^{2}  \nonumber \\
    &= \sum\limits_{i=1}^{N-1} \Vert \mathbf{Q}_{i+2}-3\mathbf{Q}_{i+1}+3\mathbf{Q}_{i}-\mathbf{Q}_{i-1} \Vert^{2},\end{aligned}$$ which can be viewed as a sum of the squared jerk of control points on the B-spline. Note here we use this formulation which is independent of the time parametrization of the trajectory instead of the traditional time integrated cost function [@MelKum1105]. Because the time duration in each span of the B-spline may be adjusted after the optimization (Sec. \[subsubsec:time\_adjust\]), Eq. \[equ:elastic\_cost\] captures the geometric shape of the B-spline regardless of the time parametrization. Besides, Eq. \[equ:elastic\_cost\] bypasses the costly evaluation of the integration and is, therefore, more numerically robust and computationally efficient in optimization.

The safety and dynamical feasibility requirements of the B-spline are enforced as soft constraints and added to the cost function. Also, the collision cost $J_c$, dynamical feasibility cost $J_v$, and $J_a$ are evaluated at only control points. The collision cost $J_c$ is formulated as the accumulated L2-penalized closest distance to obstacles along the trajectory, which is written as $$\label{equ:colli}
    J_{c} = \sum\limits_{i=p}^{N-p} F_{c}(d(\mathbf{Q}_{i})),$$ where $d(\mathbf{Q}_{i})$ is the distance between $ \mathbf{Q}_{i} $ to its closet obstacle and is recorded in the ESDF. $ F_{c} $ is defined as $$\label{equ:potential}
    F_{c}(d) = \left\{ 
    \begin{array}{ccl}
    (d-d_{0})^{2} & & d \le d_{0} \\
    0 & & d > d_{0}
    \end{array}
     \right.$$ where $ d_{0} $ is the expected path clearance. $J_v$ and $J_a$ are applied to velocities and accelerations, which exceed the physical limits. The formulations of $J_v$ and $J_a$ are similar to Eq. \[equ:colli\] and are omitted here. The overall cost function is: $$\label{equ:cost}
    J_{total} = \lambda_{1} J_{s} + \lambda_{2} J_{c} + \lambda_{3} (J_{v} + J_{a}),$$ where $\lambda_{1}, \lambda_{2}, \lambda_{3}$ are weighting coefficients. $J_{total}$ can be minimized for a local optimal solution by general optimization methods such as Gauss-Newton or Levenberg-Marquardt.

### Iterative Refinement {#subsubsec:time_adjust}

In the above-unconstrained optimization problem, although collisions and dynamical infeasibilities are penalized, there is no hard guarantee on generating a strictly feasible solution. To improve the success rate in practice, we add a post-process to refine the trajectory iteratively. In each iteration, we check collisions and feasibilities of all optimized control points. If collisions are detected, we increase the collision term $J_c$ by increasing $\lambda_2$ and solve the optimization problem (Eq. \[equ:cost\]) again.

Since we wrap the local trajectory to go around obstacles, the trajectory is always lengthened after the optimization. Consequently, using the original time parametrization will unavoidably result in a higher aggressiveness, which means the quadrotor tends to fly faster. Then its velocity and acceleration would easily exceed the predefined limits. Therefore, we adjust the time parameterization of the local trajectory to squeeze out dynamical infeasibilities. We slightly enlarge infeasible knots spans of the B-spline by the following heuristic. $$\label{equ:adj2}
\Delta t_{i}^{'} = \min \{ \alpha, \max \{\frac{v_{m}}{v_{max}}, (\frac{a_{m}}{a_{max}})^{\frac{1}{2}}\}\} \cdot \Delta t_{i},$$ where $ \alpha $ is a constant slightly larger than $ 1 $. $ v_{m}, a_{m} $ are infeasible velocity and acceleration and $ v_{max}, a_{max}$ are maximum allowed acceleration and velocity of the drone. The time duration is iteratively enlarged until obtaining a feasible solution or exceeding the maximum iteration limit. If no feasible solution exists after the time adjustment, $\lambda_3$ is increased, and the trajectory is optimized again.

Results {#sec:results}
=======

Implementation Details {#subsec:implementation_details}
----------------------

The global planning method proposed in this paper is implemented with a QP solver OOQP[^4] and a SOCP solver Mosek[^5]. The local re-planning depends on a nonlinear optimization solver NLopt[^6]. The source code of all modules in our quadrotor system, including local/global localization, mapping, and planning, are released as ros-packages[^7] for the reference of community. Readers of this paper can easily replicate all the presented results. The state estimation, pose graph optimization, local mapping, local re-planning, and the controller is running onboard on a Manifold-2C[^8] mini-computer. Other modules are running on an off-board laptop which has a GTX 1080[^9] graphics card.

Our global map is built to attach to a global pose graph. Both the map and the pose graph are saved for repeating. Before the repeating, the drone is handheld to close the loop of the current VIO frame with the saved global pose graph. The relative transformation of these two frames is used to project the control commands to the VIO frame. During the repeating, pose graph optimization is also activated to calculate the pose drift and compensate for the control command.

Simulated Flight Test {#subsec:simulation}
---------------------

We first test our global and local planning methods in simulations. The simulated environments are randomly deployed with various types of obstacles and circles for drone racing, as shown in Fig. \[fig:simulation\_rviz\]. The simulating tool we use is a light-weight simulator MockaFly[^10], which contains quadrotor dynamics model, controller, and map generator. And the simulator’s is also released as an open-source package with this paper. In the simulation, a drone is controlled by a joystick to demonstrate the teaching trajectory. The simulated drone is equipped with a depth camera whose depth measurements are real-time rendered in GPU by back-projecting the drone’s surrounding obstacles. We randomly add noise on the depth measurements to mimic a real sensor. The re-planning module is activated in the simulation and is triggered by the noise added on the depth. The teaching trajectory and the flight corridor is shown in Fig. \[fig:simu1\]. The global trajectory, local re-planned trajectory, and depth measurement are shown in Fig. \[fig:simu6\]. More details about the simulation are presented in the attached video.

Benchmark Comparisons {#subsec:benchmark}
---------------------

### Corridor Generation {#subsubsec:compare_corridor}

We test the performance of the flight corridor generation methods (Sect. \[sec:corridor\_generation\]), to show the efficacy of the proposed techniques for CPU (Sect. \[subsec:cpu\_acceleration\]) and GPU (Sect. \[subsec:gpu\_acceleration\]) accelerations. For convenience, we denote the basic process for doing *convex cluster inflation* as CPU\_raw; CPU\_raw added cube initialization as CPU+; the one with cube initialization, vertex selection and early termination as CPU++; and the parallel version of the *convex cluster inflation* as GPU. We first compare the time consumed for finding the largest flight corridor with these methods, to validate the improvements of efficiency by using our proposed CPU and GPU acceleration techniques. Then, we compare the ratio of space capturing by methods with and without the polyhedron initialization, and by our previous method [@fei2019ral]. The motivation of the latter comparison is two-fold:

1.  It serves to show superior performance by replacing cubes with polyhedrons.

2.  As discussed in Sect. \[subsec:cpu\_acceleration\], the initialization process would result in different final clustering results compared to the pure *convex cluster inflation*. This comparison also validates that the initialization process only makes neglectable harm to free space capturing.

We generate 10 random maps, with 10 $\sim$ 20 random teaching trajectories given in each map. The average length of teaching trajectories is $20 m$. Results are given in Tabs. \[tab:benchmark\_compare\_corridor\_time\] and \[tab:benchmark\_compare\_corridor\_space\].

[|c|c|c|c|c|]{} & & & &\
&&&&\
Res. = $0.25m$ & 0.031& 0.111& 0.162& 0.359\
Res. = $0.20m$ & 0.055& 0.310& 0.503& 1.309\
Res. = $0.15m$ & 0.169& 1.423& 2.803& 9.583\
Res. = $0.10m$ & 0.942 & 13.940 & 30.747 & 141.659\
Res. = $0.075m$ &3.660 & 71.862 &157.181 &927.131\

As shown in Tab. \[tab:benchmark\_compare\_corridor\_time\], as the resolution of the map being finer, the computing time of the simple *convex cluster inflation* quickly becomes unacceptable huge. In CPU, with the help of *polyhedron initialization*, the computational efficiency is improved several times. Moreover, according to Tab. \[tab:benchmark\_compare\_corridor\_time\], introducing the *voxel selection* and *early termination* can increase the speed more than one order of magnitude in a fine resolution. The efficacy of the GPU acceleration is even more significant. As shown in Tab. \[tab:benchmark\_compare\_corridor\_time\], the GPU version improves the computing speed 30 times at a fine resolution ($0.075m$), and 10 times at a coarse resolution ($0.25m$). For a finer resolution, more candidate voxels are discovered in one iteration of Alg. \[alg:parallel\_polyhedron\_inflation\], thus more computations are conducted parallelly to save time.

  ----------------- ------- -------- -------
                                     
                                     
   Res. = $0.25m$    99.22   100.00   82.28
   Res. = $0.20m$    99.56   100.00   82.92
   Res. = $0.15m$    98.93   100.00   81.82
   Res. = $0.10m$    97.06   100.00   82.78
   Res. = $0.075m$   97.14   100.00   83.03
  ----------------- ------- -------- -------

  : Comparison of Space Captured of Corridor Generation[]{data-label="tab:benchmark_compare_corridor_space"}

For the second comparison, we count the number of free voxels included in the flight corridor found by each method. At each resolution, we take the result of the method without initialization as 100% and compare others against it. Tab. \[tab:benchmark\_compare\_corridor\_space\] indicates two conclusions:

1.  Using polyhedrons instead of axis-aligned cubes can significantly increase the volume of the flight corridor.

2.  Using initialization only slightly sacrifices the volume of the flight corridor. And the sacrifice is neglectable in a medium or coarse resolution ($0.15 \sim 0.25 m$).

The first conclusion holds because a simple cube only discovers free space in $x, y, z$ directions and sacrifices much space in a highly nonconvex environment, as in Fig. \[fig:cube\_polygon\_compare\]. The second conclusion comes from the fact that in a highly nonconvex environment, a regular shaped polyhedron (a cube) does not prevent the following voxel clustering in its nearby space. It shows that the initialization plus the clustering refinement does not harm the volume of the final polyhedron, and is acceptable in practice, especially for a resolution not very fine.

### Global Planning {#subsubsec:compare_planning}

![ \[fig:benchmark\_compare\_rviz\] The comparison of trajectories optimized by different methods. The manual flight trajectory is shown as the purple curve. Blue, red, green, and yellow trajectories are generated by our proposed method, our previous method [@fei2019ral], gradient-based method [@fei2017iros] and waypoints-based method [@RicBryRoy1312]. ](rviz_benchmark_2){width="0.99\columnwidth"}

[|c|c|c|c|]{} & & &\
&&&\
**Proposed Method** &**84.607** &**55.154** &**83.350**\
Previous Method &86.723 &57.736 &89.883\
Gradient-based [@fei2017iros] &89.622 &111.398 &109.575\
Waypoint-based [@RicBryRoy1312] &97.045 &94.895 &204.267\

We compare the proposed global planning method against our previous work [@fei2019ral] and other representative optimization-based trajectory generation methods, such as the waypoint-based method [@RicBryRoy1312] and the gradient-based method [@fei2017iros]. For the latter two benchmarked methods, there is no explicit way to capture the topological structure of the teaching trajectory. Therefore, we convert the teaching trajectory to a piecewise path by recursively finding a collision-free straight line path along with it. Then we use this path to initialize the waypoint-based [@RicBryRoy1312] method and the gradient-based method [@fei2017iros]. Benchmarked methods are also integrated into the coordinate descent framework with temporal optimization. Some parameters dominate the performance of these benchmarked methods, especially for the gradient-based method [@fei2017iros] where the trade-off between collision and smoothness is essential. For a fair comparison, parameters are tuned to achieve the best performances before the test. We randomly generate 10 simulated environments with dense obstacles, as in Sect. \[subsec:simulation\], and conduct 10 teach-and-repeat trials in each map. A sample result of generated trajectories is shown in Fig. \[fig:benchmark\_compare\_rviz\].

As shown in Tab. \[tab:benchmark\_compare\], our proposed method outperforms in all length, time, and energy aspects. The waypoint-based [@RicBryRoy1312] method and the gradient-based [@fei2017iros] method both require a piecewise linear path as initialization. The waypoint-based [@RicBryRoy1312] method can only add intermediate waypoints on the initial path. Therefore, it is mostly dominated by its initialization and tends to output a solution with low quality. The gradient-based [@fei2017iros] method has no such restriction and can adjust the path automatically by utilizing gradient information. However, its optimization formulation is underlying non-convex, since the collision cost is defined on a non-convex ESDF. Therefore, the gradient-based [@fei2017iros] method always finds a locally optimal solution around its initial guess. Compared to these two methods, our method is initialization-free. Both the spatial and temporal optimization of our proposed method enjoys the convexity in its formulation. They are guaranteed to find the global energy-optimal and time-optimal solutions in the flight corridor. Naturally, a smoother trajectory also tends to generate a faster time profile. So finally, under the same coordinate descent framework, our method always outperforms [@fei2017iros] and  [@RicBryRoy1312]. Compared to [@fei2019ral], the advanced corridor generation proposed in this paper (Sect. \[sec:corridor\_generation\]) can always capture more free space than using our previous axis-aligned corridor. Naturally, it provides much more freedom for global planning and results in much better solutions.

Indoor Flight Test {#subsec:indoor_exp}
------------------

### Fast Flight in a Static Environment

Firstly, we conduct experiments in a cluttered drone racing scenario. This experiment validates the robustness of our proposed system, and also pushes the boundary of aggressive flight of quadrotors. Several different types of obstacles, including circles, arches, and tunnels, are deployed randomly to composite a complex environment, as shown in Fig. \[fig:indoor\_map\]. The smallest circle only has a diameter of $0.6 m$, which is very narrow compared to the $0.3m \times 0.3m$ tip-to-tip size of our drone. The maximum velocity and acceleration of the drone are set as $3m/s$ and $3m/s^2$, respectively. And the parameter $\rho$ in Equ. \[eq:objective\] is set as 0, which means the quadrotor is expected to fly as fast as possible as long as it respects the kinodynamic limits. A dense global consistent map is pre-built using the method stated in Sect. \[subsec:localization\_mapping\]. During the teaching phase, the quadrotor is virtually piloted by a human to move amid obstacles. Then the quadrotor autonomously converts this teaching trajectory to a global repeating trajectory and starts to track it. Snapshots of the drone in the flight are shown in Fig. \[fig:indoor\_exp\]. The teaching trajectory and the convex safe flight corridor are visualized in Fig \[fig:indoor\_teach\_rviz\_2\]. And the global repeating trajectory is in Fig. \[fig:indoor\_repeat\_rviz\_2\].

### Local Re-planning Against Unknown Obstacles  {#subsec:indoor_dynamic_exp}

Our system can deal with changing environments and moving obstacles. In this experiment, we test our system also in the drone racing site to validate our local re-planning module. Several obstacles are moved or added to change the drone racing environment significantly, and some others are dynamically added during the repeating flight, as shown in Fig. \[fig:indoor\_dyn\_exp\] In this experiment, the maximum velocity and acceleration for the quadrotor are set as $2m/s$ and $2m/s^2$. The local ESDF map is sliding with the drone using a ring-buffered updating mechanism [@usenko2017real]. The resolution of the local perception is $0.075 m$. The size of the map is decided by points observed spreading in the current frame. The horizon and frequency of the local re-planning are $3.5 s$ and $15 Hz$, respectively. Re-planning is triggered 8 times during the flight in this experiment, and local safe and dynamical feasible splines are generated on time accordingly. Local trajectories, local maps, and the overview of this experiment are shown in Fig. \[fig:indoor\_dyn\_repeat\_rviz\]. We refer readers to the attached video for more details.

![ The repeating trajectory in outdoor experiments, trial 1. Marks are interpreted as the same as in previous figures. \[fig:rviz\_outdoor\_1\]](rviz_outdoor_1_1){width="0.99\columnwidth"}

Outdoor Flight Test {#subsec:outdoor_exp}
-------------------

Finally, we conduct quadrotor flight experiments with a much higher aggressiveness in two different outdoor scenes, as in Fig. \[fig:outdoor\_exp\], to show the robustness of our system in natural environments. Although these experiments are conducted outdoor, GPS or other external positioning devices are not used. The teach-repeat-replan pipeline is as the same as before indoor experiments \[subsec:indoor\_exp\]. The maximum allowed velocity and acceleration limits for these two trials are set as $5m/s$, $6m/s^2$ and $7m/s$, $6m/s^2$, respectively. The drone’s desired and estimated positions and velocities in the second trial are given in Fig. \[fig:outdoor\_pv\_plot\], which shows acceptable tracking errors. Since the flight speed is significantly higher than indoor experiments, we set a smaller re-planning horizon as $2.0 s$. Results such as the global and local trajectory and the global map are visualized in Figs. \[fig:rviz\_outdoor\_1\] and \[fig:rviz\_outdoor\_2\]. More clearly visualizations of outdoor experiments are given in the video.

![ Profiles of the desired and estimated position and velocity. The position and velocity are estimated by our localization module VINS [@qin2018vins]. \[fig:outdoor\_pv\_plot\] ](outdoor_pv_plot){width="0.9\columnwidth"}

Conclusion {#sec:conclusion}
==========

In this paper, we propose a framework, *teach-repeat-replan* for quadrotor aggressive flights in complex environments. The main idea of this work is to find the topological equivalent free space of the user’s teaching trajectory, use spatial-temporal trajectory optimization to obtain an energy-efficient repeating trajectory, and incorporate online perception and re-planning to ensure the safety against environmental changes and moving obstacles. The teaching process is conducted by virtually controlling the drone in simulation. The generated repeating trajectory captures users’ intention and respect an expected flight aggressiveness, which enables autonomous flights much more aggressive than human’s piloting in complex environments. The online re-planning guarantees the safety of the flight and also respects the reference of the repeating trajectory.

To group large free space around the teaching trajectory, we propose a GPU-accelerated convex polyhedron clustering method to find a flight corridor. The optimal global trajectory generation problem is decoupled as spatial and temporal sub-problems. Then these two sub-problems are iteratively optimized under the coordinate descent framework. Moreover, we incorporate the local perception and local trajectory re-planning modules into our framework to deal with environmental changes, dynamic obstacles, and localization drifts.

The proposed system is complete and robust. Users of our system do not have to pilot the drone carefully to give a teaching trajectory. Instead, an arbitrarily jerky/poor trajectory can be converted to an efficient and safe global trajectory. Moreover, when the environment changes or the global localization drifts, the local perception and re-planning modules guarantee the safety of the drone while tracking the global trajectory. Our system is also flexible and easily replicable, as evidenced by various types of experiments presented in this paper, and a third-party application[^11]. We release all components of our system for the reference of the community.

[^1]: <https://github.com/HKUST-Aerial-Robotics/Teach-Repeat-Replan>

[^2]: All authors are with the Department of Electronic and Computer Engineering, Hong Kong University of Science and Technology, Hong Kong, China. [$\{$fgaoaa, lwangax, bzhouai, luxin.han, jpanaa, eeshaojie$\}$@ust.hk]{}.

[^3]: $^*$These authors contributed equally to this work

[^4]: <http://pages.cs.wisc.edu/~swright/ooqp/>

[^5]: <https://www.mosek.com>

[^6]: <https://nlopt.readthedocs.io>

[^7]: <https://github.com/HKUST-Aerial-Robotics/Teach-Repeat-Replan>

[^8]: <https://store.dji.com/product/manifold-2?vid=80932>

[^9]: <https://www.nvidia.com/en-us/geforce/20-series/>

[^10]: <https://github.com/HKUST-Aerial-Robotics/mockasimulator>

[^11]: Flight demonstration at the Electrical and Mechanical Services Department (EMSD), Hong Kong government. Video: <https://youtu.be/Ut8WT0BURrM>
