Minimum-time trajectory generation for quadrotors in constrained environments

TL;DR

This paper introduces a new method for computing minimum-time trajectories for quadrotors in environments with obstacles, using a re-formulation that simplifies constraints and fixed horizon problems, validated through simulations and real-world experiments.

Contribution

The paper proposes a novel trajectory generation strategy based on transverse coordinates and spatial parameterization, enabling efficient constrained minimum-time solutions.

Findings

Effective trajectory planning in complex environments demonstrated through simulations.

Successful real-world execution of planned trajectories with a nano-quadrotor.

Re-formulation simplifies handling of constraints and fixed horizon optimization.

Abstract

In this paper, we present a novel strategy to compute minimum-time trajectories for quadrotors in constrained environments. In particular, we consider the motion in a given flying region with obstacles and take into account the physical limitations of the vehicle. Instead of approaching the optimization problem in its standard time-parameterized formulation, the proposed strategy is based on an appealing re-formulation. Transverse coordinates, expressing the distance from a frame path, are used to parameterise the vehicle position and a spatial parameter is used as independent variable. This re-formulation allows us to (i) obtain a fixed horizon problem and (ii) easily formulate (fairly complex) position constraints. The effectiveness of the proposed strategy is proven by numerical computations on two different illustrative scenarios. Moreover, the optimal trajectory generated in the second scenario is experimentally executed with a real nano-quadrotor in order to show its feasibility.

Tables

Table 1. TABLE I: Nomenclature

$𝒑\in {ℝ}^3$	position vector expressed in ${ℱ}_i$
$\text{v}\in {ℝ}^3$	velocity vector expressed in ${ℱ}_i$
$𝚽\in {ℝ}^3$	vector of angles (yaw-pitch-roll w.r.t. current frame)
$R\left(𝚽\right)\in SO\left(3\right)$	rotation matrix to map vectors in ${ℱ}_b$ into vectors in ${ℱ}_i$
$𝝎\in {ℝ}^3$	angular rate vector expressed in ${ℱ}_b$
$\widehat{𝝎}\in so\left(3\right)$	skew-symmetric matrix associated to $𝝎$
$J\left(𝚽\right)\in {ℝ}^{3\times 3}$	matrix mapping $𝝎$ into $\dot{𝚽}$
$m\in ℝ$	vehicle mass
$I\in {ℝ}^{3\times 3}$	inertia matrix
$g\in ℝ$	gravity constant
${𝒆}_3\in {ℝ}^3$	vector defined as ${𝒆}_3:={\left[\mathrm{0\hspace{0.33em}0\hspace{0.33em}1}\right]}^T$
$F\in ℝ$	thrust
$𝜸\in {ℝ}^3$	torque vector

Equations

$$\dot{\mathbf{p}}$$

$$\dot{\text{v}}$$

$$\dot{\mathbf{\Phi }}$$

$$\dot{\mathbf{\omega }}$$

$$\dot{\mathbf{x}}(t)=f(\mathbf{x}(t),\mathbf{u}(t)),$$

$$\begin{array}{rl}\underset{\mathbf{x}(\cdot ),\mathbf{u}(\cdot ),T}{\min }& T\\ \text{subj. to}& \dot{\mathbf{x}}(t)=f(\mathbf{x}(t),\mathbf{u}(t)),\mathbf{x}(0)={\mathbf{x}}_0\text{(dynamics)}\\ & \mathbf{x}(T)\in X_T\text{(final constraint)}\\ & |p(t)|\le p_{max}\text{(roll rate)}\\ & |q(t)|\le q_{max}\text{(pitch rate)}\\ & |r(t)|\le r_{max}\text{(yaw rate)}\\ & 0<F_{min}\le F(t)\le F_{max}\text{(thrust)}\\ & |\phi (t)|\le {\phi }_{max}(t)\text{(roll angle)}\\ & |\theta (t)|\le {\theta }_{max}(t)\text{(pitch angle)}\\ & |\psi (t)|\le {\psi }_{max}(t)\text{(yaw angle)}\\ & c_{obs}(\mathbf{p}(t))\le 0\text{(position constraints)},\end{array}$$

$$\bar{\mathbf{t}}(s)$$

$$\bar{\mathbf{n}}(s)$$

$$\bar{\mathbf{b}}(s)$$

$${\bar{R}}_{SF}:=[\bar{\mathbf{t}}\bar{\mathbf{n}}\bar{\mathbf{b}}]$$

$$\bar{R}^{\prime}_{SF}(s)=\bar{R}_{SF}(s)\left[\begin{array}[]{ccc}0&-\bar{k}(s)&0\\ \bar{k}(s)&0&-\bar{\tau}(s)\\ 0&\bar{\tau}(s)&0\\ \end{array}\right],$$

$$s_f(t):=\text{arg}\underset{s\in {\mathbb{R}}_0^{+}}{\min }\Vert \mathbf{p}(t)-{\bar{\mathbf{p}}}_f(s){\Vert }^2.$$

$$\mathbf{d}$$

$$\begin{array}{rl}{w}_1& :=\bar{\mathbf{n}}(s_f^t{)}^T(\mathbf{p}-{\bar{\mathbf{p}}}_f(s_f^t)),\\ w_2& :=\bar{\mathbf{b}}(s_f^t{)}^T(\mathbf{p}-{\bar{\mathbf{p}}}_f(s_f^t)),\end{array}$$

$$\mathbf{d}=[0w_1{w}_2{]}^T.$$

$$\mathbf{p}(t)={\bar{\mathbf{p}}}_f(s_f^t)+{\bar{R}}_{SF}(s_f^t)\mathbf{d}(t).$$

$$\text{v}(t)={\bar{\mathbf{p}}}_f^{\prime }(s_f^t){\dot{s}}_f^t+{\bar{R}}_{SF}^{\prime }(s_f^t){\dot{s}}_f^t\mathbf{d}(t)+{\bar{R}}_{SF}(s_f^t)\dot{\mathbf{d}}(t).$$

$$\left[\begin{array}[]{c}0\\ \dot{w}_{1}(t)\\ \dot{w}_{2}(t)\end{array}\right]+\dot{s}_{f}^{t}\left[\begin{array}[]{ccc}1-\bar{k}(s_{f}^{t})w_{1}(t)\\ -\bar{\tau}(s_{f}^{t})w_{2}(t)\\ \bar{\tau}(s_{f}^{t})w_{1}(t)\\ \end{array}\right]-\bar{R}_{SF}^{T}(s_{f}^{t})\text{{v}}(t)=0,$$

