2.1 Preliminary definition

Let ${C}_{\textrm{W}}:\left\{{O}_{\textrm{W}};\ {X}_{\textrm{W}},{Y}_{\textrm{W}},{Z}_{\textrm{W}}\right\}$ denote the world inertial frame and ${C}_{\textrm{B}}:\left\{{O}_{\textrm{B}};\ {X}_{\textrm{B}},{Y}_{\textrm{B}},{Z}_{\textrm{B}}\right\}$ denote the moving frame attached to the quadrotor’s body at its centre of mass. We also define ${C}_{{\textrm{P}}_{i}}:\left\{{O}_{{\textrm{P}}_{i}};\ {X}_{{\textrm{P}}_{i}},{Y}_{{\textrm{P}}_{i}},{Z}_{{\textrm{P}}_{i}}\right\}$ , $i=1,\cdots ,4$ , as the frames associated with the $i$ th propeller group, with ${X}_{{\textrm{P}}_{i}}$ denote the tilting actuation axis and ${Z}_{{\textrm{P}}_{i}}$ the propeller actuated spinning axis.

We define ${\mathbf{R}}_{12}\in $ SO(3) as the rotation matrix representing the orientation of the frame $1$ with respect to the frame $2$ . Therefore, the matrix ${\mathbf{R}}_{\textrm{B}\textrm{W}}$ represents the orientation of the body frame w.r.t. the world frame, while the matrix ${\mathbf{R}}_{{\textrm{P}}_{i}\textrm{B}}$ represents the orientation of the propeller group $i$ th frame w.r.t. the body frame. By denoting ${\alpha }_{i}$ as the propeller tilt angle about the axis ${X}_{{\textrm{P}}_{i}}$ , it follows from Fig. 1 that

(1) \begin{align}{\mathbf{R}}_{{\textrm{P}}_{i}\textrm{B}}={\mathbf{R}}_{Z}\left((i-1)\frac{\pi }{2}\right){\mathbf{R}}_{X}\left({\alpha }_{i}\right).\end{align}

Among them, ${\mathbf{R}}_{X}=\left[\begin{array}{lll}1 {}&\quad 0 &\quad {} 0 \\ 0 {}&\quad \textrm{c}\textrm{o}\textrm{s}\left(\varphi \right) {}&\quad -\textrm{s}\textrm{i}\textrm{n}\left(\varphi \right)\\ 0 {}&\quad \textrm{s}\textrm{i}\textrm{n}\left(\varphi \right) {}&\quad \textrm{c}\textrm{o}\textrm{s}\left(\varphi \right)\end{array}\right]$ and ${\mathbf{R}}_{Z}=\left[\begin{array}{lll}\textrm{c}\textrm{o}\textrm{s}\left(\theta \right)&\quad {} 0 {}&\quad \textrm{s}\textrm{i}\textrm{n}\left(\theta \right) \\ 0&\quad {} 1&\quad {} 0 \\ -\textrm{s}\textrm{i}\textrm{n}\left(\theta \right) {}&\quad 0&\quad {} \textrm{c}\textrm{o}\textrm{s}\left(\theta \right)\end{array}\right]$ . They represent the rotation matrices around the $X$ -axis and $Z$ -axis, respectively. Similarly, we let

(2) \begin{align}{\mathbf{o}}_{{\textrm{P}}_{i}\textrm{B}}={\mathbf{R}}_{Z}\left((i-1)\frac{\pi }{2}\right){\left[h \quad 0 \quad 0\right]}^{{\top }}\end{align}

be the origin of the frames associated with the propeller’s ${C}_{{\textrm{P}}_{i}}$ in the body frame, with $h$ representing the distance between ${\mathbf{O}}_{{\textrm{P}}_{i}}$ from ${\mathbf{O}}_{\textrm{B}}$ .

Let $\mathbf{p}={\mathbf{o}}_{\textrm{B}\textrm{W}}\in {\mathbb{R}}^{3\times 1}$ denote the position of the body in the world frame, and let ${\mathbf{R}}_{\textrm{B}\textrm{W}}$ denote its orientation. Moreover, the propeller tilt angles ${\alpha }_{i}$ (rotations about ${X}_{{\textrm{P}}_{i}}$ ) specify the orientations of the propeller groups relative to the body frame. We omit the propeller spinning angles about ${Z}_{{\textrm{P}}_{i}}$ from the set of configuration variables, although the propeller spinning velocities ${\overline{\omega }}_{i}$ about ${Z}_{{\textrm{P}}_{i}}$ are incorporated into the system model. The detailed part will be introduced later.

2.2 Equations of dynamics

In the last subsection, some basic definitions are introduced. We focus on introducing dynamic equations in this subsection.

