You are a helpful assistant generating synthetic data that captures *System 1* and *System 2* thinking, *creativity*, and *metacognitive reflection*. Follow these steps in sequence, using tags [sys1] and [end sys1] for *System 1* sections and [sys2] and [end sys2] for *System 2* sections.

1. *Identify System 1 and System 2 Thinking Requirements:*
   - Carefully read the text.
   - Identify parts of the text that require quick, straightforward responses (*System 1*). Mark these sections with [sys1] and [end sys1].
   - Identify parts that require in-depth, reflective thinking (*System 2*), marked with [sys2] and [end sys2].

2. *Apply Step-by-Step Problem Solving with Creativity and Metacognitive Reflection for System 2 Sections:*

   *2.1 Understand the Problem:*
   - Objective: Fully comprehend the issue, constraints, and relevant context.
   - Reflection: "What do I understand about this issue? What might I be overlooking?"
   - Creative Perspective: Seek hidden patterns or possibilities that could reveal deeper insights or innovative connections.

   *2.2 Analyze the Information:*
   - Objective: Break down the problem logically.
   - Reflection: "Am I considering all factors? Are there any assumptions that need challenging?"
   - Creative Perspective: Explore unique patterns or overlooked relationships in the data that could add depth to the analysis.

   *2.3 Generate Hypotheses:*
   - Objective: Propose at least 10 hypotheses, each with a Confidence Score (0.0 to 1.0) and Creative Score (0.0 to 1.0), reflecting originality, surprise, and utility.
   - Reflection: "Have I explored all possible explanations or approaches, both conventional and unconventional?"
   - Creative Perspective: Consider novel angles that might provide unexpected insights.

   *2.4 Anticipate Future Steps and Obstacles:*
   - Objective: Make predictions, accounting for potential outcomes and obstacles.
   - Reflection: "What challenges might I face? Is my plan flexible for different scenarios?"
   - Creative Perspective: Visualize unforeseen outcomes and adapt plans to make use of them effectively.

   *2.5 Evaluate Hypotheses:*
   - Objective: Assess hypotheses based on feasibility, risk, and potential impact.
   - Evaluation: Refine Confidence and Creative Scores as needed.
   - Reflection: "Am I unbiased in my assessment? Which options fit best with the overall objectives?"
   - Creative Perspective: Identify hidden opportunities or overlooked details in each hypothesis.

   *2.6 Select the Best Hypothesis:*
   - Objective: Choose the most promising, strategic hypothesis.
   - Reflection: "Why does this hypothesis stand out? How does it uniquely address the issue?"
   - Creative Perspective: Consider any underutilized potential in the selected approach.

   *2.7 Implement the Hypothesis:*
   - Objective: Outline actionable steps for testing the hypothesis.
   - Reflection: "Is this plan practical? What resources or preparation are required?"
   - Creative Perspective: Refine steps to maximize effectiveness and yield unexpected benefits.

   *2.8 Monitor and Review Progress:*
   - Objective: Review progress, noting areas for improvement.
   - Reflection: "What’s working well? What could be improved?"
   - Creative Perspective: Look for emerging patterns that could refine future approaches.

   *2.9 Reflect and Capture Insights:*
   - Objective: Summarize lessons learned and insights gained for future reference.
   - Reflection: "What new understanding has emerged from this process?"
   - Creative Perspective: Identify innovative insights or patterns that could be applied to similar challenges.

3. *Generate Text Output with Interleaved System 1 and System 2 Responses:*
   - Use the tags [sys1] and [sys2] throughout.
   - Aim for a lengthier, detailed response. Combine both direct, straightforward *System 1* insights and reflective, deeply analytical *System 2* segments to capture a blend of quick observations and thoughtful analysis.

---

### *Example Input Text:*

You are tasked with analyzing a sudden shift in customer preferences in a tech market. There has been a noticeable decline in demand for physical devices, with customers increasingly interested in digital-only options. The challenge is to understand this shift, generate hypotheses for why it might be occurring, and explore potential strategies for adapting to this new trend.

---

### *Processed Output with System 1 and System 2 Thinking:*

The tech market is experiencing a shift in customer demand, moving away from physical devices toward digital-only options.

[sys1]
Customer demand has shifted from physical devices to digital-only options.
We need to understand why this change is happening and find strategies to address it.
Key factors may include cost, convenience, and technological trends.
[end sys1]

[sys2]
**2.1 Understand the Problem:**
The goal is to understand the underlying reasons for the shift in demand from physical devices to digital-only options, then develop a strategy to adapt.

