Differential Flatness of Lifting-Wing Quadcopters Subject to Drag and Lift for Accurate Tracking

The coordinate frames of the lifting-wing quadcopter are illustrated in Fig. 3. The earth-fixed north-east-down (NED) coordinate ${}^{\rm e}{\mathcal{F}}$, whose basis is composed of the column vectors of the identity matrix $\left[\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z}\right]$, is chosen as the inertial frame. The quadcopter-body coordinate frame ${}^{\rm b}\mathcal{F}$ is made of an orthonormal basis $\left\{\mathbf{x}_{\text{b}},\mathbf{y}_{\text{b}},\mathbf{z}_{\text{b}}\right\}$ represented in ${}^{\text{e}}\mathcal{F}$ with the origin located at the center of mass. ${}^{\rm b}\mathcal{F}$ rotates along $\mathbf{y}_{\text{b}}$ by $\kappa$ to obtain the lifting-wing coordinate frame ${}^{\rm l}\mathcal{F}$, which is made of orthonormal basis $\left\{\mathbf{x}_{\text{l}},\mathbf{y}_{\text{l}},\mathbf{z}_{\text{l}}\right\}$ represented in ${}^{\rm e}\mathcal{F}$. The vector $\mathbf{x}_{\text{l}}$ coincides with the chord line and the lifting wing symmetry plane. To handle the wing-induced forces, the wind coordinate frame ${}^{\rm w}\mathcal{F}$ is defined as $\left\{\mathbf{x}_{\text{w}},\mathbf{y}_{\text{w}},\mathbf{z}_{\text{w}}\right\}$ with the vector $\mathbf{x}_{\text{w}}$ coincides with the airspeed vector $\mathbf{v}_{\rm a}$. The rotation matrix ${\mathbf{R}}_{\text{b}}^{\text{e}}=\left[\mathbf{x}_{\text{b}}\ \mathbf{y}_{\text{b}}\ \mathbf{z}_{\text{b}}\right]$ gives the transformation from ${}^{\rm b}\mathcal{F}$ (indicated by the subscript b) to ${}^{\rm e}\mathcal{F}$ (indicated by the superscript e). Throughout this paper, the superscripts $\left({}^{\text{e,b,l,w}}\mathbf{x}\right)$ are utilized to specify in which coordinate frame a vector $\mathbf{x}$ is formulated. It should be noted that the superscripts are omitted for vectors in ${}^{\rm e}\mathcal{F}$. The position of the lifting-wing quadcopter is denoted as $\mathbf{p}$, and its derivatives, velocity, acceleration and jerk as $\mathbf{v}$, $\mathbf{a}$, $\mathbf{j}$, respectively.

The lifting-wing quadcopter is modeled as a rigid body with six degrees of freedom (6-DOF). The control inputs to the system are taken as the collective thrust $f_{z}$ and the angular velocity expressed in ${}^{\rm b}\mathcal{F}$ as ${}^{\rm{b}}{\bm{\omega}}$. It is assumed that the bandwidth of the angular velocity controller is high enough that the angular dynamics can be neglected. The motion of the aircraft is formulated by Newton and Euler’s Law of Motion [17]:

	$\displaystyle\dot{{\mathbf{p}}}$	$\displaystyle=\mathbf{v}$		(1)
	$\displaystyle\dot{{\mathbf{v}}}$	$\displaystyle=\mathbf{a}=\frac{\mathbf{f}}{m}$		(2)
	$\displaystyle{\dot{\mathbf{R}}}_{\text{b}}^{\text{e}}$	$\displaystyle={\mathbf{R}}_{\text{b}}^{\text{e}}{\left[{{}^{\text{b}}{\bm{\omega}}}\right]_{\times}}$		(3)

where $m$ is the mass of the aircraft, and ${\left[{{}^{\text{b}}{\bm{\omega}}}\right]_{\times}}$ denotes the skew-symmetric matrix. In addition, $\mathbf{f}$ represents forces acting on the airframe expressed in ${}^{\rm e}\mathcal{F}$. As shown in Fig. 3, the forces acting on the aircraft are

$${\mathbf{f}}={\mathbf{R}_{\rm{b}}^{\rm{e}}}\cdot{}^{\text{b}}{\mathbf{f}_{\text{r}}}+{\mathbf{R}_{\rm{l}}^{\rm{e}}}\cdot{}^{\rm l}{\mathbf{f}_{\text{a}}}+m{\mathbf{g}}$$

(4)

where ${}^{\text{b}}{\mathbf{f}_{\text{r}}}=\left[0\ 0\ f_{z}\right]^{\text{T}}$ with $f_{z}\leq 0$ is the collective thrust produced by rotors, ${}^{\rm l}{\mathbf{f}_{\text{a}}}$ is aerodynamic forces produced by the lifting wing, and ${\mathbf{g}}=\left[0\ 0\ g\right]^{\text{T}}$ is the gravitational acceleration vector with $g=9.81{\rm m/s}^{2}$.

