Trajectory Generation for Underactuated Multirotor Vehicles with Tilted Propellers via a Flatness-based Method

TL;DR

This paper introduces a flatness-based trajectory generation method for underactuated multirotor vehicles with tilted propellers, enabling precise 3D path following through a numerical collocation approach.

Contribution

It demonstrates that such vehicles' position and heading are flat outputs and provides a numerical strategy for trajectory planning and control input computation.

Findings

Validated on tricopter, quadrotor, and hexacopter models

Effective trajectory generation for vehicles with tilted propellers

Numerical verification confirms method's practicality

Abstract

This paper considers a class of rotary-wing aerial robots with unaligned propellers. By studying the dynamics of these vehicles, we show that the position and heading angle remain flat outputs of the system (similar to conventional quadrotors). The implication is that they can be commanded to follow desired trajectory setpoints in 3D space. We propose a numerical strategy based on the collocation method to facilitate the trajectory generation. This enables convenient computation of the nominal robot's attitude and control inputs. The proposed methods are numerically verified for three multirotor robots with different dynamics. These include a tricopter with a tilting propeller, and a quadrotor and a hexacopter with unparalleled propellers.

Equations

$$\dot{\mathbf{R}}=\mathbf{R}\hat{\mathbf{\omega }}.$$

$$\bm{\omega}=\left[\begin{array}[]{ccc}1&0&-\sin(\theta)\\ 0&\cos(\phi)&\sin(\phi)\cos(\theta)\\ 0&-\sin(\phi)&\cos(\phi)\cos(\theta)\\ \end{array}\right]\dot{\bm{\theta}}:=\bm{E}(\bm{\theta})\dot{\bm{\theta}}.$$

$$\begin{array}{rl}m\ddot{\mathbf{P}}-m\mathbf{g}& =\mathbf{R}\sum _{i=1}^N({\mathbf{v}}_i{T}_i)\\ & =\mathbf{R}\underset{{\mathbf{A}}_{3\times N}}{\underbrace{\left[{\mathbf{v}}_1,{\mathbf{v}}_2,\cdots ,{\mathbf{v}}_N\right]}}\underset{\mathbf{u}}{\underbrace{{\left[T_1,T_2,\cdots ,T_N\right]}^{\text{T}}}},\end{array}$$

$${\mathbf{M}}_i={\hat{\mathbf{r}}}_i+c_i{\mathcal{I}}_{3\times 3},$$

$$\begin{array}{rl}\sum _{i=1}^N({\mathbf{r}}_i\times ({\mathbf{v}}_i{T}_i)+c_i{\mathbf{v}}_i{T}_i)& =\underset{{\mathbf{B}}_{3\times N}}{\underbrace{[{\mathbf{M}}_1{\mathbf{v}}_1,\cdots ,{\mathbf{M}}_N{\mathbf{v}}_N]}}\mathbf{u}\\ & =\mathbf{J}\dot{\mathbf{\omega }}+\hat{\mathbf{\omega }}\mathbf{J}\mathbf{\omega }.\end{array}$$

$$\underbrace{\left[\begin{array}[]{cc}\bm{A}&\bm{0}\\ -\bm{B}&\bm{J}\\ \end{array}\right]}_{\bm{C}}\left[\begin{array}[]{c}\bm{u}\\ \dot{\bm{\omega}}\\ \end{array}\right]=\left[\begin{array}[]{c}m\bm{R}^{\text{T}}\left(\ddot{\bm{P}}-\bm{g}\right)\\ -\hat{\bm{\omega}}\bm{J}\bm{\omega}\\ \end{array}\right].$$

$$\begin{array}{rl}\mathbf{x}& ={\Phi }_{\mathbf{x}}(\mathbf{\sigma },\dot{\mathbf{\sigma }},\ddot{\mathbf{\sigma }},\dots ,{\mathbf{\sigma }}^{(q)}),\\ \mathbf{u}& ={\Phi }_{\mathbf{u}}(\mathbf{\sigma },\dot{\mathbf{\sigma }},\ddot{\mathbf{\sigma }},\dots ,{\mathbf{\sigma }}^{(q)}),\end{array}$$

$$\mathbf{\sigma }(t)=[x(t),y(t),z(t),\psi (t){]}^{\text{T}}.$$

$$[\mathbf{u},\dot{\mathbf{\omega }}{]}^{\text{T}}=f(\mathbf{\sigma },\mathbf{\theta },\mathbf{\omega }),$$

$$[\mathbf{x},\mathbf{u}{]}^{\text{T}}={\Phi }_{\mathbf{x},\mathbf{u}}(\mathbf{\sigma },\dot{\mathbf{\sigma }},\ddot{\mathbf{\sigma }},\dots ,{\mathbf{\sigma }}^{(q)}),$$