$${\dot{s}}_f^t$$

$${\dot{w}}_1(t)$$

$${\dot{w}}_2(t)$$

$$\mathbf{w}(t)$$

$$\mathbf{\omega }(t)$$

$$\begin{array}{rl}\dot{\mathbf{w}}(t)& ={\bar{\mathbf{w}}}^{\prime }(s_f^t){\dot{s}}_f^t,\dot{\text{v}}(t)={\bar{\text{v}}}^{\prime }(s_f^t){\dot{s}}_f^t,\dot{\mathbf{\Phi }}(t)={\bar{\mathbf{\Phi }}}^{\prime }(s_f^t){\dot{s}}_f^t,\end{array}$$

$$\begin{array}{rl}{\bar{w}}_1^{\prime }(s_f^t)& =\bar{\mathbf{n}}(s_f^t{)}^T\text{v}(t)\frac{1}{{\dot{s}}_f^t}+\bar{\tau }(s_f^t)w_2(t),\\ {\bar{w}}_2^{\prime }(s_f^t)& =\bar{\mathbf{b}}(s_f^t{)}^T\text{v}(t)\frac{1}{{\dot{s}}_f^t}-\bar{\tau }(s_f^t)w_1(t),\\ {\bar{\text{v}}}^{\prime }(s_f^t)& =(g{\mathbf{e}}_3-\frac{F(t)}{m}R(\mathbf{\Phi }(t)){\mathbf{e}}_3)\frac{1}{{\dot{s}}_f^t},\\ {\bar{\mathbf{\Phi }}}^{\prime }(s_f^t)& =J(\mathbf{\Phi }(t))\mathbf{\omega }(t)\frac{1}{{\dot{s}}_f^t}.\end{array}$$

$$\begin{array}{rl}{\bar{w}}_1^{\prime }& ={\bar{\mathbf{n}}}^T\bar{\text{v}}\frac{1-\bar{k}{\bar{w}}_1}{{\bar{\mathbf{t}}}^T\bar{\text{v}}}+\bar{\tau }{\bar{w}}_2,\\ {\bar{w}}_2^{\prime }& ={\bar{\mathbf{b}}}^T\bar{\text{v}}\frac{1-\bar{k}{\bar{w}}_1}{{\bar{\mathbf{t}}}^T\bar{\text{v}}}-\bar{\tau }{\bar{w}}_1,\\ {\bar{\text{v}}}^{\prime }& =(g{\mathbf{e}}_3-\frac{\bar{F}}{m}R(\bar{\mathbf{\Phi }}){\mathbf{e}}_3)\frac{1-\bar{k}{\bar{w}}_1}{{\bar{\mathbf{t}}}^T\bar{\text{v}}},\\ {\bar{\mathbf{\Phi }}}^{\prime }& =J(\bar{\mathbf{\Phi }})\bar{\mathbf{\omega }}\frac{1-\bar{k}{\bar{w}}_1}{{\bar{\mathbf{t}}}^T\bar{\text{v}}}.\end{array}$$

$${\bar{\mathbf{x}}}_w^{\prime }(s)=\bar{f}({\bar{\mathbf{x}}}_w(s),\bar{\mathbf{u}}(s)),$$

$${\int }_0^{T}1dt={\int }_{s_f(0)}^{s_f(T)}{\bar{t}}_f^{{}^{\prime }}(s)ds.$$

$${\bar{t}}_f^{{}^{\prime }}(s_f^t)=1/{\dot{s}}_f^t$$

$${\bar{t}}_f^{{}^{\prime }}(s_f^t)$$

$${\bar{t}}_f^{{}^{\prime }}(s)$$

$${\int }_0^L\frac{1-\bar{k}(s){\bar{w}}_1(s)}{\bar{\mathbf{t}}(s{)}^T\bar{\text{v}}(s)}ds.$$

$$c_f({\bar{\mathbf{x}}}_w(L))\le 0,$$

$$\displaystyle c_{f,i}(\bar{x}_{w_{i}}(L))=\Big{(}\frac{2\;\bar{x}_{w_{i}}(L)-(\bar{x}_{w_{i},max}+\bar{x}_{w_{i},min})}{(\bar{x}_{w_{i},max}-\bar{x}_{w_{i},min})}\Big{)}^{2}-1,$$

Full text

I Introduction

Numerous applications involving Unmanned Aerial Vehicles (UAVs), and in particular quadrotors, require them to move inside areas characterized by physical boundaries, obstacles and even tight space constraints (as e.g., urban environments) in order to accomplish their robotics tasks. Such applications are, for example, structural inspections, transportation tasks, surveillance and search and rescue missions. Trajectory generation, a core step for physical task realization [1], becomes extremely challenging in this scenario. A physically realizable trajectory must satisfy (i) the (nonlinear) system dynamics, (ii) the physical limits of the vehicle, such as the maximum thrust, and (iii) the position constraints. Although safety (ensured by a feasible trajectory) is the primary requirement for all applications, trajectory optimization is becoming necessary in different application domanis. The cost to minimize can be, for example, the time to execute a maneuver (in a search and rescue scenario), the energy consumption (during long endurance missions), or the “distance” from a desired unfeasible state-input curve (during inspections). The further requirement of performance optimization poses an additional challenge in the trajectory generation problem.

