Nomenclature

$\mathcal{I}$

earth centred inertial frame

$\mathcal{B}$

body fixed frame

${\textbf{q}}$

quaternion vector

$\omega $

angular velocity vector

${\textbf{J}}$

inertia matrix

$\tau $

control torque

${\textbf{e}}$

error vector

$\vartheta $

norm of attitude error

$\boldsymbol{\mathcal{A}}$

control coefficient vector

$B$

control load function

${\textbf{y}}$

null control vector

$\mathcal{P}$

null projection matrix

${\textbf{I}}$

identity matrix

$\nu $

dynamic scale factor

$\mathcal{Q}$

positive definite matrix

$\boldsymbol{\mathcal{L}}$

gain matrix

${\sigma _{max}}$

maximum singular value

$s$

sliding surface

$\boldsymbol{\mathcal{K}}$

sliding mode gain

${{\hat{\textbf{W}}}}$

approximated weight vector

${\textbf{h}}$

gaussian function

$\pi $

mean

$b$

standard deviation

$V$

lyapunov function

2.0 Spacecraft dynamics modeling

For establishing the dynamic model of the spacecraft, the rotational dynamics and kinematics of spacecraft are expressed using Euler’s equations. In spacecraft modeling, it is assumed that available thrust is continuous and the system parameters and applied disturbances are bounded. The spacecraft is assumed to be a rigid body and flexible bending modes have been neglected. Spacecraft kinematics equations are derived in Earth centred inertial frame $\mathcal{I}$ , whereas dynamics equations are derived in principal axis frame $\mathcal{B}$ .

Attitude kinematics of spacecraft in $\mathcal{B}$ with respect to $\mathcal{I}$ are described using quaternion as follows [Reference Wertz29].

(1) \begin{equation}{\textbf{q}} = {\left[ {\displaystyle\begin{array}{l} \displaystyle{{{\boldsymbol{\Upsilon }}^T}} \quad {{q_o}}\end{array}} \right]^T}\end{equation}

In above expression, ${\boldsymbol{\Upsilon }} = {[{q_1}{{\;\;\;\;}}{q_2}{{\;\;\;\;}}{q_3}]^T} \in {R^3}$ represents unit quaternion vector whereas the scalar part is represented by ${q_o} \in R$ . The constraint defined on the quaternion vector elements is given as ${{\textbf{q}}^T}{\textbf{q}} = 1$ . The rotational kinematics described in terms of differential equations are given as [Reference Wertz29]

(2) \begin{equation}{\boldsymbol{\dot{\mathrm{q}}}} = \left[ {\displaystyle\begin{array}{c}\displaystyle {{\boldsymbol{\dot{\Upsilon}}}}\\ \displaystyle {{{\dot q}_o}}\end{array}} \right] = \frac{1}{2}\left[ {\begin{array}{{c}}{\displaystyle\left( {{{\boldsymbol{\Upsilon }}^ \times } + {q_o}{{\textbf{I}}_{3 \times 3}}} \right)}\\ \displaystyle{ - {{\boldsymbol{\Upsilon }}^T}}\end{array}} \right]\omega \;\;\;\;{\textbf{q}}(0) = {{\textbf{q}}_i}\end{equation}

where $\omega = {[p{{\;\;\;\;}}q{{\;\;\;\;}}r]^T} \in {R^{3 \times 1}}$ denotes angular velocity vector expressed in frame $\mathcal{B}$ , vector ${{\textbf{I}}_{3 \times 3}}$ represents identity matrix and the term ${{\boldsymbol{\Upsilon }}^ \times }$ represents the skew-symmetric cross product matrix described as

(3) \begin{equation}{{\boldsymbol{\Upsilon }}^ \times } = \left[ {\begin{array}{ccc} \displaystyle0 &{ - {q_3}} &{{q_2}}\\[4pt] \displaystyle{{q_3}} &0 &{ - {q_1}}\\[4pt] \displaystyle { - {q_2}} &{{q_1}} &0\end{array}} \right]\end{equation}

Similarly, the dynamics of spacecraft is represented by Euler’s model, which implies

(4) \begin{equation}{\boldsymbol{\dot{\omega}}} = - {{\textbf{J}}^{ - 1}}{\omega ^ \times }{\textbf{J}}\omega + {{\textbf{J}}^{ - 1}}\tau ,{{\;\;\;\;}}\omega (0) = {\omega _0}\end{equation}

where ${\textbf{J}} = diag\left[ {{j_x}{{\;\;\;\;}}{j_y}{{\;\;\;\;}}{j_z}} \right] \in {R^{3 \times 3}}$ is the inertia matrix and vector $\tau = {[{\tau _x}{{\;\;\;\;}}{\tau _y}{{\;\;\;\;}}{\tau _z}]^T} \in {R^{3 \times 1}}$ represents the control torque generated around the body axis of spacecraft.