A special parametric expression is used to model steady aerodynamic forces ${}^{\rm s}\mathbf{f}_{\rm a}$ [18]:

$${}^{\rm s}{\mathbf{f}}_{{\text{a}}}=\frac{1}{{2}}\rho S{V_{a}^{2}}\left[{\begin{array}[]{*{20}{c}}{{C_{{d_{0}}}}}+C_{L\alpha}\sin^{2}\alpha\ \ 0\ \ C_{L\alpha}\sin\alpha\cos\alpha\end{array}}\right]^{\rm T}$$

(5)

where $\rho$ is the air density, $S$ is the surface of the wing. The airspeed vector is defined as $\mathbf{v}_{\rm a}=\mathbf{v}-\mathbf{v}_{\rm w}$, $\mathbf{v}_{\rm w}$ is wind velocity expressed in ${}^{\rm e}\mathcal{F}$, $V_{a}=\lVert\mathbf{\bf v}_{\rm a}\rVert$ is the airspeed, $\alpha$ is the angle of attack; aerodynamic parameters $C_{d_{0}}$ is the minimum drag coefficient, $C_{y_{0}}$ is the minimum side force coefficient, $C_{l\alpha}$ is a coefficient related to lift, which are determined by the airfoil. This model is further developed and improved in [19], where a $\phi-$theory method is proposed and demonstrates that the model can cover the main dynamic characteristics in the full flight envelope.

When converted to ${}^{\rm e}{\mathcal{F}}$, ${}^{\rm s}\mathbf{f}_{\rm a}$ is an expression independent of the angle of attack:

$${}^{\rm l}{\mathbf{f}}_{{\text{a}}}=\frac{1}{{2}}\rho S{V_{a}}\left[{\begin{array}[]{*{20}{c}}{{C_{{d_{0}}}}}&0&0\\ 0&{{C_{{y_{0}}}}}&0\\ 0&0&C_{d_{0}}+C_{L\alpha}\end{array}}\right]{{\mathbf{R}}_{\rm b}^{\rm l}}\left({{\mathbf{R}}_{\rm b}^{\rm e}}\right)^{\rm T}\cdot{\mathbf{v}_{\rm a}}$$

(6)

Rotation matrix ${{\mathbf{R}}_{\rm b}^{\rm l}}$ transforms a vector from ${}^{\rm b}{\mathcal{F}}$ to ${}^{\rm l}{\mathcal{F}}$ given as

$${{\mathbf{R}}_{\rm b}^{\rm l}}=\left[{\begin{array}[]{*{20}{c}}\cos\kappa&0&-\sin\kappa\\ 0&1&0\\ \sin\kappa&0&\cos\kappa\end{array}}\right]$$

(7)

where $\kappa$ is the installation angle of the lifting wing. In practice applications, different installation angles can be designed according to specific flight indexes, such as flight speed and range. Generally, the installation angle is set at about 45 degrees to ensure good wind resistance [1] [20], as shown in Fig. 1(a) and (c). In particular, when the installation angle is set at 90 degrees, it is the so-called tail-sitter aircraft, as shown in Fig. 1(b) and (d). This type of aircraft is more similar to the fixed wing in cruise flight [2][7]. In this paper, we limit $\kappa\in(15^{\circ},90^{\circ}]$, because when $\kappa$ is small, additional propellers are usually added, which is a common dual-system aircraft [21].

The control design aims to accurately track the trajectory reference ${\mathbf{p}_{\rm ref}}$. The reference ${\mathbf{p}_{\rm ref}}$, which is at least 3-order continuous, may be provided by a pre-planned trajectory or an online motion planning algorithm.

To facilitate the problem description, a brief introduction to differential flattening is given. Considering the following dynamic system:

	$\displaystyle\dot{\bm{X}}$	$\displaystyle=\mathcal{A}\left(\bm{X},\bm{U}\right)$		(8)
	$\displaystyle{\bm{Y}}$	$\displaystyle=\mathcal{B}\left(\bm{X},\bm{U}\right),$		(9)

an essential property of flat systems is that their state ${\bm{X}}$ and input ${\bm{U}}$ can be directly written as algebraic functions of the flat output ${\bm{Y}}$ and a finite number of its time-derivatives. However, the input and state variables calculated from the flat output cannot be used to directly control the dynamic system (9) because of the inevitable model uncertainty. In practice, they can serve as feedforward control inputs to improve trajectory tracking performance by reducing the phase lag. At the same time, the model uncertainty can be compensated for by the feedback control.

In this paper, we choose ${\bm{Y}}=\mathbf{p}$ as the flat output, ${\bm{U}}=\left\{f_{z},{{}^{\rm b}{\bm{\omega}}}\right\}$ as inputs, $\bm{X}=\left\{\mathbf{p,v,a,R_{\rm b}^{\rm e}}\right\}$ as states, and $\bm{X}$ and $\bm{U}$ are constrained by dynamic equations (1)-(3). Then the objective of this paper is to address the following two problems :