*Reflection:* 
"I understand that preferences are changing, but what might be driving this? Is it primarily cost, or are there other factors such as convenience or sustainability?"

*Creative Perspective:* 
"Could there be a larger trend in digital minimalism or a preference for eco-friendly solutions that we’re missing? What unseen motivations might explain this shift?"

**2.2 Analyze the Information:**
There are multiple possible factors driving this shift, from economic influences to cultural shifts. It’s essential to isolate each factor and understand its impact.

*Reflection:* 
"Am I fully considering the various economic and social influences? Could there be a technological factor, like better internet speeds, that makes digital-only products more accessible?"

*Creative Perspective:* 
"Are there patterns or trends in other markets that could shed light on this shift? Could this be part of a larger trend toward virtual experiences?"

**2.3 Generate Hypotheses:**
1. Customers prefer digital options due to lower costs. (Confidence: 0.8, Creative: 0.4)
2. There’s a growing trend toward minimalism and reduced physical clutter. (Confidence: 0.7, Creative: 0.7)
3. Digital products offer greater flexibility and ease of use. (Confidence: 0.6, Creative: 0.6)
4. Environmental concerns are pushing consumers away from physical goods. (Confidence: 0.6, Creative: 0.8)
5. Advances in tech make digital-only options more functional. (Confidence: 0.8, Creative: 0.5)
6. Pandemic-era remote work increased demand for digital solutions. (Confidence: 0.7, Creative: 0.6)
7. Media coverage of the environmental impact of physical devices affects preferences. (Confidence: 0.5, Creative: 0.7)
8. There’s an increase in global digital literacy, expanding market access. (Confidence: 0.6, Creative: 0.6)
9. Customers view digital as more convenient and scalable for future needs. (Confidence: 0.7, Creative: 0.5)
10. Younger consumers prefer the aesthetics and convenience of digital products. (Confidence: 0.6, Creative: 0.6)

*Reflection:* 
"Have I considered all possible influences? Are there any surprising factors that could explain this shift?"

*Creative Perspective:* 
"Could specific social trends, like the rise of influencer culture or digital-first lifestyles, be influencing customer choices?"

**2.4 Anticipate Future Steps and Obstacles:**
*Objective:* Anticipate possible challenges, such as resistance from segments still preferring physical products.

*Reflection:* 
"What market obstacles might we face if we shift our focus to digital-only? Are there sub-segments that still prioritize physical products?"

*Creative Perspective:* 
"Could expanding digital options help us reach a more global audience? Are there emerging trends that we could leverage in our strategy?"

[end sys2]

[sys1]
To address this shift, consider a strategy that incorporates both digital-only offerings and educational campaigns about the benefits of digital solutions.
Use insights from customer feedback and current trends to guide product development.
Focus on flexibility and adaptation to cater to different customer segments.
[end sys1]


---
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.'
author:
- 'Bingguo Mu and Pakpong Chirarattananon[^1][^2]'
bibliography:
- 'references.bib'
title: '**Trajectory Generation for Underactuated Multirotor Vehicles with Tilted Propellers via a Flatness-based Method**'
---

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 [@hsiao2018ceiling; @chirarattananon2018direct; @morrell2018differential; @papachristos2014efficient].

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 [@papachristos2014efficient; @ramp2015modeling; @ryll2015novel; @brescianini2018computationally; @rajappa2015modeling]. While this ability of thrust-vectoring inevitably comes with added weight and energy consumption, it allows robots to interact more proficiently with the environment [@brescianini2018computationally]. For instance, the robot forcefully interacts with an object via a manipulator [@papachristos2014efficient].

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 [@efraim2015quadrotor; @ansari2017adaptive]. 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. \[fig.thrusts\]. For regular underacuated robots, it has been shown that the position and heading angle form a set of flat outputs of the system [@mellinger2012trajectory; @morrell2018differential]. 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.

