2.1 General description of the transformable quadrotor

The designed quadrotor has four rotating arms as shown in Figs 1, 2 and 3. To achieve a reconfigurable geometry, these arms must be rotated independently, respecting the angle of rotation ${\psi}_i(t)|_{i=1,...,4}$ . The angle between each propeller in its horizontal state (parallel to $y_m$ -axis) and its corresponding arm (in the vertical state) is $28.70^\circ$ . In addition, for improved security and to avoid collisions of the propellers, we slightly shifted two rotors vertically so that the propellers never overlap. We emphasize that the rotation mechanism is simple, using a low-cost frame and having lower weight.

Figure 1. Design of the transformable quadrotor that can change its configuration while flying. Each rotor is connected to a rotating arm. Each rotating arm is actuated by a servomotor and can rotate independently. Each servomotor is linked to a central body. (1) Propeller. (2) Rotating arm. (3) Servo-arm junction. (4) Servomotor.

Figure 2. Schematic of the designed transformable quadrotor.

Figure 3. Designed configurations.

The position of the CoG moves when the quadrotor configuration changes. Consequently, the inertia matrix and roll and pitch moments in this case vary depend on the orientations of the arms. They are written respectively as $\mathfrak{I}_{3 \times 3}(\psi_i (t))$ , $\tau_{\varphi}(\psi_i(t))$ and $\tau_{\theta}(\psi_i(t))$ . The total thrust ${\mathcal{T}}$ and the moment around the $z_m$ -axis do not depend on the configuration, and their expressions are the same as those for the classical quadrotor. The central body and servomotors do not rotate around the $z_m$ -axis, whilst the arms and rotors rotate around the same axis. The main configuration leads to many other possible configurations by rotating the arms. Moreover, by changing the speed of the four rotors ${\Omega}_i|_{i=1,...,4}$ or the position of the arms, the quadrotor can produce different motions (roll, pitch, yaw and translation).

The designed quadrotor is considered to be a body assembly composed of a central body of mass $\mathfrak{m}_0$ , four servomotors attached to the central body of mass $\mathfrak{m}_{1,i}$ , four rotating arms attached to servomotors of mass $\mathfrak{m}_{2,i}$ and four rotors of mass $\mathfrak{m}_{3,i}$ . The four rotors are numbered 1, 2, 3 and 4, and each arm is numbered according to the rotor fixed on it. We stress that the axes of rotation of the rotors and rotating arms are parallel to the $z_m$ -axis (Fig. 2).

By actuating the rotation servomotors, the designed quadrotor can transform during flight from the well-known $``\texttt{X}"$ configuration to other particular configurations depending on the flight environment and designated tasks. Each configuration is defined by the angles of rotation of the quadrotor arms ${\psi_i(t)}$ . Figure 3 and Table 1 show the six morphologies that we have used. The choice of these latter is based mainly on their advantages and complexity.

Table 1 Characteristics of the configurations

2.2 Inertia and centre of gravity analysis

Since the majority of the considered morphologies are asymmetrical, and to validate the calculations for their use in the generic model and simulation (Sections 3 and 5), the inertia and the global CoG position are determined by two different methods:

2.2.1 Computer-Assisted Design (CAD)-based calculation

Three-Dimensional (3D) CAD modeling of the UAV using SolidWorks software (Figs 1 and 3) allows the direct calculation of the inertias $\mathfrak{I}_{xx}$ , $\mathfrak{I}_{yy}$ , $\mathfrak{I}_{zz}$ , $\mathfrak{I}_{xy}=\mathfrak{I}_{yx}$ and CoG of each configuration separately. This will be exploited in this work as a tool to confirm, validate and compare our results with those obtained by the mathematical-based development.

The variation of the CoG coordinates for each configuration calculated in the SolidWorks design environment is presented in Table 2.

Table 2 CoG coordinates of the considered configurations

We also get the variation of the inertia, as illustrated in Table 3.

Table 3 Inertias of the considered configurations

2.3 Center of gravity

Unlike a standard quadrotor, the morphology of the proposed UAV is asymmetrical in most cases. Consequently, the CoG varies and must be instantly recalculated and injected into the dynamic model when the configuration is modified. The dynamic formula that estimates the change in the CoG can be represented as follows:

(1) \begin{equation}\overrightarrow{{O}\mathbb{G}} = \frac{{\sum\limits_{i = 1}^4 {\sum\limits_{j = 1}^3 {{\mathfrak{m}_{j,i}}} {{\overrightarrow {{O}\mathbb{G}} }_{j,i}}} }}{{{\mathfrak{m}_0} + \sum\limits_{i = 1}^4 {\sum\limits_{j = 1}^3 {{\mathfrak{m}_{j,i}}} } }} \end{equation}

with

(2) \begin{equation}\ \overrightarrow {O\mathbb{G}} = \left[ {\begin{array}{*{20}{c}} {{x_G}}\\[2pt] {{y_G}}\\[2pt] {{z_G}} \end{array}} \right],{\overrightarrow {O\mathbb{G}} _0} = \left[ {\begin{array}{*{20}{c}} {{x_0}}\\[2pt] {{y_0}}\\[2pt] {{z_0}} \end{array}} \right]{\overrightarrow {,O\mathbb{G}} _{1,i}} = \left[ {\begin{array}{*{20}{c}} {{x_{1,i}}}\\[2pt] {{y_{1,i}}}\\[2pt] {{z_{1,i}}} \end{array}} \right]\end{equation}

(3) \begin{equation}{\overrightarrow {O\mathbb{G}} _{2,i}} = \left[ {\begin{array}{*{20}{c}} {{x_{2,i}}}\\[2pt] {{y_{2,i}}}\\[2pt] {{z_{2,i}}} \end{array}} \right],{\overrightarrow {O\mathbb{G}} _{3,i}} = \left[ {\begin{array}{*{20}{c}} {{x_{3,i}}}\\[2pt] {{y_{3,i}}}\\[2pt] {{z_{3,i}}} \end{array}} \right]\end{equation}

$\ {O}$ : the geometric centre (the origin of the body frame).

${\overrightarrow {O\mathbb{G}}} \in \mathbb{R}^{3 \times 1}$ : the offset between the geometric centre of the system and the global CoG.

${\overrightarrow {O\mathbb{G}} _{0}} \in \mathbb{R}^{3 \times 1}$ : the CoG coordinates of the central body.

${\overrightarrow {O\mathbb{G}} _{1,i}} \in \mathbb{R}^{3 \times 1}$ : the CoG coordinates of the servomotors.

${\overrightarrow {O\mathbb{G}} _{2,i}} \in \mathbb{R}^{3 \times 1}$ : the CoG coordinates of the rotating arms.

${\overrightarrow {O\mathbb{G}} _{3,i}} \in \mathbb{R}^{3 \times 1}$ : the CoG coordinates of the rotors.