The problem of computing optimal paths (or trajectories) for UAVs (e.g., [2] and [3]) has received significant attention and a number of algorithms for quadrotors have been proposed to accomplish complex tasks, e.g., landing on a moving target [4] and blind navigation in unknown populated environments [5]. Focusing on collision avoidance, two different approaches, namely reacting or planning, can be applied. The reactive approach is based on navigation laws preventing from possible collisions. It can be performed, e.g., modulating the velocity reference [6], selecting ad-hoc reference way-points [7] and defining an harmonic potential field [8]. On the contrary, the planned approach deals with a problem involving dynamics and state-input constraints with (possibly) a performance criterion to optimize. The majority of the planning algorithms regarding quadrotors, such as [9], [10], [11], [12], [13], takes advantage of the differential flatness property and relies on approximations via motion primitives. When dealing with obstacle dense environments, trajectory generation is often performed using a decoupled approach ([14], [15], [16]). In a first stage, a collision-free path is generated by sampling-based path planning algorithms, such as the Rapidly-exploring Random Tree (RRT) in [14, 15] or the Probabilistic Roadmap (PRM) in [16], and without the dynamics constraint. In a second stage, an optimal trajectory (satisfying the system dynamics) is generated from the collision-free path. Optimization techniques such as [9], [10], [11] can be used at this stage. In order to overcome the limitations due to the decoupled approach, a variant of the RRT algorithm is developed in [17], an approximated dynamics with an a-posteriori correction is used in [18] and a space-parameterized problem reformulation, suitable for modeling complex flight scenarios, is adopted in [12]. Differently from the previous works, in [19] the structure of the minimum-time trajectories is found by the Pontryagin’s minimum principle. Nevertheless, position constraints are not considered. Finally, in [20] a discretized simple point-mass dynamics and approximated convex constraints are considered. The approximation of non-convex constraints into convex ones is also used in [21], in which a sequential convex programming approach is used to achieve a collision free motion for dancing quadrotors.

Our main contribution is the design of an optimization framework to generate feasible minimum-time quadrotor trajectories in structured environments as, e.g., rooms, corridors, passages, or urban areas. Our strategy computes optimal trajectories that satisfy the quadrotor nonlinear dynamics. The strategy can be applied to general models, which may be more complicated than the differentially flat ones. Instead of addressing the minimum-time problem in its standard free-horizon formulation, we derive a fixed-horizon reformulation in which transverse coordinates, expressing the “transverse” distance from a frame path, are used to parameterize the vehicle position. The resulting problem, having a spatial parameter as independent variable, is easier to solve than the time-parametrized one. Position constraints can be easily added into the reformulated problem by defining the constraint boundaries as a function of the spatial parameter and shaping them according to the presence of obstacles. Approximate solutions to the infinite-dimensional optimization problem are numerically computed by combining the Projection Operator Newton method for Trajectory Optimization (PRONTO) [22] with a barrier function approach [23]. This method generates trajectories in a numerically stable manner and guarantees recursive feasibility during the algorithm evolution, i.e., at each algorithm iteration a system trajectory is available. Moreover the approximated solution always satisfies the constraints since the barrier function approach is an interior function method. As an additional contribution, we present numerical computations to show the effectiveness of the proposed strategy on two challenging scenarios. In the first one, the moving space is delimited by rooms with obstacles of different shapes. In the second scenario, the constrained environment is a tubular region delimited by hula hoops. The optimal minimum-time trajectory related to this second scenario is experimentally performed on our nano-quadrotor testbed.

Our algorithm compares to the literature in the following way. The majority of works, such as [9, 10, 11, 12, 13], uses the differential flatness to avoid the integration of nonlinear differential equations, to reduce the order of the problem and to simplify the definition of constraints [9]. On the contrary, our strategy does not rely on the differential flatness hypothesis and thus it can be applied to more complex models. In the previously cited works, the optimization problem is posed in the flat output space, where outputs are approximated using motion primitives, such as polynomial functions [9, 10, 13], B-splines [11], or “convex combinations of feasible paths” [12]. The optimization variables are thus the parameters of the motion primitives. Differently from these works, we do not rely on motion primitives: the state-input trajectory is the optimization variable in our problem formulation. Similarly to the problem formulation in [24], our reformulated minimum-time problem has a spatial parameter, instead of time, as independent variable. While in [24] the maximum velocity profile (for a given path) is computed for a motorcycle model by using a quasi-static approximation of the dynamics, we optimize the whole state-input trajectory and we consider the full nonlinear dynamics of the quadrotor. Finally, other optimization strategies using the PRONTO method are [25] and [26], which aim to compute respectively minimum-energy trajectories for two-wheeled mobile robots and minimum-distance trajectories (from an unfeasible desired maneuver) for UAVs. Differently from these works, we consider a more general three-dimensional space with position constraints and we reformulate the minimum-time problem by using the transverse coordinates.

The paper is organized as follows. In Section II we present the standard formulation of the optimization problem we aim to solve. In Section III our trajectory generation strategy, based on an appealing reformulation of the problem, is described. Finally, in Section IV, we provide numerical computations and experiments, and discuss interesting features of the computed minimum-time trajectories.

II The quadrotor minimum-time problem

We first briefly introduce the quadrotor model used in the paper and then recall the standard problem formulation.

II-A Quadrotor model

The quadrotor dynamics can be described by the so called vectored-thrust dynamical model in [27], where the gravity is the only external force and the generated torque does not influence the translational dynamics, i.e.,

[TABLE]

with $\boldsymbol{p}=[p_{1}\;p_{2}\>p_{3}]^{T}$, $\text{{v}}=[\text{v}_{1}\;\text{v}_{2}\;\text{v}_{3}]^{T}$, $\boldsymbol{\Phi}=[\varphi\;\theta\;\psi]^{T}$, where $\varphi$, $\theta$, $\psi$ are respectively the roll, pitch and yaw angles, and $\boldsymbol{\omega}=[p\;q\;r]^{T}$. The symbols in equations (1-4) are defined in the following table, where $\mathcal{F}_{i}$ and $\mathcal{F}_{b}$ respectively denote the inertial and the body frame.

For the vehicle maneuvering, we adopt a cascade control scheme with an off-board position/attitude control loop and an on-board angular rate controller. Assuming that the virtual control input $\boldsymbol{\omega}$ is tracked by the on-board angular rate controller, we restrict our trajectory generation problem on the position/attitude subsystem (1-3), which can be written in state-space form as

[TABLE]

with state $\boldsymbol{x}=[\boldsymbol{p}^{T}\;\text{{v}}^{T}\;\boldsymbol{\Phi}^{T}]^{T},$ input $\boldsymbol{u}=[\boldsymbol{\omega}^{T}\;F]^{T}$ and suitably defined $f$.

II-B Quadrotor minimum-time problem: standard formulation

We deal with the following optimal control problem:

[TABLE]

