A Gimballed Control Moment Gyroscope Cluster Design for Spacecraft Attitude Control

Abstract

This paper addresses the problem of singularity avoidance in a cluster of four Single-Gimbal Control Moment Gyroscopes (SGCMGs) in a pyramid configuration when used for the attitude control of a satellite by introducing a new gimballed control moment gyroscope (GCMG) cluster scheme. Four SGCMGs were used in a pyramid configuration, along with an additional small and simple stepper motor that was used to gimbal the full cluster around its vertical (z) axis. Contrary to the use of four variable-speed control moment gyroscopes (VSCMGs), where eight degrees of freedom are available for singularity avoidance, the proposed GCMG design uses only five degrees of freedom (DoFs), and a modified steering law was designed for the new setup. The proposed design offers the advantages of SGCMGs, such as a low weight, size, and reduced complexity, with the additional benefit of overcoming the internal elliptic singularities, which create a minor attitude error. A comparison with the four-VSCMG cluster was conducted through numerical simulations, and the results indicated that the GCMG design was considerably more efficient in terms of power while achieving a better gimbal configuration at the end of the simulation, which is essential when it is desired for different manoeuvres to be consecutively executed. Additionally, for a nano-satellite of a few kilograms, the results prove that it is feasible to manufacture the GCMG concept by using affordable and lightweight commercial off-the-shelf (COTS) stepper motors.

1. Introduction

Single-gimbal control moment gyroscopes (SGCMGs) have become popular actuators for satellites that have an agility requirement [1,2]. They are widely used in spacecraft control for manoeuvering and have been successfully employed for a wide range of space missions [3,4,5]. Agile satellites are highly manoeuvrable, and they enable more data to be gathered in a single pass compared to those employed with reaction/momentum wheel actuators. However, SGCMGs come with the inherent problem of singularities. Singularities occur when the system is unable to generate a torque output along a specific direction [6]. There are two types of singularities with respect to the ability to escape from the singular state through null motion [2]. Elliptic singularities cannot be escaped through null motion, whilst there is a null motion path for escaping hyperbolic singularities. However, in cases where the angular speed of the flywheel is equal for all SGCMGs in the cluster, it is not possible to re-orient a symmetric gimbal set to a non-symmetric set through null-space motion in order to get a singularity-free path. Such an action would result in a momentum change, and attitude error would be added into the system.

Several studies have focused on finding a measure of distance from singularities, and different indices have been proposed in the literature [7], with the most commonly used being the manipulability index [8]. Multiple SGCMG singularity avoidance/escape logics and methods have been proposed to date [5]. The singularity-robust inverse is the most common, but other methods, including null motion and path planning, have also been proposed [9,10,11]. In general, the method of planning the path of a manoeuvre also includes the optimisation of the initial configuration of the gimbals before the execution of the desired manoeuvre in order to enhance the performance of the system and provide a singularity-free trajectory [12,13,14]. A different path-planning technique that blends the pseudo-spectral and direct-shooting methods upon gimbal saturation and singularity constraints was presented in [15]. Energy-consumption-based path-planning techniques for double-gimbal CMGs have also been studied [16]. A global singularity-avoidance steering law was presented in [17], in which the time integral of the quadratic sum of the gimbal rates was minimised. A trajectory-planning approach was described in [18] to reduce the possibility of CMG saturation while the system followed a reference trajectory. Another novel deterministic allocation algorithm for computing the angular rates of the gimbal was also proposed in [19] instead of the well-known pseudo-inverse technique in order to keep the keep a four-CMG cluster in a roof configuration away from singularities. An optimisation technique for reconfiguring the gimbals back to their initial positions after completing the desired manoeuvre was also addressed.

The utilisation of more than one type of control moment exchange devices has also been studied. The combination of a reaction wheel with CMGs presented an enhanced performance compared to reaction wheel clusters in terms of agility, thus providing an affordable alternative to the conventional four-CMG cluster [20]. A control strategy that exploited both thrusters and reaction wheels was presented in [21]. Thrusters were used in an open-loop scheme to drive the spacecraft to a reference path, while reaction wheels were used for a closed-loop variable-structure control law to eliminate tracking errors due to thruster biases and initial condition errors. Another simultaneous actuation of thrusters and reaction wheels through a varying time-sharing logic was studied, and it was shown that it could be efficiently adopted for time-optimal manoeuvres [22]. Additionally, three tracking control laws that made use of momentum wheels and thrusters were explored in [23], while the implementation of an SGCMG along three orthogonal momentum wheels was tested on the Distributed Spacecraft Attitude Control System Simulator and showed promising results [24].

While a CMG cluster can easily encounter singularities, a variable-speed CMG (VSCMG) cluster deals with this problem efficiently [25]. It uses four extra degrees of freedom (DoF) compared to the SGCMG cluster by changing the angular velocity of each flywheel independently. This way, a VSCMG cluster is capable of behaving as a reaction wheel cluster near a singularity and as a conventional CMG cluster away from it [26]. Power management capabilities have also been studied for VSCMGs [27,28,29,30], whilst a novel approach of combining a double-gimbal CMG and a VSCMG was presented in [31,32]. Similarly to VSCMGs, the proposed model described in this paper makes use of more than four DoFs. The extra DoF is implemented by adding a simple, miniature, and inexpensive low-power stepper motor in order to gimbal the full SGCMG cluster. Advances in small-motor mechatronics, avionics, and pin-pullers have enabled the miniaturisation of motors by a significant amount, thus creating space for small and compact VSCMG clusters [25,33]. In addition, low-power COTS stepper motors have been proven to significantly reduce CMG power compared to RW systems [34,35]. This paper analyses the ability to use a small stepper motor to gimbal a full SGCMG cluster, making this a hybrid approach by combining the elements of SGCMGs and VSCMGs without the additional complexity of a full VSCMG system. Moreover, the design explored here does not utilise actuators with low torque capabilities, such as reaction and momentum wheels, thus preserving the agility characteristics of a simple four-CMG cluster. The details of the proposed GCMG, along with the advantages and the shortcomings of the proposed method, are presented in this paper.

2. Mathematical Modeling

2.1. VSCMG Equations of Motion

The equation of motion of a rigid spacecraft deployed with four VSCMGs is described by:

$$\mathbf{J}\dot{\omega }=-\omega \times \mathbf{J}\dot{\omega }-{\mathbf{A}}_s{\mathbf{t}}_s-{\mathbf{A}}_t{\mathbf{t}}_t-{\mathbf{A}}_g{\mathbf{t}}_g$$

(1)

where matrices ${\mathbf{A}}_s,{\mathbf{A}}_t,{\mathbf{A}}_g\in {\Re }^{3\times 4}$ contain the components of the unit direction vector of each VSCMG frame of the cluster as:

$$\begin{array}{c}{\mathbf{A}}_s=[{\widehat{\mathbf{g}}}_{s_1},{\widehat{\mathbf{g}}}_{s_2},{\widehat{\mathbf{g}}}_{s_3},{\widehat{\mathbf{g}}}_{s_4}]\end{array}$$

(2)

$$\begin{array}{c}{\mathbf{A}}_t=[{\widehat{\mathbf{g}}}_{t_1},{\widehat{\mathbf{g}}}_{t_2},{\widehat{\mathbf{g}}}_{t_3},{\widehat{\mathbf{g}}}_{t_4}]\end{array}$$

(3)

$$\begin{array}{c}{\mathbf{A}}_g=[{\widehat{\mathbf{g}}}_{g_1},{\widehat{\mathbf{g}}}_{g_2},{\widehat{\mathbf{g}}}_{g_3},{\widehat{\mathbf{g}}}_{g_4}]\end{array}$$

(4)

and $\mathbf{J}$ is the spacecraft’s inertia matrix. The spin, gimbal and transverse axes are denoted by ${\widehat{\mathbf{g}}}_{s_N},{\widehat{\mathbf{g}}}_{g_N},{\widehat{\mathbf{g}}}_{t_N}\in {\Re }^{3\times 1}$, respectively, for the Nth VSCMG. $\omega \in {\Re }^{3\times 1}$ represents the angular velocity of the spacecraft expressed and derived in terms of control axes that are fixed relative to the body of the spacecraft. In this paper, the notation $\dot{(·)}$ indicates the time derivative of $(·)$. The torques ${\mathbf{t}}_s,{\mathbf{t}}_t,{\mathbf{t}}_g\in {\Re }^{4\times 1}$ can be calculated through