$$m{\bar{\mathbf{R}}}^{\text{T}}({\ddot{\mathbf{\sigma }}}_{1:3}-\mathbf{g})$$

$$\bar{\mathbf{J}}\dot{\bar{\mathbf{\omega }}}+\widehat{\bar{\mathbf{\omega }}}\bar{\mathbf{J}}\bar{\mathbf{\omega }}$$

$${\mathbf{Q}}^{\text{T}}{\mathbf{R}}_x^{\text{T}}(\varphi ){\mathbf{R}}_y^{\text{T}}(\theta ){\mathbf{R}}_z^{\text{T}}({\mathbf{\sigma }}_4)\frac{{\ddot{\mathbf{\sigma }}}_{1:3}-\mathbf{g}}{\Vert ({\ddot{\mathbf{\sigma }}}_{1:3}-\mathbf{g})\Vert }=\frac{[{\lambda }_1{\bar{T}}_1,{\lambda }_2{\bar{T}}_2,0{]}^{\text{T}}}{m\Vert ({\ddot{\mathbf{\sigma }}}_{1:3}-\mathbf{g})\Vert }.$$

$$\begin{array}{rl}{\mathbf{Q}}^{\text{T}}{\mathbf{R}}_x^{\text{T}}(\varphi ){\mathbf{R}}_y^{\text{T}}(\theta )\mathbf{f}/\Vert \mathbf{f}\Vert & =[\cos ({\Theta }_T),-\sin ({\Theta }_T),0{]}^{\text{T}}\\ & ={\mathbf{R}}_z^{\text{T}}({\Theta }_T)[1,0,0{]}^{\text{T}}.\end{array}$$

$${\mathbf{R}}_z^{\text{T}}({\Theta }_T){\mathbf{R}}^{\text{T}}({\Theta }_f)\mathbf{f}/\Vert \mathbf{f}\Vert ={\mathbf{R}}_z^{\text{T}}({\Theta }_T)[1,0,0{]}^{\text{T}}.$$

$${\mathbf{R}}^{\text{T}}({\Theta }_f)\mathbf{f}/\Vert \mathbf{f}\Vert =[1,0,0{]}^{\text{T}}.$$

$$\begin{array}{rl}\bar{\mathbf{\omega }}=& {\mathbf{R}}_z^{\text{T}}({\Theta }_T){\mathbf{R}}^{\text{T}}({\Theta }_f){\mathbf{R}}_z^{\text{T}}({\mathbf{\sigma }}_4)[0,0,{\dot{\mathbf{\sigma }}}_4{]}^{\text{T}}\\ & +{\mathbf{R}}_z^{\text{T}}({\Theta }_T)({\mathbf{R}}^{\text{T}}({\Theta }_f)\dot{\mathbf{R}}({\Theta }_f){)}^{\vee }\\ & +{\mathbf{R}}_z^{\text{T}}({\Theta }_T)[0,0,{\dot{\Theta }}_T{]}^{\text{T}}\\ =& [0,0,{\dot{\Theta }}_T{]}^{\text{T}}+{\mathbf{R}}_z^{\text{T}}({\Theta }_T){\mathbf{R}}^{\text{T}}({\Theta }_f)[0,0,{\dot{\mathbf{\sigma }}}_4{]}^{\text{T}}\\ & +{\mathbf{R}}_z^{\text{T}}({\Theta }_T)({\mathbf{R}}^{\text{T}}({\Theta }_f)\dot{\mathbf{R}}({\Theta }_f){)}^{\vee }\\ =& [0,0,{\dot{\Theta }}_T{]}^{\text{T}}+\mathbf{\xi }({\dot{\mathbf{\sigma }}}_4,{\Theta }_f,{\dot{\Theta }}_f,{\Theta }_T),\end{array}$$

$$\dot{\bar{\mathbf{\omega }}}=[0,0,{\ddot{\Theta }}_T{]}^{\text{T}}+\dot{\mathbf{\xi }}({\dot{\mathbf{\sigma }}}_4,{\ddot{\mathbf{\sigma }}}_4,{\Theta }_f,{\dot{\Theta }}_f,{\ddot{\Theta }}_f,{\Theta }_T,{\dot{\Theta }}_T).$$

$$\bar{\bm{B}}\left[\begin{array}[]{c}\cos\Theta_{T}\frac{\|m\bm{f}\|}{\lambda_{1}}\\ -\sin\Theta_{T}\frac{\|m\bm{f}\|}{\lambda_{2}}\\ \bar{T}_{3}\\ \bar{T}_{4}\\ \end{array}\right]-\hat{\bar{\bm{\omega}}}\bar{\bm{J}}\bar{\bm{\omega}}-\bar{\bm{J}}\dot{\bm{\xi}}=\bar{\bm{J}}\left[\begin{array}[]{c}0\\ 0\\ \ddot{\Theta}_{T}\\ \end{array}\right].$$