where ${X}_{T}\subset{\mathbb{R}}^{9}$ is a desired final region, $p_{max}$, $q_{max}$ and $r_{max}$ are bounds on roll, pitch and yaw rate, respectively, $F_{min}$ and $F_{max}$ are lower and upper bounds on thrust, $\varphi_{max}(\cdot)$, $\theta_{max}(\cdot)$ and $\psi_{max}(\cdot)$ are bounds on roll-pitch-yaw angles, and $c_{obs}:{\mathbb{R}}^{3}\rightarrow{\mathbb{R}}$ represents position constraints. The $t$-dependent constraints in (6) hold for all $t\in[0,T]$, with the exception of $c_{obs}(\boldsymbol{p}(t))\leq 0$, which holds for all $t\in[0,T)$. The bounds on the angular rates avoid fast solutions. The vehicle thrust is also limited: quadrotor vehicles can only generate positive thrust and the maximum rotor speed is limited. Furthermore, constraints on roll and pitch angles are imposed into the optimization problem in order to avoid acrobatic vehicle configurations and to satisfy $\theta\neq\pm\frac{\pi}{2}$ (which makes the matrix $J(\boldsymbol{\Phi})$ in (3) always well defined). Time dependent boundaries can be used for roll and pitch constraints when the vehicle has to move through small passages. The constraint on the yaw angle may be useful in applicative scenarios in which a sensor, e.g., a camera, is provided onboard the vehicle and needs to be pointed toward a target region. The position constraints take into account physical boundaries (possibly shaped by the presence of obstacles) and may also represent GPS denied areas or spaces with limited communication.

III Minimum-Time Trajectory Generation Strategy

In this section, we describe our strategy to compute minimum-time trajectories.

Minimum-time problem (6) is difficult to solve, since it is a constrained, free-horizon problem (time $T$ is an optimization variable). For this reason, instead of directly designing an algorithm to solve problem (6), we provide a strategy to obtain an equivalent, but computationally more appealing, fixed-horizon formulation. In the following, we give an informal idea of the strategy steps to derive the new problem formulation. First, we define a frame path as a (purely geometric) curve in ${\mathbb{R}}^{3}$ used to express the quadrotor position in terms of new coordinates. That is, as depicted in Figure 1, the position is identified by the arc-length of the point on the path at minimum distance and by two transverse coordinates expressing how far the quadrotor position is from the curve. Second, we rewrite the dynamics in terms of the transverse coordinates and show that it depends on time only through the arc-length time-evolution. Thus, by using the arc-length as independent variable, we obtain a “space-dependent” transverse dynamics. Third, the time $T$ can be expressed itself as a function of the arc-length over a fixed “spatial” horizon $[0,L]$, with $L$ being the total length of the frame path. Thus, minimizing $T$ can be rewritten as minimizing an integral function over the fixed spatial interval $[0,L]$. Similarly, pointwise constraints can be written in terms of the transverse coordinates and as function of the arc-length.

The resulting fixed-horizon optimal control problem is solved by using the Projection Operator Newton method (PRONTO), [22], combined with a barrier function approach to handle the constraints, [23].

We provide a detailed and formal explanation of the strategy steps in the following subsections.

III-A Frame path

The first step of the strategy is the generation of an arc-length parameterized frame path $\bar{\boldsymbol{p}}_{f}(s)$, $\forall{s}\in[0,L]$, where $s$ is the arc-length of the path and $L$ is its total length. In the following, we denote the arc-length parameterized functions with a bar, and the derivative with respect to the arc-length with a prime, i.e., $\bar{\boldsymbol{p}}_{f}^{\prime}(s):=d\bar{\boldsymbol{p}}_{f}(s)/ds$. The frame path $\bar{\boldsymbol{p}}_{f}(\cdot)$ has to be locally a non-intersecting $C^{2}$ curve with non-vanishing $\bar{\boldsymbol{p}}^{\prime}_{f}(\cdot)$. Note that the frame path is only a geometric path and it is not required to satisfy the position constraints. A possibility is the computation of the frame path as a $C^{\infty}$ geometric curve, e.g., using arctangent functions as in our numerical computations. More details on the frame path used for our numerical computations will be given in Section IV.

The frame path is used to parameterize the inertial position of the vehicle in the new transverse coordinates, as will be clear later. In order to define the transverse coordinates, we consider the Serret-Frenet frame, whose origin has $\bar{\boldsymbol{p}}_{f}(s)$ as coordinates, and defined $\forall s\in[0,L]$. In particular, the tangent, normal and bi-normal vectors, respectively $\bar{\boldsymbol{t}}(s),\bar{\boldsymbol{n}}(s),\bar{\boldsymbol{b}}(s)$, are defined, with components in the inertial frame, as

[TABLE]

where $\bar{k}(s):=\|\bar{\boldsymbol{p}}^{\prime\prime}_{f}(s)\|_{2}$ is the curvature of $\bar{\boldsymbol{p}}_{f}(\cdot)$ at $s$. Moreover, we define the rotation matrix

[TABLE]

mapping vectors with components in the Serret-Frenet frame into vectors with components in the inertial frame. According to the Serret-Frenet formulas [28], the arc-length derivative of the Serret-Frenet rotation matrix is

[TABLE]

where $\bar{\tau}(s):=\bar{\boldsymbol{n}}(s)\;\bar{\boldsymbol{b}}^{\prime}(s)$ is the torsion of $\bar{\boldsymbol{p}}_{f}(\cdot)$ at $s$.

III-B Transverse dynamics

The second step of the strategy is the derivation of the transverse dynamics by using the transverse coordinates defined with respect to the frame path $\bar{\boldsymbol{p}}_{f}(\cdot)$. In order to rewrite the standard dynamics (1-3) into the transverse dynamics, we proceed as follows.

First, we design a change of coordinates from the inertial position $\boldsymbol{p}\in{\mathbb{R}}^{3}$ to the transverse coordinate vector $\boldsymbol{w}\in{\mathbb{R}}^{2}$, such that $\boldsymbol{w}=[w_{1}\;w_{2}]^{T}$, where $w_{1}$ and $w_{2}$ are the transverse coordinates. Let us consider the quadrotor center of mass with position $\boldsymbol{p}(t)$. As depicted in Figure 1, its orthogonal projection on the frame path identifies a point with position $\bar{\boldsymbol{p}}_{f}(s_{f}(t))$, where the function $s_{f}:{\mathbb{R}}_{0}^{+}\rightarrow{\mathbb{R}}_{0}^{+}$ is

[TABLE]

