Contribution to the Control and Command of a Quadrirotor with Six Degrees of Freedom in an Urban Environment

Ndzana Jean Calvin^1,2, Yeremou Aurelien^1,2, Mane Mane Jeannot^1,2, Assiene Mouodo Luc Vivien³, Nkotto Akono Cyrille^1,2

¹Higher Normal School of Technical Education, University of Douala, Douala, Cameroon.
²Doctoral Training Unit in Engineering Sciences, Computer and Automation Engineering Laboratory, Doctoral School of and Applied Sciences, University of Douala, Douala, Cameroon.
³Materials and Methods Modeling Laboratory of the National Higher Polytechnic School, University of Douala, Douala, Cameroon.

DOI: 10.4236/wjet.2020.84059                

Abstract

The objective of this article is to make a contribution relating to the modeling, control, simulation and stabilization of a complex system, with six degrees of freedom of a particular drone which presents many advantages and challenges. On the technological, military, political and other levels with an enormous and beneficial social contribution, it is a quadrotor which is a nonlinear, strongly coupled and unstable system. Such a structure is difficult to master, because the control is multivariable in the sense that six degrees of freedom are to be controlled simultaneously and operating in an environment subject to disturbances. Two commands, in particular Backstepping and PID, will be applied to obtain the stabilization of the quadcopter at the desired values, in attitude and in altitude. This article presents the comparative results of the performance of the quadcopter under the two controls. The effect of the parameters of each command on the response time of the system is elucidated under the Matlab/Simulink environment. For a simulation time of up to 10 seconds minimum with a less good response time of almost 2 seconds for the PID control, these results prove the robustness of the Backstepping command.

Keywords

Six-Degree Quadrirotor, Control Backstepping and PID

1. Introduction

Unmanned Aerial Vehicle (UAV) is commonly known as drones, constituting the flagship of the aerospace industry. Its applications are as visible, useful and remarkable in military projects, as regards surveillance, intelligence or combat missions, remote sensing, destruction of targets and any other task that would put a crew in danger. Drones have also proven to have a very high potential for civilian applications, such as surveillance and observation of risk areas, infrastructure management, customs, the environment, field mapping, scientific experiments, intervention in hostile sites etc. [1].

A quadrotor is a miniature rotary-wing drone with four rotors. It has been the subject of several research projects focusing on its ease of implementation, vertical take-off and landing, and its efficiency in hovering. One of its major characteristics is stabilization in attitude and altitude. This stabilization allows the quadcopter to be steered to the desired position (x, y, z) or the desired angle (ϕ, θ, ψ) [2]. The literature is full of several mathematical models describing the translational and rotational dynamics of a quadrotor, with six degrees of freedom, as well as the control algorithms. Among defined, the linear controls PID [1] [2] [3], LQR, H∞ [4] applied to a linearized model (around an operating point using the Jacobian method). The nonlinear controls take into account the important nonlinearities of the vehicle dynamics. We can cite the hierarchical control [5], the sliding mode control [6], the Backstepping algorithm [7] - [12]. Other commands exploit visual data extracted from acquired images; we will speak of visual serving or command by vision [13]. Some have also been interested in author’s flight commands based on neural networks and fuzzy logic [13]. In this article, we will focus on two approaches to controlling a quadcopter. This is the Backstepping controller first, and the PID controller second. Each of the commands will be simulated in the Matlab/Simulink environment. The results will be compared by analyzing their performances in terms of system response time, stabilization, convergence towards the desired setpoints and their behavior in the face of disturbances. The article is divided into three parts: the first describes the dynamic model of the quadrotor according to the Newton-Euler formalism. The proposed control strategies will find their applications in the second part. In the third part, the discussion of the results and the analysis of the performance of each control algorithm will be given.

2. Dynamic Model

­ The flight dynamics of a quadrotor are complex and strongly coupled. To develop a dynamic model close to the real physical model, taking into account the effects of the environment, we make these assumptions:

­ The structure of the quadrotor is assumed to be rigid and symmetrical, which implies that the inertia matrix will be assumed to be diagonal.

­ The propellers are supposed to be rigid in order to be able to neglect the effect of their deformation during rotation.

­ The center of mass and the origin of the coordinate system linked to the structure of the quadrotor coincide.

­ The lift and drag forces are proportional to the squares of the rotational speed of the rotors, which is a very close approximation of aerodynamic behavior.

1) Definition of benchmarks (Figure 1)

