Barrier Lyapunov function and adaptive backstepping-based control of a quadrotor UAV

Abstract

This paper presents backstepping control and backstepping constraint control approaches for a quadrotor unmanned aerial vehicle (UAV) control system. The proposed methods are applied to a Parrot Mambo drone model to control rotational motion along the $x$, $y$, and $z$ axes during hovering and trajectory tracking. In the backstepping control approach, each state of the system controls the previous state and is called “virtual control.” The last state is controlled by the real control input. The idea is to compute, in several steps, a control law that ensures the asymptotic stability of the system. The backstepping constraint control method, based on barrier Lyapunov functions (BLFs), is designed not only to track the desired trajectory but also to guarantee no violation of the position and angle constraints. Symmetric BLFs are introduced in the design of the controller. A nonlinear mathematical model is considered in this study. Based on Lyapunov stability theory, it can be concluded that the proposed controllers can guarantee the stability of the UAV system and the state converges asymptotically to the desired trajectory. To make the control robust, an adaptation law is applied to the backstepping control that estimates the unknown parameters and ensures their convergence to their respective values. Validation of the proposed controllers was performed by simulation on a flying UAV system.

1. Introduction

Unmanned aerial vehicles (UAVs) can fly autonomously and are useful for applications such as inspection of areas dangerous to humans, detection of toxic gases or radiation, and monitoring in sectors such as agriculture, mining, and construction. They can also be used for aerial photography and even for pizza or parcel deliveries. Great interest in these flying robots, commonly called drones, has motivated a considerable body of research to improve their movement control. Studies have shown that the quadrotor UAV configuration is superior to others, due to its structural simplicity, reduced size, low weight, and capacity for vertical takeoff and landing. A wide variety of UAVs are now available, but the newest and most popular of these are quadrotors. Due to their compact mechanical structure and low weight, they can be deployed in very complex environments. However, the flying robot remains a nonlinear, multivariable, coupled, and underactuated system. Therefore, its control system warrants further study, if it is to be stable and robust. Recently, many control strategies have been developed and applied for UAV stabilization and trajectory tracking, including feedback linearization (FBL) [1], [2], model predictive control (MPC) [3], [4], and [5], sliding mode control (SMC) [6], [7], adaptive control [8], [9], and backstepping [10], [11], to name a few.

The authors of ref. [12] used an FBL-based method for the attitude and altitude control dynamics of the nonlinear inner loop, and a linear quadratic regulator (LQR) with an integral action for the outer loop of an UAV, to follow a predefined trajectory. However, the FBL-based controller is quite sensitive to sensor noise as well as modeling uncertainty. FBL and LQR methods were applied in ref. [2] to analyze and stabilize a quadrotor UAV system. Particle swarm optimization was used to achieve the optimal weighting matrices of LQR such that the composite controller can reduce the amplitudes of system control inputs. In ref. [13], the authors applied MPC to obtain a precise trajectory in the presence of forcible wind gusts. The MPC approach is useful because of its optimal tracking characteristics [14]. In ref. [15], the authors empirically compared two control strategies, the nonlinear model predictive controller (NMPC) and the differential-flatness-based controller, by tracking a wide variety of agile trajectories. Comparisons were performed in both simulation and real-time environments to systematically evaluate both methods in terms of their tracking accuracy, robustness, and computational efficiency. The comparisons showed the superiority of the NMPC in tracking trajectories, at the cost of higher computation time and a risk of numerical convergence issues. A saturated adaptive SMC strategy was presented in ref. [16] for attitude and altitude stabilization of a quadrotor UAV, taking input saturation into consideration. Compared to conventional control strategies, this approach provided superior tracking performance and disturbance and uncertainty rejection. An adaptive fuzzy gain-scheduling SMC for attitude regulation of a quadrotor UAV with parametric uncertainties and external disturbances was considered in ref. [17]. A continuous SMC with exponential convergence, combined with a finite-time sliding mode observer estimating full states and identifying disturbances, was proposed in ref. [18]. A nonsingular fast terminal (SMC) was suggested in ref. [19] to avoid singularities and achieve tracking control of a quadrotor with actuator faults. A robust adaptive nonsingular fast terminal (SMC) was designed in ref. [9] for attitude and displacement tracking control of a quadrotor subjected to unknown modeling errors.

In ref. [20], a backstepping controller was proposed to control the altitude and attitude of quadrotor systems. Controller validity was ensured by the Lyapunov function, and simulation results showed a high precision transient and tracking response. In ref. [21], a backstepping controller for quadrotor tracking based on a neural network approach was proposed in the presence of uncertainties. A nonlinear controller combining the SMC technique for the attitude subsystem (inner loop) and the backstepping control technique for the position loop until desired attitudes are obtained was proposed in ref. [22]. In addition, an adaptive observer-based fault estimation scheme was considered for takeoff mode. Simulations were conducted to demonstrate the effectiveness of the designed robust nonlinear controller and the fault estimation scheme. In ref. [5], a finite-time adaptive integral backstepping fast terminal sliding mode control (FTSMC) was presented for the purpose of finite-time attitude/position tracking control of a quadrotor under model uncertainty and external disturbance. In ref. [23], an adaptive control strategy, based on the modified function approximation technique, was proposed for a manipulator robot with unknown dynamics. This strategy applies the backstepping control approach and uses state and output feedback to control the manipulator trajectory.

One of the issues in controlling a dynamical system is handling constraints. Various tools have been developed for this purpose, including tools for nonlinear systems. Recently, barrier Lyapunov functions (BLFs) were proposed to handle constraints, driven by the technique of tailoring the Lyapunov function according to the requirements of the problem. In order to solve the above-mentioned problems, various control approaches have been proposed. In ref. [24], a BLF-based control method taking output constraints into consideration was presented for a quadrotor, and a satisfactory tracking performance was achieved in steady state. A BLF-based finite-time backstepping control of a quadrotor with full state constraints was considered in ref. [25]. To prevent constraint violation, [26] used a BLF, which adapts to infinity when its arguments approach the desired limits. Using a BLF guarantees that these limits will never be reached. In refs. [27] and [28], a governor method was designed to satisfy the constraints of uncertain nonlinear systems, in which the input of the closed loop system is the reference signal, and modification of the reference signal does not violate the state and control constraints. In ref. [29], BLF-based robust backstepping control was applied to Multi-rotor unmanned aerial vehicles (MUAVs) to counter external disturbances. Finite-time stabilization using BLF was achieved in ref. [30] in the presence of parametric uncertainties. The integral BLF was used in conjunction with backstepping and vision-based control in ref. [31] to restrict drastic attitude motion. In ref. [32], BLF-based backstepping was implemented on a quadrotor and its performance compared to a proportional, integral and derivative controller (PID) and a SMC. Some recent developments have been made in the design of SMCs for output-constrained systems. A second-order SMC for an output-constrained system under external disturbances was developed in ref. [33]. This approach based on BLF and finite-time convergence was attained while the system output remained within the desired constraints. In this paper, we will apply backstepping control and backstepping constraint control approaches for a UAV quadrotor control system. The objective is to find a control technique that solves the problem of tracking the reference trajectory while respecting a set of constraints. Motivated by the statement and challenges presented above, the main contributions of this paper can be summarized as follows:

• • Based on the above methods, BLF is introduced into quadrotor dynamics to solve the problem of desired trajectory tracking with position and angle symmetric constraints. Using the BLF method, improved controller performance and outputs were limited within constraints.

• • Lyapunov-based analysis is employed to prove that the proposed control design makes the UAV track the desired trajectory.

• • A comparative study between the three methods is conducted, illustrating the advantages of each technique. The simulation results show that our BLF control algorithm is more effective than the classical and adaptive backstepping control algorithms.

This paper is organized as follows. In Section 2, the kinematic and dynamic model of the quadrotor is addressed. Section 3 presents the control objective. Section 4 introduces the backstepping control and considers parametric uncertainty in the model. In Section 5, we extend the designs based on symmetric BLFs. A comparison study of the different backstepping approaches is presented in Section 6. Section 6 presents and discusses the results of the numerical simulation. Finally, conclusions and possibilities for future research are presented in Section 7.

Figure 1. Quadrotor structure, forces, angles, and frames (Parrot Mambo).

2. Quadrotor dynamic

The study of how quadrotor movement evolves through space requires knowledge of both its position and orientation. Its localization requires the choice of at least two references: a global reference frame ( $R_0$ ) and a local reference frame ( $R_G$ ) attached to the quadrotor’s center of gravity (Fig. 1). The quadrotor is a flying robot with six degrees of freedom, three translations $(x, y, z)$ , and three rotations $(\phi, \theta, \psi )$ . The transition matrix from the reference frame $R_0$ to the $R_G$ frame is given by the product of three successive matrices and calculated by:

(1) \begin{align} R_{\psi \theta \phi }=&\left ( \begin{array}{c@{\quad}c@{\quad}c} c\theta c\psi & s\theta s\phi c\psi -s\psi c\phi & s\theta c\phi c\psi +s\psi s\phi \\[4pt] c\theta s\psi & s\theta s\phi s\psi +c\psi c\phi & s\theta c\phi s\psi -c\psi s\phi \\[4pt] -s\theta & c\theta s\phi & c\theta c\phi \\ \end{array} \right ) \end{align}

where $c=\cos (.)$ and $s=\sin (.)$ . For the quadrotor, a roll and pitch movement equal to $\frac{\pi }{2}$ is excluded.

Applying Newton’s law, the dynamics of the quadrotor is given by [34]:

(2) \begin{equation} \begin{array}{r@{\quad}c@{\quad}l} \ddot x&=&\dfrac{1}{m}(s\theta c\psi c\phi + s\psi s\phi )U_1\\[8pt] \ddot y&=&\dfrac{1}{m}(s\theta s\psi s\phi -c\psi s\phi )U_1\\[8pt] \ddot z&=& \dfrac{1}{m}(mg-U_1c\theta c\phi ) \\[8pt] \ddot \phi &=& \dfrac{1}{I_x}\left[(I_y-I_z)\dot{\theta }\dot{\psi }-J_r\dot{\theta }w_r+U_2\right]\\[8pt] \ddot \theta &=& \dfrac{1}{I_y}\left[(I_z-I_x)\dot{\phi }\dot{\psi }+J_r\dot{\phi }w_r+U_3\right]\\[8pt] \ddot \psi &=& \dfrac{1}{I_z}\left[(I_x-I_y)\dot{\theta }\dot{\phi }+U_4\right] \end{array} \end{equation}

where $(x, y, z)$ denotes the quadrotor’s motion along the $x$ , $y$ , and $z$ axes, respectively, $(\phi, \theta, \psi )$ denotes the roll, pitch, and yaw angles, respectively, $m$ denotes the mass of the quadrotor, and $(I_x, I_y, I_z)$ denotes the quadrotor vehicle’s moment of inertia. Quadrotors have four rotors at the ends of their four axes. These rotors are used to control the machine by acting on their rotational speeds. $U_1$ denotes the heave force of the quadrotor vehicle, and $(U_2, U_3, U_4)$ denotes the moments applied to the center of vehicle’s mass. $U_1$ , $U_2$ , $U_3$ , $U_4$ can be expressed as follows [34]:

(3) \begin{equation} \begin{cases} U_1=b \left(\omega _{1}^{2}+\omega _{2}^{2}+\omega _{3}^{2}+\omega _{4}^{2}\right)\\[4pt] U_2=bl \left(\omega _{1}^{2}-\omega _{2}^{2}-\omega _{3}^{2}+\omega _{4}^{2}\right)\\[4pt] U_3=bl \left(\omega _{1}^{2}+\omega _{2}^{2}-\omega _{3}^{2}-\omega _{4}^{2}\right)\\[4pt] U_4=d \left({-}\omega _{1}^{2}+\omega _{2}^{2}-\omega _{3}^{2}+\omega _{4}^{2}\right)\\[4pt] \omega _r=\omega _{4}+\omega _{2}-\omega _{1}-\omega _{3} \end{cases} \end{equation}

where $l$ is the distance between the center of mass and the rotation axis of the motors rotor, $b$ represents the lift coefficient, $d$ represents the torque coefficient, and $(\omega _{1},\omega _{2},\omega _{3},\omega _{4})$ are the rotational speed of the rotors. Equation (2) can be rewritten in the state space form as:

(4) \begin{equation} \dot x=f(x,u)=\left ( \begin{array}{c} x_2\\[4pt] a_1x_4x_6+a_3x_4\omega _r+b_1U_2\\[4pt] x_4\\[4pt] a_4x_2x_6+a_6x_2\omega _r+b_2U_3\\[4pt] x_6\\[4pt] a_7x_2x_4+b_3U_4\\[4pt] x_8\\[4pt] g-\dfrac{1}{m} c x_1 c x_3 U_1\\[4pt] x_{10}\\[4pt] \dfrac{1}{m} U_x U_1\\[4pt] x_{12}\\[4pt] \dfrac{1}{m} U_y U_1 \end{array} \right ) \end{equation}

where $u=[U_1,U_2,U_3,U_4]^T$ is the input torques vector, $x=[x_1,\ldots,x_{12}]^T=[\phi,\dot \phi,\theta,\dot \theta,\psi,\dot \psi,z,\dot z,x,\dot x,y,\dot y]^T$ is the state vector, and

\begin{equation*} \begin{array}{l@{\quad}l@{\quad}l} a_1=\dfrac {(I_y-I_z)}{I_x}, & a_3=\dfrac {-J_r}{I_x}, &b_1=\dfrac {1}{I_x}\\[15pt] a_4= \dfrac {(I_z-I_x)}{I_y}, & a_6=\dfrac {J_r}{I_y}, &b_2=\dfrac {1}{I_y}\\[15pt] a_7 =\dfrac {(I_x-I_y)}{I_z}, & b_3=\dfrac {1}{I_z}, &b_4=\dfrac {1}{m}\\ \end {array} \end{equation*}

\begin{equation*} U_x= cx_1 sx_3 sx_5 - sx_1 cx_5, U_y= cx_1 sx_3 cx_5 + sx_1 sx_5 \end{equation*}

3. Control objective

The control objective of this paper is to design the control input $u$ for the stratospheric UAV described by (4) in the presence of unknown parameters or state constraints such that:

• • The AUV position $x,y$ , and $z$ moves along the desired trajectory $x_d,y_d$ , and $z_d$ , and the trajectory tracking errors converge to a sufficiently small proximity to zero.

• • The output constraint is not violated.

The control strategy of the quadrotor is illustrated in Fig. 2. The quadrotor is an underactuated system. Four rotors are used to control six degrees of freedom. The control strategy is based on two loops of control. The first is the internal control loop which controls roll $\phi$ , pitch $\theta$ , yaw $\psi$ , and altitude $z$ as shown in Fig. 2. The internal control loop uses the reference value to generate the proper control signal. The second is the external control loop which controls the $x$ and $y$ position. The purpose of the external control is to calculate the desired roll $\phi _d$ and the desired pitch $\theta _d$ based on the desired position by using the positions controls $u_x, u_y$ and the desired yaw $\psi _d$ . The desired position $(x_d,y_d,z_d)$ and the desired yaw $\psi _d$ come directly from the user.

Figure 2. Control block diagram.

4. Backstepping control

In order that this paper be self-contained, the attitude control movements around the $x$ , $y$ , and $z$ axes and the attitude movement are elaborated in this section. The principle of backstepping is to divide the system into several sub-systems in a cascade form. The control laws are then made for each subsystem, in a decreasing manner, until a global control law for the whole system is generated.

4.1. Control of movement around $x$ axis

Our objective is to establish a command $U_2$ that will orient the quadrotor by the angle $\phi$ . Considering the first subsystem mentioned below:

(5) \begin{align} \dot{x}_1&=x_2\nonumber \\ \dot{x}_2&=a_1x_4x_6+a_3x_4\omega _r+b_1U_2 \end{align}

• • First iteration: To produce an output that follows a desired trajectory $\phi _d$ , we define the trajectory error

  (6) \begin{align} e_{x_1}=&\phi _d-x_1 \nonumber \\ \dot{e}_{x_1}=&\dot{\phi }_d-\dot{x}_1 \end{align}

We consider a positive definite Lyapunov candidate function $V_1(x)$ :

\begin{equation*} V_1=\frac {1}{2}e_{x_1}^{2} \end{equation*}

Its time derivative is

\begin{equation*} \dot {V}_{1}=e_{x_1}(\dot {\phi }_d-x_2) \end{equation*}

The state variable $x_2$ is then used as an intermediate command to satisfy the Lyapunov condition:

\begin{equation*} \dot {V}_{1} = -k_{1}e_{x_1}^{2} \end{equation*}

with $k_{1}\gt 0$ . For this purpose, we define a virtual command $x_{2d}$ :

(7) \begin{align} \begin{split} x_{2d} =&\dot{\phi }_d+k_{1}e_{x_1} \\ \dot{x}_{2d} =&\ddot{\phi }_d+k_{1}\dot{e}_{x_1} \end{split} \end{align}

• • Second iteration: At this stage, a new error is generated:

  (8) \begin{align} e_{x_2}=&x_{2d}-x_2=\dot{\phi }_d+k_{1}e_{x_1}-\dot \phi \end{align}
  (9) \begin{align} \dot{e}_{x2}=&\ddot{\phi }_d-\ddot \phi +k_{1}\dot{e}_{x_1} \end{align}

From Eqs. (6) and (8), we find

(10) \begin{align} \dot{e}_{x_1}=&e_{x_2}-k_{1}e_{x_1} \end{align}

To eliminate this error, the Lyapunov candidate function $V_1$ is augmented by another term, which will address the new error that was introduced earlier:

(11) \begin{align} V_2=\frac{1}{2}\left(e_{x_1}^{2}+e_{x_2}^{2} \right) \end{align}

(12) \begin{align} \dot{V}_{2}=e_{x_1}\dot{e}_{x_1}+e_{x_2}\dot{e}_{x_2} \end{align}

By replacing (9) and (10) in (12), we obtain

(13) \begin{align} \dot{V}_{2}& =e_{x_1}\left(e_{x_2}-k_{1}e_{x_1}\right)+e_{x_2}\left(\ddot{\phi }_d-\ddot \phi +k_{1}\dot{e}_{x_1}\right)\nonumber \\ &=-k_{1}e_{x_1}^{2}+e_{x_2}\left(\ddot{\phi }_d-\ddot \phi +k_{1}\left(e_{x_2}-k_{1}e_{x_1}\right) +e_{x_1}\right) \end{align}

From (4), by replacing the expression of $\ddot \phi$ in (13), we obtain

(14) \begin{align} \dot{V}_{2}=& -k_{1}e_{x_1}^{2}+e_{x_2}\left(\ddot{\phi }_d - a_1 x_4 x_6 -a_3x_4\omega _r - b_1 U_2 +k_{1}\left(e_{x_2}-k_{1}e_{x_1}\right)+e_{x_1}\right) \end{align}

In order to satisfy the Lyapunov condition ( $\dot V_2(e_{x1},e_{x2})\leq 0$ ), the following expression is used:

(15) \begin{align} \ddot{\phi }_d - & a_1x_4x_6 -a_3x_4\omega _r -b_1U_2+k_{1}\left(e_{x_2}-k_{x1}e_{x_1}\right)+e_{x1} = -k_{2}e_{x_2} \end{align}

with $k_{2}\gt 0$ . Then,

(16) \begin{align} U_2 = &\frac{1}{b_1} \left[\ddot{\phi }_d-a_1x_4x_6-a_3x_4\omega _r+k_{1}\left(e_{x_2}-k_{1}e_{x_1}\right)+e_{x_1}+k_{2}e_{x_2}\right] \end{align}

so that

(17) \begin{align} \dot{V}_{2}=-\left(k_{1}e_{x_1}^{2}+k_{2}e_{x_2}^{2}\right)\leq 0 \end{align}

By using Lyapunov theory, we conclude that the orientation error dynamics is asymptotically stable, and that the orientation angle $\phi$ of the quadrotor follows asymptotically its desired angle.

4.2. Control of movement on $y$ and $z$ axes

Following the same procedure, we obtain

\begin{eqnarray*} U_1&=&\frac{m}{cx_1 cx_3}\left [g-\ddot{z}_d-k_{x7}(e_{x_8}-k_{7}e_{x_7}) -e_{x_7}-k_{8}e_{x_8} \right ] \\ U_3&=&\frac{1}{b_2} \left[\ddot{\theta }_d-a_4x_2x_6-a_6x_2\omega _r+k_{3}(e_{x_4}-k_{3}e_{x_3})+ e_{x_3}+k_{4}e_{x_4}\right]\\ U_4&=&\frac{1}{b_3} \left[\ddot{\psi }_d-a_7x_2x_4+k_{5}(e_{x_6}-k_{5}e_{x_5})+ e_{x_5}+k_{6}e_{x_6}\right] \end{eqnarray*}

Note that $U_1$ is not singular because, as previously mentioned, a roll and pitch movement equal to $\frac{\pi }{2}$ is excluded.

5. Adaptive backstepping

In this section, a backstepping adaptive control will be applied to a quadrotor system assuming unknown parameters to ensure robustness. Adaptive backstepping is a technique where the control parameters vary in real time to adapt to changes in the system dynamics. These changes can be due to modeling errors, parametric uncertainties, or external disturbances.

