2.1 Dynamic model of the disturbed attitude and translation of the unconventional quadrotor

The foldable drone can change the angles $\sigma_i(t)\vert_{i=1,...,4}$ between its arms via the servomotors, these latter are linked to a central body as shown in Fig. 1. Each arm contains a rotor placed at its extremity. Changing the angular velocity of the rotors $\omega_i(t)\vert_{i=1,...,4}$ and rotating the arms can change the shape of the foldable drone and produce other more complex configurations from the basic X configuration.

In this paper, we have chosen six different configurations of the quadrotor with rotating arms in this order (X, Y, YI, H, T, O) as shown in Fig. 2 and Table 1.

Figure 2. Different configurations of the quadrotor with rotating arms: (a) X configuration, (b) Y configuration, (c) YI configuration, (d) H configuration, (e) T configuration, (f) O configuration.

Table 1. The arm angles of each configuration

The quadrotor with rotating arms is considered as a rigid body of which its motion in space can be represented by two coordinate frames: $B=(o,x_m,y_m,z_m)$ is the body-fixed frame attached to the centre of gravity of the aerial vehicle and $E=(o_i,x_i,y_i,z_i)$ is the earth-fixed frame as shown in Fig. 3.

Hypothesis

The following basic hypothesis are considered in quadrotor with rotating arms modeling:

• The quadrotor is a rigid body whose mass is constant.

• The propellers are rigid.

• The transition from one configuration to another is long so as not to generate abrupt phenomena.

• The angle of manoeuverability of the arms is assumed to be small.

• The rotation velocities of the arms are not too high.

We elaborate the dynamic equations of the system by using the Newton-Euler formalism.

The system can be written in the following form [Reference Derrouaoui, Bouzid, Guiatni, Kada, Dib and Moudjari4]:

(1) \begin{equation}\small{\begin{cases}\dot{\xi} = V \\[5pt]m\ddot{\xi} = R_{\eta}(F_t-F_g) + F_d + F_{ext} \\[5pt]\dot{R}_\eta=R_{\eta} \, Sk(\Omega) \\[5pt]J(\sigma(t))\dot{\Omega}=-Sk(\Omega)J(\sigma(t)) \Omega + \Gamma_f - \Gamma_a - \Gamma_{gy} + \Gamma_{ext}\end{cases},}\end{equation}

where $\xi \in \mathbb{R}^{3}$ and $V \in \mathbb{R}^{3}$ are respectively the position and velocity of the quadrotor with rotating arms. m is the total mass of the quadrotor, which is supposed to be constant and $\Omega \in \mathbb{R}^{3}$ its angular velocity expressed in the body frame.

$F_t$ , $F_g$ and $F_d$ belong to $\mathbb{R}^{3}$ are respectively the total force generated by the four rotors, the gravity force along the z-axis and the drag force along the three axes (x,y,z).

(2) \begin{equation} F_t=R \begin{bmatrix} 0\;\;\;\; &\quad 0\;\;\;\; & \quad \sum _{i=1}^4 F_i \end{bmatrix}^T, \end{equation}

where

(3) \begin{equation} F_i = b \,w_i^2, \end{equation}

(4) \begin{equation} F_g= \begin{bmatrix} 0\;\;\;\; &\quad 0\;\;\;\; & \quad -mg \end{bmatrix}^T, \end{equation}

(5) \begin{equation} F_d=-\text{diag} \begin{pmatrix} K_{dfx}\;\;\;\; &\quad K_{dfy}\;\;\;\; & \quad K_{dfz} \end{pmatrix}\dot{\xi}, \end{equation}

with b is the aerodynamic coefficients and $K_{df(x,y,z)}$ are the translation drag coefficients.

$F_{ext} \in \mathbb{R}^{3}$ represents the disturbances according to the translation movement.

Figure 3. Schematic of the quadrotor with rotating arms.

$\Gamma_f$ , $\Gamma_{gy}$ and $\Gamma_a$ belong to $\mathbb{R}^{3}$ are the moment caused by the thrust and drag forces, the gyroscopic torque of the propellers and the moment resulting from aerodynamic friction.

(6) \begin{equation}\Gamma_f= \begin{bmatrix} \Gamma_x\;\;\;\; & \quad \Gamma_y\;\;\;\; & \quad \Gamma_z \end{bmatrix}^T,\end{equation}

(7) \begin{equation} \Gamma_{gy}= \sum _{i=1}^4 \Omega \wedge J_r \begin{bmatrix} 0\;\;\;\;&\quad 0\;\;\;\;&\quad (-1)^{i+1}w_i \end{bmatrix}^T, \end{equation}

$J_r$ is the rotor inertia

(8) \begin{equation} \Gamma_a=\text{diag} \begin{pmatrix} K_{afx}\;\;\;\; & \quad K_{afy}\;\;\;\; & \quad K_{afz} \end{pmatrix}^T \Omega^2, \end{equation}

$\Gamma_{ext}\in \mathbb{R}^{3}$ is the disturbances according to the rotation movement.

$R_\eta \in \mathbb{R}^{3\times3}$ is the rotation matrix connecting the two frames of reference, it can be parameterised by Euler angles $\eta=(\phi,\theta,\psi)$ :

(9) \begin{equation}R_{\eta}=\begin{bmatrix} c\psi c\theta\;\;\;\; &\quad s\phi s\theta c\psi-s\psi c\phi\;\;\;\; &\quad c\phi s \theta c \psi +s\psi s\phi \\[5pt] s \psi c\theta\;\;\;\; & \quad s\phi s \theta s \psi + c\psi c\phi\;\;\;\; & \quad c \phi s \theta s \psi - s \phi c \psi \\[5pt] - s \theta\;\;\;\; &\quad s\phi c \theta\;\;\;\; & \quad c \phi c \theta \end{bmatrix},\end{equation}

$Sk \in \mathbb{R}^{3\times3}$ is the skew-symmetric matrix related to the angular velocity $\Omega$

