Nomenclature

2.0 Nonlinear aircraft model

Quaternion method is a more feasible method than Euler angle method, which avoids the singularity of Euler angular velocity equation when pitch angle $\theta $ is equal to ${90^ \circ }$ . This method is also called Euler four parameter method [Reference Pamadi12]. The quaternion model has been derived in detail in reference [Reference Zhuang, Zhang, Sun and Chen13]. For the convenience of readers, the derivation process of the model is described in detail as follows.

We consider the vertical takeoff and dive aircraft system described by quaternions.

(1) \begin{equation}\left[ {\begin{array}{c}{{{\dot q}_0}}\\[2pt]{{{\dot q}_1}}\\[2pt]{{{\dot q}_2}}\\[2pt]{{{\dot q}_3}}\end{array}} \right] = \frac{1}{2}\left[ {\begin{array}{r@{\quad}r@{\quad}r}{ - {q_1}} & { - {q_2}} & { - {q_3}}\\[2pt]{{q_0}} & { - {q_3}} & {{q_2}}\\[2pt]{{q_3}} & {{q_0}} & { - {q_1}}\\[2pt]{ - {q_2}} & {{q_1}} & {{q_0}}\end{array}} \right]\left[ {\begin{array}{c}\boldsymbol{p}\\[2pt]\boldsymbol{q}\\\boldsymbol{r}\end{array}} \right]\end{equation}

(2) \begin{equation}\left[ \begin{array}{c}\dot{\textbf{{p}}}\\\dot{\textbf{{q}}}\\\dot{\textbf{{r}}}\end{array} \right] = {\left(\frac{V}{{{V_0}}}\right)^2} + \left[ \begin{array}{c}{l_\beta }\beta + {l_q}q + {l_r}r + ({l_\beta }_\alpha \beta + {l_{r\alpha }})\Delta\alpha \\[3pt] + {l_p}p + {l_{{\delta _a}}}{\delta _a} + {l_{{\delta _r}}}{\delta _r}\\\hline{m_\alpha }\Delta \alpha + {m_q}q - {m_{\dot{\alpha}} }p\beta\\[3pt] + {m_v}\Delta V + {m_{{\delta _e}}}{\delta _e}\\[3pt] + {m_\alpha }({g_0}/V) \times (\cos \theta \cos \varphi - \cos {\theta _0})\\[3pt]{n_\beta }\beta + {n_r}r + {n_p}p + {n_{p\alpha }}p\Delta \alpha\\\hline + {n_q}q + {n_{{\delta _a}}}{\delta _a} + {n_{{\delta _r}}}{\delta _r}\end{array} \right] + \left[ \begin{array}{r} - {i_1}\textbf{{qr}}\\[3pt]{i_2}\textbf{{pr}}\\[3pt] - {i_3}\textbf{{pq}}\end{array} \right]\end{equation}

where ${l_a} = \dfrac{1}{2}\rho V_0^2{S_r}{l_1}\dfrac{{\partial {C_l}}}{{\partial a}}/{I_x}$ , ${m_a} = \dfrac{1}{2}\rho V_0^2{S_r}{l_2}\dfrac{{\partial {C_m}}}{{\partial a}}/{I_y}$ , and ${n_a} = \dfrac{1}{2}\rho V_0^2{S_r}{l_3}\dfrac{{\partial {C_n}}}{{\partial a}}/{I_z}$ ( $a = \left\{ {\alpha ,\beta ,q,p, \cdots } \right\})$ and so on, the relationships between these variables are detailed in Table 1 of [Reference Zhuang, Sun, Chen and Zeng10]. And ${i_1} = \left( {{I_z} - {I_y}} \right)/{I_x}$ , ${i_2} = \left( {{I_z} - {I_x}} \right)/{I_y}$ , and ${i_3} = \left( {{I_y} - {I_x}} \right)/{I_z}$ . Here we introduce the factor ${(V/{V_0})^2}$ into rotational motion on account of these parameters are proportional to the square of velocity.

The unit-quaternion is a vector defined by ${\left[ {\begin{array}{l}{{q_0}} \quad {{q_1}} \quad {{q_2}}\quad {{q_3}}\end{array}} \right]^T}$ that satisfies $q_0^2 + q_1^2 + q_2^2 + q_3^2 = 1$ . Unit quaternion provides a convenient mathematical notation for representing the orientation and rotation of 3D objects. They are easier to compose than Euler angles and avoid the problem of gimbal locking. They are numerically more stable and possibly more efficient than rotation matrices.

Assumption 1. It is can be hypothesised that we can obtain the quaternion $q$ and measure angular velocity ${\left[ {\begin{array}{*{20}{l}}\boldsymbol{p} \quad \boldsymbol{q} \quad \boldsymbol{r}\end{array}} \right]^T}$ in aircraft systems (1) and (2); and we need to measure the displacement components ${\left[ {\begin{array}{*{20}{l}}x\quad y\quad z\end{array}} \right]^T}$ and velocity $V\!\left( {V = \sqrt {V_x^2 + V_y^2 + V_z^2} } \right)$ in (1) and (2), and they are specifically contained in the following equation:

(3) \begin{equation}\left[ \begin{array}{c}{{{\dot V}_x}} \\[2pt] {{{\dot V}_y}} \\[2pt] {{{\dot V}_z}} \\[2pt] \end{array} \right] = \left[ \begin{array}{c} ({T_x} - D)/m - 2({q_0}{q_3} + {q_1}{q_2})g \\[3pt] Q/m - \left[1 - 2\!\left(q_1^2 + q_3^2\right)\right]g \\[3pt] L/m - 2({q_2}{q_3} - {q_0}{q_1})g\end{array}\right] + \left[ \begin{array}{c@{\quad}c@{\quad}c} 0 & \boldsymbol{r} & { - \boldsymbol{q}} \\[3pt] { - \boldsymbol{r}} & 0 & \boldsymbol{p} \\[3pt] \boldsymbol{q} & { - \boldsymbol{p}} & 0 \end{array}\right]\left[ \begin{array}{c} {{V_x}} \\[3pt] {{V_y}} \\[3pt] {{V_z}}\end{array} \right]\end{equation}