­ The movement of a quadrotor is described by two reference marks:

­ The Earth Inertial frame noted E:

$E=\left(O,x_e;y_e;z_e\right)$

­ The mobile mark (in English Body frame) noted B:

$B=\left(G,x_b;y_b;z_b\right)$

2) Euler angles

Euler angles are angles introduced by Leonhard Euler to describe the orientation of a solid. We have three: the roll angle ( $\frac{-\pi }{2}<\varphi <\frac{\pi }{2}$ ), pitch angle ( $\frac{-\pi }{2}<\theta <\frac{\pi }{2}$ ) and the yaw angle ( $-\pi <\psi <\pi$ ). A rotation is obtained by varying one of Euler’s three angles and a sequence of three rotations is sufficient to describe any transformation based on Euler’s theorem (Figure 2).

3) Kinematic parameters

They include four vectors of three elements describing the state of the system. The first two define the position and speed of the center of gravity of the quadrotor and the other two relate to its orientation and angular speed.

$r\left(t\right)={\left[xyz\right]}^{\text{T}}$ is the position of the center of gravity of the quadrirotor with respect to the reference {E} and defined in the reference {E}.

Figure 1. Repères.

Figure 2. Angles d’Euler.

$\eta \left(t\right)={\left[\varphi \theta \psi \right]}^{\text{T}}$ is the vector of Euler angles representing the angular position of the frame {B} linked to the quadrirotor with respect to the frame {E} and defined in the frame {E}.

$\upsilon \left(t\right)={\left[\stackrel{\dot{}}{x}\stackrel{\dot{}}{y}\stackrel{\dot{}}{z}\right]}_E^{\text{T}}={\left[uvw\right]}_B^{\text{T}}$ is the speed of the center of gravity compared to the reference vector E and described in the reference {B}. Also called linear speed. $\Omega \left(t\right)={\left[pqr\right]}_E^{\text{T}}={\left[\stackrel{\dot{}}{\varphi }\stackrel{\dot{}}{\theta }\stackrel{\dot{}}{\psi }\right]}_B^{\text{T}}$ is the instantaneous speed vector of rotation between the reference marks {E} and {B} and described in the reference {B}. Also called angular speed of rotation.

4) Kinetic parameters

They establish the link between the forces and the variation of the kinematic parameters. They are constant in our study. We can cite: the mass m of the quadrotor, the center of gravity G, the matrix of inertia J of the point G of the moving frame

$J=\left[\begin{array}{ccc}{J}_{xx}& 0& 0\\ 0& J_{yy}& 0\\ 0& 0& J_{zz}\end{array}\right]$ (1)

5) Rotation matrix

The rotation matrix is an orthogonal matrix and in fact it belongs to the subspace of orthogonal matrices, called special orthogonal group, denoted by SO (3), defined by:

$\left\{\Re \in {ℝ}^{3\times 3}/{\Re }^{\text{T}}\Re =I_{3\times 3}\text{et}\det \left(\Re \right)=1\right\}$

In the case where the mobile frame {B} linked to the quadrirotor rotates in inertial relation {E} with an angular speed Ω, the rotation matrix varies; this results in the kinematic height matrix: $Ṙ=Rsk(\; \Omega \; )$

sk(Ω) is the antisymmetric tensor associated with the angular speed Ω Î R³. The rotation matrix defines three angles of Euler rotation around the axes respectively described as:

${\Re }_x\left(\varphi \right)$ ; ${\Re }_y\left(\theta \right)$ ; ${\Re }_z\left(\psi \right)$. The product of these three elementary matrices of rotation forms the matrix of rotation R.

${\Re }_x\left(\varphi \right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \varphi & -\sin \varphi \\ 0& \sin \varphi & \cos \varphi \end{array}\right]$

${\Re }_y\left(\theta \right)=\left[\begin{array}{ccc}\cos \theta & 0& \sin \theta \\ 0& 1& 0\\ -\sin \theta & 0& \cos \theta \end{array}\right]$

${\Re }_z\left(\psi \right)=\left[\begin{array}{ccc}\cos \psi & -\sin \psi & 0\\ \sin \psi & \cos \psi & 0\\ 0& 0& 1\end{array}\right]$

$\Re ={\Re }_x\left(\varphi \right)\times {\Re }_y\left(\theta \right)\times {\Re }_z\left(\psi \right)$ (2)