(10) \begin{equation}Sk=\begin{bmatrix} 0\;\;\;\; & \quad -r\;\;\;\; & \quad q \\[5pt] r\;\;\;\; & \quad 0\;\;\;\; & \quad -p \\[5pt] -q\;\;\;\; & \quad p\;\;\;\; & \quad 0 \\[5pt] \end{bmatrix}.\end{equation}

The reconfiguration of the quadrotor, during the flight, generates a displacement of its CoG, a change of its inertia $J_{3\times3}(\sigma(t))$ and requires a readjustment of the allocation matrix. The calculation of the inertia matrix is mainly based on the theorem of Huygens-Steiner, the analytical details of calculation is presented in the works of Derrouaoui et al. [Reference Derrouaoui, Bouzid and Guiatni31, Reference Derrouaoui, Guiatni, Bouzid, Dib and Moudjari32], in which some approximations are made to simplify the development.

2.2 Developed model

Dynamic model (1) can be written in terms of position, rotation and input as follow:

(11) \begin{equation}\begin{cases}\ddot{\phi} = f_{\phi}(t)+g_{\phi}(t)U_2 \\[5pt]\ddot{\theta} =f_{\theta}(t)+g_{\theta}(t)U_3\\[5pt]\ddot{\psi} =f_{\psi}(t)+g_{\psi}(t)U_4\\[5pt]\ddot{x}=f_x(t)+g_x(t)U_x U_1\\[5pt]\ddot{y}= f_y(t)+g_y(t)U_y U_1\\[5pt]\dot{z}=f_z(t)+g_z(t)U_1\end{cases},\end{equation}

with

(12) \begin{equation}\begin{cases}f_{\phi}(t)=\dfrac{J_{yy}(\sigma(t))-J_{zz}(\sigma(t))}{J_{xx}(\sigma(t))}\dot{\theta}\dot{\psi}-\dfrac{J_r\;\omega_r}{J_{xx}(\sigma(t))}\dot{\theta}-\dfrac{K_{afx}}{J_{xx}(\sigma(t))} \dot{\phi}^2 \\[10pt] g_{\phi}(t)=\dfrac{1}{J_{xx}(\sigma(t))} \\[5pt] \\ f_{\theta}(t)=\dfrac{J_{zz}(\sigma(t))-J_{xx}(\sigma(t))}{J_{yy}(\sigma(t))}\dot{\phi}\dot{\psi}+\dfrac{J_r\;\omega_r}{J_{yy}(\sigma(t))}\dot{\phi}-\dfrac{K_{afy}}{J_{yy}(\sigma(t))}\dot{\theta}^2 \hfill \\[10pt] g_{\theta}(t)=\dfrac{1}{J_{yy}(\sigma(t))}\\[5pt] \\ f_{\psi}(t)=\dfrac{J_{xx}(\sigma(t))-J_{yy}(\sigma(t))}{J_{zz}(\sigma(t))}\dot{\phi}\dot{\theta}-\dfrac{K_{afz}}{J_{zz}(\sigma(t))}\dot{\psi}^2 \hfill\quad \\[10pt] g_{\psi}(t)=\dfrac{1}{J_{zz}(\sigma(t))}\\[5pt] \\ f_x(t)=-\dfrac{K_{dfx}}{m}\dot{x} \hfill \quad \quad \quad \\[10pt] g_x(t)=\dfrac{1}{m}\\[5pt] \\ f_y(t)=-\dfrac{K_{dfy}}{m}\dot{y} \hfill \quad \quad \quad \\[10pt] g_y(t)=\dfrac{1}{m}\\[5pt] \\ f_z(t)=-\dfrac{K_{dfz}}{m}\dot{z}-g \hfill \quad \quad \quad \\[10pt] g_z(t)=\dfrac{\text{cos} \phi\: \text{cos}\theta}{m}\end{cases},\end{equation}

and

(13) \begin{equation}\begin{cases}U_x=(c\phi s\theta c\psi + s\phi s\psi)\\[5pt]U_y=(c\phi s\theta s\psi - s\phi c\psi)\\[5pt]\omega_r=\omega_1 - \omega_2 +\omega_3 - \omega_4\end{cases}.\end{equation}

The model given by (11) can be generalised for two subsystems, one is designed for the position $\Sigma_P$ and the other is for the attitude $\Sigma_A$ in which the translational and rotational dynamics are considered as disturbed non-linear second-order systems. The following subsystems will be used in the control design:

(14) \begin{equation}\Sigma_P:\begin{cases}\dot{x}_{1} = x_{2}\\[5pt]\dot{x}_{2} = f_1(x)+g_1(x)T, \\[5pt]\end{cases}\end{equation}

(15) \begin{equation}\Sigma_A:\begin{cases}\dot{x}_{3} = x_{4}\\[5pt]\dot{x}_{4} = f_2(x)+g_2(x)U, \\[5pt]\end{cases}\end{equation}

with $f_1=\begin{bmatrix} f_x & \quad f_y & \quad f_z \end{bmatrix}^T$ is written as:

(16) \begin{equation} f_1(x)=f_{01}+ d_1, \end{equation}

and $f_2= \begin{bmatrix} f_{\phi} & \quad f_{{\theta}} & \quad f_{{\psi}} \end{bmatrix}^T$ is written as:

(17) \begin{equation}f_2(x)=f_{02}+ d_2,\end{equation}

where $x_{\xi}=\begin{bmatrix} x_1 & \quad x_2\end{bmatrix} \in \mathbb{R}^{3 \times 2} $ is the vector of states for position, $x_1 = \begin{bmatrix} x & \quad y & \quad z\end{bmatrix}^T \in \mathbb{R}^{3}$ , $ x_2 = \begin{bmatrix} \dot{x} & \quad \dot{y} & \quad \dot{z}\end{bmatrix}^T \in \mathbb{R}^{3}$ , $T = \begin{bmatrix}U_x & \quad U_y & \quad U_1 \end{bmatrix}^T$ is the vector of control inputs and the bounded disturbance $d_1= \begin{bmatrix} d_x & \quad d_y & \quad d_z \end{bmatrix} \in \mathbb{R}^{3} $