After derivation, the coordinate transformation matrix from the airframe coordinate system to the inertial coordinate system determined by quaternion is

(4) \begin{equation}\left[ \begin{array}{l}{\dot x}\\[3pt]{\dot y}\\[3pt]{\dot z}\end{array} \right] = \left[ \begin{array}{c@{\quad}c@{\quad}c}1 - 2\!\left(q_2^2 + q_3^2\right) & 2\!\left({q_1}{q_2} - {q_0}{q_3}\right) & 2\!\left({q_1}{q_3} - {q_0}{q_2}\right)\\[3pt]2\!\left({q_0}{q_3} - {q_1}{q_2}\right) & 1 - 2\!\left(q_1^2 + q_3^2\right) & 2\!\left({q_2}{q_3} - {q_0}{q_1}\right)\\[3pt]2\!\left({q_1}{q_3} - {q_0}{q_2}\right) & 2\!\left({q_2}{q_3} - {q_0}{q_1}\right) & 1 - 2\!\left(q_1^2 - q_2^2\right)\end{array} \right]\left[ {\begin{array}{c}{{V_x}}\\[3pt]{{V_y}}\\[3pt]{{V_z}}\end{array}} \right]\end{equation}

The aim of this paper is to resolve the problem of attitude tracking; it is the purpose of this paper to design a feedback controller which is bounded and satisfies $\parallel\! q\!\left( t \right)\!\parallel = 1$ that enables the state of the CLS (1) and (2) to track the expected attitude motion ${q_r}$ , which can be described in the following:

(5) \begin{align}{\mathop {{\textrm{lim}}}\limits_{t \to \infty } \xi \!\left( t \right)} & = 0\nonumber\\[5pt]{\xi \!\left( t \right)} & = q\!\left( t \right) - {q_r}\left( t \right)\end{align}

and the deflection angle vector ${\left[ {\begin{array}{*{20}{l}}{{\delta _a}}\quad {{\delta _r}}\quad {{\delta _e}}\end{array}} \right]^T}$ is the control input.

For ease of description, we re-define the aircraft model as follow:

(6) \begin{align}{{Y_1}} & = {{\left[ {\begin{array}{*{20}{l}}{{q_0}}\quad {{q_1}}\quad {{q_2}}\quad {{q_3}}\end{array}} \right]}^T}\nonumber\\[5pt]{{Y_2}} & = {{\left[ {\begin{array}{*{20}{l}}\textbf{{p}}\quad \textbf{{q}}\quad \textbf{{r}}\end{array}} \right]}^T}\nonumber\\[5pt]{{Y_3}} & = {{\left[ {\begin{array}{*{20}{l}}x\quad y\quad z\end{array}} \right]}^T}\\[5pt]{{Y_4}} & = {{\left[ {\begin{array}{*{20}{l}}{{V_x}}\quad {{V_y}}\quad {{V_z}}\end{array}} \right]}^T}\nonumber\\[5pt]U & \buildrel \Delta \over = {{\left[ {\begin{array}{*{20}{l}}{{\delta _a}}\quad {{\delta _r}}\quad {{\delta _e}}\end{array}} \right]}^T}\nonumber\end{align}

Then we rewrite the aircraft model (1) and (2) as follow:

(7) \begin{align}{{{\dot Y}_1}} & = {H_1}\!\left( {{Y_1}} \right){Y_2}\nonumber\\[5pt]{{{\dot Y}_2}} & = {H_2}\!\left( {{Y_2}} \right) + {H_3}\!\left( {{Y_1},{Y_2},{Y_3},{Y_4}} \right) + B\!\left( {{Y_1},{Y_3},{Y_4}} \right)U\end{align}

where

(8) \begin{equation}{H_1}\!\left( {{Y_1}} \right) = \frac{1}{2}\!\left[ {\begin{array}{r@{\quad}r@{\quad}r}{ - {q_1}} & { - {q_2}} & { - {q_3}}\\[5pt]{{q_0}} &{ - {q_3}} & {{q_2}}\\[5pt]{{q_3}} & {{q_0}} &{ - {q_1}}\\[5pt]{ - {q_2}} & {{q_1}} &{{q_0}}\end{array}} \right]\end{equation}

(9) \begin{equation}{H_2}\!\left( {{Y_2}} \right) = \left[ {\begin{array}{r}{ - \left( {{I_z} - {I_y}} \right)/{I_x}\textbf{{qr}}}\\[5pt]{\left( {{I_z} - {I_x}} \right)/{I_y}\textbf{{pr}}}\\[5pt]{ - \left( {{I_y} - {I_x}} \right)/{I_z}\textbf{{pq}}}\end{array}} \right]\end{equation}

(10) \begin{equation}{H_3}\!\left( {{Y_1},{Y_2},{Y_3},{Y_4}} \right) = \frac{1}{2}\rho V_0^2{S_r}{\mathcal{L}} \times \left[ \begin{array}{c} {\Big[{\dfrac{\partial {C_l}}{\partial \beta }}\beta + {\dfrac{\partial {C_l}}{\partial q}}q + {\dfrac{\partial {C_l}}{\partial r}}r} \\[15pt] { + \left( {{\dfrac{\partial {C_l}}{\partial \beta \alpha }}\beta + {\dfrac{\partial {C_l}}{\partial r\alpha }}} \right){{\Delta }}\alpha + {\dfrac{\partial {C_l}}{\partial p}}p]/{I_x}} \\[-6pt]\\\hline\\[-5pt] {\left\{ {\dfrac{\partial {C_m}}{\partial \alpha }}{{\Delta }}\alpha + {\dfrac{\partial {C_m}}{\partial q}}q - {\dfrac{\partial {C_m}}{\partial \dot \alpha }}p\beta + {\dfrac{\partial {C_m}}{\partial V}}{{\Delta }}V + \right.} \\[12pt] {{\dfrac{\partial {C_m}}{\partial \dot \alpha }}\left( {{g_0}/V} \right) \times \left( {{\textrm{cos}}\theta {\textrm{cos}}\phi - {\textrm{cos}}{\theta _0}} \right)\} /{I_y}} \\[-6pt]\\\hline\\[-5pt] {\left[ {{\dfrac{\partial {C_n}}{\partial \beta }}\beta + {\dfrac{\partial {C_n}}{\partial r}}r + {\dfrac{\partial {C_n}}{\partial p}}p + {\dfrac{\partial {C_n}}{\partial p\alpha }}p{{\Delta }}\alpha + {\dfrac{\partial {C_n}}{\partial q}}q} \right]/{I_z}} \end{array} \right]\end{equation}

