Design and Implementation of Integral Backstepping Sliding Mode Control for Quadrotor Trajectory Tracking

Abstract

A robust trajectory tracking control scheme for quadrotor unmanned aircraft vehicles under uncertainties is proposed herein. A tracking controller combined with the sliding mode and integral backstepping is performed for position and attitude tracking. The stability of the trajectory tracking controller of the quadrotor is investigated via Lyapunov stability analysis. By incorporating force and torque disturbances into numerical simulations, the results demonstrate the effectiveness of the proposed quadrotor trajectory controller. Finally, the experiments validate the feasibility of the proposed controller.

1. Introduction

Recently, the rapid development of unmanned aircraft vehicles (UAVs) has received significant attention in various applications, such as search-and-rescue, safety inspections for civil buildings, agriculture, and forest surveillance. Different types of UAVs exist, including those that can take-off and land vertically in limited spaces, as well as those that can easily hover over the desired target. Compared with a helicopter, a quadrotor is typically constructed using a symmetrical lightweight airframe, with four brushless direction current (BLDC) motors mounted on it, and four fixed-pitch propellers attached to the respective motors. A quadrotor offers better mobility in attitude change and path planning compared with helicopters. However, the design of the flight controller of a quadrotor is associated with at least two challenges. First, the quadrotor is a multiple-input multiple-output (MIMO) unstable nonlinear system. Second, quadrotors, similar to other types of UAVs, are always subjected to external and internal disturbances, model uncertainties, and parametric perturbation. Therefore, they are difficult to control accurately in terms of trajectory tracking. Hence, the design and implementation of a stable and robust quadrotor flight controller are essential [1].

Proportional–integral–derivative (PID) control is a typical control technique that is typically used in industry. PID control involves many iterative trials for adjusting the parameters of the PID controller [2]. In [3], a nonlinear PID-type controller for the total thrust and torques was applied to perform trajectory tracking and the authors derived the adjustment guidelines. Meanwhile, in terms of the design of the flight controller of a V-tail quadrotor, Castillo–Zamora et al. performed the comparisons of tracking capability among a proportional-derivative (PD) controller, PID controller, and sliding mode controller (SMC) [4]. In their study, the simulation results showed that the PD controller presented an error in the stable state that would not disappear. This can be problematic if an application with precision is required. Meanwhile, the SMC has a facile and rapid response that drives the quadrotor to converge to the desired point in space.

Model-based controllers have been shown to exhibit excellent performance in resolving nonlinear control system problems. Observer-based backstepping and SMC for the trajectory tracking of a quadrotor have been presented in many studies. In [5], the authors designed an extended state observer (ESO) to estimate unmeasurable states and external disturbances. Wind gust and actuator faults can be considered as lumped disturbances, the disturbance torque caused by lumped disturbances happen suddenly, and cannot be estimated by traditional ESO. In order to alleviate lumped disturbances, high-order ESO have been applied. Nevertheless, higher observer orders will lead to higher control gains with fixed bandwidth, which will in return excite the sensor noise and introduce them into the control loop. Shi D. et al. proposed the super-twisting based ESO (ST-ESO) to estimate and attenuate the impact of lumped disturbances in finite time. In the same study, the super-twisting based SMC (ST-SMC) is also designed as the feedback controller to drive the attitude angle and angular velocity to their desired value in finite time [6]. In [7], the nonlinear disturbance observer (NDO) was used to compensate the disturbances and model-system errors. A backstepping SMC and observer-based fault estimation for a quadrotor was proposed in [8], where a regular SMC was used for attitude control, and the backstepping SMC was utilized for position and yaw angle control. An observer-based fault estimation serves as an alarming module as well as provides an estimation of faults occurring during the take-off and hovering modes. In the same paper, a stability analysis of the overall system based on the Lyapunov stability theory was presented.

Actuator dynamics, including transient responses and control input constraints, directly affect and limit the torque inputs to the propeller rotation of a quadrotor. It is risky to design the controller when the actuator saturation situation is omitted. By incorporating the SMC and dynamic surface controller (DSC) into a backstepping-like framework, the authors of [9,10] proposed a backpropagating constraint-based controller to solve the constrained actuator dynamics problem.

The contribution of this study is the design and implementation of an integral backstepping sliding mode controller (IB-SMC) of a quadrotor for real-time trajectory tracking [11,12]. For the hardware design of the quadrotor, we used a PX4 module as the flight control board; it contains six degree-of-freedom (DOF) inertial measurement unit (IMU) sensors, i.e., a three-axis accelerometer and a three-axis gyroscope [13]. A global positioning system (GPS) was used as the position estimator of the quadrotor. A radio frequency (RF) module was utilized to receive the control commands of the quadrotor from the user. The remainder of this paper is organized as follows: Section 2 presents the kinematics and dynamical model of the quadrotor. In Section 3, the design and stability analysis of the IB-SMC for the quadrotor are discussed. Section 4 provides the simulations and experimental results of the IB-SMC of the quadrotor. The conclusions are presented in Section 5.

2. Quadrotor Model

In navigation, guidance, and control of an aircraft or rotorcraft applications, several frame systems are used in design and analysis. These frame systems are the geodetic frame system, the earth-centered earth-fixed (ECEF) frame system, the local north-east-down (NED) frame system, the vehicle-carried NED frame system, and the body-fixed frame system. The geodetic frame system is widely used in GPS-based navigation. A coordinate point near the earth’s surface is characterized by longitude, latitude, and height. The ECEF frame system rotates with the earth around its spin axis, the origin of the ECEF frame system is located at the center of the earth. The local NED frame system is known as a navigation or ground coordinate system. It is a coordinate frame fixed to the earth’s surface and the origin of the local NED frame system is arbitrarily fixed to a point on the earth’s surface. The local NED frame plays an important role in flight control and navigation. Navigation of small-scale UAVs rotorcraft is normally carried out with this frame. The vehicle-carried NED frame system is associated with the flying vehicle and the origin of the vehicle-carried NED frame system is located at the center of gravity of the flying vehicle. The body-fixed frame system is vehicle-carried and is directly defined on the body of the flying vehicle. The origin of the body-fixed frame system is located at the center of gravity of the flying vehicle. The six DOFs IMU sensor is intensively used to capturing the rigid-body dynamics of unmanned systems [14]. Due to the inherent mechanical design and power limitation of quadrotor, the applications of quadrotor are normally operated at low speeds in small regions, so the local NED frame and the body-fixed frame are considered in this article.

The structure of the quadrotor is depicted in Figure 1a; the quadrotor has six DOFs and four control inputs. The six DOFs include translational motion in three directions and rotational motion around the three axes. Four rotors driven by a motor were mounted at the endpoint of two cross-shaped frames [8].

To establish a mathematical model for the quadrotor, it is assumed that the quadrotor frame is a symmetric grid body, the center of mass is at the geometric center of the body, and the flapping dynamics of the frame are negligible. For describing the attitude and position of the quadrotor, the local NED frame, i.e., the inertial frame, and the body-fixed frame are introduced. The inertial frame denoted as $E=(X_E{,}^{}{Y}_E{,}^{}{Z}_E)$ and the body-fixed frame denoted as $B=(X_B{,}^{}{Y}_B{,}^{}{Z}_B)$. The inertial frame is based on the Earth with its origin coinciding with the origin of the body-fixed frame prior to take-off. Frame B is fixed with the quadrotor, and the origin of frame B is at the center of mass of the body.