$x_{\eta}=\begin{bmatrix} x_3 & \quad x_4\end{bmatrix} \in \mathbb{R}^{3 \times 2} $ is the vector of states for attitude, $x_3 = \begin{bmatrix} \phi & \quad \theta & \quad \psi\end{bmatrix}^T \in \mathbb{R}^{3} $ and $ x_4 = \begin{bmatrix} \dot{\phi} & \quad \dot{\theta} & \quad \dot{\psi}\end{bmatrix}^T \in \mathbb{R}^{3} $ , $U = \begin{bmatrix}U_2 & \quad U_3 & \quad U_4 \end{bmatrix}^T$ is the vector of control inputs and the bounded disturbance $d_2= \begin{bmatrix} d_{\phi} & \quad d_{\theta} & \quad d_{\psi} \end{bmatrix} \in \mathbb{R}^{3} $ . $g_1(x)$ and $g_2(x))$ belong to $\mathbb{R}^{3 \times 3}_*$ . The functions $f_{01}$ and $f_{02}$ represent the known terms of $f_1(x)$ and $f_2(x)$ respectively.

2.3 Control objective

The main objective is to design a control laws T and U such that the variable ( $x_{\xi}(t), x_{\eta}(t)$ ) track robustly the desired signals ( $x_{\xi d}(t)$ , $x_{\eta d}(t)$ ), i.e. to render the tracking error

$\begin{pmatrix}\varrho_1(t) \\[5pt] \varrho_3(t)\end{pmatrix}$ $\equiv$ $ \begin{pmatrix} x_{\xi}(t) - x_{\xi d}(t) \\[5pt] x_{\eta}(t) - x_{\eta d}(t)\end{pmatrix}$ $\equiv$ $\begin{pmatrix}0 \\[5pt] 0\end{pmatrix} $ in finite time $t_R$ .

3.1 Backstepping design

The backstepping control technique performs successive relations to make the looped system equivalent to stable cascaded first-order sub-systems in the Lyapunov sense, thus building a control law and a stabilising Lyapunov function. This method begins with the definition of the Lyapunov function. The principle is to construct at each step a gap between the current state of the system and a virtual controller that would guarantee the negativity of the derivative of the Lyapunov function. At the end of each step, the Lyapunov function compensates this gap, to ensure the convergence to zero. The process ends when the real control of the system appears, guaranteeing the convergence of the real states to the desired states.

The control laws by backstepping combined with the ST algorithm are designed as:

(18) \begin{equation}T=g_1^{-1}(x)\;(-f_{01}(x)+\dot{\upsilon}_{1}+\upsilon_{2}),\;\end{equation}

with

(19) \begin{equation} \begin{cases}\upsilon_{1}=\dot{x}_{1d}-\kappa_1\;\varrho_{1} \\[5pt]\upsilon_2=-\lambda _1\;\mathfrak{Sig}^{\frac{1}{2}}(\varrho_2)-\lambda _2\;\mathfrak{Int}(\varrho_2),\end{cases}\end{equation}

(20) \begin{equation}U=g_2^{-1}(x)\;(-f_{02}(x)+\dot{\upsilon}_{3}+\upsilon_{4}),\;\end{equation}

with

(21) \begin{equation}\begin{cases}\upsilon_3=\dot{x}_{3d}-\kappa_2\;\varrho_{3} \\[5pt]\upsilon_4=-\lambda _3\;\mathfrak{Sig}^{\frac{1}{2}}(\varrho_4)-\lambda _4\;\mathfrak{Int}(\varrho_4),\end{cases}\end{equation}

where $\left(\kappa_1=\text{diag} \begin{bmatrix} \kappa_x\;\;\;\; & \quad \kappa_y\;\;\;\; & \quad \kappa_z \end{bmatrix}, \kappa_2 =\text{diag} \begin{bmatrix} \kappa_{\phi}\;\;\;\; &\quad \kappa_{\theta}\;\;\;\; & \quad \kappa_{\psi} \end{bmatrix}\right) \in \mathbb{R}^{3 \times 3}_{+*}$

and $ \lambda_1 = \text{diag} \begin{bmatrix} \lambda_{1x}\;\;\;\; & \quad \lambda_{1y}\;\;\;\; & \quad \lambda_{1z} \end{bmatrix}, \lambda_2 = \text{diag} \begin{bmatrix} \lambda_{2x}\;\;\;\; & \lambda_{2y}\;\;\;\; & \quad \lambda_{2z} \end{bmatrix},\\[5pt] \lambda_3 = \text{diag} \begin{bmatrix} \lambda_{1\phi}\;\;\;\; & \quad \lambda_{1\theta}\;\;\;\; & \quad \lambda_{1\psi} \end{bmatrix}, \lambda_4 = \text{diag} \begin{bmatrix} \lambda_{2\phi}\;\;\;\; & \quad \lambda_{2\theta}\;\;\;\; & \quad \lambda_{2\psi} \end{bmatrix} $

The following Lemma is used for the stability analysis

Lemma 1. [Reference Bhat and Bernstein33, Reference Bhat and Bernstein34] Considering the continuous system

(22) \begin{equation}\dot{x}=f(x),\quad f(0)=0,\quad x \in \mathbb{R}^n.\end{equation}

We suppose that there is a continuous positive definite function V: $ \mathbb{R}^n \rightarrow \mathbb{R}^+,$ a real number $ \gamma > 0$ , $0 < \alpha < 1$ and an open $ U \subseteq \mathbb{R}^n$ in the neighbourhood of the origin so that the following inequalities are satisfied:

(23) \begin{equation}\dot{V}(x)+ \gamma\;V^\alpha(x)\leq 0, \quad x \in V\setminus\lbrace0\rbrace.\end{equation}

Then, the origin of the system (22) is a stable equilibrium in finite time. If $ U =\mathbb {R}^n $ , the origin is a globally stable equilibrium in finite time.

Moreover, the settling time is given by:

(24) \begin{equation}t_R \leq \frac{1}{\gamma(1-\alpha)} V^{1-\alpha}.\end{equation}

Assumption 1. The perturbation terms $d_{1}$ and $d_{2}$ are assumed to be bounded, such that:

(25) \begin{equation} \dot{ \parallel d_1 \parallel} < h_{1{max} } \quad \quad and \quad \quad \dot{\parallel d_2 \parallel} < h_{2{max} }, \end{equation}

where $h_{1{max} } $ and $h_{2{max} }$ are strictly positive constants supposed to be known.

3.1.1 First step

Consider system (14)

Let us define the tracking error $\varrho_1$ :

(26) \begin{equation}\varrho_{1}= x_{1}-x_{1d}.\end{equation}

The first Lyapunov candidate function $\mathcal{L}_1$ , associated with $\varrho_1$ is given by:

(27) \begin{equation}\mathcal{L}_1=\frac{1}{2}\varrho_1^T\varrho_1.\end{equation}

Its temporal derivative:

(28) \begin{equation}\dot{\mathcal{L}}_1=\varrho_{1}^T\;\dot{\varrho}_{1}=\varrho^T_{1}(\dot{x}_1-\dot{x}_{1d})=\varrho_{1}^T\left(\dot{x}_1-\upsilon_{1}\right)+\varrho_{1}^T\left(\upsilon_{1}-\dot{x}_{1d}\right).\end{equation}

The term $(\dot{x}_{1d})$ represents the velocity set point and $\upsilon_1$ is the stabilising controller of $\varrho_1$ .

(29) \begin{equation}\upsilon_{1}=\dot{x}_{1d}-\kappa_1\;\varrho_{1}.\end{equation}

Replacing (29) in (28), we obtain:

(30) \begin{equation}\dot{\mathcal{L}}_1=-\varrho_{1}^T\;\kappa_{1}\;\varrho_{1}+\varrho_{1}^T (\dot{x}_1-\upsilon_{1}).\end{equation}

3.1.2 Second step

The error $\varrho_2$ is defined as the difference between the real velocity and the virtual velocity control. The convergence of $\varrho_2$ will naturally lead to the convergence of $\varrho_1$ , since the velocity will tend towards the virtual velocity and the derivative of $\mathcal{L}_1$ will be strictly negative.

Let us define the error $\varrho_2$ :

(31) \begin{equation}\varrho_{2}=\dot{x}_1-\upsilon_{1}.\end{equation}

We define the second Lyapunov candidate function $\mathcal{L}_2$ , associated with $\varrho_2$ .

(32) \begin{equation}\mathcal{L}_2=\frac{1}{2}\;\varrho_2^T \varrho_2.\end{equation}

The derivative of the Lyapunov function is given by:

(33) \begin{equation}\dot{\mathcal{L}}_2=\varrho_{2}^T\;\dot{\varrho}_{2}=\varrho_{2}^T\;\left(f(x)+g(x)\;T-\dot{\upsilon}_{1}\right).\end{equation}

The control law is given by:

(34) \begin{equation}T=g^{-1}(x)\;\left(-f_0(x)+\dot{\upsilon}_{1}+\upsilon_{2}\right).\;\end{equation}

Replacing (34) in (33), we obtain

(35) \begin{equation}\dot{\mathcal{L}}_2=\varrho_2^T\;\left(d_{1} +\upsilon_2\right).\end{equation}

We take $\upsilon_2$ as

(36) \begin{equation}\upsilon_2=-\lambda _{1}\;\mathfrak{Sig}^{\frac{1}{2}}(\varrho_2)-\lambda _{2}\;\mathfrak{Int}(\varrho_2).\end{equation}

We replace (36) in (35), we find

(37) \begin{equation}\dot{\mathcal{L}}_2=\varrho_2^T\left(-\lambda _{1}\;\mathfrak{Sig}^{\frac{1}{2}}(\varrho_2)-\lambda _{2}\;\mathfrak{Int}(\varrho_2) + d_{1}\right).\end{equation}

We make the following change of variable

(38) \begin{equation}z=\begin{bmatrix}z_1 & \quad z_2\end{bmatrix}^T,\end{equation}

(39) \begin{equation} \begin{cases} z_{1}=\varrho_{2}\\[5pt]z_{2}=-\lambda _{2}\;\mathfrak{Int}(\varrho_2)+d_{1}\end{cases}.\end{equation}

Which implies

(40) \begin{equation}\begin{cases}\dot{z}_{1}=-\lambda _{1}\;\mathfrak{Sig}^{\frac{1}{2}}(z_1)+z_{2}\\[5pt]\dot{z}_{2}=-\lambda _{2}\;\dot{\mathfrak{Int}}(z_1)+\dot{d}_{1}\end{cases}.\end{equation}

The convergence of $z_1$ and $z_2$ to zero ensures the convergence of $\varrho_2$ and $\dot{\varrho}_2$ to zero.

We choose now the following Lyapunov function to analyse the stability:

(41) \begin{equation}V(z)=v^T P\;v.\end{equation}

The vector v and the symmetric matrix $P \in \mathbb{R}^{2\times2}$ are given by:

(42) \begin{equation} v=\begin{bmatrix}\;\mathfrak{Sig}^{\frac{1}{2}}(z_1)& \quad z_2\end{bmatrix}^T\quad \text{and} \quad P=\frac{1}{2}\begin{bmatrix}4 \lambda_2+\lambda_1^2 & \quad -\lambda_1\\[5pt]-\lambda_1 & \quad 2\end{bmatrix},\end{equation}