The evolution of the CoG position for each flight configuration, obtained using approximate calculation functions in MATLAB software, is shown in Fig. 4. In fact, the calculation formula given by Equation (1) takes into account all the intermediate cases; i.e. the transition from one configuration to another is included in the calculation.

Figure 4. Evolution of the CoG for the different configurations.

2.4 Inertia matrix

The calculation of the inertia for a compound system and elements with a non-regular form is not obvious. Therefore, some approximations should be made to simplify the mathematical development.

The equations that model the variation of the inertia matrix will be developed for each part and finally for the complete system. Using the Huygens–Steiner theorem, the variable inertia matrices are calculated as

(4) \begin{align}\ {I_{(0)/O}} = {m_0}diag{\left( {\frac{{{a^2} + h_0^2}}{{12}},\frac{{{a^2} + h_0^2}}{{12}},\frac{{2{a^2}}}{{12}}} \right)_{{\mathbb{G}_0}}}\\[-20pt]\nonumber\end{align}

where $I_{(0)/O}$ is the inertia matrix of the central body, a represents its length and width together and $ h_{0} $ its height.

The servomotors are assumed to be rectangular cuboids with length of $ L_{1} $ , width of $ w_{1}$ and height of $ h_{1} $ . Their inertia matrices are given as follows:

(5) \begin{equation}\ {I_{(1,i)/O}} = {I_{{\mathbb{G}_{1,i}}}} + {m_{(1,i)}}\left[ {\begin{array}{c@{\quad}c@{\quad}c} {y_{1,i}^2}&{ - {x_{1,i}}{y_{1,i}}}&0\\ { - {x_{1,i}}{y_{1,i}}}&{x_{1,i}^2}&0\\ 0&0&{x_{1,i}^2 + y_{1,i}^2} \end{array}} \right]\end{equation}

where

(6) \begin{align}\ {I_{{\mathbb{G}_{1,i}}}} = {m_{(1,i)}}diag{\left( {\frac{{L_1^2 + h_1^2}}{{12}},\frac{{w_1^2 + h_1^2}}{{12}},\frac{{L_1^2 + w_1^2}}{{12}}} \right)_{{\mathbb{G}_{1,i}}}}\\[-20pt]\nonumber\end{align}

The rotating arms are also supposed to be rectangular cuboid with length of $ L_2 $ , width of $ w_{2}$ and height of $ h_{2}$ . Their inertia matrices are determined by the following formula:

(7) \begin{align}\ {I_{(2,i)/{O_{(Rot)}}}} = {R_z}({\psi _i(t)}){I_{(2,i)/O}}{R_z}{({\psi _i(t)})^T}\\[-20pt]\nonumber\end{align}

where

(8) \begin{equation}\ {I_{(2,i)/O}} = {I_{{\mathbb{G}_{2,i}}}} + {m_{(2,i)}}\left[ {\begin{array}{c@{\quad}c@{\quad}c} {y_{2,i}^2}&{ - {x_{2,i}}{y_{2,i}}}&0\\ { - {x_{2,i}}{y_{2,i}}}&{x_{2,i}^2}&0\\ 0&0&{x_{2,i}^2 + y_{2,i}^2} \end{array}} \right]\end{equation}

(9) \begin{equation}\ {I_{{\mathbb{G}_{2,i}}}} = {m_{(2,i)}}diag{\left( {\frac{{{{L_2}^2} + h_2^2}}{{12}},\frac{{w_2^2 + h_2^2}}{{12}},\frac{{{{L_2}^2} + w_2^2}}{{12}}} \right)_{{\mathbb{G}_{2,i}}}}\end{equation}

and ${R_z}(\psi _i(t)) \in \mathbb{R}^{3 \times 3}$ is the rotation matrix on the ${z_b}$ axis.

The rotors are assumed to be cylinders with radius of $ r_3 $ and height of $ h_{3} $ . Their inertia matrices are expressed as

(10) \begin{align}\ {I_{(3,i)/{O_{(Rot)}}}} = {R_z}({\psi _i(t)}){I_{(3,i)/O}}{R_z}{({\psi _i(t)})^T}\\[-20pt]\nonumber\end{align}

where

(11) \begin{equation}\ {I_{(3,i)/O}} = {I_{{\mathbb{G}_{3,i}}}} + {m_{(3,i)}}\left[ {\begin{array}{*{20}{c}} {y_{3,i}^2}&{ - {x_{3,i}}{y_{3,i}}}&0\\ { - {x_{3,i}}{y_{3,i}}}&{x_{3,i}^2}&0\\ 0&0&{x_{3,i}^2 + y_{3,i}^2} \end{array}} \right]\end{equation}

(12) \begin{equation}\ {I_{{\mathbb{G}_{3,i}}}} = {m_{(3,i)}}diag{\left( {\frac{{{{r_3}^2}}}{4} + \frac{{h_3^2}}{{12}},\frac{{{{r_3}^2}}}{4} + \frac{{h_3^2}}{{12}},\frac{{{{r_3}^2}}}{2}} \right)_{{\mathbb{G}_{3,i}}}}\end{equation}

Finally, we obtain the global inertia matrix of the system, which depends on the angle of rotation of each arm $\psi_i(t)$ , as

(13) \begin{align}\ \begin{array}{l}{\mathfrak{I}_{3 \times 3}(\psi_i (t))} = {I_{(0)/O}} + \sum\limits_{i = 1}^4 {{I_{(1,i)/O}}} + \sum\limits_{i = 1}^4 {{I_{(2,i)/{O_{(Rot)}}}}}+ \sum\limits_{i = 1}^4 {{I_{(3,i)/{O_{(Rot)}}}}}\end{array}\\[-20pt]\nonumber\end{align}

The evolution of the inertia for the six aerial configurations is shown in Fig. 5.

Figure 5. Evolution of inertia along the $x_{m}$ , $y_{m}$ and $z_{m}$ axes as well as the inertia of the non-diagonal terms.

2.5 Results and discussion

On the basis of the results obtained by SolidWorks CAD (Table 3) and the inertias calculated by the approximate formulas (Fig. 5), the errors between the two methods according to the different axes are summarised in Table 4:

Table 4 Inertia errors

The small difference between the values given by the two methods as shown in Table 4 is justified by the assumptions imposed on the geometry of the quadrotor components.

According to Figs 3 and 4 and Table 2, we can see that, for the configurations $``\texttt{X}"$ and $``\texttt{H}"$ , the estimated CoG of the system is confused with the CoG given by SolidWorks, where its values x(m) and y(m) are equal to zero. We can explain this by the symmetry of the configurations along the two axes $x_{m}$ and $y_{m}$ . In fact, the $``\texttt{O}"$ configuration is not geometrically symmetrical about the $x_{m}$ and $y_{m}$ axes, but it is symmetrical in terms of the mass distribution in relation to these axes. For the $``\texttt{Y}"$ , $``\texttt{YI}"$ and $``\texttt{T}"$ configurations, the ordinates y(m) of the CoG are equal to zero. This is justified by the symmetry of these configurations with respect to the $x_{m}$ -axis. However, the abscissas x(m) are different from zero due to the asymmetry with respect to the $y_{m}$ -axis, where the error in this case is on the order of $10^{-3}$ (m).