2.1. Kinematics

The linear position vectors of the quadrotor expressed in the inertial frame E and body-fixed frame B are $P^E={\left[\begin{array}{ccc}{x}^E& y^E& z^E\end{array}\right]}^T$ and $P^B={\left[\begin{array}{ccc}{x}^B& y^B& z^B\end{array}\right]}^T$, respectively. The attitude representation can be described by the quaternions and the Euler angles. The advantages of the Euler angles are an intuitive way to represent the attitude and their physical meanings are quite clear. Based on Euler’s theorem, the rotation of a rigid body around one fixed point can be regarded as the composition of several finite rotations around the fixed point. Recently, quaternions are getting more and more attention due to the wide application of high-performance computer and the rapid development of the aircraft’s attitude control technologies [15].

As shown in Figure 1b, the attitude (i.e., angular position) of the quadrotor defined in frame E is described by three Euler angles ${\eta }^E={\left[\begin{array}{ccc}\varphi & \theta & \psi \end{array}\right]}^T$, where $\varphi$, $\theta$, and $\psi$ denote the angles of roll, pitch, and yaw, respectively. $\varphi$ is the angle between planes $X_E{Y}_E$ and $X_B{Y}_B$, $\theta$ is the angle between planes $Y_E{Z}_E$ and $Y_B{Z}_B$, and $\psi$ is the angle between planes $Z_E{X}_E$ and $Z_B{X}_B$. The angular velocity vector of the quadrotor is expressed as ${\dot{\eta }}^E={\left[\begin{array}{ccc}\dot{\varphi }& \dot{\theta }& \dot{\psi }\end{array}\right]}^T$. It is assumed that all attitude angles are limited between $(-\raisebox{1ex}{$\pi $}\!\left/ \!\raisebox{-1ex}{$2$}\right.,\raisebox{1ex}{$\pi $}\!\left/ \!\raisebox{-1ex}{$2$}\right.)$. The linear position transformed from frame E to frame B is expressed as $P^B=R_E^B{P}^E$, where the rotation matrix $R_E^B$ is written as follows:

$$R_E^B=\left[\begin{array}{ccc}\cos \theta \cos \psi & \cos \theta \sin \psi & -\sin \theta \\ \sin \varphi \sin \theta \cos \psi -\cos \varphi \sin \psi & \sin \varphi \sin \theta \sin \psi +\cos \varphi \cos \psi & \sin \varphi \cos \theta \\ \cos \varphi \sin \theta \cos \psi +\sin \varphi \sin \psi & \cos \varphi \sin \theta \sin \psi -\sin \varphi \cos \psi & \cos \varphi \cos \theta \end{array}\right]$$

(1)

If $R_E^B$ is nonsingular, then the rotation matrix $R_E^B$ is an orthogonal matrix, and the inverse matrix ${(R_E^B)}^{-1}$ is obtained via a simple transpose, i.e., ${(R_E^B)}^{-1}={(R_E^B)}^T$. Hence, the linear position is transformed from frame B to frame E, expressed as $P^E={(R_E^B)}^T{P}^B=R_B^E{P}^B$.

Three-axis gyroscopes are generally used in the aerospace industry to measure body rotations. The relationship between time rates of change of the Euler angles and the body axis angular rates are required. Assume that the angular velocity of the quadrotor with respect to frame B is ${\omega }^B={\left[\begin{array}{ccc}p& q& r\end{array}\right]}^T$, where $p$, $q$, and $r$ are the angular velocities of the roll, pitch, and yaw, respectively. Therefore, the transformation matrix $W_E^B$ for the angular velocity from frame E to frame B is written as

$${\omega }^B=\left[\begin{array}{c}p\\ q\\ r\end{array}\right]=W_E^B{\dot{\eta }}^E=\left[\begin{array}{ccc}1& 0& -\sin \theta \\ 0& \cos \varphi & \cos \theta \sin \varphi \\ 0& -\sin \varphi & \cos \theta \cos \varphi \end{array}\right]\left[\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right]$$

(2)

The angular velocity from the body-fixed frame to the inertial frame is expressed as [16,17,18,19].

$${\dot{\eta }}^E=\left[\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right]={\left(W_E^B\right)}^{-1}{\omega }^B=\left[\begin{array}{ccc}1& \sin \varphi \tan \theta & \cos \varphi \tan \theta \\ 0& \cos \varphi & -\sin \varphi \\ 0& \frac{\sin \varphi }{\cos \theta }& \frac{\cos \varphi }{\cos \theta }\end{array}\right]\left[\begin{array}{c}p\\ q\\ r\end{array}\right]$$

(3)

When the quadrotor executes many angular motions with low amplitude, ${\omega }^B$ can be assimilated to ${\left[\begin{array}{ccc}\dot{\varphi }& \dot{\theta }& \dot{\psi }\end{array}\right]}^T$.

2.2. Dynamics

Newton’s translational motion equation and Euler’s rotational motion equation are used to describe the dynamics of a quadrotor.

2.2.1. Translational Dynamics

As depicted in Figure 1a, the speed ${\Omega }_i$ of the rotor $R_i$ generates a force $F_i$ in the direction of the rotor axis and creates a torque $M_i$ around the rotor axis. Hence,

$$F_i=k_t{\Omega }_i^2 \mathrm{and} M_i=k_d{\Omega }_i^2,$$

(4)

where $k_t$ is the thrust coefficient, and $k_d$ is the drag force coefficient. The thrust T in the z-axis direction of the body is obtained by the combined effect of forces generated by the four motors.

$$T={\displaystyle \sum _{i=1}^4{F}_i}=k_t{\displaystyle \sum _{i=1}^4{\Omega }_i^2}$$

(5)

Subsequently, the thrust vector $T^B$ of the quadrotor with respect to frame B is written as

$$T^B=\left[\begin{array}{c}{F}_x^B\\ F_y^B\\ F_z^B\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ T\end{array}\right]$$

(6)

The moment is generated by creating a differential thrust across the two motors on the same arm of the quadrotor. $M^B$ is defined as the rotational torque provided by the four rotors with respect to frame B.

$$M^B=\left[\begin{array}{c}{M}_{\varphi }^B\\ M_{\theta }^B\\ M_{\psi }^B\end{array}\right]=\left[\begin{array}{c}d{k}_t({\Omega }_{{}_4}^2-{\Omega }_{{}_2}^2)\\ d{k}_t({\Omega }_{{}_3}^2-{\Omega }_{{}_1}^2)\\ k_d({\Omega }_2^2+{\Omega }_4^2-{\Omega }_1^2-{\Omega }_3^2)\end{array}\right],$$

(7)