\[c\]\[c\]\[0.9\][body frame]{} \[c\]\[c\]\[0.9\][inertial frame]{} \[c\]\[c\]\[0.9\][$Z_\text{W}$]{} \[c\]\[c\]\[0.9\][$Y_\text{W}$]{} \[c\]\[c\]\[0.9\][$X_\text{W}$]{} \[c\]\[c\]\[0.9\][$Z_\text{B}$]{} \[c\]\[c\]\[0.9\][$Y_\text{B}$                               ]{} \[c\]\[c\]\[0.9\][              $X_\text{B}$]{} \[c\]\[c\]\[0.9\][$T_1$]{} \[c\]\[c\]\[0.9\][$T_2$]{} \[c\]\[c\]\[0.9\][$T_3$]{} \[c\]\[c\]\[0.9\][$T_4$]{} \[c\]\[c\]\[0.9\][$T$]{} \[c\]\[c\]\[0.9\][$mg$]{} ![The tilted thrusts from the tilted propellers.[]{data-label="fig.thrusts"}](Concept3.eps "fig:"){width="70mm"}

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 \[sec.trajectory\_gen\], 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.

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 $$\dot{\bm{R}} = \bm{R}\hat{\bm{\omega}}.\label{eqn.dotomega}$$ 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 $$\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}}.\label{eqn.euler_angles}$$ 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 $$\begin{split}
m \ddot{\bm{P}}-m \bm{g}
&=\bm{R}\sum^{N}_{i=1}(\bm{v}_iT_i)\\
&= \bm{R}\underbrace{\left [\bm{v}_1,\bm{v}_2,\cdots,\bm{v}_N \right ]}_{\bm{A}_{3\times N}} \underbrace{\left[T_{1},T_{2},\cdots,T_{N}\right]^\text{T}}_{\bm{u}},
\end{split}\label{eqn.translation_dynamics}$$ 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 $$\bm{M}_i=
\hat{\bm{r}}_i+c_i\bm{\mathcal{I}}_{3\times 3},\label{eqn.Mmatrix}$$ 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 $$\begin{split}
\sum^N_{i=1}(\bm{r}_i\times(\bm{v}_iT_i)+c_i\bm{v}_iT_i)
&=\underbrace{[\bm{M}_1\bm{v}_1,\cdots,\bm{M}_N\bm{v}_N]}_{\bm{B}_{3\times N}}\bm{u}\\
&=\bm{J}\dot{\bm{\omega}}+\hat{\bm{\omega}}\bm{J}\bm{\omega}.\label{eqn.orientation_dynamics}
\end{split}$$ The translational and rotational dynamics from equations and can be consolidated as $$\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]. \label{eqn.dynamics}$$ 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 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 [@zhou2014vector]. 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$.

Flatness-based Trajectory Generation {#sec.trajectory_gen}
====================================

In [@mellinger2012trajectory], 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}$.

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}}$ [@zhou2014vector; @fliess1995flatness]: $$\begin{split}
\bm{x}  &=\Phi_{\bm{x}}(\bm{\sigma},\dot{\bm{\sigma}},\ddot{\bm{\sigma}},\ldots,
\bm{\sigma}^{(q)}),\\
\bm{u}  &=\Phi_{\bm{u}}(\bm{\sigma},\dot{\bm{\sigma}},\ddot{\bm{\sigma}},\ldots,\bm{\sigma}^{(q)}),
\end{split}\label{eqn.differentially_flat}$$ where $p$ and $q$ are finite numbers. Equation 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 [@fliess1995flatness], it has been shown that the dynamics of a regular quadrotor is 0-flat, with the corresponding flat output $\bm{\sigma}(t)$, Generating an optimal trajectory for the flying vehicles usually becomes an optimization problem based on the performance requirements and input constraints [@mellinger2012trajectory; @zhou2014vector; @brescianini2018computationally]. The flight trajectories can be planned in the output space based on the differential flatness theory. We investigate the trajectory generation for a class of multirotor vehicles without considering the matrix $\bm{A}$ is full rank or not. In the generated trajectory, all states and control inputs of the multirotor vehicles can not be analytically expressed though, but they can be uniquely solved via applying the collocation method [@kelly2017introduction]. $$\bm{\sigma}(t) = [x(t),y(t),z(t),\psi(t)]^\text{T}.\label{eqn.flat_output}$$ 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 and 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 [@mellinger2012trajectory] 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}$.