3.1 Flight dynamics

The flight dynamics of the aerial vehicle is described in two coordinate frames, as shown previously in Fig. 2. The inertial frame $\mathfrak{R}_{i}(O_{i},x_i,y_i,z_i)$ is assumed to be fixed while $\mathfrak{R}_{m}(O,x_{m},y_{m},z_{m})$ is considered to be a mobile frame. Unlike most conventional quadrotors, the CoG does not coincide with the Geometric Center (GC).

The rotation from the mobile reference system $\mathfrak{R}_{m}$ to the inertial fixed frame $\mathfrak{R}_{i}$ is achieved by three successive rotations around $X-Y-Z$ axes, meaning that the attitude is obtained first by the roll angle $\varphi$ , then by the pitch angle $\theta$ and then by the yaw angle $\psi$ .

The direction cosine matrix, denoted R^{(Reference Spong, Hutchinson and Vidyasagar63)}, from the mobile frame $\mathfrak{R}_{m}$ to the inertial fixed frame $\mathfrak{R}_{i}$ can be expressed as

(14) \begin{equation}\textrm{R}\textrm{=}\left[\begin{array}{c@{\quad}c@{\quad}c}{\textrm{c}}_{\psi }{\textrm{c}}_{\theta } & {\textrm{c}}_{\psi }{\textrm{s}}_{\theta }{\textrm{s}}_{\varphi }\textrm{-}{\textrm{s}}_{\psi }{\textrm{c}}_{\varphi } & {\textrm{c}}_{\psi }{\textrm{s}}_{\theta }{\textrm{c}}_{\varphi }\textrm{+}{\textrm{s}}_{\psi }{\textrm{s}}_{\varphi } \\[6pt] {\textrm{s}}_{\psi }{\textrm{c}}_{\theta } & {\textrm{s}}_{\psi }{\textrm{s}}_{\theta }{\textrm{s}}_{\varphi }\textrm{+}{\textrm{c}}_{\psi }{\textrm{c}}_{\varphi } & {\textrm{s}}_{\psi }{\textrm{s}}_{\theta }{\textrm{c}}_{\varphi }\textrm{-}{\textrm{c}}_{\psi }{\textrm{s}}_{\varphi } \\[6pt] \textrm{-}{\textrm{s}}_{\theta } & {\textrm{c}}_{\theta }{\textrm{s}}_{\varphi } & {\textrm{c}}_{\theta }{\textrm{c}}_{\varphi }\\[3pt] \end{array}\right]\end{equation}

where $ \textrm{R}\in SO(3)=\{\textrm{R}\in \textrm{R}^{3 \times3} | \textrm{R}^T\textrm{R}=I_{3 \times3},\ det (\textrm{R})=1\} $ and $s_{{(\cdot)}}$ and $c_{(\cdot)}$ are abbreviations for ${\textrm{sin} \left(\cdot\right)\ }$ and ${\textrm{cos} \left(\cdot\right)}$ , respectively.

The linear and angular velocity vectors of the body in the mobile frame are represented respectively as ${\Lambda}^m=(u,v,w)^T \in \mathbb{R}^3$ and $\varsigma=(p,q,r)\in \mathbb{R}^{3}$ .

The transformable quadrotor has a mass $\mathfrak{m}$ and a variable inertia $\mathfrak{I}(\psi_i (t)) \in \mathbb{R}^{3 \times 3}$ , which is calculated with respect to the mobile frame $\mathfrak{R}_{m}(o_{m},x_{m},y_{m},z_{m})$ and depends on the angular positions of the four arms $\psi_i (t)$ , where $i=1,2,3,4.$

Let $\Upsilon= (\varphi,\theta,\psi)^{T} \in \mathbb{R}^3$ describes the orientation of the mobile and $\xi$ = $(x,y,z)^{T}\in \mathbb{R}^3$ denotes its position with respect to $\mathfrak{R}_{i}$ .

The relation between the velocities and the external forces $\textit{f}^m=(\textit{f}_x^m, \textit{f}_y^m, \textit{f}_z^m)^T \in \mathbb{R}^{3}$ and moments ${\tau}^m=({\tau}_x^m, {\tau}_y^m, {\tau}_z^m)^T \in \mathbb{R}^{3}$ , applied to the CoG, can be written using the Newton–Euler formalism as

(15) \begin{equation}\begin{bmatrix}\mathfrak{m}\mathfrak{I}_{3 \times 3}(\psi_i (t)) & O_{3 \times 3}\\[3pt] O_{3 \times 3} & \mathfrak{I}(\psi_i (t))\end{bmatrix}\begin{bmatrix}\dot \varLambda^m \\[3pt] \dot \varsigma\end{bmatrix}+ \begin{bmatrix}\varsigma \times \mathfrak{m}\varLambda^m \\[3pt] \varsigma \times \mathfrak{I}(\psi_i (t)) \varsigma\end{bmatrix}= \begin{bmatrix}\textit{f}^m \\[3pt] {\tau}^m\end{bmatrix}\end{equation}

The symbol $O_{3 \times 3}$ means a $3 \times 3$ -dimensional zero matrix, and $\times$ denotes the cross product.

3.2 Forces and moments applied to the quadrotor

3.2.1 Forces

The quadrotor is affected by several forces, such as gravity, thrust, hub and drag.

Remark 1: All the parameters used in this subsection are mentioned in Table 5.

Table 5 Modeling parameters

The gravity vector $\mathcal{G}$ = $(0,0,-\mathfrak{m}\mathfrak{g})^{T}\in \mathbb{R}^3$ is applied along the $z_i$ -axis, where $\mathfrak{g}$ represents the gravity.

The rotation of the four rotors can easily generate a thrust force along the $z_{ri}$ -axis and a hub force in the hub plane ( $x_{ri},y_{ri}$ ). Thus, the total thrust force $\mathcal{T}$ is the sum of the forces produced by the rotation of each propeller, expressed in the body frame. In addition, the hub force is defined by two components $\mathcal{U}_{xi}$ and $\mathcal{U}_{yi}$ along two axes $x_{ri}$ and $y_{ri}$ in the plane corresponding to each rotor $\mathfrak{R}_{ri}(o_{ri},x_{ri},y_{ri},z_{ri})$ ^{(Reference Bouzid, Siguerdidjane and Bestaoui64)}. The $\mathcal{T}$ , $\mathcal{U}_{xi}$ and $\mathcal{U}_{yi}$ forces are expressed as

(16) \begin{align}\mathcal{T}=\sum_{i=1}^{4}\textit{f}_i,\ \mathcal{U}_x=-\sum_{i=1}^{4}\mathcal{U}_{xi},\ \mathcal{U}_y=-\sum_{i=1}^{4}\mathcal{U}_{yi}\end{align}