$\Re =\left[\begin{array}{ccc}{c}^1\theta c\psi & c\psi s^2\varphi s\theta -c\varphi s\psi & c\varphi s\theta c\psi +s\varphi s\psi \\ c\theta s\psi & s\varphi s\theta s\psi +c\varphi c\psi & c\varphi s\theta s\psi -s\varphi c\psi \\ -s\theta & s\varphi c\theta & c\varphi c\theta \end{array}\right]$

6) Aerodynamic forces

A moving quadrirotor is under the action of four (4) main forces: the weight P, the total thrust U₁, the drag in the propellers T_r (propeller) and the drag along the axes T_r(ax).

Note: $P=m{g}_{z_e}=\left[\begin{array}{c}0\\ 0\\ mg\end{array}\right]$

$U_1={\displaystyle \underset{i=1}{\overset{4}{\sum }}b{\omega }_i^2}$

$T_r\left(\text{helice}\right)=d{\omega }_i^2$

$T_r\left(\text{axe}\right)=K_{ft}V=\left[\begin{array}{c}{K}_{ftx}\\ K_{fty}\\ K_{ftz}\end{array}\right]\left[\begin{array}{c}\stackrel{\dot{}}{x}\\ \stackrel{\dot{}}{y}\\ \stackrel{\dot{}}{z}\end{array}\right]=\left[\begin{array}{c}{K}_{ftx}\stackrel{\dot{}}{x}\\ K_{fty}\stackrel{\dot{}}{y}\\ K_{ftz}\stackrel{\dot{}}{z}\end{array}\right]$

Avec:

m is the mass of the quadrirotor in kg; g is the gravitational acceleration in m/s²; $z_e={\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^{\text{T}}$ is a unit vector in the terrestrial inertial frame {E}.

b is the coefficient of lift; ω_i the speed of rotation of a motor number i (i ranging from 1 to 4).

d is the coefficient of proportionality of drag.

$K_{ftx}$ Translation drag coefficient along the x axis.

$K_{fty}$ Translation drag coefficient along the y axis.

$K_{ftz}$ Translation drag coefficient along the z axis.

7) Moments due to thrust forces of rotors along the axes

­ We have two moments due to the pushing force:

­ Moment due to thrust forces on the axis x:

$U_2=\mathcal{l}b\left({\omega }_4^2-{\omega }_2^2\right)$ (3)

$U_2$ is called roll control.

Moment due to thrust forces on the axis:

$U_3=\mathcal{l}b\left({\omega }_3^2-{\omega }_1^2\right)$ (4)

$U_3$ is called pitch control.

8) Moments due to drag forces

­ They left as follows:

­ The moment due to the drag forces in the propellers along the axis z:

$U_4=d\left({\omega }_1^2-{\omega }_2^2+{\omega }_3^2-{\omega }_4^2\right)$ (5)

$U_4$ is called lace command.

By grouping together with the expression of the total thrust $U_1$ and Equations (3), (4) and (5),

$\{\begin{array}{l}{U}_1=b\left({\omega }_1^2+{\omega }_2^2+{\omega }_3^2+{\omega }_4^2\right)\\ U_2=\mathcal{l}b\left({\omega }_4^2-{\omega }_2^2\right)\\ U_3=\mathcal{l}b\left({\omega }_3^2-{\omega }_1^2\right)\\ U_4=d\left({\omega }_1^2-{\omega }_2^2+{\omega }_3^2-{\omega }_4^2\right)\end{array}$ (6)

­ The moment resulting from aerodynamic friction:

$M_a=K_{f_a}{\Omega }^2=\left[\begin{array}{c}{K}_{f_{ax}}{\stackrel{\dot{}}{\varphi }}^2\\ K_{f_{ay}}{\stackrel{\dot{}}{\theta }}^2\\ K_{f_{az}}{\stackrel{\dot{}}{\psi }}^2\end{array}\right]$ (7)

$K_{f_{ax}}$ coefficient of aerodynamic friction along the axis x.

$K_{f_{ay}}$ coefficient of aerodynamic friction along the axis y.

$K_{f_{az}}$ coefficient of aerodynamic friction along the axis z.

9) Gyroscopic effects

The gyroscopic effect is defined as the difficulty of modifying the position or the orientation of the plane of rotation of a rotating mass. Each rotor can thus rotate around its own vertical axis with a speed of rotation ω_i. Even the vertical axis moves as the vehicle rotates. This action produces two additional moments among which, on one: the gyroscopic moment of the helices M_gh and the gyroscopic moment due to the movements of the quadrirotor M_gm. Note:

$M_{gh}={\displaystyle \underset{i=1}{\overset{4}{\sum }}{\Omega }_r}\wedge J_r{\left[\begin{array}{ccc}0& 0& {\left(-1\right)}^{i+1}{\omega }_i\end{array}\right]}^{\text{T}}$

$M_{gh}=\left[\begin{array}{c}\stackrel{\dot{}}{\varphi }\\ \stackrel{\dot{}}{\theta }\\ \stackrel{\dot{}}{\psi }\end{array}\right]\wedge J_r\left[\begin{array}{c}0\\ 0\\ {\Omega }_r\end{array}\right]=\left[\begin{array}{c}{J}_r{\Omega }_r\stackrel{\dot{}}{\theta }\\ -J_r{\Omega }_r\stackrel{\dot{}}{\varphi }\\ 0\end{array}\right]$ (8)

or $\wedge$ : vector product operator; $J_r$ is the inertia of the rotors and ${\Omega }_r={\omega }_2+{\omega }_4-{\omega }_1-{\omega }_3$

$M_{gm}=\Omega \wedge J\Omega$

$M_{gm}=\left[\begin{array}{c}\left(J_{zz}-J_{yy}\right)\stackrel{\dot{}}{\theta }\stackrel{\dot{}}{\psi }\\ \left(J_{xx}-J_{zz}\right)\stackrel{\dot{}}{\varphi }\stackrel{\dot{}}{\psi }\\ \left(J_{yy}-J_{xx}\right)\stackrel{\dot{}}{\varphi }\stackrel{\dot{}}{\theta }\end{array}\right]$ (9)

10) Translational dynamics

The equations governing the translational motion of a quadrotor are:

$\{\begin{array}{l}\stackrel{\dot{}}{r}=\upsilon \\ m\ddot{r}=\Re \times \left[\begin{array}{ccc}0& 0& 1\end{array}\right]\times U_1-T_r\left(\text{axe}\right)-P\end{array}$ (10)

11) Rotational dynamics

The equations governing the rotational motion of a quadrotor are:

$\{\begin{array}{l}\stackrel{\dot{}}{\Re }=\Re sk\left(\Omega \right)\\ J\stackrel{\dot{}}{\Omega }=-M_{gm}-M_a-M_{gh}+M_f\end{array}$ (11)

12) Complete dynamic model (12)

From Equations (10) and (11) we obtain the complete dynamic model:

$\{\begin{array}{l}\ddot{x}=-\frac{K_{ftx}}{m}\stackrel{\dot{}}{x}+\frac{1}{m}{U}_x{U}_1\\ \ddot{y}=-\frac{K_{fty}}{m}\stackrel{\dot{}}{y}+\frac{1}{m}{U}_y{U}_1\\ \ddot{z}=-\frac{K_{ftz}}{m}\stackrel{\dot{}}{z}+\frac{\cos \varphi \cos \theta }{m}{U}_1-g\\ \ddot{\varphi }=\frac{J_{yy}-J_{zz}}{J_{xx}}\stackrel{\dot{}}{\theta }\stackrel{\dot{}}{\psi }-\frac{J_r}{J_{xx}}{\Omega }_r\stackrel{\dot{}}{\theta }-\frac{K_{f_{ax}}}{J_{xx}}{\stackrel{\dot{}}{\varphi }}^2+\frac{1}{J_{xx}}{U}_2\\ \ddot{\theta }=\frac{J_{zz}-J_{xx}}{J_{yy}}\stackrel{\dot{}}{\varphi }\stackrel{\dot{}}{\psi }+\frac{J_r}{J_{yy}}{\Omega }_r\stackrel{\dot{}}{\varphi }-\frac{K_{f_{ax}}}{J_{yy}}{\stackrel{\dot{}}{\theta }}^2+\frac{1}{J_{yy}}{U}_3\\ \ddot{\psi }=\frac{J_{xx}-J_{yy}}{J_{zz}}\stackrel{\dot{}}{\varphi }\stackrel{\dot{}}{\theta }-\frac{K_{f_{ax}}}{J_{zz}}{\stackrel{\dot{}}{\psi }}^2+\frac{1}{J_{zz}}{U}_4\end{array}$

and:

$\{\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 \end{array}$ (13)

13) System state representation

For a physical system there are a multitude of state representations, in our case we choose the state vector as follows:

$X=\left[\begin{array}{ccc}{r}^{\text{T}}& {\upsilon }^{\text{T}}& \begin{array}{cc}{\eta }^{\text{T}}& {\Omega }^{\text{T}}\end{array}\end{array}\right]$ (14)

Each of the 4 vectors of relation (14) has 3 components. In total, our state vector contains (4 × 3) = 12 elements and the control vector U has 4 elements such that:

$\{\begin{array}{l}X={\left[\begin{array}{ccc}\varphi & \stackrel{\dot{}}{\varphi }& \begin{array}{cc}\theta & \begin{array}{ccc}\stackrel{\dot{}}{\theta }& \psi & \begin{array}{ccc}\stackrel{\dot{}}{\psi }& x& \begin{array}{ccc}\stackrel{\dot{}}{x}& y& \begin{array}{cc}\stackrel{\dot{}}{y}& \begin{array}{cc}z& \stackrel{\dot{}}{z}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]}^{\text{T}}\\ U={\left[\begin{array}{ccc}{U}_1& U_2& \begin{array}{cc}{U}_3& U_4\end{array}\end{array}\right]}^{\text{T}}\end{array}$

We get the state representation in the form:

$\stackrel{\dot{}}{X}=f\left(X,U\right)$

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

We ask:

$\{\begin{array}{l}{a}_1=\frac{J_{yy}-J_{zz}}{J_{xx}},a_2=-\frac{K_{f_{ax}}}{J_{xx}},a_3=-\frac{J_r}{J_{xx}}\\ b_1=\frac{1}{J_{xx}},a_4=\frac{J_{zz}-J_{xx}}{J_{yy}},a_5=-\frac{K_{f_{ax}}}{J_{yy}}\\ a_6=\frac{J_r}{J_{yy}},b_2=\frac{1}{J_{yy}},a_7=\frac{J_{xx}-J_{yy}}{J_{zz}}\\ a_8=-\frac{K_{f_{ax}}}{J_{zz}},b_3=\frac{1}{J_{zz}},a_9=-\frac{K_{ftx}}{m}\\ a_{10}=-\frac{K_{fty}}{m},a_{11}=-\frac{K_{ftz}}{m}\end{array}$ (16)

Taking into account relations (15) and (16), on Equation (17):

$\{\begin{array}{l}{\stackrel{\dot{}}{x}}_1=x_2\\ {\stackrel{\dot{}}{x}}_2=a_1{x}_4{x}_6+a_2{x}_2^2+a_3{\Omega }_r{x}_4+b_1{U}_2\\ {\stackrel{\dot{}}{x}}_3=x_4\\ {\stackrel{\dot{}}{x}}_4=a_4{x}_2{x}_6+a_5{x}_4^2+a_6{\Omega }_r{x}_2+b_2{U}_3\\ {\stackrel{\dot{}}{x}}_5=x_6\\ {\stackrel{\dot{}}{x}}_6=a_7{x}_2{x}_4+a_8{x}_6^2+b_3{U}_4\\ {\stackrel{\dot{}}{x}}_7=x_8\\ {\stackrel{\dot{}}{x}}_8=a_9{x}_8+\frac{1}{m}{U}_x{U}_1\\ {\stackrel{\dot{}}{x}}_9=x_{10}\\ {\stackrel{\dot{}}{x}}_{10}=a_{10}{x}_{10}+\frac{1}{m}{U}_y{U}_1\\ {\stackrel{\dot{}}{x}}_{11}=x_{12}\\ {\stackrel{\dot{}}{x}}_{12}=a_{11}{x}_{12}+\frac{\cos \varphi \cos \theta }{m}{U}_1-g\end{array}$

We can derive a simplified state representation of our state system as follows:

$\stackrel{\dot{}}{X}=f\left(X,U\right)=\left\{\begin{array}{c}\stackrel{\dot{}}{\varphi }\\ \stackrel{\dot{}}{\theta }\stackrel{\dot{}}{\psi }{a}_1+{\stackrel{\dot{}}{\varphi }}^2{a}_2+\stackrel{\dot{}}{\theta }{a}_3{\Omega }_r+b_1{U}_2\\ \stackrel{\dot{}}{\theta }\\ \stackrel{\dot{}}{\varphi }\stackrel{\dot{}}{\psi }{a}_4+{\stackrel{\dot{}}{\theta }}^2{a}_5+\stackrel{\dot{}}{\varphi }{a}_6{\Omega }_r+b_2{U}_3\\ \stackrel{\dot{}}{\psi }\\ \stackrel{\dot{}}{\varphi }\stackrel{\dot{}}{\theta }{a}_7+{\stackrel{\dot{}}{\psi }}^2{a}_8+b_3{U}_4\\ \stackrel{\dot{}}{x}\\ \stackrel{\dot{}}{x}{a}_9+\frac{1}{m}{U}_x{U}_1\\ \stackrel{\dot{}}{y}\\ \stackrel{\dot{}}{y}{a}_{10}+\frac{1}{m}{U}_y{U}_1\\ \stackrel{\dot{}}{z}\\ \stackrel{\dot{}}{z}{a}_{11}+\frac{\cos \varphi \cos \theta }{m}{U}_1-g\end{array}\right\}$

3. Return Control

The Backstepping control algorithm is synthesized to force the system to follow a desired trajectory.

The Backstepping command controls orientations (ϕ, θ, ψ) and positions (x, y, z). This design is based on the derived state vector (17) and relying on Lyapunov’s theory (particularly Lyapunov’s second method).

We first consider the command input for the angular rotations subsystem, then the position command input will be derived. On will deal for the subsystem of angular rotations: the roll controls U₂, pitch U₃ and yaw U₄; with regard to the translations subsystem, we will highlight the commands for position in x U_x, position in y U_y and position z U₁.

$\{\begin{array}{l}{\stackrel{\dot{}}{x}}_1=x_2\\ {\stackrel{\dot{}}{x}}_2=a_1{x}_4{x}_6+a_2{x}_2^2+a_3{\Omega }_r{x}_4+b_1{U}_2\end{array}$

The synthesis of the roll command U₂ is done in two stages:

Step 1:

We take the first equation: ${\stackrel{\dot{}}{x}}_1=x_2$

We define the error e1 between $\varphi$ et ${\varphi }_d$ by:

$e_1=\varphi -{\varphi }_d=x_1-x_{1d}$

whose dynamics can be derived as follows:

${\stackrel{\dot{}}{e}}_1={\stackrel{\dot{}}{x}}_1-{\stackrel{\dot{}}{x}}_{1d}=x_2-{\stackrel{\dot{}}{x}}_{1d}$

This definition explicitly indicates our command objective: the error e₁ must asymptotically converge towards zero.

We choose the first candidate Lyapunov function of the following form:

$V\left(e_1\right)=\frac{1}{2}{e}_1^2$

The calculation of the derivative of the Lyapunov function along the trajectory is solved as follows:

$\stackrel{\dot{}}{V}\left(e_1\right)=e_1{\stackrel{\dot{}}{e}}_1=e_1\left(x_2-{\stackrel{\dot{}}{x}}_{1d}\right)$

To ensure stability, it is necessary that $\stackrel{\dot{}}{V}\left(e_1\right)\le 0$ ; for that we take as a virtual command ${\tilde{x}}_2$ avec ${\tilde{x}}_2={\stackrel{\dot{}}{x}}_{1d}-k_1{e}_1$ avec $k_1>0$

This choice makes it possible to obtain:

$\stackrel{\dot{}}{V}\left(e_1\right)=e_1\left(-k_1{e}_1\right)=-k_1{e}_1^2$

Step 2: As the virtual control cannot instantly take its desired value, we seek in what follows to stabilize the error between the virtual control and the stabilizing function. The second equation is defined by:

${\stackrel{\dot{}}{x}}_2=a_1{x}_4{x}_6+a_2{x}_2^2+a_3{\Omega }_r{x}_4+b_1{U}_2$

The new error variable between the virtual control and the stabilizing function is given by:

$e_2=x_2-{\tilde{x}}_2$

$e_2=x_2-{\stackrel{\dot{}}{x}}_{1d}+k_1{e}_1$

By linear combination, on a:

${\stackrel{\dot{}}{e}}_1=e_2-k_1{e}_1$

And by derivation:

${\stackrel{\dot{}}{e}}_2={\stackrel{\dot{}}{x}}_2-{\ddot{x}}_{1d}+k_1\left(e_2-k_1{e}_1\right)$

The new Lyapunov function of the augmented system is ready as follows:

$V\left(e_1,e_2\right)=\frac{1}{2}{e}_1^2+\frac{1}{2}{e}_2^2$

Its derivative is given by:

$\stackrel{\dot{}}{V}\left(e_1,e_2\right)=e_1\left(e_2-k_1{e}_1\right)+e_2\left({\stackrel{\dot{}}{x}}_2-{\ddot{x}}_{1d}+k_1\left(e_2-k_1{e}_1\right)\right)$

The choice of the control law must lead to a negative derivative of the Lyapunov function $\stackrel{\dot{}}{V}\left(e_1,e_2\right)\le 0$ as following:

$\stackrel{\dot{}}{V}\left(e_1,e_2\right)=-k_1{e}_1^2-k_2{e}_2^2$ avec $\left(k_1,k_2\right)>0$

The desirable dynamic for error at this stage is defined as:

${\stackrel{\dot{}}{e}}_2=-k_2{e}_2-e_1$

We find:

$-k_2{e}_2-e_1={\stackrel{\dot{}}{x}}_2-{\ddot{x}}_{1d}+k_1\left(e_2-k_1{e}_1\right)$ (18)

By replacing with the expressions of (17) in (18):

$U_2=\frac{1}{b_1}\left[{\ddot{x}}_{1d}-e_1-k_2{e}_2-k_1\left(e_2-k_1{e}_1\right)-a_1{x}_4{x}_6-a_2{x}_2^2-a_3{\Omega }_r{x}_4\right]$

$U_3,U_4,U_1,U_x,U_y$ are calculated in the same way.

$U_3=\frac{1}{b_2}\left[{\ddot{x}}_{3d}-e_3-k_4{e}_4-k_3\left(e_4-k_3{e}_3\right)-a_4{x}_4{x}_2-a_5{x}_2^2-a_6{\Omega }_r{x}_2\right]$

$U_4=\frac{1}{b_3}\left[{\ddot{x}}_{5d}-e_5-k_6{e}_6-k_5\left(e_6-k_5{e}_5\right)-a_7{x}_2{x}_4-a_8{x}_6^2\right]$

$U_1=\frac{m}{\cos x_1\cos x_3}\left[{\ddot{x}}_{11d}+g-e_{11}-k_{12}{e}_{12}-k_{11}\left(e_{12}-k_{11}{e}_{11}\right)-a_{11}{x}_{12}\right]$

$U_x=\frac{m}{U_1}\left[{\ddot{x}}_{7d}-e_7-k_8{e}_8-k_7\left(e_8-k_7{e}_7\right)-a_9{x}_8\right]$

$U_y=\frac{m}{U_1}\left[{\ddot{x}}_{9d}-e_9-k_{10}{e}_{10}-k_9\left(e_{10}-k_9{e}_9\right)-a_{10}{x}_{10}\right]$

4. Simulation Results

The six (6) control laws acquired in section II will be implemented in the Matlab/Simulink environment with the parameters of Table 1. Through these commands, it is a question of orienting and stabilizing the quadcopter to the desired values: $\left[\begin{array}{ccc}{\varphi }_d& {\theta }_d& {\psi }_d\end{array}\right]=\left[\begin{array}{ccc}0.2& 0.2& 0.2\end{array}\right]$ (en radian) $\left[\begin{array}{ccc}{x}_d& y_d& z_d\end{array}\right]=\left[\begin{array}{ccc}20& 20& 20\end{array}\right]$ (en mètre).

Table 1. Physical parameters.

The results of the simulations are shown in Figures 3-5. They are distributed as follows:

Figure 4 presents the simulations of the Backstepping command applied to roll, pitch, yaw and the three positions (x, y, z), without disturbances. It can be seen that the Backstepping control correctly orientates and stabilizes the quadrirotor at the desired values. Ease of regaining gains favoring the stabilization of the quadrirotor.

Figure 5 shows the influence of disturbances on the Backstepping command. Using the Beaufort scale, we applied wind frequencies (frequency = 0.049 Hz) corresponding to winds of at least 22 m/s to our dynamics. There is a slight phase shift between the measurement and the desired destination for roll, pitch and yaw: the quadrirotor tends asymptotically towards the reference. On the other hand, despite this significant disturbance, the air vehicle awaits its altitude z and the x and y positions with stabilization at the desired value.

To validate the robustness of the backstepping control to a quadrirotor with six (6) degrees of freedom and to consolidate our work, we compared this control to the PID control. The PID controller was designed in the Matlab/Simulink environment and the results can be seen in Figure 5. From Figure 6, the response time of the system to the PID control is quite long compared to the Backstepping control, the difficulty of finding adequate gains allowing the quadrotor to reach the desired value and the stabilizer. This command is satisfactory at altitude.

Figure 3. Simulink model of Backstepping control and Quadrotor dynamics.

Figure 4. Backstepping control of roll, pitch, yaw, x, y, z positions without disturbances.

Figure 5. Backstepping control of roll, pitch, yaw, x, y, z positions with disturbances.

Figure 6. PID control of roll, pitch, yaw, x, y, z positions.

5. Conclusion

In this article, the stabilization in attitude and altitude of a quadrirotor by the Backstepping command has been successfully demonstrated. First, we presented the dynamic model of the quadcopter with six degrees of freedom, taking into account all the aerodynamic and gyroscopic effects. In the second part, we developed the Backstepping control strategy based on the Lyapunov function in order to guarantee the stability of the system. Then, this control algorithm was applied to the quadrirotor. Finally, the simulation results come from a high level of reliability, a very short and precise response time, stabilization with or without turbulence. All these advantages prove the robustness of the Backstepping command.

NOTES

¹c: cosines.

²s: sinus.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Salih, A.-L., Moghavvemi, M., Mohamed, H.A.F. and Gaeid, K.S. (2010) Flight PID Controller Design for a UAV Quadrotor. Scientific Research and Essays, 5, 3360-3367.
[2]	Szafranski, G. and Czyba, R. (2011) Different Approaches of PID Control UAV Type Quadrotor. International Micro Air Vehicles Conference, 70-75.
[3]	García, R.A., et al. (2012) Robust PID Control of the Quadrotor Helicopter. IFAC Proceedings Volumes, 45, 229-334. https://doi.org/10.3182/20120328-3-IT-3014.00039
[4]	Yassine, J. (2016) Commande Non Linéaire Hiérarchique d’un drone de type quadrirotor sans mesure de la vitesse linéaire. Ecole de Technologie Supérieure, Université du Québec.
[5]	Belobo Mevo, B. (2019) Contribution à la commande adaptative et robuste d’un robot mobile de type unicycle avec modèle non linéaire. Université du Québec en Abitibi Témiscamingue.
[6]	Brossard, J., Bensoussan, D., Landry, R. and Hammami, M. (2019) A New Fast Compensator Design Applied to a Quadcopter. Proceedings of the 8th International Conference on Systems and Control, Marrakech, 23-25 October 2019, 238-247. https://doi.org/10.1109/ICSC47195.2019.8950601
[7]	Basri, M.A.M., Abidin, M.S.Z. and Subha, N.A.M. (2018) Simulation of Backstepping-Based Nonlinear Control for Quadrotor Helicopter. Applications of Modelling and Simulation, 2, 34-40.
[8]	Djamel, K., Abdellah, M. and Benallegue, A. (2016) Attitude Optimal Backstepping Controller Based Quaternion for a UAV. Mathematical Problems in Engineering, 2016, Article ID: 8573235. https://doi.org/10.1155/2016/8573235
[9]	Bellahcene, Z., Bouhamida, M., Benghanem, M. and Laidani, A. (2012) La Commande Intégrale Backstepping Appliquée à un Hélicoptère à quatre hélices. 2ème CIMGLE’ 2012-19/21 Novembre 2012, Oran-ALGERIE.
[10]	Babaei, R. and Ehyaei, A.F. (2015) Robust Backstepping Control of a Quadrotor UAV Using Extended Kalman Bucy Filter. International Journal of Mechatronics, Electrical and Computer Technology, 5, 2276-2291.
[11]	Rodriguez, H.R., Vega, V.P., Orta, A.S. and Salazar, O.G. (2014) Robust Backstepping Control Based on Integral Sliding Modes for Tracking of Quadrotors. Journal of Intelligent & Robotic Systems, 73, 51-66. https://doi.org/10.1007/s10846-013-9909-4
[12]	Chen, F.Y., Jiang, R.Q., Zhang, K.K. and Jiang, B. (2016) Robust Backstepping Sliding Mode Control and Observer-Based Fault Estimation for a Quadrotor UAV. IEEE Transactions on Industrial Electronics, 63, 5044-5056.
[13]	Nguyen, X.-M. and Hong, S.K. (2019) Robust Backstepping Trajectory Tracking Control of a Quadrotor with Input Saturation via Extended State Observer. Applied Sciences, 9, 5184. https://doi.org/10.3390/app9235184