Dynamic Inversion-Enhanced U-Control of Quadrotor Trajectory Tracking

Abstract

This study proposes a universal (U)-control enhanced with dynamic inversion for a class of multiple-input multiple-output (MIMO) quadrotor flight tracking operations. In the technique, a robust compensator is integrated with the U-controller to achieve stabilisation at the equilibrium and setpoint tracking in the presence of any unmodelled uncertainties and external disturbances. A series of bench tests of simulated and real experiments on a Parrot Mambo quadrotor are conducted to show the design framework from the academic formulation to Simulink simulation and real flight tests.

1. Introduction

Over the last decade, the fast advancement of unmanned aerial vehicles (UAVs) has opened numerous possibilities in various fields such as surveillance, environmental monitoring, and delivery services. Among the many different types of UAVs, quadrotor systems have gained significant attention due to their mechanical simplicity, agile manoeuvrability, hovering capability, and their ability to vertically take off and land without requiring long runways. These factors along with the rapid development and mass production of electronics in general over the last decade have driven the prices of components down, making them much more affordable, and hence spurring a revolutionary impact in this field where high-performance drones are available to hobbyists whereas this was not the case before. Due to this wide-scale availability and adoption, the control of quadrotor systems has been a hot topic of research and development, which has brought substantial challenges, particularly due to their underactuated nature—using less inputs to drive more outputs. This underactuation requires sophisticated control strategies to achieve precise trajectory tracking and stabilisation which reduce the weight of the UAVs.

Regarding quadrotor control, classical control approaches, such as PID control, have been broadly taken as the first choice in applications due to their practical simplicity and ease of implementation [1,2]. Hybrid methods that combine the PID with other methods like Linear Quadratic Regulators (LQRs) [3] and Fuzzy control [4] have also been explored. However, these methods often fall short in handling the nonlinear dynamics and coupling effects inherent in quadrotor systems in flight tests. As a result, more advanced control techniques, including model predictive control (MPC) [5], sliding mode control (SMC) [6], and other nonlinear control approaches [7], have tried to enhance the capacity and robustness of quadrotor control systems.

Even after such tremendous effort has been devoted to UAV control, there are still some improvement domains that could enhance the control system design. Consequently, this paper expands a model-based U-control system design platform [8,9,10,11] to achieve trajectory tracking while maintaining robust stability in UAV control, which uses a transparent double dynamic inversion within two loops. In the inner loop, the plant dynamic inversion makes the control system design separate from the plant model, and the outer loop dynamic inversion makes the whole system performance specification be achieved with feedback and be unilaterally independent from the plant model. In comparison to some conventional delegate approaches, feedback linearisation control (FLC) has been widely considered as a representative approach [12] in dealing with nonlinear models by transforming into a linear one through coordinate conversions. Another well-known approach is backstepping control (BSC) [13], which is a Lyapunov-based recursive design approach and is able to accomplish the stabilisation and tracking tasks for a specific set of complex nonlinear systems. These two approaches use dynamic inversion in one way or another; however, they suffer from some shortcomings: (1) the FLC first utilises the complex coordinate transformation to design a linear control system, and then transforms the designed controller to the original controller to generate the required controller output. Unfavourably, such conversion needs strong linear algebraic conditions [14]. For BSC, it is normally applicable to a class of triangle-type models. The BSC is effective for second-order dynamics; however, the complexity would exponentially increase with the increase in the system dynamic order [15]. (2) To date, to the author’s best knowledge, almost all the existing nonlinear dynamic inversion formulations have assumed linearity with the plant input (that is, output affine plants). Accordingly, a general and concise mode-based dynamic inversion-augmented control method called model-based U-control has been established by using the expression of the U-model, which converts almost all the existing models into U-realisation without referring to the coordinate transform and backstepping process [9,10,16]. It should be mentioned that model-based U-control can dramatically reduce the complexity of designing linear, and particularly nonlinear, control systems, and the U-dynamic inversion and controller design are achieved separately. However, like all the other model-based dynamic inversion approaches, the U-dynamic inversion unavoidably depends on the model accuracy and the robustness against internal uncertainty and external disturbance. Developing an effective model-based U-control method has been the first motivation to balance the universal control flexibility and proper robustness. This study provides a new angle to enhance the objectives and to bench-test the new model-based U-control feasibility and effectiveness through simulations and real flight tests of the UAV.

The key contributions of this work to justify the research motivation and novelty include the following:

• This is the first validity bench test to model-based U-control with UAV flights in a comprehensive procedure, from analytical formulation to simulation demonstration and real flight tests. The applied control system design procedure has general guidance to other industrial applications and is a representative showcase to disseminate U-control methodology.
• In the methodology of U-model-based control, this study proposes a robust compensator to ensure trajectory accuracy and robustness against disturbances in the model-based dynamic inversion. This is supplementary to robust U-control approaches.
• We expand U-sliding mode control (USMC) to MIMO systems from single-input and single-output (SISO) formulation [17].
• We lay a practically effective foundation/platform for establishing the U-control strategy in coping with the challenges posed by underactuation and coupling in the control of quadrotor MIMO systems.
• A series of bench tests of simulated and real experiments on a Parrot Mambo quadrotor are conducted to show the design framework from academic formulation to Simulink simulation and real flight tests. This framework could be used for various UAV or related mobile robot motion control system design.

For the following presentation, Section 2 describes the quadcopter physical model as a problem background and introduces the U-model and U-model dynamic inversion (UMDI) control strategy along with a robust U-model dynamic inversion plus sliding mode control (UDI + USMC) compensator. Section 3 provides the computational simulation setup and results, validates the proposed approach with the expected performances. Finally, Section 4 summarises the study and gives some potential research directions.

2. Quadrotor Dynamics and U-Control

2.1. Quadcopter System Overview

Taken as a research background for testing the developed approach, this study chooses a Parrot^© Mambo Minidrone (Parrot S.A, Paris, France) (Figure 1), since it is cost-effective and has a manipulable small size for assessing the practical performance of the designed controllers in a closed-door environment without interference from outside factors like wind, etc. The Simulink Hardware Support Package for Parrot Minidrones, developed by the Original Equipment Manufacturer (OEM), delivers a realistic simulation environment for evaluating controller designs. After being successfully validated with the simulation demonstrations, the controller can then be embedded into the onboard processor by the Hardware Support Package provided to conduct experiments on the hardware to bench-test the outcomes in a real-world scenario.

The minidrone has two sensors on its underside: an ultrasound sensor and a low-resolution (640 × 480) camera. The ultrasound sensor measures vertical distance accurately up to 13 feet, though its performance varies depending on the floor material. The camera uses optical flow to estimate horizontal motion and speed, but its effectiveness depends on floor colour and texture. Thick carpets can interfere with both sensors. The minidrone’s affordability (around EUR 70) makes it appealing for research despite these limitations. Internally, the minidrone has a pressure sensor for Altitude and an Inertial Measurement Unit (IMU) that includes a three-axis accelerometer and gyroscope. The raw sensor data undergo bias correction, frame rotation, and low-pass filtering to remove noise. Two Kalman filters are implemented to estimate Altitude and XY position, while a complementary filter handles Roll, pitch and Yaw. The sampling time interval for the simulation and flight tests is set to 0.005 s. Actuator limitations are also taken into consideration, with controller output being capped with saturation functions to prevent exceeding the defined thrust limits. These processes are standardized by the OEM’s Hardware Support Package, ensuring consistent results. The ground effect also poses considerable challenges to this class of drones; due to their inherent low values of inertia and weight, they are easily affected by it.

The quadrotor system possesses six degrees of freedom to translate in the three Cartesian directions, i.e., along the X, Y and Z axes, and rotates about these three axes, with Roll $\varphi$, pitch $\theta$ and Yaw $\psi$ denoting the rotations about the X, Y and Z axes, respectively [18].

2.1.1. Assumptions

The dynamic model of the quadrotor system was derived under the assumption that both the structure and rotors are symmetrical and rigid. Another assumption is that the the quadrotor systems’ centre of mass is not situated at the origin of the body-fixed frame; however, the body–frame axes are aligned with the principal axes of the quadrotor. This assumption results in a diagonal inertia matrix. Lastly, it is assumed that the generated aerodynamic forces and moments are directly related to the square of the rotor’s speed.