For simplicity, in the following we use $s_{f}^{t}:=s_{f}(t)$ and $\dot{s}_{f}^{t}:=\dot{s}_{f}(t)$. Note that, the minimizing arc-length is unique provided that $\bar{\boldsymbol{p}}_{f}(\cdot)$ is locally a non-intersecting $C^{2}$ curve with non-vanishing $\bar{\boldsymbol{p}}^{\prime}_{f}(\cdot)$. By mapping $\boldsymbol{p}-\bar{\boldsymbol{p}}_{f}(s_{f}^{t})$ into a vector with components in the Serret-Frenet frame attached to $\bar{\boldsymbol{p}}_{f}(s_{f}^{t})$, we obtain

[TABLE]

Noticing that the component related to the tangent vector is always zero by construction, we define the components $w_{1}$ and $w_{2}$ of the transverse error vector $\boldsymbol{w}$ as, respectively, the second and third components of $\boldsymbol{d}$, i.e.,

[TABLE]

and thus obtaining

[TABLE]

Second, we rewrite equation (1) using $\boldsymbol{w}$ instead of $\boldsymbol{p}$. We note that, the invertible function $s_{f}(\cdot)$ provides a change of variables from the time $t$ to the arc-length $s$. A generic arc-length function $\bar{\alpha}(\cdot)$ can be expressed as the time function $\bar{\alpha}(s_{f}(\cdot))$ and its time derivative is $\frac{d\bar{\alpha}(s_{f}(t))}{dt}=\bar{\alpha}^{\prime}(s_{f}^{t})\;\dot{s}_{f}^{t}$. Let us rewrite equation (1). By using equation (13), the position of the quadrotor center of mass $\boldsymbol{p}(t)$, at time instant $t$, can be written as

[TABLE]

Differentiating (16) with respect to time, since (1) holds, we get

[TABLE]

Multiplying both sides of equation (17) by $\bar{R}^{T}_{SF}$, using (11), (15) and $\bar{\boldsymbol{p}}^{\prime}_{f}(s_{f}^{t})=\bar{R}_{SF}(s_{f}^{t})[1\;0\;0]^{T}$, we get

[TABLE]

i.e., using (10),

[TABLE]

Third and final, we rewrite equations (19), (20), (2), (3), by using the arc-length $s$ as independent variable. Let us denote by $\bar{t}_{f}:{\mathbb{R}}_{0}^{+}\mapsto{\mathbb{R}}_{0}^{+}$ the inverse function of $s_{f}:{\mathbb{R}}_{0}^{+}\mapsto{\mathbb{R}}_{0}^{+}$, satisfying $t=\bar{t}_{f}(s_{f}^{t})$. Due to the invertibility of $s_{f}(\cdot)$, a generic time function $\alpha(\cdot)$ can be expressed as the arc-length function $\alpha(\bar{t}_{f}(\cdot))$ and, defining $\bar{\alpha}:=\alpha\circ\bar{t}_{f}$, we have $\alpha(t)=\bar{\alpha}(s_{f}^{t})$. In particular,

[TABLE]

Deriving with respect to time equations (21), we get

[TABLE]

and equations (19),(20),(2),(3) become

[TABLE]

Using (18), (21) and (22), equations (23) depend on time only through the variable $s_{f}^{t}$. Thus, we can rewrite the dynamics in the arc-length, $s\in[0,L]$, domain. Formally, considering $s$ as the independent variable, we get the transverse dynamics

[TABLE]

Note that the dependence by $s$ is omitted for simplicity. Equations (24) can be written in state-space form as

[TABLE]

with state $\bar{\boldsymbol{x}}_{w}=[\bar{\boldsymbol{w}}^{T}\;\bar{\text{{v}}}^{T}\;\bar{\boldsymbol{\Phi}}^{T}]^{T}$, input $\bar{\boldsymbol{u}}=[\bar{\boldsymbol{\omega}}^{T}\;\bar{F}]^{T}$ and suitable $\bar{f}$.

Remark 1

The general theory regarding the transverse coordinates is introduced in [29] and used to design a maneuver regulation controller for a bi-dimensional case in [30]. Differently from [30], we use the transverse coordinates in a more general three-dimensional case and in order to develop a trajectory optimization strategy rather than a controller.

III-C Arc-length parameterization of cost and constraints

The third step of the strategy consists into the reformulation of cost and constraints in problem (6) by using the new (arc-length dependent) variables $\bar{\boldsymbol{x}}_{w}$ and $\bar{\boldsymbol{u}}$.

The cost functional in (6), i.e., $T=\int_{0}^{T}1\;dt$, is rewritten into an arc-length parameterization by considering the change of variable from $t$ to $s$, i.e.,

[TABLE]

Since $\frac{d\bar{t}_{f}(s_{f}(t))}{dt}=\bar{t}_{f}^{\;{}^{\prime}}(s_{f}^{t})\dot{s}_{f}^{t}$ and $\bar{t}_{f}(s_{f}^{t})=t$, we get

[TABLE]

with $\dot{s}_{f}^{t}$ as in (18). Since $w_{1}(t)=\bar{w}_{1}(s_{f}^{t})$ and $\text{{v}}(t)=\bar{\text{{v}}}(s_{f}^{t})$, as in (21), equation (26) can be written as

[TABLE]

where all the variables depend on time only through $s_{f}^{t}$. Thus, we can rewrite (27) in the arc-length, $s\in[0,L]$, domain, obtaining

[TABLE]

Finally, since $s_{f}(0)=0$, $s_{f}(T)=L$, and (28) holds, we rewrite the cost functional in (6) as

[TABLE]

Notice that, according to (29), the hypothesis $\bar{\boldsymbol{t}}(s)^{T}\bar{\text{{v}}}(s)\neq 0$, has to be satisfied $\forall s\in[0,L]$, i.e., the velocity projected on the tangent vector of the frame path has to be not null.

The constraints in (6) are rewritten into an arc-length parameterization suitable to apply the barrier function approach [23]. The constraint $\boldsymbol{x}(T)\in\boldsymbol{X}_{T}$ is written in the form

[TABLE]

with scalar components

[TABLE]