5.1. Rotation command (around $x$ axis)

Let $\widehat{a}_1$ , $\widehat{a}_3$ , and $\widehat{b}_1$ be the estimates of $a_1$ , $a_3$ , and $b_1$ , respectively. Then, the dynamics (4) becomes

\begin{align*} \dot{x}_1&= x_2 \\ \dot{x}_2&=\left(\widetilde{a}_1+\widehat{a}_1\right)x_4x_6+\left(\widetilde{a}_3+\widehat{a}_3\right)x_4\omega _r+\left(\widetilde{b}_1+\widehat{b}_1\right)U_2 \end{align*}

with: $\widetilde{a}_1=a_1-\widehat{a}_1$ , $\widetilde{a}_3=a_3-\widehat{a}_3$ , and $\widetilde{b}_1=b_1-\widehat{b}_1$ .

• • First iteration: The first iteration follows the first iteration as in the regular backstepping above, that is, Eqs. (6) and (7) are the same.

• • Second iteration: At this stage, a new error is generated:

  (18) \begin{align} e_{x_2}=x_{2d}-x_2=\dot{\phi }_d -\dot \phi +k_{1}e_{x_1} \end{align}
  (19) \begin{align} \dot{e}_{x_2}=\ddot{\phi }_d-\ddot \phi +k_{1}\dot{e}_{x_1} \end{align}

From Eqs. (6) and (18), we obtain

(20) \begin{align} \dot{e}_{x_1}=e_{x_2}-k_{1}e_{x_1} \end{align}

From Eqs. (4), (19), and (20), we obtain

(21) \begin{align} \dot{e}_{x_2}&= \ddot{\phi }_d-\ddot{\phi }+k_{x1}\left(e_{x2}-k_{1}e_{x_1}\right) \nonumber \\ &= \ddot{\phi }_d-\left(\widetilde{a}_1 + \widehat{a}_1\right)x_4x_6 - \left(\widetilde{a}_3+\widehat{a}_3\right)x_4\omega _r\nonumber \\ &\quad -\left(\widetilde{b}_1+\widehat{b}_1\right)U_2+k_{x1}\left(e_{x_2}-k_{1}e_{x_1}\right) \end{align}

To develop the adaptation law of $a_1$ and the control input $U_2$ , we make a change of variable, such that $b_1 U_2=\overline{U}_2 - \widetilde{b}_1 U_2$ . Then

(22) \begin{align} \dot{e}_{x_2}&=\ddot{\phi }_d-\left(\widetilde{a}_1+\widehat{a}_1\right)x_4x_6-\left(\widetilde{a}_3+\widehat{a}_3\right)x_4\omega _r \nonumber \\ & \quad -\overline{U}_2 + \widetilde{b}_1U_2+k_{x1}\left(e_{x_2}-k_{1}e_{x_1}\right) \end{align}

We choose $\overline{U}_2$ of the following form:

(23) \begin{align} \overline{U}_2& = \nu _2-\widehat{a}_1x_4x_6 -\widehat{a}_3x_4\omega _r +\ddot{\phi }_d + k_{x1}\left(e_{x_2}-k_{1}e_{x_1}\right) \end{align}

where $\nu _2$ is a function that depends on $e_{x_1}$ and $e_{x_2}$ which will be determined later. Substituting $\overline{ U}_2$ in (22), we obtain the following expression:

(24) \begin{align} \dot{e}_{x_2} & = -\nu _2 + \widetilde{b}_1\overline{U}_2-\widetilde{a}_1x_4x_6 -\widetilde{a}_3x_4\omega _r \end{align}

The Lyapunov candidate function $V_2$ is defined as:

(25) \begin{equation} V_{2} = \frac{1}{2}\left(e_{x_1}^{2}+e_{x_2}^{2}\right)+\frac{1}{2\gamma }\left(\widetilde{a}_{1}^{2}+\widetilde{a}_{3}^{2} + \widetilde{b}_{1}^{2}\right) \end{equation}

where $\gamma$ is a positive constant. The time derivative of $V_{2}$ is

(26) \begin{align} \dot{V}_{2}& = e_{x_1}\dot{e}_{x_1}+e_{x_2}\dot{e}_{x_2} +\frac{1}{\gamma }\left(\widetilde{a}_{1}\dot{\widetilde{a}}_1+\widetilde{a}_{3}\dot{\widetilde{a}}_3+\widetilde{b}_{1}\dot{\widetilde{b}}_1\right) \nonumber \\ & = e_{x_1}\left(e_{x_2}-k_{x1}e_{x_1}\right) - e_{x_2}\nu _2 + \widetilde{b}_1\left(\overline{U}_2e_{x_2}+\frac{1}{\gamma }\dot{\widetilde{b}}_1\right) \nonumber \\ & \quad -\widetilde{a}_1\left(e_{x_2}x_4x_6-\frac{1}{\gamma }\dot{\widetilde{a}}_1\right)-\widetilde{a}_3\left(e_{x_2}x_4\omega _r-\frac{1}{\gamma }\dot{\widetilde{a}}_3\right) \end{align}

To satisfy the Lyapunov condition $\dot{V}_2\leq 0$ , we obtain the following conditions:

(27) \begin{equation} \begin{cases} \nu _2 = e_{x_1}+k_{2}e_{x_2} \,\, \mbox{with} \,\, k_{2}\gt 0\\ \dot{\widetilde{b}}_1 = -\gamma \overline{U}_2e_{x_2}\\ \dot{\widetilde{a}}_1 = \gamma x_4x_6 e_{x_2}\\ \dot{\widetilde{a}}_3 = \gamma x_4\omega _r e_{x_2} \end{cases} \end{equation}

Replacing the expression of $\nu _2$ from (27) in (23), we obtain the following command expression:

(28) \begin{align} \overline{U}_2& = \ddot{\phi }_d -\widehat{a}_1x_4x_6-\widehat{a}_3x_4\omega _r + e_{x_1}+k_{2}e_{x_2}+k_{1}\left(e_{x_2}-k_{1}e_{x_1}\right) \end{align}

and the derivative of $V_2$ becomes

(29) \begin{align} \dot{V}_{2}=-k_{1}e_{x_1}^{2}-k_{2}e_{x_2}^{2} \end{align}

By using Lyapunov theory and LaSalle theorem, we conclude that the orientation error dynamics is asymptotically stable, and that the orientation angle $\phi$ of the quadrotor follows asymptotically its desired angle.

5.2. Control movement around $y$ and $z$ axis

Following the same procedure, the controller of attitude $\theta$ is as follows:

(30) \begin{equation} \begin{cases} \overline{U}_3 = \ddot{\theta }_d -\widehat{a}_4x_2x_6 -\widehat{a}_6x_2\omega _r + e_{x3}+k_{4}e_{x_4}+k_{3}\left(e_{x_4}-k_{3}e_{x_3}\right) \\ \dot{\widetilde{b}}_2 = -\gamma \overline{U}_3e_{x_3}\\ \dot{\widetilde{a}}_4 = \gamma x_2x_6 e_{x_4} \\ \dot{\widetilde{a}}_6 = \gamma x_2\omega _r e_{x_4} \end{cases} \end{equation}