3.0 Control design formulation

In this section, formulation of GDI-NN control architecture is illustrated to achieve the pointing accuracy in spacecraft. To formulate the asymptotically stable attitude constraint dynamics, quaternion error function in the form of quaternion products is expressed as [Reference Show and Juang30]

(5) \begin{equation}{{\textbf{e}}_q}\!\left( {{\textbf{q}},t} \right) = {{\textbf{q}}^{ - 1}} \otimes {{\textbf{q}}_d}(t),\forall {{\;\;\;\;}}t \in \!\left( {0,\infty } \right)\end{equation}

where ${{\textbf{q}}_d}(t) = {[{{\boldsymbol{\Upsilon }}_d}(t){{\;\;\;\;}}{q_{od}}(t)]^T}$ is the desired quaternion unit vector and ${{\textbf{q}}^{ - 1}}$ represents quaternion conjugate described as

(6) \begin{equation}{{\textbf{q}}^{ - 1}} = \left[ {\begin{array}{ll}\displaystyle{ - {{\boldsymbol{\Upsilon }}^T}} & \displaystyle {{q_o}}\end{array}} \right]\end{equation}

The product of quaternion established in Equation (5) can be written as, see Ref. [Reference Show and Juang30].

(7) \begin{equation}{{\textbf{e}}_q} = \left[ {\begin{array}{c}\displaystyle{{{\boldsymbol{\Upsilon }}_v}}\\ \displaystyle{{q_{eo}}}\end{array}} \right] = \left[ {\begin{array}{c}\displaystyle{{q_{do}}(t){\boldsymbol{\Upsilon }} - {q_o}{{\boldsymbol{\Upsilon }}_d}(t) - {\boldsymbol{\Upsilon }}_d^ \times (t){\boldsymbol{\Upsilon }}}\\ \displaystyle{{q_o}{q_{do}}(t) + {{\boldsymbol{\Upsilon }}^T}{{\boldsymbol{\Upsilon }}_d}(t)}\end{array}} \right]\end{equation}

The quaternion error vector given by Equation (7) must satisfy the constraint such that ${\textbf{e}}_q^T{{\textbf{e}}_q} = \!\left( {{{\boldsymbol{\Upsilon }}^T}{\boldsymbol{\Upsilon }} + q_o^2} \right)\{ {\boldsymbol{\Upsilon }}_d^T(t){{\boldsymbol{\Upsilon }}_d}(t) + {q_{do}}{(t)^2}\} = 1$ . Hence, the quaternion error kinematics in the form of differential equations are given as

(8) \begin{equation}{{\boldsymbol{\dot{\mathrm{e}}}}_q} = \left[ {\begin{array}{c}\displaystyle{{{{\boldsymbol{\dot{\Upsilon}}}}_e}}\\ \displaystyle{{{\dot q}_{eo}}}\end{array}} \right] = \frac{1}{2}\left[ {\begin{array}{c}\displaystyle{{\boldsymbol{\Upsilon }}_e^ \times + {q_{eo}}{{\textbf{I}}_{3 \times 3}}}\\ \displaystyle{ - {\boldsymbol{\Upsilon }}_e^T}\end{array}} \right]{{\textbf{e}}_\omega }\end{equation}

where ${{\textbf{e}}_\omega }$ represents the error vector of angular velocity that can be described in differential form as follows

(9) \begin{align} {{\boldsymbol{\dot{\mathrm{e}}}}_\omega } &= {\boldsymbol{\dot{\omega}}} - {{\boldsymbol{\dot{\omega}}} _d}(t)\nonumber\\&= - {{\textbf{J}}^{ - 1}}{\omega ^ \times }{\textbf{J}}\omega + {{\textbf{J}}^{ - 1}}\omega _d^ \times (t){\textbf{J}}{\omega _d}(t) + {{\textbf{J}}^{ - 1}}\tau \end{align}

where ${\omega _d}(t):\!\left( {0,\infty } \right) \to {R^3}$ represents the desired vector of angular velocity, whose dynamics are governed by the following equation

(10) \begin{equation}{{\boldsymbol{\dot{\omega}}} _d}(t) = - {{\textbf{J}}^{ - 1}}\omega _d^ \times (t){\textbf{J}}{\omega _d}(t)\end{equation}

3.1 Constraint formulation

Let $\vartheta :\left[ { - 1,1\left] \to \right[0,1} \right]$ be the norm of attitude error function, described as

(11) \begin{equation}\vartheta\!\left( {{q_{eo}}} \right) = 1 - q_{eo}^2\end{equation}

The time derivatives of $\vartheta $ are computed as [Reference Bajodah31]