$\forall i=1,...,8$, where $\bar{x}_{w_{i}}$ is the $i$-th component of $\bar{\boldsymbol{x}}_{w}$, and $\bar{x}_{w_{i},min}$ and $\bar{x}_{w_{i},max}$ are the bounds on the final states. The constraints on the angular rates, thrust and roll-pitch-yaw angles are rewritten by using equations (21), (22) and reparameterizing the time-dependent bounds $\varphi_{max}(t)$, $\theta_{max}(t)$ and $\psi_{max}(t),\;\forall t\in[0,T]$, by the arc-length $s$. Thus, we have

[TABLE]

As regards the position constraints $c_{obs}(\boldsymbol{p}(t))\leq 0$, they are written in the generic form

[TABLE]

which can be particularized according to the shape of the flying region. For environments with circular sections, the inequality (33) becomes

[TABLE]

where $\bar{r}_{obs}(s)$ identifies the radius of the circular boundary at a given arc-length $s$. For environments with rectangular sections, the inequality (33) becomes

[TABLE]

$\forall i=1,2,$ where $\bar{w}_{i,min}(s)$ and $\bar{w}_{i,max}(s)$ are the lower and upper bounds at a given arc-length $s$, defining the boundaries of the region. The constraint boundaries are arc-length functions suitable to model fairly complex regions. They represent the physical boundary of a region and they can be shaped in order to take into account the presence of obstacles attached to the boundary. As an illustrative example, let us consider the environment with rectangular sections depicted in Figure 2. An obstacle restricts the collision-free space inside the physical boundary of the region.

Let us denote by $\bar{w}^{PB}_{i,max}(s)$ and $\bar{w}^{PB}_{i,min}(s)$ the (respectively) positive and negative distance of the physical boundary from the frame path at a given arc-length $s$. We first set $\bar{w}_{i,min}(s)=\bar{w}^{PB}_{i,min}(s)$ and $\bar{w}_{i,max}(s)=\bar{w}^{PB}_{i,max}(s)$, $\forall s\in[0,L]$. Then, in order to take into account the obstacle, we suitably restrict the bounds as follows. Let us consider a point $A$ on the boundary surface of the obstacle. Let $\boldsymbol{p}_{obs}$ be the position of point $A$ with components in the inertial frame. We map $\boldsymbol{p}_{obs}$ in the transverse coordinate vector $\bar{\boldsymbol{w}}_{obs}$. First, we select the arc-length on the frame path, identifying the point at minimum distance from $A$, as

[TABLE]

Second, we map $\boldsymbol{p}_{obs}-\bar{\boldsymbol{p}}_{f}(s_{obs})$ into a vector with components in the Serret-Frenet frame attached to the point identified by $s_{obs}$, obtaining

[TABLE]

Since, according to the particular scenario, the obstacle only affects the function $\bar{w}_{1,min}(\cdot)$, we update

[TABLE]

III-D Equivalent minimum-time formulation and optimal control solver

The minimum-time problem (6) is reformulated in the new (arc-length dependent) variables $\bar{\boldsymbol{x}}_{w}$ and $\bar{\boldsymbol{u}}$, by using the cost (29), the transverse dynamics (25) and the constraints (30), (LABEL:eq:constr_s), and (33). Denoting by $c(\bar{\boldsymbol{x}}_{w}(s),\bar{\boldsymbol{u}}(s))\leq 0,\;\forall s\in[0,L],$ the constraints (LABEL:eq:constr_s) and (33) in vectorial form, the reformulated problem is

[TABLE]

Note that the fixed horizon problem (39) is equivalent to (6) since trajectories solving (39) can be mapped into trajectories solving (6).

In order to solve problem (39), we use a combination of the PRojection Operator based Newton method for Trajectory Optimization (PRONTO) [22] with a barrier function approach [23]. We relax state-input constraints by adding them in the cost functional, i.e., we consider the problem

[TABLE]

where $\epsilon$ and $\epsilon_{f}$ are positive parameters and $\beta_{\ell}(\cdot)$, $\ell\in\{\nu,\nu_{f}\}$, is a function depending on the parameter $\ell$ and defined as

[TABLE]

Let an initial trajectory for the initialization of the algorithm be given. The strategy to find an approximated solution to (39) can be summarized as follows. Problem (40) is iteratively solved by reducing the parameters $\epsilon,\nu,\epsilon_{f}$ and $\nu_{f}$ at each iteration, and thus pushing the trajectory towards the constraint boundaries. Each instance of problem (40) is solved by means of the PRONTO algorithm described in Appendix A.

III-E Summary of the strategy

A pseudo code of the whole strategy to compute minimum time trajectories is reported in the following table (Algorithm 1). We denote with $(\bar{\boldsymbol{x}}_{w}(\cdot),\bar{\boldsymbol{u}}(\cdot))^{0}$ the initial trajectory to initialize the algorithm and with PRONTO the PRojection Operator based Newton method for Trajectory Optimization routine that, given a trajectory $(\bar{\boldsymbol{x}}_{w}(\cdot),\bar{\boldsymbol{u}}(\cdot))^{i-1}$, computes the solution $(\bar{\boldsymbol{x}}_{w}(\cdot),\bar{\boldsymbol{u}}(\cdot))^{i}$ to problem (40), i.e., $(\bar{\boldsymbol{x}}_{w}(\cdot),\bar{\boldsymbol{u}}(\cdot))^{i}=\texttt{PRONTO}((\bar{\boldsymbol{x}}_{w}(\cdot),\bar{\boldsymbol{u}}(\cdot))^{i-1})$.

IV Numerical computations

In this section, we present numerical computations and experimental tests on a nano-quadrotor with mass $m=0.0325\;\text{kg}$, in order to show the effectiveness of the proposed strategy. First, we consider a scenario with two obstacles: a parallelepiped and a cylinder, as depicted in Figure 3a. Second, we consider an experimental scenario and we show the results related to the execution of the optimal trajectory using our maneuver regulation control scheme [31].

IV-A Rooms with obstacles