where d is the distance between the center of the mass and the axis of the propeller. Based on Equation (7), the roll movement is obtained by decreasing the rotor velocity of the second motor and increasing the rotor velocity of the fourth motor. The pitch movement is obtained by decreasing the first rotor velocity and increasing the third rotor velocity. The yaw movement is obtained by increasing the angular velocities of the two opposite rotors and decreasing the velocities of the other two rotors [17]. Combining Equations (6) and (7), the matrix in Equation (8) shows the relationships among the vertical thrust $F_z^B$ along the z-axis, rotational torques associated with the thrust difference of each rotor pair, i.e., $M_{\varphi }^B$, $M_{\theta }^B$ and $M_{\psi }^B$, and the angular velocities of the four propellers.

$$\left[\begin{array}{c}{F}_z^B\\ M_{\varphi }^B\\ M_{\theta }^B\\ M_{\psi }^B\end{array}\right]=\left[\begin{array}{cccc}{k}_t& k_t& k_t& k_t\\ 0& -d{k}_t& 0& d{k}_t\\ -d{k}_t& 0& d{k}_t& 0\\ -k_d& k_d& -k_d& k_d\end{array}\right]\left[\begin{array}{c}{\Omega }_{{}_1}^2\\ {\Omega }_{{}_2}^2\\ {\Omega }_3^2\\ {\Omega }_4^2\end{array}\right]$$

(8)

According to Newton’s second law of motion, the translational dynamics model of the quadrotor in frame E is expressed as follows:

$${\dot{V}}^E=\frac{1}{m}(R_B^E{T}^B-G^E-F_d),$$

(9)

where $V^E={\left[\begin{array}{ccc}{\dot{x}}^E& {\dot{y}}^E& {\dot{z}}^E\end{array}\right]}^T$ is the linear velocity along the three axes in the inertial frame; $m$ is the total mass of the quadrotor; $R_B^E$ is the rotation matrix derived from Equation (1); $T^B$ is the lift force from Equation (6); $G^E={\left[\begin{array}{ccc}0& 0& mg\end{array}\right]}^T$ is gravity, where $g$ is the acceleration due to gravity; $F_d={\left[\begin{array}{ccc}{k}_1\dot{x}& k_2\dot{y}& k_3\dot{z}\end{array}\right]}^T$ is the drag force caused by translational motions of the quadrotor, where, $k_1$, $k_2$ and $k_3$ are the air drag coefficients. Substituting Equations (1) and (6) into Equation (9) yields

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{\ddot{x}}^E\\ {\ddot{y}}^E\\ {\ddot{z}}^E\end{array}\right]& =\frac{1}{m}(R_B^E{T}^B)-g\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]-\frac{1}{m}\left[\begin{array}{c}{k}_1\dot{x}\\ k_2\dot{y}\\ k_3\dot{z}\end{array}\right]\hfill \\ & =\frac{T}{m}\left[\begin{array}{c}\cos \varphi \sin \theta \cos \psi +\sin \varphi \sin \psi \\ \cos \varphi \sin \theta \sin \psi -\sin \varphi \cos \psi \\ \cos \varphi \cos \theta \end{array}\right]-g\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]-\frac{1}{m}\left[\begin{array}{c}{k}_1\dot{x}\\ k_2\dot{y}\\ k_3\dot{z}\end{array}\right]\hfill \end{array}$$

(10)

Using Equations (6) and (10), we can extract the expressions for the pitch angle, roll angle, and thrust T, as follows:

$${\varphi }^E={\sin }^{-1}\left(\frac{({\ddot{x}}^E+\frac{k_1}{m}{\dot{x}}^E)\sin \psi -({\ddot{y}}^E+\frac{k_2}{m}{\dot{y}}^E)\cos \psi }{\sqrt{{\left({\ddot{x}}^E+\frac{k_1}{m}{\dot{x}}^E\right)}^2+{\left({\ddot{y}}^E+\frac{k_2}{m}{\dot{y}}^E\right)}^2+{\left({\ddot{z}}^E+g+\frac{k_3}{m}{\dot{z}}^E\right)}^2}}\right)$$

(11)

$${\theta }^E={\tan }^{-1}\left(\frac{({\ddot{x}}^E+\frac{k_1}{m}{\dot{x}}^E)\cos \psi +({\ddot{y}}^E+\frac{k_2}{m}{\dot{y}}^E)\sin \psi }{({\ddot{z}}^E+g+\frac{k_3}{m}{\dot{z}}^E)}\right)$$

(12)

$$\begin{array}{cc}\hfill T& =m[({\ddot{x}}^E+\frac{k_1}{m}{\dot{x}}^E)(\cos \varphi \sin \theta \cos \psi +\sin \varphi \sin \psi )\hfill \\ & +({\ddot{y}}^E+\frac{k_2}{m}{\dot{y}}^E)(\cos \varphi \sin \theta \sin \psi -\sin \varphi \cos \psi )\hfill \\ & +({\ddot{z}}^E+g+\frac{k_3}{m}{\dot{z}}^E)(\cos \varphi \cos \theta )]\hfill \end{array}$$

(13)

2.2.2. Rotational Dynamics

The quadrotor is a symmetric structure with four arms aligned with the body’s x- and y-axes. Based on Euler’s second law of motion, the rotational dynamics model of the quadrotor in frame B is expressed as follows:

$${\dot{\omega }}^B={(I^B)}^{-1}\left[(M^B-{\omega }^B\times (I^B{\omega }^B)-M_g-M_d\right]$$

(14)

where $M^B$ is the torque in Equation (7); $I^B$ is the symmetric positive definite constant matrix expressed in Equation (15); $I^B\dot{\omega }$ is the angular acceleration of the inertia; $\omega \times \left(I^B\omega \right)$ is the centripetal force, where the notation $\times$ denotes the cross product; $M_g$ is the resultant torque due to gyroscopic effects; $M_d$ is the resultant aerodynamic friction torque. The expressions for $M_g$ and $M_d$ are shown in Equation (16) and Equation (17), respectively [8,17,18,20,21].

$$I^B=\left[\begin{array}{ccc}{I}_x& 0& 0\\ 0& I_y& 0\\ 0& 0& I_z\end{array}\right]$$

(15)

$$M_g={{\displaystyle \sum _{i=1}^4\omega \times I_r\left[\begin{array}{ccc}0& 0& {(-1)}^{i+1}{\Omega }_i\end{array}\right]}}^T$$

(16)

$$M_d=\left[\begin{array}{ccc}{k}_4& 0& 0\\ 0& k_5& 0\\ 0& 0& k_6\end{array}\right]\left[\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right]=\left[\begin{array}{c}{k}_4\dot{\varphi }\\ k_5\dot{\theta }\\ k_6\dot{\psi }\end{array}\right]$$

(17)

where $I_x$, $I_y$, and $I_z$ in Equation (15) are the rotary inertia for the x-, y-, and z-axes of the body, respectively; $I_r$ in Equation (16) is the moment of inertia of each rotor; $k_4$, $k_5$, and $k_6$ in Equation (17) are the aerodynamic drag coefficients.

Considering that $I^B$ is nonsingular and substituting Equations (16) and (17) into Equation (14), we obtain

$$\begin{array}{rr}{\dot{\omega }}^B& =\left[\begin{array}{c}\dot{p}\\ \dot{q}\\ \dot{r}\end{array}\right]={(I^B)}^{-1}(M^B-{\omega }^B\times (I^B{\omega }^B)-M_g-M_d)\hfill \\ & ={(I^B)}^{-1}\cdot (M^B-\left[\begin{array}{c}p\\ q\\ r\end{array}\right]\times \left[\begin{array}{c}{I}_{x}p\\ I_{y}q\\ I_{z}r\end{array}\right]-I_r\left[\begin{array}{c}p\\ q\\ r\end{array}\right]\times \left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]{\Omega }_{\Gamma }-\left[\begin{array}{c}{k}_4\dot{\varphi }\\ k_5\dot{\theta }\\ k_6\dot{\psi }\end{array}\right]),\hfill \\ & =\left[\begin{array}{c}\frac{M_{\varphi }^B}{I_x}\\ \frac{M_{\theta }^B}{I_y}\\ \frac{M_{\psi }^B}{I_z}\end{array}\right]+\left[\begin{array}{c}\frac{(I_y-I_z)qr}{I_x}\\ \frac{(I_z-I_x)pr}{I_y}\\ \frac{(I_x-I_y)pq}{I_z}\end{array}\right]-I_r\left[\begin{array}{c}\frac{q}{I_x}\\ \frac{-p}{I_y}\\ 0\end{array}\right]{\Omega }_{\Gamma }-\left[\begin{array}{c}\frac{k_4\dot{\varphi }}{I_x}\\ \frac{k_5\dot{\theta }}{I_y}\\ \frac{k_6\dot{\psi }}{I_z}\end{array}\right]\hfill \end{array}$$