where $\textit{f}_{i}=C_{\textit{f}}{\rho}A({\Omega}_iR_r)^2$ and $\mathcal{U}_{i}=C_{\mathcal{U}}{\rho}A({\Omega}_iR_r)^2$ .

The quadrotor can also produce a drag force (friction force) $\mathcal{D}\in \mathbb{R}^3$ expressed in $\mathfrak{R}_{i}(o_{i},x_i,y_i,z_i)$ . This force is due to the movement of the quadrotor body affected by the wind, and can be defined as

(17) \begin{align}\nonumber\ {\mathcal{D}} &=({\mathcal{D}x}\ {\mathcal{D}y}\ \,{\mathcal{D}z})^{T} \\ &=\textrm{diag}\left( { - {\mathcal{K}_{\mathcal{D}x}}}\ { - {\mathcal{K}_{\mathcal{D}y}}}\ { - {\mathcal{K}_{\mathcal{D}z}}} \right)\mathop \xi \limits^. \end{align}

Once the modeling of the different forces is established, the global force vector in $\mathfrak{R}_{m}(O,x_{m},y_{m},z_{m})$ can be expressed as

(18) \begin{align}\nonumber\textit{f}^m &= \begin{bmatrix}\mathcal{U}_x\\[3pt]\mathcal{U}_y \\[3pt]\mathcal{T}\end{bmatrix}+\textrm{R}^T\begin{bmatrix}{\mathcal{D}x}\\[3pt]{\mathcal{D}y} \\[3pt]{\mathcal{D}z}\end{bmatrix}+\textrm{R}^T\begin{bmatrix}0\\[3pt]0 \\[3pt]\mathcal{G}\end{bmatrix} \\[3pt]&= \begin{bmatrix}-\sum_{i=1}^{4}\mathcal{U}_{xi}\\[3pt]-\sum_{i=1}^{4}\mathcal{U}_{yi} \\[3pt]\sum_{i=1}^{4}\textit{f}_i\end{bmatrix}+\textrm{R}^T\begin{bmatrix}{\mathcal{-K}_{\mathcal{D}x}}{\dot{x}}\\[3pt]{\mathcal{-K}_{\mathcal{D}y}}{\dot{y}} \\[3pt]{\mathcal{-K}_{\mathcal{D}z}}{\dot{z}}\end{bmatrix}+\textrm{R}^T \begin{bmatrix}0\\[3pt]0 \\[3pt]-\mathfrak{m}\mathfrak{g}\end{bmatrix}\end{align}

where $\textrm{R}^{-1}=\textrm{R}^T$ .

3.2.2 Moments

Herein, we seek to present the different moments around the three axes $x_m$ , $y_m$ and $z_m$ . The quadrotor can turn about these axes by applying the moments ${\tau}_{\varphi}$ , ${\tau}_{\theta}$ and ${\tau}_{\psi}$ , respectively. The moments of roll and pitch ( ${\tau}_{\varphi}$ , ${\tau}_{\theta}$ ) are induced by the thrust forces $f_i$ generated by each rotor. In our case, they depend on the folding angle of the arms as

(19) \begin{align}\nonumber\ {{\tau}_{(\varphi,\theta)}} &= \sum\limits_{i = 1}^4 {{{\overrightarrow {\mathbb{GG}} }_{(3,i)}} \times \overrightarrow{\textit{f}_i}}\\ &=\begin{bmatrix}\sum_{i=1}^{4}(y_{3,i}(t)-y_G(t)){\textit{f}_i}\\ \\[-7pt] -\sum_{i=1}^{4}(x_{3,i}(t)-x_G(t)){\textit{f}_i}\end{bmatrix}\end{align}

The moving coordinates of the centre of gravity of the rotors are given as