The first scenario is as follows. The vehicle has to move from one room to another through a narrow corridor. There is a parallelepiped in the first room and a cylinder in the second room. As an additional requirement, the quadrotor must reach a neighborhood of $\boldsymbol{x}_{w0}$ at the end of its motion. In order to fulfil this objective, we consider the final constraint (30) with $c_{f,i}(\bar{x}_{w_{i}}(L))$ as in (31), where $\bar{x}_{w_{i},min}=x_{w_{i},0}-\text{tol}_{i}$, $\bar{x}_{w_{i},max}=x_{w_{i},0}+\text{tol}_{i}$, $\text{tol}_{i}$ is a given tolerance and $x_{w_{i},0}$ is the $i$-th component of $\boldsymbol{x}_{w0}$. Results are depicted in Figures 3, 4, 5. The initial trajectory is depicted in dot-dashed green, intermediate trajectories in dotted black and the minimum-time trajectory in solid blue. Collision-free boundaries are depicted in grey and remaining state-input constraints in dashed red.

We choose as frame path a $C^{\infty}$ curve on the $\bar{p}_{1}-\bar{p}_{2}$ plane with constant binormal vector $\bar{\boldsymbol{b}}=[0\;0\;1]^{T}$ and curvature

[TABLE]

The collision free region is defined by constraint (35) where obstacle boundaries $\bar{w}_{i,min}(\cdot)$ and $\bar{w}_{i,max}(\cdot)$, $i=1,2$, are chosen as follows. Functions $\bar{w}_{1,max}$ and $\bar{w}_{2,max}$ are not affected by obstacles. As depicted in Figures 3b and 3c, $\bar{w}_{1,max}(\cdot)$ and $\bar{w}_{2,max}(\cdot)$ are obtained using sigmoid functions with values varying from $2\;\text{m}$ to $0.25\;\text{m}$. Functions $\bar{w}_{1,min}$ and $\bar{w}_{2,min}$ are affected by obstacles. In order to model the obstacles, we consider the position of the obstacle boundary as a function of its arc-length. We choose $10^{-3}$ as discretization step for the arc-length and for every value of the boundary position $\boldsymbol{p}_{obs}$, we compute $s_{obs}$ and $\bar{w}_{1,obs},\bar{w}_{2,obs}$ by using equations (36), (37) and (38), respectively. Thus, in order to define $\bar{w}_{1,min}$ and $\bar{w}_{2,min}$, we first set $\bar{w}_{1,min}(s)=-\bar{w}_{1,max}(s)$ and $\bar{w}_{2,min}(s)=-\bar{w}_{2,max}(s)$, $\forall s\in[0,L]$. Second, for each $s_{obs}^{R}$ and $\bar{w}_{1,obs}^{R}$ related to a point $R$ on the parallelepiped, we update

[TABLE]

Third and final, for each $s_{obs}^{C}$ and $\bar{w}_{2,obs}^{C}$ related to a point $C$ on the cylinder, we update

[TABLE]

We choose the initial trajectory as follows. We set the frame path as the position, a velocity module of $0.5\;\text{m/s}$ along the curve and a zero yaw angle. The remaining initial states and inputs are computed by using the differential flatness of the quadrotor dynamics [10]. It is worth noting that the position part of the initial trajectory does not have to necessarily match the frame path. Also, the initial trajectory could be alternatively computed through the projection of a state-input curve by using the projection operator (41) described in Appendix A, instead of using the differential flatness.

Having the initial trajectory in hand, we run the algorithm to numerically compute solutions. Note that the PRONTO method (described in Appendix A) is designed considering an $s$-dependent continuous dynamics. In order to implement it by using a numerical toolbox (Matlab), we consider a suitable tolerance. We choose $10^{-3}$ as discretization step on $s$ and we use the tolerance of the Matlab solver to integrate the differential equations. Each intermediate optimal trajectory is computed by solving the optimization problem (40) with constant values of the parameters $\epsilon$, $\nu$, $\epsilon_{f}$ and $\nu_{f}$. We start with $\epsilon=1$, $\nu=1$, $\epsilon_{f}=1$, $\nu_{f}=1$ and, following a suitable heuristic, we decrease them at each iteration. Since the algorithm operates in an interior point fashion, intermediate trajectories are all feasible and are pushed to the constraint boundaries when $\epsilon,\nu,\epsilon_{f},\nu_{f}$ are decreased.

As regards the minimum-time trajectory, the maneuver is performed in $3.57\;\text{s}$ and the path touches the constraint boundaries when the vehicle is inside the corridor and in the proximity of obstacles (Figures 3b and 3c). The velocity $\bar{\boldsymbol{t}}^{T}\bar{\text{{v}}}$ (Figure 4a) reaches a peak of about $8.5\;\text{m/s}$ in the middle of the path and approaches the final desired value at the end. Velocities $\bar{\boldsymbol{n}}^{T}\bar{\text{{v}}}$ and $\bar{\boldsymbol{b}}^{T}\bar{\text{{v}}}$ (Figures 4c and 4e, respectively) are between $-2.0\;\text{m/s}$ and $2.0\;\text{m/s}$. Roll and pitch angles (Figures 4b, 4d, respectively) do not touch constraint boundaries and alternate positive and negative values between $-50\;\text{deg}$ and $50\;\text{deg}$. The yaw angle has values between $-20\;\text{deg}$ and $50\;\text{deg}$ (Figure 4f). As regards the inputs, while constraints on roll and pitch rates (Figures 5a, 5b, respectively) are always active, yaw rate and thrust (Figures 5c and 5d, respectively) alternate intervals with active and inactive constraints. Furthermore, note that the final state reaches a neighborhood of the initial state, satisfying $||\bar{\boldsymbol{x}}_{w}(L)-\boldsymbol{x}_{w0}||<0.07$.

IV-B Tubular passage