(11) \begin{equation}\begin{array}{*{20}{l}}{B\!\left( {{Y_1},{Y_3},{Y_4}} \right) = \dfrac{1}{2}\rho V_0^2{S_r}{\mathcal{L}} \times \left[ {\begin{array}{c@{\quad}c@{\quad}c}{\dfrac{{\partial {C_l}}}{{\partial {\delta _a}}}/{I_x}} & {\dfrac{{\partial {C_l}}}{{\partial {\delta _r}}}/{I_x}} & 0\\[5pt]0 & 0 & {\dfrac{{\partial {C_m}}}{{\partial {\delta _e}}}/{I_y}}\\[5pt]{\dfrac{{\partial {C_n}}}{{\partial {\delta _a}}}/{I_z}} & {\dfrac{{\partial {C_n}}}{{\partial {\delta _r}}}/{I_z}} & 0\end{array}} \right]}\end{array}\end{equation}

It can be clearly seen from the above formula that many uncertain factors are involved in the whole system. For example, the atmospheric moment coefficient $\dfrac{{\partial {C_l}}}{{\partial {*}}}$ , $\dfrac{{\partial {C_m}}}{{\partial {\textrm{*}}}}$ and $\dfrac{{\partial {C_n}}}{{\partial {\textrm{*}}}}$ depends on the change of Mach number ${M_a}$ . The Mach number ${M_a}$ is related to the state variable. During the actual flight of the aircraft, it is difficult to determine these atmospheric moment coefficients, which leads to the uncertainty of the aircraft model. Therefore, because of the uncertainty of the atmospheric moment coefficient $\dfrac{{\partial {C_l}}}{{\partial {\textrm{*}}}}$ , $\dfrac{{\partial {C_m}}}{{\partial {\textrm{*}}}}$ and $\dfrac{{\partial {C_n}}}{{\partial {\textrm{*}}}}$ , ${H_3}$ and $B$ are unknown, which makes the control design more complex, and the system structure (7) is particularly difficult [Reference Xia, Lu, Zhu and Fu15].

In order to settle this puzzle, we introduce a new variable $F\!\left( t \right)$ as follow

(12) \begin{equation}F\!\left( t \right) = {H_3}\!\left( {{Y_1},{Y_2},{Y_3},{Y_4}} \right) + B\!\left( {{Y_1},{Y_3},{Y_4}} \right)U - {B_0}U\end{equation}

where ${B_0}$ is defined as

(13) \begin{equation}{B_0} = \frac{1}{2}\rho V_0^2{S_r}{\mathcal{L}}{{\Phi }}{|_{{M_a} = const}}\end{equation}

where ${{\Phi }}{|_{{M_a} = const}}$ is defined as

\begin{equation*}{{\Phi }}{\Big|_{{M_a} = const}} = \left[ {\begin{array}{c@{\quad}c@{\quad}c}{\dfrac{{\partial {C_l}}}{{\partial {\delta _a}}}/{I_x}} & {\dfrac{{\partial {C_l}}}{{\partial {\delta _r}}}/{I_x}} & 0\\[9pt]0 & 0 & {\dfrac{{\partial {C_m}}}{{\partial {\delta _e}}}/{I_y}}\\[9pt]{\dfrac{{\partial {C_n}}}{{\partial {\delta _a}}}/{I_z}} & {\dfrac{{\partial {C_n}}}{{\partial {\delta _r}}}/{I_z}} & 0\end{array}} \right]{\Big|_{{M_a} = const}}\end{equation*}

We can obtain the coefficients $\dfrac{{\partial {C_l}}}{{\partial {\delta _a}}}$ , $\dfrac{{\partial {C_l}}}{{\partial {\delta _r}}}$ , $\dfrac{{\partial {C_m}}}{{\partial {\delta _e}}}$ , $\dfrac{{\partial {C_n}}}{{\partial {\delta _a}}}$ , and $\dfrac{{\partial {C_n}}}{{\partial {\delta _r}}}$ by selecting the appropriate constant Maher number ${M_a}$ . These coefficients are part of ${B_0}$ , and decompose the uncertain part into variable $F\!\left( t \right)$ . Thus, the dynamic uncertainty existing in (7) can be reduced to the total uncertainty $F\!\left( t \right)$ , hence, the complicacy of control design is also reduced. Therefore, we can rewrite the system (7) as follow:

(14a) \begin{equation}{\dot Y_1} = {H_1}\!\left( {{Y_1}} \right){Y_2}\end{equation}

(14b) \begin{equation}{\dot Y_2} = H_2\!\left( {{Y_2}} \right) + F\!\left( t \right) + {B_0}U\!\left( t \right)\end{equation}

To proceed the controller design and stability proof, the following lemmas are given.

Lemma 1. [Reference Ioannou and Sun16] Let $f,{{\;\;}}V:\left[ {0,\infty } \right) \to \mathbb{R}$ . Then

\begin{equation*}\dot V \le - \alpha V + f,{{\;\;\;\;\;\;\;\;}}\forall t \geqslant {t_0} \geqslant 0\end{equation*}

implies that

\begin{equation*} {V\!\left( t \right) \le } {e^{ - \alpha \!\left( {t - {t_0}} \right)}}V\!\left( {{t_0}} \right) + \mathop \smallint \nolimits_{{t_0}}^t {e^{ - \alpha \!\left( {t - \tau } \right)}}f\!\left( \tau \right)d\tau ,\quad \forall t \ge {t_0} \ge 0\end{equation*}

