Spacecraft attitude control based on generalised dynamic inversion with adaptive neural network

Abstract

This paper proposes a robust generalised dynamic inversion (GDI) control system design with adaptive neural network (NN) estimation for spacecraft attitude tracking under the absence of knowledge of the spacecraft inertia parameters. The robust GDI control system works to enforce attitude tracking, and the adaptive NN augmentation compensates for the lack of knowledge of the spacecraft inertia parameters. The baseline GDI control law consists of a particular part and an auxiliary part. The particular part of the GDI control law works to realise a desired attitude dynamics of the spacecraft, and the auxiliary part works for finite-time stabilisation of the spacecraft angular velocity. Robustness against modeling uncertainties and external disturbances is provided by augmenting a siding mode control element within the particular part of the GDI control law. The singularity that accompanies GDI control is avoided by modifying the Moore-Penrose generalised inverse by means of a dynamic scaling factor. The NN weighting matrices are updated adaptively through a control Lyapunov function. A detailed stability analysis shows that the closed loop system is semi-global practically stable. For performance assessment, a spacecraft model is developed, and GDI-NN control is investigated for its attitude control problem through numerical simulations. Simulation results reveal the efficacy, robustness and adaptive attributes of proposed GDI-NN control for its application to spacecraft attitude control.

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

1.0 Introduction

Various parts of a spacecraft are required to be pointing in certain directions for a successful accomplishment of the mission. The solar panels should be pointed towards the sun, the cameras should be pointed towards the targeted areas and antennas should be pointed towards the ground stations, etc. To achieve the desired pointing capabilities, highly agile attitude tracking is required that adjusts spacecraft orientation in the three-dimensional space.

Many practical spacecraft experience large variations in their inertia parameters during the course of the mission. For example a spacecraft that has to service another satellite in space will require getting hold of the target satellite. As a result, the inertia of the combined satellite will change significantly, requiring its control system to adapt accordingly to maintain the desired orientation [1]. LightSail 2 employs a partially deployed solar sail boom [2] and solar panels [3] that may cause uncertainties in its inertia parameters. Similarly the missions requiring a big change in velocity [4], also experience change in mass and hence inertia. Factors external to the spacecraft may also affect the dynamics of a spacecraft such as the solar nutation [5].

The subject of spacecraft stability and control has been studied extensively in recent decades. Several control methodologies have been proposed by researchers to address accurate spacecraft tracking, including robust control [6], back-stepping control [7], intelligent control [3], neural networks-based control [8–10], fault tolerant control [11, 12], sliding mode control [13–16], trajectory linearlisation control [17], optimal control [18] and adaptive control [19–21].

Nonlinear dynamic inversion (NDI) [22, 23] is a control design methodology that is used extensively in spacecraft control applications. The NDI control methodology relies on linearising the plant dynamics through feedback and applying linear control consequently. Nevertheless, NDI has several limitations, including the need of exact knowledge of system dynamics, square dimensionality constraints, the simplifying approximations made to obtain the inverse plant model and cancellation of some useful non-linearities. To overcome these challenges, the generalised dynamic inversion (GDI) control design methodology has been introduced in Ref. [24], and has been applied in solving several control problems for aerospace and robotic applications [25, 26].

GDI-based control works on the principle of defining prime objectives of control system through constraint differential equations that are asymptotically stable and then establishing the control law by inverting these dynamical equations employing dynamically scaled Moore-Penrose generalised inversion (MPGI) [25]. The architecture of this control methodology is composed of two components that are orthogonally complementary and are defined as the particular and auxiliary control components. Asymptotic stabilisation of the constraint dynamics is realised by the particular part, whereas auxiliary control offers additional control design freedom to meet secondary control objectives without interfering with the particular part.

In this work, a robust GDI control methodology aided with adaptive NN is implemented to address attitude control problem of spacecraft. Initially, GDI control law is constructed and contains a baseline GDI control augmented with sliding mode control (SMC) based switching control [27]. The switching control action is introduced to strengthen robustness against modeling uncertainties, non-linearities and external disturbances.

Additionally, radial basis function (RBF) NN is incorporated within the GDI framework to limit the dependencies of proposed control on spacecraft dynamics [28]. Lyapunov stability criteria have been employed to design adaptive laws for updating the weight matrices of the RBF-NN to estimate the spacecraft inertial parameters. Furthermore, comprehensive stability analysis of GDI-NN control is established using the concept of Lyapunov stability to ensure the attainment of reliable asymptotic tracking. It is noteworthy to mention here that the primary objective of incorporating RBF-NN with RGDI framework is not to ascertain the precise inertia parameters of the satellite, but rather to identify the parameters that facilitate the convergence of tracking error to zero.

The effectiveness of proposed GDI-NN control is demonstrated through computer simulations performed on spacecraft simulator developed in Simulink/Matlab environment.

The paper is organised in the following manner. The rotational kinematics and dynamics of a rigid body spacecraft are described in Section 2. The detailed formulation of GDI-NN control approach is presented in Section 3. Lyapunov based stability analysis of proposed control law is discussed in Section 4. Simulation results along with necessary discussions are presented in Section 5 followed by conclusion in Section 6.

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 [29].

(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 [29]

(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 [30]

(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. [30].

(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 [31]

(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 [32]

(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$ [33].

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 [31].

(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 [34]. 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 [35]. 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 [26].

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).

6.0 Conclusion

This paper presents the GDI-NN control system for attitude control of spacecraft. In proposed controller, constraint dynamics have been prescribed for stable attitude tracking, whereas switching control is introduced to ensure robustness. Furthermore RBF-NN are engaged for real-time estimation of inertial parameters which relaxes the requirement of these parameters for the controller. Stability analysis of proposed GDI-NN control is presented using Lyapunov stability criteria. Simulation results reveal that the attitude and angular body rate response curves rapidly converge towards the desired commands by generating the required torque commands proving the feasibility of proposed controller for its application to spacecraft attitude control.