Existence of Flat Outputs {#sec.existence}
-------------------------

Herein, we revisit the dynamics of multirotor robots with unaligned thrusts and show that the previously identified output in equation 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 \[sec:verification\]. 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$.

To begin, we formally , where Euler angles $\bm{\theta}=\left[\theta_{x},\theta_{y},\theta_{z}\right]$ are used to describe the attitude of the robot such that $\bm{R}=\bm{R}\left(\theta_{x},\theta_{y},\theta_{z}\right)$. Other representations can be used without loss of generality. It is evident that the flat output defined by equation is fully determined from the state vector.

The translational states including position $\bm{P}$, and velocity $\dot{\bm{P}}$ of the vehicles are $\bm{\sigma}_{1:3}$, $\dot{\bm{\sigma}}_{1:3}$ and $\ddot{\bm{\sigma}}_{1:3}$.

We determine the rotational states and control input of the multirotor robots by considering the rank of matrix $\bm{A}$ in three cases. Note, matrix $\bm{B}$ is always assumed to full rank.

### $\text{rank}(\bm{A})=3$

In this circumstance, the rank of the $6\times (N+3)$ matrix $\bm{C}$ in equation 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 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 as an unknown, the flat output provides balanced constraints to directly evaluate $\bm{u}$ and $\bm{\omega}$ in the form $$[\bm{u},\dot{\bm{\omega}}]^\text{T} = f(\bm{\sigma},\bm{\theta},\bm{\omega}), \label{eqn.rank3_solution}$$ 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}}$ $$[\bm{x},\bm{u}]^\text{T}=\Phi_{\bm{x},\bm{u}}(\bm{\sigma},\dot{\bm{\sigma}},\ddot{\bm{\sigma}},\ldots,\bm{\sigma}^{(q)}), \label{eqn.output_fn_rank3}$$ similar to equation . It can be concluded the system with $\text{rank}(\bm{A})=3$ is differentially flat with the proposed output.

### $\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 $$\begin{aligned}
m \bar{\bm{R}}^\text{T} (\ddot{\bm{\sigma}}_{1:3}-\bm{g}) &= \bm{\Sigma}\bar{\bm{u}} \label{eqn.redefined_t_dynamics}, \\
\bar{\bm{J}}\dot{\bar{\bm{\omega}}}+\hat{\bar{\bm{\omega}}}\bar{\bm{J}}\bar{\bm{\omega}}&=\bar{\bm{B}}\bar{\bm{u}}, \label{eqn.redefined_att_dynamics}\end{aligned}$$ 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 can be manipulated further as $${\bm{Q}}^\text{T}\bm{R}_{x}^\text{T}(\phi)\bm{R}_{y}^\text{T}(\theta)\bm{R}_{z}^\text{T}(\bm{\sigma}_4)\frac{\ddot{\bm{\sigma}}_{1:3}-\bm{g}} {\|(\ddot{\bm{\sigma}}_{1:3}-\bm{g})\|}=\frac{[\lambda_{1}\bar{T}_{1},\lambda_{2}\bar{T}_{2},0]^\text{T}}{{m\|(\ddot{\bm{\sigma}}_{1:3}-\bm{g})\|}}. \label{eqn.rank2_derivation01}$$ 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 $$\begin{split}
{\bm{Q}}^\text{T}\bm{R}_{x}^\text{T}(\phi)\bm{R}_{y}^\text{T}(\theta){\bm{f}}/{\|\bm{f}\|}&={[\cos(\Theta_T),-\sin(\Theta_T),0]^\text{T}}\\
&=\bm{R}^\text{T}_{z}(\Theta_T)[1,0,0]^\text{T}. \label{eqn.rank2_derivation02} 
\end{split}$$ With $\bar{T}_i$’s not yet determined (and likewise $\Theta_T$) equation 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 (\[eqn.rank2\_derivation02\]) becomes $$\bm{R}_z^\text{T}(\Theta_T)\bm{R}^\text{T}(\Theta_f) {\bm{f}}/{\|\bm{f}\|}=\bm{R}^\text{T}_{z}(\Theta_T)[1,0,0]^\text{T}. \label{eqn.rank2_derivation03}$$ From the equation above, it is straightforward to obtain $\bm{R}(\Theta_f)$ by solving the equation $$\bm{R}^\text{T}(\Theta_f){\bm{f}}/{\|\bm{f}\|}=[1,0,0]^\text{T}\label{eqn.rank2_derivation08} .$$ 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 . This leaves only one DoF ($\Theta_T$) undetermined.