• •

  Find the function of $\mathcal{H}$ and $\mathcal{G}$ that satisfy

  	$\displaystyle{\bm{X}}$	$\displaystyle=\mathcal{H}\left(\bm{Y},\dot{\bm{Y}},\ldots,\bm{Y}^{(n)}\right)$		(10)
  	$\displaystyle{\bm{U}}$	$\displaystyle=\mathcal{G}\left(\bm{Y},\dot{\bm{Y}},\ldots,\bm{Y}^{(m)}\right).$		(11)

• •

  Given a reference trajectory $\mathbf{p}_{\rm ref}(t)$, design a controller for the lifting-wing quadcopter with uncertainty such that $\mathop{\lim}\limits_{t\to\infty}\left\|{\mathbf{p}(t)-\mathbf{p}_{\rm ref}(t)}\right\|=0$ by using the differential flatness property given by (10) and (11).

The two problems will be solved in Section III and Section IV, respectively.

During a coordinated turn, there is no lateral acceleration expressed in ${{}^{\rm w}\mathcal{F}}$ of the aircraft, that is ${}^{\rm w}{\mathbf{a}}={{}^{\rm w}{\mathbf{f}}}/m=\left[a_{x_{\rm w}}\ 0\ a_{z_{\rm w}}\right]^{\rm T}$, as shown in Fig. 4. Then, the dynamic equation (2) is reformulated as

$${\mathbf{x}}_{\rm w}{a_{x_{\rm w}}}+{\mathbf{z}}_{\rm w}{a_{z_{\rm w}}}+{\mathbf{g}}-\mathbf{a}=0.$$

(12)

Since the vector $\mathbf{x}_{\rm w}$ coincides with the airspeed vector, $\mathbf{x}_{\rm w}$ is calculated as

$\displaystyle{\mathbf{x}}_{\rm w}=\mathbf{v}/\lVert\mathbf{v}\rVert$

(13)

when the wind velocity $\mathbf{v}_{\rm wind}=\bm{0}$ and the norm of the reference velocity $\lVert\mathbf{v}_{\rm{ref}}\rVert\neq 0$. In fact, if the wind velocity can be estimated online, the results of differential flatness hold true. For the convenience of the following derivation, the wind velocity is assumed to be 0. The case where $\lVert\mathbf{v}_{\rm{ref}}\rVert=0$ will be discussed later in section III-C.

Left-multiplying (12) by ${\mathbf{x}}_{\rm w}^{\rm T}$ yields

$$a_{x_{\rm w}}={\mathbf{x}}_{\rm{w}}^{\rm{T}}\left({{\mathbf{a}}-{\mathbf{g}}}\right).$$

(14)

In ${}^{\rm w}{\mathcal{F}}$, the reference acceleration is constrained in ${\mathbf{x}_{\rm w}}-{\mathbf{z}_{\rm w}}$ plane. Then, the component of the reference acceleration along $\mathbf{z}_{\rm w}$ can be derived from (12) by substituting (14) and applying the Euclidean norm

$\displaystyle a_{z_{\rm w}}$

$\displaystyle=-\lVert\left({{\mathbf{I}}-{{\mathbf{x}}_{\rm w}}{\mathbf{x}}_{\rm w}^{\rm T}}\right)\left({{\mathbf{a}}-{\mathbf{g}}}\right)\rVert.$

(15)

Further, considering aerodynamic forces and rotor thrust, the transitional dynamics of the aircraft expressed in ${}^{\rm w}\mathcal{F}$ are

$$m\cdot{{}^{\rm w}{\bf a}}={\mathbf{R}_{\rm b}^{\rm w}}\cdot{{}^{\rm b}{{\bf f}_{r}}}+{\mathbf{R}_{\rm l}^{\rm w}}\cdot{{}^{\rm l}{\bf f}_{a}}$$

(16)

where ${\mathbf{R}_{\rm b}^{\rm w}}={\mathbf{R}_{\rm l}^{\rm w}}{\mathbf{R}_{\rm b}^{\rm l}}$ and ${\mathbf{R}_{\rm l}^{\rm w}}$ is the transformation from ${}^{\rm l}{\mathcal{F}}$ to ${}^{\rm w}{\mathcal{F}}$ given by

$${\mathbf{R}_{\rm l}^{\rm w}}=\left[{\begin{array}[]{*{20}{c}}\cos\alpha&0&\sin\alpha\\ 0&1&0\\ -\sin\alpha&0&\cos\alpha\end{array}}\right].$$

(17)