$$\begin{array}{c}{\mathbf{t}}_s={[t_{s_1},t_{s_2},t_{s_3},t_{s_4}]}^T\end{array}$$

(5)

$$\begin{array}{c}{\mathbf{t}}_t={[t_{t_1},t_{t_2},t_{t_3},t_{t_4}]}^T\end{array}$$

(6)

$$\begin{array}{c}{\mathbf{t}}_g={[t_{g_1},t_{g_2},t_{g_3},t_{g_4}]}^T\end{array}$$

(7)

and

$$\begin{array}{c}{t}_{s_N}=J_s({\dot{\Omega }}_N+{\dot{\delta }}_N{\omega }_{t_N})-(J_{t_N}-J_{g_N}){\dot{\delta }}_N{\omega }_{t_N}\end{array}$$

(8)

$$\begin{array}{c}{t}_{t_N}=J_s({\Omega }_N+{\omega }_{s_N}){\dot{\delta }}_N-(J_{t_N}-J_{g_N}){\omega }_{s_N}{\dot{\delta }}_N+J_s{\Omega }_N{\omega }_{g_N}\end{array}$$

(9)

$$\begin{array}{c}{t}_{g_N}=J_g{\ddot{\delta }}_N-J_s{\Omega }_N{\omega }_{t_N}\end{array}$$

(10)

where the inertia about the spin, gimbal and transverse axis is given by $J_s,J_g$ and $J_t$. For simplicity, these inertias are considered to be the same for all VSCMGs. Let $\delta ={[{\delta }_1,{\delta }_2,{\delta }_3,{\delta }_4]}^T\in {\Re }^{4\times 1}$ be the gimbal angles vector for a spacecraft employed with 4 CMGs. The vectors $\dot{\delta },\ddot{\delta }\in {\Re }^{N\times 1}$ and $\dot{\Omega }\in {\Re }^{4\times 1}$ contain the gimbal rate, gimbal acceleration and flywheel angular velocity, respectively. The stability condition as given in [36,37] is described by:

$${\mathbf{C}}_0\dot{\Omega }+\mathbf{B}\ddot{\delta }+({\mathbf{C}}_1-{\mathbf{C}}_2)\dot{\delta }={\mathbf{T}}_c$$

(11)

where

$$\begin{array}{c}\mathbf{B}=[{\widehat{\mathbf{g}}}_{g_1}{J}_g,{\widehat{\mathbf{g}}}_{g_2}{J}_g,{\widehat{\mathbf{g}}}_{g_3}{J}_g,{\widehat{\mathbf{g}}}_{g_4}{J}_g]\end{array}$$

(12)

$$\begin{array}{c}{\mathbf{C}}_0=[{\widehat{\mathbf{g}}}_{s_1}{J}_s,{\widehat{\mathbf{g}}}_{s_2}{J}_s,{\widehat{\mathbf{g}}}_{s_3}{J}_s,{\widehat{\mathbf{g}}}_{s_4}{J}_s]\end{array}$$

(13)

$$\begin{array}{c}{\mathbf{C}}_1=[{\widehat{\mathbf{g}}}_{t_1}{J}_s({\Omega }_1+{\omega }_{s_1}),{\widehat{\mathbf{g}}}_{t_2}{J}_s({\Omega }_2+{\omega }_{s_2}),{\widehat{\mathbf{g}}}_{t_3}{J}_s({\Omega }_3+{\omega }_{s_3}),{\widehat{\mathbf{g}}}_{t_4}{J}_s({\Omega }_4+{\omega }_{s_4})]\end{array}$$

(14)

$$\begin{array}{c}{\mathbf{C}}_2=[{\widehat{\mathbf{g}}}_{t_1}{J}_t{\omega }_{s_1},{\widehat{\mathbf{g}}}_{t_2}{J}_t{\omega }_{s_2},{\widehat{\mathbf{g}}}_{t_3}{J}_t{\omega }_{s_3},{\widehat{\mathbf{g}}}_{t_4}{J}_t{\omega }_{s_4}]\end{array}$$

(15)

To calculate the control torque ${\mathbf{T}}_c\in {\Re }^{3\times 1}$, the satellite’s kinematics equation of motion is expressed in quaternion form:

$$\dot{\mathbf{q}}=\frac{1}{2}\mathbf{q}\odot {\omega }_{\mathit{q}}$$

(16)

where ⊙ denotes the quaternion multiplication and $\mathbf{q}={[q_0,q_1,q_2,q_3]}^T\in {\Re }^{4\times 1}$ is the attitude quaternion. ${\omega }_{\mathit{q}}$ is the angular velocity vector of the satellite, $\omega$, given in a quaternion form as ${\omega }_{\mathit{q}}={[0,{\omega }^T]}^T\in {\Re }^{4\times 1}$. For the implementation of the attitude control, the control torque ${\mathbf{T}}_c$ applied on the satellite’s body is a function of the vector part of the error quaternion ${\mathbf{q}}_{err}$ and $\omega$ as described by the following equation:

$${\mathbf{T}}_c=K_p{\mathbf{q}}_{err}^v+K_{\omega }\omega$$

(17)

where $K_p$ and $K_w$ are fixed gains. The error quaternion between the current attitude quaternion and a desired quaternion ${\mathbf{q}}_{des}$ as:

$${\mathbf{q}}_{err}=\left[\begin{array}{c}{q}_{err}^s\\ {\mathbf{q}}_{err}^v\end{array}\right]={\mathbf{q}}_{des}^{*}\odot \mathbf{q}$$

(18)

where $q_{err}^s$ and ${\mathbf{q}}_{err}^v={[q_{err}^{roll},q_{err}^{pitch},q_{err}^{yaw}]}^T$ are the scalar and the vector part, respectively, and ${\mathbf{q}}_{des}^{*}$ expresses the conjugate quaternion of ${\mathbf{q}}_{des}$. The normalization of the quaternions is required before evaluating the ${\mathbf{q}}_{err}$. The “3-2-1” sequence is used to convert the quaternion error to the corresponding Euler angles error.

Assuming that the inertia of each VSCMG about the spin axis is much larger compared to the inertia in gimbal and transverse axis, the stability condition is

$$\mathbf{R}\dot{\mathit{p}}={\mathbf{T}}_c$$

(19)

where $\mathbf{R}\in {\Re }^{3\times 8}$ and $\dot{\mathit{p}}\in {\Re }^{8\times 1}$

$$\mathbf{R}=\left[\begin{array}{cc}{\mathbf{C}}_0& {\mathbf{C}}_1\end{array}\right]$$

(20)

$$\dot{\mathit{p}}=\left[\begin{array}{c}\dot{\Omega }\\ \dot{\delta }\end{array}\right]$$

(21)

Thus, the steering law in compact form is given by:

$$\dot{\mathit{p}}={\mathbf{R}}^{\#}{\mathbf{T}}_c$$

(22)

The matrix $\mathbf{R}$ is the Jacobian matrix of the system and except for the gimbal angles, the Jacobian matrix depends on geometric characteristics of the CMGs like the skew angle. Since the Jacobian matrix is not rectangular and the inverse of the matrix cannot be calculated, several definitions have been proposed for the inverse of the Jacobian matrix ${\mathbf{R}}^{\#}$. Using a weighting matrix $\mathbf{W}=diag\left([{\mathbf{W}}_{RW},{\mathbf{W}}_{CMG}]\right)\in {\Re }^{8\times 8}$ where

$$\begin{array}{c}{\mathbf{W}}_{RW}=W_{{RW}_0}{e}^{-\zeta m}\hfill \end{array}$$

(23)

and $\zeta ,W_{{RW}_0}$ are constant gains and ${\mathbf{W}}_{RW},{\mathbf{W}}_{CMG}\in {\Re }^{1\times 4}$, the matrix ${\mathbf{R}}^{\#}$ is

$${\mathbf{R}}^{\#}=\mathbf{W}{\mathbf{R}}^{T}inv\left(\mathbf{R}\mathbf{W}{\mathbf{R}}^T\right)$$

(24)

m is a measure of the linear independence of the columns of the matrix ${\mathbf{C}}_1$, usually referred as manipulability index:

$$m=\sqrt{det\left({\mathbf{C}}_1{{\mathbf{C}}_1}^T\right)}$$

(25)

This way, the VSCMG cluster is capable of performing as the conventional SGCMG cluster away from the singularity but as the manipulability index value decreases, the variable speed capability is exploited. The null space motion can be added upon the given steering law, assuming that ${\delta }_{des}$ and ${\Omega }_{des}$ are the vectors that contain the desired values for the gimbal angles and the flywheels’ angular velocity. The steering law is described below:

$$\dot{\mathit{p}}={\mathbf{R}}^{\#}{\mathbf{T}}_c+\rho ({\mathbf{R}}^{\#}\mathbf{R}-{\mathbf{I}}_{8\times 8})\mathbf{G}\left[\begin{array}{c}{\Omega }_{err}\\ {\delta }_{err}\end{array}\right]$$

(26)

where ${\mathbf{I}}_{8\times 8}$ is the identity matrix and $\rho$ a constant gain. The diagonal matrix $\mathbf{G}\in {\Re }^{8\times 8}$ is given by

$$\mathbf{G}=\left[\begin{array}{cc}{G}_{RW}{\mathbf{I}}_{4\times 4}& {\mathbf{0}}_{4\times 4}\\ {\mathbf{0}}_{4\times 4}& G_{CMG}{\mathbf{I}}_{4\times 4}\end{array}\right]$$

(27)

where ${\mathbf{0}}_{4\times 4}$ is the 4 by 4 zero matrix and the sets ${\Omega }_{err},{\delta }_{err}\in {\Re }^{4\times 1}$ represent the error with respect to the desired values as

$$\begin{array}{c}{\Omega }_{err}=\Omega -{\Omega }_{des}\end{array}$$

(28)

$$\begin{array}{c}{\delta }_{err}=\delta -{\delta }_{des}\end{array}$$

(29)

The $G_{RW}$ and $G_{CMG}$ values are set to 1 if it is desired the system to follow the ${\Omega }_{des}$ and ${\delta }_{des}$, respectively. Otherwise, they are set to 0. The total momentum of the VSCMG cluster $\mathbf{h}\in {\Re }^{3\times 1}$ is calculated by the sum of the individual momentum from each VSCMG:

$$\mathbf{h}={\mathbf{h}}_1+{\mathbf{h}}_2+{\mathbf{h}}_3+{\mathbf{h}}_4=\left[\begin{array}{c}-c\beta s{\delta }_1{h}_1-c{\delta }_2{h}_2+c\beta s{\delta }_3{h}_3+c{\delta }_4{h}_4\\ c{\delta }_1{h}_1-c\beta s{\delta }_2{h}_2-c{\delta }_3{h}_3+c\beta s{\delta }_4{h}_4\\ s\beta s{\delta }_1{h}_1+s\beta s{\delta }_2{h}_2+s\beta s{\delta }_3{h}_3+s\beta s{\delta }_4{h}_4\end{array}\right]$$

(30)

where s, c are the abbreviations for sin and cos, respectively. The parameter $\beta$ denotes the skew angle of the 4-CMG cluster in pyramid configuration and it is chosen in such a way that the momentum envelope is nearly tree-axis symmetric and spherical [6]. $h_i,i=1,\dots ,4$ are given by:

$$h_i=J_s{\Omega }_i$$

(31)

2.2. GCMG Equations of Motion

Adding an extra motor in the conventional SGCMG cluster it is possible to rotate the whole pyramid structure along the z direction, for ${\theta }_z$ degrees, w.r.t the body frame as shown in Figure 1. For a symmetric SGCMG cluster, its products of inertia are equal to zero and the angular momentum of the cluster, here referred as ${\mathbf{h}}_{GCMG}$ is given by Equation (32).

$${\mathbf{h}}_{GCMG}={\mathbf{h}}_{cmg}+{\mathbf{h}}_{stepper}=\left[\begin{array}{ccc}cos\left({\theta }_z\right)& -sin\left({\theta }_z\right)& 0\\ sin\left({\theta }_z\right)& cos\left({\theta }_z\right)& 0\\ 0& 0& 1\end{array}\right]\begin{array}{c}\mathbf{h}\end{array}+J_{cmg}^{zz}\left[\begin{array}{c}0\\ 0\\ {\dot{\theta }}_z\end{array}\right]$$

(32)

where $J_{cmg}^{zz}$ is the moment of inertia about the z axis of the SGCMG cluster. The total angular momentum $\mathbf{H}$ of the system is

$$\mathbf{H}=\mathbf{J}\omega +{\mathbf{h}}_{GCMG}$$

(33)

Let the moment of inertia of the cluster be much smaller than this of the satellite. It applies that changing the angle ${\theta }_z$ negligibly affects the orientation of the satellite. As a result, the stepper motor only rotates the actuator and the rate of the total angular momentum $\mathbf{H}$ can be calculated by:

$$\dot{\underset{¯}{\mathbf{H}}}=\left[\begin{array}{cc}\underset{¯}{\mathbf{A}}& \underset{¯}{\mathbf{D}}\end{array}\right]\left[\begin{array}{c}\dot{\delta }\\ \dot{{\theta }_z}\end{array}\right]+J_{cmg}^{zz}\left[\begin{array}{c}0\\ 0\\ {\ddot{\theta }}_z\end{array}\right]+\mathbf{J}\dot{\omega }+\omega \times \mathbf{H}$$

(34)

Because the disturbance torques are neglected and the total angular momentum is preserved, Equation (34) becomes:

$$\mathbf{0}=\left[\begin{array}{cc}\underset{¯}{\mathbf{A}}& \underset{¯}{\mathbf{D}}\end{array}\right]\left[\begin{array}{c}\dot{\delta }\\ \dot{{\theta }_z}\end{array}\right]+J_{cmg}^{zz}\left[\begin{array}{c}0\\ 0\\ {\ddot{\theta }}_z\end{array}\right]+\mathbf{J}\dot{\omega }+\omega \times {\mathbf{h}}_{cmg}+\omega \times J_{cmg}^{zz}\left[\begin{array}{c}0\\ 0\\ {\dot{\theta }}_z\end{array}\right]+\omega \times \mathbf{J}\omega$$

(35)

Let $h_1=h_2=h_3=h_4=h_f$, for identical flywheels running at the same angular speed. Then,

$$\begin{array}{cc}\hfill \underset{¯}{\mathbf{A}}& =h_f\left[\begin{array}{cc}s{\delta }_{1}s{\theta }_z-c\beta c{\delta }_{1}c{\theta }_z& s{\delta }_{2}c{\theta }_z+c\beta c{\delta }_{2}s{\theta }_z\\ -s{\delta }_{1}c{\theta }_z-c\beta c{\delta }_{1}s{\theta }_z& s{\delta }_{2}s{\theta }_z-c\beta c{\delta }_{2}c{\theta }_z\\ c{\delta }_{1}s\beta & c{\delta }_{2}s\beta \end{array}\right.\hfill \\ & \left.\begin{array}{cc}c\beta c{\delta }_{3}c{\theta }_z-s{\delta }_{3}s{\theta }_z& -s{\delta }_{4}c{\theta }_z-c\beta c{\delta }_{4}s{\theta }_z\\ s{\delta }_{3}c{\theta }_z+c\beta c{\delta }_{3}s{\theta }_z& c\beta c{\delta }_{4}c{\theta }_z-s{\delta }_{4}s{\theta }_z\\ c{\delta }_{3}s\beta & c{\delta }_{4}s\beta \end{array}\right]\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill \underset{¯}{\mathbf{D}}& =h_f\left[\begin{array}{c}c{\delta }_{2}s{\theta }_z-c{\delta }_{4}s{\theta }_z-c{\delta }_{1}c{\theta }_z+c{\delta }_{3}c{\theta }_z+\\ c{\delta }_{3}s{\theta }_z-c{\delta }_{1}s{\theta }_z-c{\delta }_{2}c{\theta }_z+c{\delta }_{4}c{\theta }_z-\\ 0\end{array}\right.\hfill \\ \hfill & \left.\begin{array}{c}+c\beta s{\delta }_{2}c{\theta }_z-c\beta s{\delta }_{4}c{\theta }_z+c\beta s{\delta }_{1}s{\theta }_z-c\beta s{\delta }_{3}s{\theta }_z\\ -c\beta s{\delta }_{1}c{\theta }_z+c\beta s{\delta }_{3}c{\theta }_z+c\beta s{\delta }_{2}s{\theta }_z-c\beta s{\delta }_{4}s{\theta }_z\\ 0\end{array}\right]\hfill \end{array}$$

(37)

Using the Euler formula, the Jacobian matrix of the system $\mathbf{R}\in {\Re }^{3\times 5}$ is formulated as:

$$\mathbf{R}=\left[\begin{array}{cc}\underset{¯}{\mathbf{A}}& \underset{¯}{\mathbf{D}}+\mathbf{M}\end{array}\right]$$

(38)

where

$$\mathbf{M}=J_{cmg}^{zz}\left[\begin{array}{c}{\omega }_y\\ -{\omega }_x\\ 0\end{array}\right]$$

(39)

The Equation (35) is given in a compact form by

$$\mathbf{0}=\mathbf{R}\left[\begin{array}{c}\dot{\delta }\\ \dot{{\theta }_z}\end{array}\right]+J_{cmg}^{zz}\left[\begin{array}{c}0\\ 0\\ {\ddot{\theta }}_z\end{array}\right]+\omega \times {\mathbf{h}}_{cmg}+\mathbf{J}\dot{\omega }+\omega \times \mathbf{J}\omega$$

(40)

The control vector ${\mathbf{T}}_c$ is selected as

$${\mathbf{T}}_c=\mathbf{R}\left[\begin{array}{c}\dot{\delta }\\ \dot{{\theta }_z}\end{array}\right]+J_{cmg}^{zz}\left[\begin{array}{c}0\\ 0\\ {\ddot{\theta }}_z\end{array}\right]+\omega \times {\mathbf{h}}_{cmg}$$

(41)

The gimbal angles as well as the rate of ${\theta }_z$ can be found by

$$\left[\begin{array}{c}\dot{\delta }\\ {\dot{\theta }}_z\end{array}\right]={\mathbf{R}}^{\#}({\mathbf{T}}_c-J_{cmg}^{zz}\left[\begin{array}{c}0\\ 0\\ {\ddot{\theta }}_z\end{array}\right]-\omega \times {\mathbf{h}}_{cmg})$$

(42)

where

$${\mathbf{R}}^{\#}={\mathbf{R}}^T{\left(\mathbf{R}{\mathbf{R}}^T\right)}^{-1}$$

(43)

is the standard pseudo-inverse, or commonly referred as Moore-Penrose steering law, used due to its simplicity and the real-time implementation it offers. To guarantee the attitude stability of the spacecraft, consider the following the candidate Lyapunov function [38,39,40]

$$V=\frac{1}{2}{\omega }^{\mathrm{T}}\mathbf{J}\omega +K_p\mathbf{q}{{}_{err}^v}^{\mathrm{T}}\mathbf{q}{}_{err}^v+K_p{\left(1-q_0\right)}^2$$

(44)

The first time derivative of V is given by

$$\dot{V}={\omega }^T\mathbf{J}\dot{\omega }+K_p\dot{\mathbf{q}}{{}_{err}^v}^T\mathbf{q}{}_{err}^v+K_p\mathbf{q}{{}_{err}^v}^T\dot{\mathbf{q}}{}_{err}^v-2K_p\left(1-q_0\right){\dot{q}}_0$$

(45)

Because $\mathbf{q}{{}_{err}^v}^{\mathbf{T}}\dot{\mathbf{q}}{}_{err}^v$ is a scalar, it can be shown that

$$\mathbf{q}{{}_{err}^v}^T\dot{\mathbf{q}}{}_{err}^v={\left(\dot{\mathbf{q}}{{}_{err}^v}^T\mathbf{q}{}_{err}^v\right)}^T=\dot{\mathbf{q}}{{}_{err}^v}^T\mathbf{q}{}_{err}^v$$

(46)

and the Equation (45) becomes

$$\dot{V}={\omega }^T\mathbf{J}\dot{\omega }+2K_p\dot{\mathbf{q}}{{}_{err}^v}^T\mathbf{q}{}_{err}^v-2K_p\left(1-q_0\right){\dot{q}}_0$$

(47)

Equation (47) becomes

$$\dot{V}={\omega }^T\mathbf{J}\dot{\omega }+2K_p{(\frac{1}{2}S\left(\omega \right){\mathbf{q}}_{err}^v+\frac{1}{2}{q}_0\omega )}^T\mathbf{q}{}_{err}^v+2K_p\left(1-q_0\right)\frac{1}{2}{\omega }^T{\mathbf{q}}_{err}^v$$

(48)

which can be simplified to

$$\dot{V}={\omega }^T\mathbf{J}\dot{\omega }+K_p({{\mathbf{q}}_{err}^v}^T{S\left(\omega \right)}^T+q_0{\omega }^T)\mathbf{q}{}_{err}^v+K_p\left(1-q_0\right){\omega }^T{\mathbf{q}}_{err}^v$$

(49)

Because $\mathbf{q}{{}_{err}^v}^{T}S\left(\omega \right)\mathbf{q}{}_{err}^v=0$, where $S\left(\omega \right)$ is the skew matrix of $\omega$, the Equation (49) becomes

$$\dot{V}={\omega }^T(\mathbf{J}\dot{\omega }+K_p\mathbf{q}{}_{err}^v)$$

(50)

Using Equations (17), (40) and (41), the Equation (50) can be simplified to

$$\dot{V}={\omega }^T(-\omega \times \mathbf{J}\omega -K_w\omega )$$

(51)

Note that ${\omega }^T(\omega \times \mathbf{J}\omega )=0$. Finally,

$$\dot{V}=-{\omega }^T{K}_w\omega$$

(52)

and the global stability is guaranteed for $K_w>0$ [41].

The redundancy of the system allows to exploit motion in the null space and Equation (42) can be modified to contain the null motion term as

$$\left[\begin{array}{c}\dot{\delta }\\ {\dot{\theta }}_z\end{array}\right]={\mathbf{R}}^{\#}({\mathbf{T}}_c-J_{cmg}^{zz}\left[\begin{array}{c}0\\ 0\\ {\ddot{\theta }}_z\end{array}\right]-\omega \times {\mathbf{h}}_{cmg})+\kappa ({\mathbf{R}}^{\#}\mathbf{R}-{\mathbf{I}}_{5\times 5})\Gamma \left[\begin{array}{c}{\delta }_{err}\\ {{\theta }_z}_{err}\end{array}\right]$$

(53)

where $\kappa$ is a constant gain. $\Gamma$ is a diagonal matrix the values of which denote which elements of ${\left[\begin{array}{c}{\delta }_{err},{{\theta }_z}_{err}\end{array}\right]}^T$ are selected for the null motion. The possible values of each element of $\Gamma$ are 0 and 1 and ${{\theta }_z}_{err}$ is given by

$${{\theta }_z}_{err}={\theta }_z-{{\theta }_z}_{des}$$

(54)

Linearized Model

Using the Equations (17) and (42), the input vector $\mathbf{u}=\left[\begin{array}{c}\dot{\delta }\\ {\dot{\theta }}_z\end{array}\right]$ can be calculated by

$$\mathbf{u}={\mathbf{R}}^{\#}(K_p{\mathbf{q}}_{err}^v+K_{\omega }\omega -J_{cmg}^{zz}\left[\begin{array}{c}0\\ 0\\ {\ddot{\theta }}_z\end{array}\right]-\omega \times {\mathbf{h}}_{cmg})$$

(55)

which can be embedded in the Equation (40) resulting in

$$\mathbf{0}=K_p{\mathbf{q}}^v+K_{\omega }\omega +\mathbf{J}\dot{\omega }+\omega \times \mathbf{J}\omega$$

(56)

where ${\mathbf{q}}^v={[q_1,q_2,q_3]}^T\in {\Re }^{3\times 1}$ because ${\mathbf{q}}_{des}=[1,0,0,0]$ for simplicity. Then, the closed loop system is described by the Equation (56) and the reduced kinematics equation

$${\dot{\mathbf{q}}}^v=\frac{1}{2}\left[\begin{array}{ccc}\sqrt{1-q_1^2-q_2^2-q_3^2}& -q_3& q_2\\ q_3& \sqrt{1-q_1^2-q_2^2-q_3^2}& -q_1\\ -q_2& q_1& \sqrt{1-q_1^2-q_2^2-q_3^2}\end{array}\right]\omega$$

(57)

Let the spacecraft’s products of inertia be zero for simplicity, i.e., $\mathbf{J}=\mathrm{diag}\left([J_{11},J_{22},J_{33}]\right)$, in accordance to the simulation parameters presented in Section 3. Defining the state vector ${\mathbf{x}}_s\triangleq {[\omega ,{\mathbf{q}}^v]}^T\in {\Re }^{6\times 1}$, the linearization of the nonlinear model results in

$${\dot{\mathbf{x}}}_s=\left[\begin{array}{cccc}& {\mathbf{M}}_{K_{\omega }}& & \multicolumn{1}{c}{{\mathbf{M}}_{K_p}}\\ \frac{1}{2}\sqrt{1-q_1^2-q_2^2-q_3^2}& -q_3/2& q_2/2& \\ q_3/2& \frac{1}{2}\sqrt{1-q_1^2-q_2^2-q_3^2}& -q_1/2& {\mathbf{M}}_{\omega }\\ -q_2/2& q_1/2& \frac{1}{2}\sqrt{1-q_1^2-q_2^2-q_3^2}& \end{array}\right]{\mathbf{x}}_s$$

(58)

where ${\mathbf{M}}_{K_{\omega }},{\mathbf{M}}_{K_p},{\mathbf{M}}_{\omega }\in {\Re }^{3\times 3}$ are defined by

$$\begin{array}{c}{\mathbf{M}}_{K_{\omega }}=\left[\begin{array}{ccc}\frac{-K_{\omega }}{J_{11}}& \frac{{\omega }_3(J_{22}-J_{33})}{J_{11}}& \frac{{\omega }_2(J_{22}-J_{33})}{J_{11}}\\ -\frac{{\omega }_3(J_{11}-J_{33})}{J_{22}}& \frac{-K_{\omega }}{J_{22}}& -\frac{{\omega }_1(J_{11}-J_{33})}{J_{22}}\\ \frac{{\omega }_2(J_{11}-J_{22})}{J_{33}}& \frac{{\omega }_1(J_{11}-J_{22})}{J_{33}}& \frac{-K_{\omega }}{J_{33}}\end{array}\right]\end{array}$$

(59)

$$\begin{array}{c}{\mathbf{M}}_{K_p}=\left[\begin{array}{ccc}\frac{-K_p}{J_{11}}& 0& 0\\ 0& \frac{-K_p}{J_{22}}& 0\\ 0& 0& \frac{-K_p}{J_{33}}\end{array}\right]\end{array}$$

(60)

$$\begin{array}{c}{\mathbf{M}}_{\omega }=\left[\begin{array}{ccc}0& \frac{{\omega }_3}{2}& -\frac{{\omega }_2}{2}\\ \frac{-{\omega }_3}{2}& 0& \frac{{\omega }_1}{2}\\ \frac{{\omega }_2}{2}& -\frac{{\omega }_1}{2}& 0\end{array}\right]\end{array}$$

(61)

For the set point ${\mathbf{x}}_s={\mathbf{x}}_s^e={[0,0,0,0,0,0]}^T$, the Equation (58) becomes

$${\dot{\mathbf{x}}}_s=\mathbf{Q}{\mathbf{x}}_s=\left[\begin{array}{cccccc}\frac{-K_p}{J_{11}}& 0& 0& & & \\ 0& \frac{-K_p}{J_{22}}& 0& & {\mathbf{M}}_{K_p}& \\ 0& 0& \frac{-K_p}{J_{33}}& & & \\ \frac{1}{2}& 0& 0& 0& 0& 0\\ 0& \frac{1}{2}& 0& 0& 0& 0\\ 0& 0& \frac{1}{2}& 0& 0& 0\end{array}\right]{\mathbf{x}}_s$$

(62)

The eigenvalues of $\mathbf{Q}$ are

$$eig\left(\mathbf{Q}\right)=\left[\begin{array}{c}-\frac{K_{\omega }-\sqrt{K_{\omega }^2-2J_{11}{K}_p}}{2J_{11}}\\ -\frac{K_{\omega }+\sqrt{K_{\omega }^2-2J_{11}{K}_p}}{2J_{11}}\\ -\frac{K_{\omega }-\sqrt{K_{\omega }^2-2J_{22}{K}_p}}{2J_{22}}\\ -\frac{K_{\omega }+\sqrt{K_{\omega }^2-2J_{22}{K}_p}}{2J_{22}}\\ -\frac{K_{\omega }-\sqrt{K_{\omega }^2-2J_{33}{K}_p}}{2J_{33}}\\ -\frac{K_{\omega }+\sqrt{K_{\omega }^2-2J_{33}{K}_p}}{2J_{33}}\end{array}\right]$$

(63)

The asymptotic stability is guaranteed for the parameters presented in Section 3 since every eigenvalue of $\mathbf{Q}$ has a negative real part.

3. Singularities—Preliminaries

As initially stated, the main problem when SGCMGs are utilized for spacecraft attitude control, is singularities. A brief overview of the singularities is presented, in order to better explain the selected task as presented in the next section. Near a singularity, the SGCMG cluster is incapable of producing torque along at least one direction and due to the ill-conditioned/singular Jacobian of the system in this state, the gimbal rates approach infinity. Even though it is known that only three SGCMG are necessary for three-axis control, the utilization of a redundant SGCMG provides the opportunity to avoid some types of singularities, increasing the momentum and the torque capacity of the cluster in a mass-efficient way.

Not all singularities are of the same type. They are divided into two categories. If a singular state can be escaped by null motion, it is passable or hyperbolic. Otherwise, the singularity is elliptic. The null motion of a SGCMG cluster is defined as the SGCMG gimbal motions that produce no net torque, i.e., the actuators’ generated torque does not change the attitude of the satellite. Most of the times, the null motion is exploited to avoid singularities by internally configuring the gimbal angles to better performance configurations without affecting the execution of the manoeuvre towards to commanded attitude.

In addition, singularities are sub-categorized into external and internal. External singularities correspond to the maximum angular momentum in any single direction. Such singularities, also known as saturation singularities, are accommodated in the sizing of the control system and they are not considered for singularity avoidance. In contrast, internal singularities occur when the system is in a singularity and the cluster momentum is far from the limit of the cluster’s capability. Depending on the number of the momentum axes aligned, the singularities of a 4-SGCMG cluster are classified to 4 h, 2 h and 0 h (Figure 2).

Caution must be taken when interpreting the singular surfaces since the 4 h surfaces do not always have a momentum magnitude of 4 h [6]. When two of the SGCMGs point opposite the other two, then it is a 0 h singularity. Every momentum state corresponding to a singular gimbal set is a singular surface. When these 4 h, 2 h and 0 h surfaces are plotted together, they form a continuous connection between all possible singular momentum states as illustrated in Figure 3.

The detailed understanding on the location and inner working of singularities can be found in the literature [2,11,12,39,42,43,44] as it is not the main focus of the work presented. However, the mathematical framework for examining the null motion as well as for determining the singularity types is extensively detailed in [6,42,45,46,47] and is briefly presented below for the typical 4-SGCMG pyramid cluster.

The null motion condition for the null motion momentum vector $\delta {\mathit{h}}_N\in {\Re }^{3\times 1}$ is defined by

$$\delta \mathit{H}=\sum _{N=1}^4\delta {\mathit{h}}_N=0$$

(64)

where $\mathit{H}$ is the total angular momentum of the 4-SGCMG cluster and N corresponds to the Nth SGCMG of the cluster. Using the Taylor series, the null motion momentum vector $\delta {\mathit{h}}_N$ can be expanded to

$$\begin{array}{c}\delta \mathit{H}=\sum _{N=1}^4\delta {\mathit{h}}_N=\sum _{N=1}^4\left(\frac{d{\mathit{h}}_N}{d{x}_N}\delta x_N+\frac{1}{2!}\frac{d^2{\mathit{h}}_N}{d{x}_N^2}\delta x_N^2+\frac{1}{3!}\frac{d^3{\mathit{h}}_N}{d{x}_N^3}\delta x_N^3+\cdots \right)\end{array}$$

(65)

$$\begin{array}{c}=\sum _{N=1}^4\left({\mathit{f}}_N\delta x_N-\frac{1}{2!}{\mathit{h}}_N\delta x_N^2-\frac{1}{3!}{\mathit{f}}_N\delta x_N^3+\cdots \right)=0\end{array}$$

(66)

where ${\mathit{f}}_N=d{\mathit{h}}_N/d{x}_N$ and $\delta x_N$ are the null motion gimbal angle displacements. The first-order necessary condition for null motion is

$$\sum _{N=1}^4{\mathbf{f}}_N\delta x_N=\mathbf{A}\delta \mathbf{x}=0$$

(67)

where $\delta \mathbf{x}=\left(\delta x_1,\dots ,\delta x_n\right)$ and $\mathbf{A}$ is the Jacobian matrix of the system. The null vector $\delta \mathbf{x}$ is the null-space vector of $\mathbf{A};$ that is, $\delta \mathbf{x}=\mathbf{n}=null\left(\mathbf{A}\right)$ and $\mathbf{A}\mathbf{n}=0$. In order to define the type of a singular state, let the null motion constraint

$$\delta \mathit{H}=\mathit{H}(\mathbf{x}+\delta \mathbf{x})-\mathit{H}\left(\mathbf{x}\right)=\sum _{N=1}^4\left({\mathbf{f}}_N\delta x_N-\frac{1}{2!}{\mathbf{h}}_N\delta x_N^2-\frac{1}{3!}{\mathbf{f}}_N\delta x_N^3+\cdots \right)=0$$

(68)

which is an alternative form of Equation (65). Let $\mathbf{u}$ be an arbitrary vector. The inner product of $\delta \mathbf{H}$ with $\mathbf{u}$ results in

$${\mathbf{u}}^T\delta \mathit{H}={\mathbf{u}}^T\left\{\sum _{N=1}^4\left[{\mathbf{f}}_N\delta x_N-\frac{1}{2!}{\mathbf{h}}_i\left(x_N\right)\delta x_N^2-\frac{1}{3!}{\mathbf{f}}_N\delta x_N^3+\cdots \right]\right\}=0$$

(69)

Since ${\mathbf{u}}^T{\mathbf{f}}_N=0$ when $\mathbf{u}$ is along the singular direction, it applies that

$$\begin{array}{c}0={\mathbf{u}}^T\left[\sum _{N=1}^4\left(-\frac{1}{2!}{\mathbf{h}}_N\delta x_N^2+\frac{1}{4!}{\mathbf{h}}_N\delta x_N^n+\cdots \right)\right]\end{array}$$

(70)

$$\begin{array}{c}=\sum _{N=1}^4\left({\mathbf{u}}^T{\mathbf{h}}_N\right)\left(-\frac{1}{2!}\delta x_N^2+\frac{1}{4!}\delta x_N^n-\cdots \right)\end{array}$$

(71)

$$\begin{array}{c}=\sum _{N=1}^4{e}_N\left(cos\delta x_N-1\right)\end{array}$$

(72)

where $e_N={\mathbf{u}}^T{\mathbf{h}}_N$. Keeping only the second-order terms, the Equation (70) becomes

$$\sum _{N=1}^4{e}_N\delta x_N^2=\delta {\mathbf{x}}^T\mathbf{E}\delta \mathbf{x}=0$$

(73)

where $\mathbf{E}=diag\left(e_N\right).$ When $\mathbf{E}$ is a sign-definite matrix, the only solution to Equation (73) is $\delta \mathbf{x}=0$, and there is not null motion. The null motion of gimbal angles can be described by the first-order condition

$$\delta \mathbf{x}=\sum _{N=1}^2{c}_N{\mathbf{n}}_N=\mathbf{Z}\mathbf{c}$$

(74)

where $\mathbf{c}={[c_1,c_2]}^T$ and $c_N$ is the Nth weighting coefficient. ${\mathbf{n}}_N$ are the null-space basis vectors of the Jacobian matrix

$$\mathbf{A}=\left[\begin{array}{cccc}-c\beta c{\delta }_1& s{\delta }_2& c\beta c{\delta }_3& -s{\delta }_4\\ -s{\delta }_1& -c\beta c{\delta }_2& s{\delta }_3& c\beta c{\delta }_4\\ s\beta c{\delta }_1& s\beta c{\delta }_2& s\beta c{\delta }_3& s\beta c{\delta }_4\end{array}\right]$$

(75)

such that $\mathbf{A}{\mathbf{n}}_N=0$ or $\mathbf{Z}=null\left(\mathbf{A}\right).$ From Equations (73) and (74), it applies that

$${\mathbf{c}}^T\mathbf{V}\mathbf{c}=0$$

(76)

where $\mathbf{V}={\mathbf{Z}}^T\mathbf{EZ}$. If matrix $\mathbf{V}$ is sign-definite, then the singularity is elliptic and null motion is not possible. Otherwise, the singularity is hyperbolic and it is possible to exploit null motion to escape the singular state. For the 4-SGCMG cluster in pyramid configuration, it is easy to prove that the system is in an elliptic singular state when the gimbal configuration is $\delta ={[-90,0,90,0]}^T$ deg as follows

$$\begin{array}{c}\mathbf{A}(-90,0,90,0)=\left[\begin{array}{cccc}0& 0& 0& 0\\ 1& -0.5774& 1& 0.5774\\ 0& 0.8164& 0& 0.8164\end{array}\right]\end{array}$$

(77)

$$\begin{array}{c}\mathbf{Z}=\mathrm{null}\left(\mathbf{A}\right)=\left[\begin{array}{cc}-0.7071& -0.3536\\ 0.0000& -0.6124\\ 0.7071& -0.3536\\ -0.0000& 0.6124\end{array}\right]\end{array}$$

(78)

$$\begin{array}{c}\mathbf{u}=\mathrm{null}\left({\mathbf{A}}^T\right)={\left[\begin{array}{c}1,0,0\end{array}\right]}^T\end{array}$$

(79)

$$\begin{array}{c}\mathbf{E}=\mathrm{diag}\left({\mathbf{u}}^T{\mathbf{h}}_N\right)=\mathrm{diag}(0.5774,-1,0.5774,1)\end{array}$$

(80)

$$\begin{array}{c}\mathbf{V}={\mathbf{Z}}^T\mathbf{E}\mathbf{Z}=\left[\begin{array}{cc}0.5774& 0\\ 0& 0.1444\end{array}\right]\end{array}$$

(81)

The $\mathbf{V}$ matrix is sign-definite because it has a positive determinant. Thus, the gimbal configuration is $\delta ={[-90,0,90,0]}^T$ deg corresponds to an elliptic singularity. A visualization of the gimbals in this state is presented in Figure 4.

In terms of controllability, a singularity can also be examined considering the dynamics equation

$$\dot{\omega }={\mathbf{J}}^{-1}(-\omega \times \left(\mathbf{J}\omega \right)-\dot{\mathbf{h}}-\omega \times \mathbf{h})$$

(82)

and the Equation (57) where the state vector is defined as ${\mathbf{x}}_s\triangleq {[\omega ,{\mathbf{q}}^v]}^T\in {\Re }^{6\times 1}$ and ${\mathbf{u}}_s\triangleq \dot{\delta }\in {\Re }^{4\times 1}$. Without loss of generality, let $\mathbf{J}={\mathbf{I}}_{3\times 3}$ for simplicity. Then, for the set point ${\mathbf{x}}_s={\mathbf{x}}_s^e={[0,0,0,0,0,0]}^T$ and ${\mathbf{u}}_s={\mathbf{u}}_s^e={[0,0,0,0]}^T$, the open loop system can be described in the linearized form of ${\dot{\mathbf{x}}}_s={\mathbf{A}}_l{\mathbf{x}}_s+{\mathbf{B}}_l{\mathbf{u}}_s$ where

$$\begin{array}{c}{\mathbf{A}}_l=\left[\begin{array}{cccccc}0& -s\beta s{\delta }_1+s{\delta }_2+s{\delta }_3+s{\delta }_4& c{\delta }_1-c{\delta }_3-c\beta s{\delta }_2+c\beta s{\delta }_4& 0& 0& 0\\ s\beta s{\delta }_1+s{\delta }_2+s{\delta }_3+s{\delta }_4& 0& c{\delta }_2-c{\delta }_4+c\beta s{\delta }_1-c\beta s{\delta }_3& 0& 0& 0\\ c{\delta }_3-c{\delta }_1+c\beta s{\delta }_2-c\beta s{\delta }_4& c{\delta }_4-c{\delta }_2-c\beta s{\delta }_1+c\beta s{\delta }_3& 0& 0& 0& 0\\ 0.5& 0& 0& 0& 0& 0\\ 0& 0.5& 0& 0& 0& 0\\ 0& 0& 0.5& 0& 0& 0\end{array}\right]\end{array}$$

(83)

$$\begin{array}{c}{\mathbf{B}}_l=\left[\begin{array}{cccc}c\beta c{\delta }_1& -s{\delta }_2& -c\beta c{\delta }_3& s{\delta }_4\\ s{\delta }_1& c\beta c{\delta }_2& -s{\delta }_3& -c\beta c{\delta }_4\\ -c{\delta }_{1}s\beta & -c{\delta }_{2}s\beta & -c{\delta }_{3}s\beta & -c{\delta }_{4}s\beta \\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]\end{array}$$

(84)

Then, the controllability matrix can be calculated by

$${\mathbf{C}}_l=[{\mathbf{B}}_l,{\mathbf{A}}_l{\mathbf{B}}_l,{\mathbf{A}}_l^2{\mathbf{B}}_l,\dots ,{\mathbf{A}}_l^5{\mathbf{B}}_l]$$

(85)

and the $\mathrm{rank}\left({\mathbf{C}}_l\right)$ indicates whether the system is controllable or not. E.g., at the configuration $\delta ={[0,0,0,0]}^T$ deg, the $\mathrm{rank}\left({\mathbf{C}}_l\right)=6$ and the system is fully controllable in this state. In contrast, the well-known elliptic singularity at $\delta ={[-90,0,90,0]}^T$ deg results in $\mathrm{rank}\left({\mathbf{C}}_l\right)=4<6$ indicating the uncontrollability at the certain gimbal configuration.

4. Simulation Set-Up

Following the description of the singular states and the problems they introduce to the spacecraft attitude control, it is evident that exploring the behaviour of a system employed with SGCMGs near these states is essential. Thus, the task selected to study the behaviour of each model near an internal elliptic singularity is a manoeuvre about the roll axis by −90 deg because the system encounters the $\delta ={[-90,0,90,0]}^T$ deg singularity [46]. The initial non-zero gimbal configuration ${\delta }_{ini}$ near the singular configuration is selected to guarantee that the system will meet the singularity while the initial attitude quaternion is ${\mathbf{q}}_{ini}={[1,0,0,0]}^T$. Similar to [12,34,48], studying this type of singular state in addition to the initial gimbal values that are extremely near the singularity is a useful simulation strategy which enables to prove that the proposed GCMG design is effective even in challenging conditions. The total simulation time is 20 s. The first 2 s the systems is at rest, then the system is commanded to follow the desired orientation given by the quaternion ${\mathbf{q}}_{des}$. Since the simulation describes a discrete time system, it is required to integrate the quaternion rate $\dot{\mathbf{q}}$ at the ith iteration as given by Equation (16) using the equation:

$${\mathbf{q}}_i={\mathbf{q}}_{i-1}\odot [cos(\left|\right|\omega \left|\right|\frac{dt}{2}),\frac{\omega }{\left|\right|\omega \left|\right|}sin\left(\right|\left|\omega \right||\frac{dt}{2}{\left)\right]}^T$$

(86)

where $dt$ is the time between two consecutive iterations. The SGCMG rates are saturated when they exceed the specific threshold ${\dot{\delta }}_{th}$, using the following formula:

$${\dot{\delta }}_{sat}=\dot{\delta }\frac{{\dot{\delta }}_{th}}{max\left(\right|{\dot{\delta }}_1|,|{\dot{\delta }}_2|,|{\dot{\delta }}_3|,|{\dot{\delta }}_4\left|\right)}$$

(87)

It is preferred to saturate the gimbal angle rates using Equation (87) over applying a boundary value to every gimbal angle rate that exceeds the required threshold because the characteristics of the motion are conserved. For a fair comparison, it is desired to keep the number of the unrestricted DoF of the two models the same. Thus, every DoF apart from the motion of the four gimbals in the cluster is selected to follow certain criteria, i.e., every flywheel of the VSCMG cluster as well as the stepper motor of the GCMG cluster are used for null space motion following a preferred value. Moreover, the rotation angle and rate of the stepper motor are limited in order to better depict a real-life application.

A 2.7 GHz Inter Core i5 with 8 GB of RAM computer was used and all the simulations were run on the operating system MacOS Mojave in Matlab^® R2016b. The exact simulation parameters used, are shown in

Table 1. Simulation Parameters.

Parameter	VSCMG	GCMG
Manoeuvre	0 to −90 deg	0 to −90 deg
Moment of Inertia J	$\mathrm{diag}\left([1,1,1]\right){\mathrm{kgm}}^2$	$\mathrm{diag}\left([1,1,1]\right){\mathrm{kgm}}^2$
$J_s,J_t,J_g$	$9.5\times {10}^{-5}{\mathrm{kgm}}^2$	-----
$J_{cmg}^{zz}$	-----	38 $\times {10}^{-4}{\mathrm{kgm}}^2$
Time-step dt	0.01 s	0.01 s
Simulation time	20 s	20 s
PD, $K_p$, $K_{\omega }$	1.6, 3	1.6, 3
Skew angle $\beta$	54.73 deg	54.73 deg
${\dot{\delta }}_{th}$	50 deg/s	50 deg/s
${\delta }_N^{ini}$	${[-70,0,75,0]}^T$ deg	${[-70,0,75,0]}^T$ deg
${\dot{\theta }}_z^{sat}$	-----	50 deg/s
$[{\theta }_z^{-sat},{\theta }_z^{+sat}]$	-----	$[-180,180]$ deg
$\Omega$	$1047.2{[1,1,1,1]}^T$ rad/s	$1047.2{[1,1,1,1]}^T$ rad/s

. The exact steering law parameters used for each model are presented in

Table 2. VSCMG steering law parameters.

Parameter	Value
$\zeta$	${10}^{-4}$
$G_{RW},G_{CMG}$	1, 0
$\rho$	1
${\mathbf{W}}_{R{W}_0}$	$6\times {10}^5[1,1,1,1]$
${\mathbf{W}}_{CMG}$	$10[1,1,1,1]$
${\Omega }_{des}$	$1047.2{[1,1,1,1]}^T$ rad/s
${\delta }_{des}$	-----
$G_{RW}$	1
$G_{CMG}$	0

and

Table 3. GCMG steering law parameters.

Parameter	Value
$\kappa$	1
$\Gamma$	$\mathrm{diag}\left(\right[0,0,0,0,1\left]\right)$
${{\theta }_z}_{des}$	−90 deg

.

5. Simulation Results

Figure 5 illustrates the results derived by using the VSCMG cluster in avoiding the singularity encountered. The steering law of Equation (26) is applied and the null motion guides the system to follow only the desired ${\Omega }_{des}$ set because the $G_{CMG}$ equals zero.

The attitude angle in x axis is shown in Figure 5a. After t = 2 s the attitude begins to follow the commanded quaternion. The system reaches the steady state at approximately t = 18 s where the attitude error about the roll axis is slightly less than 1 deg. Figure 5b shows that the manoeuvre is executed with a minor attitude error about the pitch and yaw axis because the the system exploits the variable speed part to handle the approach to singularity. In particular, the maximum deviation of the attitude angle about pitch and yaw is −0.18 deg and 0.9 deg, respectively. The angular velocity of the satellite is presented in Figure 5c and the maximum absolute value is 18.9 deg/s in the negative of the roll axis. It is possible to accomplish a lower maximum velocity value by decreasing the value of the gain $K_p$. The momentum profile of the VSCMG cluster, is presented in Figure 5d. The momentum in x direction increases approaching the value of 1.15 Nms where the gimbals are in singular configuration [47]. Figure 5e shows that the angle of the first and the third gimbal is heading towards −90 deg and +90 deg, respectively, when the system approaches the singularity but in the steady state these values change to −74.3 deg and 71.1 deg. The rest of the gimbal angles remain almost zero through the whole simulation time. Figure 5f presents the gimbal rates. It is observed that mainly the first and the third gimbal contribute to the rotation whilst the rest gimbal angle rates are significantly lower. The rate of the first gimbal reaches the saturation for a slight period of time but in the steady state, no rates are generated by the steering law. The control torque as presented in Figure 5g is about the roll axis. From t = 2 s to t = 2.9 s the torque spans to the negative of roll axis but after t = 2.9 s the sign changes to decelerate the satellite until the steady state. The manipulability index drops from 0.446 (t = 0 s) to 0.149 (t = 2.9 s). It is known that the minimum value is obtained once the system reaches closer to the singular configuration. After that, the index increases till the steady state where the obtain value is equal to 0.47 (Figure 5h). The angular velocity of each flywheel in the cluster is presented in Figure 6. A small change in the angular velocity of each flywheels is capable of solving the singularity problem. The maximum deviation from the nominal velocity is 0.31% at t = 2 s but in the steady state all flywheels are following the ${\Omega }_{des}$ set.

Overall, the 4-VSCMG cluster is capable of avoiding the singular state due to the redundancy in the DoF it offers. The system avoids the elliptic singularity because it is capable to change the flywheel’s angular velocity while configuring the four gimbals. When the null space motion is applied, the system can efficiently track the desired set ${\Omega }_{des}$ and follow the commanded manoeuvre while maintaining a high performance configuration.

The same manoeuvre is executed using the GCMG model and the results are presented in Figure 7. Figure 7a shows that after t = 2 s the system follows the desired quaternion. The attitude angle about x direction is starting to decrease reaching the steady state at t = 18 s with a steady state error similar to this of the VSCMG. The GCMG exploits the fifth DoF to handle the encountering singularity completing the manoeuvre with slight attitude error (Figure 7b). Compared to the VSCMG, the maximum deviation is −1 deg and 0.53 deg about the pitch and yaw axis, respectively. This deviation is also depicted in the angular velocity of the satellite mainly in the negative of the y direction with an absolute maximum velocity of 4.4 deg/s (Figure 7c). The main velocity component is about the x axis and the maximum absolute value is 18.7 deg/s. The angular momentum of the GCMG cluster is presented in Figure 7d. When the momentum in x direction reaches the value of 1.15 Nms where is known that the system encounters the singular state, a change in the momentum about y is presented. The gimbal angles profile is shown in Figure 7e and in the steady state the gimbal angle vector $\delta$ is ${[2.3,4.9,128.5,-54]}^T$ deg. The rate of the third gimbal becomes saturated from t = 2 s to t = 2.37 s (Figure 7f) while every other gimbal rate is scaled down according to Equation (87) for this time period. The control torque can be seen in Figure 7g and the profile in x direction is similar to this of the VSCMG. It spans in the negative of the axis till t = 3.2 s where the overshoot in the positive of the axis indicates again the deceleration of the satellite. It can be seen in Figure 7h that the manipulability index of the CMG cluster, i.e., use of $\underset{¯}{\mathbf{A}}$—drops to 0.166 at t = 2.4 s but in the steady state it converges to 1.22. In contrast, the lowest value of the manipulability index of the GCMG cluster—use of $\mathbf{R}$—is 0.466 at t = 2.15 s. Even though the manipulability index in these two cases is not comparable because of the unequal DoF are used in each case, the profile of the index is. The significant difference is noticed from t = 2.16 s to t = 3 s. In this period of time the manipulability index of the GCMG increases quickly because of the action of the stepper motor, in contrast to this of the CMG which remains near the singularity. In Figure 8a,b the angle ${\theta }_z$ and rate ${\dot{\theta }}_z$ are presented. The stepper motor rotates the cluster −90 deg following the commanded ${{\theta }_z}_{des}$ set and the rate ${\dot{\theta }}_z$ remains saturated from t = 2 s to t = 2.5 s. Then, it decreases again because of the motion in the null space and in the steady state, it drops to zero.

Table 4. Percentage differences.

Parameter	Value
Max attitude error in y	455%
Max attitude error in z	−41%
Max velocity in x	−1%
Min manipulability	11%
Steady state manipulability	160%

presents the percentage differences for the main parameters among the two simulations.

As the results indicate, both the VSCMG and GCMG cluster are capable of avoiding the internal singularity exploiting the redundant DoF while maintaining a high manipulability index through the whole manoeuvre and similar settling time. Contrary to the results derived by using the GCMG, the maximum deviation of the attitude angle and satellite’s angular velocity is lower about pitch and higher about yaw axis, respectively, in the case of VSCMG. The same applies for the control torque. The manipulability index profile for the VSCMG is significantly different from this of the GCMG. It presents a wide dip and slowly increases until the end of the simulation. The GCMG presents a higher index value both in singularity and the steady state compared to the VSCMG. Thus, after completing the desired manoeuvre the GCMG cluster comes to a better configuration which is significant for the beginning of the next manoeuvre without re-configuring its gimbals, e.g., through null space.

No major difference is presented in the maximum absolute value of the angular velocity of the satellite about the roll axis for the two models. Both models are capable of avoiding the elliptic singular state with attitude error lower than 1 deg which makes the GCMG an efficient alternative to the VSCMG when it is desired to explicitly follow a given trajectory. The memory allocation is almost identical for the two models (approx. 2% less in the case of GCMG). The difference in execution time is less than 1% as measured by 100 consecutive measurements. However the mathematical complexity of the GCMG is highly reduced. In applications where only the final attitude is of interest while high attitude error is allowed during the manoeuvre, the extra motor of the GCMG design is adding weight and complexity without providing any extra value compared to the conventional SGCMG. Such a task can be efficiently accomplished by a SGCMG cluster combined with the modified singularity robust inverse steering law.

Caution must be taken in selecting the moment of inertia matrix of the CMG cluster. The used value is negligible compared to this of the satellite and as a result a minor rotation is applied to the system. A large CMG inertia may result in rotating the satellite in the opposite direction resulting in higher attitude errors.

Figure 9 illustrates an important difference between the total power required by the VSCMG compared to the GCMG to execute the same manoeuvre. The power component that is related to the initial spin up and the stabilization of the flywheels to their nominal velocities is the same in both cases and as a result, it has been omitted. The power profiles presented consist of two power components. The first is related to power used for the gimbal motion whilst the second is related to the power needed to change the velocity of the flywheels and control the stepper motor in the VSCMG and GCMG design, respectively. The power difference shown in Figure 9 is mainly due to the second component because the order of magnitude of the power of the first component is much smaller. There is an 88.1% and a 74.1% decrease in the maximum and the mean power needed when the GCMG is exploited. The maximum and the mean values for both cases are presented in

Table 5. Power values.

Power	GCMG	VSCMG
Maximum	0.153 W	1.28 W
Mean	0.006 W	0.023 W

. The values indicate that is feasible to use a commercial off-the-shelf stepper motor to manufacture the GCMG design with a required motor torque capability of less than 0.34 Nm. Such demands can be easily fulfilled by low mass and affordable stepper motors allowing the application of the GCMG actuator to micro and nano-satellites missions.

A more detailed power comparison of COTS based CMGs using stepper motors will be presented in future work showing experimental results which as shown previously in [34,35] show the power advantage of using stepper motors for CMG systems. The proposed GCMG uses less power than a conventional VSCMG cluster as the added gimbal motor/extra degree of freedom brings savings in power compared to the increased speed/power needed when using VSCMGs.

6. Conclusions

This paper addresses the problem of singularity avoidance in a 4-single gimbal control moment gyroscope (SGCMG) cluster in pyramid configuration as used for the attitude control of a satellite. A novel SGCMG cluster design is proposed that uses an extra stepper motor capable of rotating the full cluster and a comparison is held among the design proposed and the VSCMG cluster. Both are capable of avoiding the encountering singularity even though the GCMG presents a lower and a higher attitude error about the yaw and the pitch, respectively, for the exact same tuning. However, the GCMG shows a clear advantage over the VSCMG in the gimbal configuration in the end of the manoeuvre. The gimbals end in a configuration of higher performance that is essential for agile and small satellites that perform multiple consecutive manoeuvres in one pass. Moreover, the results indicate that GCMG provides a superior performance in terms of the power needed to execute the manoeuvre while it has been shown that a simple COTS stepper motor can be utilized in practice for manufacturing the GCMG concept. With three less DoF and reduced mathematical complexity, the GCMG design may be an effective alternative to the VSCMG in applications that is desired to explicitly follow a commanded trajectory. Even though the present study only investigates the GCMG cluster using the Moore-Penrose steering law, an extensive analysis is kept for future work to understand how the GCMG behaves to different steering laws.