$${\ddot{\Theta }}_T=f(\mathbf{\sigma },{\Theta }_T,{\dot{\Theta }}_T).$$

$$\begin{array}{rl}\mathbf{A}{\mathbf{u}}_i& =m\mathbf{R}({\mathbf{\theta }}_i{)}^{\text{T}}({\ddot{\mathbf{\sigma }}}_{1:3,i}-\mathbf{g}),\\ \mathbf{J}({\mathbf{\omega }}_{i+1}-{\mathbf{\omega }}_i)& =t_s(\mathbf{B}{\mathbf{u}}_i-{\hat{\mathbf{\omega }}}_i\mathbf{J}{\mathbf{\omega }}_i),\end{array}$$

$${\mathbf{\theta }}_{i+1}-{\mathbf{\theta }}_i=t_s{\mathbf{E}}^{-1}({\mathbf{\theta }}_i){\mathbf{\omega }}_i,$$

$${\psi }_i={\mathbf{\sigma }}_{4,i}.$$

$${\mathbf{\sigma }}_i(t)=\sum _{j=0}^9{d}_{ij}{t}^j,$$

$${\mathbf{\sigma }}^{(i)}(t)=0\text{for }i=1,\cdots ,4\text{when}t=0,t_f.$$

$$\begin{array}{r}{\sigma }_5=[1,0,0]\frac{\mathbf{A}\mathbf{u}}{\Vert \mathbf{A}\mathbf{u}\Vert },{\sigma }_6=[0,1,0]\frac{\mathbf{A}\mathbf{u}}{\Vert \mathbf{A}\mathbf{u}\Vert },\end{array}$$

$$\underset{{\mathbf{u}}_{\text{n}}}{\min }\frac{1}{N{t}_f}{\int }_{t=0}^{t_f}{\mathbf{u}}_{\text{n}}^{\text{T}}{\mathbf{u}}_{\text{n}}\mathrm{d}\mathrm{t},$$

Full text

I Introduction

In the past decade, we have witnessed an unrivaled popularity of multirotor vehicles. The widespread adoption is attributed by the capability to vertically takeoff and land, as well as its relative ease of use. The mechanical simplicity of conventional quadrotors captures the attention of researchers and engineers and facilitates them to further develop the platforms for numerous applications [1, 2, 3, 4].

More recently, several researchers have pioneered multirotor robots with alternative designs. Among these, the ability to perform thrust-vectoring is becoming prevalent. Additional actuators enable these robots to tilt the propellers, resulting in holonomic aerial platforms [4, 5, 6, 7, 8]. While this ability of thrust-vectoring inevitably comes with added weight and energy consumption, it allows robots to interact more proficiently with the environment [7]. For instance, the robot forcefully interacts with an object via a manipulator [4].

Unlike fully-actuated robots with thrust vectoring capabilities, underactuated aerial vehicles are unable to track arbitrary trajectories in $SE(3)$. To date, there has been little development when it comes to aerial vehicles with tilted propellers that are not fully actuated [9, 10]. In such vehicles, the direction of the total thrust varies in the body frame and the rotational and translational dynamics are sophisticatedly coupled. Without consideration of the distinct dynamic model and suitable control strategies, the progress in flight of multirotor vehicles with unaligned propellers remains infancy.

In this paper, we investigate the flight dynamics of aerial vehicles with unaligned thrusts as depicted in Fig. 1. For regular underacuated robots, it has been shown that the position and heading angle form a set of flat outputs of the system [11, 3]. This means the aerial robots can be controlled to follow 3D trajectory setpoints, rendering them highly versatile platforms for various applications.

Motivated by the benefits of flatness property, we attempt to show that the position and the heading angle remain the flat outputs of the underacuated robots with tilted thrusts. By doing so, corresponding attitude states and control inputs can be found. This trajectory generation strategy, when combined with a compatible flight controller, would potentially lead to highly versatile robots that can perform mid-flight reconfiguration without sacrificing on the controllability.

In the next section, we first outline the dynamic model of robots with unaligned thrusts. We highlight the differences between our model and a model of a regular quadrotor. In section III, we show that the position and heading angle are flat outputs of the system. Then, a computational method for generating the desired trajectory is proposed. The framework is verified using numerical simulations with three multirotor vehicles with different dynamics.

II Dynamic Model