(18)

where ${\Omega }_{\Gamma }={\Omega }_4+{\Omega }_3-{\Omega }_2-{\Omega }_1$ represents the overall speed of the quadrotor rotor. As shown in Equation (3), the angular acceleration vector in frame E is obtained from the angular acceleration in frame B comprising the transformation matrix, ${(W_E^B)}^{-1}$, and its time derivative is expressed as

$${\ddot{\eta }}^E=\left[\begin{array}{c}\ddot{\varphi }\\ \ddot{\theta }\\ \ddot{\psi }\end{array}\right]=\frac{d}{dt}\left[{(W_E^B)}^{-1}{\omega }^B\right]=\left[\frac{d}{dt}{(W_E^B)}^{-1}\right]{\omega }^B+{(W_E^B)}^{-1}{\dot{\omega }}^B$$

(19)

By substituting Equations (3) and (18) into Equation (19), we obtain

$${\ddot{\eta }}^E=\left[\begin{array}{c}\ddot{\varphi }\\ \ddot{\theta }\\ \ddot{\psi }\end{array}\right]=\left[\begin{array}{c}\frac{\left(I_y-I_z\right)\dot{\theta }\dot{\psi }}{I_x}\\ \frac{\left(I_z-I_x\right)\dot{\varphi }\dot{\psi }}{I_y}\\ \frac{\left(I_x-I_y\right)\dot{\varphi }\dot{\theta }}{I_z}\end{array}\right]-I_r\left[\begin{array}{c}\frac{\dot{\theta }}{I_x}\\ \frac{-\dot{\varphi }}{I_y}\\ \underset{}{\overset{}{0}}\end{array}\right]{\Omega }_{{\rm T}}+\left[\begin{array}{c}\frac{M_{\varphi }^B}{I_x}\\ \frac{M_{\theta }^B}{I_y}\\ \frac{M_{\psi }^B}{I_z}\end{array}\right]-\left[\begin{array}{c}\frac{k_4}{I_x}\dot{\varphi }\\ \frac{k_5}{I_y}\dot{\theta }\\ \frac{k_6}{I_z}\dot{\psi }\end{array}\right]$$

(20)

2.3. Model Simplification

The motions of the quadrotor are generated by combining a series of forces and moments from different physical effects. The model developed using Equations (10) and (20) describes the differential equation of the quadrotor system. In order to implement the flight controller in the resource limited microcontroller and consider the requirements of the dynamic response for the attitude controller which is 4 to 10 times faster than that of the position controller, it is assumed that the quadrotor is operated in a hovering state, so the models of quadrotor are linearized and the high-order effects of the force and moment models of the quadrotor can be simplified. Hence, the drag force $F_d$ in Equation (9) and the aerodynamic friction torque $M_d$ in Equation (14) were disregarded, and the thrust coefficient $k_t$ and drag force coefficient $k_d$ in Equation (4) were assumed to be constants [20]. The quadrotor system can be rewritten in the state-space form $\dot{X}=f(X{,}_{}U)$. The following input vector $U$ and state vector $X$ were selected:

$$\begin{array}{l}X={\left[\begin{array}{llllllllllll}{x}_1\hfill & x_2\hfill & x_3\hfill & x_4\hfill & x_5\hfill & x_6\hfill & x_7\hfill & x_8\hfill & x_9\hfill & x_{10}\hfill & x_{11}\hfill & x_{12}\hfill \end{array}\right]}^T\\ ={\left[\begin{array}{llllllllllll}\varphi \hfill & \dot{\varphi }\hfill & \theta \hfill & \dot{\theta }\hfill & \psi \hfill & \dot{\psi }\hfill & x\hfill & \dot{x}\hfill & y\hfill & \dot{y}\hfill & z\hfill & \dot{z}\hfill \end{array}\right]}^T\end{array}$$

(21)

$$U=\left[\begin{array}{c}{U}_1\\ U_2\\ U_3\\ U_4\end{array}\right]=\left[\begin{array}{c}{F}_z^B\\ M_{\varphi }^B\\ M_{\theta }^B\\ M_{\psi }^B\end{array}\right]$$

(22)

The state representations of the system expressed in Equations (10) and (20) were simplified to the following