for any finite constant $\alpha $ .

Lemma 2. [Reference Guo and Zhao17] The linear tracking differentiator is considered as follow:

(15) \begin{equation}\left\{ {\begin{array}{*{20}{l}}{{{\dot \chi }_1} = {\chi _2}}\\[5pt]{{{\dot \chi }_2} = - {\kappa _1}{R^2}\left( {{\chi _1} - {Y_r}} \right) - {\kappa _2}R{\chi _2}}\end{array}} \right.\end{equation}

where ${\kappa _1} \gt 0$ and ${\kappa _2} \gt 0$ are constants, and $R \gt 0$ is the tuning parameter.

Suppose that ${\kappa _1} \gt 0,{\kappa _2} \gt 0$ and ${Y_r}\,:\,\left[ {0,\infty } \right) \to R$ is a function satisfying $\mathop {{\textrm{sup}}}\limits_{t \in \left[ {0,\infty } \right)} \left( {\left| {{Y_r}\left| + \right|{{\dot Y}_r}} \right|} \right) = {M_1} \lt \infty $ for constant ${M_1} \gt 0$ . Then the linear tracking differentiator (15) is convergent in the sense that, for $\forall a \gt 0$ ,

\begin{equation*}\left\{ {\begin{array}{*{20}{l}}{\mathop {{\textrm{lim}}}\limits_{R \to \infty } \left| {{\chi _1} - {Y_r}} \right| = 0}\\[9pt]{\mathop {{\textrm{lim}}}\limits_{R \to \infty } \left| {{\chi _2} - {{\dot Y}_r}} \right| = 0}\end{array}} \right.\end{equation*}

uniformly for $t \in \left[ {a,\infty } \right)$ .

Lemma 3. [Reference Buckholtz18] The slope of the saturation function, $sat\!\left( {\dfrac{s}{{{{\Phi }}\left| \eta \right|}}} \right)$ , increases as the approach angle magnitude, $\left| \eta \right|$ , decreases; hence,

\begin{equation*}{\mathop {{\textrm{lim}}}\limits_{\left| \eta \right| \to 0} sat\!\left( {\frac{s}{{{{\Phi }}\left| \eta \right|}}} \right) = sign\!\left( s \right)}\end{equation*}

3.0 Back-stepping adaptive smc based on rbf network approximation design

3.1. Back-stepping procedure

Next, we will design the BSM of the system (14): first, let’s define some variables, as shown in Table 2. Next, we start with the definition of state error ${z_1}$

\begin{equation*}{{z_1} = {Y_1} - {Y_{r,1}}}\end{equation*}

where the reference value of ${Y_r}$ is ${Y_{r,1}}$ , next, take the derivative of ${z_1}$ , and we get

(16) \begin{equation}{{{\dot z}_1} = {{\dot Y}_1} - {{\dot Y}_r} = {H_1}\!\left( {{Y_1}} \right){Y_2} - {{\dot Y}_r}}\end{equation}

Based on the principle of BSM, ${Y_2}$ is regarded as the virtual control input, used for applying the following dynamic expectations

(17) \begin{equation}{\dot z_1} = - {{\Lambda }}{z_1} = - {\textrm{diag}}\!\left[ {\begin{array}{*{20}{l}}{{\lambda _1}}\quad {{\lambda _2}}\quad {{\lambda _3}} \quad{{\lambda _4}}\end{array}} \right]{z_1}\end{equation}

To ensure the asymptotic stability of (17), the design matrix ${{\Lambda }}$ is chosen as ${\lambda _i} \gt 0$ , $i = 1,2,3,4$ . Therefore, combining (16) with (17), we can get the solution

(18) \begin{equation}{Y_{r,2}} = H_{1L}^{ - 1}\left( {{Y_1}} \right)\left( {{{\dot Y}_r} - {{\Lambda }}{z_1}} \right)\end{equation}

where the left inverse matrix of ${H_1}\!\left( {{Y_1}} \right)$ is denoted by $H_{1L}^{ - 1}\left( {{Y_1}} \right)$ .

Table 2. The variable definition in BSM

Remark 3. In (8), ${H_1}\!\left( {{Y_1}} \right)$ has four ${\boldsymbol{R}^{3 \times 3}}$ subdeterminants. ${q_0}$ , ${q_1}$ , ${q_2}$ and ${q_3}$ are respectively the values of these subdeterminants. The column rank of ${H_1}\!\left( {{Y_1}} \right)$ is reduced if and only if ${q_0}{ = q_1} = {q_2} = {q_3} = 0$ , otherwise, if any ${q_i}$ satisfies ${q_i} \ne 0$ , then the column rank of ${H_1}\!\left( {{Y_1}} \right)$ is full. Note $q_0^2 + q_1^2 + q_2^2 + q_3^2 = 1$ , ${q_0} = {q_1} = {q_2} = {q_3} = 0$ cannot happen at the same time. Hence ${H_1}\!\left( {{Y_1}} \right)$ is full column rank, which means the existence of the left inverse matrix of ${H_1}\!\left( {{Y_1}} \right)$ [Reference Xia, Shi, Liu, Rees and Han5].

3.2. RBF network approximation

For the sake of overcoming the interference inside the system and obtain the accurate modeling information $F(t)$ of the controlled object, fractional order SMC and BSM based on RBF network approximation should be combined.

Then, we design the sliding mode surface (SMS) as follow:

(19) \begin{equation}s = {z_2}\end{equation}

Next, consider the following reaching law

(20a) \begin{equation}\dot s = - \Pi s - \Sigma {\textrm{sig}}{{\textrm{n}}^\varsigma }\left( s \right)\end{equation}

(20b) \begin{equation}\Pi = {\textrm{diag}}\!\left[ {\begin{array}{*{20}{l}}{{\pi _1}}\quad {{\pi _2}}\quad {{\pi _3}}\end{array}} \right],{\pi _i} \gt 0\end{equation}

(20c) \begin{equation}\Sigma = {\textrm{diag}}\!\left[ {\begin{array}{*{20}{l}}{{\sigma _1}}\quad {{\sigma _2}}\quad {{\sigma _3}}\end{array}} \right],{\sigma _i} \gt 0\end{equation}

where ${\textrm{sig}}{{\textrm{n}}^\varsigma }\left( s \right) = [{\textrm{sign}}\left( {{s_1}} \right)\left| {{s_1}{|^\varsigma },{\textrm{sign}}\left( {{s_2}} \right)} \right|{s_2}{|^\varsigma },{\textrm{sign}}\left( {{s_3}} \right)|{s_3}{|^\varsigma }]$ with $0 \lt \varsigma \lt 1$ .

Taking the derivative of $s$ , and combining the reaching law (20a), we get

(21) \begin{align}{\dot s} & = {H_2}\!\left( {{Y_2}} \right) + F\!\left( t \right) + {B_0}U\!\left( t \right) - {{\dot Y}_{r,2}}\nonumber\\[5pt] & = - {{\Pi }}s - {{\Sigma sig}}{{\textrm{n}}^\varsigma }\left( s \right)\end{align}

we can acquire the control law by solving for $U\!\left( t \right)$ in equality (21)

(22) \begin{equation}{U\!\left( t \right) = B_0^{ - 1}\left( { - {{\Pi }}s - {{\Sigma sig}}{{\textrm{n}}^\varsigma }\left( s \right) - {H_2}\!\left( {{Y_2}} \right) - F\!\left( t \right) + {{\dot Y}_{r,2}}} \right)}\end{equation}

Note that, the controller (25) contains the uncertainty $F\!\left( t \right)$ , which is partially unknown to us; Therefore, it cannot be directly applied to the actual system before obtaining $F\!\left( t \right)$ . Next, we need to utilise the RBF NN to approximate $F\!\left( t \right)$ .

The uncertainty $F\!\left( t \right)$ is approximated adaptively by RBF NN. The RBF NN algorithm is:

(23) \begin{align}{R_j} & = {exp\!\left( { - \dfrac{{\parallel {\textbf{t}} - {c_j}{\parallel ^2}}}{{2b_j^2}}} \right),{{\;\;\;\;\;\;}}j = 1,2, \cdots ,m}\nonumber\\[8pt]F\!\left( t \right) & = {{\textbf{WR}}\!\left( {\textbf{t}} \right) + \varepsilon }\end{align}

where t is input signal of network; $j$ is the number of hidden layer nodes in the network; ${c_j}$ is the center of the basis function; ${b_j}$ is the width of the Gaussian function of the $j$ th hidden layer unit; ${\textbf{R}} = {[{r_1},{r_2}, \cdots ,{r_m}]^T}$ is the output of the Gaussian function; ${\textbf{W}} \in {R^{3 \times m}}$ is the ideal neural network weights; $\varepsilon $ is the approximation error of neural network, $\varepsilon = {[{\varepsilon _1},{\varepsilon _2},{\varepsilon _3}]^T}$ , $\left| {{\varepsilon _i}} \right| \le {\varepsilon _N}$ .

$F\!\left( t \right)$ is approximated by the RBF NN, and the output of RBF NN is

(24) \begin{equation}\hat F\!\left( {\textbf{t}} \right) = \hat{\textbf{W}}{{\textbf{R}}\left( {\textbf{t}} \right)}\end{equation}

where, $\hat F\!\left( t \right)$ is the approximation of the RBF NN. We adopt the NN minimum parameter learning method, let $\varphi = \parallel {\textbf{W}}{\parallel ^2}$ , $\varphi \gt 0$ , the estimation of $\varphi $ is $\hat \varphi $ , $\tilde \varphi = \hat \varphi - \varphi $ .

Remark 4. In the actual control system design, for guaranteeing that the input value of the NN is a Gaussian basis function (GBF) within the valid range, the coordinate vector ${c_j}$ of the center point of the GBF should be determined according to the actual range of the network input value. In order to ensure the effective mapping of the GBF, the width ${b_j}$ of the GBF needs to be taken as an appropriate value. By proving the stability of the closed-loop Lyapunov function, the adjustment of $\hat W$ is designed.

Using adaptive RBF NN to approximate the total uncertainty of $F\!\left( t \right)$ , the controller (22) is rewritten as

(25) \begin{equation}{U\!\left( t \right) = B_0^{ - 1}\left( { - {{\Pi }}s - {{\Sigma sig}}{{\textrm{n}}^\varsigma }\left( s \right) - {H_2}\!\left( {{Y_2}} \right) - \frac{1}{2}s\hat \varphi {{\textbf{R}}^T}{\textbf{R}} + {{\dot Y}_{r,2}}} \right)}\end{equation}

4.0 Main results

Theorem 1. Consider the subsystem (14a), and SMS given by (19). If the RBF neural network algorithm and control law are respectively established as (23) and (25), then the plant trajectories asymptotically converge into a small region of SMS, and the global asymptotic stability of the CLS (14b) is guaranteed within a finite reaching time by defining the adaptive control laws:

(26) \begin{equation} {\dot{\hat{\varphi}} = \frac{\rho}{2}\parallel\!{s}{\parallel ^2}{\boldsymbol{R}^T}\boldsymbol{R} - \kappa \rho \hat \varphi }\end{equation}

where $\kappa \gt 0$ .

Proof. A Lyapunov function candidate is defined

(27) \begin{equation}{L = {{\frac{1}{2}}^2} + \frac{1}{{2\rho }}{{\tilde \varphi }^2}}\end{equation}

where $\rho \gt 0$ .

With (21) and (25), the derivative of $L$ can be calculated by

(28) \begin{align} {\dot L} & = {s\dot s + {1 \over \rho }\tilde \varphi {\dot{\hat{\varphi}}} } \nonumber\\[5pt] & = {s\!\left( {{\textbf{WR}} + \varepsilon - \Pi s - {{\Sigma sig}}{{\textrm{n}}^\varsigma }\left( s \right) - \frac{1}{2}s\hat \varphi {{\textbf{R}}^T}{\textbf{R}}} \right) + {1 \over \rho }\tilde \varphi {\dot{\hat{\varphi}}} } \nonumber\\[5pt] & \le {\frac{1}{2}\parallel\!{s}{\parallel ^2}\varphi {{\textbf{R}}^T}{\textbf{R}} + \frac{1}{2} - \frac{1}{2}\parallel\!{s}{\parallel ^2}\hat \varphi {{\textbf{R}}^T}{\textbf{R}} + \parallel\!{\varepsilon}\!\parallel \parallel\!{s}\!\parallel - \parallel\!{\Sigma}\!\parallel \parallel\!{s}\!\parallel + {1 \over \rho }\tilde \varphi {\dot{\hat{\varphi}}} - \parallel\!{\Pi}\!\parallel \parallel\!{s}{\parallel ^2}} \nonumber\\[5pt]& = { - \frac{1}{2}\parallel\!{s}{\parallel ^2}\tilde \varphi {{\textbf{R}}^T}{\textbf{R}} + \frac{1}{2} + \parallel\!{\varepsilon}\!\parallel \parallel\!{s}\!\parallel - \parallel\!{\Sigma}\!\parallel \parallel\!{s}\!\parallel + {1 \over \rho }\tilde \varphi {\dot{\hat{\varphi}}} - \Pi {s^2}} \nonumber\\[5pt]& = {\tilde \varphi \!\left( { - \frac{1}{2}\parallel\!{s}{\parallel ^2}{{\textbf{R}}^T}{\textbf{R}} + {1 \over \rho }{\dot{\hat{\varphi}}} } \right) + \frac{1}{2} + \parallel\!{\varepsilon}\!\parallel \parallel\!{s}\!\parallel } \nonumber\\[5pt] & { - \parallel\!{\Sigma}\!\parallel \parallel\!{s}\!\parallel - \parallel\!{\Pi}\!\parallel \parallel\!{s}{\parallel ^2}} \nonumber\\[5pt] & \le {\tilde \varphi \!\left( { - \frac{1}{2}\parallel\!{s}{\parallel ^2}{{\textbf{R}}^T}{\textbf{R}} + {1 \over \rho }{\dot{\hat{\varphi}}} } \right) + \frac{1}{2} - \parallel\!{\Pi}\!\parallel \parallel\!{s}{\parallel ^2}}\end{align}

Substituting the adaptive control law (26) into equation (28), it can be acquired

(29) \begin{align}{\dot L \le } & { - \kappa \tilde \varphi \hat \varphi + \frac{1}{2} - \Pi {s^2}} \nonumber\\[5pt] \le & { - {\kappa \over 2}\left( {{{\tilde \varphi }^2} - {\varphi ^2}} \right) + \frac{1}{2} - \Pi {s^2}} \\[5pt] = & { - {\kappa \over 2}{{\tilde \varphi }^2} - \Pi {s^2} + \frac{1}{2}\left( {\kappa {\varphi ^2} + 1} \right)} \nonumber\end{align}

Let $\kappa = \dfrac{{2\parallel\!{\Pi}\!\parallel }}{\rho }$ , next, we get

(30) \begin{align} {\dot L \le } & { - {{\parallel {{\Pi }}\parallel } \over \rho }{{\tilde \varphi }^2} - \parallel {{\Pi }}\parallel {s^2} + \left( {{\kappa \over 2}{\varphi ^2} + \frac{1}{2}} \right)} \nonumber\\[5pt] = & { - 2\parallel {{\Pi }}\parallel \left( {{1 \over {2\rho }}{{\tilde \varphi }^2} + \frac{1}{2}{s^2}} \right) + \left( {{\kappa \over 2}{\varphi ^2} + \frac{1}{2}} \right)} \\[5pt] = & { - 2\parallel {{\Pi }}\parallel L + Q}\nonumber\end{align}

where $Q = \dfrac{\kappa }{2}{\varphi ^2} + \dfrac{1}{2}$ .

According to Lemma 1, solve inequality (30), it can be obtained

(31) \begin{equation}{L \le \frac{Q}{{2\parallel\!{\Pi}\!\parallel }} + \left( {L\!\left( 0 \right) - \frac{Q}{{2\parallel\!{\Pi}\!\parallel }}} \right){e^{ - 2\parallel\!{\Pi}\!\parallel t}}}\end{equation}

namely,

(32) \begin{align}{\mathop {{\textrm{lim}}}\limits_{t \to \infty } L }& = {\frac{Q}{{2\parallel\!{\Pi}\!\parallel }} = \frac{{\frac{\kappa }{2}{\varphi ^2} + \frac{1}{2}}}{{2\parallel\!{\Pi}\!\parallel }} = \frac{{\kappa {\varphi ^2} + 1}}{{4\parallel\!{\Pi}\!\parallel }}}\nonumber\\[5pt]& = {\frac{{\frac{{2\parallel\!{\Pi}\!\parallel }}{\rho }{\varphi ^2} + 1}}{{4\parallel\!{\Pi}\!\parallel }} = \frac{{{\varphi ^2}}}{{2\rho }} + \frac{1}{{4\parallel\!{\Pi}\!\parallel }}}\end{align}

Substituting (32) into Equation (30), it can be getted

(33) \begin{equation}\begin{array}{*{20}{l}}{\dot L \le } { - 2\parallel {{\Pi }}\parallel L + Q = 0}\end{array}\end{equation}

Hence, the state trajectory of the CLS (14b) is asymptotically stable under the application of the controller (25).

Proof. The proof process is divided into two steps.

Step 1: According to (17), the selection of ${{\Lambda }}$ is based on the principle of back-stepping, Lyapunov functional is chosen as follows:

(34) \begin{equation}{V_1} = \frac{1}{2}z_1^2\end{equation}

Then

(35) \begin{equation}{\dot V_1} = {z_1}{\dot z_1} = - \mathop \sum \limits_{i = 1}^4 {\lambda _i}z_{1i}^2 - z_1^T{H_1}\left( {{Y_1}} \right){z_2}\end{equation}

by selecting positive ${\lambda _i}$ large enough, we obtain ${\dot V_1} \lt 0$ when ${V_1}$ is out of a certain bounded region.

Step 2: Choosing the Lyapunov function as follows

(36) \begin{equation}{V_2} = \frac{1}{2}z_2^2\end{equation}

The time derivative of Lyapunov function (36) is

(37) \begin{equation}{\dot V}_2 = {z_2^T{{\dot z}_2}}\end{equation}

With (21) and (25), the derivative of ${V_2}$ can be calculated by

(38) \begin{align}{{\dot V}_2} & = {s\!\left( {WR + \varepsilon - \Pi s - \Sigma {\textrm{sig}}{{\textrm{n}}^\varsigma }\left( s \right) - \frac{1}{2}s\hat \varphi {{\textbf{R}}^T}{\textbf{R}}} \right)} \nonumber\\[5pt]& \le {\frac{1}{2}\parallel\!{s}{\parallel ^2}\varphi {{\textbf{R}}^T}{\textbf{R}} + \frac{1}{2} - \frac{1}{2}\parallel\!{s}{\parallel ^2}\hat \varphi {R^T}R + \parallel\!{\varepsilon}\!\parallel \parallel\!{s}\!\parallel - \parallel\!{\Sigma}\!\parallel \parallel\!{s}\!\parallel - \parallel\!{\Pi}\!\parallel \parallel\!{s}{\parallel ^2}} \nonumber\\[5pt]& = { - \frac{1}{2}\parallel\!{s}{\parallel ^2}\tilde \varphi {{\textbf{R}}^T}{\textbf{R}} + \frac{1}{2} + \parallel\!{\varepsilon}\!\parallel \parallel\!{s}\!\parallel - \parallel\!{\Sigma}\!\parallel \parallel\!{s}\!\parallel - {{\Pi }}{s^2}} \\[5pt]& \le { - \frac{1}{2}\parallel\!{s}{\parallel ^2}\tilde \varphi {{\textbf{R}}^T}{\textbf{R}} + \frac{1}{2} - \parallel\!{\Pi}\!\parallel \parallel\!{s}{\parallel ^2}}\nonumber\end{align}

Next, scaling down ${\dot V_2}$ , it can be acquired

(39) \begin{align}{{\dot V}_2} \lt & - \frac{1}{2}\parallel\!{s}{\parallel ^2}\tilde \varphi {{\textbf{R}}^T}{\textbf{R}} + {1 \over \rho }\tilde \varphi {\dot{\hat{\varphi}}} + \frac{1}{2} - \parallel\!{\Pi}\!\parallel \parallel\!{s}{\parallel ^2} \nonumber\\[5pt] & \le {\tilde \varphi \!\left( { - \frac{1}{2}\parallel\!{s}{\parallel ^2}{{\textbf{R}}^T}{\textbf{R}} + {1 \over \rho }{\dot{\hat{\varphi}}} } \right) + \frac{1}{2} - \parallel\!{\Pi}\!\parallel \parallel\!{s}{\parallel ^2}}\end{align}

According to Theorem 1, we know ${\dot V_2} \lt \dot L \le 0$ , hence, it can be obtained

(40) \begin{equation}{{{\dot V}_2} \lt 0}\end{equation}

In conclusion, define Lyapunov function $V = V_1 + {V_2}$ , $\dot V = {\dot V_1} + {\dot V_2} \lt 0$ .

Hence, the tracking errors ${z_1}$ and ${z_2}$ are driven by the controller (25) to converge to the neighbourhood of the origin and remain there for all subsequent times.

Remark 6. We can acquire the virtual control input ${Y_{r,2}}$ by calculating (18), but ${\dot Y_r}$ may not be acquired without difficulty because of the complex structure of $H_{1L}^{ - 1}\left( {{Y_1}} \right)$ . Here, the tracking differentiator (TD) [Reference Xia, Shi, Liu, Rees and Han5] is introduced so as to obtain the ${\dot Y_r}$ . Here’s a brief design of TD which is applied to track reference signal.

(41) \begin{equation}\left\{ {\begin{array}{*{20}{l}}{{{\dot \chi }_1} = {\chi _2}}\\[5pt]{{{\dot \chi }_2} = - {\kappa _1}{R^2}\left( {{\chi _1} - {Y_r}} \right) - {\kappa _2}R{\chi _2}}\end{array}} \right.\end{equation}

where ${\chi _1}$ and ${\chi _2}$ are the state variables of TD, ${\kappa _1} \gt 0,{\kappa _2} \gt 0$ are constants that denote the maximum actuation available in the system; and $R \gt 0$ is the tuning parameter.

In the light of Lemma 2, the TD states ${\chi _1}$ , ${\chi _2}$ will be approximate to ${Y_r}$ , ${\dot Y_r}$ , respectively as fast as possible. Xia, Y.Q. etc. have proved the advantages of this tracking method, for more details, please refer to Refs [Reference Xia, Shi, Liu, Rees and Han5, Reference Xia, Zhu and Fu6]. See Ref. [Reference Han19] for further explanation of TD. Thus ${\dot Y_r}$ can be gotten with the help of designing of TD for ${Y_r}$ .

Remark 7. In theory, we give an analysis to prove the stability of the CLS, which at least shows that this method is theoretically feasible. However, SMC has a disadvantage, that is, chattering caused by switching functions cannot be eliminated. For this reason, we use the saturation function ${\textbf{sat}}\!\left( {\dfrac{s}{{\Phi \!\left| \eta \right|}}} \right)$ instead of the switching function ${\textbf{sign}}\!\left( s \right)$ to weaken the chattering generated by the switching function. According to Lemma 3, we know that $\mathop {lim}\limits_{\left| \eta \right| \to 0} {\textbf{sat}}\left( {\dfrac{s}{{\Phi \!\left| \eta \right|}}} \right) = {\textbf{sign}}\left( s \right)$ . The specific theoretical proof and analysis are detailed in the literature [Reference Buckholtz18], we won’t repeat it here. In fact, it can be seen from Fig. 1 that as $\eta $ gets smaller, the saturation function gets closer to the switching function. Hence, the controller can be rewritten as

(42) \begin{equation}\hat U \left( t \right) = B_0^{ - 1}\left( - {{\Pi }}s - \Sigma {\textbf{sat}}\left( {s \over {\Phi \left| \eta \right|}}\right)^{\varsigma } - {H_2}\!\left( {{Y_2}} \right) - \hat F\!\left( t \right) + {{\dot Y}_{r,2}} \right)\end{equation}

Figure 1. Comparison of saturation function and switching function.

5.0 Simulation results and discussion

To demonstrate the validity and effectiveness of the proposed method in this paper, it was applied the Airbus-300 aircraft model [Reference Liu20]. See Table 3 for specific flight conditions and data.

Table 3. The aerodynamic derivatives computed at the trim condition.

Datas from: Rudolf B., Wolfgang A., Robert L., Flugregelung, Springer, Berlin, 2011.

It can be seen from the Sections 3 and 4 that the convergence speed of the state trajectory can be adjusted by these parameters ${{\Pi }}$ , ${{\Sigma }}$ , $\zeta $ , which can be tuned to weaken the chattering on the sliding surface according to Ref. [Reference Zhuang, Sun, Chen and Zeng10]. Table 4 shows the simulation parameters, which are selected by basis of the theorem and remark in the paper. That is to say, the controller (42) can make the system track converge to the predetermined SMS with fast speed and accurate response in the presence of uncertainty and disturbance, and maintained there for all subsequence time.

Table 4. Parameter value of simulation

It is hypothesised that the desired quaternion is as follows

(43) \begin{equation}{q_r} = {\left[ {\begin{array}{*{20}{l}}{{q_{0,r}}} \quad {{q_{1,r}}} \quad{{q_{2,r}}} \quad{{q_{3,r}}}\end{array}} \right]^T} = {\left[ {\begin{array}{*{20}{l}}{0.7{\textrm{sin}}\left( t \right)} \quad{{\textrm{cos}}\left( t \right)}\quad {2{\textrm{sin}}\left( t \right)}\quad {{\textrm{cos}}\left( t \right)}\end{array}} \right]^T}\end{equation}

and the initial attitude orientation of the aircraft is the same as in Ref. [Reference Chen and Huang21], which is $q\!\left( 0 \right) = {\left[ {\begin{array}{*{20}{l}}{0.3}\quad {0.2}\quad { - 0.3}\quad {0.8832}\end{array}} \right]^T}$ .

The attitude history for state $_0$ is shown in the top graph of Fig. 2. Similar curves for ${q_1},{q_2}$ and ${q_3}$ are achieved as well. The attitude quaternion tracking errors converges to zero as shown in the Fig. 3, which shows that the controller (42) achieves high-precision performance on the attitude stabilisation in the presence of uncertainties and disturbances. It is clear that the proposed aircraft control system can ensure that the attitude quaternions effectively track the commanded quaternions by choosing a reasonable positive parameter ${{\Lambda }} = 0.6{I_3}$ , where ${I_3}$ is the ${\mathbb{R}^{3 \times 3}}$ identity matrix. The control deflection angles is illustrated in the Fig. 4.

Figure 2. Profiles of the attitude.

Figure 3. The quaternion errors.

Figure 4. F(t) and approximation.

The sliding surface function parameters and corresponding approximated parameters are chosen in Table 4. As shown in Fig. 5, the sliding surface function reaches the reachable state of ${s_i} = 0,i = 1,2,3$ in 4s, and stays there all the time, which proves that the sliding surface selected in this paper is very ideal. It is clear that quickly approximates to $F\!\left( t \right)$ in Fig. 6.

Figure 5. The control deflection angles (in degrees).

Figure 6. The sliding surface.

We introduce the actuator fault scenario (AFS) as follow in order to further illustrate the performance of the designed controller.

(44a) \begin{equation}{\dot Y_1} = {H_1}\!\left( {{Y_1}} \right){Y_2},\end{equation}

(44b) \begin{equation}{\dot Y_2} = {H_2}\!\left( {{Y_2}} \right) + F\!\left( t \right) + {B_0}{{\Delta \hat U}}\left( t \right).\end{equation}

where ${{\Delta }} = {\textrm{diag}}\!\left( {{\delta _1},{\delta _2},{\delta _3}} \right)$ is the actuator effectiveness, which satisfies $0 \lt {\delta _i} \le 1,i = 1,2,3$ . Note that the case ${\delta _i}\left( t \right)$ means that the $i$ th actuator works normally, and $0 \lt {\delta _i} \lt 1$ corresponds the case in which the $i$ th actuator has partially lost its effectiveness, but it still works.

We define an AFS ${{\Delta }} = {\textrm{diag}}\left( {{\delta _1},{\delta _2},{\delta _3}} \right)$ as follow so as to further illustrate the performance of designed controller (42)

(45) \begin{equation}{\delta _i}\left( t \right) = \left\{ {\begin{array}{l@{\quad}l}{1,} & {{\textrm{if}}\ t \lt 5s}\\[6pt]{0.3 + 0.1{\textrm{sin}}\left( {0.5t + i\pi /3} \right),} & {{\textrm{if}}\ t \geqslant 5s}\end{array}} \right.\end{equation}

As shown in Fig. 7, the simulation experiment of the tracking errors in the AFS. It is set at 5s that the actuator fails, and there is a slight overshoot and chattering at 5s. The state response trajectory also fluctuates slightly, but it stabilises within 3s. It proves that the controller not only has good stability, but also has strong robustness and anti-interference.

Figure 7. The quaternion errors.

Compared with Fig. 4, the control signal under the AFS has a slight step at 5s and can be driven into the region of the SMS, which is shown in Fig. 8

1. (i) From the dynamic convergence process of the system, the robustness and stability of the proposed controller are obvious;

2. (ii) From the convergence rate and convergence trend of the state trajectory, the convergence speed is fast;

3. (iii) From Figs 7 and 8, even in the case of AFS, the controller still has strong stability and robustness. Although the state trajectory fluctuates slightly near zero at 3s, it becomes stable instantly.

Figure 8. The control deflection angles (in degrees).