(20) \begin{align}\ \left\{ {\begin{array}{l} {{x_{3,1}(t)} = \frac{a}{2} + \left( {{L_2}+{r_3}} \right)\sin {\psi _1(t)}}\\[5pt] {{x_{3,2}(t)} = \frac{a}{2} + \left( {{L_2}+{r_3}} \right)\cos {\psi _2(t)}}\\[5pt] {{x_{3,3}(t)} = - \frac{a}{2} - \left( {{L_2}+{r_3}} \right)\sin {\psi _3(t)}}\\[5pt] {{x_{3,4}(t)} = - \frac{a}{2} - \left( {{L_2}+{r_3}} \right)\cos {\psi _4(t)}}\\[5pt] {{y_{3,1}(t)} = - \frac{a}{2} - \left( {{L_2}+{r_3}} \right)\cos {\psi _1(t)}}\\[5pt] {{y_{3,2}(t)} = \frac{a}{2} + \left( {{L_2}+{r_3}} \right)\sin {\psi _2(t)}}\\[5pt] {{y_{3,3}(t)} = \frac{a}{2} + \left( {{L_2}+{r_3}} \right)\cos {\psi _3(t)}}\\[5pt] {{y_{3,4}(t)} = - \frac{a}{2} - \left( {{L_2}+{r_3}} \right)\sin {\psi _4(t)}} \end{array}} \right.\end{align}

In addition, the total yaw moment $\tau_{\psi}$ along $z_m$ -axis is considered to be

(21) \begin{align}\ \tau_{\psi}=\sum_{i=1}^{4}(-1)^{i+1}\mathcal{Q}_{i}\\[-20pt]\nonumber\end{align}

where $\mathcal{Q}_{i}=C_{\mathcal{Q}}{\rho}AR_r({\Omega}_iR_r)^2$ .

Wind effects can generate an aerodynamic friction torque ${\mathcal{M}_{\mathcal{A}}}\in \mathbb{R}^3$ , expressed as

(22) \begin{align}\nonumber\ {\mathcal{M}_{\mathcal{A}}} &={\begin{array}{*{20}{c}} ({\mathcal{M}_{\mathcal{A}x}}&{\mathcal{M}_{\mathcal{A}y}}&{\mathcal{M}_{\mathcal{A}z}})^{T} \end{array}}\\ &=\textrm{diag}\left({\begin{array}{*{20}{c}} {\mathcal{K}_{\mathcal{A}x}}&{\mathcal{K}_{\mathcal{A}y}}&{\mathcal{K}_{\mathcal{A}z}} \end{array}} \right){\dot{\varsigma} ^2}\end{align}

The flapping moment $\mathcal{L}_i$ is created due to the difference in lift between the blade that moves back and the blade that advances. It is given by its two components as

(23) \begin{align}\mathcal{L}_x=\sum_{i=1}^{4}(-1)^{i+1}\mathcal{L}_{xi},\ \mathcal{L}_y&=\sum_{i=1}^{4}(-1)^{i+1}\mathcal{L}_{yi}\end{align}

where $\mathcal{L}_{i}=C_{\mathcal{L}}AR_r({\Omega}_iR_r)^2$ .

In the case when the quadrotor rotates around its $x_m$ -axis, a gyroscopic moment $\mathcal{O}_y$ occurs on each rotating rotor along the $y_m$ -axis. However, when the quadrotor rotates around its $y_m-$ axis, the rotating rotors have also a gyroscopic moment $\mathcal{O}_x$ along the $x_m$ -axis. The expression of the global gyroscopic moment is

(24) \begin{align}\mathcal{O}_x = \mathfrak{J}_rq\sum_{i=1}^{4}(-1)^{i+1}\Omega_i,\ \mathcal{O}_y = -\mathfrak{J}_rp\sum_{i=1}^{4}(-1)^{i+1}\Omega_i\end{align}

where $\mathfrak{J}_r$ is the inertia of the rotor.

Another moment is produced from the hub forces along the $x_m$ -axis and $y_m$ -axis, called the hub moment. It is given as

(25) \begin{align}\nonumber\mathcal{H} &= \sum_{i=1}^{4}(x_{3,i}(t)-x_G(t))({\textrm{s}}_{\psi_i (t) }\mathcal{U}_{xi}+{\textrm{c}}_{\psi_i (t) }\mathcal{U}_{yi})\\ &-\sum_{i=1}^{4}(y_{3,i}(t)-y_G(t))({\textrm{c}}_{\psi_i (t) }\mathcal{U}_{xi}-{\textrm{s}}_{\psi_i (t) }\mathcal{U}_{yi})\end{align}

Finally, the global moment vector $\mathcal{\tau}^m$ , which includes various moments, is expressed as

(26) \begin{align} \nonumber \mathcal{\tau}^m &= \begin{bmatrix} \mathcal{\tau}_\varphi+\mathcal{O}_x+\mathcal{L}_x+{\mathcal{M}_{\mathcal{A}x}}\\[2pt] \mathcal{\tau}_\theta +\mathcal{O}_y+\mathcal{L}_y+{\mathcal{M}_{\mathcal{A}y}}\\[2pt] \mathcal{\tau}_\psi+\mathcal{H}+{\mathcal{M}_{\mathcal{A}z}} \end{bmatrix}\\[2pt] &= \begin{bmatrix} \sum_{i = 1}^4[(y_{3,i}(t)-y_G(t))\textit{f}_i+\mathfrak{J}_rq{(-1)}^{i+1}\Omega_i+{(-1)}^{i+1}\mathcal{L}_{xi}]+{\mathcal{K}_{\mathcal{A}x}}p^2\\ \\ \sum_{i = 1}^4[(x_G(t)-x_{3,i}(t))\textit{f}_i-\mathfrak{J}_rp{(-1)}^{i+1}\Omega_i+{(-1)}^{i+1}\mathcal{L}_{yi}]+{\mathcal{K}_{\mathcal{A}y}}q^2\\ \\ \sum_{i = 1}^4[{(-1)}^{i+1}Q_i +(x_{3,i}(t)-x_G(t))({\textrm{s}}_{\psi_i (t) }\mathcal{U}_{xi}+{\textrm{c}}_{\psi_i (t) }\mathcal{U}_{yi})\\ \quad -(y_{3,i}(t)-y_G(t))({\textrm{c}}_{\psi_i (t) }\mathcal{U}_{xi}-{\textrm{s}}_{\psi_i (t) }\mathcal{U}_{yi})] +{\mathcal{K}_{\mathcal{A}z}}r^2 \end{bmatrix} \end{align}

Equation (26) clearly shows the major difference in the dynamic model with respect to the conventional quadrotor.

3.3 Servomotor modeling

The proposed dynamic model for the servomotors is based on the physical principles for modern Dynamixel servomotors^{(Reference Maximo, Ribeiro and Afonso65)} and includes several phenomena such as viscous friction and saturation. These servomotors are used to rotate the arms of the proposed quadrotor. They are controlled by an internal Proportional–Integral–Derivative (PID) regulator. After some mathematical operations and using the Laplace transformation, the mechanical open-loop transfer function as well as the electrical part are given respectively by the following relations:

(27) \begin{equation}\dfrac{\mathcal{\theta}_l(s)}{\mathcal{I}(s)}=\dfrac{N\textit{n}K_t}{(J_l+J_mN^2\textit{n})s^2+(b_mN^2\textit{n})s}\end{equation}

(28) \begin{equation}{\mathcal{I}(s)}=\dfrac{U(s)}{R_\mathcal{I}+L_\mathcal{I}s}-\dfrac{sN\theta_l(s)}{k_w(R_\mathcal{I}+L_\mathcal{I}s)}\end{equation}

The closed-loop transfer function including the PID regulator is given as follows:

(29) \begin{equation}{{\mathcal{F}(s)}=\dfrac{A_{\mathcal{F}}(k_Ds^2+k_Ps+k_I)}{B_{\mathcal{F}}s^3+C_{\mathcal{F}}s^2+D_{\mathcal{F}}s+E_{\mathcal{F}}}}\end{equation}

The servomotor parameters are described in Table 6, where $A_{\mathcal{F}}$ , $B_{\mathcal{F}}$ , $C_{\mathcal{F}}$ , $D_{\mathcal{F}}$ , $E_{\mathcal{F}}$ are constants specific to the servomotor.

Table 6 Dynamixel actuator and rotor parameters^{(Reference Maximo, Ribeiro and Afonso65,Reference Derafa, Madani and Benallegue66)}

3.4 Rotor modeling

Knowledge of the rotor dynamics allows low-level rotor control. This control presents the energy required by the different rotors in terms of the current and voltage needed to turn them with a desired speed. Studies on brushless motors have shown that rotors are driven by Direct-Current (DC) motors, they can be modeled as^{(Reference Derafa, Madani and Benallegue66)}

(30) \begin{align}\ \left\{ {\begin{array}{*{20}{c}} {\mathcal{V}= R_rI +L \dfrac{dI}{dt}+E}\\[8pt] {\mathcal{T}_m =\mathfrak{J}_r \dfrac{d\Omega}{dt}+C_s+\mathcal{T}_r}\\ \end{array}} \right.\end{align}

Since this type of motors usually have a negligible inductance, the rotor dynamics can be approximated as

(31) \begin{align}{\mathcal{V}=\dfrac{\dot{\Omega}+q_2{\Omega}^2+q_1{\Omega}+q_0}{q_3}}\\[-18pt]\nonumber\end{align}

where $E=k_e\Omega$ , $\mathcal{T}_m=k_mI$ , $\mathcal{T}_r=k_r{\Omega}^2$ , ${q_0}=\dfrac{C_s}{\mathfrak{J}_r}$ , ${q_1}=\dfrac{k_ek_m}{R_r\mathfrak{J}_r}$ , ${q_2}=\dfrac{k_r}{\mathfrak{J}_r}$ and ${q_3}=\dfrac{k_m}{R_r\mathfrak{J}_r}$ .

3.5 Control allocation matrix

We assume that the quadrotor is controlled to follow paths with small speed and maneuvers. Consequently, the Euler angular velocities are considered as $\dot\varphi\approx p, \ \ \dot\theta\approx q, \ \ \dot\psi\approx r$ and the aerodynamic coefficients are taken with their simplified formulas^{(Reference Bangura and Mahony67)} as follows:

(32) \begin{align} \textit{f}_{i}=b\Omega_i^2, \mathcal{Q}_{i}=d\Omega_i^2\\[-18pt]\nonumber\end{align}

where b and d are the thrust and drag coefficients, respectively.

The transformable quadrotor has eight control inputs. The control vector of its altitude and attitude is defined as $u=(u_1, u_2, u_3, u_4)^T\in \mathbb{R}^{4}$ , while the control vector of the servomotor positions is given as $s=(u_5, u_6, u_7, u_8)^T\in \mathbb{R}^{4}$ , where

(33) \begin{align} u_1=\mathcal{T}, \ \ u_2={{\tau}_{\varphi}}, \ \ u_3={{\tau}_{\theta}}, \ \ u_4={{\tau}_{\psi}}\\[-18pt]\nonumber \end{align}

The relation between the total thrust force, the moments applied to the CoG and the propeller square velocities can be expressed in matrix notation as

(34) \begin{align} u={\varDelta}(\psi_i (t))\eta \\[-18pt]\nonumber\end{align}

with $\eta \in \mathbb{R}^{4}$ is the vector including the squared propeller velocities:

(35) \begin{align} \eta=\left[\Omega^2_1, ...,\Omega^2_{4}\right]^T \\[-18pt]\nonumber\end{align}

and $\varDelta(\psi_i (t)) \in \mathbb{R}^{4*4}$ is the control allocation matrix. The latter represents a crucial element in the architecture of the quadrotor control as it transforms the square of the rotational speeds of the propellers ${\Omega}^2_i|_{i=1,...,4}$ into a total thrust force ${\mathcal{T}}$ and moments ${\tau}_{\varphi}$ , ${\tau}_{\theta}$ , ${\tau}_{\psi}$ as

(36) \begin{equation} \ \left[ {\begin{array}{*{20}{c}} u_1\\[3pt] {u_2}\\[3pt] {u_3}\\[3pt] {u_4} \end{array}} \right] = \varDelta(\psi_i (t))\left[ {\begin{array}{*{20}{c}} {\Omega^2_1}\\[3pt] {\Omega^2_2}\\[3pt] {\Omega^2_3}\\[3pt] {\Omega^2_4} \end{array}} \right]\end{equation}

Considering the new expressions of ( ${\mathcal{T}},{\tau}_{\varphi}$ , ${\tau}_{\theta}$ , ${\tau}_{\psi}$ ) generated by the quadrotor, given by (19) and (32), the dynamic allocation matrix $\varDelta(\psi_i (t))$ becomes

(37) \begin{align} \ \varDelta(\psi_i (t)) = {\left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} {{b}}&{{b}\left [{{y_{3,1}(t)} - {y_G}(t)} \right]}&{{b}\left [{ {x_G}(t) - {x_{3,1}}(t)} \right] }&{{d}}\\[6pt] {{b}}&{{b}\left [{{y_{3,2}(t)} - {y_G}(t)} \right] }&{{b}\left [{ {x_G}(t) - {x_{3,2}(t)}} \right] }&{ - {d}}\\[6pt] {{b}}&{{b}\left [{{y_{3,3}(t)} - {y_G}(t)} \right] }&{{b}\left [{ {x_G}(t) - {x_{3,3}(t)}} \right] }&{{d}}\\[6pt] {{b}}&{{b}\left [{{y_{3,4}(t)} - {y_G}(t)} \right] }&{{b}\left [{ {x_G}(t) - {x_{3,4}(t)}} \right] }&{ - {d}}\\[6pt] \end{array}} \right]^T}\\[-18pt]\nonumber\end{align}