(12) \begin{equation}\dot \vartheta = {q_{eo}}{\boldsymbol{\Upsilon }}_e^T{{\textbf{e}}_\omega }\end{equation}

(13) \begin{equation}\ddot \vartheta = \frac{1}{2}{\textbf{e}}_\omega ^T\!\left( {q_{eo}^2{{\textbf{I}}_{3 \times 3}} - {{\boldsymbol{\Upsilon }}_e}{\boldsymbol{\Upsilon }}_e^T} \right){{\textbf{e}}_\omega } + {q_{eo}}{\boldsymbol{\Upsilon }}_e^T{{\textbf{J}}^{ - 1}}\boldsymbol{\tau} \end{equation}

Based on the time derivatives of $\vartheta $ , spacecraft constraint dynamics can be represented through following differential equation

(14) \begin{equation}\ddot \vartheta + {c_1}\dot \vartheta + {c_2}\vartheta = 0\end{equation}

Here ${c_1}$ and ${c_2}$ are positive gain constants. The constraint differential equation given by Equation (14) is made asymptotically stable through suitable selection of the gain constants ${c_1}$ and ${c_2}$ . By obtaining the time derivatives of $\vartheta $ in Equation (14), following constraint differential equation in algebraic form is obtained

(15) \begin{equation}\boldsymbol{\mathcal{A}}\!\left( {{{\textbf{e}}_q}} \right)\tau = B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right)\end{equation}

where $\boldsymbol{\mathcal{A}}:{{\textbf{R}}^n} \times \left[ {{t_0},\infty } \right) \to {{\textbf{R}}^{1 \times 3}}$ is the control coefficient vector function described as

(16) \begin{equation}\boldsymbol{\mathcal{A}}\!\left( {{{\textbf{e}}_q}} \right) = {q_{eo}}{\boldsymbol{\Upsilon }}_e^T{{\textbf{J}}^{ - 1}}\end{equation}

and $B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) \in {\textbf{R}}$ symbolise the control load function prescribed as

(17) \begin{equation}B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) = - \frac{1}{2}{\textbf{e}}_\omega ^T\!\left( {q_{eo}^2{{\textbf{I}}_{3 \times 3}} - {{\boldsymbol{\Upsilon }}_e}{\boldsymbol{\Upsilon }}_e^T} \right){{\textbf{e}}_\omega } - {c_1}{q_{eo}}{\boldsymbol{\Upsilon }}_e^T{{\textbf{e}}_\omega } - {c_2}\!\left( {1 - q_{eo}^2} \right)\end{equation}

3.2 Design of robust GDI control

The objective of GDI control is to guarantee asymptotic realisation of prescribed constraint dynamics for stable attitude tracking such that $\vartheta \!\left( {{q_{eo}}} \right) = 0$ or ${{\boldsymbol{\Upsilon }}_e} = {0_{3 \times 1}}$ . Control law is obtained through inversion of prescribed constraint dynamics given by Equation (15) employing dynamically scaled MPGI, implies

(18) \begin{equation}\boldsymbol{\tau ^{\textrm{*}}} = \boldsymbol{\boldsymbol{\mathcal{A}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) + P\!\left( {{{\textbf{e}}_q}} \right){\textbf{y}}\end{equation}

where ${\textbf{y}}$ represents null control vector and $\mathcal{P}\!\left( {{{\textbf{e}}_q}} \right)$ represents null projection matrix described as

(19) \begin{equation}\boldsymbol{\mathcal{P}}\!\left( {{{\textbf{e}}_q}} \right) = {{\textbf{I}}_{3 \times 3}} - \boldsymbol{\boldsymbol{\mathcal{A}}^ + }\!\left( {{{\textbf{e}}_q}} \right)\boldsymbol{\mathcal{A}}\!\left( {{{\textbf{e}}_q}} \right)\end{equation}

The term ${\boldsymbol{\mathcal{A}}^{\textrm{*}}}$ in the control law defined in Equation (18) indicates dynamically scaled MPGI and is given by

(20) \begin{equation}\boldsymbol{\boldsymbol{\mathcal{A}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right) = \boldsymbol{\boldsymbol{\mathcal{A}}^T}\!\left( {{{\textbf{e}}_q}} \right){\{ \boldsymbol{\mathcal{A}}\!\left( {{{\textbf{e}}_q}} \right)\boldsymbol{\boldsymbol{\mathcal{A}}^T}\!\left( {{{\textbf{e}}_q}} \right) + \nu (t)\} ^{ - 1}}\end{equation}

where $\nu $ represents dynamic scale factor, its differential form is describes as [Reference Ansari and Bajodah32]

(21) \begin{equation}\dot \nu (t) = - \nu (t) + \frac{\gamma }{{\|{{\textbf{e}}_\omega }^2\|}},{\textrm{\;\;}}\nu (0) \gt 0\end{equation}

here $\gamma $ is a positive constant and $\parallel .{\parallel ^2}$ is the Euclidean norm. Hence elements of $\boldsymbol{\boldsymbol{\mathcal{A}}^{\textrm{*}}}$ remain bounded for all $t \gt 0$ [Reference Ansari, Bajodah and Kada33].

Similarly closed-loop stability of angular body rates is achieved by utilising null control vector ${\textbf{y}} = \boldsymbol{\mathcal{L}}{{\textbf{e}}_\omega }$ , its gain matrix $\boldsymbol{\mathcal{L}} \in {R^{3 \times 3}}$ is established by Lyapunov principle [Reference Bajodah31].

(22) \begin{equation}\boldsymbol{\mathcal{L}} = - \dot{\mathcal{P}}\!\left( {{{\textbf{e}}_q}} \right){{\textbf{e}}_\omega }) - {\sigma _{max}}\!\left( {\dot{\mathcal{P}}\!\left( {{{\textbf{e}}_q}} \right){{\textbf{e}}_\omega }} \right)) - \mathcal{Q}\end{equation}