The controller of attitude $\psi$ is as follows:

(31) \begin{equation} \begin{cases} \overline{U}_4 = \ddot{\psi }_d -\widehat{a}_7x_2x_4 + e_{x_5} +k_{5}\left(e_{x_6}-k_{5}e_{x_5}\right)+k_{6}e_{x_6} \\ \dot{\widetilde{b}}_3 = -\gamma \overline{U}_4e_{x_6}\\ \dot{\widetilde{a}}_7 = \gamma x_2x_4 e_{x_6} \end{cases} \end{equation}

The controller of altitude $z$ is as follows:

(32) \begin{equation} \begin{cases} \overline{U}_1 = e_{x_7}-k_{x8}e_{x_{8}}+\ddot{z}_d+g + k_{7}\left(e_{x_{8}}-k_{7}e_{x_7}\right) \\ \dot{\widetilde{b}}_z = \gamma \overline{U}_xe_{x_{8}}\\ b_z =b_4 cx_1 cx_3 \end{cases} \end{equation}

The controller in the $x$ direction can be expressed as follows:

(33) \begin{equation} \begin{cases} \overline{U}_x = e_{x_9}-k_{10}e_{x_{10}}+\ddot{x}_d + k_{9}\left(e_{x_{10}}-k_{9}e_{x_9}\right) \\ \dot{\widetilde{b}}_x = \gamma \overline{U}_xe_{x_{10}}\\ b_x = b_4 U_1 \end{cases} \end{equation}

The controller in the $y$ direction can be expressed as follows:

(34) \begin{equation} \begin{cases} \overline{U}_y = e_{x_{11}}-k_{12}e_{x_{12}}+\ddot{y}_d + k_{11}\left(e_{x_{12}}-k_{11}e_{x_{11}}\right)\\ \dot{\widetilde{b}}_y = \gamma \overline{U}_ye_{x_{12}}\\ b_y = b_4 U_1 \end{cases} \end{equation}

6. Control design of BLF

In this section, control design and analysis are presented for the quadrotor model. We consider the case where the system model is known. The control design is based on backstepping, employing a BLF candidate in the first step and quadratic Lyapunov function candidates in the remaining steps.

6.1. BLF PRELIMINARIES:

To better establish constraint compensation and performance boundaries, and to make the paper more self-contained, we recall here the following definitions, assumptions, and lemma.

Definition [26]: A BLF is a scalar function $V(x)$ , defined with respect to the system $\dot x = f (x, t)$ on an open region $\Omega$ containing the origin, that is, continuous, positive definite, has continuous first-order partial derivatives at every point of $\Omega$ , has the property $V(x)\rightarrow +\infty$ as $x$ approaches the boundary of $\Omega$ , and satisfies $V(x(t)) \leq b, \forall t \geq 0$ , along the solution of $\dot x = f (x, t)$ for $x(0)\in \Omega$ and some positive constant $b$ .

A BLF should be symmetric or asymmetric according to the boundary character, which is shown in Fig. 3.

Figure 3. Symmetric (left) and asymmetric (right) barrier Lyapunov function.

Assumption 1. [35]: For any $K_{c_1}\gt 0,$ there exists positive constants $Y_1, Y_2$ , $\underline{Y}_0$ and $ \overline{Y}_0$ satisfying $\max \{\underline{Y}_0, \overline{Y}_0\} \lt k_{c_1}$ , such that the desired trajectory $y_d(t)$ and its time derivatives satisfy $-\underline{Y}_0\leq y_d(t)\leq \overline{Y}_0$ and $|\dot{y}_d(t)|\leq Y_1,\,\, |\ddot{y}_d(t)|\leq Y_2,\,\, \forall t\geq 0$ .

Lemma 1. [36]: For any positive constants $k_{a_1}, k_{b_1}$ let $\chi _1=\{x_1\in \mathbf{R},\,\, k_{a_1}\lt x_1\lt k_{b_1}\}\subset \mathbf{R}$ and $\mathbf{N}:=\mathbf{R}^l\times \chi _1\subset \mathbf{R}^{l+1}$ be open sets. Consider the system

\begin{equation*}\dot \eta =h(\eta,t)\end{equation*}

where $\eta =[\omega,x_1]\in \mathbf{N}$ , and $h\,:\, \mathbf{R}_+ \times \mathbf{N} \rightarrow \mathbf{R}^{l+1}$ is piecewise continuous with respect to $t$ and locally Lipschitz with respect to $\eta$ , uniformly with respect to $t$ , on $\mathbf{R}_+ \times \mathbf{N}$ . Suppose that there exist functions $U: \mathbf{R}^l \times \mathbf{R}_+$ and $V_1: \chi _1 \rightarrow \mathbf{R}_+$ continuously differentiable and positive definite in their respective domains, such that

\begin{equation*}V_1(x_1)\rightarrow \infty \,\, as \,\, x_1\rightarrow -k_{a_1}\,\, or \,\, x_1\rightarrow k_{b_1},\end{equation*}

\begin{equation*}\gamma _1(\parallel \omega \parallel )\leq U(\omega )\leq \gamma _2(\parallel \omega \parallel )\end{equation*}

where $\gamma _1$ and $\gamma _2$ are class $\mathbb{K}^\infty$ functions. Let

\begin{equation*}V(\eta )=V_1(x_1)+U(\omega )\end{equation*}

and $x_1(0) \in \chi _1$ . If the inequality

\begin{equation*}\dot V=\frac {\partial V}{\partial \eta }h \leq -cV+b\end{equation*}