this function is defined positive if $\lambda_2$ is defined positive.

(43) \begin{equation}\delta_{{min} } (P)\Vert v \Vert^2 \;\leq\; V(z) \;\leq \; \delta_{{max} }(P)\Vert v \Vert^2,\end{equation}

where v is given by:

(44) \begin{equation} \Vert v \Vert^2 =\Vert z_1 \Vert+z_2^2, \end{equation}

$\delta_{{min}} (P)$ and $\delta_{{max}}(P)$ are respectively the minimal and maximal eigenvalues of P.

Knowing that $\dot{\parallel d_1 \parallel} \leq h_{1{max} }$ (see Assumption 1), the temporal derivative of the candidate Lyapunov function verify that:

(45) \begin{equation}\dot{V} \leq -\frac{1}{\Vert z_1 \Vert^\frac{1}{2}} \;\delta_{{min} } (Q)\;\Vert v \Vert ^2\; + h_{1{max} } \;\Vert \xi \Vert \Vert v \Vert,\end{equation}

where

(46) \begin{equation}Q=\frac{\lambda_1}{2}\begin{bmatrix}2 \lambda_2+\lambda_1^2 & \quad -\lambda_1\\[5pt]-\lambda_1 & \quad 1\end{bmatrix}\quad \quad ;\quad \xi=\begin{bmatrix}-\lambda_1 & \quad 2\end{bmatrix}.\end{equation}

Using (44), we can deduce that

(47) \begin{equation}\frac{1}{\Vert z_1 \Vert^\frac{1}{2}} \geq \frac{1}{\Vert v \Vert}.\end{equation}

We find

(48) \begin{equation}\dot{V} \leq (-\delta_{{min} }(Q)+h_{1{max} } \Vert \xi \Vert)\;\Vert v \Vert.\end{equation}

Using (43), we can see that:

(49) \begin{equation}\frac{V}{\delta_{{max} }(P)}\leq \Vert v \Vert ^2 \leq\frac{V}{\delta_{{min} }(P)}.\end{equation}

We find

(50) \begin{equation}\dot{V} \leq - (\delta_{{min} }(Q)-h_{1{max} } \Vert \xi \Vert)\;\frac{V^\frac{1}{2}}{\sqrt{\delta_{{max} }(P)}}.\end{equation}

To conclude on the negativity of $\dot{V}(z)$ , it is necessary that:

(51) \begin{equation}\delta_{{min} }(Q) \geq h_{1{max} } \Vert \xi \Vert.\end{equation}

The matrix Q will be defined positive with a minimal eigenvalue $ \delta_{{min} }(Q) \geq h_{1{max} } || \xi || $ , $ \lambda _1$ defined positive and $h_{1{max} } > 0$ , hence:

(52) \begin{equation}\lambda_2 \geq \frac{2 \, h_{1{max} }^2 \Vert \xi \Vert^2}{\lambda _1^2 }.\end{equation}

Therefore, we get

(53) \begin{equation}\dot{V}\leq - \gamma\; V ^{\frac{1}{2}}.\end{equation}

with

(54) \begin{equation} \gamma= \frac{\delta_{{min} }(Q)- h_{1{max} } \Vert \xi \Vert}{\sqrt{\delta_{{max} }(P)}}.\end{equation}

The solution is

(55) \begin{equation}\int_0^{t_R} \frac{dV}{V^\frac{1}{2}}\leq -\gamma \int_0^{t_R} dt \Longrightarrow t_R= \frac{2 V^\frac{1}{2}z(0)}{\gamma}.\end{equation}

So, we have demonstrated the stability in finite time. i.e. $z_1 $ and $z_2 $ tend to zero in a finite time.

(56) \begin{equation}z_2=0\;\Longrightarrow \quad d_1=\lambda _2 \;\mathfrak{Int}(z_1),\end{equation}

where $ \lambda _2 \;\mathfrak{Int}(z_1) $ represents the perturbation estimation term.

As soon as $z_1=z_2=0 \in \mathbb{R}_3 $ and according to Equation (56), then:

(57) \begin{equation}\begin{split}\quad\Longrightarrow \dot{\mathcal{L}}_2=\varrho_2^T\left(-\lambda _{1}\;\mathfrak{Sig}^{\frac{1}{2}}(\varrho_2)\right) \\[5pt]\quad \quad \quad \Longrightarrow \dot{\mathcal{L}}_2=- \sum^3_{i=1} \lambda_{1i} \mid \varrho_{2i} \mid \mid \varrho_{2i} \mid ^{\frac{1}{2}} \\[5pt]\quad \Longrightarrow \dot{\mathcal{L}}_2=- \sum^3_{i=1} \lambda_{1i} \mid \varrho_{2i} \mid ^{\frac{3}{2}},\end{split}\end{equation}

we define $\tilde{\lambda}_1={min} \lbrace \lambda_{11}, \lambda_{12}, \lambda_{13}\rbrace $ , we get

(58) \begin{equation}\begin{split}\quad \Longrightarrow \dot{\mathcal{L}}_2 \leq - \tilde{\lambda}_1 \sum^3_{i=1} \mid \varrho_{2i} \mid ^{\frac{3}{2}} \\[5pt]\quad \quad \quad \Longrightarrow \dot{\mathcal{L}}_2 \leq - \tilde{\lambda}_1 \sum^3_{i=1} \left( \frac{1}{2} 2 \mid \varrho_{2i} \mid ^{\frac{3}{2}} \right) \\[5pt]\quad \quad \quad \Longrightarrow \dot{\mathcal{L}}_2 \leq - \tilde{\lambda}_1 2^{\frac{3}{4} } \sum^3_{i=1} \left( \frac{1}{2} \mid \varrho_{2i} \mid ^{\frac{3}{2}}\right).\end{split}\end{equation}

Let $\tilde{\kappa}_1=2^{\frac{3}{4}}\tilde{\lambda}_1$