here $\mathcal{Q} \in {R^{3 \times 3}}$ is a positive definite matrix and variable ${\sigma _{max}}$ symbolise the maximum singular value. By incorporating the null control vector, the baseline control law established in Equation (18) is described as

(23) \begin{equation}\boldsymbol{\tau ^{\textrm{*}}} = \boldsymbol{\boldsymbol{\mathcal{A}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) + P\!\left( {{{\textbf{e}}_q}} \right)\!\boldsymbol{\mathcal{L}}{{\textbf{e}}_\omega }\end{equation}

By placing the control expression given by Equation (23) in Equation (9), the spacecraft dynamics are expressed as

(24) \begin{equation}{\boldsymbol{\dot{\omega}}} = - {{\textbf{J}}^{ - 1}}{\omega ^ \times }{\textbf{J}}\omega + {{\textbf{J}}^{ - 1}}\left\{ {\boldsymbol{\boldsymbol{\mathcal{A}}^ + }\!\left( {{{\textbf{e}}_q}} \right)B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) + \mathcal{P}\!\left( {{{\textbf{e}}_q}} \right)\boldsymbol{\mathcal{L}}{{\textbf{e}}_\omega }} \right\}\end{equation}

Additionally GDI control given by Equation (23) is augmented with a switching term, which is synthesised to incorporate robustness against variations in system parameters and nonlinearities as well as to provide improved tracking performance. The resultant robust GDI control expression is written as