The thrust $f_{z}$ and angle of attack $\alpha$ are now obtained by solving the equation (16) as

	$\displaystyle f_{z}$	$\displaystyle=-\frac{1}{k_{1}}\left|k_{f_{xs}}^{2}+{C_{L\alpha}}QSk_{f_{xw}}+{m^{2}}a_{z_{\rm w}}^{2}\right|,$		(18)
	$\displaystyle\alpha$	$\displaystyle=2{\rm atan}\left(\frac{{\left(C_{L\alpha}QS+k_{f_{xw}}\right)\sin\kappa-m{a_{z_{\rm w}}}\cos\kappa-{k_{1}}}}{{k_{f_{xw}}\cos\kappa+m{a_{z_{\rm w}}}\sin\kappa}}\right)$		(19)

where

	$\displaystyle k_{f_{xw}}$	$\displaystyle=C_{d0}QS+m{a_{x_{\rm w}}},$
	$\displaystyle{k_{1}}$	$\displaystyle=\sqrt{\Big{\{}\begin{aligned} {k_{f_{xs}}^{2}}{{\cos}^{2}}\kappa&+{{(k_{f_{xs}}+{C_{L\alpha}}QS)}^{2}}{{\sin}^{2}}\kappa+{m^{2}}a_{z_{\rm w}}^{2}\\ &-m{C_{L\alpha}}QS{a_{z_{\rm w}}}\sin 2\kappa\end{aligned}\Big{\}}}$

Substituting (4) and (6) into (2), the position dynamic equation is reformulated as

	$\displaystyle{c_{z}}{{\mathbf{z}}_{\rm b}}$	$\displaystyle-\left({{Q_{a}}{C_{{d_{x}}}}{\mathbf{x}}_{\rm b}^{\rm T}{{\mathbf{v}}}+{Q_{a}}{C_{{d_{xz}}}}{\mathbf{z}}_{\rm b}^{\rm T}{\mathbf{v}}}\right){{\mathbf{x}}_{\rm b}}-\left({{Q_{a}}{C_{{y_{0}}}}{\mathbf{y}}_{\rm b}^{\rm T}{\mathbf{v}}}\right){{\mathbf{y}}_{\rm b}}$
		$\displaystyle-\left({{Q_{a}}{C_{{d_{z}}}}{\mathbf{z}}_{\rm b}^{\rm T}{\mathbf{v}}+{Q_{a}}{C_{{d_{xz}}}}{\mathbf{x}}_{\rm b}^{\rm T}{\mathbf{v}}}\right){{\mathbf{z}}_{\rm b}}+g{{\mathbf{z}}_{\rm e}}-{\mathbf{a}}=0$		(20)

where ${C_{{d_{x}}}}={C_{{d_{0}}}}{\cos^{2}}\kappa+({C_{{L_{\alpha}}}}+{C_{{d_{0}}}}){\sin^{2}}\kappa$, ${C_{{d_{z}}}}={C_{{d_{0}}}}{\sin^{2}}\kappa+({C_{{L_{\alpha}}}}+{C_{{d_{0}}}}){\cos^{2}}\kappa$, ${C_{{d_{xz}}}}={C_{{L_{\alpha}}}}\sin\kappa\cos\kappa$, ${c_{z}}={{{f_{z}}}\mathord{\left/{\vphantom{{{f_{z}}}m}}\right.\kern-1.2pt}m}$, ${Q_{a}}={{\rho S{V_{a}}}\mathord{\left/{\vphantom{{\rho S{V_{a}}}{2m}}}\right.\kern-1.2pt}{2m}}$.

Left-multiplying (20) by ${\mathbf{y}}_{\rm b}^{\rm T}$ yields

$${\mathbf{y}}_{\rm b}^{\rm T}{{\mathbf{y}}_{\rm{perp}}}=0$$

(21)

with

$${{\mathbf{y}}_{\rm perp}}={Q_{a}}{C_{{y_{0}}}}{\mathbf{v}}-g{{\mathbf{z}}_{e}}+{\mathbf{a}}.$$

(22)

Under the coordinated turn condition, $\mathbf{x}_{\rm w}$ can be obtained by rotating $\mathbf{x}_{\rm b}$ along $\mathbf{y}_{\rm b}$ by $\alpha-\kappa$, which implies that $\mathbf{y}_{\rm b}$ is perpendicular to $\mathbf{x}_{\rm w}$. With this and (21), $\mathbf{y}_{\rm b}$ is obtained as

$${{\mathbf{y}}_{\rm b}}=\frac{{{\mathbf{x}}_{\rm w}}\times{{\mathbf{y}}_{\rm perp}}}{{\left\|{{{\mathbf{x}}_{\rm w}}\times{{\mathbf{y}}_{\rm perp}}}\right\|}}.$$

(23)

According to Rodrigues’ formula, $\mathbf{x}_{\rm b}$ is obtained as

	$\displaystyle{{\mathbf{x}}_{\rm b}}=\big{(}$	$\displaystyle\cos\left({\alpha-\kappa}\right){\mathbf{I}}+\left({1-\cos\left({\alpha-\kappa}\right)}\right){\mathbf{y}}_{\rm b}{\mathbf{y}}_{\rm b}^{\rm T}$
		$\displaystyle+\sin\left({\alpha-\kappa}\right){\left[{\mathbf{y}}_{\rm b}\right]_{\times}}\big{)}{{\mathbf{x}}_{\rm w}}.$		(24)