(59) \begin{equation}\begin{split}\Longrightarrow \dot{\mathcal{L}}_2\leq - \tilde{\kappa}_1 \mathcal{L}_2^{\frac{3}{4}} \\[5pt]\Longrightarrow \dot{\mathcal{L}}_2\leq - \tilde{\kappa}_1 \mathcal{L}_2^{\alpha} \end{split},\end{equation}

with $\alpha = \frac{3}{4} \in\;[0 , 1] $

According to Lemma 1, $\mathcal{L}_2$ is stable, so the error $\varrho_2 \rightarrow 0$ .

This implies that $\dot{\mathcal{L}}_1= - \kappa_1\parallel\varrho_1\parallel^2$ is stable. Then, the error $\varrho_1 \rightarrow 0$ .

Desired Euler angles

The unconventional quadrotor can be separated into two connected subsystems, based on the hierarchical backstepping. According to (13), we can calculate the desired Euler angles $\phi_d$ and $\theta_d$ , which will introduce the next error variable:

(60) \begin{equation}\begin{cases}\phi _d= \sin ^{-1}(s \psi _d U_x - c \psi _d U_y) \\[5pt]\theta _d= \sin ^{-1}\left(\dfrac{c \psi _d U_x +s \psi _d U_y}{c \phi _d}\right)\end{cases}.\end{equation}

3.1.3 Third step

Consider system (15), a new error variable is defined as:

(61) \begin{equation}\varrho_{3}= x_{3}-x_{3d}.\end{equation}

Consider the Lyapunov candidate function $\mathcal{L}_3$ , associated with the error $\varrho_3$ :

(62) \begin{equation}\mathcal{L}_3=\frac{1}{2}\varrho_3^T \varrho_3.\end{equation}

Its derivative with respect to time yields:

(63) \begin{equation}\dot{\mathcal{L}}_3=\varrho_{3}^T\;\dot{\varrho}_{3}=\varrho_{3}^T(\dot{x}_3-\dot{x}_{3d})=\varrho_{3}^T(\dot{x}_3-\upsilon_{3})+\varrho_{3}^T(\upsilon_{3}-\dot{x}_{3d}).\end{equation}

The term $(\dot{x}_{3d})$ represents the velocity set point and $\upsilon_3$ is the stabilising controller of $\varrho_3$ .

(64) \begin{equation}\upsilon_{3}=\dot{x}_{3d}-\kappa_2\;\varrho_{3}.\end{equation}

Replacing (64) in (63), we obtain:

(65) \begin{equation}\dot{\mathcal{L}}_3=-\varrho_{3}^T\;\kappa_{2}\;\varrho_{3}+\varrho_{3}^T(\dot{x}_3-\upsilon_{3}),\end{equation}

the second term will be eliminated in the next step.

3.1.4 Fourth step

Consider the second term $\varrho_{3}(\dot{x}_3-\upsilon_{3})$ to continue the process of backstepping. The error $\varrho_4$ is defined as the difference between the real velocity and the virtual velocity control. The convergence of $\varrho_4$ will naturally lead to the convergence of $\varrho_3$ , since the velocity will tend towards the virtual velocity and the derivative of $\mathcal{L}_3$ will be strictly negative.

Let define the error $\varrho_4$ :

(66) \begin{equation}\varrho_{4}=\dot{x}_3-\upsilon_{3}.\end{equation}

The Lyapunov candidate function $\mathcal{L}_4$ associated with $\varrho_4$ is defined as:

(67) \begin{equation}\mathcal{L}_4=\frac{1}{2}\;\varrho_4^T \varrho_4.\end{equation}

Its temporel derivative is given by:

(68) \begin{equation}\dot{\mathcal{L}}_{4}=\varrho_{4}^T\;\dot{\varrho}_{4}=\varrho_{4}^T\;(f(x)+g(x)\;U-\dot{\upsilon}_{3}).\end{equation}

The controller is given by:

(69) \begin{equation}U=g^{-1}(x)\;(-f_0(x)+\dot{\upsilon}_{3}+\upsilon_{4}).\;\end{equation}

Replacing (69) in (68), we obtain

(70) \begin{equation}\dot{\mathcal{L}}_4=\varrho_4^T\;(d_2 +\upsilon_4).\end{equation}

We take $\upsilon_4$ as

(71) \begin{equation}\upsilon_2=-\lambda _{3}\;\mathfrak{Sig}^{\frac{1}{2}}(\varrho_4)-\lambda _{4}\;\mathfrak{Int}(\varrho_4).\end{equation}

We replace (71) in (70), we find

(72) \begin{equation}\dot{\mathcal{L}}_4=\varrho_4^T\left(-\lambda _{3}\;\mathfrak{Sig}^{\frac{1}{2}}(\varrho_4)-\lambda _{4}\;\mathfrak{Int}(\varrho_4) + d_2\right).\end{equation}

Consider the following change of variable

(73) \begin{equation}\bar{z}=\begin{bmatrix}z_3 & \quad z_4\end{bmatrix}^T,\end{equation}

(74) \begin{equation} \begin{cases} z_3=\varrho_4\\[5pt]z_4=-\lambda _{4}\;\mathfrak{Int}(\varrho_4)+d_2\end{cases}.\end{equation}

Which implies

(75) \begin{equation}\begin{cases}\dot{z}_{1}=-\lambda _{3}\;\mathfrak{Sig}^{\frac{1}{2}}(z_3)+z_{4}\\[5pt]\dot{z}_{2}=-\lambda _{4}\;\dot{\mathfrak{Int}}(z_3)+\dot{d}_{2}\end{cases}.\end{equation}

The convergence of $z_3$ and $z_4$ to zero ensures the convergence of $\varrho_4$ and $\dot{\varrho}_4$ to zero. Therefore, the disturbance $d_2$ will be estimated in finite time.

As soon as $z_3$ = $z_4$ =0 and the disturbance is estimated, then:

(76) \begin{equation}\Longrightarrow \dot{\mathcal{L}}_4\leq - \tilde{\kappa}_2 \mathcal{L}_4^{\alpha},\end{equation}

where $\alpha = \frac{3}{4} $ , $\tilde{\kappa}_2 = 2^{\frac{3}{4}}\tilde{\lambda}_3$ and $\tilde{\lambda}_3={min} \lbrace \lambda_{31}, \lambda_{32}, \lambda_{33}\rbrace $

According to Lemma 1, $\mathcal{L}_4$ is stable $\Rightarrow$ the error $\varrho_4 \rightarrow 0$ .

which implies that $\dot{\mathcal{L}}_3= - \kappa_2\parallel\varrho_3\parallel^2$ is stable. Then, the error $\varrho_3 \rightarrow 0$

3.2 Control architecture

Hierarchical control is known in aeronautics as guidance and piloting control. The system is decomposed into two cascaded subsystems: one for translation corresponding to the slow dynamics called the high level (outer loop) and the other for rotation corresponding to the fast dynamics called the low level (inner loop) (see Fig. 4).

Figure 4. Control architecture.

The reference trajectories $(x_d, y_d, z_d$ and $\psi_d)$ are generated by the guidance system block. $U_1$ controls the attitude while the two other translation variables x and y are controlled indirectly by controlling the two rotation angles $(\phi$ and $\theta)$ with the virtual controllers $U_x$ and $U_y$ . The correction block will generate the desired roll and pitch angles. The foldable drone is controlled by the velocity of the motors. These velocities are obtained from the control matrix. The CoG and the inertia matrix are instantly calculated according to the variation of the angles $\sigma_i$ . The different configurations of the quadrotor are generated by the switching block.

The main difficulty of this control structure is to demonstrate the stability of the global system in a closed loop and to guarantee good performance as well as robust behaviour towards disturbances. Servomotors are controlled by a classical PID to turn the arms, where:

(77) \begin{equation}U_s=k_pe_{\sigma_i}(t)+k_I\int e_{\sigma_i}(t)dt+k_D \dot{e}_{\sigma_i}(t),\end{equation}

with $e_{\sigma_i}(t)$ is the tracking error. $k_p, k_I$ and $k_D$ are the controller’s gains.

The studied quadrotor is controlled by a backstepping controller with a correction term based on the super-twisting algorithm, which is designed to follow the desired trajectory $(x_d, y_d, z_d)$ and the angle $\psi_d$ in the presence of disturbances.

3.3 Control inputs

Applying the control laws (34) and (69) on the systems (14) and (15) respectively, we obtain the final control inputs as:

Altitude control (z)

$\begin{cases}\varrho_{11}=z-z_d\\[5pt]\varrho_{21}=\dot{z}-\upsilon_{11}\\[5pt]\upsilon_{11}=\dot{z}_d-\kappa_z\;\varrho_{11}\\[5pt]\upsilon_{21}=-\lambda _{1z}\;|\varrho_{21}|^{\frac{1}{2}}\; \text{sign} (\varrho_{21})-\lambda _{2z}\;\int \text{sign}(\varrho_{21}(z))dz\end{cases}$

(78) \begin{equation}U_1=\frac{1}{g_z(x)}\left(-f_z(x)+\dot{\upsilon}_{11}+\upsilon_{21}\right).\end{equation}

Roll control ( $\phi$ )

$\begin{cases}\varrho_{12}=\phi -\phi _d\\[5pt]\varrho_{22}=\dot{\phi}-\upsilon_{12}\\[5pt]\upsilon_{12}=\dot{\phi}_d-\kappa_{\phi}\;\varrho_{12}\\[5pt]\upsilon_{22}=-\lambda _{1\phi}\;|\varrho_{22}|^{\frac{1}{2}}\; \text{sign}(\varrho_{22})-\lambda _{2\phi}\;\int \text{sign}(\varrho_{22})(z)dz\end{cases}$

(79) \begin{equation}U_2=\frac{1}{g_{\phi}(x)}\left(-f_{\phi}(x)+\dot{\upsilon}_{12}+\upsilon_{22}\right).\end{equation}

Pitch control ( $\theta$ )

$\begin{cases}\varrho_{13}=\theta -\theta _d\\[5pt]\varrho_{23}=\dot{\theta}-\upsilon_{13}\\[5pt]\upsilon_{13}=\dot{\theta}_d-\kappa_{\theta}\;\varrho_{13}\\[5pt]\upsilon_{23}=-\lambda _{1\theta}\;|\varrho_{23}|^{\frac{1}{2}}\; \text{sign}(\varrho_{23})-\lambda _{2\theta}\;\int \text{sign}(\varrho_{23}(z))dz\end{cases}$

(80) \begin{equation}U_3=\frac{1}{g_{\theta}(x)}\left(-f_{\theta}(x)+\dot{\upsilon}_{13}+\upsilon_{23}\right).\end{equation}

Yaw control ( $\psi$ )

$\begin{cases}\varrho_{14}=\psi -\psi _d\\[5pt]\varrho_{24}=\dot{\psi}-\upsilon_{14}\\[5pt]\upsilon_{14}=\dot{\psi}_d-\kappa_{\psi}\;\varrho_{14}\\[5pt]\upsilon_{24}=-\lambda _{1\psi}\;|\varrho_{24}|^{\frac{1}{2}}\; \text{sign}(\varrho_{24})-\lambda _{2\psi}\;\int \text{sign}(\varrho_{24}(z))dz\end{cases}$

(81) \begin{equation}U_4=\frac{1}{g_{\psi}(x)}(-f_{\psi}(x)+\dot{\upsilon}_{14}+\upsilon_{24}).\end{equation}

Virtual controllers

Control of (x)