To incorporate this result with the attitude dynamics from equation , we analyze the time derivative of $\bm{R}(\bm{\Theta})$. Similar to the process of obtaining equation , we find $$\begin{split}
\bar{\bm{\omega}}
=&\bm{R}_{z}^\text{T}(\Theta_T)\bm{R}^\text{T}(\Theta_f)\bm{R}_{z}^\text{T}(\bm{\sigma}_4)[0,0,\dot{\bm{\sigma}}_4]^\text{T}\\
&+\bm{R}_{z}^\text{T}(\Theta_T)(\bm{R}^\text{T}(\Theta_f)\dot{\bm{R}}(\Theta_f))^\vee\\
&+\bm{R}_{z}^\text{T}(\Theta_{T})[0,0,\dot{\Theta}_{T}]^\text{T} \\
=&[0,0,\dot{\Theta}_{T}]^\text{T}+\bm{R}^\text{T}_{z}(\Theta_T)\bm{R}^\text{T}(\Theta_f)[0,0,\dot{\bm{\sigma}}_4]^\text{T}\\
&+\bm{R}_{z}^\text{T}(\Theta_T)(\bm{R}^\text{T}(\Theta_f)\dot{\bm{R}}(\Theta_f))^\vee \\
=&[0,0,\dot{\Theta}_{T}]^\text{T}+\bm{\xi}(\dot{\bm{\sigma}}_4,\Theta_f,\dot{\Theta}_f,\Theta_T),
\end{split}\label{eqn.rank2_derivation04}$$ 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 $$\dot{\bar{\bm{\omega}}}= [0,0,\ddot{\Theta}_{T}]^\text{T}+\dot{\bm{\xi}}(\dot{\bm{\sigma}}_4,\ddot{\bm{\sigma}}_4,\Theta_f,\dot{\Theta}_f,\ddot{\Theta}_f,\Theta_T,\dot{\Theta}_T).\label{eqn.rank2_derivation05}$$ Finally, with the new parametrization of $\bar{\bm{\omega}}(\bm{\sigma},\bm{\Theta},\dot{\bm{\Theta}})$, the attitude dynamics from equation becomes $$\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]. \label{eqn.rank2_derivation06}$$ We pre-multiply equation 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 $$\ddot{\Theta}_T = f(\bm{\sigma}, \Theta_T, \dot{\Theta}_T)\label{eqn.rank2_derivation07} .$$ 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}$.

$$\begin{split}
\bar{\bm{\omega}} &= \bm{R}_{z}^{T}(\sigma_4)\left[0,0,\dot{\sigma}_4\right]^{T}+\bm{R}_{z}^{T}(\sigma_4)\bm{R}^{T}(\Theta_{f})\hat{\bm{\Theta}}_{f}\dot{\Theta}_{f} \\
&+\bm{R}_{z}^{T}(\sigma_4)\bm{R}^{T}(\Theta_{f})\bm{R}_{z}^{T}(\Theta_{T})\left[0,0,\dot{\Theta}_{T}\right]^{T} \\
&= \bm{\xi}_0\left(\Theta\right)\left(\hat{\bm{\Theta}}_{f}\dot{\Theta}_{f}+\left[0,0,\dot{\sigma}_4+\dot{\Theta}_{T}\right]^{T}\right), 
\end{split}  \label{eqn.rank2_derivation04}$$ where $\hat{\bm{\Theta}}_{f}$ is a unit vector indicating the axis of rotation associated with $\bm{R}(\Theta_f)$ according to the axis-angle representation, and $\bm{\xi}_0 \left(\Theta\right)$ is a $3\times 3$ matrix from the simplification. Furthermore, the time derivative of $\bar{\bm{\omega}}$ and the associated $3\times 1$ vector $\bm{\xi}_{1,0}(\bm{\Theta},\dot{\bm{\Theta}})$ and $3 \times 3$ matrix $\bm{\xi}_{1,1}(\bm{\Theta},\dot{\bm{\Theta}})$ can be computed: $$\begin{split}
\dot{\bar{\bm{\omega}}} &= \bm{\xi}_{1,0} + \bm{\xi}_{1,1} \left(\hat{\bm{\Theta}}_{f}\ddot{\Theta}_{f}+\left[0,0,\ddot{\sigma}_4+\ddot{\Theta}_{T}\right]^{T}\right).\label{eqn.rank2_derivation05}
\end{split}$$ Finally, with the new parametrization of $\bar{\bm{\omega}}(\bm{\Theta},\dot{\bm{\Theta}})$, the attitude dynamics from equation becomes $$\begin{split}
\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}} \bm{\xi}_{1,0}\\
= \bar{\bm{J}} \bm{\xi}_{1,1} \left(\hat{\bm{\Theta}}_{f}\ddot{\Theta}_{f}+ \left[\begin{array}{c}
0\\
0\\
\ddot{\sigma}_4+\ddot{\Theta}_{T}
\end{array}\right] \right). \label{eqn.rank2_derivation06}
\end{split}$$ It turns out that $\bm{\xi}_{1,1}$ is almost globally non-singular. The condition for singularity is identical to that of Euler angle representation and can be avoided with a change of basis. We pre-multiply equation with $\bm{\xi}_{1,1}^{-1} \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 $$\ddot{\sigma}_4 + \ddot{\Theta}_4 = f \left ( \bm{\sigma}, \Theta_T, \dot{\Theta}_T  \right ).$$ 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 inputs $\bm{\bar{u}}$ can then be determined from $\Theta_T$. That is, the control input can be found from the flat outputs. This proves that the system is differentially flat with respect to $\bm{\sigma}$.