Then, utilizing the orthogonality of rotation matrix, $\mathbf{R}_{\rm b}^{\rm e}$ is constructed as

$$\mathbf{R}_{\rm b}^{\rm e}=\left[{\mathbf{x}}_{\rm b}\ {{\mathbf{y}}_{\rm b}}\ {{\mathbf{x}}_{\rm b}}\times{{\mathbf{y}}_{\rm b}}\right].$$

(25)

Taking the derivative of (20) yields

	$\displaystyle{\mathbf{j}}=$	$\displaystyle{\dot{c}_{z}}{{\mathbf{z}}_{\rm b}}+{c_{z}}{{\mathbf{\dot{z}}}_{\rm b}}-{Q_{a}}\mathbf{R}_{\rm b}^{\rm e}\left({{\left[{{}^{\text{b}}{\bm{\omega}}}\right]_{\times}}{\mathbf{D}}-{\mathbf{D}}{\left[{{}^{\text{b}}{\bm{\omega}}}\right]_{\times}}}\right){\left({\mathbf{R}_{\rm b}^{\rm e}}\right)^{\rm T}}{\mathbf{v}}$
		$\displaystyle-2{Q_{a}}\mathbf{R}_{\rm b}^{\rm e}{\mathbf{D}}{\left({\mathbf{R}_{\rm b}^{\rm e}}\right)^{\rm T}}{\mathbf{a}}$		(26)

where

$${\mathbf{D}}=\left[{\begin{array}[]{*{20}{c}}{{C_{{d_{x}}}}}&0&{{C_{{d_{xz}}}}}\\ 0&{{C_{{y_{0}}}}}&0\\ {{C_{{d_{xz}}}}}&0&{{C_{{d_{z}}}}}\end{array}}\right].$$

Left-multiplying (26) by ${\mathbf{x}}_{\rm b}^{\rm T}$ yields

$$\begin{gathered}{\mathbf{x}}_{\rm b}^{\rm T}{\mathbf{j}}+{Q_{a}}{{\mathbf{v}^{\rm T}\mathbf{a}}/V_{a}^{2}}\left({{d_{x}}{\mathbf{x}}_{\rm b}^{\rm T}{\mathbf{a}}+{d_{xz}}{\mathbf{z}}_{\rm b}^{\rm T}{\mathbf{a}}}\right)\hfill\\ ={\omega_{y_{\rm b}}}\left({{c_{z}}-{Q_{a}}\left({\left({{d_{z}}-{d_{x}}}\right){\mathbf{z}}_{\rm b}^{\rm T}{\mathbf{v}}+2{d_{xz}}{\mathbf{x}}_{\rm b}^{\rm T}{\mathbf{v}}}\right)}\right)\\ +{\omega_{x_{\rm b}}}{Q_{a}}{d_{xz}}{\mathbf{y}}_{\rm b}^{\rm T}{\mathbf{v}}-{\omega_{z_{\rm b}}}{Q_{a}}\left({{d_{x}}-{C_{{y_{0}}}}}\right){\mathbf{y}}_{\rm b}^{\rm T}{\mathbf{v}}.\end{gathered}$$

(30)

Left-multiplying (26) by ${\mathbf{y}}_{\rm b}^{\rm T}$ yields

	$\displaystyle{\mathbf{y}}_{\rm b}^{\rm T}{\mathbf{j}}+{Q_{a}}\left({{\mathbf{v}^{\rm T}\mathbf{a}}/V_{a}^{2}}\right){C_{{y_{0}}}}{\mathbf{y}}_{\rm b}^{\rm T}{\mathbf{a}}$
	$\displaystyle={\omega_{x_{\rm b}}}\left({-{c_{z}}+{Q_{a}}\left({{d_{xz}}{\mathbf{x}}_{\rm b}^{\rm T}{\mathbf{v}}+\left({{d_{z}}-{C_{{y_{0}}}}}\right){\mathbf{z}}_{\rm b}^{\rm T}{\mathbf{v}}}\right)}\right)$
	$\displaystyle-{\omega_{z_{\rm b}}}{Q_{a}}\left({\left({{d_{x}}-{C_{{y_{0}}}}}\right){\mathbf{x}}_{\rm b}^{\rm T}{\mathbf{v}}+{d_{xz}}{\mathbf{z}}_{\rm b}^{\rm T}{\mathbf{v}}}\right).$		(31)

Under the coordinated turn condition, there is also no lateral acceleration expressed in ${{}^{\rm b}\mathcal{F}}$ of the aircraft, that is

$${\mathbf{y}}_{\rm e}^{\rm T}\cdot{}^{\rm b}{\mathbf{\dot{v}}}=0.$$