As a second test, we consider a region delimited by hula hoops as constrained environment. First, we compute a minimum-time trajectory through our optimization strategy and second, we experimentally execute the minimum-time trajectory on the CrazyFlie nano-quadrotor (https://www.bitcraze.io/crazyflie/), by using a suitable controller. We invite the reader to watch the attached video related to this experiment.

We set-up the optimization algorithm as follows. We approximate the collision free region as a tube with circular section. We choose as frame path a curve on the $\bar{p}_{2}-\bar{p}_{3}$ plane with constant binormal vector $\bar{\boldsymbol{b}}=[1\;0\;0]^{T}$ and curvature

[TABLE]

Moreover, we consider the constraint (34) with constant $\bar{r}_{obs}=r_{hh}-l-e_{p}$, where $r_{hh}=0.33\;\text{m}$ is the hula hoop radius, $l=0.04\;\text{m}$ is the distance between the quadrotor center of mass and its propellers and $e_{p}=0.01\;\text{m}$ is the estimated position error arising during control.

As regards input constraints, we impose, for safety reasons, more severe bounds than the ones required by the physical vehicle limitations. In this way, we also ensure that the “experimental” trajectory remains feasible although the imperfect tracking of desired inputs by actual values (naturally arising during control). We choose $p_{min}=-15\;\text{deg/s}$ and $p_{max}=15\;\text{deg/s}$ for the roll rate and $F_{min}=0.1779\;\text{N}$ and $F_{max}=0.3411\;\text{N}$ for thrust.

By using our minimum-time strategy, we obtain the following result. The optimal trajectory, performed in $2.38$ s, is depicted in Figure 6 (solid blue). Constraint boundaries are depicted in dashed red and the hula hoops are depicted in solid green. The optimal path (blue line in Figure 6a) first takes negative values of $p_{2}$ until changing direction toward positive $p_{2}$ values, touching the constraint boundary in the proximity of the maximum curvature of the tube, and staying in the middle of the feasibility region at the end. The roll angle (blue line in Figure 6b) decreases in order to push the vehicle to negative $p_{2}$ values and then it monotonically increases during the remaining time interval. The velocity module (blue line in Figure 6c) always increases, as we expect for a minimum-time trajectory. As regards the inputs, in the beginning, the angular rate $p$ (blue line in Figure 6d) stays on the lower bound and then it switches to the upper bound. The thrust $F$ (blue line in Figure 6e) always takes the upper bound.

We execute the computed minimum-time trajectory on the CrazyFlie nano-quadrotor by using the closed-loop, maneuver regulation controller developed in [31], in which the minimum-time trajectory is used for the desired maneuver. The maneuver regulation controller computes thrust and angular rate virtual inputs, which are tracked by the standard off-the-shelf angular rate controller provided on board the CrazyFlie. The actual (experimental) trajectory performed using our maneuver regulation controller is depicted in Figure 6 in solid magenta. Snapshots of the experiment are depicted in Figure 6f. As expected, the quadrotor passes close to the second hula hoop maintaining the distance imposed by the restrictive constraints in the optimization problem. The actual velocity does not perfectly match (at higher velocities) the desired one, due to the unmodeled drag effect. Since the vehicle is asked to follow the desired thrust, the actual velocity becomes lower than the desired one because of the opposing aerodynamic force. The experiment shows the actual feasibility of the optimal trajectory and also reveals that a more accurate model including aerodynamic effects would improve the control performance.

V Conclusion

In this paper, we have presented a strategy to address the minimum-time problem for quadrotors in constrained environments. Our approach consists of: (i) generating a frame path, (ii) expressing the quadrotor dynamics in a new set of coordinates “transverse” with respect to that path, and (iii) redefining cost and constraints in the new coordinates. Thus, we obtain a reformulation of the problem, which we solve by combining the PRONTO algorithm with a barrier function approach. Numerical computations on two challenging scenarios prove the effectiveness of the strategy and allow us to show interesting dynamic capabilities of the vehicle. Moreover, the experimental test of the second scenario shows the feasibility of the computed trajectory. As a future work, we aim at extending our strategy to a scenario with moving obstacles. Challenges to be addressed include how to combine trajectory generation and control, and how to take into account a fast integration of the dynamics for realtime computation.

Appendix A Projection Operator Newton Method

Here, we provide a brief description of the PRONTO algorithm [22]. The PRONTO algorithm is based on a properly designed projection operator $\mathcal{P}:\xi_{c}\rightarrow\xi$, mapping a state-control curve $\xi_{c}=(\bar{\boldsymbol{x}}_{w,c}(\cdot),\bar{\boldsymbol{u}}_{c}(\cdot))$ into a system trajectory $\xi=(\bar{\boldsymbol{x}}_{w}(\cdot),\bar{\boldsymbol{u}}(\cdot))$, by the nonlinear feedback system

[TABLE]

where the feedback gain $\bar{K}(\cdot)$ is designed by solving a suitable linear quadratic optimal control problem on the linearized dynamics of (25) about the trajectory $\xi$. Note that the feedback gain $\bar{K}(\cdot)$ is only used to define the projection operator and it is not related to the controller used to execute the optimal trajectory in our experimental test. The projection operator is used to convert the dynamically constrained optimization problem (40) into the unconstrained problem

[TABLE]

where $g(\xi;\bar{k})=h(\mathcal{P}(\xi);\bar{k})$, and $h(\xi;\bar{k}):=\int_{0}^{L}\!(\frac{1-\bar{k}(s)\bar{w}_{1}(s)}{\bar{\boldsymbol{t}}^{T}(s)\bar{\text{{v}}}(s)}+\epsilon\sum_{j}\beta_{\nu}(-c_{j}(\bar{\boldsymbol{x}}_{w}(s),\bar{\boldsymbol{u}}(s))))ds+\;\epsilon_{f}\sum_{i}\beta_{\nu_{f}}(-c_{f,i}(\bar{\boldsymbol{x}}_{w}(L)))$. Then, using an (infinite dimensional) Newton descent method, a local minimizer of (42) is computed iteratively. Given the current trajectory iterate $\xi_{i}$, the search direction $\zeta_{i}$ is obtained by solving a linear quadratic optimal control problem with cost $Dg(\xi_{i};\bar{k})\cdot\zeta+\frac{1}{2}D^{2}g(\xi_{i};\bar{k})(\zeta,\zeta)$, where $\zeta\mapsto Dg(\xi_{i};\bar{k})\cdot\zeta$ and $\zeta\mapsto D^{2}g(\xi_{i};\bar{k})(\zeta,\zeta)$ are respectively the first and second Fréchet differentials of the functional $g(\xi,\bar{k})$ at $\xi_{i}$. Then, the curve $\xi_{i}+\gamma_{i}\zeta_{i}$, where $\gamma_{i}$ is a step size obtained through a standard backtracking line search, is projected, by means of the projection operator, in order to get a new trajectory $\xi_{i+1}$.

The strength of this approach is that the local minimizer of (42) is obtained as the limit of a sequence of trajectories, i.e., curves satisfying the dynamics. Furthermore, the feedback system (41), defining the projection operator, allows us to generate trajectories in a numerically stable manner.

Remark 2

An elegant extension of the PRONTO method to Lie groups is developed in [32] and could be alternatively used in our strategy.