(25) \begin{equation}\boldsymbol{\tau ^{\textrm{*}}} = \boldsymbol{\boldsymbol{\mathcal{A}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) + {\textbf{P}}\!\left( {{{\textbf{e}}_q}} \right)\boldsymbol{\boldsymbol{\mathcal{L}}}{{\textbf{e}}_\omega } - \boldsymbol{\boldsymbol{\mathcal{K}}}\boldsymbol{\boldsymbol{\mathcal{A}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)\frac{s}{{\parallel s\parallel }}\end{equation}

where $\boldsymbol{\mathcal{K}}$ denotes sliding mode gain constant and $s$ is the sliding surface defined as

(26) \begin{equation}s = \dot \vartheta + {c_1}\vartheta + {c_2}\int \vartheta dt\end{equation}

3.3 Real-time parameter estimation

For real-time estimation of inertia matrix ${\textbf{J}}$ in control coefficient vector $\boldsymbol{\mathcal{A}}$ given by Equation (16), RBF-NNs have been engaged and is described as

(27) \begin{equation}\hat f = {{{\hat{\textbf{W}}}}^T}{\textbf{h}}\!\left( {{{\textbf{e}}_q}} \right)\end{equation}

where ${{\hat{\textbf{W}}}} \in {R^{1 \times 5}}$ is the approximated weight vector and ${\textbf{h}} \in {R^{1 \times 5}}$ is the Gaussian function described as

(28) \begin{equation}{h_j}\!\left( {\textbf{u}} \right) = exp\left\{ {\frac{{\parallel\!{\textbf{u}} - {\pi _{ij}}{\parallel ^2}}}{{2b_j^2}}} \right\}\end{equation}

where number of inputs and hidden layers are denoted by subscript $i$ and $j$ respectively, $\pi $ is the mean and $b$ is the standard deviation of Gaussian function [Reference Ansari and Bajodah34]. Based on this, inertial parameters of the spacecraft assuming axial symmetry are approximated by

(29) \begin{equation}{\textbf{J}} = diag\left\{\, {{{\hat f}_{{j_x}}}{{\;\;\;\;}}{{\hat f}_{{j_y}}}{{\;\;\;\;}}{{\hat f}_{{j_z}}}} \right\}\end{equation}

where

(30) \begin{equation}{\hat f_{{j_x}}}\!\left( {{{\textbf{e}}_q}} \right) = {{{\hat{\textbf{W}}}}_{\textbf{x}}}^T{{\textbf{h}}_{\textbf{x}}}\!\left( {{{\textbf{e}}_q}} \right){{\;\;\;\;}}{\hat f_{{j_y}}}\!\left( {{{\textbf{e}}_q}} \right) = {{{\hat{\textbf{W}}}}_{\textbf{y}}}^T{{\textbf{h}}_{\textbf{y}}}\!\left( {{{\textbf{e}}_q}} \right){{\;\;\;\;}}{\hat f_{{j_z}}}\!\left( {{{\textbf{e}}_q}} \right) = {f_{{j_x}}}\end{equation}

In Equation (30), weight vectors ${{{\hat{\textbf{W}}}}_{\textbf{x}}}$ and ${{{\hat{\textbf{W}}}}_{\textbf{y}}}$ are updated in real-time using Lyapunov approach, whose derivation is explained in the following section. Since the control coefficient vector $\boldsymbol{\mathcal{A}}$ given by Equation (16) depends on the inertia matrix ${\textbf{J}}$ which is estimated using RBF-NN, therefore the resultant control expression given by Equation (25) is re-defined as

(31) \begin{equation}\boldsymbol{\tau ^{\textrm{*}}} = \boldsymbol{\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)\textbf{B}\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) + {\textbf{P}}\!\left( {{{\textbf{e}}_q}} \right)\boldsymbol{\boldsymbol{\mathcal{L}}}{{\textbf{e}}_\omega } - \boldsymbol{\boldsymbol{\mathcal{K}}}\boldsymbol{\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)\frac{s}{{\parallel s\parallel }}\end{equation}

where ${\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}$ contains the elements of the inertia matrix ${\textbf{J}}$ which are estimated using RBF-NN. The overall block diagram of GDI-NN control loop is depicted in Fig. 1.

Figure 1. Block diagram of proposed control architecture.

4.0 Stability analysis

In this section, the proof of asymptotic convergence of attitude error and angular body rates of spacecraft has been established through following theorem.

Theorem 1. By considering the closed loop system defined in Equation (24), GDI-NN control expression given by Equation (31) will guarantee that the controlled system is globally uniformly bounded and the error function of body attitudes and angular velocities are semi-global practically stable.

Proof. Computation of the derivative of sliding surface $s$ is given as

(32) \begin{equation}\dot s = \ddot \vartheta + {c_1}\dot \vartheta + {c_2}\vartheta \end{equation}

Asymptotic convergence of $s$ towards zero infers asymptotic realisation of the constraint dynamics established in Equation (15), which implies

(33) \begin{equation}\dot s = \boldsymbol{\mathcal{A}}\!\left( {{{\textbf{e}}_q}} \right)\boldsymbol{\tau^{\!\textrm{*}}} - B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right)\end{equation}

By updating the control law $\boldsymbol{\tau ^{\textrm{*}}}$ established in Equation (31) and replacing ${\textbf{P}}$ with ${{\textbf{I}}_{3 \times 3}} - \boldsymbol{\boldsymbol{\mathcal{A}}^ + }\!\left( {{{\textbf{e}}_q}} \right)\boldsymbol{\mathcal{A}}\!\left( {{{\textbf{e}}_q}} \right)$ , sliding mode dynamics given by Equation (33) can be written as

(34) \begin{align}\dot s &= \boldsymbol{\mathcal{A}}\!\left( {{{\textbf{e}}_q}} \right)\{ {\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) + \left\{ {{{\textbf{I}}_{3 \times 3}} - \boldsymbol{\boldsymbol{\mathcal{A}}^ + }\!\left( {{{\textbf{e}}_q}} \right)\boldsymbol{\mathcal{A}}\!\left( {{{\textbf{e}}_q}} \right)} \right\}\nonumber\\ &\qquad\qquad \qquad \qquad\boldsymbol{\mathcal{L}}{{\textbf{e}}_\omega } - \boldsymbol{\mathcal{K}}{\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)\frac{s}{{\parallel s\parallel }}\} - B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right)\end{align}

With respect to the property of pseudo inverse i.e., $\boldsymbol{\mathcal{A}}\!\left( {{{\textbf{e}}_q}} \right)\boldsymbol{\boldsymbol{\mathcal{A}}^ + }\!\left( {{{\textbf{e}}_q}} \right) = 1$ for all $\boldsymbol{\mathcal{A}}\!\left( {{{\textbf{e}}_q}} \right) \ne {0_{1 \times 3}}$ , and replacing $\boldsymbol{\mathcal{A}}$ with $\skew{10}\hat{\boldsymbol{\mathcal{A}}} + \skew{10}\tilde{\boldsymbol{\mathcal{A}}}$ , sliding mode dynamics described in Equation (34) is written in simplified form as