(32)

With (32) and (2), another constraint on angular velocity is obtained

$${\omega_{x_{\rm b}}}{v_{{z_{\rm b}}}}-{\omega_{z_{\rm b}}}{v_{{x_{\rm b}}}}=-{g_{{y_{\rm b}}}}.$$

(33)

The angular velocity can now be obtained by solving the linear equations composed of (30), (31) and (33).

In this section, we will introduce some practical solutions to solve the special case that the proposed method is invalid.

Case 1 - $\lVert\mathbf{v}\rVert=0$: For the lifting-wing quadcopter, the case often occurs. For example, when the aircraft is commanded to hover, the reference velocity $\mathbf{v}_{\rm{ref}}$ is desired to be $\mathbf{0}$. In fact, when ${\mathbf{v}}\to{\mathbf{0}}$, ${{\mathbf{x}}_{\rm w}}$ and $\alpha$ are not defined, but this does not mean that ${{\mathbf{y}}_{\rm b}}$ and ${{\mathbf{x}}_{\rm b}}$ cannot be solved. In Section III-A, we use the wind speed vector to determine the nose orientation of the aircraft, but when the wind speed vector is 0, the nose orientation of the aircraft cannot be determined accordingly. In this case, $\mathbf{y}_{\rm b}$, originally determined by (23), can be any vector on the plane perpendicular to $\mathbf{y}_{\rm perp}$. To determine the attitude $\mathbf{R}_{\rm b}^{\rm e}$ of the aircraft, the yaw angle $\psi$ should be given when $\mathbf{v}_{\rm{ref}}=0$, then $\mathbf{y}_{\rm b}$ is obtained as

$${{\mathbf{y}}_{\rm b}}=\frac{{{\mathbf{x}}_{\rm c}}\times{{\mathbf{y}}_{\rm perp}}}{{\left\|{{{\mathbf{x}}_{\rm c}}\times{{\mathbf{y}}_{\rm perp}}}\right\|}}$$

(34)

where ${{\mathbf{x}}_{\rm{c}}}={[\begin{array}[]{*{20}{c}}{\cos\psi}&{\sin\psi}&0\end{array}]^{\rm{T}}}$.

Left-multiplying (20) by ${\mathbf{x}}_{\rm b}^{\rm T}$ yields

$${\mathbf{x}}_{\rm b}^{\rm T}\left({Q_{a}}{C_{{d_{x}}}}{\mathbf{v}}-g{{\mathbf{z}}_{\rm e}}+{\mathbf{a}}\right)={Q_{a}}{C_{{d_{xz}}}}{\mathbf{z}}_{\rm b}^{\rm T}{\mathbf{v}}.$$

(36)

That is, when ${\mathbf{v}}={\mathbf{0}}$, ${\mathbf{x}}_{\rm b}^{\rm T}{{\mathbf{y}}_{\rm perp}}={\mathbf{0}}$. With this and (21), $\mathbf{x}_{\rm b}$ is obtained as

$${{\mathbf{x}}_{\rm b}}=\frac{{{\mathbf{y}}_{\rm perp}}\times{{\mathbf{y}}_{\rm b}}}{{\left\|{{{\mathbf{y}}_{\rm perp}}\times{{\mathbf{y}}_{\rm b}}}\right\|}}.$$

(37)

However, state mutation may occur by directly giving an arbitrary yaw angle. Another simple and effective method is to keep $\mathbf{x}_{\rm w}$ to the last state when $\mathbf{v}$ is not equal to 0.

Case 2 - ${{\mathbf{x}}_{\rm w}}\times{{\mathbf{y}}_{\rm perp}}=0$: This case occurs if ${{\mathbf{y}}_{\rm perp}}$ is aligned with ${{\mathbf{x}}_{\rm w}}$, for example when the aircraft needs to execute the vertical takeoff or landing command. In this case, $\mathbf{y}_{\rm b}$ can be any vector on the plane perpendicular to $\mathbf{y}_{\rm perp}$. To overcome this singularity, $\mathbf{y}_{\rm b}$ can be obtained using the same method as in Case 1.

Cascaded PID controllers for position (1) and velocity (2) are designed as

$$\mathbf{v}_{\rm d}={\mathbf{K}_{\mathbf{p}p}}\left(\mathbf{p}_{\rm{ref}}-\mathbf{p}\right)+\mathbf{v}_{\rm{ref}},$$

(38)

$$\mathbf{a}_{\rm d}={\mathbf{K}_{\mathbf{v}p}}\left(\mathbf{v}_{\rm{d}}-\mathbf{v}\right)+{\mathbf{K}_{\mathbf{v}i}}\int\left(\mathbf{v}_{\rm{d}}-\mathbf{v}\right){\rm d}t+\mathbf{a}_{\rm{ref}}$$

(39)