holds in the set $\eta \in \mathbf{N}$ and $c$ , $b$ are positive constants, then $x_1(t)$ remains in the open set $\chi _1, \forall t\in [0,\infty [$ .

6.2. Control design with a BLF

$\bullet$ First iteration: Denote $e_{x_1}=x_1-\phi _d$ , $e_{x_3}=x_3-\theta _d$ , $e_{x_5}=x_5-\psi _d$ , $e_{x_7}=x_7-z_d$ , $e_{x_9}=x_9-x_d$ , $e_{x_{11}}=x_{11}-y_d$ as the tracking errors, and $e_{x_{2i}}=x_{2i}-\alpha _{2i}$ , $\forall i=1,\ldots,6$ as virtual errors, where $\alpha _{2i}$ is a virtual control function. By selecting the symmetric BLF

(35) \begin{align} V_1& = \frac{1}{2}\log \left(\frac{\varepsilon _{1}^{2}}{\varepsilon _{1}^{2}-e_{x_1}^{2}}\right) +\frac{1}{2}\log \left(\frac{\varepsilon _{3}^{2}}{\varepsilon _{3}^{2}-e_{x_3}^{2}}\right) \nonumber \\ & \quad + \frac{1}{2}\log \left(\frac{\varepsilon _{5}^{2}}{\varepsilon _{5}^{2}-e_{x_5}^{2}}\right) +\frac{1}{2}\log \left(\frac{\varepsilon _{7}^{2}}{\varepsilon _{7}^{2}-e_{x_7}^{2}}\right) \nonumber \\ & \quad + \frac{1}{2}\log \left(\frac{\varepsilon _{9}^{2}}{\varepsilon _{9}^{2}-e_{x_9}^{2}}\right) +\frac{1}{2}\log \left(\frac{\varepsilon _{11}^{2}}{\varepsilon _{11}^{2}-e_{x_{11}}^{2}}\right) \end{align}

where $-\varepsilon _{2i-1}$ and $\varepsilon _{2i-1}$ are the symmetric constraints on $e_{x_{2i-1}}$ , and each of these can be set independently, depending on the upper and lower bounds of $\phi _d,\theta _d,\psi _d,z_d,x_d$ , and $y_d$ . By taking the derivative of Eq. (35), we get the following expression:

(36) \begin{align} \dot{V}_{1}& = \frac{e_{x_1}\dot{e}_{x_1}}{\varepsilon _{1}^{2}-e_{x1}^{2}}+\frac{e_{x_3}\dot{e}_{x_3}}{\varepsilon _{3}^{2}-e_{x_3}^{2}}+\frac{z_{5}\dot{z}_{5}}{\varepsilon _{5}^{2}-e_{x_5}^{2}}\nonumber \\ &\quad +\frac{z_{7}\dot{e}_{x_7}}{\varepsilon _{7}^{2}-e_{x_7}^{2}}+\frac{e_{x_9}\dot{e}_{x_9}}{\varepsilon _{9}^{2}-e_{x_9}^{2}}+\frac{e_{x_{11}}\dot{e}_{x_{11}}}{\varepsilon _{11}^{2}-e_{x11}^{2}} \end{align}

So, we can have

(37) \begin{align} \dot{V}_{1}& = \frac{e_{x_1}(e_{x_2}+\alpha _{2}-\dot{\phi }_{d})}{\varepsilon _{1}^{2}-e_{x_1}^{2}} +\frac{e_{x_3}(e_{x_4}+\alpha _{4}-\dot{\theta }_{d})}{\varepsilon _{3}^{2}-e_{x_3}^{2}}\nonumber \\ & \quad + \frac{e_{x_5}(e_{x_6}+\alpha _{6}-\dot{\psi }_{d})}{\varepsilon _{5}^{2}-e_{x_5}^{2}}+\frac{e_{x_7}(e_{x_8}+\alpha _{8}-\dot{z}_{d})}{\varepsilon _{7}^{2}-e_{x_7}^{2}}\nonumber \\ & \quad + \frac{e_{x_9}(e_{x_{10}}+\alpha _{10}-\dot{x}_{d})}{\varepsilon _{9}^{2}-e_{x_9}^{2}}+\frac{e_{x_{11}}(e_{x_{12}}+\alpha _{12}-\dot{y}_{d})}{\varepsilon _{11}^{2}-e_{x_{11}}^{2}} \end{align}

Design the virtual controllers $\alpha _{2i}$ as

(38) \begin{align} \alpha _{2}&=-k_{1}\left(\varepsilon _{1}^{2}-e_{x_1}^{2}\right)e_{x_1}+\dot{\phi }_{d} \nonumber \\ \alpha _{4}&=-k_{3}\left(\varepsilon _{3}^{2}-e_{x_3}^{2}\right)e_{x_3}+\dot{\theta }_{d} \nonumber \\ \alpha _{6}&=-k_{5}\left(\varepsilon _{5}^{2}-e_{x_5}^{2}\right)e_{x_5}+\dot{\psi }_{d} \nonumber \\ \alpha _{8}&=-k_7\left(\varepsilon _{7}^{2}-e_{x_7}^{2}\right)e_{x_7}+\dot{z}_{d} \\ \alpha _{10}&=-k_9\left(\varepsilon _{9}^{2}-e_{x_9}^{2}\right)e_{x_9}+\dot{x}_{d} \nonumber \\ \alpha _{12}&=-k_{11}\left(\varepsilon _{11}^{2}-e_{x_{11}}^{2}\right)e_{11}+\dot{y}_{d} \nonumber \end{align}

where $k_{2i-1}$ are positive constants. Equation (38) is substituted into (37) to obtain

(39) \begin{align} \dot{V}_{1}& = -k_{1}e_{x_1}^{2}-k_{3}e_{x_3}^{2}-k_{5}e_{x_5}^{2}-k_{7}e_{x_7}^{2}-k_{9}e_{x_9}^{2}-k_{11}e_{x_{11}}^{2}\nonumber \\ & \quad +\frac{e_{x_1}e_{x_2}}{\varepsilon _{1}^{2}-e_{x_1}^{2}}+\frac{e_{x_3}e_{x_4}}{\varepsilon _{3}^{2}-e_{x_3}^{2}} +\frac{e_{x_5}e_{x_6}}{\varepsilon _{5}^{2}-e_{x_5}^{2}} \nonumber \\ & \quad + \frac{e_{x_7}e_{x_8}}{\varepsilon _{7}^{2}-e_{x_7}^{2}} +\frac{e_{x_9}e_{x_{10}}}{\varepsilon _{9}^{2}-e_{x_9}^{2}}+\frac{e_{11}e_{x_{12}}}{\varepsilon _{11}^{2}-e_{x_{11}}^{2}} \end{align}

When $e_{x_{2i}}\rightarrow 0$ , there is

(40) \begin{align} \dot{V}_1& = -\sum _{i=1}^{6}k_{2i-1}e_{x_{2i-1}}^{2}\leq 0,\,\,\, \forall t\geq 0 \end{align}

Then, $V_1\rightarrow 0$ when $t\rightarrow \infty$ . The second term of $\dot{V}_1$ can be canceled in a later step.

$\bullet$ Second iteration: Define the Lyapunov function candidate as

(41) \begin{align} V_{2}& = V_{1}+\frac{1}{2}\sum _{i=1}^{6}e_{x_{2i}}^{2} \end{align}

Differentiating Eq. (41), we can obtain

(42) \begin{align} \dot{V}_{2}& = \dot{V}_1+e_{x_2}\left (b_1U_2+a_1(e_{x_4}+\alpha _{4})(e_{x_6}+\alpha _{6}) \right. \nonumber \\ & \quad \left. +a_3(e_{x_4}+\alpha _{4})\omega _r - \dot{\alpha }_{2}\right ) + e_{x_4}\left (b_2U_3 \right. \nonumber \\ & \quad\left. + a_4(e_{x_2}+\alpha _{2})(e_{x_6}+\alpha _{6}) +a_6(e_{x_2}+\alpha _{2})\omega _r-\dot{\alpha }_{4}\right ) \nonumber \\ & \quad+e_{x_6}(b_3U_4+a_7(e_{x_2}+\alpha _{2})(e_{x_4}+\alpha _{4})-\dot{\alpha }_{6})\nonumber \\ &\quad +e_{x_8}\left(g-\frac{U_1}{m} c(e_{x_1}+\phi _d)c(e_{x_3}+\theta _d)-\dot{\alpha }_{8}\right)\nonumber \\ &\quad +e_{x_{10}}\left(\frac{U_1}{m}U_x-\dot{\alpha }_{10}\right) +e_{x_{12}}\left(\frac{U_1}{m}U_y-\dot{\alpha }_{12}\right) \end{align}

The control law can be written as

(43) \begin{align} U_1& = \frac{m}{c x_1 cx_3}\left[k_{8}e_{x_8}-\dot{\alpha }_{8}+\frac{e_{x_7}}{\varepsilon _{7}^{2}-e_{x_7}^{2}}\right]\nonumber \\ U_2& = \frac{1}{b_1}\left[-k_{2}e_{x_2}+\dot{\alpha }_{2}-\frac{e_{x_1}}{\varepsilon _{1}^{2}-e_{x_1}^{2}}-a_1(e_{x_4}+\alpha _{4})(e_{x_6}+\alpha _{6}) -a_3(e_{x_4}+\alpha _{4})\omega _r\right]\nonumber \\ U_3& = \frac{1}{b_2}\left[-k_{4}e_{x_4}+\dot{\alpha }_{4}-\frac{e_{x_3}}{\varepsilon _{3}^{2}-e_{x_3}^{2}} -a_4(e_{x_2}+\alpha _{2})(e_{x_6}+\alpha _{6}) -a_6(e_{x_2}+\alpha _{2})\omega _r\right] \nonumber \\ U_4& = \frac{1}{b_3}\left[-k_{6}e_{x_6}+\dot{\alpha }_{6}-\frac{e_{x_5}}{\varepsilon _{5}^{2}-e_{x_5}^{2}}-a_7(e_{x_2}+\alpha _{2})(e_{x_4}+\alpha _{4} )\right] \nonumber \\ U_x& = \frac{m}{U_1}\left[-k_{10}e_{x_{10}}+\dot{\alpha }_{10}-\frac{e_{x_9}}{\varepsilon _{9}^{2}-e_{x_9}^{2}}\right]\nonumber \\ U_y& = \frac{m}{U_1}\left[-k_{12}e_{x_{12}}+\dot{\alpha }_{12}-\frac{e_{11}}{\varepsilon _{11}^{2}-e_{x_{11}}^{2}}\right] \end{align}

where $k_{2i}$ are positive constants. Equation (43) is substituted into (42) to obtain

(44) \begin{align} \dot{V}_{2}& = -\sum _{i=1}^{12}k_ie_{x_i}^{2}\leq 0 \end{align}

By using Barbalat’s lemma, we conclude that the error $e_{x_i}$ is asymptotically stable, and that the trajectory of the quadrotor follows asymptotically its desired trajectory.

Remark 1.

1. 1. $U_1=\sum _{i=1}^{4}f_i$ with $f_i=b\omega _{i}^{2}.$ $b\gt 0$ is a given constant and $w_i$ the angular speed resulting of motor $i$ . The system input is such $U_1\gt 0.$

2. 2. The axillary control inputs $U_x$ and $U_y$ are used to calculate the desired roll $\phi _d$ and desired pitch $\theta _d$ , then, the roll and pitch will be controlled in $U_2,U_3$ . The desired roll and pitch are found as [42], [43]:

   (45) \begin{align} \phi _d& = \arcsin (U_x\sin \psi _d-U_y\cos \psi _d)\nonumber \\ \theta _d& = \arcsin \left(\frac{U_x\sin \psi _d-U_y\cos \psi _d}{\cos \phi _d}\right) \end{align}

Proposition 1. Consider the closed loop system ( 4 ), ( 38 ) and ( 43 ) under Assumption 1. If the initial conditions are such that $e_{x_{i}}(0)\in \left\{ e_{x_{i}}\in \mathbf{R}^{12};\, |e_{x_{2i-1}}|\lt \varepsilon _{2i-1}\right\}$ , then the following properties hold.

1. 1. The signals $e_{x_i}(t)$ remain in the compact set defined by

   (46) \begin{align} \Omega:=\left\{e_{x_{i}}\in \mathbf{R}^{12}, |e_{x_{2i-1}}| \lt D_{e_{x}} ; \,\,\|e_{x_{2i}}\|\leq \sqrt{2V_2(0)}\right\} \end{align}
   (47) \begin{align} D_{e_{x}}=\varepsilon _{2i-1}\sqrt{1-e^{-2V_2(0)}} \end{align}
   where $V_2$ is the overall Lyapunov function candidate obtained from ( 35 ) and ( 41 ).

2. 2. The output tracking error $e_{x_{2i-1}}(t)$ converges to zero asymptotically,

3. 3. All closed loop signals are bounded.

Proof.

• • From $\dot{V}_2\leq 0$ , it follows that $V_2(t)\leq V_2(0)$

  \begin{equation*}\frac {1}{2}\log \left(\frac {\varepsilon _{2i-1}^{2}}{\varepsilon _{2i-1}^{2}-e_{x_{2i-1}}^{2}}\right)\leq V_2(0)\end{equation*}
  Applying the exponential function, the following inequality can be obtained:
  \begin{equation*}|e_{x_{2i-1}}(t)|\leq \varepsilon _{2i-1}\sqrt {1-e^{-2V_2(0)}}\end{equation*}
  Similarly, from the fact that $\frac{1}{2}\sum _{i=1}^{6}e_{x_{2i}}^{2}\leq V_2(0)$ , we can show that $\|e_{x_{2i}}\|\leq \sqrt{V_2(0)}$ .

• • As can be seen in Eq. (42), $\dot{V}_2$ is negative definite when the control law is (43). Thus, the system is asymptotically stable. When $t\rightarrow \infty$ , $e_{x_i} \rightarrow 0$ . Based on Lemma 1, $e_{x_{2i-1}}$ remains in $]-\varepsilon _{2i-1},\varepsilon _{2i-1}[, \forall t\in [0,\infty )$ . Therefore, the constraints will never be violated.

• • Since $x_{2i-1}(t)=e_{x_{2i-1}}(t)+\eta _d(t)$ , where $\eta _d=(\phi _d,\theta _d,\psi _d,x_d,y_d,z_d)^T,$ $ |e_{x_{2i-1}}(t)|\leq \varepsilon _{2i-1}$ and $\underline{Y}_0(t)\leq \eta _d(t)\leq \overline{Y}_0(t)$ , then,

  \begin{equation*}|x_{2i-1}| \leq \varepsilon _{2i-1}+\overline {Y}_0(t)\end{equation*}
  The boundedness of $e_{x_{2i-1}}$ and $\eta _d$ implies that the state $x_{2i-1}$ is bounded. Together with the fact that $\dot{\eta }_d(t)$ is bounded from Assumption 1, it is clear, from (38), that the stabilizing functions $\alpha _{2i}$ are also bounded.

7. Simulation results

In this subsection, numerical simulations are performed to illustrate the theoretical results of the previous sections. The system’s parameter values (IS units) used for simulations are given in Table I.

Table I. System’s parameters.

In the simulation, the desired trajectories for the quadrotor are

(48) \begin{equation} \begin{cases} x_d=0\,\,\,\,\,\,\,\, if \,\, t\lt 3 \,\, s\\ x_d=0.5\sin (0.1t)\,\,\,\,\,\,\,\, if \,\, 3\lt t\lt 30 \,\, s\\ \end{cases} \end{equation}

(49) \begin{equation} \begin{cases} y_d=0\,\,\,\,\,\,\,\, if \,\, t\lt 3 \,\, s\\ y_d=0.5\sin (0.1t)+0.5\,\,\,\,\,\,\,\, if \,\, 3\lt t\lt 30 \,\, s\\ \end{cases} \end{equation}

\begin{equation*}z_d=1,\,\,\, \psi _d=0\end{equation*}

and all the initial values are set to zero.

7.1. Simulation results: Backstepping

In this section, we will present the simulation results of the backstepping method. In this implementation, the controllers parameters are chosen by trial and error as follows:

\begin{equation*}\begin {array}{ll} k_1=55, k_2=45, k_3=52& k_4=47, k_5=2, k_6=1.9\\ k_7=2.7, k_8=1.4, k_9=2.8&k_{10}=1.3, k_{11}=2.5, k_{12}=2.5\\ \end {array}\end{equation*}

Figure 4 shows the real and desired positions for the backstepping method. Figure 5 displays the control signals. Figure 6 represents the real and reference trajectories of the quadrotor in 3D. Figures 7 and 8 show the time evolution of the tracking positions and attitude errors.

Figure 4. Real and reference positions (backstepping).

Figure 5. Control heave force, roll, pitch, and yaw torque (backstepping).

Figure 6. Real and reference trajectories of quadrotor in 3D (backstepping).

Figure 7. Position tracking errors (backstepping).

Figure 8. Attitude tracking errors (backstepping).

7.2. Simulation results: Adaptive backstepping

In this section, we will present the simulation results of the adaptive backstepping method. The parameters are chosen as follows:

\begin{equation*}\begin {array}{ll} k_1=211, k_2=179,& k_3=211, k_4=179\\ k_5=201, k_6=171,& k_7=2.7, k_8=1.4\\ k_9=291, k_{10}=131,& k_{11}=57, k_{12}=45\\ \end {array}\end{equation*}

Similarly, Figs. 9 –13 show the real and desired positions, control signals, 3D trajectories, position errors, and attitude errors, respectively, for the adaptive backstepping method.

Figure 9. Real and reference position (adaptive).

Figure 10. Control heave force, roll, pitch, and yaw torques (adaptive).

Figure 11. Real and reference trajectories of quadrotor in 3D (adaptive).

Figure 12. Position tracking errors (adaptive).

Figure 13. Attitude tracking errors (adaptive).

7.3. Simulation results: BLF control

In this section, we will present the simulation results of the BLF backstepping method. In the implementation, the parameters are chosen as follows:

\begin{equation*}\begin {array}{l@{\quad}l@{\quad}l} \varepsilon _1=\dfrac {\pi }{10}&\varepsilon _2=\dfrac {\pi }{10}&\varepsilon _3=\pi \\[3pt] \varepsilon _4=3&\varepsilon _5=5&\varepsilon _6=5\\ \end {array}\end{equation*}

\begin{equation*}\begin {array}{l@{\quad}l} k_1=k_3=k_5=0.1&k_7=k_9=k_{11}=1\\ [3pt]k_2=k_4=k_6=10&k_8=k_{10}=k_{12}=10\\ \end {array}\end{equation*}

Figures 14 –18 show the time evolution of the real and desired positions, control signals, 3D trajectories, position errors, and attitude errors, respectively, for the BLF method.

Figure 14. Real and reference positions (BLF).

Figure 15. Control heave force, roll, pitch, and yaw torques (BLF).

Figure 16. Real and reference trajectories of quadrotor in 3D (BLF).

Figure 17. Position errors (BLF).

7.4. Results and comments

Figures 6, 11, and 16 show that the quadrotor is able to reach the desired position using the different methods and strategies. Figures 7, 8, 12, 13, 17, and 18 show the errors between the real and reference state variables. The position and attitude errors converge within a small proximity to zero.

Figures 5, 10, and 15 show the control forces $U_1$ and control torques $U_2$ , $U_3$ , and $U_4$ needed for tracking stabilization of these different methods. Figure 10 indicates that these control inputs perform sufficiently in driving the quadrotor along the desired trajectory under internal parameter uncertainties. The tracking performance of the classical backstepping-based controller is quite similar to that of the adaptive controller. However, less force is required in the BLF control method to ensure trajectory tracking. For example, for heave force to stabilize the quadrotor, 0.6 $N$ is needed using BLF control, as compared to 0.8 $N$ for the adaptive controller and 0.9 $N$ for the classical backstepping controller. Figures 4, 9, and 14 show the time evolution of the state variables $x$ , $y$ , and $z$ . We can see that the inertial position converges to the desired positions $x_d$ , $y_d$ , and $z_d$ . Using the adaptive or BLF method, we can eliminate the error at position $x$ in the classical backstepping method. From the simulation results, it is evident that asymptotic tracking performance was achieved. Figures 4, 9, and 14 show that the quadrotor can track the desired trajectory using backstepping or BLF methods. Using BLF, the trajectory of the quadrotor converges to the desired value within 5 seconds.

Figures 17 and 18 illustrate that the BLF controller offered good tracking performance, with the actual positions and angles of the quadrotor always within the constraint range (0 0.5), and tracking error falling strictly within the error constraints. Figures 17 and 18 indicate that the constraints are never violated. The response time of the step response (Figs. 4, 9, and 14) is less than 5 seconds, and the error approaches zero asymptotically. Table II show the values of the performance criteria of the root mean square errors between real and reference states. For example, in the 2nd column of the table, the quadratic error in the $x$ -position is $78e^{-4}m$ for the backstepping method, $8.6405e^{-4}m$ for the adaptive method, and $8.4716e^{-4}m$ for the BLF method. At the 5th column, the root mean square errors of the roll angle $\phi$ for the classical backstepping method are $0.3445 deg$ , and for the adaptive method the error is $1.6405e^{-4}deg$ , while the error for the BLF method is $5.227e^{-28} deg$ . The BLF control strategy achieves the requirement: it has a smaller root mean square errors and no boundary violations. However, with the backstepping and adaptive controllers, the boundary of the position and angle tracking errors are violated.

Table II. Root mean square errors between real and reference states.

Figure 18. Attitude errors (BLF).

For example, in contrast to the control effects of backstepping controller (Fig. 8) and adaptive controller (Fig. 13), the BLF controller (Fig. 18) will not cause oscillation of the pitch angle.

Due to the constraint boundary, the BLF controller is more complex as well as computationally more complicated than traditional backstepping or adaptive methods. However, when the constraint boundary is taken into account, the motion position and angle range of the quadrotor are limited to a reasonable range.

8. Conclusion

Backstepping techniques are employed in this paper to develop a tracking controller for an underactuated quadrotor UAV. Two backstepping control approaches, a classical and an adaptative controller, are presented, such that the tracking error converges to zero. A BLF-based backstepping control scheme is proposed for quadrotors with full-state constraints. The BLF is constructed in each step of the backstepping design to guarantee that the full-state variables of the quadrotor system are constrained. Lyapunov-based analysis is utilized to prove that asymptotic tracking is achieved with minimal assumption of UAV roll and pitch angles. Simulation results are provided to demonstrate the effectiveness of the proposed method. The simulation results show that the BLF control algorithm is more effective than the classical and adaptive backstepping control algorithms. For future research, we will consider the implementation of backstepping methods on a real quadrotor system and design of a robust controller to counteract parametric uncertainties and disturbances.