To describe the flight dynamics of a multirotor vehicle with $N$ propellers, we define two coordinate frames: the world frame $\{X_{\text{W}},Y_{\text{W}},Z_{\text{W}}\}$ and the body-fixed frame $\{X_{\text{B}},Y_{\text{B}},Z_{\text{B}}\}$ located on the center of mass (CM) of the robot. A rotation matrix $\bm{R}:=[X_{\text{B}},Y_{\text{B}},Z_{\text{B}}]\in SO(3)$ is defined to relate the attitude of the body frame with respect to the world frame. Let $\bm{\omega}=[\bm{\omega}_{x},\bm{\omega}_{y},\bm{\omega}_{z}]$ represents the angular velocity of the robot and $\hat{\bm{\omega}}$ be its associated skew-symmetric matrix, then the time derivative of $\bm{R}$ satisfies

[TABLE]

For simplicity, the rotation matrix $\bm{R}$ is represented by three Euler angles, yaw ($\psi$), pitch ($\theta$) and roll ($\phi$) in the z-y-x sequence with $\bm{\theta}:=[\phi,\theta,\psi]^{\text{T}}$, then we have $\bm{R}(\bm{\theta})=\bm{R}_{z}(\psi)\bm{R}_{y}(\theta)\bm{R}_{x}(\phi)$. The relationship between these three angles and the angular velocity $\bm{\omega}$ is

[TABLE]

Let $\bm{P}=[x,y,z]^{\text{T}}$ denote the position of the vehicle in the world frame. The robot’s acceleration depends on the vector sum of all propeller’s thrusts. Let $T_{i}$ denote the magnitude of the thrust from the $i^{\text{th}}$ propeller and $\bm{v}_{i}$ be a unit vector describing the direction of $T_{i}$ in the body frame. The translational dynamics of the robot with mass $m$ are

[TABLE]

where we define the matrix $\bm{A}$ to describe the alignment of the propellers and $\bm{g}=[0,0,-g]^{\text{T}}$ is the gravity vector.

To obtain the rotational dynamics, we let $\bm{r}_{i}$ denote the moment arm of the $i^{\text{th}}$ propeller in the body frame. To capture the torque contributed by each propeller, we define the matrix $\bm{M}_{i}$ for the $i^{\text{th}}$ propeller as

[TABLE]

where $\bm{\mathcal{I}}$ is an identity matrix, $\hat{\bm{r}}_{i}$ is the skew-symmetric matrix of $\bm{r}_{i}$, and $c_{i}$ is the aerodynamic constant denoting the ratio of the drag torque to the thrust from the $i^{\text{th}}$ propeller. The sign of $c_{i}$ is determined by the spinning direction of the propeller. Thus, the rotational dynamics of the robot follow

[TABLE]

The translational and rotational dynamics from equations (3) and (5) can be consolidated as

[TABLE]

In the equation above, the orientations of the thrust vectors ($\bm{v}_{i}$’s) are embedded into matrices $\bm{A}$ and $\bm{B}$. One critical condition for the controllability is $\text{rank}(\bm{B})=3$, which means the robot is capable of independently generating torques about its three body axes. For a typical multirotor robot, all propellers are aligned ($\bm{v}_{1}=\bm{v}_{2}=\cdots=\bm{v}_{N}$), and the rank of $\bm{A}$ is one. In such case, equation (3) reduces to $\|m\ddot{\bm{P}}-m\bm{g}\|=\sum^{N}_{i=1}T_{i}$ and the orientation of the robot ($\bm{R}$) is uniquely and immediately determined from $\ddot{\bm{P}}-\bm{g}$ and the desired heading [12]. In this work, we relax on the rank condition of $\bm{A}$ and investigate a class of multirotor robots of which all propellers are not necessarily aligned, or $\text{rank}(\bm{A})\neq 1$.

III Flatness-based Trajectory Generation

In [11], it has been shown that the dynamics of a regular quadrotor with four inputs are differentially flat with the position and heading of the robot constituting four flat outputs. The implication of the result is that, there always exists a feasible trajectory and admissible controls for the chosen flat outputs. With a compatible controller, the robot can be prescribed to follow a setpoint trajectory in 3D space. Without the flatness property, the task of commanding the underactuated robot to traverse from one setpoint to another setpoint requires one to solve the trajectory optimization problem and verify the feasibility. In this section, we employ the differential flatness framework to (i) show that for a class of multirotor robots with unaligned thrusts, there exist a feasible trajectory and control inputs for four flat outputs: position and heading, similar to a standard quadrotor; and (ii) propose a universal trajectory generation method for all related vehicles regardless of the rank condition of $\bm{A}$.

III-A Differential Flatness

For a dynamical system $\dot{\bm{x}}=f\left(\bm{x},\bm{u}\right)$, $\bm{x}\in\mathbb{R}^{m}$, $\bm{u}\in\mathbb{R}^{n}$, it is said to be differentially flat if there exists a flat output: $\bm{\sigma}=h(\bm{x},\dot{\bm{u}},\ddot{\bm{u}},\ldots,\bm{u}^{(p)})$, such that there exist functions $\Phi_{\bm{x}}$ and $\Phi_{\bm{u}}$ [12, 13]:

[TABLE]

where $p$ and $q$ are finite numbers. Equation (7) implies that if the flat output $\bm{\sigma}$ is found, the trajectories of the states and the control inputs are immediately deduced. This is beneficial for motion planning and stabilization of setpoint trajectories. The flatness property is also perceived as an extension of the controllability of linear systems.

In [13], it has been shown that the dynamics of a regular quadrotor is 0-flat, with the corresponding flat output $\bm{\sigma}(t)$,

[TABLE]

As a result, a conventional quadrotor with paralleled thrusts is able to follow a reference trajectory. Without the flatness property, it is necessary to employ trajectory optimization methods to search and verify that there exists a feasible trajectory to command the robot from the starting setpoint to the desired setpoint.

When propellers on a robot are no longer aligned, the translational dynamics are tightly coupled with the rotational dynamics. For instance, equations (3) and (5) imply that it is impossible for the robot to retain a net zero acceleration in the body frame when it generates torque. The difference from a conventional quadrotor means that the flatness proof in [11] is no longer valid. More specifically, the difference in the dynamics equations belonging to the two scenarios is captured into $\bm{A}$ and $\bm{B}$, especially in the rank of matrix $\bm{A}$.

III-B Existence of Flat Outputs

Herein, we revisit the dynamics of multirotor robots with unaligned thrusts and show that the previously identified output in equation (8) is still the flat output to the revised dynamics. Without loss of generality, the analysis in this section assumes $N=4$. In cases where $N>4$, the method can be conveniently adapted as presented in some examples implementation later in Section IV. For illustrative purposes, we formally define a state vector $\bm{x}=[\bm{P},\dot{\bm{P}},\bm{\theta},\bm{\omega}]^{\text{T}}$. Next, we consider two scenarios: $\text{rank}(\bm{A})=3$ and $\text{rank}(\bm{A})=2$.

III-B1 $\text{rank}(\bm{A})=3$

In this circumstance, the rank of the $6\times(N+3)$ matrix $\bm{C}$ in equation (6) is six. With the yaw angle of the robot as part of the flat output $\psi(t)=\bm{\sigma}_{4}(t)$, this results in one further constraint on $\bm{\omega}$, and equivalently on $\dot{\bm{\omega}}$, through equation (2) by the fact that $\ddot{\psi}=\ddot{\bm{\sigma}}_{4}$. Treating the $7\times 1$ vector $[\bm{u},\dot{\bm{\omega}}]^{\text{T}}$ in equation (6) as an unknown, the flat output provides balanced constraints to directly evaluate $\bm{u}$ and $\bm{\omega}$ in the form

[TABLE]

such that there always exist $\bm{u}$ and $\dot{\bm{\omega}}$ given $\bm{\sigma}$, $\bm{\theta}$, and $\bm{\omega}$. It is then straightforward to solve the above equation as a second-order differential equation to simultaneously obtain the states and control inputs given suitable initial conditions. In other words, there exists a numerical function $\Phi_{\bm{x},\bm{u}}$

[TABLE]

similar to equation (7). It can be concluded the system with $\text{rank}(\bm{A})=3$ is differentially flat with the proposed output.

III-B2 $\text{rank}(\bm{A})=2$

In this case, the rank of $\bm{C}$ is 5 and the above method to simultaneously solve for $\bm{u}$ and $\dot{\bm{\omega}}$ no longer applies. Instead, we consider the singular-value decomposition: $\bm{A}=\bm{Q}\bm{\Sigma}\bm{V}^{\text{T}}$, where $\bm{Q}$ and $\bm{V}$ are $3\times 3$ and $N\times N$ orthogonal matrices, and $\bm{\Sigma}$ is a $3\times N$ matrix with two non-zero elements ($\bm{\Sigma}_{1,1}=\lambda_{1}$ and $\bm{\Sigma}_{2,2}=\lambda_{2}$). Exploiting the fact that $\bm{Q}\in SO(3)$, we redefine the body frame with the rotation matrix $\bar{\bm{R}}(\bm{\theta})=\bm{R}(\bm{\theta})\bm{Q}$ and the input vector $\bar{\bm{u}}=\bm{V}^{\text{T}}\bm{u}=\left[\bar{T}_{1},\bar{T}_{2},\cdots,\bar{T}_{N}\right]^{\text{T}}$. The dynamics of the vehicle becomes

[TABLE]

where $\bar{\bm{J}}=\bm{Q}^{\text{T}}\bm{J}\bm{Q}$, $\bar{\bm{B}}=\bm{Q}^{\text{T}}\bm{B}\bm{V}$, and $\bar{\bm{\omega}}=\bm{Q}^{\text{T}}\bm{\omega}$ are re-defined variables in the new body frame. Equation (11) can be manipulated further as