where ${{\mathbf{K}}_{{\mathbf{p}}{\text{p}}}}$, ${{\mathbf{K}}_{{\mathbf{v}}{\text{p}}}}$, ${{\mathbf{K}}_{{\mathbf{v}}{i}}}\in{\mathbb{R}^{3\times 3}}$ are diagonal matrix acting as the feedback control gain. The integral term is added to compensate for model errors and external disturbances. According to the results of Section III, the desired collective thrust $f_{z,{\rm d}}$ and desired attitude $\mathbf{R}_{\rm d}$ can be obtained from $\bf a_{\rm d}$. However, in the actual flight, there is the inevitable deviation between the aerodynamic force calculated by (6) and acting on the aircraft. Multiplying the aerodynamic force by an attenuated feedforward coefficient can reduce the adverse effects caused by the inaccurate model. So the desired mass-normalized collective thrust expressed in ${}^{\rm e}{\mathcal{F}}$ is obtained as

$$\mathbf{a}_{\rm r,d}=\mathbf{a}_{\rm{d}}-{\mathbf{g}}-{\mathbf{K}_{\rm ff}}\mathbf{a}_{\rm{af}}$$

(40)

where ${\mathbf{K}_{\rm ff}}\in{\mathbb{R}^{3\times 3}}$ are diagonal matrix acting as the feedforward control gain, and $\mathbf{a}_{\rm{af}}$ is the acceleration due to aerodynamic forces, obtained as

$${{\mathbf{a}}_{{\rm{af}}}}=\frac{1}{{2m}}\rho S{V_{a}}{{\mathbf{R}}_{{\rm{ref}}}}\mathbf{D}{\mathbf{R}}_{{\rm{ref}}}^{\rm T}{\mathbf{v}_{\rm ref}}.$$

(41)

The desired mass-normalized collective thrust $\mathbf{a}_{\rm{r,d}}$ is completely produced by rotor thrust, so the desired attitude $\mathbf{R}_{\rm d}$ is obtained as

$$\begin{gathered}{{\mathbf{z}}_{\rm{b,d}}}=-\frac{{{{\mathbf{a}}_{\rm{r,d}}}}}{{\left\|{{{\mathbf{a}}_{\rm{d}}}}\right\|}}\hfill\\ {{\mathbf{y}}_{\rm{b,d}}}=\frac{{{{\mathbf{z}}_{\rm{b,d}}}\times{{\mathbf{x}}_{\rm{w}}}}}{{\left\|{{{\mathbf{z}}_{\rm{b,d}}}\times{{\mathbf{x}}_{\rm{w}}}}\right\|}}\hfill\\ {{\mathbf{x}}_{\rm{b,d}}}={{\mathbf{y}}_{\rm{b,d}}}\times{{\mathbf{z}}_{\rm{b,d}}}\hfill\\ \end{gathered}$$

(46)

where $\mathbf{x}_{\rm w}$ is obtained by (13). By projecting the desired acceleration onto the quadcopter-body $z$-axis, the collective thrust input is computed as

$$f_{z,{\rm d}}=m\left(\mathbf{z}_{\rm b,d}^{\rm T}\cdot\mathbf{a}_{\rm d}\right).$$

(47)

The attitude error is presented in the form of quaternion based on which the corresponding controller is designed:

$${\mathbf{q}_{\rm err}}={\mathbf{q}_{\rm d}^{*}}\otimes{\mathbf{q_{\rm b}^{\rm e}}}={\left[q_{\rm{e}_{0}}\ \ {q_{\rm{e}_{1}}}\ \ {q_{\rm{e}_{2}}}\ \ {q_{\rm{e}_{3}}}\right]^{\rm{T}}}$$

(48)

where ${\mathbf{q}}_{\rm{d}}$ is transformed from $\mathbf{R}_{\rm d}$ , $(\cdot)^{*}$ is the conjugate of a quaternion, and $\otimes$ is the quaternion product. Then, ${\mathbf{q}}_{\rm{err}}$ is transformed into the axis-angle form ${\mathbf{q}}_{\rm{err}}=\left[\cos\frac{\vartheta}{2}\ \ {\bm{\xi}}^{\rm T}_{\rm e}\sin\frac{\vartheta}{2}\right]^{\rm T}$ by

$$\begin{gathered}\vartheta={\rm{wrap}}_{\pi}\left({2{{\arccos}}\left(q_{\rm{e}_{0}}\right)}\right)\hfill\\ {\bm{\xi}}_{\rm{e}}={\rm{sign}}(q_{\rm{e}_{0}})\frac{\vartheta}{{\sin{\vartheta\mathord{\left/{\vphantom{\vartheta 2}}\right.\kern-1.2pt}2}}}{{\left[{q_{\rm{e}_{1}}\ \ q_{\rm{e}_{2}}\ \ q_{\rm{e}_{3}}}\right]}^{\rm{T}}}.\end{gathered}$$

(51)