\begin{equation*}\begin{cases}\varrho_{1x}=x-x_d\\[5pt]\varrho_{2x}=\dot{x}-v_{1x}\\[5pt]\upsilon_{1x}=\dot{x}_d-\kappa_x\;\varrho_{1x}\\[5pt]\upsilon_{2x}=-\lambda _{1x}\;|\varrho_{2x}|^{\frac{1}{2}}\; \text{sign}(\varrho_{2x})-\lambda _{2x}\;\int \text{sign}(\varrho_{2x}(z))dz\end{cases}\end{equation*}

(82) \begin{equation}U_x=\frac{m}{U_1}(\dot{\upsilon}_{1x}+\upsilon_{2x}).\end{equation}

Control of (y)

\begin{equation*}\begin{cases}\varrho_{1y}=y-y_d\\[5pt]\varrho_{2y}=\dot{y}-v_{1y}\\[5pt]\upsilon_{1y}=\dot{y}_d-\kappa_y\;\varrho_{1y}\\[5pt]\upsilon_{2y}=-\lambda _{1y}\;|\varrho_{2y}|^{\frac{1}{2}}\; \text{sign} (\varrho_{2y})-\lambda _{2y}\;\int \text{sign}(\varrho_{2y}(z))dz\end{cases}\end{equation*}

(83) \begin{equation}U_y=\frac{m}{U_1}(\dot{\upsilon}_{1y}+\upsilon_{2y}).\end{equation}

4.1 Scenario 1

In this scenario, we make a square trajectory where the foldable drone changes its configuration during the flight. Simulation results are presented in Figs 5 to 11. Figure 5 represents the absolute position of the quadrotor during its flight. The translation and attitude responses are displayed in Figs 6 and 7, while the control signals are illustrated in Fig. 10.

Figure 5. Absolute position of the quadrotor under disturbances.

Figure 6. Trajectory of positions.

Figure 7. Euler angles.

Figure 8. Quadrotor tracking errors.

Figure 9. Quadrotor attitude errors.

According to Figs 6 and 7, the quadrotor follows the reference trajectories. It is observed in Figs 8 and 9 that the tracking error converges to zero due to the effectiveness of the proposed controller. Figure 10 shows the evolution of the control signals. We notice that $U_2$ and $U_3$ contain chattering. This is due to the discontinuous function sign in control laws, this switching in control excites the unmodeled dynamics, which leads to oscillations of the state vector. From Fig. 11, the servomotors’ output follows the desired angles with some error in the descent and ascent, which is generated by the rotation of their arms.

Figure 10. Evolution of control signals.

Figure 11. Servomotors outputs.

The efficiency of the proposed control method is illustrated in Figs 5 to 10 via the comparison between the results obtained using the proposed controller and the classical backstepping presented in [Reference Derrouaoui, Bouzid, Guiatni, Kada, Dib and Moudjari4] in the presence of external disturbances. In comparison to the classical backstepping controller provided, the simulation results clearly indicate better performance, quicker convergence, high precision tracking, and resilience.

Figure 12. Absolute position of the quadrotor under disturbances.

Figure 13. Trajectory of positions.

Figure 14. Euler angles.

In the light of these results, we have demonstrated the efficiency of this control approach to follow the desired trajectories and simulate them in the presence of disturbances.

Figure 15. Quadrotor tracking errors.

Figure 16. Quadrotor attitude errors.

Figure 17. Evolution of control signals.

4.2 Scenario 2

In this scenario, we make a circular trajectory, which is more complex. The obtained results are displayed in Figs 12 to 17. Figure 12 shows the absolute position of the quadrotor during the flight in the circular trajectory. Figures 13 and 14 represent the evolution of the position and orientation of the vehicle. The actuator force controls are displayed in Fig. 17.

On the basis of these results, we can clearly see from Figs 12, 13 and 14 a good tracking of the desired trajectories. We can also see that the outputs of the system converge towards the setpoint trajectories quickly. We notice the presence of chattering in the evolution of the control signals due to the switching in control laws.

The comparison between the results obtained using the proposed controller and the classical backstepping is illustrated in Figs 12 to 17 in the presence of external disturbances. When compared to the classical backstepping controller presented in Ref.(Reference Derrouaoui, Bouzid, Guiatni, Kada, Dib and Moudjari4), the simulation results clearly demonstrate satisfactory performance, faster convergence, high accuracy tracking and robustness.

Using the simulation results, we have demonstrated the efficiency of this control approach to follow the desired trajectories and to adapt to the change in the flight configuration. Furthermore, simulations have shown the robustness of the proposed control in the presence of disturbances. We interpret this by the introduction of the correction term based on the super-twisting algorithm.

To better clarify the comparison and strengthen our results, we have opted for a quantitative comparison with respect to the trajectory error and the control inputs.

The Mean Squared Error MSE and the Mean Squared Input MSI are calculated for N points as follows:

(84) \begin{equation}\begin{cases}MSE_{\text{position}}=\dfrac{1}{N}\left(\varrho_x^T\varrho_x+\varrho_y^T\varrho_y+\varrho_z^T\varrho_z\right)\\[9pt]MSE_{\text{attitude}}=\dfrac{1}{N}\left(\varrho_{\phi}^T\varrho_{\phi}+\varrho_{\theta}^T\varrho_{\theta}+\varrho_{\psi}^T\varrho_{\psi}\right)\\[9pt]MSI=\dfrac{1}{N}\left(U_1^TU_1+U_2^TU_2+U_3^TU_3+U_4^TU_4\right)\end{cases}.\end{equation}

In addition to the results found previously, the comparison between the two control laws in terms of error and input shows the efficiency of our proposed controller, in which the calculated values of the controller based on the super-twisting algorithm are lower than the classical backstepping controller, and this is proven in both scenarios (see Tables 5 and 6). The two values, MSE and MSI are considered, respectively, as the tracking error and the energy consumed during the flight.

Table 5. Quantitative comparison for the $1^{st}$ scenario

Table 6. Quantitative comparison for the $2^{nd}$ scenario