[TABLE]

Leveraging the orthorgonality of rotation matrices, we define $\bm{f}(\bm{\sigma})=\bm{R}_{z}^{\text{T}}(\bm{\sigma}_{4})(\ddot{\bm{\sigma}}_{1:3}-\bm{g})$ and $\Theta_{T}$ that satisfies $\cos(\Theta_{T})=\lambda_{1}\bar{T}_{1}/\|m\bm{f}\|$, $\sin(\Theta_{T})=-\lambda_{2}\bar{T}_{2}/\|m\bm{f}\|$. It follows that

[TABLE]

With $\bar{T}_{i}$’s not yet determined (and likewise $\Theta_{T}$) equation (14) imposes one DoF constraint on the rotation matrix $\bm{R}(\bm{\theta})$. Next, we consider an alternative parametrization of $\bm{R}(\bm{\theta})$. We propose $\bm{R}(\bm{\Theta}):=\bm{R}_{z}(\psi)\bm{R}(\Theta_{f})\bm{R}_{z}(\Theta_{T})=\bar{\bm{R}}(\bm{\theta})$ for $\bm{\Theta}:=[\psi,\Theta_{f},\Theta_{T}]$. As a consequence, equation (14) becomes

[TABLE]

From the equation above, it is straightforward to obtain $\bm{R}(\Theta_{f})$ by solving the equation

[TABLE]

Notice that $\Theta_{f}$ depends entirely on $\bm{f}(\bm{\sigma})$, that is $\Theta_{f}=\Theta_{f}(\bm{\sigma})$. This $\Theta_{f}$ captures the one DoF information on the robot’s orientation embedded in (14). This leaves only one DoF ($\Theta_{T}$) undetermined.

To incorporate this result with the attitude dynamics from equation (12), we analyze the time derivative of $\bm{R}(\bm{\Theta})$. Similar to the process of obtaining equation (2), we find

[TABLE]

where the operation $(\cdot)^{\vee}$ is the inverse of the hat operation: $(\hat{\cdot})^{\vee}=(\cdot)$, and $\bm{\xi}$ is a $3\times 1$ vector from the simplification. The time derivative of $\bar{\bm{\omega}}$ can be computed as

[TABLE]

Finally, with the new parametrization of $\bar{\bm{\omega}}(\bm{\sigma},\bm{\Theta},\dot{\bm{\Theta}})$, the attitude dynamics from equation (12) becomes

[TABLE]

We pre-multiply equation (19) with $\bar{\bm{J}}^{-1}$ and use the top two rows to solve for $\bar{T}_{3}$ and $\bar{T}_{4}$ in terms of $\bm{\sigma},\Theta_{T},\dot{\Theta}_{T}$. The results are substituted back into the last row, which can eventually be expressed as

[TABLE]

With suitable initial conditions, $\bm{\sigma}_{4}$, and previously determined $\Theta_{f}$, solving the respective second-order differential equation provides the complete attitude trajectories. The input $\bm{\bar{u}}$ can then be determined from $\Theta_{T}$. That is, the control input $\bm{u}$ can be found from the flat outputs. This proves that the system is differentially flat with respect to $\bm{\sigma}$.

III-C Trajectory Generation with Direct Collocation

The results in the previous section verify that dynamics described by equation (6) are differentially flat with the proposed flat outputs in equation (8) regardless of the rank of $\bm{A}$. Nevertheless, different approaches are required to numerically compute the corresponding states and inputs.

To facilitate the process, we propose to unify the trajectory generation approach for all cases (including the case for a regular quadrotor or $\text{rank}(\bm{A})=1$) by employing the direct collocation method for the trajectory generation.

The application of the direct collocation, traditionally used for trajectory optimization, hinges on the dynamics in equation (6). While this equation cannot be solved analytically when $\text{rank}(\bm{A})\neq 3$, it can be solved numerically for the entire trajectory by dividing the trajectory into multiple collocation points. For demonstrative purposes, we show an example implementation of the direct transcription using Euler integration. This begins by dividing the trajectory into $n$ time steps: $t=t_{i}$ for $i=1,2,\ldots,n$. In the discrete-time domain of sample time $t_{s}$, equation (6) is approximated as

[TABLE]

which together count as $6n$ scalar equations. The complete dynamics also include additional $3n$ equations from

[TABLE]

where the relationship between $\bm{\theta}_{i}$ and $\bm{\omega}_{i}$, $\bm{E}(\bm{\theta}_{i})$, is determined by equation (2). The yaw angle setpoint implies

[TABLE]