$$\{\begin{array}{l}{x}_1=\varphi \\ x_2={\dot{x}}_1=\dot{\varphi }\\ x_3=\theta \\ x_4={\dot{x}}_3=\dot{\theta },\begin{array}{cc}& \end{array}\\ x_5=\psi \\ x_6={\dot{x}}_5=\dot{\psi }\end{array}\{\begin{array}{l}{x}_7=x\\ x_8={\dot{x}}_7=\dot{x}\\ x_9=y\\ x_{10}={\dot{x}}_9=\dot{y}\\ x_{11}=z\\ x_{12}={\dot{x}}_{11}=\dot{z}\end{array}$$

(23)

$$f(X{,}_{}U)=\left[\begin{array}{l}\dot{\varphi }\hfill \\ \ddot{\varphi }\hfill \\ \dot{\theta }\hfill \\ \ddot{\theta }\hfill \\ \dot{\psi }\hfill \\ \ddot{\psi }\hfill \\ \dot{x}\hfill \\ \ddot{x}\hfill \\ \dot{y}\hfill \\ \ddot{y}\hfill \\ \dot{z}\hfill \\ \ddot{z}\hfill \end{array}\right]=\left[\begin{array}{l}\dot{\varphi }\hfill \\ \dot{\theta }\dot{\psi }{a}_1-\dot{\theta }{a}_2{\Omega }_{{\rm T}}+b_1{U}_2\hfill \\ \dot{\theta }\hfill \\ \dot{\varphi }\dot{\psi }{a}_3+\dot{\varphi }{a}_4{\Omega }_{{\rm T}}+b_2{U}_3\hfill \\ \dot{\psi }\hfill \\ \dot{\varphi }\dot{\theta }{a}_5+b_3{U}_4\hfill \\ \dot{x}\hfill \\ u_x\frac{1}{m}{U}_1\hfill \\ \dot{y}\hfill \\ u_y\frac{1}{m}{U}_1\hfill \\ \dot{z}\hfill \\ u_z\frac{1}{m}{U}_1-g\hfill \end{array}\right]$$

(24)

$$\{\begin{array}{l}{a}_1=\frac{\left(I_y-I_z\right)}{I_x}\\ a_2=\frac{I_r}{I_x}\\ a_3=\frac{\left(I_z-I_x\right)}{I_y}\\ a_4=\frac{I_r}{I_y}\\ a_5=\frac{\left(I_x-I_y\right)}{I_z}\end{array},\begin{array}{cc}& \end{array}\{\begin{array}{l}{b}_1=\frac{1}{I_x}\\ b_2=\frac{1}{I_y}\\ b_3=\frac{1}{I_z}\end{array},\begin{array}{cc}& \end{array}\{\begin{array}{l}{u}_x=\cos \varphi \sin \theta \cos \psi +\sin \varphi \sin \psi \\ u_y=\cos \varphi \sin \theta \sin \psi -\sin \varphi \cos \psi \\ u_z=\cos \varphi \cos \theta \end{array}$$

(25)

Subsequently, Equations (11) and (13) can be rewritten as

$${\varphi }^E={\sin }^{-1}\left(\frac{{\ddot{x}}^E\sin \psi -{\ddot{y}}^E\cos \psi }{\sqrt{{\left({\ddot{x}}^E\right)}^2+{\left({\ddot{y}}^E\right)}^2+{\left({\ddot{z}}^E+g\right)}^2}}\right)$$

(26)

$${\theta }^E={\tan }^{-1}\left(\frac{{\ddot{x}}^E\cos \psi +{\ddot{y}}^E\sin \psi }{{\ddot{z}}^E+g}\right)$$

(27)

$$T=m\left[{\ddot{x}}^E{u}_x+{\ddot{y}}^{E}u{}_y+({\ddot{z}}^E+g)u_z\right]$$

(28)

The revolutions per minute (rpm) to throttle conversion between the angular velocity of the rotors and the inputs is written as

$$\left[\begin{array}{c}{\Omega }_{{}_1}^2\\ {\Omega }_{{}_2}^2\\ {\Omega }_3^2\\ {\Omega }_4^2\end{array}\right]=\left[\begin{array}{llll}\frac{1}{4k_t}\hfill & {}^{}0\hfill & \frac{-1}{2d{k}_t}\hfill & \frac{-1}{4k_d}\hfill \\ \frac{1}{4k_t}\hfill & \frac{-1}{2d{k}_t}\hfill & {}^{}0\hfill & \frac{1}{4k_d}\hfill \\ \frac{1}{4k_t}\hfill & {}^{}0\hfill & \frac{1}{2d{k}_t}\hfill & \frac{-1}{4k_d}\hfill \\ \frac{1}{4k_t}\hfill & \frac{1}{2d{k}_t}\hfill & {}^{}0\hfill & \frac{1}{4k_d}\hfill \end{array}\right]\left[\begin{array}{c}{U}_1\\ U_2\\ U_3\\ U_4\end{array}\right]$$

(29)

3. Quadrotor Control

The SMC switches between two structures and causes the system states to approach the sliding surface. This sliding surface is defined over the state space or the error state space of the system. During the design of the SMC, a Lyapunov function was used to derive the stability of the system. The backstepping control approach sets a Lyapunov function and results in a stabilizing control law. The backstepping controller allows the entire system to be segregated into small subsystems. An error or regulatory variable is defined based on the designed Lyapunov function for each subsystem. Subsequently, a virtual control law is applied such that the subsystem is equivalent to a linear system and the derivative of the Lyapunov function is negative definite to ensure the stability of the system. However, the backstepping controller is not robust to parametric variations. To design a robust controller against matched and mismatched uncertainties, an integrator backstepping controller and an SMC are combined to achieve a continuous controlled signal, thereby eliminating the chattering effect from the control input [11].

3.1. Control Architecture

The system structure of a quadrotor is shown in Figure 2. The external uncertainties $D={\left[\begin{array}{cc}{d}_1& d_2\end{array}\right]}^T$ are considered, where $d_1$ and $d_2$ represent the force and torque disturbances, respectively. The general form of the IB-SMC that reflects the dynamics of the quadrotor is shown in Figure 3. The backstepping controller is used to control the first subsystem and the SMC is utilized to control the second subsystem.

3.2. Attitude Tracking Control

As shown in Figure 2, the rotational controller is designed to ensure the trajectory tracking capability of the current pitch angle ${\varphi }_c^E(t)$, roll angle, ${\theta }_c^E(t)$ and yaw angle ${\psi }_c^E(t)$ of the body along the desired pitch angle ${\varphi }_d^E(t)$, roll angle, ${\theta }_d^E(t)$ and yaw angle ${\psi }_d^E(t)$, respectively. Consider the first subsystem of the rotational controller with virtual control input ${\alpha }_1$ as follows:

$${\dot{\eta }}^E={\alpha }_1$$

(30)

The angle tracking error $e_1$ of the first subsystem is defined as

$$e_1={\eta }_d^E-{\eta }^E+J_1{\displaystyle \int ({\eta }_d^E-{\eta }^E)}dt,$$

(31)

where $e_1\in R^{3\times 1}$ and $J_1\in R^{3\times 3}$ are positive definite diagonal integral coefficient matrices. A backstepping controller was used to control the first subsystem. A Lyapunov function $V_1(t)=\frac{1}{2}{e}_1^T{e}_1$ was considered for the first subsystem, and the time derivative of $V_1(t)$ is expressed as

$$\begin{array}{cc}\hfill {\dot{V}}_1(t)& =e_1^T{\dot{e}}_1\hfill \\ & =e_1^T\left[{\dot{\eta }}_d^E-{\dot{\eta }}^E+J_1({\eta }_d^E-{\eta }^E)\right]\hfill \\ & =e_1^T\left[{\dot{\eta }}_d^E-{\alpha }_1+J_1({\eta }_d^E-{\eta }^E)\right]\hfill \end{array}$$

(32)

The stabilization of the first subsystem can be obtained by designing the virtual control signal ${\alpha }_1$ as follows:

$${\alpha }_1={\dot{\eta }}_d^E+A_1{e}_1+J_1({\eta }_d^E-{\eta }^E),$$

(33)

where $A_1\in R^{3\times 3}$ is a positive-definite diagonal matrix. Substituting Equation (33) into Equation (32), we obtain

$${\dot{V}}_1(t)=-e_1^T{A}_1{e}_1<0{,}^{}\forall e_1\ne 0$$

(34)

Considering the simplified version of Equation (14) discussed in Section 2.3 and the torque disturbance $d_2$, the second subsystem with virtual control input ${\alpha }_2$ can be defined as

$${\dot{\omega }}^B={(I^B)}^{-1}{\alpha }_2-f({\omega }^B)+{(I^B)}^{-1}{d}_2$$

(35)

where ${\alpha }_2\in R^{3\times 1}$ and $f({\omega }^B)$ are expressed as

$$\{\begin{array}{l}f({\omega }^B)={(I^B)}^{-1}\left[{\omega }^B\times \left(I^B{\omega }^B\right)\right]\\ {\alpha }_2=M^B-M_g\end{array}$$

(36)

The SMC was used to control the second subsystem, and the angle-tracking error $e_2$ of the second subsystem is defined as

$$e_2={\alpha }_1-{\dot{\eta }}^E$$

(37)

Next, the sliding surface $S_1$ is expressed as

$$\begin{array}{cc}\hfill S_1& ={\dot{e}}_1+A_1{e}_1\hfill \\ & =({\dot{\eta }}_d^E-{\dot{\eta }}^E)+J_1({\eta }_d^E-{\eta }^E)+A_1{e}_1\hfill \\ & ={\alpha }_1-{\dot{\eta }}^E\hfill \\ & =e_2\hfill \end{array}$$

(38)

A Lyapunov function $V_2(t)=V_1(t)+\frac{1}{2}{S}_1^T{S}_1$ was considered for the second subsystem. Combining Equations (34) and (38),the time derivative of $V_2(t)$ is written as

$$\begin{array}{cc}{\dot{V}}_2(t)\hfill & ={\dot{V}}_1(t)+S_1^T{\dot{S}}_1\hfill \\ & =(-e_1^T{A}_1{e}_1)+S_1^T\left({\ddot{e}}_1+A_1{\dot{e}}_1\right)\hfill \\ & =(-e_1^T{A}_1{e}_1)+S_1^T\left[{\ddot{\eta }}_d^E+J_1({\dot{\eta }}_d^E-{\dot{\eta }}^E)+A_1(S_1-A_1{e}_1)-{\ddot{\eta }}^E\right]\hfill \\ & =(-e_1^T{A}_1{e}_1)+S_1^T\left[{\ddot{\eta }}_d^E+J_1({\dot{\eta }}_d^E-{\dot{\eta }}^E)+A_1(S_1-A_1{e}_1)-{(I^B)}^{-1}{\alpha }_2+f({\omega }^B)-{(I^B)}^{-1}{d}_2\right]\hfill \end{array}$$

(39)

The second subsystem can be stabilized by designing the virtual control signal ${\alpha }_2$ as follows:

$${\alpha }_2=I^B\left[{\ddot{\eta }}_d^E+A_1(S_1-A_1{e}_1)+J_1({\dot{\eta }}_d^E+S_1-{\alpha }_1)+f({\omega }^B)+{\epsilon }_{1}sat(S_1)+Q_1{S}_1\right]$$

(40)

where ${\epsilon }_1\in R^{+}$ and $Q_1\in R^{3\times 3}$ are the positive real number and positive definite diagonal matrix, respectively; $sat(S)$ is the saturation function of the sliding surface, defined as [12]

$$sat(S)=\{\begin{array}{l}{{}^{}}^{}\frac{S}{\Vert S\Vert }\begin{array}{ccc}& & S>\Delta \end{array}\\ {{}^{}}^{}{}^{}0\begin{array}{ccc}& & {}^{}S=0\end{array}\\ -\frac{S}{\Vert S\Vert }\begin{array}{ccc}& & S<\Delta \end{array}\\ (\frac{S}{\Delta \Vert S\Vert })S\begin{array}{ccc}& otherwise& \end{array}\end{array},$$

(41)

where $\Delta$ is a positive constant. To improve the robustness of the rotational controller, it was assumed that ${\epsilon }_1\ge \Vert {(I^B)}^{-1}{d}_2\Vert$. By substituting Equation (40) into Equation (39), we obtain

$$\begin{array}{cc}\hfill {\dot{V}}_2(t)& =(-e_1^T{A}_1{e}_1)-S_1^T{Q}_1{S}_1-S_1^T{\epsilon }_{1}sat(S_1)-S_1^T{(I^B)}^{-1}{d}_2\hfill \\ & \le (-e_1^T{A}_1{e}_1)-S_1^T{Q}_1{S}_1-{\epsilon }_1\Vert S_1\Vert +\Vert {(I^B)}^{-1}{d}_2\Vert \Vert S_1\Vert \hfill \\ & =(-e_1^T{A}_1{e}_1)-S_1^T{Q}_1{S}_1-({\epsilon }_1-\Vert {(I^B)}^{-1}{d}_2\Vert )\Vert S_1\Vert \hfill \\ & <0{,}^{}\forall e_1\ne 0{,}^{}{S}_1\ne 0\hfill \end{array}$$

(42)

Based on Equations (34) and (42), the stability of the rotational controller of the quadrotor is globally asymptotically stable; hence, the attitude trajectory tracking capability is confirmed [12].

3.3. Position Tracking Control

As shown in Figure 2, the objective of the translational controller is to ensure the excellent tracking performance of the position and yaw angle along the desired reference path $P_d^E={[\begin{array}{ccc}{x}_d^E(t)& y_d^E(t)& z_d^E(t)\end{array}]}^T$ and ${\psi }_d^E(t)$. The first subsystem of the translational controller with a virtual control input ${\alpha }_3$ is defined as

$${\dot{P}}^E={\alpha }_3$$

(43)

The second subsystem of the translational controller with a virtual control input ${\alpha }_4$ and force disturbance $d_1$ is defined as

$${\dot{V}}^E={\alpha }_4+\frac{1}{m}{d}_1$$

(44)

where ${\alpha }_4=\frac{1}{m}(R_B^E{T}^B-G^E)$. The position tracking error $e_3$ of the first subsystem is defined as

$$e_3=P_d^E-P^E+J_2{\displaystyle \int (P_d^E-P^E)}dt$$

(45)

where $J_2\in R^{3\times 3}$ is a positive-definite diagonal matrix. A backstepping controller was used to control the first subsystem. A Lyapunov function $V_3(t)=\frac{1}{2}{e}_3^T{e}_3$ was considered for the first subsystem, and the time derivative of $V_3(t)$ is written as

$$\begin{array}{cc}\hfill {\dot{V}}_3(t)& =e_3^T{\dot{e}}_3\hfill \\ & =e_3^T\left[{\dot{P}}_d^E-{\dot{P}}^E+J_2(P_d^E-P^E)\right]\hfill \\ & =e_3^T\left[{\dot{P}}_d^E-{\alpha }_3+J_2(P_d^E-P^E)\right]\hfill \end{array}$$

(46)

The first subsystem can be stabilized by designing the virtual control signal ${\alpha }_3$ as

$${\alpha }_3={\dot{P}}_d^E+A_2{e}_3+J_2(P_d^E-P^E)$$

(47)

where $A_2\in R^{3\times 3}$ is a positive-definite diagonal matrix. Substituting Equation (47) into Equation (46) yields

$${\dot{V}}_3(t)=-e_3^T{A}_2{e}_3<0{,}^{}\forall e_3\ne 0$$

(48)

The SMC was used to control the second subsystem, as shown in Equation (44), and the position tracking error $e_4$ of the second subsystem is defined as

$$e_4={\alpha }_3-{\dot{P}}^E$$

(49)

The sliding surface $S_2$ is defined as

$$\begin{array}{cc}\hfill S_2& ={\dot{e}}_3+A_2{e}_3\hfill \\ & ={\dot{P}}_d^E-{\dot{P}}^E+J_2(P_d^E-P^E)+A_2{e}_3\hfill \\ & ={\alpha }_3-{\dot{P}}^E\hfill \\ & =e_4\hfill \end{array}$$

(50)

A Lyapunov function $V_4(t)=V_3(t)+\frac{1}{2}{S}_2^T{S}_2$ was considered for the second subsystem. Combining Equations (48) and (50), the time derivative of $V_4(t)$ is expressed as

$$\begin{array}{cc}\hfill {\dot{V}}_4(t)& ={\dot{V}}_3(t)+S_2^T{\dot{S}}_2\hfill \\ & =(-e_3^T{A}_2{e}_3)+S_2^T({\dot{\alpha }}_3-{\ddot{P}}^E)\hfill \\ & =(-e_3^T{A}_2{e}_3)+S_2^T\left[{\ddot{P}}_d^E+A_2(S_2-A_2{e}_3)+J_2({\dot{P}}_d^E-{\dot{P}}^E)-{\alpha }_2-\frac{1}{m}{d}_1\right]\hfill \end{array}$$

(51)

The second subsystem can be stabilized by designing the virtual control signal ${\alpha }_4$ as

$${\alpha }_4={\ddot{P}}_d^E+A_2(S_2-A_2{e}_3)+J_2({\dot{P}}_d^E+S_2-{\alpha }_3)+Q_2{S}_2+{\epsilon }_{2}sat(S_2)$$

(52)

where ${\epsilon }_2\in R^{+}$ and $Q_2\in R^{3\times 3}$ are the positive real number and positive definite diagonal matrix, respectively; $sat(S)$ is the saturation function of the sliding surface shown in Equation. (41). To improve the robustness of the controller, it was assumed that ${\epsilon }_2\ge \frac{1}{m}\Vert d_1\Vert$. By substituting Equation (52) into Equation (51), we obtain

$$\begin{array}{cc}{\dot{V}}_4(t)\hfill & =(-e_3^T{A}_2{e}_3)+S_2^T\left[-Q_2{S}_2-{\epsilon }_{2}sat(S_2)-\frac{1}{m}{d}_1\right]\hfill \\ & \le (-e_3^T{A}_2{e}_3)-S_2^T{Q}_{2}S-{\epsilon }_2\Vert S_2\Vert +\frac{1}{m}\Vert d_1\Vert \Vert S_2\Vert \hfill \\ & =(-e_3^T{A}_2{e}_3)-S_2^T{Q}_{2}S-({\epsilon }_2-\frac{1}{m}\Vert d_1\Vert )\Vert S_2\Vert \hfill \\ & <0{,}^{}\forall e_3\ne 0{,}^{}{S}_2\ne 0\hfill \end{array}$$

(53)

The combination of Equations (48) and (53) indicates that the stability of the translational controller of the quadrotor is globally asymptotically stable, and that positional trajectory tracking capability is confirmed [12].

4. Simulations and Experimental Results

In this section, numerical simulations and experiments for trajectory tracking were implemented to validate the merit and effectiveness of the IB-SMC of the quadrotor. By considering external disturbances, the trajectory tracking results, including those for take-off, hovering, and landing, are shown in Section 4.1. The simulations were conducted using the MATLAB & Simulink software tool, and the computations were performed on a personal computer with an Intel Core i5-8265 CPU, 16 GB of RAM, and an Nvidia Geforce MX110 graphical card. The trajectory tracking experiments for the IB-SMC of the quadrotor are explained in Section 4.2.

4.1. Simulations

The designs of the IB-SMC of the quadrotor for attitude control and position control described in Section 3 were simulated. For simplicity, the physical phenomenon presented in the real world, i.e., fluid dynamics, the friction of the four rotors, the error range from the GPS modules which are mounted on the quadrotor are ignored, therefore, the conceptual behavior for the trajectory tracking of the tracking controller in the real world are simulated. The parameters of the quadrotor model are listed in

Table 1. Model parameters of quadrotor.

Scheme	Value (Units)	Symbol	Value (Unit)
m	0.5 (kg)	d	0.2 (m)
$I_x$	$2.35\times {10}^{-3}{}^{}(\mathrm{kg}\cdot {\mathrm{m}}^2)$	g	${9.81}^{}(\mathrm{m}/{\mathrm{s}}^2)$
$I_y$	$2.35\times {10}^{-3}{}^{}(\mathrm{kg}\cdot {\mathrm{m}}^2)$	$k_t$	$6.13\times {10}^{-5}{}^{}(\mathrm{N}\cdot {\mathrm{s}}^2)$
$I_z$	$5.26\times {10}^{-2}{}^{}(\mathrm{kg}\cdot {\mathrm{m}}^2)$	$k_d$	$2.5\times {10}^{-6}{}^{}(\mathrm{N}\cdot \mathrm{m}\cdot {\mathrm{s}}^2)$

. The parameters of the quadrotor controllers are listed in

Table 2. Controller parameters.

Symbol	Value (Unit)	Symbol	Value (Unit)
$A_1$	$1.5I_{3\times 3}$	$A_2$	$9I_{3\times 3}$
${\epsilon }_1$	10	${\epsilon }_2$	2
$Q_1$	$5I_{3\times 3}$	$Q_2$	$50I_{3\times 3}$
$J_1$	$3I_{3\times 3}$	$J_2$	$0.9I_{3\times 3}$
$\Delta$	5

. The initial states of the quadrotor were ${\eta }^E(0)=0$, $V^E(0)=0$, and ${\omega }^B(0)=0$.

To demonstrate the physical disturbance, the force $d_1$ and torque disturbance $d_2$ were assumed as the following:

$$\{\begin{array}{l}{d}_1=d_{1,static}+d_{1,random}\begin{array}{ccc}& & 0 \mathrm{s}\le t<5 \mathrm{s}\end{array}\\ d_2=d_{2,static}+d_{2,random}\begin{array}{ccc}& & 5 \mathrm{s}\le t<10 \mathrm{s}\end{array}\end{array}$$

(54)

$$\begin{array}{l}{d}_{1,static}=0.25{\left[\begin{array}{ccc}1& 1& 1\end{array}\right]}^T,{\begin{array}{l}\hfill \end{array}}^{}-0.05{\left[\begin{array}{ccc}1& 1& 1\end{array}\right]}^T\le d_{1,random}\le 0.05{\left[\begin{array}{ccc}1& 1& 1\end{array}\right]}^T\\ d_{2,static}=0.1{\left[\begin{array}{ccc}1& 1& 1\end{array}\right]}^T,{\begin{array}{l}\hfill \end{array}}^{}-0.05{\left[\begin{array}{ccc}1& 1& 1\end{array}\right]}^T\le d_{2,random}\le 0.05{\left[\begin{array}{ccc}1& 1& 1\end{array}\right]}^T\end{array}$$

(55)

Here, $d_{i,static}$ ($i=1,2$) represents a slow change, and the unknown constant disturbances along with time $d_{i,static}$ can be regarded as static wind or unmodeled dynamic errors; $d_{i,random}$ ($i=1,2$) represents the fast varying and unknown stochastic disturbance with the assumption of a uniform distribution, and $d_{i,random}$ can be regarded as stochastic wind or uncertain measurement noise; $t$ is the simulation time [12].

To emulate the tracking performance of the quadrotor in the experiment, five trajectory paths were planned, as follows:

$$(x,y,z)=\{\begin{array}{l}(0, 0,-5+2.5t)\\ (1.5(t-2), 2(t-2), 0)\\ (3-1.5(t-4), 4+2(t-4), 0)\\ (-1.5(t-6), 8-2(t-6), 0)\\ (-3+1.5(t-8), 4-2(t-8), 0)\end{array}{\begin{array}{l}\hfill \end{array}}^{}\begin{array}{l}0 \mathrm{s}\le t\le 2 \mathrm{s}\hfill \\ 2 \mathrm{s}<t\le 4 \mathrm{s}\hfill \\ 4 \mathrm{s}<t\le 6 \mathrm{s}\hfill \\ 6 \mathrm{s}<t\le 8 \mathrm{s}\hfill \\ 8 \mathrm{s}<t\le 10 \mathrm{s}\hfill \end{array}$$

(56)

The initial posture of the quadrotor was set to $P_0^E(t)=P_0^E(0)=(0,0,-5)$ at $t=0 \mathrm{s}$. Five target positions were set: the first position was located at $P_1^E(2)=(0,0,0)$ at $t=2 \mathrm{s}$, the second position was $P_2^E(4)=(3,4,0)$, the third position was $P_3^E(6)=(0,8,0)$, the fourth position was $P_4^E(8)=(-3,4,0)$; Subsequently, the quadrotor returns to the final target position $P_5^E(10)=(0,0,0)$ at $t=10 \mathrm{s}$. Figure 4 shows the tracking for the desired, IB-SMC trajectories and IB-SMC tracking velocity vector.

Figure 5a,c,e show the position response of the trajectory tracking presented in Figure 4. The position error between the desired and real trajectories with disturbance is shown in Figure 5b,d,f. At $t=5 \mathrm{s}$, the position errors in the x-, y-, and z-axes were −0.0009372, −0.001187, and −0.001044, respectively.

Figure 6a,c,e show the Euler’s angle response for the trajectory tracking. The Euler’s angle error between the desired and real trajectories under disturbance is shown in Figure 6b,d,f. As shown in these figures, when disturbance was considered, the developed IB-SMC enabled the attitude-tracking capability of the quadrotor.

As shown in Figure 7a, the value of the thrust T was $3.142\times {10}^4$. However, owing to the physical limitation of the rotor, the speed of the rotor was limited; therefore, the real thrust was less than $3.142\times {10}^4$. Figure 7b–d show the rotational torques $M_{\varphi },$ $M_{\theta }$ and $M_{\psi }$, and the variation range of these rotational torques was from $-0.05 \mathrm{N}\cdot \mathrm{m}$ to $-0.15 \mathrm{N}\cdot \mathrm{m}$. It was observed that the developed IB-SMC of the quadrotor was robust against disturbance.

4.2. Experiments

Figure 8 shows the prototype of the proposed quadrotor. The hardware components of the quadrotor are listed in

Table 3. Hardware components quadrotor model parameters.

Component	Part Name	Specifications
Airframe	Tarot FY450	Diameter: 450 mm Weight: 282 g X-configuration skeleton
$4$ BLDC motors	Tarot2214–920 KV	920 rpm/V Maximum current: 10.3 A $\mathrm{Resistance}: 180 \Omega$
Electric battery	ATC Li–Po 3 Cells	Voltage: 11.1 V Capacity: 5200 mAh
GPS module	NEO–M8N	Frequency: 1575.42 MHz, Transmission rate: 9600 bps
$4$ propellers	$10$4.5	2 blades (CC and CCW pair)
RF module	RFD900	Frequency: 902–928 MHz
Electronic speed controller	ORun 30 A	PWM to PPM
Flight control module	Pixhawk	Version: 2.4.8

.

In order to validate the path tracking performance of the proposed quadrotor, the original proportional-PID (P-PID) translational controller implemented in Pixhawk was replaced by the IB-SMC translational controller. Figure 9 shows the trajectory tracking plan in the experiment, where the tracking targets were set using the QGround control (QGC) software tool in the ground station. The experiments were performed at the stadium of the university. The tracking targets were denoted as $P_1$, $P_2$, $P_3$ and $P_4$ on the Google map. The original position was $P_1$, which was located at $({24}^{\circ }{58}^{\prime }{05.9}^{″} \mathrm{N}, {121}^{\circ }{11}^{\prime }{25.2}^{″} \mathrm{E})$; $P_2$ was at $({24}^{\circ }{58}^{\prime }{06.9}^{″} \mathrm{N}, {121}^{\circ }{11}^{\prime }{24.4}^{″} \mathrm{E})$, $P_3$ was at $({24}^{\circ }{58}^{\prime }{06.1}^{″} \mathrm{N}, {121}^{\circ }{11}^{\prime }{23.8}^{″} \mathrm{E})$, and $P_4$ was at $({24}^{\circ }{58}^{\prime }{05.5}^{″} \mathrm{N}, {121}^{\circ }{11}^{\prime }{24.5}^{″} \mathrm{E})$. The tracking path was from $P_1$ to $P_4$, followed a return to $P_1$ by passing through $P_2$ and $P_3$. The flight height of the quadrotor was set to 5 m.

The trajectory tracking processes of the quadrotor are shown in Figure 10a–e. The link to the video is provided in the Supplementary Material. The flight data of the quadrotor were obtained using the Pixhawk module and then transmitted via the RF module. Subsequently, these flight data were received by the ground station. Figure 11 shows the position response of the trajectory tracking in the experiment. As shown, the developed IB-SMC of the quadrotor can provide effectiveness trajectory tracking in the experiment.

5. Conclusions and Future Works

In this study, the design of the IB-SMC of a quadrotor for trajectory position tracking was investigated. A GPS and six-axis IMU sensors, which constituted the Pixhawk module, were utilized to measure the position and attitude of the quadrotor. Theoretical analysis and numerical simulations for designing the IB-SMC were provided. Finally, experiments were performed to validate the feasibility and effectiveness of the proposed IB-SMC for a quadrotor.

For designing the next generation of the quadrotor, the quadrotor with the capabilities of obstacle avoidance, slug-payload are two main research topics. For obstacle detection and collision avoidance approaches, the readouts of the ultrasonic sensor, millimeter-Wave (mmWave) radar sensor and Light-detection and ranging (LIDAR) sensor can be integrated by the multi-sensor fusion algorithm and then pattern recognition technologies are used to identify objects by utilizing the machine learning algorithms [22,23,24]. For aerial transportation applications, delivering cargo via quadrotor is a new trend. However, a quadrotor with a slug-payload system is an underactuated, strong coupling, and nonlinear system. How to achieve the quadrotor with precise positioning, maximum payload, and minimum oscillation of suspended load simultaneously needs be investigated in the future [25]. Moreover, for extending the operating time of quadrotor, power management strategies need to be carefully planned further.