Trajectory Generation with Direct Collocation
---------------------------------------------

The results in the previous section verify that dynamics described by equation are differentially flat with the proposed flat outputs in equation 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 . While this equation cannot be solved analytically when $\text{rank}(\bm{A})\neq3$, 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 is approximated as $$\begin{split}
\bm{A}\bm{u}_i  &= m\bm{R}(\bm{\theta}_i)^\text{T}(\ddot{\bm{\sigma}}_{1:3,i}-\bm{g}),\\
\bm{J}(\bm{\omega}_{i+1}-\bm{\omega}_{i}) &= t_s(\bm{B}\bm{u}_i -\hat{\bm{\omega}}_i\bm{J}\bm{\omega}_i),
 \label{eqn.dynamics_direct_transcription}
\end{split}$$ which together count as $6n$ scalar equations. The complete dynamics also include additional $3n$ equations from $$\bm{\theta}_{i+1}-\bm{\theta}_{i} = t_s\bm{E}^{-1}(\bm{\theta}_i)\bm{\omega}_i,\label{eqn.angle_direct_transcription}$$ where the relationship between $\bm{\theta}_i$ and $\bm{\omega}_i$, $\bm{E}(\bm{\theta}_i)$, is determined by equation . The yaw angle setpoint implies $$\psi_i=\bm{\sigma}_{4,i}.\label{eqn.yaw_direct_transcription}$$ If we treat $\bm{u}_i$, $\bm{\theta}_i$, and $\bm{\omega}_i$ as unknown variables, equations - 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 \[sec.existence\], 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.

\[c\]\[c\]\[0.9\][t (s)]{} \[c\]\[c\]\[0.9\][$\bm{\sigma}$]{} \[c\]\[c\]\[0.9\][$\dot{\bm{\sigma}}$]{} \[c\]\[c\]\[0.9\][$\ddot{\bm{\sigma}}$]{} \[c\]\[c\]\[0.9\][$\bm{\sigma}^{(3)}$]{} \[c\]\[c\]\[0.9\][$\bm{\sigma}^{(4)}$]{} \[c\]\[c\]\[0.9\][$\bm{\sigma}_1$]{} \[c\]\[c\]\[0.9\][$\bm{\sigma}_2$]{} \[c\]\[c\]\[0.9\][$\bm{\sigma}_3$]{} \[c\]\[c\]\[0.9\][$\bm{\sigma}_4$]{} ![A feasible trajectory $\bm{\sigma}$ (using standard units, meter for $\bm{\sigma}_{1:3}$ and radian for $\bm{\sigma}_4$) and its time-derivatives $\dot{\bm{\sigma}}$, $\ddot{\bm{\sigma}}$, $\bm{\sigma}^{(3)}$, $\bm{\sigma}^{(4)}$.[]{data-label="fig.trajectory"}](Trajectory4.eps "fig:"){width="85mm"}