(35) \begin{align}\dot s &= \left\{ {\skew{10}\hat{\boldsymbol{\mathcal{A}}}\!\left( {{{\textbf{e}}_q}} \right) + \skew{10}\tilde{\boldsymbol{\mathcal{A}}}\!\left( {{{\textbf{e}}_q}} \right)} \right\}\left\{ {{{\skew{10}\hat{\boldsymbol{\mathcal{A}}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) - \boldsymbol{\mathcal{K}}{{\skew{10}\hat{\boldsymbol{\mathcal{A}}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)\frac{s}{{\parallel s\parallel }}} \right\}\nonumber\\ &\qquad - B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right)\nonumber\\ &= \skew{10}\hat{\boldsymbol{\mathcal{A}}}\!\left( {{{\textbf{e}}_q}} \right){\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) - B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) - \boldsymbol{\mathcal{K}}\skew{10}\hat{\boldsymbol{\mathcal{A}}}\!\left( {{{\textbf{e}}_q}} \right){\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)\nonumber\\&\qquad\frac{s}{{\parallel s\parallel }} + \skew{10}\tilde{\boldsymbol{\mathcal{A}}}\!\left( {{{\textbf{e}}_q}} \right){\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) - \boldsymbol{\mathcal{K}}\skew{10}\tilde{\boldsymbol{\mathcal{A}}}\!\left( {{{\textbf{e}}_q}} \right){\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)\frac{s}{{\parallel s\parallel }}\nonumber\\&= \left\{ {\rho - 1} \right\}B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) - \boldsymbol{\mathcal{K}}\skew{10}\hat{\boldsymbol{\mathcal{A}}}\!\left( {{{\textbf{e}}_q}} \right){\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)\frac{s}{{\parallel s\parallel }}\nonumber\\&\qquad+ \skew{10}\tilde{\boldsymbol{\mathcal{A}}}\!\left( {{{\textbf{e}}_q}} \right){\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)B\!\left( {{{\textbf{e}}_q},{{\textbf{e}}_\omega }} \right) - \boldsymbol{\mathcal{K}}\skew{10}\tilde{\boldsymbol{\mathcal{A}}}\!\left( {{{\textbf{e}}_q}} \right){\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)\frac{s}{{\parallel s\parallel }}\end{align}

where $\rho = \skew{10}\hat{\boldsymbol{\mathcal{A}}}\!\left( {{{\textbf{e}}_q}} \right){\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}\!\left( {{{\textbf{e}}_q},\nu } \right)$ . Defining a Lyapunov function $V$ that is positive definite

(36) \begin{equation}V = \frac{1}{2}{s^2} + \frac{1}{2}{\gamma _1}{{{\tilde{\textbf{W}}}}_{\textbf{x}}}{{{\tilde{\textbf{W}}}}_{\textbf{x}}}^T + \frac{1}{2}{\gamma _2}{{{\tilde{\textbf{W}}}}_{\textbf{y}}}{{{\tilde{\textbf{W}}}}_{\textbf{y}}}^T\end{equation}

By differentiating Lyapunov function $V$ with respect to time, we get the following

(37) \begin{equation}\dot V = s\dot s + {\gamma _1}{{{\tilde{\textbf{W}}}}_{\textbf{x}}}{\dot{\hat{\textbf{W}}}}_{\textbf{x}}^{\textbf{T}} + {\gamma _2}{{{\tilde{\textbf{W}}}}_{\textbf{y}}}{\dot{\hat{\textbf{W}}}}_{\textbf{y}}^{\textbf{T}}\end{equation}

Substituting the expression of $\dot s$ given by Equation (35) in Equation (37), implies

(38) \begin{align} \dot V &= \left\{ {\rho - 1} \right\}Bs - \boldsymbol{\mathcal{K}}\rho \parallel s\parallel + \skew{10}\tilde{\boldsymbol{\mathcal{A}}}\left\{ {{{\skew{10}\hat{\boldsymbol{\mathcal{A}}}}^{\textrm{*}}}B - \boldsymbol{\mathcal{K}}{{\skew{10}\hat{\boldsymbol{\mathcal{A}}}}^{\textrm{*}}}\parallel s\parallel } \right\}\nonumber\\ &\qquad+ {\gamma _1}{{{\tilde{\textbf{W}}}}_{\textbf{x}}}{\dot{\hat{\textbf{W}}}}_{\textbf{x}}^{\textbf{T}} + {\gamma _2}{{{\tilde{\textbf{W}}}}_{\textbf{y}}}{\dot{\hat{\textbf{W}}}}_{\textbf{y}}^{\textbf{T}}\end{align}

Let $\delta = {\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}B - \boldsymbol{\mathcal{K}}{\skew{10}\hat{\boldsymbol{\mathcal{A}}}^{\textrm{*}}}\parallel s\parallel $ and replace $\skew{10}\hat{\boldsymbol{\mathcal{A}}}$ with ${{{\tilde{\textbf{W}}}}_{\textbf{x}}}^T{{\textbf{h}}_{\textbf{x}}}\!\left( {{{\textbf{e}}_q}} \right) + {{{\tilde{\textbf{W}}}}_{\textbf{y}}}^T{{\textbf{h}}_{\textbf{y}}}\!\left( {{{\textbf{e}}_q}} \right)$ , then expression of $\dot V$ implies

(39) \begin{align}\dot V &= \left\{ {\rho - 1} \right\}Bs - \boldsymbol{\mathcal{K}}\rho \parallel s\parallel + {{{\tilde{\textbf{W}}}}_{\textbf{x}}}^T{{\textbf{h}}_{\textbf{x}}}\!\left( {{{\textbf{e}}_q}} \right)\delta + {{{\tilde{\textbf{W}}}}_{\textbf{y}}}^T{{\textbf{h}}_{\textbf{y}}}\!\left( {{{\textbf{e}}_q}} \right)\delta \nonumber\\ &\qquad+ {\gamma _1}{{{\tilde{\textbf{W}}}}_{\textbf{x}}}{\dot{\hat{\textbf{W}}}}_{\textbf{x}}^{\textbf{T}} + {\gamma _2}{{{\tilde{\textbf{W}}}}_{\textbf{y}}}{\dot{\hat{\textbf{W}}}}_{\textbf{y}}^{\textbf{T}}\nonumber\\ &= \left\{ {\rho - 1} \right\}Bs - \boldsymbol{\mathcal{K}}\rho \parallel s\parallel + {{{\tilde{\textbf{W}}}}_{\textbf{x}}}^T\left\{ {{{\textbf{h}}_{\textbf{x}}}\!\left( {{{\textbf{e}}_q}} \right)\delta + {\gamma _1}{\dot{\hat{\textbf{W}}}}_{\textbf{x}}^{\textbf{T}}} \right\}\nonumber\\ &\qquad + {{{\tilde{\textbf{W}}}}_{\textbf{y}}}^T\left\{ {{{\textbf{h}}_{\textbf{y}}}\!\left( {{{\textbf{e}}_q}} \right)\delta + {\gamma _2}{\dot{\hat{\textbf{W}}}}_{\textbf{y}}^{\textbf{T}}} \right\}\end{align}

The updating principle of the weight vectors are computed as

(40) \begin{equation}\dot{\hat{\textbf{W}}}_{\textbf{x}}^{\textbf{T}} = - \frac{1}{{{\gamma _1}}}{{\textbf{h}}_{\textbf{x}}}\!\left( {{{\textbf{e}}_q}} \right)\delta ,{{\;\;\;\;}}{\dot{\hat{\textbf{W}}}}_{\textbf{y}}^{\textbf{T}} = - \frac{1}{{{\gamma _2}}}{{\textbf{h}}_{\textbf{y}}}\!\left( {{{\textbf{e}}_q}} \right)\delta \end{equation}

By placing weight vectors in Equation (39), the simplified expression of $\dot V$ is written as

(41) \begin{equation}\dot V = \left\{ {\rho - 1} \right\}Bs - \boldsymbol{\mathcal{K}}\rho \parallel s\parallel \end{equation}

Since $\nu \in \!\left( {0,\infty } \right)$ , therefore $0 \lt \rho \!\left( {{{\textbf{e}}_q},\nu } \right) \lt 1$ for all $\boldsymbol{\mathcal{A}}\!\left( {{{\textbf{e}}_q}} \right) \ne {0_{1 \times 3}}$ , which follows that

(42) \begin{equation}\mathop {{\textrm{lim}}}\limits_{t \to \infty } \rho \!\left( {{{\textbf{e}}_q},\nu } \right) = 0{\textrm{}} \Leftrightarrow {\textrm{}}\mathop {{\textrm{lim}}}\limits_{t \to \infty } \boldsymbol{\mathcal{A}}\!\left( {{{\textbf{e}}_q}} \right) = {0_{1 \times 3}}\end{equation}

Therefore, a sliding mode gain $\boldsymbol{\mathcal{K}}$ that ensures

(43) \begin{equation}\boldsymbol{\mathcal{K}} \gt \frac{{\rho - 1}}{\rho }\left| B \right|\end{equation}

will also make sure that $\dot V$ is negative definite and consequently $s = 0$ , as established in Lyapunov’s direct method [Reference Khalil35]. By achieving the condition $s = 0$ , $\nu \in \!\left( {0,\infty } \right)$ also converges towards zero, which implies

(44) \begin{equation}\mathop {{\textrm{lim}}}\limits_{{{\textbf{e}}_q} \to 0} \frac{{\rho \!\left( {{{\textbf{e}}_q},\nu } \right) - 1}}{{\rho \!\left( {{{\textbf{e}}_q},\nu } \right)}} = - \infty \end{equation}

This implies that an infinite value of sliding mode gain $\boldsymbol{\mathcal{K}}$ is required to achieve negative definiteness of $\dot V$ which is not practical. Therefore to address weaker stability condition, sliding mode gain is designed to satisfy semi-global practical stability condition [Reference Ansari, Bajodah and Hamayun26].

Figure 2. Quaternion vs. time (stabilisation).

Figure 3. Body rates vs. time (stabilisation).

5.0 Numerical simulations and discussion

Performance of GDI-NN control for spacecraft attitude control is analysed and discussed through numerical simulations. As indicated in Section 2, the simulations model is developed with the assumptions that the spacecraft is a rigid body, the control action is continuous, disturbances are bounded and co-ordinate system of spacecraft coincides with the orbital co-ordinate system. For numerical simulations, the initial values of inertia matrix is assumed to be ${\textbf{J}} = diag\!\left( {8,8,8} \right)$ , whereas width and mean values of the RBF-NN functions are initialised as

(45) \begin{equation}{b_j} = {[1.5{{\;\;\;\;}}1.5{{\;\;\;\;}}1.5{{\;\;\;\;}}1.5{{\;\;\;\;}}1.5]^T}\end{equation}

(46) \begin{equation}{\pi _{ij}} = \begin{array}{c@{\,\,\,}c@{\,\,\,}c@{\,\,\,}c@{\,\,\,}c}0 &{25} &{50} &{75} &{100}\\0 &{25} &{50} &{75} &{100}\\0 &{25} &{50} &{75} &{100}\\0 &{25} &{50} &{75} &{100}\end{array}\end{equation}

In first scenario, stabilisation effect with GDI-NN is analysed by re-orienting the spacecraft from initial attitude ${{\textbf{q}}_0} = {[0.3{{\;\;\;\;}}0.4{{\;\;\;\;}}0.2{{\;\;\;\;}}0.85]^T}$ expressed in the form of quaternion to a desired quaternion, i.e. ${{\textbf{q}}_d} = {[0.7{{\;\;\;\;}}0.5{{\;\;\;\;}} - 0.4{{\;\;\;\;}}0.35]^T}$ . For this manoeuver, the initial body rates ${\omega _0}$ are set to be $10$ deg/sec and the simulation is executed for $5$ seconds. The actual and the commanded attitude profiles in terms of quaternion are shown in Fig. 2. Similarly the angular body rates are depicted in Fig. 3. It is apparent from results that the proposed controller smoothly stabilises the spacecraft and enables it to make a transition from its initial orientation to the desired orientation within $1.5$ seconds. The generated input torque commands ${\tau _x},{\tau _y},{\tau _z}$ are shown in Fig. 4, which are well within the specified saturation limits of $10$ Nm.

Figure 4. Control torque vs. time (stabilisation).

Figure 5. Quaternion vs. time (tracking).

Figure 6. Body rates vs. time (tracking).

The controller is further investigated by commanding the following body rate profile

(47) \begin{equation}{\omega _d} = {[10{\textrm{sin}}\!\left( {0.1t} \right){{\;\;\;\;}}15{\textrm{cos}}\!\left( {0.2t} \right){{\;\;\;\;}}20{\textrm{cos}}\!\left( {0.3t} \right)]^T},{\textrm{deg}}/{\textrm{sec}}\end{equation}

In this case, it is assumed that the spacecraft has an initial angular speed of ${\omega _0} = {[10{{\;\;\;\;}}20{{\;\;\;\;}}30]^T}$ deg/sec. Spacecraft attitude and angular body rates plots are depicted in Figs 5 and 6, respectively. These plots show that the designed control system enables the spacecraft to seamlessly follow the desired body rate profile from its initial angular speed. In this case, the simulation is executed for $100$ seconds. The time histories of the corresponding control torques commands required for subject manoeuvering are shown in Fig. 7. The figure shows that the control effort required to achieve the desired manoeuver is well within the saturation limit of $10$ Nm and once the spacecraft starts tracking the desired body rate profile, the control effort is minimal.

Figure 7. Control torque vs. time (tracking).