2.1.2. Frames of Reference

The quadrotor’s dynamics are described by employing the following pair of reference frames: the inertial reference frame (IRF) $\mathbb{E}(X_E,Y_E,Z_E)$ and the Body Reference Frame (BRF) $\mathbb{B}(X_B,Y_B,Z_B)$, as illustrated in Figure 2. The quadrotor orientation in the coordinate space is described by $\mathrm{\Theta }={\left[\begin{array}{ccc}\varphi & \theta & \psi \end{array}\right]}^T$, where the Euler angle $\varphi$ represents the Roll (rotation about the X axis), $\theta$ represents the pitch (rotation about the Y axis) and $\psi$ represents the Yaw (rotation about the Z axis). These attitude angles are bounded as $\varphi \in (-{\displaystyle \frac{\pi }{2}},{\displaystyle \frac{\pi }{2}}),\theta \in (-{\displaystyle \frac{\pi }{2}},{\displaystyle \frac{\pi }{2}})$ and $\psi \in (-\pi ,\pi )$ because of the various acrobatic flying manoeuvres, which are not admissible. ${\omega }^B={\left[\begin{array}{ccc}p& q& r\end{array}\right]}^T$ represents the quadrotor’s angular velocity in the BRF.

To model the quadrotor UAV, a coordinate transformation matrix between these two reference frames is necessary as the translational equations of motion are described in the IRF while the rotational equations and angular velocities are defined in the BRF. The coordinate transformation matrix is therefore defined as follows $\left(R_{\varphi \theta \psi }\right)$ in Equation (1) [19].

$$\begin{array}{c}\hfill R_B^E=\left[\begin{array}{ccc}{C}_{\theta }{C}_{\psi }& S_{\theta }{C}_{\psi }{S}_{\varphi }-S_{\psi }{C}_{\varphi }& S_{\theta }{C}_{\psi }{C}_{\varphi }+S_{\psi }{S}_{\varphi }\\ C_{\theta }{S}_{\psi }& S_{\theta }{S}_{\psi }{S}_{\varphi }+C_{\psi }{C}_{\varphi }& S_{\theta }{S}_{\psi }{C}_{\varphi }-C_{\psi }{S}_{\varphi }\\ -S_{\theta }& C_{\theta }{S}_{\varphi }& C_{\theta }{C}_{\varphi }\end{array}\right]\end{array}$$

(1)

The transformational relationship between angular velocities calculated in the two different frames of reference, i.e., the IRF and BRF, can be described with the help of the vector $\mathbb{M}$, which is used to transform the angular velocities ${\omega }^B$ measured in the BRF to Euler rates $\dot{\mathrm{\Theta }}$ measured in the IRF and vice versa. Equations (2) and (3) show this relationship and how the transformations can be carried out from one frame of reference to another using the matrix $\mathbb{M}$ and its inverse ${\mathbb{M}}^{-1}$ as follows when ${\omega }^B=\mathbb{M}\dot{\mathrm{\Theta }}$:

$$\left[\begin{array}{c}p\\ q\\ r\end{array}\right]=\left[\begin{array}{ccc}1& 0& -S_{\theta }\\ 0& C_{\varphi }& S_{\varphi }{C}_{\theta }\\ 0& -S_{\varphi }& C_{\varphi }{C}_{\theta }\end{array}\right]\left[\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right]$$

(2)

and $\dot{\mathrm{\Theta }}={\mathbb{M}}^{-1}{\omega }^B$

$$\left[\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right]=\left[\begin{array}{ccc}1& S_{\varphi }{T}_{\theta }& C_{\varphi }{T}_{\theta }\\ 0& C_{\varphi }& -S_{\varphi }\\ 0& {\displaystyle \frac{S_{\varphi }}{C_{\theta }}}& {\displaystyle \frac{C_{\varphi }}{C_{\theta }}}\end{array}\right]\left[\begin{array}{c}p\\ q\\ r\end{array}\right]$$

(3)

where “S”, “C” and “T” denote “sin”, “cos” and “tan” trigonometric functions, respectively.

2.1.3. Rotor Dynamics

The thrust T and moment M generated by the $i^{th}$ rotor are calculated as follows:

$$\begin{array}{c}\hfill T_i=\left(C_T\rho A{R}^2\right){\mathrm{\Omega }}_i^2=\left(k_{Ti}\right){\mathrm{\Omega }}_i^2\\ \hfill M_i=\left(C_Q\rho A{R}^3\right){\mathrm{\Omega }}_i^2=\left(k_{Mi}\right){\mathrm{\Omega }}_i^2\end{array}$$

(4)

where ${\mathrm{\Omega }}_i$ represents the $i^{th}$ rotor speed. $k_{Ti}$ and $k_{Mi}$ are aerodynamic constants that are associated with the thrust and torque generated by the rotor depending on its shape and configuration. These coefficients can be determined experimentally (

Table 1. Parrot Mambo parameters.

Parameter	Value	Units
Thrust Coefficient, $C_T$	$0.0107$
Torque Coefficient, $C_Q$	$7.8264\times {10}^{-4}$
Air density, $\rho$	$1.184$	$\mathrm{kg}/{\mathrm{m}}^3$
Rotor Radius, R	$0.033$	m
Rotor Disk Area, $A=\pi R^2$	$0.0034212$	${\mathrm{m}}^2$
L	$0.0624$	m
Rotor inertia, $J_r$	$1.0209375\times {10}^{-7}$	$\mathrm{kg}.{\mathrm{m}}^2$
Inertia moment along X axis, $J_{xx}$	$0.0000582857$	$\mathrm{kg}.{\mathrm{m}}^2$
Inertia moment along Y axis, $J_{yy}$	$0.0000716914$	$\mathrm{kg}.{\mathrm{m}}^2$
Inertia moment along Z axis, $J_{zz}$	$0.0001$	$\mathrm{kg}.{\mathrm{m}}^2$

Values taken from Parrot Minidrone Simulink hardware package.

).

The thrust and moment generated by the rotor rely on the rotor’s physical properties, the air density, and the speed at which the rotor is spinning. The term ${\mathrm{\Omega }}^2$ indicates that the thrust and moment are in proportion to the square of the angular velocity of the rotor, which is a key aspect in understanding how changes in rotor speed affect the generated thrust and moment.

The total moments are formulated below:

$$\begin{array}{cc}\hfill {\tau }_x& ={\displaystyle \frac{\sqrt{2}}{2}}L(T_1-T_2-T_3+T_4)\hfill \\ \hfill {\tau }_y& ={\displaystyle \frac{\sqrt{2}}{2}}L(-T_1-T_2+T_3+T_4)\hfill \\ \hfill {\tau }_z& =(M_1-M_2+M_3-M_4)\hfill \end{array}$$

(5)

where L denotes the distance from the quadrotor mass centre to the rotor rotation axis. The control inputs $U_i$ where $i\in 1,\dots ,4$ are described as follows:

$$\left[\begin{array}{c}{U}_1\\ U_2\\ U_3\\ U_4\end{array}\right]=\left[\begin{array}{cccc}1& 1& 1& 1\\ {\displaystyle \frac{\sqrt{2}}{2}}L& -{\displaystyle \frac{\sqrt{2}}{2}}L& -{\displaystyle \frac{\sqrt{2}}{2}}L& {\displaystyle \frac{\sqrt{2}}{2}}L\\ -{\displaystyle \frac{\sqrt{2}}{2}}L& -{\displaystyle \frac{\sqrt{2}}{2}}L& {\displaystyle \frac{\sqrt{2}}{2}}L& {\displaystyle \frac{\sqrt{2}}{2}}L\\ {\displaystyle \frac{k_{M1}}{K_{T1}}}& {\displaystyle \frac{k_{M2}}{K_{T2}}}& {\displaystyle \frac{k_{M3}}{K_{T3}}}& {\displaystyle \frac{k_{M4}}{K_{T4}}}\end{array}\right]\left[\begin{array}{c}{T}_1\\ T_2\\ T_3\\ T_4\end{array}\right]$$

(6)

where $U_1$ is for the whole thrust added on the quadrotor, $U_2$ represents the Roll torque ${\tau }_x$, $U_3$ represents the pitch torque ${\tau }_y$ and $U_4$ represents the Yaw torque ${\tau }_z$.