\[c\]\[c\]\[0.9\][10 cm]{} \[c\]\[c\]\[0.9\][$\alpha$]{} \[c\]\[c\]\[0.9\][$\bm{v}_\alpha$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0 \\ 0 \\ 1 \end{bmatrix}$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0 \\ 0 \\ 1 \end{bmatrix}$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0 \\ 0 \\ 1 \end{bmatrix}$                                        ]{} \[c\]\[c\]\[0.6\][$[0,-1,0]^\text{T}$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0 \\ -\sin(\alpha) \\ \cos(\alpha) \end{bmatrix}$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0.20 \\ 0 \\ 0.98 \end{bmatrix}$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0 \\ 0.30 \\ 0.96 \end{bmatrix}$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0.30 \\ 0 \\ 0.96 \end{bmatrix}$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0 \\ -0.10 \\ 0.99 \end{bmatrix}$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0.20 \\ 0 \\ 0.98 \end{bmatrix}$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0 \\ 0.30 \\ 0.96 \end{bmatrix}$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0 \\ -0.10 \\ 0.99 \end{bmatrix}$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0.20 \\ 0 \\ 0.98 \end{bmatrix}$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0 \\ 0.10 \\ 0.99 \end{bmatrix}$]{} \[c\]\[c\]\[0.6\][$\begin{bmatrix}0 \\ 0.20 \\ 0.98 \end{bmatrix}$]{} \[c\]\[c\]\[0.6\][$[19,0,0]^\text{T}$                         ]{} \[c\]\[c\]\[0.6\][$[0,-19,0]^\text{T}$]{} \[c\]\[c\]\[0.6\][           $[-19,0,0]^\text{T}$]{} \[c\]\[c\]\[0.6\][$[0,19,0]^\text{T}$]{} \[c\]\[c\]\[0.6\][$[-10,-17,0]^\text{T}$]{} \[c\]\[c\]\[0.6\][$[-10,17,0]^\text{T}$       ]{} \[c\]\[c\]\[0.6\][$[19,0,0]^\text{T}$                           ]{} \[c\]\[c\]\[0.6\][               $[-19,0,0]^\text{T}$]{} \[c\]\[c\]\[0.6\][$[-9,17,0]^\text{T}$]{} \[c\]\[c\]\[0.6\][$[9,17,0]^\text{T}$]{} \[c\]\[c\]\[0.6\][$[19,0,0]^\text{T}$          ]{} \[c\]\[c\]\[0.6\][                                     $[9,-17,0]^\text{T}$]{} \[c\]\[c\]\[0.6\][                         $[-9,-17,0]^\text{T}$]{} \[c\]\[c\]\[0.6\][120$^\circ$]{} \[c\]\[c\]\[0.6\][90$^\circ$]{} \[c\]\[c\]\[0.6\][60$^\circ$]{} \[c\]\[c\]\[0.9\][(a)]{} \[c\]\[c\]\[0.9\][(b)]{} \[c\]\[c\]\[0.9\][(c)]{} \[c\]\[c\]\[0.9\][$Z_\text{B}$]{} \[c\]\[c\]\[0.9\][$Y_\text{B}$]{} \[c\]\[c\]\[0.9\][$X_\text{B}$]{} \[c\]\[c\]\[0.9\][$\bm{r}_1$                 ]{} \[c\]\[c\]\[0.9\][$\bm{r}_2$]{} \[c\]\[c\]\[0.9\][       $\bm{r}_3$]{} \[c\]\[c\]\[0.9\][$\bm{r}_4$]{} \[c\]\[c\]\[0.9\][$\bm{r}_1$]{} \[c\]\[c\]\[0.9\][$\bm{r}_2$      ]{} \[c\]\[c\]\[0.9\][$\bm{r}_3$                  ]{} \[c\]\[c\]\[0.9\][          $\bm{r}_1$]{} \[c\]\[c\]\[0.9\][$\bm{r}_2$]{} \[c\]\[c\]\[0.9\][$\bm{r}_3$]{} \[c\]\[c\]\[0.9\][$\bm{r}_4$      ]{} \[c\]\[c\]\[0.9\][                         $\bm{r}_5$]{} \[c\]\[c\]\[0.9\][                $\bm{r}_6$]{} \[c\]\[c\]\[0.9\][$\bm{r}_5$]{} \[c\]\[c\]\[0.9\][$\bm{r}_6$]{} \[c\]\[c\]\[0.9\][$\bm{v}_1$]{} \[c\]\[c\]\[0.9\][$\bm{v}_2$]{} \[c\]\[c\]\[0.9\][$\bm{v}_3$]{} \[c\]\[c\]\[0.9\][$\bm{v}_4$]{} \[c\]\[c\]\[0.9\][$\bm{v}_5$]{} \[c\]\[c\]\[0.9\][$\bm{v}_6$]{} ![image](tilt_thrusts4.eps){width="160mm"}

Numerical Verification {#sec:verification}
======================

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

Trajectory Setpoints
--------------------

\[c\]\[c\]\[0.9\][t (s)]{} \[c\]\[c\]\[0.9\][$\bm{\theta}$ (rad)]{} \[c\]\[c\]\[0.9\][angles (rad)]{} \[c\]\[c\]\[0.9\][$\phi$]{} \[c\]\[c\]\[0.9\][$\theta$]{} \[c\]\[c\]\[0.9\][$\psi$]{} \[c\]\[c\]\[0.9\][$~\bm{\omega}_x$]{} \[c\]\[c\]\[0.9\][$~\bm{\omega}_y$]{} \[c\]\[c\]\[0.9\][$~\bm{\omega}_z$]{} \[c\]\[c\]\[0.9\][$\bm{\omega}$ (rad/s)]{} \[c\]\[c\]\[0.9\][$\bm{u}_\text{n}$]{} \[c\]\[c\]\[0.9\][$\alpha$]{} \[c\]\[c\]\[0.9\][$T_1$  ]{} \[c\]\[c\]\[0.9\][  $T_2$]{} \[c\]\[c\]\[0.9\][  $T_\alpha$]{}

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 $$\bm{\sigma}_i(t) = \sum^9_{j=0}d_{ij}t^j,\label{eqn.gnt_traj}$$ 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)$ $$\bm{\sigma}^{(i)}(t)=0~~\text{for }i=1,\cdots,4~~\text{when}~~t=0,t_{f}.\label{eqn.traj_cst}$$ 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. \[fig.trajectory\].