The elements of the matrix $\varDelta(\psi_i (t))$ depend strongly on the geometrical variations. Therefore, they must be instantly updated to be taken into account when the configuration changes.

To establish the complete dynamic model, translational velocities and accelerations are expressed in the frame $\mathfrak{R}_{i}(o_{i},x_i,y_i,z_i)$ . However, the angular velocities and accelerations are expressed in the frame $\mathfrak{R}_{m}(O,x_{m},y_{m},z_{m})$ .

Based on Equations (15), (18) and (26), the simplified control model can be arranged as (40).

The choice of the state vector is

(38) \begin{align}\ {X} = {\left[ {\begin{array}{*{20}{c}} \varphi, \dot{\varphi}, \theta, \dot{\theta}, \psi, \dot{\psi}, z, \dot{z}, x, \dot{x}, y, \dot{y} \end{array}}\right]^T}\\[-18pt]\nonumber\end{align}

such as

(39) \begin{equation}\ {X}= {\left[ {\begin{array}{*{20}{c}} x_1, x_{2}, x_{3}, x_{4}, x_{5}, x_{6}, x_{7}, x_{8}, x_{9}, x_{10}, x_{11}, x_{12} \end{array}}\right]^T}\end{equation}

(40) \begin{equation} \left\{\begin{array}{rl}\displaystyle \dot{x}_1 &= x_{2} \\[5pt] \dot{x}_2 &= \beta_{1}(t)x_{4}x_{6}+\beta_{2}(t)x_{4}{\Omega}_r+\beta_{3}(t)u_{2}+\beta_{4}(t){x_{2}}^2 \\[5pt] \dot{x}_3 &= x_{4} \\[5pt] \dot{x}_4 &= \beta_{5}(t)x_{2}x_{6}+\beta_{6}(t)x_{2}{\Omega}_r+\beta_{7}(t)u_{3}+\beta_{8}(t){x_{4}}^2 \\[5pt] \dot{x}_5 &= x_{6} \\[5pt] \dot{x}_6 &=\beta_{9}(t)x_{2}x_{4}+\beta_{10}(t)u_{4}+\beta_{11}(t){x_{6}}^2 \\[5pt] \dot{x}_7 &= x_{8} \\[5pt] \dot{x}_8 &=-g+u_{1}\dfrac{cos(x_{1})cos(x_{3})}{m}+\beta_{12}x_{8} \\[5pt] \dot{x}_9 &= x_{10} \\[5pt] \dot{x}_{10} &= u_{1}\frac{u_{x}}{m}+\beta_{13}x_{10} \\[5pt] \dot{x}_{11} &= x_{12} \\[5pt] \dot{x}_{12} &= u_{1}\frac{u_{y}}{m}+\beta_{14}x_{12} \end{array} \right. \end{equation}

where

$\beta_{1}(t)$ = $\dfrac{I_{yy}(\psi_i(t))-I_{zz}(\psi_i(t))}{I_{xx}(\psi_i(t))}$ , $\beta_{2}(t)$ = $\dfrac{-\mathfrak{J}_r}{I_{xx}(\psi_i(t))}$ $\beta_{3}(t)$ = $\dfrac{1}{I_{xx}(\psi_i(t))}$ , $\beta_{5}(t)$ = $\dfrac{I_{zz}(\psi_i(t))-I_{xx}(\psi_i(t))}{I_{yy}(\psi_i(t))}$ $\beta_{4}(t)$ = $\dfrac{-{\mathcal{K}_{\mathcal{A}x}}}{I_{xx}(\psi_i(t))}$ , $\beta_{9}(t)$ = $\dfrac{I_{xx}(\psi_i(t))-I_{yy}(\psi_i(t))}{I_{zz}(\psi_i(t))}$ $\beta_{6}(t)$ = $\dfrac{\mathfrak{J}_r}{I_{yy}(\psi_i(t))}$ , $\beta_{7}(t)$ = $\dfrac{1}{I_{yy}(\psi_i(t))}$ , $\beta_{8}(t)$ = $\dfrac{-{\mathcal{K}_{\mathcal{A}y}}}{I_{yy}(\psi_i(t))}$

$\beta_{10}(t)$ = $\dfrac{1}{I_{zz}(\psi_i(t))}$ , $\beta_{11}(t)$ = $\dfrac{-{\mathcal{K}_{\mathcal{A}z}}}{I_{zz}(\psi_i(t))}$ , $\beta_{12}$ = $\dfrac{ - {\mathcal{K}_{\mathcal{D}z}}}{\mathfrak{m}}$

$\beta_{13}$ = $\dfrac{ - {\mathcal{K}_{\mathcal{D}x}}}{\mathfrak{m}}$ , $\beta_{14}$ = $\dfrac{ - {\mathcal{K}_{\mathcal{D}y}}}{\mathfrak{m}}$ , ${\Omega}_r= \sum_{i=1}^{4}(-1)^{i+1}\Omega^2_i$

4.1 Control review

According to recent literature, a few strategies have been developed to control UAVs subject with variable geometric properties.

In Ref.^{(Reference Barbaraci62)}, the authors discussed the use of Linear Quadratic Regulation (LQR) for a quadrotor with variable geometry arms, where the angles between the arms change but the structure always remains symmetrical. The same controller is adopted in Refs^{(Reference Wallace13,Reference Avant, Lee, Katona and Morgansen26)} to stabilize a morphing quadrotor with rotating and extendable arms. Adaptive LQR was exploited by Ref.^{(Reference Falanga, Kleber, Mintchev, Floreano and Scaramuzza10)} to control a foldable quadrotor with the goal of guaranteeing stable flight. Considering the rotation of the quadrotor arms, Xiong et al.^{(Reference Xiong, Hu and Diao29)} were interested in optimizing energy consumption during the flight, where the attitude and trajectory of the prototype were controlled using two PID controllers. An X-Morf quadrotor, with two independent rotating arms, was designed in Ref.^{(Reference Desbiez, Expert, Boyron, Diperi, Viollet and Ruffier35)}, where Model Reference Adaptive Control (MRAC) is implemented to deal with the changes of the CoG and the inertia and to control the angular velocities. In addition, a classic PID controller combined with a dynamic reversal process was exploited for attitude control. Riviere et al.^{(Reference Riviere, Manecy and Viollet18)} used an adaptive PID controller to control a quad-morphing robot, taking into account the variation of the inertia matrix. Zhao et al.^{(Reference Zhao, Kawasaki, Chen, Noda, Okada and Inaba43)} developed a controller for a transformable multirotor based on the Linear Quadratic Integral (LQI) for attitude and altitude control. However, position control in the horizontal plane was achieved by a conventional PID. In addition, the developed controller is based on a simplified model and assumes that the CoG is always close to the geometric centre. The authors in Ref.^{(Reference Bucki and Mueller15)} presented an LQR attitude controller for an aerial morphing robot that can fold and unfold its arms vertically. A quadrotor with extendable arms was controlled in position and attitude using a conventional PID regulator in Ref.^{(Reference Kamil, Hazry, Wan, Razlan and Abu Bakar50)}. The same control strategy was exploited to control the position and attitude of the transformable UAV in Refs.^{(Reference Sakaguchi, Takimoto and Ushio16,Reference D’Sa42,Reference Kumar, Sridhar, Cazaurang, Cohen and Kumar53)} . Recently, Raj et al. proposed a dynamic feedback linearisation-based strategy in Ref.^{(Reference Raj, Banavar and Kothari68)} for attitude control of a transformable quadrotor UAV.

Analyzing the above-cited works, it can be concluded that the most widely used techniques are PID, LQR, LQI and feedback linearisation, or their adaptive versions. Furthermore, due to the complexity of the control strategy, they did not take into account the variation of the CoG, the inertia and the allocation matrix at the same time.

In the current work, we propose the design of a new approach for the control of the designed transformable quadrotor. This approach was chosen for its effectiveness with respect to geometric variations, as well as the limitations of the control strategies proposed in literature. In addition, it is very suitable for the hierarchical architecture of the studied system, and it is also approached as a system control tool with various unknown parameters.

4.2 Architecture and controller design

In this paper, the controller is designed to ensure the tracking of the desired trajectory ( ${x_d}$ , ${y_d}$ , ${z_d}$ ) along the three axes in addition to the yaw angle, with the different quadrotor morphologies previously given according to the mission. These reference trajectories are provided on-line by a Trajectory Generator Block (TGB), as illustrated in Fig. 6. The quadrotor is controlled by the speeds of the four motors ${\Omega}_i$ , which are deduced from the allocation matrix.

Figure 6. Control architecture of the transformable quadrotor.

The transformable quadrotor is considered as a series of fixed-configuration quadrotors, where each configuration is defined by a range of angles and each range has its own parameters (gains $k_i$ ). Consequently, the adaptation block (Fig. 6) contains the gains for each configuration. When the angles correspond to a configuration, the adaptation block assigns the corresponding parameters to the different controllers to calculate the appropriate control laws.

The CoG, the inertia and the allocation matrix vary instantaneously and have to be recomputed on-line and incorporated into the generic model of the system when the configuration is changed (Fig. 6).

The angular position of the four servomotors is controller using PID controllers ( ${u_5}$ , ${u_6}$ , ${u_7}$ , ${u_8}$ ), as displayed in Fig. 6. Their controllers are given by the following equation:

(41) \begin{equation}{u_j=k_{P}e_{\psi_i (t)}+k_{I}\int_{0}^{t}e_{\psi_i (t)} dt+k_{D}\dot{e}_{\psi_i (t)}}, j,i=5,..,8\end{equation}

where $e_{\psi_i (t)}=\psi_{id}-\psi_{i}$ is the tracking error. $k_{P}, k_{I}$ and ${k}_{D}$ denote the usual proportional, integral and derivative tuning gains. They are positive constants.

The servomotors control the position of the arms. These positions are used in the adaptation block to calculate the gains of the quadrotor controllers and by other blocks to update the current inertia matrix, the allocation matrix as well as the CoG for each configuration (Fig. 6). Then, the inertia terms will also be used by the various controller blocks.

The two controls $u_x$ and $u_y$ and the desired yaw angle ${\psi}_d$ are used to calculate the roll and pitch angles, which are required to achieve the desired position.

The different control laws will be implemented and tested in simulations in the next section.

4.3 Controller synthesis

The first subsystem is given as

(42) \begin{equation}\left\{\begin{array}{rl}\dot{x}_1 &= x_{2}\\[4pt]\dot{x}_2 &=\beta_{1}(t)x_{4}x_{6}+\beta_{2}(t)x_{4}{\Omega}_r+\beta_{3}(t)U_{2}+\beta_{4}(t){x_{2}}^2\end{array}\right.\end{equation}

For clarity and conciseness, steps 1 and 2 are described with detailed explanations, while the other steps are simplified.

The detailed design procedure is described in the two following steps.

Firstly:

The tracking error of the first subsystem $e_{1}$ and its corresponding derivative $\dot{e_{1}}$ are expressed as follows:

(43) \begin{align}e_{1}=x_1-x_{1d}, \dot{e_{1}}=x_2-\dot{x}_{1d}\\[-18pt]\nonumber\end{align}

Let $V_1(e_1)$ be a candidate positive-definite Lyapunov function,

(44) \begin{align}V_{1}(e_1)=\dfrac{1}{2} {e_1}^{2}\\[-18pt]\nonumber\end{align}

The dynamics of $\dot{V}_{1}(e_1) $ is negative definite as

(45) \begin{align}\dot{V}_{1}(e_1)=-k_1{e_1}^{2}, k_1>0.\\[-18pt]\nonumber\end{align}

Now, asymptotic Lyapunov stability is guaranteed. For this, the virtual control law of the first subsystem could be chosen as

(46) \begin{align}x_{2d}=\dot{x}_{1d}-k_1{e_1}.\\[-18pt]\nonumber\end{align}

Secondly:

The tracking error of the second subsystem and its derivative are given as

(47) \begin{align}e_{2}=x_2-x_{2d}, \dot{e_{2}}=\dot{x}_2-\dot{x}_{2d}\\[-18pt]\nonumber\end{align}

Then, the error dynamics ( $e_{1}$ , $e_{2}$ ) are written as follows:

(48) \begin{equation}\left\{\begin{array}{rl}\dot{e}_1 &= e_{1}+{x}_{2d}-\dot{x}_{1d}\\[4pt]\dot{e}_2 &=\beta_{1}(t)x_{4}x_{6}+\beta_{4}(t){x_{2}}^2+\beta_{2}(t)x_{4}{\Omega}_r+\beta_{3}(t)u_{2}+k_1(x_2-\dot{x}_{1d})-\ddot{x}_{1d}\end{array}\right.\end{equation}

To obtain the control law $u_2$ , we increase the already existing Lyapunov function (44) by adding a quadratic term.

Consider the new Lyapunov function candidate,

(49) \begin{equation}V_{2}(e_1,e_2)=V_{1}(e_1)+\dfrac{1}{2}{e_2}^{2}.\end{equation}

Its derivative is negative definite as

(50) \begin{equation}\dot{V}_{2}(e_1,e_2)=-k_1{e_1}^{2}-k_2{e_2}^{2}<0, k_2>0.\end{equation}

To achieve objective (50) and ensure Lyapunov stability, the control law is designed as

(51) \begin{equation}u_2 =\dfrac{1}{\beta_{3}(t)}[-\beta_{1}(t)x_{4}x_{6}-\beta_{4}(t){x_{2}}^2-\beta_{2}(t)x_{4}{\Omega}_r+\ddot{x}_{1d}-k_1(-k_1e_1+e_2)-e_1-k_2e_2]\end{equation}

To extract the other controllers, we follow the same steps as above.

5.1 Flight scenario

In this scenario, we want to make a circular trajectory followed by the quadrotor where it changes its configuration every 25s. To achieve this, the quadrotor starts in the classic $``\texttt{X}"$ configuration and then changes shape to other special morphologies such as $``\texttt{H}"$ , $``\texttt{O}"$ , $``\texttt{Y}"$ , $``\texttt{YI}"$ and $``\texttt{T}"$ (Fig. 3).

Figure 7. Evolution of servomotor angles.

The simulation results are shown in Figs 7, 8, 9, 10, 11, 12 and 13.

Figure 8. Evolution of servomotor position errors.

Figure 9. Evolution of 3D trajectory and positions.

Figure 10. Evolution of position errors.

5.1.1 Results and discussion for the flight scenario

Figure 7 shows the evolution of the real and desired positions of the servomotors used to rotate the four arms of the quadrotor. It is observed that the outputs of the four servomotors ( $\psi_1,\psi_2,\psi_3,\psi_4$ ) achieve their desired positions ( $\psi_{1d},\psi_{2d},\psi_{3d},\psi_{4d}$ ) at about 3s, without any overshoot. The steady-state errors are approximatively zero. To illustrate this result, the servomotor position errors are plotted in Fig. 8.

Figure 11. Evolution of Euler angles.

Figure 12. Evolution of attitude errors.

Figure 13. Evolution of control signals.

The evolution of the 3D trajectory (Fig. 9) shows that the vehicle moves in the horizontal plane (x,y) with a slow variation in altitude. Moreover, the altitude z and the positions (x,y) track their desired trajectories $z_d$ , $x_d$ , $y_d$ .

For the y-axis, the error varies approximately between $-0.14$ and $0.12$ m, while for the x-axis, it reaches 1m (approximatively equivalent to 50%) at the beginning of the trajectory tracking. However, this error decreases rapidly to reach about $0.18$ m (approximatively equivalent to 9%). The error on the z-axis is equal to 1m at $t=0$ s and then rapidly tends to zero after about $1.8$ s (Fig. 10).

To reach the desired altitude $z_d$ , the corresponding control increases rapidly at the beginning, and as soon as the rates of variation of the desired and actual altitudes are equal, the control takes a constant value of approximatively 11N (Fig. 13).

From Figs 11 and 12, it is observed that the yaw response follows its desired state, where the tracking error converges towards zero. The yaw dynamics is well illustrated by the evolution of the $u_4$ control (Fig. 13), increasing at the start and then approaching zero.

For the roll $\phi$ and pitch $\theta$ dynamics, the drone remains stable, with some oscillations appearing in the curves of the angles in Fig. 11. These oscillations produce errors along the two axes (x,y) during the whole flight period, where their values are around zero (Fig. 12).

The evolution of the $u_2$ and $u_3$ controls is shown in Fig. 13, where $u_2$ shows a negligible variation over time while $u_3$ takes a maximum value of $4.1$ N at $t = 0s$ , then for the rest of the time varies sinusoidally around zero.

In this scenario, we change the configuration of the quadrotor every 25s. The simulation results are satisfactory in terms of tracking, speed and accuracy.