(3) \begin{align}\left\{\begin{array}{l}\ddot{\mathbf{p}}=-{\mathbf{e}}_{3}g+{m}^{-1}{\mathbf{R}}_{\mathrm{B}\mathrm{W}}\sum _{i=1}^{4} {\mathbf{R}}_{{\mathrm{P}}_{i}\mathrm{B}}{\mathbf{t}}_{{\mathrm{P}}_{i}}+{m}^{-1}{\mathbf{f}}_{1}\\ \\ {\dot{\mathit{\omega }}}_{\mathrm{B}}={\mathbf{I}}_{\mathrm{B}}^{-1}\sum _{i=1}^{4} \left({\mathbf{o}}_{{\mathrm{P}}_{i}\mathrm{B}}\times {\mathbf{R}}_{{\mathrm{P}}_{i}\mathrm{B}}{\mathbf{t}}_{{\mathrm{P}}_{i}}\right)+{\mathbf{I}}_{\textrm{B}}^{-1}{\mathbf{f}}_{2}\end{array}\right.\end{align}

Equation (3) represents the dynamic model of the vector quadrotor, where $\mathbf{p}={\left[x\textrm{}\textrm{}\textrm{}\textrm{}y\textrm{}\textrm{}\textrm{}\textrm{}z\right]}^{{\top }}$ and ${\mathit{\omega }}_{\textrm{B}}={\left[\varphi \textrm{}\textrm{}\textrm{}\textrm{}\theta \textrm{}\textrm{}\textrm{}\textrm{}\psi \right]}^{{\top }}$ denote the position and attitude in an inertial coordinate system. The total mass of the vector quadrotor, is $m$ , the gravitational acceleration is $g$ and the vector ${\mathbf{e}}_{3}=[0\textrm{}\textrm{}\textrm{}\textrm{}0\textrm{}\textrm{}\textrm{}\textrm{}1{]}^{{\top }}$ . The rotation matrix ${\mathbf{R}}_{\textrm{B}\textrm{W}}$ represents the body coordinate system to the world coordinate system. The rotation matrix ${\mathbf{R}}_{{\textrm{P}}_{i}\textrm{B}}$ represents the rotor coordinate system to the body coordinate system. The thrust ${\mathbf{t}}_{{\textrm{P}}_{i}}$ is generated by a single motor, and ${\mathbf{I}}_{\textrm{B}}$ is the moment of inertia.

2.2.1 Analysis of dynamic equations

By applying standard techniques, a comprehensive description of the quadrotor’s dynamic model can be derived by accounting for the forces and torques produced by the propeller’s motion, along with any cross-coupling effects resulting from gyroscopic and inertial interactions arising from the relative motion of the five bodies that compose the quadrotor. Regarding aerodynamics and torque, we only consider the factors that drive the quadrotor, and the sum of the known and unknown disturbances in the system represented by ${\mathbf{f}}_{1}$ and ${\mathbf{f}}_{2}$ . We now outline the key conceptual steps required to derive the quadrotor’s dynamic model. To this end, let ${\mathit{\omega }}_{\textrm{B}}\in {\mathbb{R}}^{3\times 1}$ denote the angular velocity of the quadrotor body $\textrm{B}$ expressed in its body frame $\beta $ , and consider the $i$ th propeller group ${\textrm{P}}_{i}$ . The angular velocity of the propeller, expressed in its frame, is defined as

(4) \begin{align}{\mathit{\omega }}_{{\textrm{P}}_{i}}={\mathbf{R}}_{{\textrm{P}}_{i}\textrm{B}}^{{\top }}{\mathit{\omega }}_{\textrm{B}}+{\left[{\dot{\alpha }}_{i} \quad 0 \quad {\overline{\omega }}_{i}\right]}^{{\top }}\end{align}

where ${\overline{\omega }}_{i}$ is the spinning velocity of the propeller about ${Z}_{{\textrm{P}}_{i}}$ . This produces an angular acceleration

(5) \begin{align}{\dot{\mathit{\omega }}}_{{\textrm{P}}_{i}}={\mathbf{R}}_{{\textrm{P}}_{i}\textrm{B}}^{{\textrm{}}^{{\top }}}{\dot{\mathit{\omega }}}_{\textrm{B}}+{\dot{\mathbf{R}}}_{{\textrm{P}}_{i}\textrm{B}}^{{\top }}{\mathit{\omega }}_{\textrm{B}}+{\left[{\ddot{\alpha }}_{i} \quad 0 \quad {\overline{\omega }}_{i}\right]}^{{\top }}.\end{align}

From the application of Euler’s equations of motion, it follows that

(6) \begin{align}{\mathit{\tau }}_{{\textrm{P}}_{i}}={\mathbf{I}}_{{\textrm{P}}_{i}}{\dot{\mathit{\omega }}}_{{\textrm{P}}_{i}}+{\mathit{\omega }}_{{\textrm{P}}_{i}}\times {\mathbf{I}}_{{\textrm{P}}_{i}}{\mathit{\omega }}_{{\textrm{P}}_{i}},\end{align}

where the torque ${\mathit{\tau }}_{{\textrm{P}}_{i}}$ represents the torque generated when a single rotor rotates, and matrix ${\mathbf{I}}_{{\textrm{P}}_{i}}\in {\mathbb{R}}^{3\times 3}$ is the constant symmetric and positive definite inertia matrix of the propeller group. We assume the presence of a counter-rotating torque about the ${Z}_{{\textrm{P}}_{i}}$ axis, caused by air drag and modelled as

(7) \begin{align}{\mathit{\tau }}_{\textrm{e}\textrm{x}{\textrm{t}}_{i}}={\left[0 \quad 0 \quad -{k}_{m}{\mathit{\omega }}_{{\textrm{P}}_{{i}_{Z}}}\left|{\mathit{\omega }}_{{\textrm{P}}_{{i}_{Z}}}\right|\right]}^{{\top }}\end{align}

where ${\mathit{\tau }}_{\textrm{e}\textrm{x}{\textrm{t}}_{i}}$ denotes any external torque applied to the propeller and ${k}_{m}$ is the drag coefficient. Considering the quadrotor body $\textrm{B}$ and its four propellers ${\textrm{P}}_{i}$ , we obtain

(8) \begin{align}{\mathit{\tau }}_{\textrm{B}}={\mathbf{I}}_{\textrm{B}}{\dot{\mathit{\omega }}}_{\textrm{B}}+{\mathit{\omega }}_{\textrm{B}}\times {\mathbf{I}}_{\textrm{B}}{\mathit{\omega }}_{\textrm{B}}+\sum _{i=1}^{4} {\mathbf{R}}_{{\textrm{P}}_{i}\textrm{B}}{\mathit{\tau }}_{{\textrm{P}}_{i}},\end{align}

where ${\mathbf{I}}_{\textrm{B}}\in {\mathbb{R}}^{3\times 3}$ is the constant, symmetric and positive definite inertia matrix of $\textrm{B}$ . The torques ${\mathit{\tau }}_{\textrm{B}}\in {\mathbb{R}}^{3\times 1}$ acting on $\textrm{B}$ are generated by the moments produced by the four propeller forces along the ${Z}_{{\textrm{P}}_{i}}$ axes,

(9) \begin{align}{\mathbf{t}}_{{\textrm{P}}_{i}}={\left[0 \quad 0 \quad {k}_{f}{\overline{\omega }}_{i}\left|{\overline{\omega }}_{i}\right|\right]}^{{\top }}\end{align}

acting at ${\mathbf{o}}_{{\textrm{P}}_{i}\textrm{B}}$ in ${C}_{\textrm{B}}$ (the propeller drag torques are already accounted for in the terms ${\mathit{\tau }}_{{\textrm{P}}_{i}}$ ), and ${k}_{f}$ is the thrust coefficient. According to (9), we obtain

(10) \begin{align}{\mathit{\tau }}_{\textrm{B}}=\sum _{i=1}^{4} \left({\mathbf{o}}_{{\textrm{P}}_{i}\textrm{B}}\times {\mathbf{R}}_{{\textrm{P}}_{i}\textrm{B}}{\mathbf{t}}_{{\textrm{P}}_{i}}\right).\end{align}

Regarding the translational dynamics, we assume for simplicity that the barycentre of each propeller group ${\textrm{P}}_{i}$ coincides with ${\mathbf{O}}_{{\textrm{P}}_{i}}$ , allowing us to neglect the inertial effects on the propeller groups caused by the quadrotor body’s acceleration in space. Therefore, recalling that $\mathbf{p}$ represents the quadrotor body’s position in the world frame, we obtain

(11) \begin{align}m\ddot{\mathbf{p}}=-mg{\mathbf{e}}_{3}+{\mathbf{R}}_{\textrm{B}\textrm{W}}\sum _{i=1}^{4} {\mathbf{R}}_{{\textrm{P}}_{i}\textrm{B}}{\mathbf{t}}_{{\textrm{P}}_{i}}.\end{align}

In summary, Equations (4)-(11) characterise the rotational and translational dynamics of the quadrotor body and propeller groups. Note that the inputs of the model are the rotor torques that actuate the propeller tilting axes ${X}_{{\textrm{P}}_{i}}$ and spinning axes ${Z}_{{\textrm{P}}_{i}}$ . The propeller spinning velocities ${\overline{\omega }}_{i},$ which generate the forces and torques (6)-(8), which influence the translational and rotational motion of the quadrotor body $\textrm{B}$ as a function of its current configuration, particularly the tilting angles ${\alpha }_{i}$ which are actuated by ${\mathit{\tau }}_{{\alpha }_{i}}$ .

2.2.2 Summarising of dynamic equations

We now define servo angle $\mathit{\alpha }={\left[{\alpha }_{1},\cdots ,{\alpha }_{4}\right]}^{{\top }}\in {\mathbb{R}}^{4\times 1}$ and rotor speed $\mathit{\omega }={\left[{\overline{\omega }}_{1}\left|{\omega }_{1}\right|,\cdots ,{\overline{\omega }}_{4}\left|{\omega }_{4}\right|\right]}^{{\top }}\in {\mathbb{R}}^{4\times 1}$ . Note that the elements of the vector $\mathit{\omega }$ are the signed squares of the spinning velocities ${\overline{\omega }}_{i}$ , since the torques and forces in (8)–(11) are approximately linear functions of these quantities. Therefore, in the following analysis, ${\omega }_{i}={\overline{\omega }}_{i}\left|{\overline{\omega }}_{i}\right|$ will be considered as the spinning velocity of the $i$ th propeller, with the understanding that the actual speed can always be recovered as $\overline{\omega }=\textrm{s}\textrm{i}\textrm{g}\textrm{n}\left({\omega }_{i}\right)\sqrt{\left|{\omega }_{i}\right|}$ , where $\textrm{s}\textrm{i}\textrm{g}\textrm{n}\left({\omega }_{i}\right)$ is signum function. Under these assumptions and neglecting external disturbances, the vector quadrotor dynamic model simplifies to:

(12) \begin{align}2pt\left\{\begin{array}{ll}\ddot{\mathbf{p}} &= -{\mathbf{e}}_{3}g+{m}^{-1}\textrm{}{\mathbf{R}}_{\textrm{B}\textrm{W}}\mathbf{F}\left(\alpha \right)\mathit{\omega }\\[3pt] {\dot{\mathit{\omega }}}_{\textrm{B}} &={\mathbf{I}}_{\textrm{B}}^{-1}\mathit{\tau }\left(\alpha \right)\mathit{\omega }\\[3pt] {\dot{\mathbf{R}}}_{\textrm{B}\textrm{W}} &={\mathbf{R}}_{\textrm{B}\textrm{W}}\textrm{}[{\mathit{\omega }}_{\textrm{B}}{]}_{\wedge }\end{array}\right.,\end{align}

where $[\cdot {]}_{\wedge }$ denotes the standard mapping from a vector $\alpha \in {\mathbb{R}}^{4\times 1}$ to its corresponding skew-symmetric matrix $\left[\alpha \right]\in \textrm{S}\textrm{O}\left(3\right)$ , and

(13) \begin{align}\mathbf{F}\left(\alpha \right)=\left[\begin{array}{cccc} 0 {}&\quad -{k}_{f}{\textrm{s}}_{2} {}&\quad 0 {}&\quad {k}_{f}{\textrm{s}}_{4} \\ -{k}_{f}{\textrm{s}}_{1} {}&\quad 0 &\quad {} {k}_{f}{\textrm{s}}_{3} &\quad 0 \\ {k}_{f}{\textrm{c}}_{1} &\quad -{k}_{f}{\textrm{c}}_{2} &\quad {} {k}_{f}{\textrm{c}}_{3} &\quad {} -{k}_{f}{\textrm{c}}_{4}\\ \end{array}\right],\end{align}

(14) \begin{align}\begin{array}{l}\mathit{\tau }\left(\alpha \right)=\left[\begin{array}{cccc}0 &\quad {} h{k}_{f}{\textrm{c}}_{2}-{k}_{m}{\textrm{s}}_{2} &\quad {} 0 {}&\quad h{k}_{f}{\textrm{c}}_{4}+{k}_{m}{\textrm{s}}_{4}\\ -h{k}_{f}{\textrm{c}}_{1}+{k}_{m}{\textrm{s}}_{1} {} &\quad 0 {}&\quad h{k}_{f}{\textrm{c}}_{3}-{k}_{m}{\textrm{s}}_{3} &\quad {} 0\\ -h{k}_{f}{\textrm{s}}_{1}-{k}_{m}{\textrm{c}}_{1} {} &\quad h{k}_{f}{\textrm{s}}_{2}-{k}_{m}{\textrm{c}}_{2} {}&\quad -h{k}_{f}{\textrm{s}}_{3}-{k}_{m}{\textrm{c}}_{3} &\quad {} h{k}_{f}{\textrm{s}}_{4}-{k}_{m}{\textrm{c}}_{4}\\ {} {} {} \end{array}\right] \end{array}, \\[3pt] \nonumber \end{align}

matrices $\mathbf{F}\!\left(\alpha \right)\in {\mathbb{R}}^{3\times 4}$ , $\mathit{\tau }\left(\alpha \right)\in {\mathbb{R}}^{3\times 4}$ representing the allocation matrix of vector quadrotor, where ${\textrm{s}}_{i}=\textrm{s}\textrm{i}\textrm{n}\left({\alpha }_{i}\right)$ and ${\textrm{c}}_{i}=\textrm{c}\textrm{o}\textrm{s}\left({\alpha }_{i}\right)$ . It is noteworthy that the input $\mathit{\omega }$ appears linearly in Equation (12), as anticipated. The subsequent controller design is conducted based on the simplified model given by (12)–(14).

2.3 Problem formulation

This paper aims to design a controller for vector quadrotors that regulates position and attitude, enabling the quadrotor to track a desired trajectory while compensating for external disturbances and nonlinear system dynamics. Specifically, the control problem is defined as follows: using the available control inputs, track an arbitrary desired trajectory characterised by ${\mathbf{p}}_{\textrm{d}}\left(t\right)\in {\mathbb{R}}^{3\times 1}$ (desired position) and ${\mathbf{r}}_{\textrm{d}}\left(t\right)\in {\mathbb{R}}^{3\times 1}$ (desired attitude). Here, $\mathbf{p}\left(t\right)\in {\mathbb{R}}^{3\times 1}$ (actual position) and $\mathbf{r}\left(t\right)\in {\mathbb{R}}^{3\times 1}$ (actual attitude) denote the position and attitude of the vector quadrotor’s main body, which serve as the system’s output functions. It should be emphasised that this trajectory-tracking problem is inherently ill-posed for standard quadrotors with fixed propellers, due to their underactuated nature. In contrast, the overactuated design of the vector quadrotor–featuring eight independent control inputs–guarantees complete controllability of both the position and orientation of its main body. This controllability will be demonstrated in subsequent analyses.

The primary objective of the controller design is high-precision trajectory tracking, which enables the vector quadrotor’s position and attitude to converge rapidly and robustly to their commanded reference values. The controller should minimise rise time, overshoot and steady-state error while maintaining stability under model uncertainties, external disturbances and sensor noise. To achieve this goal, a tightly coordinated inner-loop attitude regulator and outer-loop position tracker are required. These components should be supplemented with adaptive or robust mechanisms that sustain performance throughout the vector quadrotor’s full flight envelope.

In summary, requirements for trajectory tracking are as follows

(15) \begin{align}\left\{\begin{array}{lll}\underset{t\to {\infty }}{\textrm{l}\textrm{i}\textrm{m}}{\mathbf{p}}_{\textrm{d}}\left(t\right)-\mathbf{p}\left(t\right)=0 {} {} \\[5pt] \underset{t\to {\infty }}{\textrm{l}\textrm{i}\textrm{m}}{\mathbf{r}}_{\textrm{d}}\left(t\right)-\mathbf{r}\left(t\right)=0 {} {} \\ {} {} \end{array}\right..\end{align}

Subsequently, corresponding controllers for attitude and position is designed, along with methods for calculating the spinning velocities ${\omega }_{i}$ and tilt angular ${\alpha }_{i}$ , in the next section.

3.1 Model compensating control

In this subsection, we will introduce model compensation control, which includes high-order differentiator and a compensation function observer.

A second-order differentiable single-input single-output (SISO) controlled system is defined as

(16) \begin{align}\ddot{y}=bu+{f}_{k}+{f}_{u},\end{align}

where $y$ denotes the system output, and $b$ represents the control input coefficient. $u$ denotes the system input, ${f}_{k}$ represents the known model information, and ${f}_{u}$ accounts for the unknown disturbances in the model.

In highly manoeuverable flights, disturbances mainly originate from strong winds, gyroscopic torques, aerodynamic forces and complex mechanical structures, all of which are bounded and continuous. The block diagram illustrating the vector quadrotor’s control structure is shown in Fig. 2. HOD estimates the derivatives of the desired signal. CFO estimates the unknown disturbances, providing replacements for the unknown perturbations ${\widehat{\mathbf{f}}}_{1}$ , ${\widehat{\mathbf{f}}}_{2}$ in (3). CFO and HOD will be discussed in the following.

Figure 2.

Block diagram of vector quadrotor control that SO(3)-based MCC with SQP.

For (16), the control rate is

(17) \begin{align}u=\frac{1}{b}(\!\!-{k}_{p}(y-{y}_{\textrm{d}})-{k}_{d}(\dot{y}-{\widehat{\dot{y}}}_{\textrm{d}})+{\widehat{\ddot{y}}}_{\textrm{d}}-{\widehat{f}}_{u}-{\widehat{f}}_{k}),\end{align}

where ${k}_{\textrm{p}}$ , ${k}_{\textrm{d}}$ represents the error gain factors, ${\widehat{\dot{y}}}_{\textrm{d}}$ , ${\widehat{\ddot{y}}}_{\textrm{d}}$ are estimated by HOD, and ${\widehat{f}}_{u}$ is observed by CFO.

MCC fully utilises the known model information ${\widehat{f}}_{k}$ and compensates for the unknown disturbance ${\widehat{f}}_{u}$ in real time, thereby enhancing its anti-disturbance performance. Additionally, the differential information in MCC provides a foundation for high-manoeuverability control in servo tracking. Therefore, MCC compensates for system nonlinearities, enabling the closed-loop system to achieve predefined dynamics similar to those of linear systems.

Certain issues arise when applying it to attitude control. Specifically, the formulation in (17) is only applicable within Euler’s framework. This approach decomposes attitude into pitch, roll and yaw, simplifying attitude control. However, it disrupts the continuity of the control process. In the $\textrm{S}\textrm{O}\left(3\right)$ space, the attitude control trajectories lose their smoothness. When the vector quadrotor undergoes large-angle flips, the decoupled control significantly restricts its manoeuverability. In extreme cases, the control system may encounter singularities. Control formulated in $\textrm{S}\textrm{O}\left(3\right)$ is free from this issue. It is grounded in mathematical theories, including attitude dynamics, directional representation and rotation vectors. It’s goal is ${\textrm{l}\textrm{i}\textrm{m}}_{t\to {\infty }}\mathbf{r}\left(t\right)={\mathbf{r}}_{\textrm{d}}$ . Essentially, it represents a continuous geometric control approach. It avoids the issue of universal joint locking present in Euler angles and the ambiguity inherent in quaternion representations. The control approach exhibits near-global stability. It maintains good performance even under drastic attitude changes. Following, we extends MCC framework to the $\textrm{S}\textrm{O}\left(3\right)$ space.

3.2 Attitude errors and controller

Assume that an arbitrary, smooth, desired attitude trajectory ${\mathbf{r}}_{\textrm{d}}\left(t\right)\in \textrm{S}\textrm{O}\left(3\right)$ is provided. The corresponding angular velocity of command can be derived using the attitude dynamic equation, ${\widehat{\mathit{\omega }}}_{\textrm{B}}={\mathbf{r}}_{\textrm{d}}^{{\top }}{\dot{\mathbf{r}}}_{\textrm{d}}$ . We define the error terms related to the attitude dynamics of the quadrotor. The attitude error and angular velocity error function investigated in (18) and some of its key properties are summarised as follows. For a given command ${\mathbf{r}}_{\textrm{d}},{\mathit{\omega }}_{\textrm{d}}$ , and the current attitude and angular velocity of the vector quadrotor, denoted by $\mathbf{r},{\mathit{\omega }}_{\textrm{B}}$ . We define the attitude error vector as ${\mathbf{e}}_{\mathbf{r}}\in {\mathbb{R}}^{3\times 1}$ , and we define the angular velocity error vector ${\mathbf{e}}_{\mathit{\omega }}\in {\mathbb{R}}^{3\times 1}$ as follows

(18) \begin{align}\left\{\begin{array}{l}{\mathbf{e}}_{\mathbf{r}}=\frac{1}{2}{\left({\mathbf{r}}_{\textrm{d}}^{{\top }}-{\mathbf{r}}^{{\top }}{\mathbf{r}}_{\textrm{d}}\right)}^{\vee }\\[3pt] {\mathbf{e}}_{\mathit{\omega }}={\mathit{\omega }}_{\textrm{B}}-{\mathbf{r}}^{{\top }}{\mathbf{r}}_{\textrm{d}}{\mathit{\omega }}_{\textrm{d}}\end{array}\right..\end{align}

The authors in Ref. (Reference Lee, Leok and McClamroch35) introduce a nonlinear controller for the attitude-controlled flight mode, defined by the moment vector expression

(19) \begin{align}\begin{array}{ll}{\mathit{\tau }}_{\textrm{d}}= {} -{k}_{\mathbf{r}}{\mathbf{e}}_{\mathbf{r}}-{k}_{\mathit{\omega }}{\mathbf{e}}_{\mathit{\omega }}+{\mathit{\omega }}_{\textrm{B}}\times {\mathbf{I}}_{\textrm{B}}{\mathit{\omega }}_{\textrm{B}}-{\mathbf{I}}_{\textrm{B}}\left({\widehat{\mathit{\omega }}}_{\textrm{B}}{\mathbf{r}}^{{\top }}{\mathbf{r}}_{\textrm{d}}{\mathit{\omega }}_{\textrm{d}}-{\mathbf{r}}^{{\top }}{\mathbf{r}}_{\textrm{d}}{\dot{\mathit{\omega }}}_{\textrm{d}}\right)\end{array},\end{align}

where ${k}_{\mathbf{r}}$ , ${k}_{\mathit{\omega }}\in {\mathbb{R}}^{3\times 1}$ represent the error gain factors, ${\mathit{\omega }}_{\textrm{d}}$ and ${\dot{\mathit{\omega }}}_{\textrm{d}}$ are computed using HOD. Afterwards, we introduced CFO to calculate external disturbances ${\widehat{\mathbf{f}}}_{1}$ and obtain a new control rate

(20) \begin{align}\begin{array}{ll}{\mathit{\tau }}_{\textrm{d}}= {} -{k}_{\mathbf{r}}{\mathbf{e}}_{\mathbf{r}}-{k}_{\mathit{\omega }}{\mathbf{e}}_{\mathit{\omega }}+{\mathit{\omega }}_{\textrm{B}}\times {\mathbf{I}}_{\textrm{B}}{\mathit{\omega }}_{\textrm{B}}-{\mathbf{I}}_{\textrm{B}}\left({\widehat{\mathit{\omega }}}_{\textrm{B}}{\mathbf{r}}^{{\top }}{\mathbf{r}}_{\textrm{d}}{\mathit{\omega }}_{\textrm{d}}-{\mathbf{r}}^{{\top }}{\mathbf{r}}_{\textrm{d}}{\dot{\mathit{\omega }}}_{\textrm{d}}-{\widehat{\mathbf{f}}}_{1}\right)\end{array}.\end{align}

3.3 Position errors and controller

We introduce a nonlinear controller for the position-controlled flight mode. This flight mode necessitates analysing the coupled translational and rotational equations of motion. Therefore, we utilise the notation and analysis from the previous subsection to characterise the properties of the closed-loop system in this mode.

We define the position error ${\mathbf{e}}_{\mathbf{p}}$ and velocity error ${\mathbf{e}}_{\dot{\mathbf{p}}}$ . An arbitrary smooth position tracking reference ${\mathbf{p}}_{\textrm{d}}\left(t\right)\in {\mathbb{R}}^{3\times 1}$ is selected. The position and velocity tracking errors are defined as follows

(21) \begin{align}\left\{\begin{array}{l}{\mathbf{e}}_{\mathbf{p}}=\mathbf{p}-{\mathbf{p}}_{\textrm{d}}\\[3pt] {\mathbf{e}}_{\mathbf{v}}=\mathbf{v}-{\dot{\mathbf{p}}}_{\textrm{d}}\end{array}\right..\end{align}

Next, we introduce the position controller, which is defined by control expressions for the thrust magnitude vector, as follows

(22) \begin{align}\mathbf{f}=\left({k}_{\mathbf{p}}{\mathbf{e}}_{\mathbf{p}}+{k}_{\dot{\mathbf{p}}}{\mathbf{e}}_{\dot{\mathbf{p}}}+mg{\mathbf{e}}_{3}-m{\ddot{\mathbf{p}}}_{\textrm{d}}-{\widehat{\mathbf{f}}}_{2}\right){\mathbf{R}}_{\textrm{B}\textrm{W}}{\mathbf{e}}_{3},\end{align}

where ${k}_{\mathbf{p}}$ , ${k}_{\dot{\mathbf{p}}}$ are the gain factors, ${\dot{\mathbf{p}}}_{\textrm{d}}$ , ${\ddot{\mathbf{p}}}_{\textrm{d}}$ are estimated by HOD and introduced CFO to calculate external disturbances ${\widehat{\mathbf{f}}}_{u2}$ .

3.4 Compensation function observer

The term ${f}_{u}$ in (16) is crucial for the vector quadrotor’s anti-disturbance capability. To estimate $n$ th order CFO is designed in Ref. (Reference Qi, Hu, Li and Li25). The state equation is

(23) \begin{align}\left\{\begin{array}{l}{\dot{z}}_{\textrm{c}1}={z}_{\textrm{c}2}\\[3pt] {\dot{z}}_{\textrm{c}2}={z}_{\textrm{c}3}\\[3pt] \textrm{}\textrm{}\textrm{}\textrm{} \vdots \\[3pt] {\dot{z}}_{\textrm{c}n}=\mathbf{l}{\mathbf{e}}_{\textrm{c}}+{z}_{\textrm{c}(n+1)}+bu\\[3pt] {\dot{z}}_{\textrm{c}(n+1)}=\lambda \mathbf{l}{\mathbf{e}}_{\textrm{c}}\\[3pt] \widehat{f}=\mathbf{l}{\mathbf{e}}_{\textrm{c}}+{z}_{\textrm{c}(n+1)}\end{array}\right.,\end{align}

where ${z}_{\textrm{c}1},\cdots ,{z}_{\textrm{c}n},{z}_{\textrm{c}(n+1)}$ represent the internal state variables of the observer, ${\mathbf{e}}_{\textrm{c}}={\left[{e}_{\textrm{c}1},\cdots ,{e}_{\textrm{c}n}\right]}^{{\top }}$ is the error vector, and ${e}_{\textrm{c}i}={x}_{i}-{z}_{\textrm{c}i},i=\textrm{1,2},\cdots ,n$ . The matrix $\mathbf{l}=[{l}_{\textrm{c}1},\cdots ,{l}_{\textrm{c}n}]$ contains adjustable coefficients, $\lambda $ is the filter factor, $\widehat{f}$ is the estimate of disturbance.

3.5 High-order differentiator

Given the desired signals ${y}_{\textrm{d}}$ , the differential and high-order derivative information can be extracted. The designed $n$ th order HOD in Ref. (Reference Li, Qi, Guo, Chen and Zhao27) is

(24) \begin{align}\left\{\begin{array}{l}{\dot{z}}_{\textrm{h}1}={z}_{\textrm{h}2}+{l}_{\textrm{h}1}({y}_{\textrm{d}}-{z}_{\textrm{h}1})\\[3pt] {\dot{z}}_{\textrm{h}2}={z}_{\textrm{h}3}+{l}_{\textrm{h}2}({y}_{\textrm{d}}-{z}_{\textrm{h}1})\\[3pt] \textrm{}\textrm{}\textrm{}\textrm{} \vdots \\[3pt] {\dot{z}}_{\textrm{h}(n+1)}={z}_{\textrm{h}n}+{l}_{\textrm{h}n}({y}_{\textrm{d}}-{z}_{\textrm{h}1})\end{array}\right.\end{align}

where ${z}_{\textrm{h}1},{z}_{\textrm{h}2},\cdots ,{z}_{\textrm{h}1(n+1)}$ represent the internal states of the HOD, ${l}_{\textrm{h}1},{l}_{\textrm{h}2},\cdots ,{l}_{\textrm{h}n}$ the parameters to be adjusted. At time ${t}_{0}$ , ${y}_{\textrm{d}}({t}_{0}+{\Delta }t)$ is expressed using its Taylor series expansion

(25) \begin{align}\begin{array}{l}{y}_{\textrm{d}}({t}_{0}+{\Delta }t)={y}_{\textrm{d}}\left({t}_{0}\right)+{\dot{y}}_{\textrm{d}}\left({t}_{0}\right){\Delta }t+\cdots +\frac{{y}_{\textrm{d}}^{\left(n\right)}\left({t}_{0}\right){\Delta }{t}^{\left(n\right)}}{n!}+o\left({\Delta }{t}^{(n+1)}\right)\end{array}\end{align}

where $o\left({\Delta }{t}^{(n+1)}\right)$ is the infinitesimal term.

In practice, the signal sampling time in vector quadrotors is short. Thus, the higher-order terms $o\left({\Delta }{t}^{(n+1)}\right)$ can be neglected. Assumption 2 includes constant, first-order, and second-order signals, as well as their combinations, encompassing a broad range of cases.

3.6 Closed-loop system analysis

In practice, ${f}_{1}$ is replaced by ${\widehat{f}}_{1}$ , which is observed by first-order CFO and ${\widehat{f}}_{2}$ is observed by second-order CFO. Therefore, for attitude and position control, we have

(26) \begin{align}\left\{\begin{array}{ll}{\mathit{\tau }}_{\textrm{d}} {} &=-{k}_{\mathbf{r}}{\mathbf{e}}_{\mathbf{r}}-{k}_{\mathit{\omega }}{\mathbf{e}}_{\mathit{\omega }}+{\mathit{\omega }}_{\textrm{B}}\times {\mathbf{I}}_{\textrm{B}}{\mathit{\omega }}_{\textrm{B}}-{\mathbf{I}}_{\textrm{B}}\left({\widehat{\mathit{\omega }}}_{\textrm{B}}{\mathbf{r}}^{{\top }}{\mathbf{r}}_{\textrm{d}}{\mathit{\omega }}_{\textrm{d}}-{\mathbf{r}}^{{\top }}{\mathbf{r}}_{\textrm{d}}{\dot{\mathit{\omega }}}_{\textrm{d}}-{\widehat{\mathbf{f}}}_{1}\right)\\[6pt] \mathbf{f} {} &=\left({k}_{\mathbf{p}}{\mathbf{e}}_{\mathbf{p}}+{k}_{\dot{\mathbf{p}}}{\mathbf{e}}_{\dot{\mathbf{p}}}+mg{\mathbf{e}}_{3}-m{\ddot{\mathbf{p}}}_{\textrm{d}}-{\widehat{\mathbf{f}}}_{2}\right){\mathbf{R}}_{\textrm{B}\textrm{W}}{\mathbf{e}}_{3}\end{array}\right.\end{align}

Usually, when the requirements for control are not stringent, ${\widehat{f}}_{1}$ and ${\widehat{f}}_{2}$ can be ignored. But as UAV velocity increases, the effect of ${\widehat{f}}_{1}$ and ${\widehat{f}}_{2}$ on the UAV will become more severe. Estimation and compensation by observers are simple and effective.

Taylor expansion of f based on Assumption 1, we have

(27) \begin{align}\left\{\begin{array}{ll}f({t}_{0}+{\Delta }t) {} &=f\left({t}_{0}\right)+\dot{f}\left({t}_{0}\right){\Delta }t+\frac{1}{2!}\ddot{f}\left({t}_{0}\right)({\Delta }t{)}^{2}+o\left(({\Delta }t{)}^{3}\right)\\[6pt] {y}_{d}({t}_{0}+{\Delta }t) {} &={y}_{d}\left({t}_{0}\right)+{\dot{y}}_{d}\left({t}_{0}\right){\Delta }t+\frac{1}{2!}{\ddot{y}}_{d}\left({t}_{0}\right)({\Delta }t{)}^{2}+o\left(({\Delta }t{)}^{3}\right)\end{array}\right.\end{align}

Equation (16) is a second-order controlled system. Its main problems are to solve the effect of second-order disturbances. For ordinary UAV, (16) can be small enough to be ignored due to the short sampling time and boundary $f$ . On the other hand, when higher control precisions are required, higher-order CFO also can be chosen. Here, use of second-order CFO is based on practicality and simplicity of analysis. For ${y}_{d}({t}_{0}+{\Delta }t)$ , it is usually given by the remote. When the sampling period of the UAV is short, it is reasonable to ignore $o\left(({\Delta }t{)}^{3}\right)$ .

Based on Lemma 1 and Assumption 3, second-order CFO achieves unstable estimates for any $f$ . That is ${\textrm{l}\textrm{i}\textrm{m}}_{x\to 0}\widehat{f}=f$ . Based on Lemma 2 and Assumption 3, second-order HOD achieves no steady-state error estimates for any ${y}_{d}$ .

Proof. We aim to prove the stability of the following control law:

(28) \begin{align}\mathbf{f}=\left({k}_{\mathbf{p}}{\mathbf{e}}_{\mathbf{p}}+{k}_{\dot{\mathbf{p}}}{\mathbf{e}}_{\dot{\mathbf{p}}}+mg{\mathbf{e}}_{3}-m{\ddot{\mathbf{p}}}_{\textrm{d}}-{\widehat{\mathbf{f}}}_{2}\right){\mathbf{R}}_{\textrm{B}\textrm{W}}{\mathbf{e}}_{3}\end{align}

using Lyapunov stability theory, combined with the analysis of the system’s error dynamics. Here, ${\mathbf{e}}_{\mathbf{p}}=\mathbf{p}-{\mathbf{p}}_{\textrm{d}}$ denotes the position error, ${\mathbf{e}}_{\dot{\mathbf{p}}}=\dot{\mathbf{p}}-{\dot{\mathbf{p}}}_{\textrm{d}}$ denotes the velocity error, ${\widehat{\mathbf{f}}}_{2}$ is the estimated disturbance, and ${\mathbf{R}}_{\textrm{B}\textrm{W}}{\mathbf{e}}_{3}$ represents the thrust direction in the world frame.

Substituting the control law into the dynamic equation of the system (assuming the dynamic equation is $m\ddot{\mathbf{p}}={\mathbf{f}}_{2}+{\widehat{\mathbf{f}}}_{2}-mg{\mathbf{e}}_{3}$ , where $mg{\mathbf{e}}_{3}$ is the gravity term):

(29) \begin{align}m\ddot{\mathbf{p}}=\left({k}_{\mathbf{p}}{\mathbf{e}}_{\mathbf{p}}+{k}_{\dot{\mathbf{p}}}{\mathbf{e}}_{\dot{\mathbf{p}}}+mg{\mathbf{e}}_{3}-m{\ddot{\mathbf{p}}}_{\textrm{d}}-{\widehat{\mathbf{f}}}_{2}\right){\mathbf{R}}_{\textrm{B}\textrm{W}}{\mathbf{e}}_{3}+{\mathbf{f}}_{2}-mg{\mathbf{e}}_{3}\end{align}

Rearranging for the second derivative of the position error ( ${\ddot{\mathbf{e}}}_{\textrm{p}}=\ddot{\mathbf{p}}-{\ddot{\mathbf{p}}}_{\textrm{d}}$ ):

(30) \begin{align}m{\ddot{\mathbf{e}}}_{\textrm{p}}=m\left(\ddot{\mathbf{p}}-{\ddot{\mathbf{p}}}_{\textrm{d}}\right)=\left({k}_{\textrm{p}}{\mathbf{e}}_{\textrm{p}}+{k}_{\textrm{d}}{\dot{\mathbf{e}}}_{\textrm{p}}+mg{\mathbf{e}}_{3}-{\widehat{\mathbf{f}}}_{2}\right){\mathbf{R}}_{\textrm{B}\textrm{W}}{\mathbf{e}}_{3}+{\mathbf{f}}_{2}-mg{\mathbf{e}}_{3}\end{align}

For simplified analysis, assume ${\mathbf{R}}_{\textrm{B}\textrm{W}}{\mathbf{e}}_{3}={\mathbf{e}}_{3}$ (i.e. the thrust direction aligns with the world z-axis), then:

(31) \begin{align}m{\ddot{\mathbf{e}}}_{\textrm{p}}={k}_{\textrm{p}}{\mathbf{e}}_{\textrm{p}}+{k}_{\textrm{d}}{\dot{\mathbf{e}}}_{\textrm{p}}+mg{\mathbf{e}}_{3}-{\widehat{\mathbf{f}}}_{2}+{\mathbf{f}}_{2}-mg{\mathbf{e}}_{3}\end{align}

The gravity terms cancel, giving:

(32) \begin{align}m{\ddot{\mathbf{e}}}_{\textrm{p}}={k}_{\textrm{p}}{\mathbf{e}}_{\textrm{p}}+{k}_{\textrm{d}}{\dot{\mathbf{e}}}_{\textrm{p}}+\left({\mathbf{f}}_{2}-{\widehat{\mathbf{f}}}_{2}\right)\end{align}

If the disturbance estimation is accurate ( ${\widehat{\mathbf{f}}}_{2}\to {\mathbf{f}}_{2}$ ), then:

(33) \begin{align}m{\ddot{\mathbf{e}}}_{\textrm{p}}+{k}_{\textrm{d}}{\dot{\mathbf{e}}}_{\textrm{p}}+{k}_{\textrm{p}}{\mathbf{e}}_{\textrm{p}}=0\end{align}

We choose the Lyapunov function as:

(34) \begin{align}V=\frac{1}{2}m{\dot{\mathbf{e}}}_{\textrm{p}}^{{\top }}{\dot{\mathbf{e}}}_{\textrm{p}}+\frac{1}{2}{k}_{\textrm{p}}{\mathbf{e}}_{\textrm{p}}^{{\top }}{\mathbf{e}}_{\textrm{p}}\end{align}

The time derivative is:

(35) \begin{align}\dot{V}=m{\dot{\mathbf{e}}}_{\textrm{p}}^{{\top }}{\ddot{\mathbf{e}}}_{\textrm{p}}+{k}_{\textrm{p}}{\mathbf{e}}_{\textrm{p}}^{{\top }}{\dot{\mathbf{e}}}_{\textrm{p}}\end{align}

Substituting error dynamics:

(36) \begin{align}m{\ddot{\mathbf{e}}}_{\textrm{p}}=-{k}_{\textrm{d}}{\dot{\mathbf{e}}}_{\textrm{p}}-{k}_{\textrm{p}}{\mathbf{e}}_{\textrm{p}}\end{align}

gives:

(37) \begin{align}\dot{V}={\dot{\mathbf{e}}}_{\textrm{p}}^{{\top }}(\!\!-{k}_{\textrm{d}}{\dot{\mathbf{e}}}_{\textrm{p}}-{k}_{\textrm{p}}{\mathbf{e}}_{\textrm{p}})+{k}_{\textrm{p}}{\mathbf{e}}_{\textrm{p}}^{{\top }}{\dot{\mathbf{e}}}_{\textrm{p}}\end{align}

Simplifying:

(38) \begin{align}\dot{V}=-{k}_{\textrm{d}}{\dot{\mathbf{e}}}_{\textrm{p}}^{{\top }}{\dot{\mathbf{e}}}_{\textrm{p}}\end{align}

Since ${k}_{\textrm{d}} \gt 0$ , we have $\dot{V}\le 0$ . Therefore, the system is globally stable in the sense of Lyapunov, and the error dynamics are asymptotically stable.

Proof. The attitude dynamics are

(39) \begin{align}{\mathbf{I}}_{\textrm{B}}{\dot{\mathit{\omega }}}_{\textrm{B}}=-{\mathit{\omega }}_{\textrm{B}}\times \left({\mathbf{I}}_{\textrm{B}}{\mathit{\omega }}_{\textrm{B}}\right)+{\mathit{\tau }}_{\textrm{d}}+{\mathbf{f}}_{1}.\end{align}

Define the attitude and angular velocity errors by

(40) \begin{align}{\mathbf{e}}_{\mathbf{r}}=\frac{1}{2}{\left({\mathbf{R}}_{\textrm{d}}^{{\top }}\mathbf{R}-{\mathbf{R}}^{{\top }}{\mathbf{R}}_{\textrm{d}}\right)}^{\vee }, \end{align}

(41) \begin{align}{\mathbf{e}}_{\mathit{\omega }}={\mathit{\omega }}_{\textrm{B}}-{\mathbf{R}}^{{\top }}{\mathbf{R}}_{\textrm{d}}{\mathit{\omega }}_{\textrm{d}}. \\[6pt] \nonumber \end{align}

They satisfy

(42) \begin{align} {\dot{\mathbf{e}}}_{\mathbf{r}}=\mathbf{E}(\mathbf{R},{\mathbf{R}}_{\textrm{d}})\textrm{}{\mathbf{e}}_{\mathit{\omega }}, \\[-8pt] \nonumber \end{align}

(43) \begin{align} \parallel \mathbf{E}(\mathbf{R},{\mathbf{R}}_{\textrm{d}})\parallel \le 1. \\[9pt] \nonumber \end{align}

Consider the control law

(44) \begin{align}{\mathit{\tau }}_{\textrm{d}}=-{k}_{\mathbf{r}}\textrm{}{\mathbf{e}}_{\mathbf{r}}-{k}_{\mathit{\omega }}\textrm{}{\mathbf{e}}_{\mathit{\omega }}+{\mathit{\omega }}_{\textrm{B}}\times \left({\mathbf{I}}_{\textrm{B}}{\mathit{\omega }}_{\textrm{B}}\right)-{\mathbf{I}}_{\textrm{B}}\left({\widehat{\mathit{\omega }}}_{\textrm{B}}{\mathbf{R}}^{{\top }}{\mathbf{R}}_{\textrm{d}}{\mathit{\omega }}_{\textrm{d}}-{\mathbf{R}}^{{\top }}{\mathbf{R}}_{\textrm{d}}{\dot{\mathit{\omega }}}_{\textrm{d}}\right)-{\widehat{\mathbf{f}}}_{1},\end{align}

where ${k}_{\mathbf{r}},{k}_{\mathit{\omega }}\in {\mathbb{R}}^{3\times 3}$ are positive definite gains and ${\widehat{\mathbf{f}}}_{1}$ is the disturbance estimate.

Let the estimation error be

(45) \begin{align}{\stackrel{\sim }{\mathbf{f}}}_{1}={\mathbf{f}}_{1}-{\widehat{\mathbf{f}}}_{1}.\end{align}

Substituting the control into the dynamics yields the error system

(46) \begin{align}\left\{\begin{array}{l}{\dot{\mathbf{e}}}_{\mathbf{r}}=\mathbf{E}\textrm{}{\mathbf{e}}_{\mathit{\omega }},\\ {\dot{\mathbf{e}}}_{\mathit{\omega }}=-{\mathbf{I}}_{\textrm{B}}^{-1}{k}_{\mathbf{r}}\textrm{}{\mathbf{e}}_{\mathbf{r}}-{\mathbf{I}}_{\textrm{B}}^{-1}{k}_{\mathit{\omega }}\textrm{}{\mathbf{e}}_{\mathit{\omega }}+{\Delta }+{\mathbf{I}}_{\textrm{B}}^{-1}{\stackrel{\sim }{\mathbf{f}}}_{1},\end{array}\right.\end{align}

with

(47) \begin{align} {\Delta }\triangleq -{\widehat{\mathit{\omega }}}_{\textrm{B}}{\mathbf{R}}^{{\top }}{\mathbf{R}}_{\textrm{d}}{\mathit{\omega }}_{\textrm{d}} &+ {\mathbf{R}}^{{\top }}{\mathbf{R}}_{\textrm{d}}{\dot{\mathit{\omega }}}_{\textrm{d}}-{\dot{\mathit{\omega }}}_{\textrm{d}}, \\[-8pt] \nonumber \end{align}

(48) \begin{align} \parallel {\Delta }\parallel &\le \delta . \\[9pt] \nonumber\end{align}

Choose the Lyapunov function

(49) \begin{align} V={k}_{\mathbf{r}}\textrm{}{\Psi }(\mathbf{R},{\mathbf{R}}_{\textrm{d}})+\frac{1}{2}{\mathbf{e}}_{\mathit{\omega }}^{{\top }}{\mathbf{I}}_{\textrm{B}}{\mathbf{e}}_{\mathit{\omega }}+\frac{1}{2}{\stackrel{\sim }{\mathbf{f}}}_{1}^{{\top }}{{\Gamma }}^{-1}{\stackrel{\sim }{\mathbf{f}}}_{1}, \\[-10pt] \nonumber \end{align}

(50) \begin{align} {\Psi }=\frac{1}{2}\textrm{t}\textrm{r}({\mathbf{I}}_{3}-{\mathbf{R}}_{\textrm{d}}^{{\top }}\mathbf{R}), \\[9pt] \nonumber\end{align}

with ${\Gamma } \gt 0$ . Using $\dot{{\Psi }}={\mathbf{e}}_{\mathbf{r}}^{{\top }}{\mathbf{e}}_{\mathit{\omega }}$ gives

(51) \begin{align}\begin{array}{ll}\dot{V} {} &= {k}_{\mathbf{r}}\textrm{}{\mathbf{e}}_{\mathbf{r}}^{{\top }}{\mathbf{e}}_{\mathit{\omega }}+{\mathbf{e}}_{\mathit{\omega }}^{{\top }}{\mathbf{I}}_{\textrm{B}}{\dot{\mathbf{e}}}_{\mathit{\omega }}+{\stackrel{\sim }{\mathbf{f}}}_{1}^{{\top }}{{\Gamma }}^{-1}{\dot{\stackrel{\sim }{\mathbf{f}}}}_{1}\\[3pt] &={k}_{\mathbf{r}}\textrm{}{\mathbf{e}}_{\mathbf{r}}^{{\top }}{\mathbf{e}}_{\mathit{\omega }}+{\mathbf{e}}_{\mathit{\omega }}^{{\top }}\left(-{k}_{\mathbf{r}}{\mathbf{e}}_{\mathbf{r}}-{k}_{\mathit{\omega }}{\mathbf{e}}_{\mathit{\omega }}+{\mathbf{I}}_{\textrm{B}}{\Delta }+{\stackrel{\sim }{\mathbf{f}}}_{1}\right)+{\stackrel{\sim }{\mathbf{f}}}_{1}^{{\top }}{{\Gamma }}^{-1}({\dot{\mathbf{f}}}_{1}-{\dot{\widehat{\mathbf{f}}}}_{1}).\end{array}\end{align}

Adopt the adaptation law

(52) \begin{align}{\dot{\widehat{\mathbf{f}}}}_{1}={\Gamma }\textrm{}{\mathbf{e}}_{\mathit{\omega }}.\end{align}

Then

(53) \begin{align}\dot{V}=-{\mathbf{e}}_{\mathit{\omega }}^{{\top }}{k}_{\mathit{\omega }}{\mathbf{e}}_{\mathit{\omega }}+{\mathbf{e}}_{\mathit{\omega }}^{{\top }}{\mathbf{I}}_{\textrm{B}}{\Delta }+{\stackrel{\sim }{\mathbf{f}}}_{1}^{{\top }}{{\Gamma }}^{-1}{\dot{\mathbf{f}}}_{1}.\end{align}

By Cauchy–Schwarz and Young inequalities, for any $\eta \gt 0$ ,

(54) \begin{align} {\mathbf{e}}_{\mathit{\omega }}^{{\top }}{\mathbf{I}}_{\textrm{B}}{\Delta }\le \frac{{\lambda }_{\textrm{m}\textrm{a}\textrm{x}}\left({\mathbf{I}}_{\textrm{B}}\right)}{2\eta }\parallel {\Delta }{\parallel }^{2}+\frac{\eta }{2}\parallel {\mathbf{e}}_{\mathit{\omega }}{\parallel }^{2}, \\[-8pt] \nonumber\end{align}

(55) \begin{align} {\stackrel{\sim }{\mathbf{f}}}_{1}^{{\top }}{{\Gamma }}^{-1}{\dot{\mathbf{f}}}_{1}\le \frac{1}{2}\parallel {\stackrel{\sim }{\mathbf{f}}}_{1}{\parallel }_{{{\Gamma }}^{-1}}^{2}+\frac{1}{2}\parallel {\dot{\mathbf{f}}}_{1}{\parallel }_{{{\Gamma }}^{-1}}^{2}. \\[9pt] \nonumber\end{align}

Hence

(56) \begin{align}\dot{V}\le -{\lambda }_{\textrm{m}\textrm{i}\textrm{n}}\left({k}_{\mathit{\omega }}\right)\parallel {\mathbf{e}}_{\mathit{\omega }}{\parallel }^{2}+\frac{\eta }{2}\parallel {\mathbf{e}}_{\mathit{\omega }}{\parallel }^{2}+\frac{{\lambda }_{\textrm{m}\textrm{a}\textrm{x}}\left({\mathbf{I}}_{\textrm{B}}\right)}{2\eta }{\delta }^{2}+\frac{1}{2}\parallel {\stackrel{\sim }{\mathbf{f}}}_{1}{\parallel }_{{{\Gamma }}^{-1}}^{2}+\frac{1}{2}\parallel {\dot{\mathbf{f}}}_{1}{\parallel }_{{{\Gamma }}^{-1}}^{2}.\end{align}

Choose

(57) \begin{align}\eta \lt 2{\lambda }_{\textrm{m}\textrm{i}\textrm{n}}\left({k}_{\mathit{\omega }}\right).\end{align}

Then

(58) \begin{align}\dot{V}\le -{c}_{\mathit{\omega }}\parallel {\mathbf{e}}_{\mathit{\omega }}{\parallel }^{2}+{c}_{{\Delta }}{\delta }^{2}+{c}_{f},\end{align}

for some constants ${c}_{\mathit{\omega }},{c}_{{\Delta }},{c}_{f} \gt 0$ . Together with

(59) \begin{align}{\dot{\mathbf{e}}}_{\mathbf{r}}=\mathbf{E}\textrm{}{\mathbf{e}}_{\mathit{\omega }},\end{align}

and $\parallel \mathbf{E}\parallel \le 1$ , it follows that $({\mathbf{e}}_{\mathbf{r}},{\mathbf{e}}_{\mathit{\omega }},{\stackrel{\sim }{\mathbf{f}}}_{1})$ is uniformly ultimately bounded.

If

(60) \begin{align}{\Delta }\equiv 0 \quad \textrm{a}\textrm{n}\textrm{d} \quad {\dot{\mathbf{f}}}_{1}\equiv 0,\end{align}

then

(61) \begin{align}\dot{V}\le -{c}_{\mathit{\omega }}\parallel {\mathbf{e}}_{\mathit{\omega }}{\parallel }^{2},\end{align}

which implies almost global asymptotic stability:

(62) \begin{align}\underset{t\to {\infty }}{\textrm{l}\textrm{i}\textrm{m}}{\mathbf{e}}_{\mathit{\omega }}=0, \quad \textrm{}\textrm{}\textrm{}\textrm{}\textrm{}\textrm{}\textrm{}\textrm{}\underset{t\to {\infty }}{\textrm{l}\textrm{i}\textrm{m}}{\mathbf{e}}_{\mathbf{r}}=0.\end{align}

3.7 Control allocation

This subsection presents an optimisation method for generating feasible tilt angles ( ${\alpha }_{i}$ ) and rotor speeds ( ${\omega }_{i}$ ) for the vector quadrotor system. As the vector quadrotor is an overactuated system, its control input allocation requires a systematic optimisation approach. The complete optimisation framework integrates cost functions and constraints, ultimately forming a quadratic programming problem with linear constraints (SQP with linear constraints), defined as follows:

(63) \begin{align}\begin{array}{ll}\underset{{\omega }_{i}}{\textrm{m}\textrm{i}\textrm{n}} \quad E={k}_{1}{\left|\left|{v}_{\textrm{B}}\right|\right|}^{2}+{k}_{2}{\left|\left|{\mathbf{e}}_{\mathbf{p}}\right|\right|}^{2}+{k}_{3}{\left|\left|{\mathbf{e}}_{\mathbf{r}}\right|\right|}^{2}\\[3pt] s.t. \quad mg \, {\le} \, {\left({\mathbf{R}}_{{\textrm{P}}_{i}\textrm{B}}{\mathbf{t}}_{{\textrm{p}}_{i}}\right)}_{3}\\[3pt] {} \left[\begin{array}{l}\ddot{\mathbf{p}}\\ {\dot{\mathit{\omega }}}_{\textrm{B}}\end{array}\right]={\mathbf{e}}_{6}g+\mathbf{J}\left(\alpha \right)\left[\begin{array}{l}\mathit{\omega }\\[3pt] {\mathit{\omega }}_{\alpha }\end{array}\right]\\[3pt] {} {\omega }_{\textrm{m}\textrm{i}\textrm{n}}\le {\omega }_{i}\le {\omega }_{\textrm{m}\textrm{a}\textrm{x}}\\[3pt] {} {\alpha }_{\textrm{m}\textrm{i}\textrm{n}}\le {\alpha }_{i}\le {\alpha }_{\textrm{m}\textrm{a}\textrm{x}}\\ {} \end{array}\end{align}

where ${v}_{\textrm{B}}=\sqrt{{v}_{x}^{2}+{v}_{y}^{2}+{v}_{z}^{2}}$ denotes the simplified kinetic energy of the vector quadrotor. Here, ${e}_{\mathbf{p}}$ and ${e}_{\mathbf{r}}$ represent the position error and attitude error, respectively, while ${k}_{1},{k}_{2},{k}_{3}$ are weight coefficients. The inequality constraint $mg{\le}{\left({\mathbf{R}}_{{\textrm{P}}_{i}\textrm{B}}{\mathbf{t}}_{{\textrm{p}}_{i}}\right)}_{3}$ is designed to ensure that the thrust generated by the rotor speed maintains the vector quadrotor in hovering equilibrium; the symbol ${\left({\mathbf{R}}_{{\textrm{P}}_{i}\textrm{B}}{\mathbf{t}}_{{\textrm{p}}_{i}}\right)}_{3}$ denotes the third element of the vector ${\mathbf{R}}_{{\textrm{P}}_{i}\textrm{B}}{\mathbf{t}}_{{\textrm{p}}_{i}}$ . The equality constraint is given by: $\left[\begin{array}{l}\ddot{\mathbf{p}}\\ {\dot{\mathit{\omega }}}_{\textrm{B}}\end{array}\right]={\mathbf{e}}_{6}g+\mathbf{J}\left(\alpha \right)\left[\begin{array}{l}\mathit{\omega }\\ {\mathit{\omega }}_{\alpha }\end{array}\right]$ where ${\mathbf{e}}_{6}={[\begin{array}{l@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}r}0&0&0&0&0&1\end{array}]}^{{\top }}$ , and the matrix $\mathbf{J}\left(\alpha \right)$ (referred to as the output Jacobian, with $\mathbf{J}\left(\alpha \right)\in {\mathbb{R}}^{6\times 8}$ ) is defined as $\mathbf{J}\left(\alpha \right)={\mathbf{J}}_{\textrm{R}}{\mathbf{J}}_{\alpha }\left(\alpha \right)$ . Here, ${\mathbf{J}}_{\textrm{R}}\in {\mathbb{R}}^{6\times 6}$ (also an output Jacobian) and ${\mathbf{J}}_{\alpha }\left(\alpha \right)$ are expressed as: ${\mathbf{J}}_{\textrm{R}}=\left[\begin{array}{l@{\quad}r}{m}^{-1}{\mathbf{R}}_{\textrm{BW}} {}& {0}^{3\times 3}\\ {0}^{3\times 3} {}& {\mathbf{I}}_{\textrm{B}}^{-1}\end{array}\right],\textrm{}\textrm{}\textrm{}\textrm{}{\mathbf{J}}_{\alpha }\left(\alpha \right)=\left[\begin{array}{l@{\quad}r}\mathbf{F}(\alpha)&{0}^{3\times 4}\\ \mathit{\tau}(\alpha)&{0}^{3\times 4}\end{array}\right]$ Note that ${e}_{\mathbf{r}}$ and ${e}_{\mathbf{p}}$ are defined in (18) and (21), respectively. For the vector quadrotor’s end-effector: ${\mathbf{r}}_{\textrm{d}}\left(t\right)$ and $\mathbf{r}\left(t\right)$ represent the desired and actual attitudes, while ${\mathbf{p}}_{\textrm{d}}\left(t\right)$ and $\mathbf{p}\left(t\right)$ represent the desired and actual positions. Additional constraints are introduced to ensure accurate computation of the tilt angles ${\alpha }_{i}$ and rotor speeds ${\omega }_{i}$ –proper constraints are critical for this purpose. Specifically, the constraints ${\omega }_{\textrm{m}\textrm{i}\textrm{n}}\le {\omega }_{i}\le {\omega }_{\textrm{m}\textrm{a}\textrm{x}}$ and ${\alpha }_{\textrm{m}\textrm{i}\textrm{n}}\le {\alpha }_{i}\le {\alpha }_{\textrm{m}\textrm{a}\textrm{x}}$ define the permissible ranges for ${\omega }_{i}$ (rotor speed) and ${\alpha }_{i}$ (tilt angle), respectively. It is important to note that the trajectory generation algorithm can be computed either offline or online. Given that quadratic optimisation problems can be efficiently solved using state-of-the-art algorithms, this approach is particularly suitable for online computation—even on hardware with limited computational resources. Online iteration allows for dynamic adjustment of optimisation parameters during specific mission phases, which facilitates the implementation of different redundancy resolution strategies by tuning cost function weights based on the mission phase or external triggers. This capability is demonstrated in the experimental section through a vector quadrotor mission.

3.8 Theoretical analysis of allocation algorithm

The core foundation of global convergence lies in the intrinsic properties of the optimisation problem, which determine the existence, uniqueness and reachability of the optimal solution.

The objective function is defined as

(64) \begin{align}E={k}_{1}{||{v}_{\textrm{B}}||}^{2}+{k}_{2}{||{\mathbf{e}}_{\mathbf{p}}||}^{2}+{k}_{3}{||{\mathbf{e}}_{\mathbf{r}}||}^{2}\end{align}

where the weight coefficients ${k}_{1},{k}_{2},{k}_{3} \gt 0$ . By definition, squared norm terms are strictly convex functions, and the weighted sum of strictly convex functions remain strictly convex. This implies that the objective function has a unique minimum point in its domain.

The feasible region ${\Omega }$ is defined by linear constraints, including equality constraints describing system dynamics, and inequality constraints such as rotor speed limits, tilt angle bounds and thrust requirements. Linear constraints inherently form convex sets, and the intersection of convex sets is also a convex set. Thus, ${\Omega }$ is a convex set.

Due to physical limits on rotor speeds and tilt angles:

(65) \begin{align} {\omega }_{\textrm{m}\textrm{i}\textrm{n}}\le {\omega }_{i}\le {\omega }_{\textrm{m}\textrm{a}\textrm{x}},\textrm{}\textrm{}\textrm{}\textrm{}{\alpha }_{\textrm{m}\textrm{i}\textrm{n}}\le {\alpha }_{i}\le {\alpha }_{\textrm{m}\textrm{a}\textrm{x}}\end{align}

the feasible region ${\Omega }$ is bounded. Meanwhile, sets defined by linear constraints are necessarily closed. A ‘bounded and closed’ set in Euclidean space is a compact set.

For the problem of minimising a strictly convex function over a non-empty convex set, there exists at most one global optimal solution. Furthermore, according to the Weierstrass extreme value theorem, a continuous function (here, the strictly convex objective function $E$ ) must attain its minimum on a compact set. Combining these two points, the optimisation problem admits a unique global optimal solution.

The proposed problem is essentially a convex SQP problem, which can be efficiently solved using state-of-the-art SQP solvers (e.g. interior-point methods or active-set methods). The global convergence of these solvers for convex SQP problems has been fully verified.

Interior-point methods transform constrained problems into a series of unconstrained subproblems using barrier functions. For convex SQP, each subproblem remains strictly convex, and the sequence of solutions generated by interior-point methods converges to the global optimal solution with a quadratic convergence rate, independent of the initial point selection.

Active-set methods iteratively identify ‘active constraints’ (i.e. constraints that hold with equality at the optimal solution) and reduce the original problem to a low-dimensional unconstrained SQP problem. For convex SQP with a compact feasible region, the number of active-set updates is finite, ensuring that the algorithm converges to solution in a finite number of steps.

When the algorithm is used for online optimisation (e.g. real-time adjustment of weight coefficients), global convergence remains valid:

• • Weight adjustments only alter the ‘shape’ of the objective function without violating its strict convexity (the Hessian matrix of the new objective function remains positive definite);

• • The feasible region ${\Omega }$ remains unchanged (still a compact set), so the new problem still admits a unique global optimal solution;

• • Smooth adjustments of weights lead to continuous changes in the optimal solution, and the iterative sequence can stably track these changes, ensuring convergence in online scenarios.

The optimisation framework for the vector quadrotor system is a convex quadratic programming problem with a compact feasible region and a unique global optimal solution. When using standard SQP algorithms such as interior-point methods or active-set methods, the iterative sequence converges to this global optimal solution regardless of the initial point. Convergence is maintained even in scenarios with online parameter adjustments. Therefore, the algorithm theoretically possesses global convergence.

4.1 Effect of control algorithm

To verify the feasibility of the control algorithm, MCC controller combine with SQP based on SO(3). In the subsection, we calculate the trajectory and attitude of tracking a stationary target using the Lyapunov vector field navigation method [Reference Muslimov and Munasypov36]. The simulation considers a vector quadrotor equipped with a fixed-angle camera to track a stationary target. We limit the maximum yaw rate of the vector quadrotor to ${\dot{\psi }}_{\textrm{m}\textrm{a}\textrm{x}}=0.1\textrm{}\textrm{}\textrm{r}\textrm{a}\textrm{d}/\textrm{s}$ and a nominal speed of ${v}_{0}=15\textrm{}\textrm{}\textrm{m}/\textrm{s}$ . The desired angle ${\theta }_{\textrm{d}}$ is determined by the distance between the vector quadrotor and the tracked target, given by ${\theta }_{\textrm{d}}=\textrm{a}\textrm{t}\textrm{a}\textrm{n}\left(\frac{z}{x}\right)$ assuming no background wind, where $z,x$ representing the altitude and horizontal displacement of the vector quadrotor, respectively. While tracking the target, the vector quadrotor obtains the target’s velocity and position using the onboard camera. To ensure effective camera operation and keep the target centred in its field of view, the distance between the camera and the target should be maintained within a specified range.

Considering that the vector quadrotor’s flight altitude is constant and simplify the problem, this paper requires the vector quadrotor to maintain an approximately constant standoff distance from the target on the horizontal plane while tracking the target. The vector quadrotor starts from the position $x=600\textrm{}\textrm{}\textrm{m},y=600\textrm{}\textrm{}\textrm{m},z=100\textrm{}\textrm{}\textrm{m}$ with an initial speed of ${v}_{x}=15\textrm{}\textrm{}\textrm{m}/\textrm{s},{v}_{y}=-15\textrm{}\textrm{}\textrm{m}/\textrm{s},{v}_{z}=0\textrm{}\textrm{}\textrm{m}/\textrm{s}$ and tracks a stationary target located at $x=100\textrm{}\textrm{}\textrm{m},y=100\textrm{}\textrm{}\textrm{m},z=0\textrm{}\textrm{}\textrm{m}$ . The standoff distance is set to $150\textrm{m}$ .

Figure 3 shows the trajectory of the vector quadrotor in the $\left(x,y\right)$ coordinate plane, which represents the variables in $x$ and $y$ , respectively. Vector quadrotor starting from $x=600\textrm{}\textrm{}\textrm{m},y=600\textrm{}\textrm{}\textrm{m}$ , vector quadrotor’s $x$ and $y$ values oscillate between $x=-106\textrm{}\textrm{}\textrm{m},y=306\textrm{}\textrm{}\textrm{m}$ after entering the standoff range. Adjusting the vector quadrotor’s required flight speed, a single vector quadrotor can successfully track a stationary ground target under heading control.

Figure 3.

The changes in position and attitude of vector quadrotor when tracking stationary targets.

Figure 3(a) indicates that the vector quadrotor can correctly enter the standoff range and remain stable within this range. Figure 3(b) shows that the yaw angle $\psi $ of the vector quadrotor continuously increases while the roll angle $\varphi $ increases during the initial stage. Once the vector quadrotor enters the standoff range, it maintains stability. Clearly, vector quadrotor can accomplish the task of tracking stationary targets.

Figure 4(a) illustrates the temporal evolution of the four tilt angles ${\alpha }_{i}$ , while Figure 4(b) presents the corresponding rotor speeds ${\omega }_{i}$ , both obtained by numerically solving (63) using the proposed control allocation framework. As the system initiates motion, significant fluctuations can be observed in both ${\alpha }_{i}$ and ${\omega }_{i}$ , reflecting the system’s dynamic response to the initial conditions and external disturbances. As the vector quadrotor approaches and enters the standoff range, these fluctuations progressively diminish, indicating convergence and stabilisation of the control variables. This behaviour demonstrates the system’s ability to self-stabilise under the designed control strategy.

Figure 4.

The tilt angle ${\alpha }_{i}$ of the servo and the rotor speed ${\omega }_{i}$ are calculated using quadratic programming algorithms.

The above results collectively confirm the feasibility and effectiveness of using a quadratic programming method to simultaneously determine the tilt angles and rotor speeds. Moreover, they validate the precision of the dynamic model in predicting the system’s physical outputs, thereby laying a solid foundation for real-time implementation.

4.2 Comparison of observer effects

In this subsection, we will compare the effectiveness of different observers. We consider four methods: MCC, $\textrm{S}\textrm{O}\left(3\right)$ -based MCC, ADRC and $\textrm{S}\textrm{O}\left(3\right)$ -based ADRC. The $\textrm{S}\textrm{O}\left(3\right)$ -based MCC, which is applied to vector quadrotor, is proposed in this paper. The key difference between the $\textrm{S}\textrm{O}\left(3\right)$ -based ADRC method and the $\textrm{S}\textrm{O}\left(3\right)$ -based MCC lies in the observer and differentiator. The former employs LESO and TD, while MCC utilises HOD and CFO. LESO and CFO in the controller share the same characteristic equation in both experiments and simulations. The parameter selection is determined based on the analysis of characteristic equations. The dynamic performance of CFO and HOD is determined by the poles of the characteristic equations. The principle of pole selection is to balance overshoot, delay time and settling time. In the simulations, we took these factors into account. The characteristic equation of the first-order HOD is given by ${s}^{2}+{l}_{\textrm{h}1}s+{l}_{\textrm{h}2}=0$ $({l}_{\textrm{h}1}=15,{l}_{\textrm{h}2}=75)$ . This ensures that the poles of HOD are located at 5 and 15. Similarly, the characteristic equation of the second-order HOD is given by ${s}^{3}+{l}_{\textrm{h}1}{s}^{2}+{l}_{\textrm{h}2}s+{l}_{\textrm{h}3}=0$ $({l}_{\textrm{h}1}=8.25,{l}_{\textrm{h}2}=22.695,{l}_{\textrm{h}3}=20.80)$ . In a similar way, the characteristic equations of the LESO and the first-order CFO are given by ${s}^{2}+\lambda s+{l}_{\textrm{c}1}=0$ $({l}_{\textrm{c}1}=27,\lambda =12)$ . The characteristic equation of the second-order CFO is given by ${s}^{3}+{l}_{\textrm{c}1}{s}^{2}+(\lambda {l}_{\textrm{c}1}+{l}_{\textrm{c}2}s+\lambda {l}_{\textrm{c}2}=0)$ . For all second-order CFOs in simulations, we set ${l}_{\textrm{c}1}=63.48,{l}_{\textrm{c}2}=20.7,\lambda =3.07$ . Next, the controller parameters should be tuned based on the observer’s performance. The parameters ${k}_{\mathbf{p}},{k}_{\dot{\mathbf{p}}},{k}_{\mathbf{r}}$ and ${k}_{\mathit{\omega }}$ offer large adjustment margins, ensuring improved performance over a wide range of conditions.

The initial state at $t=0$ is given by $\mathbf{p}=\dot{\mathbf{p}}=\mathbf{r}={\mathit{\omega }}_{\textrm{B}}={0}^{3\times 1}$ . Given a disturbance ${f}_{\theta }={f}_{x}=0.5$ applied at $10\textrm{s}$ . Given persistent Gaussian noise in the interval $15\textrm{s}-25\textrm{s}$ , with a mean of $\mu =0$ and a variance of $\sigma =0.01$ . CFO and LESO independently estimate the unknown disturbance in the pitch channel. The comparative performance of the observers is illustrated in Fig. 5(a). For ${f}_{\theta }=0.5$ , the maximum estimation error of LESO is 0.134, while that of CFO is 0.0532. For Gaussian noise, the mean absolute error (MAE) of LESO is 0.0450, while that of CFO is 0.0294. CFO exhibits superior performance. Moreover, increasing the error gain of LESO can reduce both delay and steady-state error. However, this results in high-frequency oscillations. The convergence of LESO is inherently constrained by its structural design. Lemma 1 states that the second-order CFO can achieve zero error estimation for both slope and acceleration signals. However, LESO provides a bounded estimation of the slope signal but diverges when estimating higher-order signals (Qi et al., 2024). In practice, first-order or second-order CFOs are preferred due to their sufficiently high estimation accuracy. Higher-order CFOs require additional sensor information, which is often impractical for vector quadrotors.

Figure 5.

Performance comparison of different controllers under disturbance without expectation signals.

The quality of observation directly impacts controller performance. The variation in pitch angle is illustrated in Fig. 5(b). For ${f}_{\theta }=0.5$ , the maximum pitch angles for ADRC, $\textrm{S}\textrm{O}\left(3\right)$ -based ADRC, MCC and $\textrm{S}\textrm{O}\left(3\right)$ -based MCC are 0.0661 rad, 0.0517 rad, 0.027 rad and 0.0106 rad, respectively. For Gaussian noise, the MAEs for ADRC, $\textrm{S}\textrm{O}\left(3\right)$ -based ADRC, MCC and $\textrm{S}\textrm{O}\left(3\right)$ -based MCC are 0.0262rad, 0.0241rad, 0.0055rad and 0.0026rad, respectively. Therefore, $\textrm{S}\textrm{O}\left(3\right)$ -based MCC exhibits superior anti-disturbance performance against high-frequency and step disturbances. The variation of $x$ is illustrated in Fig. 5(c). For $x$ , the MAEs for ADRC, $\textrm{S}\textrm{O}\left(3\right)$ -based ADRC, MCC, and $\textrm{S}\textrm{O}\left(3\right)$ -based MCC are 0.1167rad, 0.0784rad, 0.0559rad and 0.0334rad, respectively. On the other hand, for pitch angle and $x$ , Fig. 5(b) and (c) further demonstrate that $\textrm{S}\textrm{O}\left(3\right)$ -based MCC has a faster response time compared to other methods. Overall, the $\textrm{S}\textrm{O}\left(3\right)$ -based MCC method demonstrates superior anti-disturbance performance in both position and attitude control.

4.3 Comparison of tracking effects between four controllers

The vector quadrotor is subjected to a disturbance given by ${f}_{x}={f}_{\varphi }=0.6\textrm{s}\textrm{i}\textrm{n}(t+0.5)+\textrm{s}\textrm{i}\textrm{n}(0.8t+0.5)+\textrm{c}\textrm{o}\textrm{s}\left(1.2t\right)$ and the initial state at $t=0s$ is $\mathbf{p}=\mathbf{v}=\mathbf{r}={\mathit{\omega }}_{\textrm{B}}={0}^{3\times 1}$ . At the same time, the desired position for $x$ is specified, as shown in Fig. 6(b). Figure 6(a) illustrates the performance of LESO and CFO in estimating the disturbance. The MAE of LESO is 0.6834, while that of CFO is 0.0455. The estimation performance of CFO is approximately 13.48 times better than that of LESO. Figure 6(b) illustrates the performance of different controllers in tracking the same desired signal. Intuitively, control methods based on LESO exhibit higher delays and larger errors when tracking sine wave signals. MCC method tracks sine signals with low delay and minimal error but exhibits larger overshoots and a longer settling time when responding to step signals. For step signal tracking, the settling times of ADRC, $\textrm{S}\textrm{O}\left(3\right)$ -based LESO, MCC, and $\textrm{S}\textrm{O}\left(3\right)$ -based MCC are 3.46, 4.23, 2 and 0.32 seconds, respectively. Meanwhile, compared to MCC, the overshoot of $\textrm{S}\textrm{O}\left(3\right)$ -based MCC is reduced by approximately 35%.

Figure 6.

Observation errors of CFO and LESO and trajectory tracking effects of four different controllers.

We evaluate the performance of MCC and $\textrm{S}\textrm{O}\left(3\right)$ -based MCC methods in tracking step signals with amplitudes of 2m, 3m, 4m and 5m, each lasting for 3 seconds. The tracking performance of $x$ is illustrated in Fig. 7(a) and (b). For the two methods, the settling time does not differ significantly. For a step signal of 2m/s, the overshoot of $\textrm{S}\textrm{O}\left(3\right)$ -based MCC was 6.9%, while that of MCC was 22%. The overshoot of MCC method is eight times greater than that of the $\textrm{S}\textrm{O}\left(3\right)$ -based MCC. For a step signal of 3m/s, the overshoot of $\textrm{S}\textrm{O}\left(3\right)$ -based MCC was 6.9%, while that of MCC was 21.2%. The overshoot of MCC method is 11.5 times greater than that of the $\textrm{S}\textrm{O}\left(3\right)$ -based MCC. For a step signal of 4m/s, the overshoot of $\textrm{S}\textrm{O}\left(3\right)$ -based MCC was 6.9%, while that of MCC was 21.8%. The overshoot of MCC method is 14.2 times greater than that of the $\textrm{S}\textrm{O}\left(3\right)$ -based MCC. For a step signal of 5m/s, the overshoot of $\textrm{S}\textrm{O}\left(3\right)$ -based MCC was 6.9%, while that of MCC was 21.9%. The overshoot of MCC method is five times greater than that of the $\textrm{S}\textrm{O}\left(3\right)$ -based MCC. Therefore, the tracking precision of the $\textrm{S}\textrm{O}\left(3\right)$ -based MCC method is significantly higher.

Figure 7.

Position comparison of different controllers under disturbance with expectation signals.

4.4 Comparison of optimisation algorithm performance

The experimental environment is built on the Windows 11 operating system. In terms of hardware configuration, it is equipped with an Intel Core i5-13600KF processor and an NVIDIA GeForce RTX 3080 graphics card (12 GB memory), which can provide strong computing and graphics processing capabilities for the experiment. The development and operation of the experiment are mainly completed relying on the Simulink platform, which can effectively support experimental procedures such as model construction and simulation analysis.

4.4.1 Optimisation of objective functions by different algorithms

In the evaluation of algorithm performance, convergence time is one of the key indicators for measuring the efficiency and effectiveness of an algorithm. In particular, the convergence time to reach the final stable value directly reflects the time consumed by the algorithm to produce reliable and unchanging results. This study focuses on three algorithms, namely SQP, MPC and PINV, and conducts an in-depth comparison of their convergence times to reach the final stable value, aiming to comprehensively understand the performance characteristics and applicable scenarios of each algorithm.

In the Fig. 8, from these values, it can be intuitively observed that the SQP algorithm has the shortest convergence time to reach the final stable value, which is only 1.03. The convergence time of the MPC algorithm is 1.23, placing it in an intermediate position. The PINV algorithm has the longest convergence time, which is 1.70. To more clearly illustrate the differences in convergence speeds among the algorithms, we calculated the ratio of each algorithm’s convergence time to that of the fastest algorithm, SQP. The convergence time of MPC is approximately 1.194 times that of SQP, and that of PINV is approximately 1.650 times. This multiplicative relationship further emphasises the advantage of the SQP algorithm in terms of convergence speed. In the application of vector quadrotor, this advantage enables the reduction of decision-making waiting time, allowing the vector quadrotor to seize fleeting opportunities.

Figure 8.

Optimisation of objective functions by different algorithms.

To conduct an in-depth comparison of the differences in objective function values between different algorithms and the benchmark algorithm SQP, we calculated the deviation degree of the MPC and PINV algorithms relative to the final objective function value of the SQP algorithm. Specifically, by obtaining the final objective function value of each algorithm, we computed the percentage of the absolute value of the difference between this value and the final objective function value of the SQP algorithm relative to the final value of the SQP algorithm, so as to measure the degree of deviation.

The calculation results show that the deviation degree of MPC relative to SQP is 0.12%, and that of PINV relative to SQP is 0.22%. From this data, it can be inferred that in terms of the final value of the objective function, the result of the MPC algorithm is closer to that of the SQP algorithm, and its deviation degree is significantly smaller than that of the PINV algorithm.

4.4.2 Comparison of utilisation rates of different algorithms

From the Fig. 9, the average CPU utilisation rates of SQP, MPC and PINV are 16.53%, 23.51% and 11.00%, respectively. Taking the average CPU utilisation rate of SQP (16.53%) as the benchmark, the average CPU utilisation rate of MPC (23.51%) is significantly higher than that of SQP, with an increase of approximately 42.2%. In contrast, the average CPU utilisation rate of PINV (11.00%) is lower than that of SQP, with a decrease of approximately 33.4%.

Figure 9.

Comparison of utilisation rates of different algorithms.

In terms of memory usage, the average values of SQP, MPC and PINV are 9.02%, 20.15% and 13.75%, respectively. The average memory usage rate of MPC is approximately 2.23 times that of SQP. The average memory usage rate of PINV, which is 13.75%, is also higher than that of SQP, with an increase of approximately 52.4%.

In terms of energy consumption, the average values of SQP, MPC and PINV are 20.02W, 28.01W and 15.00W, respectively. The average energy consumption of MPC is approximately 1.4 times that of SQP. PINV has the lowest average energy consumption, which is 75% of that of SQP.

Overall, the MPC algorithm offers faster speed; it consumes significant resources, including higher CPU utilisation, memory usage and energy consumption. The PINV algorithm is relatively lightweight, with low and stable resource requirements and efficient energy utilisation, but it is inferior in terms of computational accuracy and speed. The SQP algorithm, on the other hand, has indicators that lie between the two, making it a more balanced choice. In practical applications, during the high-speed movement of tilt-rotor UAVs, the SQP algorithm can achieve a better balance between performance and resource consumption.