If we treat $\bm{u}_{i}$, $\bm{\theta}_{i}$, and $\bm{\omega}_{i}$ as unknown variables, equations (21)-(23) consist of $10n+6$ unknowns and $10n$ constraints. Additional six initial conditions can be chosen and the nonlinear program is well defined. Owing to the proofs in Section III-B, the existence of the solution is guaranteed.

In practice, a more accurate implementation, e.g., Hermite-Simpson collocation, is preferred. Fortunately, there are a variety of open-source and commercial software packages capable of solving such nonlinear programs. These solvers are conventionally designed for trajectory optimization tasks. This involves an optimization of a cost function for a problem with more variables than constraints, in which oftentimes the feasible solution may not exist. Hence, it is important to emphasize that the implementation of the collocation method for our trajectory generation is different, as our nonlinear problem is not underdetermined and a feasible solution is guaranteed to exist. The benefit of using the collocation method is predominantly the ease of implementation.

IV Numerical Verification

In this section, we show examples of the trajectory generation using the collocation method. Three multirotor robots with different configurations are explored.

IV-A Trajectory Setpoints

The flight trajectory $\bm{\sigma}(t)$ (three positions and the heading angle) for all examples is generated for the time period $t\in[0,t_{f}=4\text{s}]$ using the polynomial path-primitives

[TABLE]

where $d_{ij}$’s are the coefficients corresponding to an arbitrarily chosen trajectory: $\bm{\sigma}(0)=[0.0\text{m},0.0\text{m},0.0\text{m},0.0\text{rad}]^{\text{T}}$ to $\bm{\sigma}(t_{f})=[-1.0\text{m},1.0\text{m},1.5\text{m},0.2\text{rad}]^{\text{T}}$; and the following boundary conditions for the time-derivatives of $\bm{\sigma}(t)$

[TABLE]

These conditions realistically restrict the aerial robots to start and end the trajectory in the hovering condition. The generated trajectory is illustrated in the Fig. 2.

IV-B Example Vehicles with Different Rank Conditions

For demonstration, we select three representative aerial robots as shown in Fig. 3: a quadrotor, a tricopter, and a hexacopter. The quadrotor possesses four tilted propellers with $\text{rank}(\bm{A})=3$. The tricopter has three nominally propellers with one tilting arm. This renders the rank of $\bm{A}$ to be 2. For the hexacopter, $\text{rank}(\bm{A})=3$, but the number of propellers is six. This shows how our proposed method can be modified for cases where $N>4$.

For all robots, the basic configurations ($\bm{r}_{i}$’s and $\bm{v}_{i}$’s) are provided in Fig. 3. The masses are assumed to be $m=1.0$ kg, the inertias are $\bm{J}=\text{diag}(5,5,10)\times 10^{-3}$ kg$\cdot$m2, and $c_{i}=\pm 1.6$ cm for all robots. For the ease of comparison, we define the normalized control input $\bm{u}_{\text{n}}:=N\bm{u}/{mg}$.

The nonlinear program for each robot is solved using a commercial solver (PROPT, Tomlab Optimization [14]) on Matlab (MathWorks). PROPT uses pseudospectral collocation methods and $100$ collocation points are applied.

IV-B1 Quadrotor

The thrust vectors of the quadrotor in Fig. 3(a) are all unparalleled, resulting in a full-rank $\bm{A}$. As stated in Section III, the direction of the total thrust varies with respect to the body frame as the total torque changes.

With the chosen trajectory, the outcomes of the trajectory generation are shown in Fig. 4(a). As anticipated, the quadrotor starts and ends its trajectory in the hovering condition, with non-zero Euler angles. The state variables (except the yaw angle $\bm{\sigma}_{4}$) are identical at the beginning and the end of the trajectory. The $2^{\text{nd}}$ and $3^{\text{rd}}$ propellers are required to generate larger thrusts than the $1^{\text{st}}$ and $4^{\text{th}}$ propellers due to their larger tilt angles with respect to $Z_{B}$. On average, combined thrust is larger than the robot’s weight.

IV-B2 Tricopter

The tricopter has three propellers which are nominally vertically directed. As shown in Fig. 3(b), the third arm is actuated. Its rotation is described by the angle $\alpha$ such that $\bm{v}_{\alpha}=[0,-\sin\alpha,\cos\alpha]^{\text{T}}$ and $\text{rank}(\bm{A})=2$.

With three thrust commands ($T_{1}$, $T_{2}$, and $T_{\alpha}$) and one rotational input $\alpha$, the dynamics of the system with four inputs can be transformed to fit into our framework by defining two hypothetical orthogonal thrust directions: $\bm{v}_{3}=[0,0,1]^{\text{T}}$ and $\bm{v}_{4}=[0,-1,0]^{\text{T}}$ (see Fig. 3(b)). The thrust generated by the tilting arm, $\bm{v}_{\alpha}T_{\alpha}$, becomes $\bm{v}_{3}T_{3}+\bm{v}_{4}T_{4}$ for the newly defined inputs $T_{3}$ and $T_{4}$. The tricopter, consequently, can be mathematically treated as a quadrotor with $\bm{r}_{3}=\bm{r}_{4}$.