2.1.4. Translational Subsystem

The equations of motion that describe the system behaviour along the three Cartesian axes in the inertial reference frame (IRF) are

$$m\ddot{X}=mg+R_{\varphi \theta \psi }^T{u}_1{Z}_E$$

(7)

where m is the quadrotor mass, g is the acceleration due to the gravity of Earth, $R_{\varphi \theta \psi }^T$ is the coordinate transformation matrix from Equation (1) and $z_E={\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^T$ is a unit vector in the inertial frame of reference (IRF).

Equation (7) is then expanded and rewritten as follows:

$$\left[\begin{array}{c}\ddot{x}\\ \ddot{y}\\ \ddot{z}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ -g\end{array}\right]+\left[\begin{array}{c}{C}_{\varphi }{S}_{\theta }{C}_{\psi }+S_{\varphi }{S}_{\psi }\\ C_{\varphi }{S}_{\theta }{S}_{\psi }-S_{\varphi }{C}_{\psi }\\ C_{\varphi }{C}_{\theta }\end{array}\right]{\displaystyle \frac{U_1}{m}}$$

(8)

2.1.5. Rotational Subsystem

The rotational dynamics of a quadrotor are derived by utilizing the Newton–Euler equations as follows:

$$J{\dot{\omega }}^B+{\omega }^B\times J{\omega }^B+M_g=\tau$$

(9)

where J is the quadrotor inertia matrix, $M_g$ is the gyroscopic moment generated by the rotors and $\tau ={\left[\begin{array}{ccc}{\tau }_x& {\tau }_y& {\tau }_z\end{array}\right]}^T$ represents the cumulative effect of moments exerted on the quadrotor. The gyroscopic moment can be characterized as

$$M_g={\omega }^B\times {\left[\begin{array}{ccc}0& 0& J_r{\omega }_r\end{array}\right]}^T$$

(10)

where $J_r$ denotes the propeller’s inertia and ${\omega }_r$ represents the relative speed of the propeller, which is defined as

$${\mathrm{\Omega }}_r=\sum _{i=1}^4{(-1)}^{i+1}{\omega }_i$$

(11)

Consequently, the attitude dynamics of the system can be summarised by using Equations (6) and (9) as follows:

$$\left[\begin{array}{c}\ddot{\varphi }\\ \ddot{\theta }\\ \ddot{\psi }\end{array}\right]=\left[\begin{array}{c}\left({\displaystyle \frac{J_{yy}-J_{zz}}{J_{xx}}}\right)\dot{\theta }\dot{\psi }\\ \left({\displaystyle \frac{J_{zz}-J_{xx}}{J_{yy}}}\right)\dot{\varphi }\dot{\psi }\\ \left({\displaystyle \frac{J_{xx}-J_{yy}}{J_{zz}}}\right)\dot{\varphi }\dot{\theta }\end{array}\right]-J_r{\mathrm{\Omega }}_r\left[\begin{array}{c}{\displaystyle \frac{\dot{\theta }}{J_{xx}}}\\ {\displaystyle \frac{-\dot{\varphi }}{J_{yy}}}\\ 0\end{array}\right]+\left[\begin{array}{c}{\displaystyle \frac{U_2}{J_{xx}}}\\ {\displaystyle \frac{U_3}{J_{yy}}}\\ {\displaystyle \frac{U_4}{J_{zz}}}\end{array}\right]$$

(12)

Subsystems (8) and (12) together represent the full dynamic description of a quadrotor system. As can be observed, the equations that describe the rotational motion of the quadrotor are completely independent of the translational elements X, Y and Z and their derivatives. However, the quadrotor translational motion is dependent on the rotational components $\varphi ,\theta$ and $\psi$. Therefore, the overall system is made up of two subsystems with the rotational system in the inner loop, while the translational subsystem in the outer loop generates the required reference signals for the inner loop in the form of desired Roll and pitch angles to track the desired location. Effective attitude control, therefore, is critical to ensure that the UAVs can perform a wide range of tasks that can range from simple hovering to complex aerial manoeuvrers while maintaining stability and responding appropriately to external disturbances like wind since the UAV’s position is controlled indirectly by manipulating the pitch, Roll and Yaw.

2.2. Model-Based U-Control Platform

Figure 3 shows the model-based U-control framework [8,20]; the controller consists of two parts, with the first being the universal controller ($G_{cl}$) which is to be designed completely independently from the plant, according to the desired indices, and the second part being the plant inverse ($G_p^{-1}$). Consequently, the primary objective in the entire designing process boils down to the determination of the inverse $G_p^{-1}$ of the plant $G_p$ that is to be regulated and the combination of the plant inverse with the specified universal controller $G_{cl}$.

The following is a brief overview highlighting the up-to-date research of the U-control design framework approach. In [20], the authors propose a pole placement control strategy for complex dynamic systems by leveraging U-models, demonstrating improved system stability and performance. In [8], a general U-block model-based design procedure is proposed that aims to simplify the control of a class of complex polynomial difference equations. The methodology effectively transforms nonlinear control problems into linear ones using U-models, facilitating easier controller design and implementation. Further developing on the previous work, [9] extends the U-model methodology to rational nonlinear systems, offering insights into controlling more complex dynamics, highlighting the flexibility of U-models in addressing the challenges posed by rational nonlinearities in control systems. In [21], the authors address the issue of input delay in nonlinear processes, using a U-model-based predictive control approach that significantly enhances control accuracy and system responsiveness, showcasing the practical applicability of U-models in real-world scenarios. More recently, [22] explores the application of U-model-augmented control in the context of underactuated unmanned marine robotics. Ref. [10] presents solutions for reversing dynamics based on the U-model and [11] outlines a comprehensive U-control framework for continuous-time systems. Ref. [16] consolidates the U-model and U-control methodologies, providing a detailed analysis of their application in complex dynamic systems.

Remark 1. 

It is worth noting that robust flat control (RFC) [23,24] has been a very effective approach using dynamic inversion to enhance the control system design in a large class of wheeled mobile robots, UAVs and other motion systems. In comparison with model-based U-control, RFC needs coordinate transform to achieve the dynamic inversion to determine the control input and the model-based U-control directly takes dynamic inversion from the plant U-model realisation, which is straightforward in the reformulation of the plant model into a control-oriented time-varying model [16]. Regarding the robustness of the two approaches, RFC could be more effective than model-based U-control. Surely, RFC provides a very useful reference to improve U-control robustness and we will provide simulation comparison as well in future publications.

2.2.1. Design of U-Controller $G_{cl}$

To begin the process of designing the universal controller, assign the desired system performance in advance with the selection of an ideal closed-loop transfer function as

$${\displaystyle \frac{Y\left(s\right)}{R\left(s\right)}}=G\left(s\right)={\displaystyle \frac{{\omega }_n^2}{s^2+2\zeta {\omega }_{n}s+{\omega }_n^2}}={\displaystyle \frac{4}{s^2+4s+4}}$$

(13)

where $\zeta =1$ is the selected damping ratio and ${\omega }_n$ = 2 is the chosen undamped natural frequency, so that the system output response is critically damped to prevent any overshoot, ensuring that the system returns to the desired state in the shortest possible time without oscillating [25,26].

The determination of the universal controller treats the pair of the plant and its inverse as a unit in the feedback control system, and then takes the inverse of the closed-loop transfer function in Equation (13) as follows:

$$G_{cl}={\displaystyle \frac{G}{1-G}}={\displaystyle \frac{4}{s^2+4s}}$$

(14)

2.2.2. U-Polynomial Realisation

Before calculating the plant inverse, the system needs to be converted to its U-realisation (either U-Polynomial or U-state space) depending on the system. A brief overview of the realisation of the system to its appropriate U-realisation is given in this section before proceeding to controller design.

Consider a general SISO continuous-time polynomial system defined with $\left(y\right(t),\lambda (t)$ and $u\left(t\right))$ where $y\left(t\right)\in \mathbb{R}$ and $u\left(t\right)\in \mathbb{R}$ are the system output and system input, respectively, at time $t\in {\mathbb{R}}^{+}$, which is described as follows:

$$\mathrm{{\rm Y}}{\left(t\right)}^{\left(M\right)}=\sum _{j=0}^J{\lambda }_j({\mathrm{{\rm Y}}}_{M-1},({\mathrm{\Xi }}_{N-1},\mu ){\left(u{\left(t\right)}^{\left(N\right)}\right)}^j,M\ge N$$

(15)

where the term ${y\left(t\right)}^{\left(M\right)}$ represents the $M^{th}$-order derivative of the output $y\left(t\right)$ and the term ${u\left(t\right)}^{\left(N\right)}$ represents the $N^{th}$-order derivative of the system input $u\left(t\right)$, while the term ${\lambda }_j(\ast )\in \mathbb{R}$ represents a time-varying parameter that absorbs all the other inputs and outputs.

$$\begin{array}{cc}\hfill {\mathrm{{\rm Y}}}_{M-1}& =\left[y{\left(t\right)}^{(M-1)},y{\left(t\right)}^{(M-2)},\dots y\left(t\right)\right]\in {\mathbb{R}}^M\hfill \\ \hfill {\mathrm{\Xi }}_{N-1}& =\left[u{\left(t\right)}^{(N-1)},u{\left(t\right)}^{(N-2)},\dots u\left(t\right)\right]\in {\mathbb{R}}^M\hfill \\ \hfill \mu & ={\left(u{\left(t\right)}^{\left(N\right)}\right)}^j\hfill \end{array}$$

(16)

Consider an example showcasing the U-realisation of a classical nonlinear polynomial representation [27]:

$$\begin{array}{c}\hfill \ddot{z}=z+(z-2)\dot{z}+(4+z^2)v+(z+{\dot{z}}^2)v^2\end{array}$$

(17)

The U-model realisation, then, is

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \ddot{z}& ={\lambda }_0+{\lambda }_{1}v+{\lambda }_2{v}^2\hfill \\ \hfill {\lambda }_0& =z+(z-2)\dot{z}\hfill \\ \hfill {\lambda }_1& =(4+z^2)\hfill \\ \hfill {\lambda }_2& =(z+{\dot{z}}^2)\hfill \end{array}\hfill \end{array}\right.$$

(18)

Comparing Equations (17) and (18), the model’s U-realisation is quite forthright. The U-realisation is the same as the original polynomial in attributes but is expressed in a control-oriented structure.

2.2.3. U-State Space Realisation

Consider a general SISO state space description as follows:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \dot{Q}& =G(Q,u)\hfill \\ \hfill y& =H\left(Q\right)\hfill \end{array}\hfill \end{array}\right.$$

(19)

where $u,y,\in \mathbb{R},Q\in {\mathbb{R}}^n,$$G\in {\mathbb{R}}^n$→ and $H\in \mathbb{R}$ are the smooth mappings $u\to Q\in \mathbb{R}\to {\mathbb{R}}^n$ and $Q\to y\in {\mathbb{R}}^n\to \mathbb{R}$, respectively.

This state space representation of the system is expressed as a multilayer polynomial representation as follows:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill {\dot{q}}_1& =G_1(q_1,q_2,\dots ,q_n),\hfill \\ \hfill {\dot{q}}_2& =G_2(q_1,q_2,\dots ,q_n),\hfill \\ \hfill ⋮\\ \hfill {\dot{q}}_n& =G_n(q_1,q_2,\dots ,q_n,u),\hfill \\ \hfill y& =H(q_1,q_2,\dots ,q_n)\hfill \end{array}\hfill \end{array}\right.$$

(20)

Convert this multilayer polynomial expression into a multilayer U-realisation as

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill {\dot{q}}_1& =\sum _{i=0}^n{\lambda }_{1i}{g}_{1i}\left(q_2\right),\hfill \\ \hfill {\dot{q}}_2& =\sum _{i=0}^n{\lambda }_{2i}{g}_{2i}\left(q_3\right),\hfill \\ \hfill ⋮\\ \hfill {\dot{q}}_n& =\sum _{i=0}^n{\lambda }_{ni}{g}_{ni}\left(u\right),\hfill \\ \hfill y& =\sum _{i=0}^n{h}_i(q_1,q_2,\dots ,q_n)\hfill \end{array}\hfill \end{array}\right.$$

(21)

where $\lambda (\ast )$ and $g(\ast )$ are the time-varying parameters absorbing the remaining terms in the state space model. It is worth mentioning that the U-state space realisation is a form of expanded multilayer U-Polynomial expression because each line of the U-state space realisation can be considered a U-Polynomial by itself.

In case Equation (21) is unable to exhibit the explicit association between the system output term y and the input term u, differentiate y and substitute to $x_i$ items in system (21) with the polynomial expression for each $x_i$ in Equation (20). This iterative process should be carried out until an explicit expression between $y^{\left(m\right)}$ and u is exhibited, where $y^{\left(m\right)}$ denotes the $m^{th}$ derivative of the model output y.

2.3. Determining the Dynamic Inverse

After establishing a direct association between $y^{\left(m\right)}$ and u in system (21), the system equation can be restructured as follows:

$$y^{\left(m\right)}={\lambda }_1(q_1,q_2,\dots ,q_{n-1},q_n)+{\lambda }_2(q_1,q_2,\dots ,q_{n-1},q_n)u$$

(22)

where ${\lambda }_i$ ($i\in 1,2,\dots$) is a time-varying gain vector soaking up all the state variables. Equation (22) can be identified as an alternative U-model polynomial described by Equation (15). The output of the plant inverter $u_i$ can be calculated by reorganising and isolating the input term u in Equation (22):

$$u={\displaystyle \frac{y^{\left(m\right)}-{\lambda }_1(q_1,q_2,\dots ,q_{n-1},q_n)}{{\lambda }_2(q_1,q_2,\dots ,q_{n-1},q_n)}}$$

(23)

In addition, the state variables $q_i$ are worked out by the integral of ${\dot{q}}_i$ associated with the plant. This relationship is expressed as $q_i={\dot{q}}_i·{\displaystyle \frac{1}{s}}$, where ${\displaystyle \frac{1}{s}}$ is the integral operation of the Laplace transform.

For illustration purposes, consider a simple nonlinear plant modelled by

$$\left\{\begin{array}{cc}{\dot{q}}_1\hfill & =q_2\hfill \\ {\dot{q}}_2\hfill & =-q_1-sin\left(q_2\right)+u\hfill \\ y\hfill & =q_1\hfill \end{array}\right.$$

(24)

Since no direct relation between the input u and output y can be observed here, y is differentiated repeatedly till an expression linking the two explicitly is obtained.

$$\begin{array}{cc}\hfill \ddot{y}& =-q_1-sin\left(q_2\right)+u\hfill \end{array}$$

(25)

Equation (25) establishes a direct link between the input and output, making it a viable option for dynamic inversion.

The associated U-realisation and the dynamic inverse is then

$$\left\{\begin{array}{cc}\ddot{y}\hfill & ={\lambda }_0+{\lambda }_{1}u\hfill \\ {\lambda }_0\hfill & =-q_1-sin\left(q_2\right);{\lambda }_1=1\hfill \\ u_i\hfill & ={\displaystyle \frac{\ddot{y}-{\lambda }_0}{{\lambda }_1}}=\ddot{y}+q_1+sin\left(q_2\right)\hfill \end{array}\right.$$

(26)

With the aim of assigning the plant output y equivalent to the output (v) of the controller $\left(G_{cl}\right)$, the plant dynamic inversion $G_p^{-1}$ can be computed. This involves substituting the highest order of the y term with v of the same order. Following this, the expression is integrated twice to yield y:

$$y={\displaystyle \frac{{\lambda }_0}{s^2}}+{\displaystyle \frac{{\lambda }_1}{s^2}}u$$

(27)

To directly determine the second-order derivative $\ddot{v}$, the universal controller can be modified directly to ensure that the same dynamic orders of v and y are in the designed control system.

Remark 2. 

For the system represented by Equation (24), the method to accomplish this goal is through multiplication of the second-order Laplace operator $s^2$ with the universal controller $G_{cl}$ to yield $\ddot{v}=s^{2}v$. Consequently, the output y is substituted by the second-order derivative of the invariant controller $G_{cl}$ output, denoted as $\ddot{v}$. Therefore, $u=q_1+sin{q}_2+\ddot{v}$ results in $G_p{G}_p^{-1}=1.$ Furthermore, the universal controller $G_{cl}$ from Equation (14) can be rewritten as follows:

$$G_{cl}={\displaystyle \frac{4s^2}{s^2+4s}}$$

(28)

2.4. Compensated Plant Inverter $G_c^{-1}$ Design

As shown in Figure 4, expanded from Figure 3, this study introduces a sliding mode control (SMC) into model-based U-control for MIMO systems to enhance their robustness through the addition of a robust SMC to accommodate errors induced from the plant inverter $G_p^{-1}$ due to the plant uncertainties and system disturbances. To determine $G_c^{-1}$, for a general second-order system that is defined as

$$\left\{\begin{array}{cc}{\dot{q}}_1\hfill & =g_1(q_1,q_2)\hfill \\ {\dot{q}}_2\hfill & =g_2(q_1,q_2)+h(q_1,q_2)u+d\hfill \\ y\hfill & =q_1\hfill \end{array}\right.$$

(29)

where $g_i\in {\mathbb{R}}^n$ is a smooth function absorbing all states, h represents the relationship between state variables and control outputs u and d are the disturbance of where $\left|d\right|<\overline{d}$.

The state error vector is defined as follows:

$$\left\{\begin{array}{c}{\tilde{q}}_1=q_1-q_d\hfill \\ {\tilde{q}}_2={\dot{q}}_1-{\dot{q}}_d\hfill \end{array}\right.$$

(30)

where ${\dot{q}}_1$ and ${\dot{q}}_d$ are the derivatives of the plant output $q_1$ and desired output $q_d$, respectively. Then, we can define a sliding mode manifold $\sigma$ by [28]

$$\sigma =c_1{\tilde{q}}_1+{\tilde{q}}_2$$

(31)

And its derivative $\dot{\sigma }$ is calculated as follows:

$$\dot{\sigma }=c_1\left({\dot{q}}_1-{\dot{q}}_d\right)+\left({\dot{q}}_2-{\ddot{q}}_d\right)$$

(32)

where $c_i\in {\mathbb{R}}^{+}andi\in {\mathbb{R}}^{+}$ are the coefficient specifying a stable Hurwitz polynomial and the asymptotically convergent speed [28].

Similar to the state space expression, replace the highest-order derivative of the desired output ${\ddot{q}}_d$ with $vs.$ Consequently, the lower-order derivatives of the desired output $q_d$ are

$$\begin{array}{cc}\hfill {\dot{q}}_d& =vs.\left({\displaystyle \frac{1}{s}}\right)\hfill \\ \hfill q_d& =vs.{\left({\displaystyle \frac{1}{s}}\right)}^2\hfill \end{array}$$

(33)

where the term ${\displaystyle \frac{1}{s}}$ is the Laplace integration operator. Accordingly, convert error Equations (30) and (32) into

$$\left\{\begin{array}{c}{\tilde{q}}_1=q_1-v{\left({\displaystyle \frac{1}{s}}\right)}^2\hfill \\ {\tilde{q}}_2={\dot{q}}_1-v\left({\displaystyle \frac{1}{s}}\right)\hfill \end{array}\right.$$

(34)

$$\dot{\sigma }=c_1\left({\dot{q}}_1-v\left({\displaystyle \frac{1}{s}}\right)\right)+\left({\ddot{q}}_1-v\right)$$

(35)

The calculation process for the $m^{th}$-order derivative of state $q_1$, i.e., $\stackrel{m}{q_1}$, is found by continuously differentiating its lower-order derivatives; the detailed process is described in [17]. In the case of a second-order system, as described here, m = 2, so ${\ddot{q}}_1$ is calculated as

$$\left\{\begin{array}{c}{\dot{q}}_1=G_1(q_1,q_2)={\mathcal{G}}_1(q_1,q_2)\hfill \\ {\ddot{q}}_1={\displaystyle \frac{\partial {\mathcal{G}}_1}{\partial q_1}}{G}_1+{\displaystyle \frac{\partial {\mathcal{G}}_1}{\partial q_2}}(G_2+h{u}_s+d)\hfill \end{array}\right.$$

(36)

Substituting this value of ${\ddot{q}}_1$ in Equation (35), it is then rewritten as follows:

$$\dot{\sigma }=c_1\left({\dot{q}}_1-v\left({\displaystyle \frac{1}{s}}\right)\right)+\left(\left({\displaystyle \frac{\partial {\mathcal{G}}_1}{\partial q_1}}{G}_1+{\displaystyle \frac{\partial {\mathcal{G}}_1}{\partial q_2}}(G_2+hu+d)\right)-v\right)$$

(37)

The robust sliding mode controller output is made of two parts, with the first being the equivalent control $u_{eq}$ and the second part being the switching control $u_{sw}$. To obtain the equivalent controller output, let $\dot{\sigma }=0$ [28].

Similarly, the USMC final controller output consists of two parts, with the first being from the ideal plant inverse $G_p^{-1}$, i.e., ideal inversion output $u_i$ and the compensated inversion output $u_c$. Equating them and reorganising the terms gives the following:

$$\begin{array}{cc}\hfill u_i+u_c& =u_{eq}+u_{sw}\hfill \\ \hfill u_c& =u_{sw}+u_{eq}-u_i\hfill \end{array}$$

(38)

Differentiating the system output y from Equation (29) until you obtain a direct relationship between y and u shows

$$\begin{array}{cc}\hfill \dot{y}& ={\dot{q}}_1=q_2\hfill \\ \hfill \ddot{y}& ={\ddot{q}}_1={\dot{q}}_2=g(q_1,q_2)+u={\lambda }_0+{\lambda }_{1}u\hfill \end{array}$$

(39)

In Equation (39), the term ${\dot{q}}_2$ is rewritten into its U-model equivalent form, i.e., ${\lambda }_0+{\lambda }_{1}u$. Referring to the relationship between the plant input u and the system output y, the ideal inverted output, denoted as $u_i$, can be determined as follows:

$$u_i={\displaystyle \frac{\ddot{y}-{\lambda }_0}{{\lambda }_1}}={\displaystyle \frac{v-{\lambda }_0}{{\lambda }_1}}$$

(40)

Replacing the term ${\dot{q}}_2$ with its U-model equivalent from Equation (35) gives

$$\dot{\sigma }=c\left({\dot{q}}_1-v{\displaystyle \frac{1}{s}}\right)+({\lambda }_0+{\lambda }_{1}u-v)$$

(41)

Setting $\dot{\sigma }=0$ to calculate $u_{eq}$ gives

$$u_{eq}={\displaystyle \frac{v-c\left({\dot{q}}_1-v\frac{1}{s}\right)-{\lambda }_0}{{\lambda }_1}}$$

(42)

$$u_{sw}=-(\epsilon sign\left(\sigma \right)+k\sigma )$$

(43)

where $\epsilon =(\overline{d}+\rho )$ with $\rho >0$ and $k>0$. Then, from Equation (38),

$$\begin{array}{cc}\hfill u_c& =u_{sw}+u_{eq}-u_i\hfill \\ \hfill & =-(\epsilon sign\left(\sigma \right)+k\sigma )+{\displaystyle \frac{v-c\left(q_2-v\frac{1}{s}\right)-{\lambda }_0}{{\lambda }_1}}-{\displaystyle \frac{(v-{\lambda }_0)}{{\lambda }_1}}\hfill \\ \hfill & =-(\epsilon sign\left(\sigma \right)+k\sigma )+{\displaystyle \frac{\cancel{v}-c\left(q_2-v\frac{1}{s}\right)-\cancel{{\lambda }_0}-\cancel{v}+\cancel{{\lambda }_0}}{{\lambda }_1}}\hfill \\ \hfill & =-(\epsilon sign\left(\sigma \right)+k\sigma )+{\displaystyle \frac{-c\left(q_2-v\frac{1}{s}\right)}{{\lambda }_1}}\hfill \end{array}$$

(44)

where $c\in {\mathbb{R}}^{+},k\in {\mathbb{R}}^{+}$ and $\epsilon =\in {\mathbb{R}}^{+}$ are design parameters associated with the Sliding Mode Compensator. This $u_c$ represents the compensator output and together with the ideal inverter output $u_i$ forms the final overall controller output $u_s$, as can be shown in Figure 4.

Remark 3. 

To reduce the effects of the chattering effect, replace the discontinuous function $sign\left(\sigma \right)$ with a smooth sigmoid function [28] as $sign\left(\sigma \right)\approx {\displaystyle \frac{\sigma }{\left|\sigma \right|+\xi }}$, where ξ is a small positive scalar and it can be observed that, point-wise,

$$\underset{\xi \to 0}{lim}{\displaystyle \frac{\sigma }{\left|\sigma \right|+\xi }}=sign\left(\sigma \right)$$

Equation (44) can therefore be rewritten by substituting the $sign(\ast )$ function as follows:

$$u_c=\left(-(\epsilon {\displaystyle \frac{\sigma }{\left|\sigma \right|+\xi }}+k\sigma )+{\displaystyle \frac{-c\left(q_2-v\frac{1}{s}\right)}{{\lambda }_1}}\right)$$

(45)

Taking the final USMC controller output, $u_s=u_c+u_i$, and substituting in the value of $u_s$ in Equation (37), we end up with the following:

$$\dot{\sigma }=-(\epsilon sign\left(\sigma \right)+k\sigma )+{\displaystyle \frac{\partial {\mathcal{G}}_1}{\partial q_2}}d$$

(46)

To analyse the stability of this approach, we use Lyapunov theory, where a typical Lyapunov candidate function and its time derivative are defined as

$$\begin{array}{cc}\hfill V\left(\sigma \right)& ={\displaystyle \frac{1}{2}}{\sigma }^2\ge 0\hfill \\ \hfill \dot{V}\left(\sigma \right)& =\sigma \dot{\sigma }\hfill \end{array}$$

(47)

using control input u that stabilizes the system satisfying $\dot{\sigma }=-(\epsilon sign\left(\sigma \right)+k\sigma )+{\displaystyle \frac{\partial {\mathcal{G}}_1}{\partial q_2}}d$, where k is a positive gain and $\epsilon =(\overline{d}+\rho )$, with $\rho >0$. Then, substituting the control law into the time derivative,

$$\begin{array}{cc}\hfill \dot{V}\left(\sigma \right)=\sigma \dot{\sigma }& =\sigma \left(-\epsilon sign\left(\sigma \right)-k\sigma +{\displaystyle \frac{\partial {\mathcal{G}}_1}{\partial q_2}}d\right)\hfill \\ \hfill & =-\epsilon \left|\sigma \right|-k{\sigma }^2+\sigma \left({\displaystyle \frac{\partial {\mathcal{G}}_1}{\partial q_2}}d\right)\hfill \\ \hfill & =-(\overline{d}+\rho )\left|\sigma \right|-k{\sigma }^2+\sigma \left({\displaystyle \frac{\partial {\mathcal{G}}_1}{\partial q_2}}d\right)\le -\rho \left|\sigma \right|-k{\sigma }^2\le 0\hfill \end{array}$$

(48)

and by selecting gains where $k>0$ and $\epsilon =(\overline{d}+\rho )$, where $\rho$ is a positive number that is greater than $\left|d\right|$ to cancel out the effect of disturbances, it follows that

$$\dot{V}\left(\sigma \right)\le 0$$

(49)

This demonstrates that when appropriate positive gains are used in the control law, $\dot{V}\left(\sigma \right)$ is negative semi-definite, leading to stability.

Secondly, the controller $G_{cl}\left(s\right)$ is used to specify the desired system performance in Equation (13), where $\zeta =1$ represents the selected damping ratio and ${\omega }_n=2$ is the undamped natural frequency so that the closed loop control system generates a fast monotonic response and zero steady state error to level reference. The corresponding characteristic equation $D\left(s\right)$ and its roots $s_1$ and $s_2$ are

$$\begin{array}{cc}\hfill D\left(s\right)& =s^2+4s+4=0\hfill \\ \hfill s_1,s_2& =-2\hfill \end{array}$$

(50)

Since the first column of the Routh array is $[1,4,4]$ and all its elements are positive, indicating no change of sign, and all the poles are in the left half of the s-plane, the system as a whole is Routh-stable while the inner loop dynamic inversion results in $G_p{G}_p^{-1}\to 1$.

3. Main Results

3.1. Model-Based U-Control of the Rotational Subsystem

The Roll of the quadrotor system is represented by the following subsystem:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =\dot{\varphi }=x_2\hfill \\ \hfill {\dot{x}}_2& =\ddot{\varphi }={\displaystyle \frac{J_{yy}-J_{zz}}{J_{xx}}}\dot{\theta }\dot{\psi }-J_r{\mathrm{\Omega }}_r{\displaystyle \frac{\dot{\theta }}{J_{xx}}}+{\displaystyle \frac{U_2}{J_{xx}}}\hfill \end{array}$$

(51)

The U-realisation of the Roll subsystem is

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \ddot{\varphi }& ={\lambda }_0+{\lambda }_1{U}_2\hfill \\ \hfill {\lambda }_0& ={\displaystyle \frac{J_{yy}-J_{zz}}{J_{xx}}}\dot{\theta }\dot{\psi }-J_r{\mathrm{\Omega }}_r{\displaystyle \frac{\dot{\theta }}{J_{xx}}}\hfill \\ \hfill {\lambda }_1& ={\displaystyle \frac{1}{J_{xx}}}\hfill \end{array}\hfill \end{array}\right.$$

(52)

The control input term $U_2$ is determined by

$$U_2={\displaystyle \frac{\ddot{\varphi }-{\lambda }_0}{{\lambda }_1}}$$

(53)

Similarly, for the pitch channel, the U-controller is designed by first converting the state space subsystem into its U-control expression and then extracting the control input term as above.

$$\begin{array}{cc}\hfill {\dot{x}}_3& =\dot{\theta }=x_4\hfill \\ \hfill {\dot{x}}_4& =\ddot{\theta }={\displaystyle \frac{J_{zz}-J_{xx}}{J_{yy}}}\dot{\varphi }\dot{\psi }-J_r{\mathrm{\Omega }}_r{\displaystyle \frac{-\dot{\varphi }}{J_{yy}}}+{\displaystyle \frac{U_3}{J_{yy}}}\hfill \end{array}$$

(54)

The U-model realisation of the pitch subsystem is

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \ddot{\theta }& ={\lambda }_0+{\lambda }_1{U}_3\hfill \\ \hfill {\lambda }_0& ={\displaystyle \frac{J_{zz}-J_{xx}}{J_{yy}}}\dot{\varphi }\dot{\psi }-J_r{\mathrm{\Omega }}_r{\displaystyle \frac{\dot{\varphi }}{J_{yy}}}\hfill \\ \hfill {\lambda }_1& ={\displaystyle \frac{1}{J_{yy}}}\hfill \end{array}\hfill \end{array}\right.$$

(55)

We can work out the control input term $U_3$ using

$$U_3={\displaystyle \frac{\ddot{\theta }-{\lambda }_0}{{\lambda }_1}}$$

(56)

Again, repeating the same process as before by converting the state space subsystem into its U-control realisation, the Yaw subsystem can therefore be written as follows:

$$\begin{array}{cc}\hfill {\dot{x}}_5& =\dot{\psi }=x_6\hfill \\ \hfill {\dot{x}}_6& =\ddot{\psi }={\displaystyle \frac{J_{xx}-J_{yy}}{J_{zz}}}\dot{\varphi }\dot{\theta }+{\displaystyle \frac{U_4}{J_{zz}}}\hfill \end{array}$$

(57)

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \ddot{\psi }& ={\lambda }_0+{\lambda }_1{U}_4\hfill \\ \hfill {\lambda }_0& ={\displaystyle \frac{J_{xx}-J_{yy}}{J_{zz}}}\dot{\varphi }\dot{\theta }\hfill \\ \hfill {\lambda }_1& ={\displaystyle \frac{1}{J_{zz}}}\hfill \end{array}\hfill \end{array}\right.$$

(58)

The control input term $U_4$ can now be calculated as

$$U_4={\displaystyle \frac{\ddot{\psi }-{\lambda }_0}{{\lambda }_1}}$$

(59)

3.2. Model-Based U-Control of Translational Subsystem

For the Altitude channel, since there is a direct control input $U_1$ for its regulation, we similarly designed the channel to follow that of the rotational subsystems where they are processed as a series of multiple SISO systems. For the position control along the X and Y directions, the control is indirect in the sense that a decoupling algorithm needs to calculate the specified Roll $\left({\varphi }_d\right)$ and pitch $\left({\theta }_d\right)$ that are passed onto the attitude controller to successfully track the desired trajectories.

3.2.1. Altitude Controller

The Altitude subsystem has the following principle model:

$$\begin{array}{cc}\hfill {\dot{x}}_7& =\dot{z}=x_8\hfill \\ \hfill {\dot{x}}_8& =\ddot{z}=-g+\left(C_{\varphi }{C}_{\theta }\right){\displaystyle \frac{U_1}{m}}\hfill \end{array}$$

(60)

Converting this into the U-control realisation, the Altitude subsystem is

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \ddot{z}& ={\lambda }_0+{\lambda }_1{U}_1\hfill \\ \hfill {\lambda }_0& =-g\hfill \\ \hfill {\lambda }_1& ={\displaystyle \frac{\left(C_{\varphi }{C}_{\theta }\right)}{m}}\hfill \end{array}\hfill \end{array}\right.$$

(61)

The control input term $U_4$ can now be calculated as follows:

$$U_1={\displaystyle \frac{\ddot{z}-{\lambda }_0}{{\lambda }_1}}$$

(62)

3.2.2. Model-Based Position U-Controller

From Figure 5, it is clear that the outer loop deals with the positioning by generating the desired setpoints for Roll $\left({\varphi }_d\right)$ and pitch $\left({\theta }_d\right)$, which are then passed onto the rotational subsystem as reference co‘mmands. The desired trajectories $x_d$ and $y_d$ are passed on to the position controller along with the Yaw output and thrust as the change in Yaw affects the heading, and, therefore, trajectories need to be calculated accordingly and the vertical and horizontal thrust components from Roll and pitch manoeuvres will affect the Altitude thrust. In order to address these issues, virtual inputs $u_x$, $u_y$ and $u_z$ are used as they simplify and decouple the complex underactuated dynamics of the quadrotor, making it easier to design and implement control strategies. The decoupling algorithm simplifies the design of controllers because each virtual input can be associated with a specific axis or motion:

• $u_x$ controls the Roll and lateral motion (sideways movement).
• $u_y$ controls the pitch and longitudinal motion (forward/backward movement).
• $u_z$ controls the thrust and vertical motion (up/down movement).

One can formulate the required pitch ${\theta }_d$ and $Roll{\varphi }_d$ angles using the following decoupling algorithm [29]:

$$\left\{\begin{array}{c}\begin{array}{c}\hfill {\varphi }_d=arcsin\left(m{\displaystyle \frac{u_{x}sin\left(\psi \right)-u_{y}cos\left(psi\right)}{u_1}}\right)\\ \hfill {\theta }_d=arctan\left({\displaystyle \frac{u_{x}cos\left(\psi \right)+u_{y}sin\left(psi\right)}{u_z}}\right)\end{array}\hfill \end{array}\right.$$

(63)

The position of the system is determined by the following subsystem:

$$\begin{array}{ccc}\hfill \left[\begin{array}{c}\ddot{x}\\ \ddot{y}\end{array}\right]& =\left[\begin{array}{c}(C_{\varphi }{S}_{\theta }{C}_{\psi }+S_{\varphi }{S}_{\psi })\\ (C_{\varphi }{S}_{\theta }{S}_{\psi }-S_{\varphi }{C}_{\psi })\end{array}\right]{\displaystyle \frac{U_1}{m}}\hfill & \hfill =\left[\begin{array}{c}{u}_x\\ u_y\end{array}\right]\end{array}$$

(64)

Converting these into the U-control realisation, the position control subsystem is as follows:

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \ddot{x}& ={\lambda }_0+{\lambda }_1{u}_x;\mathrm{where}{\lambda }_0=0;{\lambda }_1=1\hfill \\ \hfill \ddot{y}& ={\lambda }_0+{\lambda }_1{u}_y;\mathrm{where}{\lambda }_0=0;{\lambda }_1=1\hfill \end{array}\hfill \end{array}\right.$$

(65)

The virtual control inputs $u_y$ and $u_y$ can now be calculated as

$$\begin{array}{ccc}\hfill u_x& ={\displaystyle \frac{v-{\lambda }_0}{{\lambda }_1}};u_y\hfill & \hfill ={\displaystyle \frac{v-{\lambda }_0}{{\lambda }_1}}\end{array}$$

(66)

Figure 5. Schematic of the proposed control approach.

Figure 5. Schematic of the proposed control approach.

Implementing the derived control laws, the simulated system performance can be seen in Figure 6. It is noted that the model-based U-dynamic inversion approach is sensitive to model inaccuracies and its performance deteriorates in the presence of unmodelled dynamics or disturbances, resulting in model–inverse mismatch (i.e., $G_p\times G_p^{-1}\ne 1$). To accommodate the problem, a robust compensator with a Quasi-Sliding Mode [28] is added to strengthen the stability and eliminate convergent error.

3.3. Robust Compensator Design Procedure

Following from Section 2.4, a robust Sliding Mode-based Compensator is implemented to improve the performance of the UMDI due to its sensitivity to model inaccuracies. Considering the Yaw subsystem from earlier, defined by its state space Equation (57), U-state space realisation (58) and finally the ideal inversion defined by Equation (59),

$$\begin{array}{cc}\hfill {\dot{x}}_5& =\dot{\psi }=x_6\hfill \\ \hfill {\dot{x}}_6& =\ddot{\psi }={\displaystyle \frac{J_{xx}-J_{yy}}{J_{zz}}}\dot{\varphi }\dot{\theta }+{\displaystyle \frac{U_4}{J_{zz}}}\hfill \end{array}$$

(67)

$$\left\{\begin{array}{c}\begin{array}{cc}\hfill \ddot{\psi }& ={\lambda }_0+{\lambda }_1{U}_4\hfill \\ \hfill {\lambda }_0& ={\displaystyle \frac{J_{xx}-J_{yy}}{J_{zz}}}\dot{\varphi }\dot{\theta }\hfill \\ \hfill {\lambda }_1& ={\displaystyle \frac{1}{J_{zz}}}\hfill \end{array}\hfill \end{array}\right.$$

(68)

The ideal inversion output term $U_4$ is calculated directly as there is a direct relationship between input u and output $\psi$:

$$U_4={\displaystyle \frac{\ddot{\psi }-{\lambda }_0}{{\lambda }_1}}={\displaystyle \frac{v-{\lambda }_0}{{\lambda }_1}}$$

(69)

Based on Equations (34) and (35), the error equations for $m=2$ (i.e., a second-order system) are as follows:

$$\left\{\begin{array}{c}{\tilde{q}}_1=\psi -v{\left({\displaystyle \frac{1}{s}}\right)}^2\hfill \\ {\tilde{q}}_2={\dot{\tilde{q}}}_1=\dot{\psi }-v\left({\displaystyle \frac{1}{s}}\right)\hfill \end{array}\right.$$

(70)

The sliding mode manifold $\sigma$ and its derivative $\dot{\sigma }$ are given by

$$\begin{array}{cc}\hfill \sigma & =c{\tilde{q}}_1+{\tilde{q}}_2\hfill \\ \hfill \dot{\sigma }& =c{\dot{\tilde{q}}}_1+{\dot{\tilde{q}}}_2=c\left(\dot{\psi }-v{\displaystyle \frac{1}{s}}\right)+(\ddot{\psi }-v)\hfill \end{array}$$

(71)

Replacing the term $\ddot{\psi }$ with its U-model equivalent from Equation (58) gives

$$\dot{\sigma }=c\left(\dot{\psi }-v{\displaystyle \frac{1}{s}}\right)+({\lambda }_0+{\lambda }_{1}u-v)$$

(72)

Setting $\dot{\sigma }=0$ to calculate $u_{eq}$ gives

$$u_{eq}={\displaystyle \frac{v-c\left(\dot{\psi }-v\frac{1}{s}\right)-{\lambda }_0}{{\lambda }_1}}$$

(73)

$$u_{sw}=-(\epsilon sign\left(\sigma \right)+k\sigma )$$

(74)

From Equation (38),

$$\begin{array}{cc}\hfill u_c& =u_{sw}+u_{eq}-u_i\hfill \\ \hfill & =-(\epsilon sign\left(\sigma \right)+k\sigma )+{\displaystyle \frac{v-c\left(\dot{\psi }-v\frac{1}{s}\right)-{\lambda }_0}{{\lambda }_1}}-{\displaystyle \frac{(v-{\lambda }_0)}{{\lambda }_1}}\hfill \\ \hfill & =-(\epsilon sign\left(\sigma \right)+k\sigma )+{\displaystyle \frac{\cancel{v}-c\left(\dot{\psi }-v\frac{1}{s}\right)-\cancel{{\lambda }_0}-\cancel{v}+\cancel{{\lambda }_0}}{{\lambda }_1}}\hfill \\ \hfill & =-(\epsilon sign\left(\sigma \right)+k\sigma )+{\displaystyle \frac{-c\left(\dot{\psi }-v\frac{1}{s}\right)}{{\lambda }_1}}=-\left(\epsilon {\displaystyle \frac{\sigma }{\left|\sigma \right|+\xi }}+k\sigma \right)+{\displaystyle \frac{-c\left(\dot{\psi }-v\frac{1}{s}\right)}{{\lambda }_1}}\hfill \end{array}$$

(75)

where $c\in {\mathbb{R}}^{+},k\in {\mathbb{R}}^{+}$ and $\epsilon \in {\mathbb{R}}^{+}$ are design parameters associated with the Sliding Mode Compensator. The results after the implementation of the compensator can be observed in Figure 7 showcasing elimination of the steady state error and improved response time.

Following a similar approach, robust compensators for the Roll, pitch and Altitude channels can be designed using the design parameters in

Table 2. Compensator design parameters.

	c	k	$\mathit{\xi }$
Altitude Z (m)	10	1	0.1
Yaw ($\psi$) (radians)	10	5	0.1
Roll ($\varphi$) (radians)	10	10	0.03
Pitch ($\theta$) (radians)	10	10	0.03

.

Remark 4. 

Equation (75) gives the compensated controller output $u_c$, where $u_{sw}$ represents the switching controller for driving the system to the sliding surface (or sliding boundary in this specific case), $u_{eq}$ represents the equivalent control that ensures the system state stays inside the sliding boundary and $u_i$ represents the ideal plant inversion calculated using the UMDI controller.

To compare the simulated performance of the designed controllers quantitatively, the Root Mean Square values of the simulated tracking errors are calculated in

Table 3. Root Mean Square Errors (RMSEs) in trajectory following.

	PID	UMDI	USMC
X (m)	0.1026	0.2749	0.1113
Y (m)	0.1589	0.2885	0.0860
Z (m)	0.0980	0.1159	0.1244
Roll ($\varphi$) (radians)	0.0025	0.0304	0.0107
Pitch ($\theta$) (radians)	0.0021	0.0327	0.0105
Yaw ($\psi$) (radians)	0.0543	0.0809	0.0715

. Figure 8 shows the flight trajectories in 3D and 2D coordinates. Figure 9 shows the time sequence of the trajectories in X, Y, Z directions, and Yaw channel respectively. It can be observed that while both PID and the USMC controllers can drive the quadrotor to the desired trajectory, the U-dynamic inversion-based controller suffers from deteriorated performance due to unmodelled dynamics/parameters and therefore is unable to accurately track the desired trajectory.

Flight tests to verify the simulation results were performed on the Parrot Mambo Minidrone using the same controller parameters. While the fully actuated subsystem (Altitude and Yaw) showed good tracking (Figure 10), the position tracking was unstable and exhibited very high chattering, suggesting that the measures taken while designing the controller, e.g., the introduction of the sigmoid function to replace the discontinuous function $sign\left(\sigma \right)$ to reduce system chattering, were inadequate. The outer loop that generates the desired pitch and Roll angles for position control (Equation (63)) uses the Altitude controller output terms $u_1,u_z$ and Yaw $\psi$; chattering in these terms generates unstable trajectories, causing the minidrone to become unstable. Since the minidrone is a very light weight and has a small test platform with very low inertia, it cannot handle such mechanical chattering and therefore becomes unstable. However, the designed controller did show improved performance in both Yaw and Altitude channels, reducing the steady state error and achieving an overall smoother hover performance, as can be seen in Figure 11.

Table 4. Root Mean Square Errors (RMSEs) in trajectory following for real flight tests.

	PID	USMC
Z (m)	0.0501	0.0138
Yaw ($\psi$) (radians)	0.0284	0.0177

shows the effectiveness of these experimental results and compares them to the simulated results from earlier by using the RMS errors. It can further be observed from Figure 12 that the experimental results follow the simulated results more closely using the proposed approach when compared to a PID controller. This can be partially explained by the fact that the PID gains that were tuned for simulations need further tuning for flight tests, i.e., the performance of PID controllers varies in simulated and real flight tests. In contrast, the proposed approach after simulations is quite straightforward and does not need any further modifications for successful implementation during flight tests.

4. Discussions

Comparing the results of a tuned PID controller, a UMDI controller and the compensated robust USMC, it can be concluded that the UMDI controller results in deteriorated performance in the presence of unmodelled disturbances, but the addition of a robust compensator improves the performance, as can be seen in Table 3 where the RMSE is taken to statistically analyse the performance on controllers. It can be concluded that if a system is modelled accurately without outside disturbance, the UMDI-based controller is effective since there is no model–inversion mismatch and, therefore, $G_p.{G_p}^{-1}=1$. However, real dynamic systems always have some outside interference and unmodelled dynamics; therefore, these disturbances and model inaccuracies need to be taken into account while designing a controller. Therefore, some sort of a compensator is a must while using this approach to impart some robustness to the system. The Sliding Mode-based Compensator, while ensuring robustness, introduces chattering. This is partly fixed by using a Quasi-Sliding Mode control to approximate the $sign(\ast )$ function; however, that does not drive the dynamics to the sliding surface but rather to a sliding boundary, therefore possibly affecting system performance. Other possible ways to avoid chattering effects need to be studied and integrated with U-control possibly with a higher-order sliding mode approach like super twisting sliding mode. The resultant simulation plots show that the new procedure successfully addresses the challenges involved, like system underactuation, coupling, unmodelled dynamics and disturbances.

There are some warnings that one needs to be aware of in the implementation of the model-based U-control for UAV:

• The actuator saturation could limit the SMC for large control inputs.
• Comparing the model-free trial-and-error PID (which is time-consuming in the gain tuning and difficult to achieve desired performance) embedded in the Simulink, the model-based U-control requires a nominal model for each different type of UAV, so it requires more design and implementation expertise/cost in calibration.
• More effective model-free control with merits in concise implementation and easy robust tuning could be considered in the next stage of this research.

5. Conclusions

This paper integrates a U-model-based dynamic inversion with an SMC compensator to control the trajectory of a quadrotor system flight operation. The designed controller shows strong stability and accurate setpoint tracking, even facing the model uncertainties.

On the whole, after comparing the results and relaxing implementation constraints, the built-in PID control is still preferred over advanced approaches due to its simplicity, even if tuning it involves some trial-and-error and guesstimating. The unmodelled dynamics and especially the external disturbances in the Roll and pitch channels destabilize the system and deteriorate its performance. This can affect the system performance, greatly preventing its implementation on physical systems. While the designed compensator improves the performance of the controller greatly, it might not be a feasible solution when running hardware-in-loop tests due to the chattering effects introduced by sliding mode control. Therefore, to ensure robustness along with practical applicability there is a need to integrate a chattering-free (or reduced chattering) alternative controller that can provide robustness and improve performance. Possible future work should be focused on this where the bounded disturbances/model inaccuracies can be lumped together and studied. Another alternative would be to implement a chattering-free robust controller like a super twisting [28] controller design in future work. Model-free control approaches [30,31,32] can also provide a subject for exploration and their integration with U-control can improve overall performance and system design simplicity. Furthermore, comprehensive tests of the different UAVs are on the cards for future studies.