The function ${\rm{wrap}_{\pi}}(\vartheta)$ constrains the $\vartheta$ in $\left[{-\pi\ \ \pi}\right]$ to ensure the shortest rotation path. To eliminate the attitude error, the attitude controller is designed as

$${}^{\rm{b}}{{\bm{\omega}}_{\rm{d}}}={{\mathbf{K}}_{{\mathbf{\Theta}}{\rm{p}}}}{{\bm{\xi}}_{\rm e}}+{{\bm{\omega}_{\rm{ref}}}}$$

(52)

where ${{\mathbf{K}}_{{\mathbf{\Theta}}{\rm{p}}}}\in{\mathbb{R}^{3\times 3}}$ is diagonal matrix acting as the control gain.

We prepared a lifting-wing quadcopter with a special installation angle for experiments, as shown in Fig. 5. The lifting wing is made of foam and connected to the quadcopter through a 3D-printing connector which controls the installation angle. The power system involves a 3300mAh, four-cell lithium polymer (LiPo) battery, EMAX RS2205 2600KV brushless DC motors with HQ5045 BN propeller, and 45A electronic speed control (ESC). A Pixhawk flight controller¹¹1https://docs.px4.io/master/en/assembly/quick_start_pixhawk4.html is rigidly attached to the center of the quadcopter frame, which has an on-board inertial measurement unit (IMU), barometer and SD card for flight data logging. A separate NEO-M8N based GPS module with a magnetometer is used to read the positioning and heading. A Benewake TFmini Lidar with a 40 m range is used for altitude estimation. A 900MHz telemetry link performs communication between the ground control station and the Pixhawk. The flight control algorithm runs on the Pixhawk at 250 Hz.

To validate our theoretical results, we present two trajectories, that is, circle and Lemniscate trajectories, and use the root mean square position error (RMSE) $E_{p}$ to evaluate the trajectory tracking performance, which is defined as

$${E_{p}}=\sqrt{\frac{1}{N}\sum\limits_{k=1}^{N}{\left({{{\bf{p}}_{{\rm{ref}}}}\left(k\right)-{\bf{p}}\left(k\right)}\right)}}$$

(53)

where $N$ is the control cycles of the position controller required to execute a given trajectory.

The circle trajectory is expressed as

$$\mathbf{p}_{\rm ref}=\mathbf{p}_{0}+\left[r(1-\cos(0.5\omega t^{2}))\ r\sin(0.5\omega t^{2})\ 0\right]^{\rm T}.$$

(54)

where $\mathbf{p}_{0}$ is the local position of the vehicle when the task starts. Here, we set $r=15$ m, $\omega=0.06$ rad/s, and when the magnitude of the reference velocity is accelerated to $10$ m/s, the system will switch to a circle trajectory with a constant speed.

We compare the position error of the lifting-wing quadcopter flying the circle trajectory described above for four conditions: (i) PID velocity controller expressed in (39) with feedforward references computed from the differential flatness transform considering aerodynamic forces (PID+DFAF), (ii) PD velocity controller with feedforward references considering aerodynamic forces (PD+DFAF), (iii) PID velocity controller with feedforward references without considering aerodynamic forces (PID+DF), and (iv) PD velocity controller with feedforward references without considering aerodynamic forces (PD+DF). We have also tried to remove the angular velocity feedforward. As a result, the aircraft deviated significantly from the trajectory when the speed is increased to about 9m/s. Therefore, this flight data is not included in the comparison.

Fig. 6 shows experimental results for tracking a circular trajectory reference on conditions (ii) and (iv). The velocity has a low-frequency sine wave, caused by the state estimator, because the update frequency of GPS is only 10Hz, which brings a considerable delay to the position estimation. The altitude tracking error is far less than that of the horizontal position, because, on the one hand, the Lidar provides higher measurement accuracy than GPS, and on the other hand, the altitude command is constant. The tracking performance statistics for four conditions are summarized in Fig. 7. From these statistics, it can be observed that the trajectory tracking performance has been improved significantly when considering aerodynamic forces and the integral can reduce tracking error, but its effect is very limited. The $E_{p}$ of conditions (i) and (ii) are about half of that of (iii) and (iv). Fig. 8 shows states of the vehicle agree well with the feedforward commands, demonstrating that the simplified aerodynamic model (6) has fairly good fitting accuracy.

The Lemniscate trajectory is expressed as

$$\mathbf{p}_{\rm ref}=\mathbf{p}_{0}+\left[r(1-\cos(0.5\omega t))\ r\sin(\omega t)\ 0\right]^{\rm T}.$$

(55)

Here, we set $r=20$ m, $\omega=0.33$ rad/s. Fig. 9 shows experimental results for tracking a Lemniscate of Gerono with an increasing speed from 0 m/s to 8 m/s on conditions (ii) and (iv) with the same aerodynamic coefficients in the circle trajectory. It also can be observed that the trajectory tracking performance has been improved significantly when considering aerodynamic forces, with the $E_{p}$ reduced by half.