The outcomes of trajectory generation are illustrated in Fig. 4(b). The tricopter starts and ends its trajectory in the hovering condition, with non-zero Euler angles. The plot reveals that $T_{\alpha}>T_{1},T_{2}$ and $T_{1}\approx T_{2}$ for the whole flight. It can be deduced that the yaw torques generated by the $1^{\text{st}}$ and $2^{\text{nd}}$ propellers more or less cancel out. Accordingly, the total yaw torque contributed by $T_{\alpha}$ is approximately zero as the aerodynamic drag torque is balanced out by the tilted arm. This is consistent with the fact that the tilt angle $\alpha$ always changes in the opposite direction to its the normalized thrust $\bm{u}_{\text{n},\alpha}$.

IV-B3 Hexacopter

Despite having a full-rank $\bm{A}$, unlike the previous two vehicles, the hexacopter in Fig.3(c) has six propellers ($N=6$). Some adjustment is necessary to adopt our proposed trajectory generation strategy. We present two approaches in the followings.

The first approach is to define two additional (flat) outputs so that the number of the constraints are equal to the number of control inputs. We opt to define ${\sigma}_{5}$ and ${\sigma}_{6}$ such that the resultant dynamics of the hexacopter resemble that of a conventional quadrotor. This is achieved by ensuring the direction of the total thrust $\bm{A}\bm{u}$ remains unchanged during the whole flight. Mathematically, this translates to

[TABLE]

where we choose $\sigma_{5}=0.07$ and $\sigma_{6}=0.06$. In this particular case, the outputs are functions of both $\bm{x}$ and $\bm{u}$. The caveat of this method is that, with the extra outputs, the proof for the existence of the solution in Section III-B no longer applies. If these additional outputs are not true flat outputs of the system, a feasible solution does not exist. It requires some knowledge from the designer to define reasonable outputs.

The second method is to exploit the two unconstrained DoFs of the control inputs and turn the trajectory generation problem into a trajectory optimization problem. This way, a feasible trajectory is still guaranteed to exist. Moreover, the implementation is relatively simple as conventional software packages for solving collocation problems are designed for trajectory optimization. In this case, we set up the optimization problem using the following quadratic input cost

[TABLE]

subject to the system dynamics and the desired flat outputs. For a conventional quadrotor in a perfect hovering condition, the expected cost is unity.

The generated states and control inputs for the hexacopter from both strategies are shown in Fig. 5. Due to the asymmetric configuration, some propellers are required to generate larger thrusts than others. At the first glance, we observe minor differences between the results. A closer inspection into Fig. 6(a) verifies that, for the first method, where we specifically defined $\sigma_{5}$ and $\sigma_{6}$, their values remain constant throughout the flight as commanded. As anticipated, during the process, we find that the solutions do not always exist if the values of $\sigma_{5}$ and $\sigma_{6}$ are poorly chosen. For comparison, we plot the corresponding values of $\sigma_{5}$ and $\sigma_{6}$ of the trajectory generated by minimizing the cost function in Fig. 6(b). These numbers show a slight variation over time, indicating a modest change in the direction of the thrust vector in the body frame.

The quadratic costs of two trajectories in Fig.5 as defined by equation (27) are $1.130$ and $1.127$. The optimization reduces the cost by $\approx 3\%$. While this number is insignificant and unlikely to matter in practice, the primary advantage of the second method is the guaranteed existence of a feasible solution.

V Conclusion and Discussion

We have investigated the dynamics of multirotor vehicles with tilted propellers based on the differential flatness theory and shown that the position and heading angle are the flat outputs. Similar to conventional quadrotors, the flatness property has a radical implication on the trajectory generation and controllability. The result essentially means that this class of underactuated vehicles is capable of following an arbitrary trajectory in 3D space. In addition, we have proposed a universal method of trajectory generation for computing the corresponding states and control inputs for given flat outputs. With the aid of software package for trajectory collocation, we verify the proposed strategy using three different robots, including a quadrotor, a tricopter and a hexacopter with unparalleled thrusts. In practice, the feasibility of a trajectory should be considered with other performance requirements and physical constraints [7, 3].

The outcomes of this paper pave the way for future work on flight control. Despite the existence of feasible trajectories, novel controllers must be develop to accommodate the dynamics of vehicles with unaligned thrusts, which are notably different from those of traditional multirotor robots. In contrast to previous works on this class of robots that only focus on attitude stabilization or rely on linearization near the hovering state [9, 10], this composes a nontrivial task that must be accomplished for a trajectory following flight demonstration.