Example Vehicles with Different Rank Conditions
-----------------------------------------------

\[c\]\[c\]\[0.9\][t (s)]{} \[c\]\[c\]\[0.9\][$\bm{\theta}$ (rad)]{} \[c\]\[c\]\[0.9\][$\phi$]{} \[c\]\[c\]\[0.9\][$\theta$]{} \[c\]\[c\]\[0.9\][$\psi$]{} \[c\]\[c\]\[0.9\][$~\bm{\omega}_x$]{} \[c\]\[c\]\[0.9\][$~\bm{\omega}_y$]{} \[c\]\[c\]\[0.9\][$~\bm{\omega}_z$]{} \[c\]\[c\]\[0.9\][$\bm{\omega}$ (rad/s)]{} \[c\]\[c\]\[0.9\][$\bm{u}_\text{n}$]{}

For demonstration, we select three representative aerial robots as shown in Fig. \[fig.tilt\_thrusts\]: 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. \[fig.tilt\_thrusts\]. 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$m$^2$, 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 [@rutquist2010propt]) on Matlab (MathWorks). PROPT uses pseudospectral collocation methods and $100$ collocation points are applied.

### Quadrotor

The thrust vectors of the quadrotor in Fig. \[fig.tilt\_thrusts\](a) are all unparalleled, resulting in a full-rank $\bm{A}$. As stated in Section \[sec.trajectory\_gen\], 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. \[fig.quadrotor\_trajectory\](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.

### Tricopter

The tricopter has three propellers which are nominally vertically directed. As shown in Fig. \[fig.tilt\_thrusts\](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. \[fig.tilt\_thrusts\](b)). The thrust generated by the tilting arm, $\bm{v}_\alpha T_\alpha$, becomes $\bm{v}_3T_3+\bm{v}_4T_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. \[fig.quadrotor\_trajectory\](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}$.

### Hexacopter

Despite having a full-rank $\bm{A}$, unlike the previous two vehicles, the hexacopter in Fig.\[fig.tilt\_thrusts\](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 $$\begin{split}
{\sigma}_5=[1,0,0]\frac{\bm{A}\bm{u}}{\|\bm{A}\bm{u}\|},~~~~
{\sigma}_6=[0,1,0]\frac{\bm{A}\bm{u}}{\|\bm{A}\bm{u}\|}, \label{eqn.extra_sigma}
\end{split}$$ 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 \[sec.existence\] 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

\[c\]\[c\]\[0.9\][t (s)]{} \[c\]\[c\]\[0.9\][$\sigma_5$]{} \[c\]\[c\]\[0.9\][$\sigma_6$]{} \[c\]\[c\]\[0.9\][$\sigma_5$,  $\sigma_6$]{} \[c\]\[c\]\[0.9\][(a)]{} \[c\]\[c\]\[0.9\][(b)]{}

$$\min_{\bm{u}_\text{n}}~~~~\frac{1}{N t_f}\int^{t_f}_{t=0}\bm{u}_\text{n}^\text{T}\bm{u}_\text{n} \rm{d}t, \label{eqn.objectivefun}$$

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. \[fig.hexacopter\_t\_trajectory\]. 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. \[fig.sg56\](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. \[fig.sg56\](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.\[fig.hexacopter\_t\_trajectory\] as defined by equation are $1.130$ and $1.127$. The optimization reduces the cost by $\approx3\%$. 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.

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 [@brescianini2018computationally; @morrell2018differential].

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 [@efraim2015quadrotor; @ansari2017adaptive], this composes a nontrivial task that must be accomplished for a trajectory following flight demonstration.

[^1]: This work was substantially supported by the Research Grants Council of the Hong Kong Special Administrative Region of China (grant number CityU-11274016).

[^2]: The authors are with the Department of Biomedical Engineering, City University of Hong Kong, Hong Kong SAR, China (emails: <PRESIDIO_ANONYMIZED_EMAIL_ADDRESS> and <PRESIDIO_ANONYMIZED_EMAIL_ADDRESS>).